Diffusion processes involving multiple conserved charges: a first study from kinetic theory and implications to the fluid-dynamical modeling of heavy ion collisions

The bulk nuclear matter produced in heavy ion collisions carries a multitude of conserved quantum numbers: electric charge, baryon number, and strangeness. Therefore, the diffusion processes associated to these conserved charges cannot occur independently and must be described in terms of a set of coupled diffusion equations. This physics is implemented by replacing the traditional diffusion coefficients for each conserved charge by a diffusion coefficient matrix, which quantifies the coupling between the conserved quantum numbers. The diagonal coefficients of this matrix are the usual charge diffusion coefficients, while the off-diagonal entries describe the diffusive coupling of the charge currents. In this paper, we show how to calculate this diffusion coefficient matrix from kinetic theory and provide results for a hadron resonance gas and a gas of partons. We further find that the off-diagonal entries can reach similar magnitudes compared to the diagonal entries. In order to provide some insight on the influence that the coupling between the net charge diffusion currents can have on heavy ion observables, we present first results for the diffusive evolution of a hadronic system in a simple (1+1)D-fluid dynamics approach, and study different configurations of the diffusion matrix.


I. INTRODUCTION
The main motivation for studying nuclear collisions at relativistic energies is to understand the properties of strongly interacting matter. Especially the possibility of observing the transition from hadronic matter to quark-gluon plasma (QGP), as predicted by Quantum Chromodynamics (QCD), has been the driving force behind the active experimental heavy-ion programs at the Brookhaven National Laboratory (BNL) and the Conseil Européen Pour La Recherche Nucléaire (CERN). During the last couple of decades, the high-energy nuclear collision experiments performed in the Relativistic Heavy Ion Collider (RHIC), at BNL, and in the Large Hadron Collider (LHC), at CERN, have shown that a considerable amount of QCD matter is produced in these collisions and that it is possible to infer the properties of such matter from the experimental data.
A dissipative process that is usually neglected, is the diffusion of conserved charges due to temperature or density gradients. Diffusion is a dissipative process which occurs as soon as inhomogeneities arise in a conserved quantity. In the simplest non-relativistic case, the diffusion current j q of charge q is described through Fick's first law [29,30], where the current is generated by gradients in net charge density n q (t, x) and the diffusion coefficient D q characterizes the reaction strength of this thermal force. In the highest-energy nuclear collisions, the created matter has almost zero net baryon density at midrapidity, and the effects of diffusion are expected to be small in this region [31]. However, diffusion is expected to play an increasingly important role as the net baryon density increases with the decreasing collision energy.
Recently, the Beam Energy Scan (BES) program was initiated at RHIC. In this program, nuclear collisions were systematically performed at lower energies in order to investigate the phase diagram and transport properties of nuclear matter at finite net baryon (and net electric charge) densities [32][33][34]. At beam energies down to, e.g., √ s NN = 7.7 GeV in the RHIC BES, the baryon chemical potential can reach values up to µ B ∼ 400 MeV, which is significant when compared to the temperatures that are reached in this system [35,36]. Furthermore, the Facility for Antiproton and Ion Research (FAIR) at the Gesellschaft für Schwerionenforschung (GSI) in Darmstadt, Germany, and the Nuclotron-based Ion Collider facility (NICA) in Dubna, Russia, aim to generate and study compressed hadronic matter at large baryon densities [37]. The theoretical description of those collisions could rely on diffusion dynamics.
The constituents of strongly interacting matter carry a multitude of conserved quantum numbers: baryon number, strangeness, electric charge, among others. As a result, the diffusion currents of the conserved charges must be coupled with each other. This multicomponent nature of diffusion in strongly interacting matter was first fully embraced in Ref.
[1], where the full matrix of diffusion coefficients was computed, and it was subsequently found that the diffusion coefficients, describing the cross-coupling between the diffusion currents, are of the same magnitude as the "normal" (diagonal) diffusion coefficients. The purpose of this study is to complement Ref.
[1] and provide more details on the computation of the diffusion matrix for strongly interacting matter, as well as to provide an initial hydrodynamic calculation that illustrates the influence of the cross-couplings in relativistic nuclear collisions. As we will show, a novel phenomenon emerging from the coupling is a generation of regions of non-zero net strangeness from initially net strangeness neutral matter.
This work is organized in two parts. In the first part we discuss the diffusion coefficients and in the second part we present a first investigation with fluid dynamics. In Section II we define the most important notations and expressions used in the paper. Section III provides a short review of diffusion in a relativistic gas with multiple conserved charges and introduces the diffusion coefficient matrix, which characterizes the coupling of the diffusion currents. We present the derivation of the diffusion coefficient matrix within a linear response approach from relativistic kinetic theory in Section IV and we further discuss its properties and results in relaxation time approximation (RTA) in Section V. The first part of this work is concluded with detailed discussions of the results for the coefficient matrix for a hadronic and a massless partonic system in Sections VI and VII. In Section VIII we provide a short overview of the fluid dynamic approach used and also present our first results for the longitudinal diffusive evolution of a hadronic system. A summarizing conclusion and an outlook is provided in Section IX. We use natural units, = c = k B = 1, and greek indices run from 0 to 3.

