Baryon stopping as a relativistic Markov process in phase space

We reconsider baryon stopping in relativistic heavy-ion collisions in a nonequilibrium-statistical framework. The approach combines earlier formulations based on quantum chromodynamics with a relativistic diffusion model through a suitably derived fluctuation-dissipation relation, thus allowing for a fully time-dependent theory that is consistent with QCD. We use an existing framework for relativistic stochastic processes in spacetime that are Markovian in phase space, and adapt it to derive a Fokker-Planck equation in rapidity space, which is solved numerically. The time evolution of the net-proton distribution function in rapidity space agrees with stopping data from the CERN Super Proton Synchrotron and the BNL Relativistic Heavy Ion Collider.


I. INTRODUCTION
In relativistic heavy-ion collisions at the CERN Super Proton Synchrotron (SPS), the BNL Relativistic Heavy Ion Collider (RHIC), or the CERN Large Hadron Collider (LHC), the incoming baryons are being slowed down ("stopped") as they interpenetrate each other, while in the spatial region between the receding, highly Lorentzcontracted fragments [1] a hot fireball is formed, which cools during its expansion and eventually hadronizes in a parton-hadron crossover.Of particular interest is the initial stage of such a collision with the local thermalization of quarks and gluons, and the simultaneous stopping of the baryons.The latter occurs essentially through collisions of the incoming valence quarks with soft gluons in the respective other nucleus.
Various models to account for the stopping process and its energy dependence have been developed, for example, in Refs.[2][3][4] and related works, which are relying on the appropriate parton distribution functions and hence on quantum chromodynamics (QCD).These models yield agreement with the available stopping data at SPS and RHIC, such as the distributions of net protons (protons minus produced antiprotons) in longitudinal rapidity space, and also provide predictions at LHC energies, where stopping data at forward rapidities are not yet available.They do not, however, provide the time development from the initial distribution at the instant of the collision to the final, measured one.
Complementary time-dependent approaches to stopping and local equilibration have relied on phenomenological, nonequilibrium-statistical approaches: A linear Fokker-Planck equation for the net-baryon rapidity distribution function had been proposed in Ref. [5], which accounts for the time evolution of the net-baryon or netproton rapidity distributions in a two-source relativistic diffusion model (RDM).Variants of the model with a nonlinearity in the diffusion term have subsequently been suggested in Refs.[6][7][8] and related works, which assume, however, the debatable validity [9] of nonextensive statistics.A linear diffusion model for particle production was also put forward in Ref. [10] and subsequent works.For produced particles at RHIC and LHC energies, a third (midrapidity) source is essential to cover pair production processes in the central fireball that provide the bulk of charged-hadron generation at sufficiently high energy [9].When considering net baryons or protons, however, the midrapidity source cancels out because it is equally composed of particles and antiparticles.The highly nonlinear local thermalization of quarks and gluons in the initial stages of the collision can be modeled through quantum Boltzmann-like collision terms, which require numerical solutions [11,12], but also a schematic model has been developed that accounts for the fast local equilibration through analytical solutions of a nonlinear partial differential equation [13].
Diffusion models are being used in many areas of physics, chemistry, and biology [14,15].They have originally been developed by Einstein and Smoluchowski to provide a mesoscopic theory of Brownian motion [16][17][18] as linear differential equations for the Brownian particles' single-particle distribution function.Alternative microscopic approaches treat the Brownian particles' trajectories as stochastic processes in position space, for example in the form of a Wiener [19] or Ornstein-Uhlenbeck process [20].Stochastic processes and stochastic differential equations have subsequently been considered in more generality in a newly established branch of mathematics, stochastic calculus, with notable contributions by Itô [21,22], Stratonovich [23], Fisk [24], and Klimontovich [25] who introduced various concepts of a stochastic integral, each with different mathematical properties and physical interpretations.In some cases, connections between the micro-and mesoscopic formulation can be established through a Kramers-Moyal expansion [26,27] or the Feynman-Kac formula [28].
Following the discovery of special relativity, it became clear that statistical physics in general, and diffusion models in particular, would have to be adapted [29] to meet the requirements imposed by a limited velocity of light.Especially, nontrivial Lorentz-invariant stochastic processes for spacetime coordinates are necessarily non-Markovian [30][31][32], which makes a straightforward generalization of the nonrelativistic diffusion equation to special relativity impossible.While a general stochastic process may depend on any finite or infinite number of its prior realizations, Markov processes [33] have no memory of the past apart from their current state, which greatly simplifies their mathematical treatment.Consequently, stochastic processes used in physical models often have the Markov property, in spite of being mathematically the exception rather than the rule.
As long as a process' memory is finite, i. e., only a finite number of its previous realizations affect the next value, it can be reformulated as a coupled system of multiple Markov processes through the introduction of additional variables [34].This has been used in Refs.[35][36][37][38][39][40] to formulate relativistic phase-space diffusion processes based on a generalized Ornstein-Uhlenbeck process for the Brownian particle's momentum.These processes are Markovian in phase space but lose the Markov property when expressed solely in spacetime coordinates.The authors also deduced associated relativistic Kramers and Fokker-Planck equations for the particles' phasespace and momentum distribution functions, and derived fluctuation-dissipation relations suitable for an isotropic thermal background.
In the present work, we aim to derive a nonequilibriumstatistical diffusion model for baryon stopping in rapidity space that is based on the key premises of the phenomenological RDM, but is constructed from a consistent approach with relativistic Markov processes in phase space and incorporates the QCD-based theory through a suitably adapted fluctuation-dissipation relation.The corresponding Fokker-Planck equation will enable us to account for the time evolution of the initial distribution functions from the onset of a relativistic heavy-ion collision to the final, measured distributions of net baryons or net protons in agreement with the available SPS and RHIC data.
The key assumptions for our nonequilibrium-statistical approach to stopping in relativistic heavy-ion collisions are presented in the next section, followed by the equations of motion in Langevin and Fokker-Planck formulation in Sec.III.Drift and diffusion terms are discussed in Sec.IV, as well as the expected stationary state derived from the earlier QCD formulation that allows us to formulate an appropriate fluctuation-dissipation relation that determines the course of the time evolution from the initial to the final net-proton rapidity distribution functions.The latter are compared with available stopping data from SPS and RHIC experiments in Sec.V.The conclusions are drawn in Sec.VI.