A. Basic definitions
Throughout this paper, we will express the momentum as k µ and the coordinates as x µ . We denote the metric as g µν and impose the (+, −, −, −)-signature. It is convenient to express all tensors in terms of irreducible tensors regarding the local fluid velocity, u µ ≡ u µ (x). Therefore, we introduce the orthogonal projectors ∆ µ ν ≡ g µ ν −u µ u ν and ∆ µν αβ ≡ 1 2 ∆ µ α ∆ ν β + ∆ µ β ∆ ν α − 1 3 ∆ µν ∆ αβ . The projectors are symmetric (∆ µν = ∆ νµ = ∆ (µν) and ∆ µν αβ = ∆ (µν) (αβ) ) and are orthogonal to the fluid velocity (u ν ∆ µ ν = 0 and u α ∆ µν αβ = 0). More details can be found in Refs. [38,39]. We denote the projected tensors as A µ ≡ ∆ µ α A α and A µν ≡ ∆ µν αβ A αβ . The four-derivative can then be decomposed into the comoving derivative D = u ν ∂ ν and the projected derivative or gradient ∇ µ = ∆ ν µ ∂ ν : For later use, we define the particle energy in the local rest frame (LRF) as where the index i refers to the particle's species. The state of the system is characterized by the single-particle distribution function of each particle species, f i (x, p). It can be decomposed into an equilibrium part, f (0) i,k , and an off-equilibrium part, δf i,k , as f i (x, p) = f (0) i,k + δf i,k . We introduce the following notation for the integration measure: The momentum integrals over the distribution functions will be expressed using the following notation: B. Kinetic theory The evolution of f i (x, k) =: f i,k is given by the Boltzmann equation, where C ij is the collision term. The energy-momentum tensor T µν and the net charge currents N µ q are expressed as the following momentum integrals of the single-particle distribution function [38] and they fulfill the local conservation laws: ∂ ν T µν = 0 and ∂ µ N µ q = 0. It is convenient to decompose T µν and N µ q in terms of the fluid velocity field, u µ . Without loss of generality, we use Landau's definition of the fluid velocity [40], where u µ is an eigenvector of T µν with an eigenvalue given by the energy density in the local rest frame of the fluid (LRF), . That is, T µν u ν = u µ . The decompositions read where we introduced the local isotropic equilibrium pressure P 0 in the LRF, the bulk viscous pressure Π in the LRF, the shear stress tensor π µν , the net charge densities n q with q ∈ {B, Q, S} in the LRF, and the corresponding net charge diffusion currents j µ q . The bulk viscous pressure, the shear-stress tensor and the diffusion currents represent the dissipative corrections in the energy-momentum tensor and the four-currents of the charges. The diffusion currents of the net charges q are the main objects of our investigation and represent the charges diffusing orthogonally to the flow of the fluid. In this scheme, each introduced quantity can also be expressed as a contraction of the currents, T µν and N µ q , with u µ and ∆ µν , By specifying an equation of state, we can define the temperature and the chemical potentials for this system using the traditional matching conditions [40], where eq and n q,eq are the energy density and net charge densities of the system in local thermodynamic equilibrium, respectively. These quantities are calculated in kinetic theory by introducing the local equilibrium distribution function. In this work we want to restrict ourselves to classical statistics, and therefore the equilibrium distribution is given by the Maxwell-Juettner function where µ i = B i µ B + Q i µ Q + S i µ S is the chemical potential and g i is the degeneracy of the i-th species. Furthermore, the local equilibrium pressure is determined by the temperature and chemical potentials,

III. NET CHARGE DIFFUSION
In order to describe diffusion processes in relativistic fluids, a relativistic version of Fick's law must be employed. For a fluid with only one conserved charge, q, the relativistic Fick's law reads [40,41]: where the diffusion current is generalized to be generated by a gradient in the corresponding thermal potential of the charge α q ≡ µ q /T = β 0 µ q , and β 0 = 1/T is the inverse temperature. Note that (in flat Minkowski space) in the local rest frame ∇ µ ≡ (0, − ∇), and because of the sign, diffusion currents dissipate the existing inhomogeneities that originally generated the current. Often, instead of the charge diffusion coefficient, κ q , the corresponding charge conductivity, σ q ≡ κ q /T , is used. We can relate κ q to D q (introduced in Eq. (1)) by evaluating and imposing that the temperature is homogeneous, ∇ µ β 0 = 0, leading to As already stated, there are multiple conserved charges in nuclear matter: the baryon number, strangeness and electric charge. Moreover, the constituents of quark and hadronic matter carry multiple types of these charges, e.g. the proton carries baryon number and electric charge while the hyperons carry strangeness, baryon number and electric charge. Therefore, these constituents must react to multiple types of gradients in charge chemical potentials, in such a way that a gradient in baryon number does not only generate a baryon current, but can also produce currents in strangeness and electric charge (depending on the chemistry of the system). In order to account for this coupling, we introduced the diffusion coefficient matrix in Ref.
[1], which relates the charge diffusion currents to gradients in all thermal potentials, α q , as The objective of this paper is the evaluation and a first look at the possible dynamic implications of the complete diffusion matrix. In the first part of this work we will present a method of computation from relativistic kinetic theory.

IV. LINEAR RESPONSE THEORY: FIRST-ORDER CHAPMAN-ENSKOG EXPANSION
In this chapter, we present a method of evaluating the full diffusion coefficient matrix in relativistic kinetic theory. Here, we follow Refs. [1,25] and consider a dilute gas consisting of N species particle species, with the i-th particle species having degeneracy g i , electric charge Q i , strangeness S i and baryon number B i . The system shall be under the influence of spatial gradients in baryon, strangeness and electric chemical potentials over temperature ∇ µ (µ q /T ) (with q ∈ {B, Q, S}), but no other external forces, as assumed in Ref. [25]. The gradients are assumed to be small, such that the distortions from (local) equilibrium are small and linear response theory is applicable.

A. The linearized Boltzmann equation
We consider the system to be initially in global equilibrium. Next, we apply small gradients in the chemical potentials that are instantly switched on and cause a small perturbation δf i p of the single-particle distribution function (of the i-th particle species) from equilibrium. This perturbation generates a diffusion current in the corresponding charges. The aim of this section is to set up the Boltzmann equation for this situation.
The magnitude of gradients can be quantified by introducing the so-called Knudsen numbers, Kn, which are constructed as ratios of the characteristic microscopic and macroscopic length scales, Kn = micro / macro . Thus, Knudsen numbers are small if the corresponding macroscopic length scales, macro , are large in comparison to the microscopic length scales, micro . The later is often taken to be the mean free-path of a particle in the gas, and macro is related to the gradients in the system. If the gradients that generate the perturbation of the single-particle distribution are small, it may be possible to expand the distribution in terms of the Knudsen number and truncate such an expansion at lower order: This expansion is also referred to as Chapman-Enskog expansion [42]. If the Knudsen number is sufficiently small and the series defined above converges, it is possible to neglect contributions that are of second order or higher and the perturbed single-particle distribution function can be approximated solely in terms of its first order terms, Applying the Chapman-Enskog expansion to the Boltzmann equation (7) and only retaining the terms that are of first order in the Knudsen number leads to the following equation: where we introduced the linearized collision term and σ ij→ab (s, Ω) is the differential cross section for the binary interaction of incoming particles of species i and j, with outgoing particles of species a and b (denoted as ij → ab), at the center of mass collision energy √ s in a solid angle Ω. Further, we introduced the symmetry factor γ ij = 1 − 1 2 δ ij and the notationf i,k , where a = 1 for fermions, a = −1 for bosons or a = 0 for classical particles. In this paper, we limit our discussion to classical statistics and thereforef (0) i,k = 1, and to binary elastic processes, where only processes ij → ij are considered (σ ij→ij ≡ σ ij ). The achieved equation is typical for perturbation theory: the perturbed quantity on the right-hand side is determined by unperturbed quantities on the left-hand side of the equation. The left-hand side in the linearized Boltzmann equation (19) can be evaluated by first decomposing the four-derivative into comoving time derivative and projected derivative, and then substituting the comoving time derivatives of the primary fields , n q (or β 0 ≡ 1/T and α q ≡ µ q /T correspondingly) and u µ using the explicit local conservation laws from ideal fluid dynamics, Above, we introduced the expansion scalar θ ≡ ∇ µ u µ . Further, using the Euler relation, and the Gibbs-Duhem relation, in the form, we find the following equivalent form to the momentum conservation equation: Following this procedure, we derive the following source term [43] for a system with multiple conserved charges (terms related to the shear-stress tensor and bulk viscous pressure are omitted in the last line) where we defined the shear tensor σ µν ≡ ∂ µ u ν . The diffusion current is then calculated as, Thus, the correction f (1) i,k related to net-charge diffusion can be calculated from the following linear equation, by inverting the linearized collision term, The source term can be understood as a force term that generates the perturbation of the single-particle distribution due to the gradients in the thermal potentials, ∇ µ α q , which will eventually give rise to diffusion currents in the conserved charges q, according to Eq. (25).
Applying this to the linearized collision term (20) in the classical limit and for elastic scatterings only results in We can rewrite the linearized Boltzmann equation into a matrix equation by multiplying Eq. (26) with E n−1 i,k k ν i and then integrating over the momentum k i , such that we evaluate orthogonal moments of the Boltzmann equation. This becomes even more apparent when we realize that k µ i fulfills the following orthogonality relation for arbitrary scalar functions F in energy. Furthermore, because the gradients in the thermal potentials, ∇ µ α q , are arbitrary, we can split the linearized Boltzmann equation into separate equations for each charge q. Using the abovementioned evaluation of moments, the orthogonality properties of the momentum basis (29), and the separation by charge, we arrive at the following set of linear equations [1]: where we introduced the abbreviations and we used the dispersion relation k i, α k n are the entries of an (N species ·M )-dimensional source vector and λ (j) q,m are the entries of an (N species ·M )-dimensional vector of the expansion coefficients from Eq. (27), which are the solutions of the linear set of equations. In order to make matrix M quadratic, we set the parameter n to run from 0 to M . Furthermore, there are as many sets of such matrix equations (30) as there are considered charge types. In this paper, we limit ourselves to baryon number, strangeness and electric charge, and therefore there are three sets of linear equations to solve.
Up until this point, all steps were done without imposing the definition of the local rest frame. As already stated, in this work, we use the Landau definition of the four-velocity [40], in which all orthogonal energy-momentum flow vanishes: Applying the expansion (27) i,k to Eq. (32) gives us an additional constraint for the expansion coefficients, Together with the matrix equations (30), Eq. (33) forms a set of linear equations of which the expansion coefficients λ q,m are the solutions. By applying the expansion in Eq. (27) to the diffusion current (25) and directly comparing to its Navier-Stokes form, we arrive at an explicit form for the entries of the diffusion coefficient matrix: In Ref. [25], the quick convergence of the series in Eq. (35) was demonstrated and thus, we restrict the linearized calculations to the truncation order M = 1 in this paper.

V. RELAXATION TIME APPROXIMATION
In order to verify the results of the above calculations, it is useful to compare the results with a simple, analytic estimate. We apply the relaxation time approximation (RTA), in which the collision term is assumed to take a simple form where τ is the relaxation time. The relaxation time can be interpreted as a global mean free-time between collisions of particles, and is an input parameter.

A. Diffusion coefficients in RTA
Applying the RTA to the linearized Boltzmann equation (26) allows us to directly identify its analytical solution for the perturbation The diffusion currents then take the form and by direct comparison we arrive at the RTA expression for the diffusion coefficients: This expression can also be written as where the last integral gives the partial equilibrium pressure P 0i of particle species i that, in the Boltzmann gas, can be written as P 0i = n i T , where n i is the total number density of the particle species. Thus, the expression for the diffusion coefficients becomes Even if derived in the RTA, this expression allows us to identify the main features of the diffusion coefficients, in particular, its temperature dependence. We first note that in Eq. (41) the symmetry of κ qq [44,45] with respect to charge q ↔ q is explicit. Moreover, we note that the charge combination q i q i is the same for a particle and its corresponding anti-particle, so that the first term increases as the total density of charge carriers increases, with the largest contribution coming from the lightest particle species that carries both charges q and q . The last term is proportional to the net charge densities (particle minus antiparticle), and it becomes important when the net charge densities are comparable to the total charge density. Furthermore, the relaxation time is related to the inverse of the scattering rate, τ ∼ 1/Γ scatt. ∼ 1/(n tot σ tot ), and thus we deduce that the diffusion coefficients are suppressed by the scattering rate of the charged particles, which is strongly related to the total particle density of the medium with which they scatter. The dependence of κ qq on temperature and µ B is discussed below in more detail when we show the results for hadron gas.