II. NET-PROTON RAPIDITY SPECTRA
We model baryon stopping as a diffusive process in rapidity space of the participating nucleons whose dynamics are governed by a-not necessarily thermalized-fluctuating background representing the quarks and gluons of the fragments.In this process, interactions between the partons and the background are assumed to prevail such that the nucleon distribution function reduces to a superposition of single-particle probability density functions.As the distribution functions of participant baryons are experimentally inaccessible, we consider only participant protons in our model and compare our result to the measured net-proton number density in rapidity space, i. e., the difference between the distributions of protons and antiprotons produced in the collision, which we expect to be reasonably close to the participant-proton distribution function.
To incorporate the spatial separation of the two nuclear fragments, we use a two-source ansatz [5] and completely disconnect the time evolution of particles originating from the forward-and backward-moving fragment through separate probability densities and fluctuationdissipation relations.Taking advantage of the symmetry of the system with respect to its center of momentum, we then write the net-proton number density in rapidity space dN p− p/dy in the system's center-of-momentum frame F as the superposition where N p− p denotes the net-proton number and ψ(t; ±y) dy the probability to find a participant proton from the forward-or backward-moving fragment, respectively, at time t with rapidity in [y, y + dy].

A. Initial state
Prior to the collision of the nuclei at some time t i , we assume the system to be in an initial state where each nucleus can be approximated by a zero-temperature Fermi gas with appropriate Fermi momentum p F .Then, the protons of each nucleus are distributed in the nucleus's rest frame F * according to the momentum-space probability density function which is given by a Heaviside step function Θ scaled by a normalizing factor.We determine the Fermi momentum through a simple potential well model, [41] Here, Z denotes the nucleus's proton number, V * the nuclear charge volume, and r * the nuclear charge radius, which we take from Ref. [42].Choosing the orientation of F * such that p 3 * is parallel to the beam axis, we define the cylindrical coordinates with the proton mass m.Boosting to F leaves the transverse degrees of freedom unaffected (γ ⊥ = γ ⊥ * , ϕ = ϕ * ) while the longitudinal rapidity coordinate is shifted by the 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 with the Fermi rapidity y F = asinh(p F /m). Fig. 1 shows ψ i as a function of y * for the isotope 208 Pb.The numerical values of y F for 197 Au and 208 Pb are 0.3134 and 0.3136, respectively, and differ only slightly in the fourth decimal place.
A more realistic description of the initial state is principally desirable (for example, with a finite temperature); however, the exact form of the initial probability density function hardly influences the later stages of the time evolution.Hence, normal or delta distributions are often used as convenient approximations for the initial state [43].