B. Ultrarelativistic limit
In the ultrarelativistic limit, all coefficients for fixed chemical potentials, µ q = const., have the same asymptotic limit. In order to show this, we first make use of the fact that in the case of massless, classical particles, the thermodynamic integrals simplify which leaves us with the expression for the diffusion coefficients in the massless case in the RTA For fixed chemical potentials, all mass and chemical scales can be neglected in the ultrarelativistic limit, since m i /T → 0 and α i = µ i /T → 0 for all particle species. Because all thermal potentials vanish, all net charge densities n q also vanish in this limit. The high temperature limit follows directly from the massless limit expression (43) with exp(α i ) → 1, and thus reads The relaxation time can be related to the scattering rate as where n tot = i E i,k i,0 is the total particle density, σ tot the total averaged cross section for the interaction between the constituents of the gas, and C is a constant. For the massless case, we can then write that and in the limit introduced above, this simplifies to With this, the diffusion coefficients in the ultrarelativistic limit read: which is independent of any chemical potential. Furthermore, it becomes a constant if the total cross section, σ tot , is constant, while, in the conformal limit (where the cross section must scale with the temperature as σ tot ∼ 1/T 2 ) κ qq /T 2 becomes constant. These are properties that we also found to be true in the full linearized numerical evaluation.
In recent publications [46,47], the authors took the baryon diffusion coefficient of a massless QGP in RTA to be , where the relaxation time was assumed to be τ = C B /T . This relation is also a special case of Eq. (43) for the case of a massless gas with a particle and a corresponding anti-particle species with baryon charge B = ±1: C. Validity of the relaxation time approximation In this section, we show that the relaxation time approximation retains the correct scaling behavior in temperature and baryon chemical potential for constant cross sections. In order to investigate when the RTA is applicable, we compute the baryon diffusion κ BB for the lightest 19 hadron species (listed in Appendix X A) with a constant isotropic cross section (10 mb), using both the linearized collision term (later denoted as "full" in Fig. 1), Eq. (28), and its relaxation time approximation (RTA), Eq. (36). To this end, we employ the transport relaxation time τ tr ((45) with C = 3 2 ), This form originates from the assumption of a constant differential cross section, dσ(ϕ,ϑ) dϕdϑ = const., which is weighted at large scattering angles: where The comparison is shown in Fig. 1 for several values of baryon chemical potential µ B . We note that the temperature dependence of the full calculation is reproduced well by the RTA evaluation. Additionally, we remark that the µ Bdependence vanishes at high temperatures and that a (1/T 2 )-scaling of κ BB /T 2 is achieved at very high temperatures (which is not shown in Fig. 1), which we discussed in the last section for the case of constant cross sections. The full calculation deviates by a factor of 1 − 3 from the RTA and improves with larger temperatures accordingly. Finally, we conclude that the RTA becomes more reliable at higher temperatures, but any quantitative study should retain the full collision term as we have done; especially if non-constant cross sections (which in general introduce additional dependencies on temperature and chemical potential) are present.

VI. DIFFUSION COEFFICIENT MATRIX OF A HADRONIC GAS
We provide results for the diffusion coefficient matrix computed for a gas of hadrons and characterized by elastic binary hadron-hadron collision cross sections. We model the hadron gas using the most dominant mesons and baryons in a hot gas close to the QGP phase transition, that is, pions, kaons, nucleons, as well as lambda-and sigma-baryons (for particle properties see Appendix X A). From the particle data group [48], we use all available elastic, isotropic cross sections and complement other theoretically-described resonant cross section from GiBUU [49] and SMASH [50], as shown in Fig. 2. All missing cross sections are approximated by the constant values taken from UrQMD [51,52], or (approximated from) GiBUU [49], as given in table II in Appendix X A.
In the following, we present and discuss results for calculations completed in this particular example of a hadronic system, where we assumed a temperature range of T = 60 MeV to 180 MeV and a baryon chemical potential range of µ B = 0 to 600 MeV. The electric chemical potential is set to zero, µ Q = 0, for simplicity. Furthermore, in this section, when we show the transport coefficients we always set the net strangeness density to be zero, n S = 0 (as is expected to occur in the initial stages of heavy-ion collisions). This condition results in a strangeness chemical potential that cannot be zero, but must be a function of temperature and baryon chemical potential, µ S = µ S (T, µ B ). However, as we will see later, the cross-coupling between the diffusion currents can dynamically generate regions of non-zero net strangeness during the (fluid) dynamical evolution, even if it is initially zero. Therefore, in order to perform fluid dynamical simulations where the diffusion is taken fully into account, it is necessary to compute the full table of diffusion coefficients with arbitrary combinations of temperature and chemical potentials. Here, we show the positive baryon chemical potential range, however, we emphasize that there is in general no symmetry along the individual µ q axes. The coefficients are only symmetric under point reflections: if the sign of all chemical potentials are changed simultaneously. Further, we emphasize that there is no phase transition included in this model because this is a evaluation from (weakly-coupled) kinetic theory.
Due to the systematic uncertainty in the cross sections, we vary all approximated constant values in Tab. II in Appendix X A by multiplying them by a factor k = 0.5, 1, 2. We show this uncertainty of the diffusion coefficients by transparent bands in Figs. 3a -8a. Furthermore, in Figs. 3b -8b we show the full T and µ B dependence of the diffusion coefficients in 3D plots.
As discussed in Sec. V, the magnitude of the diffusion coefficients depends on both: the total and net density of the corresponding charge carriers, as well as their scattering rate. In turn, the scattering rate depends on the total particle density (with pion as the lightest hadron giving the largest contribution) and on the scattering cross section of the charge carriers. Many of the cross sections show quite a strong dependence on the particle collision energy (see Fig. 2) and this also is somewhat reflected in the temperature dependence of the scattering rates. Thus, the temperature and chemical potential dependence of the diffusion coefficients is a result of the interplay between the energy and hadron type dependence of the scattering cross sections, as well as the temperature and chemical potential dependence of the hadron densities.

A. Baryon diffusion
First, we begin with the diffusion of baryon number, which concerns the diffusion coefficients κ BB , κ BS and κ BQ in the coefficient matrix. They measure the response of baryon number due to gradients in baryon-, strangeness-and electric-chemical thermal potential π +/0/-+π -/±/+ π 0 +p, π 0 +n π + +p π -+p π -+n p+p p+p n+p n+p K + +p K -+p K -+n K+π Figure 2: The isotropic resonance cross sections from the particle data book [48], GiBUU [49] and SMASH [50] that we use for the computation of the diffusion coefficients. All combinations of species not listed are approximated by constant values (see Appendix X A). The grey bars denote the minimal √ s for the corresponding interaction. This plot was taken from Ref. [1].
All the diffusion coefficients in the baryon sector, which are shown in Figs. 3a, 3b, 4a, 4b, 5a, and 5b, display a rather strong dependence on temperature and baryon chemical potential. In particular, at µ B = 0, the magnitude of the coefficients increases by a factor of ∼ 10 4 , within the studied temperature range (note that the diffusion coefficients are divided by T 2 in the plots), and a similar increase is observed at lowest temperatures, within the studied µ B range. In both cases, the fast increase can be attributed to the rapid increase of the total baryon density. With nucleons being the lightest baryon number carriers, the baryon density is strongly suppressed by the Boltzmann mass factor exp(−m/T ) and the relative increase of the baryon density with temperature is much faster than the increase of the total particle density, which is mainly determined by pions. Therefore, at µ B = 0, the increase of baryon density clearly wins over the increase in the scattering rate as temperature increases, and this results in the strong temperature dependence of the diffusion coefficients at low temperatures. At non-zero µ B , the rapid increase of density with temperature is tamed by the fugacity factor exp(µ B /T ) and the temperature dependence becomes considerably milder at µ B = 600 MeV.
Similarly, the increase in µ B mainly affects the baryon density, with pion density being unaffected, and thus results in the fast increase of the diffusion coefficients with increasing µ B . When µ B becomes sufficiently large, the baryon density is comparable to the total hadron density, and thus begins to affect the scattering rate, and the increasing scattering rate also limits the growth of the diffusion coefficients. At the same time, the net density of baryons also becomes comparable to the total density and this further limits the diffusion coefficient, see Eq. (41).
All the diffusion coefficients of the baryon sector (κ BB , κ BS and κ BQ ) display very similar behavior. This is due to the fact that the lightest hadron that contributes to these coefficients are nucleons for κ BB and κ BQ , and hyperons for κ BS , with similar masses, and thus very similar behavior of densities. We further note that the qualitative behavior of κ BB is very similar in the test case with constant cross section (see Fig. 1). In this case, the energy and hadron type dependence of the cross sections only play a small role. The most visible difference is that with constant cross section, κ BB /T 2 at µ B = 600 MeV actually decreases with increasing temperature, whereas this does not happen with more realistic cross sections.
In the case of the baryon diffusion coefficient, κ BB , the shown bands demonstrate that the multiplicative factor in front of the constant cross sections changes the results more strongly at high temperatures, where the constant cross sections dominate the interactions in our study. Contrary to this, in the results for κ BS we see that these bands have a large width over the whole temperature range. This is because most of the assumed interactions of hyperons were modeled with constant cross sections and these are the only charge carriers contributing to this particular coefficient due to the charge combination. We also note that κ BS is negative due to the definition of the strangeness: hyperons with positive baryon number carry negative strangeness. This indicates a possible anti-correlation of baryon number and strangeness in dynamic simulations. In Fig. 5a we find for κ BQ that throughout most of the temperature range, the given bands are narrow.     3D-plot of the same coefficient over temperature and baryon chemical potential.

B. Strangeness diffusion
The diffusion of strangeness is characterized by the coefficients κ SS , κ SB and κ SQ via The κ SB -coefficient was already discussed because the diffusion matrix is symmetric and κ BS = κ SB [44,45]. For both of the remaining coefficients, κ SS and κ SQ , shown in Figs. 6a, 6b, 7a, and 7b, the lightest hadron that carries strangeness, or both strangeness and electric charge, is the kaon. Similarly to the baryons, the increase in temperature leads to an increase in total strangeness density compared to the total density determining the scattering rate, but in this case the effect on the diffusion coefficient is much weaker since kaons are significantly lighter than the lightest baryons. Further, there is no significant dependence on the baryon chemical potential since the kaons do not carry any baryon charge.

C. Electric diffusion
The response of electric charge due to gradients in all thermal charge potentials, α q , is measured by the coefficients κ QQ , κ QB and κ QS The only coefficient left to discuss is κ QQ , which is shown in the same manner as before in Figs. 8a and 8b. Contrary to κ QB , κ QQ again shows no significant µ B -dependence, since the most dominant charge carriers, the pions, do not carry any baryon charge. The fact that the lightest electric-charge carriers are also the lightest hadrons results in the fact that the total density of charge carriers grows at the same rate as the total hadron density, which in turn, determines the scattering rate. In this case, κ QQ depends very weakly on temperature, and the shown temperature dependency in Fig. 8a is from the 1/T 2 scaling.
We remark that for µ B = 0, the electric diffusion coefficient, κ QQ , coincides with the electric conductivity calculated in Ref. [25], κ QQ (µ B = 0) = T σ el (µ B = 0). The similarity of electric conductivity and diffusion (or in the Eckart frame [41] the heat conductivity) is the manifestation of the Wiedemann-Franz law [53]. Similar to the electric conductivity, the electric diffusion coefficient also decreases strongly with temperature and only shows a mediocre dependence on µ B [25]. This is because the dominant electric charge carriers are the pions, but a significant amount of baryonic species also contribute to the electric diffusion current. Similar to Fig. 3a, the band widths vary strongly with increasing temperatures where the constant cross sections dominate the interactions. The computation of the hadronic diffusion coefficients presented above can only be extended in temperature up to T ∼ 160 MeV where the transition to QGP is expected to happen. In order to extend the computation to higher temperatures, we need to complement the hadronic part by calculating the diffusion matrix also for the QGP. Here, we consider a toy model for the QGP where it is described as a massless gas of quarks and gluons undergoing isotropic elastic binary collisions. To this end, we take gluons and the three lightest quark flavors u, d, and s, together with their antiparticles. The degeneracy factors accounting for the spin and color degrees of freedom are g = 6 for quarks and g = 16 for gluons.
The magnitude of the diffusion coefficients is then determined by the collision cross section, which we must specify, and we do so in a very simplistic manner in order to get a first estimate. We can either fix the total cross section to a constant value, e.g. σ tot = 10 mb, or set the shear viscosity to entropy ratio η/s to a fixed value. The latter choice is more suitable, since shear viscosity of the QGP is often extracted from experiment making this assumption. The shear viscosity is given by η = 2 /(5n 0 σ tot ) and the entropy density in chemical equilibrium is s = 4n 0 [54,55]. Using the theoretical minimum η/s = 1/(4π) [56], the total isotropic cross section can be fixed to σ tot ≈ 0.716/T 2 [54,55].
The quarks carry baryon number, strangeness and electric charge 1 , and the gluons contribute to the diffusion coefficients, mainly through the scattering rate. In the case of s-independent, isotropic cross sections, the diffusion coefficients scale with the total cross section σ tot . As an example, we give the massless limit of the complete diffusion coefficient matrix at vanishing chemical potential, µ q = 0 for q ∈ {B, Q, S}: As discussed in Section V B, in the case of constant cross sections and µ i = 0, the diffusion coefficients κ qq are also constant in the ultrarelativistic limit and the scaled coefficients κ qq /T 2 therefore scale with inverse temperature squared. Contrary to this, in the conformal limit where η/s = const. and σ tot ∼ T −2 the scaled coefficients are constant over temperature. At non-zero µ q , all the coefficient acquire temperature dependence through the Boltzmann factors exp(µ q /T ) in the densities. However, at fixed µ q the temperature dependence vanishes at large temperatures. In Fig. 9 we plot the full diffusion coefficient matrix, where we show all acquired results for the hadronic gas already discussed in Section VI and also the results for the massless simple QGP model. We show results by fixing µ B = 0, 300 and 600 MeV, the electric chemical potential to zero, µ Q = 0, and also the net strangeness density to zero n S = 0 as in the hadronic case. We then simply compare and present our results for the two models in one summarizing plot and switch the model at 160 MeV temperature 2 . We already note that on the QGP side there is very little dependence on temperature and baryon chemical potential, especially at large temperatures, as expected.
Surprisingly, for µ B = 0 the coefficients for the two different models almost match at T = 160 MeV, where the phase transition or the smooth crossover would normally occur. The only exception seems to be the κ BQ , where there is a large discrepancy between the hadronic model and the QGP model. In the latter case, the coefficient vanishes at µ B = 0 since the generated currents of the quarks exactly cancel out the net flow of electric charge due to symmetry. However, we speculate that the hadronic results for κ BQ will decrease in magnitude if more hadronic particles are included in the computation, and will thus reduce the discrepancy between both models.
We further compare to the holographic results for the diagonal entries of the diffusion coefficient matrix from Ref. [24] (grey dashed and dash-dot-dotted lines). These results approach the conformal limit at high temperatures and only show moderate µ B dependence, but the overall shape and magnitude seems to be qualitatively consistent with our results. However, the results for the simple QGP model with fixed η/s = 1/4π coincides with the conformal limit at very high temperatures, as it should. We note that the κ SB -coefficient has the same magnitude as the baryon diffusion coefficient, κ BB , but is negative. We further remark that the strangeness diffusion coefficient is the largest coefficient in magnitude.
We found that the off-diagonal entries of the diffusion coefficient matrix can reach similar magnitudes to the diagonal coefficients (which are usually considered). We therefore would expect significant corrections to the diffusion currents due to the mixing of charge types compared to approaches when parts of the diffusion coefficient matrix are neglected. Nevertheless, the phenomenological consequences are still not known. In the following, we take a first step in this direction and investigate the influence of the full diffusion matrix in a one-dimensional fluid dynamics approach. -κ SB /T 2 Figure 9: Complete diffusion coefficient matrix plotted over temperature and for the baryon chemical potentials µ B = 0, 300 and 600 MeV. We show results for the assumed hadronic system in the temperature range T = 50 to 160 MeV and for the simple QGP model for fixed η/s = 1/4π for temperatures above 160 MeV. We compare to the holographic results achieved in Ref. [24]. This plot was taken from Ref. [1].