B. Final state
For t > t i , the system evolves in time, driven by the fluctuating background, until it reaches a final state at some time t f when the partonic interactions between the nuclei effectively cease due to their increasing spatial distance.Consequently, for comparison with experimental data 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 not important at this stage of the model, since it does only appear in products with other a priori unknown quantities; it is not an observable.

III. TIME EVOLUTION
The time evolution of the system will ultimately be governed by a Fokker-Planck equation for the singleparticle probability density function ψ, whose drift and diffusion coefficient functions will be derived in Sec.IV from the expected mesoscopic behavior.A similar evolution equation had previously been determined on a phenomenological level in the RDM [9] from comparisons with available SPS and RHIC data.In this work, we will derive it from the underlying particle dynamics to shed light on the different assumptions entering our model.

A. Langevin formulation
We describe the trajectories of the individual protons as relativistic stochastic processes in spacetime that are Markovian when expressed in phase-space coordinates [34,44].In the following, these stochastic processes will be designated with uppercase letters, while lowercase letters denote the corresponding coordinates in the system's center-of-momentum frame F. The equations of motion for the spacetime position X α (t) and energy and momentum P α (t) as a function of time t follow from a relativistic generalization [35][36][37][38][39][40] of the Ornstein-Uhlenbeck process [20] in phase space, dX α = P α /P 0 dt , (6a) where the Greek index α runs from 0 to 3 and the Latin indices i, k from 1 to 3. The momenta P i are driven by three independent standard Wiener processes W k [19] representing the fluctuating background, while the particle energy P 0 is fixed by the mass-shell condition, P 0 = m 2 + P 2 , where m denotes the proton mass.The interaction between particles and background is governed by the drift coefficients µ p i and diffusion coefficients σ p i ,k : The former represent directed, deterministic effects (for example, friction or pressure gradients) and determine the mean value of the stochastic process; the latter are connected to undirected, stochastic interactions (such as random particle collisions) and its variance.In general, they can be functions of all involved stochastic processes X α and P i .Here, however, we will assume that they depend on the momentum processes only.
If we let the 3-direction of the coordinate system coincide with the beam axis, we can replace P 3 with the stochastic process Y = atanh(P 3 /P 0 ) , (7) which corresponds to the (longitudinal) rapidity y.
The associated drift coefficient µ y and diffusion coefficients σ y,k can be related to µ p i and σ p i ,k through differential calculus.
To simplify the following computations, we will assume that the longitudinal drift and diffusion coefficients' dependence on the transverse degrees of freedom is negligible, i. e., ∂ p i µ y = ∂ p i σ y,3 = 0 and σ p i ,3 = σ y,k = 0 for i, k = 1, 2. Then, the Langevin equations for Y decouple from those of P 1 and P 2 , Further, we want to treat σ y,3 as a constant with respect to rapidity for now, since the nonconstant case entails some technical subtleties regarding discretization and interpretation of the Langevin equations [21][22][23][24][25]45].We intend, however, to address this issue in a forthcoming publication.
The choice of a constant diffusion coefficient is permissible here because the rapidity Y may assume any real value, and hence arbitrarily large changes by the driving Wiener process still result in a physically permissible state of Y .By contrast, if we were to formulate a stochastic process for the particle's velocity or Lorentz factor, the diffusion coefficient would necessarily have to be nonconstant to prevent superluminal motion by suppressing fluctuations that would lead the stochastic process to an unphysical state [46].