VIII. HYDRODYNAMIC EVOLUTION
In the last sections we evaluated and discussed the diffusion coefficient matrix for a simple hadronic and (massless) partonic system and showed that in the chosen basis of charge definitions there are non-vanishing off-diagonal contributions arising from the fact that hadrons and partons can carry several different charges. The goal of this chapter is to provide initial investigations of its implications with the help of relativistic fluid dynamics. After providing a short review of our framework, we present the first results for the dynamic evolution of a system with multiple conserved charges. Here, we assume the same hadronic system as presented in Chapter VI as an example. More sophisticated studies will follow in the future.

A. Transient dissipative relativistic fluid dynamics
The foundation of fluid dynamics is the exact conservation of energy, momentum and the net quantum numbers (or charges) q. In the same fashion as in the last sections, we assume conserved baryon number B, strangeness S and electric charge Q. The local conservation equations of energy, momentum and net charge q can then be expressed in general (curved) spacetime as where we introduced the covariant derivative, (· · · ) ;µ , and the Christoffel symbols of the second kind, Γ µ αβ . The central assumption of fluid dynamics is that the evolution of the fluid is taking place close to local equilibrium. This holds as long the characteristic microscopic scales of the system -e.g. the mean free-path of the particles -are sufficiently small compared to the (dominating) characteristic macroscopic scales of the system. Secondly, the dissipative corrections of the fluid dynamic tensors must be small in comparison to the equilibrium quantities. Both requirements are quantified by introducing the Knudsen numbers, Kn, which are defined as the ratios of the microscopic and macroscopic scales, and the inverse Reynolds numbers Rn −1 , which are defined as the ratios of the magnitude of the dissipative quantities (e.g. the diffusion currents), as well as the corresponding primary hydrodynamic fields (e.g. the local net charge densities). It is often argued that if both measures are small, fluid dynamics is applicable. However, it was recently shown that in some situations the applicability of fluid dynamics extends even up to Kn ∼ 1 [57]. In order to close the set of fluid dynamic equations, one needs to introduce additional equations of motion for the dissipative quantities. Following the approach of transient fluid dynamics in DNMR (Denicol-Niemi-Molnár-Rischke) theory [58], the equations of motion for the bulk viscous pressure, the diffusion currents, and the shear-stress tensor are introduced. The source terms responsible for the generation of dissipation in these equations are expanded in orders of Knudsen numbers and inverse Reynolds numbers under the assumption that they are sufficiently small so that the higher order contributions can be neglected.
For these first investigations where we want to examine the impact of the off-diagonal terms in the diffusion coefficient matrix, we only expand the source term to first order in the Knudsen numbers. The transient equations of motion then read [58]: where we introduced the bulk viscosity ζ, the shear viscosity η and accounted for the diffusion coefficient matrix is the comoving time derivative. The first order source terms correspond with the source terms from Navier-Stokes-Fourier theory [40,41]. We note that the Navier-Stokes terms do not contain any direct cross-couplings between the dissipative fields. In the following, we neglect bulk and shear, Π = π µν = 0, and focus on the diffusion without viscous corrections. Thus, the only dissipative equations of motion we consider in this work are the equations for the net diffusion currents, In order to solve these fluid dynamic equations of motion, we rewrite the set of equations in an appropriate manner ( see Appendix X B) and use the numerical solver SHASTA [59,60]. For the sake of simplicity, we only assume longitudinal dynamics in a hyperbolic (1+1)D-geometry characterized by the proper time, τ ≡ √ t 2 − z 2 , and the spacetime rapidity, η s ≡ arctanh (z/t). We then solve Equations (78), (79), (80) and (81) numerically, and use Eqs. (75) and (76) to infer the LRF quantities. Furthermore, we solve Eq. (77) with Newtons secant algorithm in order to find the velocity. Please note that all equations in Appendix X B are already given without any viscous corrections.

B. Equation of state
In order to close the set of fluid dynamics equations, we need to impose an equation of state, P 0 (T, µ B , µ Q , µ S ). In the non-interacting hadron gas it is straightforward to compute thermodynamic quantities as a function of T and µ q , eq ≡ eq (T, µ B , µ Q , µ S ), n q,eq ≡ n q,eq (T, µ B , µ Q , µ S ).
However, in fluid dynamics the natural variables are energy-and net charge densities, and we need to invert these relations numerically in order to obtain the pressure, temperature and the chemical potentials as a function of and n q . Here we assume the same classical hadronic system as presented in Chapter VI, and thus the single-particle distribution function is of Maxwell-Juettner type (12), and the thermodynamic quantities can be expressed as, for q ∈ {B, Q, S}, P 0,eq ≡ 1 3 We note that the equation of state constructed in this way is consistent with the equilibrium state in the computation of the diffusion matrix.

C. Results
In order to obtain some understanding of the diffusive interplay between the multiple conserved charges and the importance of the diffusion coefficient matrix, we simulate the dynamics of the hadronic system presented in Chapter VI.

a. Case study
For the sake of simplicity, we consider only two conserved charges in the system, the net baryon number, and net strangeness by setting the electric chemical potential to zero, µ Q = 0. We set simple initial conditions and consider four different configurations of the diffusion coefficient matrix of the system: • Case 1 : No diffusion; all the diffusion coefficients are set to zero, κ qq = 0.
• Case 2 : Baryon diffusion only; the only non-vanishing coefficient is κ BB , which is taken from the abovementioned evaluation. This case is usually assumed in other works [46,47]. κ BB is computed with the linear response method as described in the first part of this paper in the relevant range of temperature and chemical potentials. The only restriction implied is that the electric chemical potential vanishes, µ Q = 0.
• Case 3 : Off-diagonal entries neglected; all the off-diagonal entries of the coefficient matrix are artificially set to zero. Note that the only off-diagonal coefficient is κ SB = 0. All the diagonal coefficients are again taken from the above-mentioned calculation.
• Case 4 : Full diffusion matrix; the complete diffusion coefficient matrix of the system is considered.
We assume simple transversally homogeneous initial conditions for a heavy ion collision at small collisional energies, which is entirely in the hadronic phase and suffers large longitudinal gradients in net baryon number. In all of the above-introduced cases, the system is initialized at proper time τ 0 = 2 fm/c, with a homogeneous temperature of 160 MeV and a double-gaussian profile in initial net baryon number density, where η s,0 = 1.0, n B,max = 0.5 fm −3 and R 0 = 0.5. Furthermore, we set the initial net strangeness density to zero everywhere, n S = 0, and as usual, the initial fluid velocity is zero, u µ = 0 3 . From these specifications, the energy density is calculated from the equation of state, which results in a non-homogeneous profile. This implies that besides the diffusion, the dynamics of the system are also determined by gradients in pressure. Therefore, even in the non-diffusive case (Case 1), the baryon number is transported with the flow of the system.

b. Description
We show our results for the evolution of the system for each of the four assumed cases in Fig. 10 for proper times starting at τ = τ 0 = 2 fm/c until τ = 7 fm/c. The evolution of the net baryon number (left side of the plot) and the net strangeness density (right side of the plot) is presented for the four cases: no diffusion (top row), full diffusion matrix (second row), no off-diagonal entries (third row) and baryon diffusion only (bottom row). Various colored and dashed curves are plotted over spacetime rapidity, η s , representing the state at four different proper times: initial state at τ = 2 fm/c (black solid curve), at τ = 3 fm/c (blue dashed line), at τ = 5 fm/c (orange dotted line), and finally at τ = 7 fm/c (red mixed dashed line). We emphasize that two plots in one row belong to the same case. We show the densities multiplied by the Bjorken factor τ /τ 0 in order to account for the longitudinal expansion [61].
First, we note that in the non-diffusive case (Case 1), only a small amount of the baryon number is transported towards the mid-and outwards rapidities from the regions of high baryon densities due to the motion of the fluid (convection generated by pressure gradients). Furthermore, there is no transportation of net strangeness. Accordingly, in Case 2, there is also no extra transported net strangeness and the distribution of net strangeness remains flat at zero. However, in contrast to the non-diffusive case, there is significant diffusive transport of the net baryon number. All diffusive cases (Cases 2 to 4) show a very similar evolution of the net baryon number, but the evolution of the net strangeness is sensitive to the assumed configuration of the diffusion coefficient matrix. Contrary to Case 2, a wave-like profile in net strangeness density builds up over time in Cases 3 and 4, while the total net strangeness is conserved globally. This profile is more pronounced if the off-diagonal entry is neglected. To assess the magnitude of this effect, we compare the net strangeness density to the total number density, n tot . For Case 3, we find ratios up to |n S /n tot | ∼ 6%, and in the consideration of the full diffusion matrix (Case 4), only ratios up to ∼ 3% are reached during the evolution in this example. Aside from the differences in magnitude, there are also differences in the wave-like profile that appears, depending on the assumed case.