B. Fokker-Planck formulation
To obtain an equation for the time evolution of the single-particle probability density function associated with the particle trajectories discussed in the preceding section, we perform a Kramers-Moyal expansion [26,27] with respect to the longitudinal stochastic processes defined in Eqs.(8) and the transverse stochastic processes from Eqs. ( 6).As we have decoupled X 3 and Y from the other processes, we can immediately integrate out the transverse coordinates x 1 , x 2 and p 1 , p 2 , which leaves us with the Kramers equation for the marginal probability density function f of longitudinal position x 3 and rapidity y, with f (t; x 3 , y) dx 3 dy giving the probability to find a participant proton at time t with X 3 ∈ [x 3 , x 3 + dx 3 ] and Y ∈ [y, y + dy] in F. To ease notation, we drop the subscripts of the longitudinal drift and diffusion coefficients from now on as they are the only coefficient functions left.
Given an appropriate initial condition, we could in principle solve Eq. ( 9).However, since the position coordinate x 3 is unobservable, we integrate it out, thus reducing Eq. ( 9) to a Fokker-Planck equation for the marginal rapidity probability density ψ, where we have used that f must vanish at the boundaries and that µ and σ were assumed to be independent of x 3 .Alternatively, Eq. ( 10) can be rewritten as a continuity equation with the probability current density which can be decomposed into an advective (j a ) and a diffusive part (j d ), j a (t; y) = µ(y)ψ(t; y) , (14a) In this context, the prefactor σ 2 /2 can be understood as the protons' diffusivity D in rapidity space.When defining nonlocal observables as in Eq. ( 11) in relativistic statistical physics, care has to be taken [47] because the involved integral introduces a dependence on the chosen hypersurface.In our case, integration is done with respect to isochronous hypersurfaces in F; we expect this to give a reasonable representation of the dN p− p/dy measuring process.A more accurate treatment would require precise knowledge of the particle positions and detector 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 and diffusion coefficients, apart from setting the latter constant with respect to rapidity.Instead of deriving them from microscopic considerations, we will set the coefficients in a way that the solutions of the Fokker-Planck equation ( 10) reproduce a certain expected mesoscopic behavior of the physical system to be modeled, as proposed in Refs.[35,40].Possible choices include presetting the system's stationary state or specifying the time evolution of some macroscopic observable.Generally, two such criteria are needed to uniquely determine both coefficients [39], however, having set σ 2 /2 = D to a constant that can be numerically deduced by fitting the model to experimental data, one constraint will suffice in our case.In an earlier version of the RDM, a linear approximation was used for the drift coefficient function that enabled an analytical solution of the Fokker-Planck equation [5].

A. Expected stationary state
The stochastic process defined in Eq. (8b) would approach a stationary state if its time evolution continued past t = t f .We can estimate this state by assuming the formation of a color-glass condensate (CGC) [48][49][50][51], a coherent state based on the saturation of the gluon density below a characteristic momentum scale Q s .In the CGC framework, the post-collision distribution of the forward-moving participant protons is given by [52][53][54] where x is the longitudinal momentum fraction carried by the protons' valence quarks and q v denotes the valencequark distribution function for which we use the NNLO results from [55].C is a normalizing constant that sets the integral of ψ CGC to unity.To determine the distribution function g of the soft gluons from the backwardmoving fragment, we choose the Golec-Biernat-Wüsthoff model [56] in which g reduces to a simple function of the scaling variable The gluon-saturation-scale exponent λ determines the while the constant Q 2 0 sets its dimension and the mass number A its scaling with the nucleus's size.Together with the center-of-mass energy per nucleon pair √ s NN , the same three parameters determine the rapidity dependence of ψ CGC through the dimensionless function More details on the subject can be found in Refs.[2,3], where similar distribution functions were fitted directly to proton-stopping data without considering a time evolution of the system.
The integral in Eq. ( 15) has no analytic solution, and hence, we solve it numerically with adaptive Gauss-Kronrod quadrature in all our computations.For rapidities far from zero, however, analytical approximate solutions exist, which we will discuss briefly below.
For large positive rapidities, the argument of g becomes very small such that we can approximate g(ζ) ≈ 4πζ and separate the x-and y-dependent terms, With the integral yielding a constant numerical factor, the distribution function thus decays exponentially for large positive y, ψ CGC (y) = O exp(α + y) , with decay constant α + = −2(1 + λ).
For large negative rapidities, only small x values contribute to the integral due to the exponential damping with τ (y).If the low-x behavior of the valence-quark distribution is given by x q v (x) ∼ ax b , Eq. ( 15) reduces to the definition of the gamma function times an exponential function of τ (y), Accordingly, the distribution function exhibits an exponential tail also for large negative values of y, where

B. Fluctuation-dissipation relation
A Fokker-Planck equation of the form Eq. ( 10) possesses a unique stationary solution ψ s .It can be easily calculated by using the fact that its time derivative vanishes, ∂ t ψ s (y) = 0, resulting in [57,58] where the lower integration limit is chosen such that the integral exists.All solutions of Eq. ( 10) would converge against this state for t → ∞, lim t→∞ ψ(t; y) = ψ s (y), if we continued their time evolution past t f , which is, however, physically impossible since the fragments separate.Hence, fixing the drift coefficient and diffusivity determines ψ s and vice versa: Inverting Eq. ( 21) yields the fluctuation-dissipation relation associated with a given stationary state ψ s [35,39,40], If we then identify ψ s ≡ ψ CGC with the CGC distribution from Eq. ( 15), the drift coefficient µ can thus be fixed as a function of D and y.Like ψ CGC , the resulting expression for µ is not analytic, but can be evaluated numerically as shown in Fig. 2. The graph is roughly S-shaped and converges toward constant values for y → ±∞ due to the exponential tails of ψ CGC , Its zero crossing marks the peak position of ψ CGC ; the maximum close to y ≈ 0 indicates an inflection point of the logarithm of ψ CGC .

V. RESULTS
We obtain a dimensionless form of the Fokker-Planck equation (10) by substituting the time t with the evolution parameter s(t) = (t − t i )/(t f − t i ) and reordering some terms.The transformed equation reads with ∆t = t f − t i .While s, y, and ψ are dimensionless by definition, we have arranged the remaining quantities such that they form the composite dimensionless factors D∆t and µ(y)/D.The latter is given by the fluctuation-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 with the diffusivity D, while ∆t is defined as the time span during which the system is subject to the associated forces, the compound variable D∆t can be interpreted as the net impact of the partonic interactions between the nuclei.Appearing only on the right-hand side of Eq. ( 24), it can be completely absorbed into the evolution parameter 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 hence indicate that the system remains close to its initial state, while larger values drive it closer toward the stationary state imposed by the fluctuation-dissipation relation.The transformed Fokker-Planck equation ( 24) is solved numerically for 0 < s ≤ 1 by discretizing the rapidity derivative operators and solving the resulting system of ordinary differential equations with an additive Runge-Kutta method [59].
We compare the results of our calculations to SPS and RHIC data from the NA49 and BRAHMS Collaboration, respectively [60][61][62][63].The gluon-saturation-scale exponent λ and prefactor Q 2 0 , the net-proton number N p− p, and the diffusivity times elapsed time D∆t are free parameters 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 experimental data [60][61][62][63].For √ sNN = 200 GeV, final and stationary state were so close to each other that no meaningful fit result could be determined for D∆t (see text); accordingly, no numerical value is given at this energy.λ and Q 2 0 are shared parameters and hence take the same numerical values for all collisions.Reduced sums of squared residuals χ 2 /ndf (excluding shared parameters) are given for each setting as measures for the individual goodness-of-fit.[61] and 2008 [62], respectively; uncertainties are depicted as bars.