c. Interpretation
The reason for this separation of strangeness is that the Navier-Stokes terms of the corresponding diffusion currents introduce a coupling between the charge currents via the diffusion coefficient matrix (see Eq. (17)) and the assumed equation of state.
In Case 3, positive baryon number and positive strangeness is transported to the mid-and outward rapidity region. Due to charge conservation, less baryon number and negative net strangeness stays behind in the regions of the baryon source. From this case we see that even though we did not assume any explicit coupling through κ SB in the fluid dynamic equations, we can still witness a correlation between the conserved charges. The origin of this intrinsic correlation of charges introduced by the equation of state alone is the same as the correlation introduced by the diffusion matrix: the particles carry a multitude of conserved charge types. This in turn results in the fact that chemical potentials are generally dependent on all assumed charge densities and vice versa. Thus, this chemistry of "mixed" charges already encodes charge-correlation into the equation of state, see e.g. [62]. However, in order to achieve physically correct results for the charge-correlation during the dynamic evolution, it is important to ensure that the same chemistry is assumed for the calculation of the diffusion matrix as well.
This is demonstrated with Case 4, where we assumed the full diffusion coefficient matrix. We find a similar picture as in Case 3. However, because κ SB is negative (as shown in Section VI), and the gradients in α q have the same sign, the influences of both gradients on the diffusion currents inhibit or even cancel each other out in this configuration, which leads to a different evolution of the net strangeness in comparison to Case 3. The exact effects of this off-diagonal entry in the diffusion matrix clearly depend on the profiles of temperature and chemical potentials in a complicated manner, since the diffusion coefficients are also a function of these quantities.

d. Summary
In this section, we presented first results which imply that choosing an equation of state and diffusion coefficients in an inconsistent way could make a difference in the evolution of a system that consists of particles carrying a multitude of conserved quantum numbers. However, these investigations were done in a simple manner and results from the evolution were not transferred to particle spectra. The question of whether the influence of the full diffusion coefficient matrix is significant is therefore left for more sophisticated and detailed future works.

IX. CONCLUSION AND OUTLOOK
In the first part of this paper, we introduced the diffusion coefficient matrix in order to account for the fact that, especially in hadronic gases, particles generally carry multiple types of conserved quantum numbers. We find that the mixed chemistry results in a coupling of all diffusion currents that correspond to the conservation of these conserved quantities. In order to describe heavy ion collision, we propose that this coupling must be accounted for in dynamic simulations. We evaluated the complete diffusion coefficient matrix for two examples: a hadron gas and a simple model for a massless QGP, both containing conserved baryon number, strangeness and electric charge. This was done by using a semi-analytical linear response approach in relativistic kinetic theory, and we compared our results for the diagonal coefficients to Ref. [24] in the case of the massless, conformal QGP. We find that the off-diagonal coefficients κ BQ , κ SB and κ SQ , which describe the mixing between the diffusion currents, can reach similar magnitudes to the diagonal coefficients, κ BB , κ QQ and κ SS , which are usually evaluated in other approaches e.g. Ref. [24,63].
Dynamic simulations or other model descriptions of high density heavy ion collisions in experiments like RHIC BES, NICA or FAIR will become increasingly important. We used the evaluated diffusion coefficient matrix and presented a first study of the influence of the matrix in a simple (1+1)D-fluid dynamic simulation of a hadronic system. In addition, signals of strangeness separation and significant baryon diffusion were found and discussed. The results imply that inconsistently choosing the equation of state and the diffusion coefficient matrix of the system results in false dynamics of the conserved charge, which could mislead the physical interpretation. We therefore advise that the mixing between the diffusion currents should not be neglected in simulations of high density heavy ion collisions. However, the relevance of these effects for experimental observables has not yet been investigated. Furthermore, significant effects from e.g. transverse dynamics, shear viscosity and second order contributions to the diffusion currents are expected, but were neglected in this first investigation. This, as well as other aspects remain open for more sophisticated works in the future. Moreover, a comparison of our results to lattice QCD, other transport models or dynamic approaches are also desirable for future research.  Figure 10: Longitudinal fluid dynamic evolution of the net baryon number (left plots) and net strangeness (right plots) multiplied by the Bjorken factor τ /τ 0 of a classical, hadronic system, with 19 assumed particle species (see Appendix X A) for different configurations of its diffusion matrix (see upper left corner of left plots). The state of the net densities is shown at different evolution times (various colored and dashed lines) beginning at the initial proper time τ 0 = 2 fm/c (black, solid curve) and plotted over the spacetime rapidity η s . The system is prepared at initial, homogeneous temperature T 0 = 160 MeV, a double-gaussian profile in net baryon number with maxima at η s = ±1.0 with value n B,max = 0.5 fm −3 and initial vanishing local net strangeness density n S .  π + π − π 0 K + K − K 0K 0 p npn Λ 0Λ0 Σ 0Σ0 Σ +Σ+ Σ −Σ− π + 10 res res 10 10 res 10 res 10 10 res 23.1 23.1 5 5 5 5 5 5 π − 10 res res 10 Fig. 2. We use constant cross sections where no resonance cross section was available from UrQMD [51,52] .

B. Fluid dynamic equations: (1+1)-dimensional longitudinal system in hyperbolic coordinates without bulk and shear viscosity
The transformation law between cartesian and hyperbolic coordinates (t, x, y, z) ↔ (τ, x, y, η s ) reads t = τ cosh(η s ), x = x, y = y, z = τ sinh(η s ), where τ is the proper time and η s is the spacetime rapidity. The metric in hyperbolic coordinates reads g µν = diag 1, −1, −1, −τ 2 and the only non-vanishing Christoffel symbols of second kind are: Γ τ ηη = τ and Γ η τ η = Γ η ητ = 1 τ . The fluid velocity simplifies to u µ = γ η (1, 0, 0, v η ), with the Lorentz factor γ η = 1 √ 1−τ 2 (v η ) 2 , and the four-derivative reading ∂ µ = (∂ τ , 0, 0, ∂ η ). The expansion scalar can then be expressed as We set all initial values of the dissipative fields to zero. Therefore, the non-vanishing fluid dynamic fields read: T τ η = (T τ τ + P 0 )v η (69) T ηη = T τ η v η + P 0 τ 2 (71) N τ q = n q γ η + j τ q (72) Due to the orthogonality of the diffusion current, j η q is the only independent component, and therefore Using the explicit form of the fluid dynamic tensor, we can express local rest frame quantities in terms of the lab frame quantities as Further, the fluid velocity can also be connected to the lab frame quantities, and is evaluated using Newtons secant algorithm because the components of the energy-momentum tensor are dependent on the fluid velocity itself. In order to evaluate the primary fluid dynamic fields -, n q and v η -the lab frame quantities -T τ τ , T τ η , N τ q and j η q -must be calculated by solving the fluid dynamic equations of motion. The explicit form of these read: and the equation of motion for the diffusion current in η s -direction reads: and Further, we propose simple estimates for the relaxation times motivated from Ref. [58] τ q ≡ 12κ qq n tot .