Nuclei
simultaneous weighted least-squares fit of the final netproton distribution functions to the experimental data, where minimization of the fit objective is done numerically with a quasi-Newton method [64,65].We restrict N p− p to deviate not more than 10 % from the respective Glauber result, while λ and Q 2 0 are treated as common parameters that take the same numerical values for all collisions in the SPS to RHIC energy region.Our results are given in Tab.I; the combined sum of squared residuals divided by the total number of degrees of freedom is χ 2 /ndf ≈ 0.89.The estimates for λ and Q 2 0 compare well to literature results, where λ ≈ 0.288 and Q 2 0 ≈ 0.097 GeV 2 were obtained in a fit to deep-inelasticscattering data from the DESY Hadron-Electron Ring Accelerator (HERA) [56].
The time evolution of the net-proton distribution functions for the two collisions with lower energy is shown in Fig. 3.As expected [66], the distribution functions converge exponentially in time toward the stationary state, which is indicated in the plot by a logarithmic spacing of the intermediate time steps.While the final distribution functions appear to differ only slightly from the stationary ones, the systems are still far from their stationary states in a temporal sense, as further convergence slows down exponentially.
At the lower center-of-mass energies √ s NN = 17.3 and 62.4 GeV, the final and stationary state differ enough for a reasonable estimate of D∆t, which takes values between 3 and 4. At the higher energy √ s NN = 200 GeV, however, final and stationary state are too close compared to the experimental errors.As a consequence of the exponential convergence in time, the uncertainty in the determination of D∆t becomes orders of magnitude larger than the actual value.Therefore, a meaningful estimate of D∆t is not possible and no values are given in Tab.I at this center-of-mass energy.Fig. 4 therefore shows only the stationary net-proton distribution functions for the two collisions at √ s NN = 200 GeV, which are nearly indistinguishable from the final distribution functions.A time evolution from the initial to the final state cannot be given due to the indeterminate value of D∆t.The net-proton numbers differ for the two centralities, being higher for 0-5 % and lower for 0-10 %, which is consistent with the latter data containing additional events with fewer participants.All required numerical routines were implemented with the Julia programming language [67]; functionalities for the solution of differential equations and parameter optimization were provided by the packages Differen-tialEquations.jl [68] and Optim.jl[69], respectively.

VI. CONCLUSIONS
A relativistic phase-space diffusion model for the time evolution of net-proton distribution functions in rapidity space was presented to account for the transition process from the initial to the final state in baryon stopping.Inspired by the phenomenological RDM, the model uses similar key assumptions, but is based on stochastic particle trajectories constructed from relativistic Markov processes in phase space that are equivalent to non-Markovian spacetime processes.The drift and diffusion coefficient, which carry over from the Langevin to the Fokker-Planck formulation of the system's time evolution, were determined by assuming a constant diffusivity in rapidity space and setting the stationary solution of the Fokker-Planck equation to a QCD-inspired distribution function.Due to the latter's exponential tails, the associated fluctuation-dissipation relation was found to be virtually constant for large absolute rapidities.Analytic expressions for the limiting values were derived.
A simultaneous least-squares fit was used to determine the free model parameters for four data sets recorded at SPS and RHIC by the NA49 and BRAHMS Collaboration, respectively.In the fit, the net-proton number N p− p was restricted to a neighborhood of the corresponding Glauber result for each collision.The gluon-saturationscale exponent λ and prefactor Q 2 0 were treated as common parameters taking the same value in all comparisons with experiment.No constraints, apart from positivity, were placed on the dimensionless factor D∆t composed of the diffusivity and the elapsed time between initial and final state.For 208 Pb and 197 Au collisions at √ s NN = 17.3 and 62.4 GeV, respectively, agreement with the data could be reached and an estimate of the time evolution from the initial to the final state was given.At 200 GeV, the latter was not possible, since the final and stationary distribution functions were found to be too close compared to the experimental uncertainties.
The phase-space diffusion framework adopted in this article is easily adaptable to different physical systems and allows to construct the drift and diffusion coefficient functions, which can be difficult to access theoretically, from mesoscopic considerations.In a forthcoming work, we will examine a possible application to charged-particle production in relativistic heavy-ion collisions.

FIG. 1 .
FIG. 1. Marginal rapidity probability density function ψi (solid, red) of the participant protons in a 208 Pb nucleus prior to the initial collision.For comparison, a normal distribution with zero mean and standard deviation yF/2 is depicted as dashed blue curve.Dotted vertical lines indicate the Fermi rapidity ±yF.