Role of bulk viscosity in deuteron production in ultrarelativistic nuclear collisions

We use a Bayesian-calibrated multistage viscous hydrodynamic model to explore deuteron yield, mean transverse momentum and flow observables in LHC Pb-Pb collisions. We explore theoretical uncertainty in the production of deuterons, including (i) the contribution of thermal deuterons, (ii) models for the subsequent formation of deuterons (hadronic transport vs coalescence) and (iii) the overall sensitivity of the results to the hydrodynamic model -- in particular to bulk viscosity, which is often neglected in studies of deuteron production. Using physical parameters set by a comparison to only light hadron observables, we find good agreement with measurements of the mean transverse momentum $\langle p_T \rangle$ and elliptic flow $v_2$ of deuterons; however, tension is observed with experimental data for the deuteron multiplicity in central collisions. The results are found to be sensitive to each of the mentioned theoretical uncertainties, with a particular sensitivity to bulk viscosity, indicating that the latter is an important ingredient for an accurate treatment of deuteron production.

hadrons, and light nuclei. The recent measurements of deuterons, 3 He, 4 He, and 3 Λ H at the Large Hadron Collider (LHC) in Pb-Pb collisions at √ s NN = 2.76 TeV and 5.02 TeV [1][2][3][4] have renewed interest in the mechanisms of light nucleus production. The improved precision of measurements provides an opportunity to revisit and test current models of deuteron production. Three broad types of models are typically used to calculate the production of deuterons in such collisions. The thermal model assumes that light nuclei reach equilibrium with hadrons, until a point of chemical freeze-out, where hadron and nuclear abundances are frozen, which happens at approximately the same temperature across the entire system [5]. To obtain a momentum distribution for the light nuclei, one can combine the thermal model with a blast wave model at a later kinetic freeze-out [6]. Another approach is the coalescence model, which predicts that the number of produced light nuclei is a convolution of (i) a source function characterizing the nucleons' distribution in phase space, and (ii) the light nucleus' Wigner function. 1 Realistic source functions (space-timemomentum distributions) for the final hadrons can be obtained from hydrodynamics and transport models of heavy-ion collisions. The main assumption of the coalescence model is that light nuclei are formed when hadronic interactions cease. Finally, one can model the production of light nuclei with a transport approach which implements specific production and dissociation reactionsfor example N pn ↔ N d, πpn ↔ πd, or πd ↔ pn [12][13][14][15][16] (N denotes either p or n) -or binding of nucleons in transport by potentials [17,18].
Ultra-relativistic heavy-ion collisions are frequently simulated by multistage approaches, where the evolution of the plasma is modeled by relativistic viscous hydrodynamics, and the subsequent hadronic rescattering at a more dilute stage is described by kinetic theory. Conveniently, in a multistage simulation, one can test all three types of light nucleus production models. One can sample light nuclei from a near-equilibrium distribution at the transition from hydrodynamics to transport, similar to a thermal model at chemical freeze out; one can include the reactions involving deuterons into the transport phase; or one can use the final state nucleons as a source function for coalescence.
In all of these models of light nucleus production, there is a close connection between the distribution of nucleons and the light nucleus bound states that they form. Therefore, it is evident that a successful description of 1 Although the quantum mechanical foundations of the coalescence model have been studied for many years [7][8][9][10], modern implementations can still vary substantially both in methods and results, especially for nuclei heavier than the deuteron. Different types of coalescence models are listed, for example, in Ref. [11]. We will use a model where nucleons from transport, free-streamed to the larger of their last interaction times, are convoluted with a Gaussian deuteron Wigner function in their center of mass frame; see Section II C.
light nuclei relies on a good description of the underlying system evolution. To test the above light nucleus production models, we employ a hybrid hydrodynamic and hadronic transport model that was calibrated to a wide set of hadron observables [19,20]. Importantly, this multistage model includes bulk viscosity, unlike most previous studies of light nuclei. While uncertainties remain in the magnitude of bulk viscosity of QCD, and in particular the modeling of its effect on particlization, it is generically expected to have a larger effect on heavier particles, and therefore should be important for the production of light nuclei. We show in this work that bulk viscosity indeed has a large effect on deuteron production, regardless of the underlying deuteron production model. This work is organized as follows. In Section II, we first briefly describe our multistage hydrodynamic model and deuteron production models used in this work. We then compare the multistage model with measured deuteron observables -yield dN/dy, mean transverse momenta p T , and azimuthal angular anisotropy v 2 -in Pb-Pb at the LHC in Sections III A-III B. We further quantify via Bayesian inference the additional constraints provided by deuteron observables on the initial conditions and bulk viscosity of the plasma in Section III C, and summarize our results in Section IV.

A. Multistage model of heavy ion collisions
A detailed description of the hybrid hydrodynamictransport model used throughout this work can be found in Ref. [20]. Briefly, the T R ENTo model [21] is used as a parametric initialization for the energy density shortly after the impact of the nuclei. This profile is then freestreamed for a short proper time, and is then used as initialization for second-order relativistic viscous hydrodynamics (MUSIC [22][23][24]). We use the shear and bulk viscosity parametrizations described in Refs. [19,20] and an equation of state which matches the trace anomaly of lattice calculations [25] and the hadron resonance gas. On a surface of constant temperature T sw , the Cooper-Frye prescription [26] is employed to switch description ("particlization") from fluid to a kinetic theory of hadrons, which then scatter, decay and form resonances in the SMASH hadronic afterburner [27][28][29].
In the ideal hydrodynamic limit, as in the case of thermal models, the hadrons in the fluid are in local equilibrium. In kinetic theory, it implies that their momentum distribution is given by the Bose-Einstein or Fermi-Dirac distribution. In the more general case where there is dissipation in hydrodynamics, there is an off-equilibrium correction to the hadron distribution, and this viscous correction is model dependent [30][31][32]. In this work, we use the "Grad" model [33] for viscous corrections to the equilibrium distribution function, the one out of several studied in Refs. [19,20] that gave the best agreement with light hadron measurements (listed below).
The model parameters were calibrated by Bayesian parameter estimation against observables for Pb-Pb collisions at √ s NN = 2.76 TeV as well as Au-Au collisions at √ s NN = 0.2 TeV. The Pb-Pb observables include the yields and mean transverse momenta of pions, kaons and protons [34], the charged particle multiplicity and transverse energy [35,36], the charged particle elliptic, triangular and quadrangular flow [37], and the charged particle mean transverse momentum fluctuations [38]. The calibration observables for Au-Au collisions include the yields and mean transverse momenta of pions and kaons [39], and the charged particle elliptic and triangular flow [40,41]. In particular, no deuteron or light nuclei measurements were used in the calibration of these parameters. The parameters used to generate the predictions in the next section are listed in Table I. 2 We explore several different models of deuteron production, including sampling on the switching hypersurface, production via three-body scattering in the hadronic afterburner, and coalescence on the kinetic freezeout surface.

B. Deuteron production with transport
We investigate two models of deuteron production with transport, which are (i): "transport only", which assumes that no deuterons are present at the transition from hydrodynamics to hadronic transport (i.e., at "particlization"); rather, all deuterons are created during the hadronic rescattering phase by reactionsπd ↔ πnp, N d ↔ N np,N d ↔N np, πd ↔ N N and all of their CPT-conjugates, with elastic πd, N d andN d also changing the deuterons' momenta; (ii): "Cooper-Frye + transport" which assumes that deuterons are nearly equilibrated with the hadron resonance gas at the transition from hydrodynamics, so that they are sampled according to near-equilibrium distribution functions and subsequently allowed to be formed and/or destroyed in the hadronic rescattering phase. Notice that the yield of deuterons at the Cooper-Frye sampling in the model (ii) is closely related to the results of a thermal model (discussed in the introduction) with freeze-out temperature T sw and volume V = dσ µ u µ , where dσ µ are the normal 4-vectors to the hypersurface, and u µ are the collective velocities of the hypersurface elements. These quantities enter the Cooper-Frye formula [26] with f (P ν u ν /T ) for deuterons being the Bose-Einstein distribution, δf being a correction due to viscosity, and g their degeneracy factor. Typically, thermal models assume local equilibrium, δf = 0.
In this work we use the implementation of these 3 → 2 reactions via an intermediate fictitious d resonance [14]: Recently a possibility to simulate 3 → 2 reactions directly, without d , via stochastic rates was implemented in SMASH [42], but it is not employed in this work. In Ref. [42] it was shown that the main difference between direct 3 → 2 reactions and reactions with the intermediate d resonance is that 3 → 2 reactions bring the deuteron yield to equilibrium faster.

C. Deuteron production with coalescence
Here we employ a Wigner-function coalescence model [9] that has no free parameters. The implementation of coalescence is based on the equation where the factor 3/8 originates from spin and isospin averaging. The deuteron Wigner function is where the deuteron wavefunction ϕ d is assumed to be a Gaussian with d = 3.2 fm [43]. The function W is the probability to find a pair of nucleons at positions r d ±r/2 with momenta p d /2 ± q. In our model the distribution W is represented by the nucleon pairs in the hadronic afterburner. The integrals in Eq.
(2) are computed as a sum over all nucleon pairs in the simulation. Every pair is transported in time to the latest of their last collision times. At this moment r 1 , r 2 are coordinates and p 1 and p 2 are momenta of the nucleons. Then and the pair contributes to deuteron spectra with the weight 3 8 D(r, q).

III. RESULTS
A. Comparing transport, coalescence, and thermal-like deuteron production To investigate the different mechanisms of deuteron production, we compare the three specific scenarios presented in Sections II B and II C within our multistage model of heavy ion collisions: (i) Transport only: deuterons not present at particlization, followed by hadronic transport with deuteron reactions; (ii) Cooper-Frye + Transport: deuterons present at particlization and allowed to react in the hadronic transport; (iii) Coalescence only: deuterons not present at particlization, hadronic transport without deuteron reactions, followed by coalescence at kinetic freezeout.
While these do not represent all possible variations, they provide sufficient information for understanding the effects of each mechanism.
In Fig. 1 one can see that these three scenarios result in rather similar yield and mean transverse momentum of deuterons. Let us concentrate on the comparison of the transport scenarios first. As in previous studies [14,15], deuteron yields and average transverse momenta are similar whether or not deuterons are sampled at particlization. To understand this effect, note that (i) the reaction πd → πpn has a large cross section, (ii) the implementation of the reverse reaction obeys the detailed balance principle, and (iii) a large number of pions is produced in heavy ion collisions. Given these conditions, as long as the πd ↔ πpn reaction rate is larger than the expansion rate of the plasma, deuterons rapidly approach relative equilibrium with protons. When the deuteron is not sampled at particlization, the deuteron yield approaches but does not completely reach equilibrium [14,15], because the fireball's expansion freezes out the πd ↔ πpn reactions; this is the reason for the slightly smaller number of deuterons found in the "Transport only" scenario. Note that, to avoid confusion, one should think about deuterons being in equilibrium in a statistical sense, averaging over a large number of heavy-ion collisions, and not in a single event -notice that on average only around 0.1 deuterons are produced per event per unit of rapidity, as shown in Fig. 1.
In Refs. [14,15] both scenarios (i) and (ii) were found to be compatible with the experimental data, although sampling deuterons at particlization was slightly preferred. In this work, the situation is different: the experimentally measured deuteron yield and transverse momentum are better described by omitting the deuteron from particlization (see Fig. 1). While there are multiple differences between the multistage models used in the present work and Refs. [14,15], the likely reason for this difference is bulk viscosity, which was not included in previous works. We discuss this effect in detail in Section III C. We note that the (mostly systematic) uncertainties are still significant on the deuteron measurements [44] shown in Fig. 1. Reduction of these systematic uncertainties would increase the discrimination power of deuterons even further.
The creation of deuterons by reactions has a certain similarity with coalescence. With reactions, most of finalstate deuterons are produced rather late; earlier produced deuterons tend to get destroyed by subsequent collisions. Moreover, by the kinematics of the πpn → πd reaction, the incoming proton and neutron have to be close in phase space, as coalescence assumes. In Fig. 1 one can indeed see that the coalescence model provides a very similar deuteron yield and transverse momentum as the models with deuteron-producing reactions. Therefore, it is the underlying proton phase space distribution (which is the same in all three models) that influences the deuteron observables here, while the deuteron production mechanism is less important. 3 As a consequence, the model parameters that influence proton production will also influence deuteron production. This means that by combining proton and deuteron observables one could potentially constrain these parameters better than by using only proton observables.
Based on the agreement of model predictions with measured proton data, one might expect a better agreement for deuteron yields. There are several reasons why the deuteron yield is not as well described as expected. First, although the model is tuned to describe integrated proton yield and p T precisely, the proton p T -spectrum in fact exhibits deviations from experiment [20, Fig. 17]. Second, previous studies suggesting that a good agreement of proton data implies good agreement with deuteron data [14,15] did not take into account bulk viscosity, unlike in the present work. Importantly, the bulk viscous corrections to the proton p T -spectrum are substantial (see Fig. 4 below). In this section we provide a prediction for the multiplicity and mean p T of deuterons for Pb-Pb √ s NN = 5.02 TeV collisions using the "Transport only" approach. We further compare our calculation for p T -differential v 2 with ALICE measurements for deuterons in both Pb-Pb √ s NN = 2.76 TeV and √ s NN = 5.02 TeV.
The initial conditions, transport coefficients and other parameters of the multistage model given in Table I [19,20]. We assume that all model parameters remain the same except for two T R ENTo initial condition parameters that are expected to be center-of-mass energy dependent: (i) the nucleonnucleon inelastic cross section and (ii) the overall normalization of the initial energy density. For the inelastic nucleon-nucleon cross section at √ s NN = 5.02 TeV, we used 70 mb. The normalization of the initial energy density is typically a parameter tuned to heavy ion measurements, mostly the hadronic multiplicities. In our approach, instead of re-tuning it to measurements, we simply estimated it from a previous Bayesian inference that also used T R ENTo initial condition [45]. In that work, 3 This may be different in small systems where the span of the deuteron wave function is comparable to the system size. prediction for the deuteron multiplicity and mean p T at √ s NN = 5.02 TeV is shown on the same figure. As is the case for light hadrons, it is natural to expect our predictions for deuterons at 5.02 TeV to have very similar agreement as for the 2.76 GeV results (see the "Transport only" curve in Fig. 1) -that is, generally good agreement except for an overestimated yield in central collisions.
The p T -differential v 2 of deuterons in Pb-Pb collisions at √ s NN = 2.76 TeV and 5.02 TeV is described very well for different collision centralities, as shown in Fig. 3. The √ s NN = 5.02 TeV v 2 result was shown as a prediction in the ALICE publication [4]. We evaluate the p Tdifferential deuteron momentum anisotropy v 2 {2} using the Q-cumulant method [49].

C. Sensitivity to medium properties
In this section, we explore the sensitivity of the deuteron yield and mean transverse momentum to prop- Deuterons are produced at particlization and allowed to dynamically form and be destroyed, corresponding to the "Cooper-Frye+Transport" scenario discussed in Section III A. Note that the effect of the viscous correction on pions, protons and other hadrons propagate to deuterons through the transport phase. ALICE measurements [34,44] are plotted as black triangles. erties of the hydrodynamic medium.
As discussed in Section II, deviations of the plasma from local thermal equilibrium lead to modifications in the corresponding hadronic momentum distribution from Bose-Einstein or Fermi-Dirac. This "viscous correction" to the equilibrium distribution function is related to the magnitude of the bulk pressure. Its dependence on the hadron mass depends on the model used to calculate the viscous corrections. For the Grad model used in this work, this viscous correction increases with the hadron mass. 5 As a result, heavy particles such as protons, neutrons, and especially deuterons might be expected to have a higher sensitivity to bulk viscosity, compared to the majority of produced hadrons. Despite this, there has been parameter min. max (ζ/s)max 0.03 0.15 k 0.3 2 w [fm] 0.5 1. 5   TABLE II. Range of model parameters used to produce the model predictions in this section. Note that the highest bulk viscosity to entropy ζ/s attained in the fluid before switching (at temperature Tsw = 0.136 GeV) ranges from 0.029 (ζ/s)(Tsw) 0.143; this is slightly smaller than the nominal maximum value (ζ/s)max defined at temperature T ζ,c = 0.12 GeV (see Table I for the value of the model parameters).
FIG. 5. The allowed range of specific bulk viscosity permitted by the chosen Bayesian prior. Only the magnitude was varied, while the shape parameters were held fixed by the values in Table I. no systematic study of the role of bulk viscosity in the production of deuterons until now.
In Fig. 4, we investigate the relative importance of bulk viscosity by comparing the yield of each particle (solid line) to the case where the bulk viscous correction at particlization has been set by hand to zero for all particles (dotted line). One can see that the importance of the viscous correction indeed increases significantly with mass, and that the yield of deuterons is affected much more than that of lighter hadrons; generically the identified hadron multiplicity gets enhanced whereas the mean p T gets reduced by the bulk viscous correction. Note that, while the distribution of deuterons at particlization does not have a strong effect on the final deuteron distribution (see discussion in Section III A), any change in the distribution of protons and pions subsequently feeds down to the deuteron.
To better understand the role of bulk viscosity on deuteron production, we proceed with a Bayesian analysis with three parameters of interest. Two parameters, k and w enter the initial conditions via the T R ENTo model. The parameter k ≡ 1/σ 2 k controls the magnitude of the fluctuations of the deposited energy in each nucleonnucleon collision. The parameter w, referred to as the "nucleon width", controls the transverse radius of the nu-cleons in T R ENTo, and defines the transverse size of deposited energy given a nucleon-nucleon collision. Both of these parameters largely control the homogeneity of the initially-deposited energy density. Finally, we vary the magnitude of the specific bulk viscosity at its peak value, (ζ/s) max . The temperature dependence of bulk viscosity in this work is assumed to have a skewed-Cauchy form as in Ref. [20, Fig. 1]. The priors for each parameter are assumed to be uniform within the ranges listed in Table II. The range of temperature-dependent specific bulk viscosities spanned by the prior for (ζ/s) max is shown in Fig. 5. For deuteron production, we use the "Coalescence only" model described in Sections II and III A.
Forty-five design points were sampled from the prior using a Latin hypercube design [50]. Combined with the fixed values of all remaining model parameters, the model's prediction (given by the model's prior predictive distribution) for these 45 samples of k, w and (ζ/s) max are shown in Fig. 6 for the distributions of pion, proton and deuteron yields and mean p T in Pb-Pb √ s NN = 2.76 TeV collisions. We perform the Bayesian inference along the lines of previous works [19,20,45,[51][52][53][54][55][56]: For each observable, an emulator is used to interpolate the model's results between the parameter point samples. The major difference with previous work is that we use a more sophisticated Gaussian process emulator. 6 To illustrate the sensitivity of deuteron observables to the magnitude of the bulk viscosity, we fix the T R ENTo fluctuation k and nucleon width w to the midpoints of their prior, and plot the emulated response of the deuteron yield to changes in the specific bulk viscosity (ζ/s) max . This is shown in Fig. 7, and we see that the deuteron yield indeed shows strong sensitivity to the magnitude of bulk viscosity.
In Fig. 8 we show the effect of adding deuteron observables on constraining the nucleon width w, the fluctuation parameter k and the maximum of bulk viscosity (ζ/s) max . We see that the deuteron's dependence on bulk viscosity modifies the value of (ζ/s)(T ) that is in best agreement with measurements. Moreover, Fig. 8 shows how this change in bulk viscosity correlates with changes in the preferred value of the initial condition parameters. For example, we see that including deuteron observables would favor a slightly smaller (ζ/s) max and a larger w; the former can be understood from the fact that FIG. 6. Prior predictive distributions for the particle yields (top) and mean pT (bottom), given by model predictions at the 45 points sampled uniformly in the prior volume. These points form the Gaussian process design. ALICE data from Refs [34,44].
simulations with parameters calibrated without deuteron observables overestimate the deuteron multiplicity (see Fig. 1a) whereas decreasing (ζ/s) max helps to reduce the tension (see Fig. 7), and the latter can be understood by its anticorrelation with w. We do not intend these results to represent accurate calibrations of the three model parameters considered here, in part because we only explored a small subspace of the much larger parameter space considered in Ref. [20], and in part because we only calibrated the three model parameters against a small set of observables. Nevertheless, Fig. 8 illustrates the sensitivity of deuterons to bulk viscosity and their potential for improving constraints on the properties of quark-gluon plasma.

IV. SUMMARY
We explored deuteron production in ultra-relativistic Pb+Pb collisions at √ s NN of 2.76 TeV and 5.02 TeV, two collision systems where recent deuteron measurements are available. For this purpose we employed a multistage approach (hydrodynamics + hadronic afterburner) tuned to reproduce the yields, mean transverse momenta, and flow of pions, kaons, and protons in Au-Au collisions at √ s NN = 0.2 TeV and Pb-Pb collisions at √ s NN = 2.76 TeV; no deuteron observables were used for tuning. Three different models of deuteron production were tested -"Transport only", "Cooper-Frye + Transport", "Coalescence only" (Section III A). Overall, all three models produce rather similar results. At 2.76 TeV they reproduce the centrality dependence of deuteron p T and v 2 (p T ) within error bars, while the deuteron yields are overestimated in central collisions but reproduced well in more peripheral ones. It is possible that a more realistic deuteron wave function might affect the centrality dependence and lead to improvement in central collisions but we have not checked this.
Our predictions for deuteron flow at 5.02 TeV had been confronted with experimental data in Ref. [4], with good overall agreement. This was expected, since the deuteron flow is not much different at 2.76 and 5.02 TeV and the model reproduced the deuteron v 2 (p T ) at 2.76 TeV very well. The data for deuteron yields at 5.02 TeV are not yet published by ALICE -although analogously to 2.76 TeV, it would not be surprising that our prediction overestimates the yield in central collisions, reproduces it in peripheral collisions, and reproduces the p T precisely. Despite the above tension with 2.76 TeV measurements in central collisions, we expect our prediction for the ratio of measured yields at 5.02 TeV and 2.76 TeV to be more robust.
The main conclusion of this study is that deuteron observables are particularly sensitive to bulk viscosity. We have seen in Fig. 4 that when bulk viscous corrections change the proton and neutron yield by 20-25%, the deuteron yield can be changed by as much as 50%. While the quantitative values for the bulk viscous corrections quoted above are quite large -and might be even pushing the multistage model to its limits [58][59][60][61] -the stronger dependence on bulk viscosity of deuterons compared to protons should be robust.
The fact that deuterons are sensitive to the bulk viscous corrections has an interesting implication: proton femtoscopic radii should also be sensitive to the bulk viscosity. Indeed, a relation between proton femtoscopic radii and coalescence has been explicitly demonstrated recently [10].
The overall dependence of light nuclear observables on bulk viscosity could be used to improve constraints on this transport coefficient, as discussed in Section III C. We have provided a preliminary constraint in Fig. 8; a more robust constraint will require a better understanding of the bulk viscosity in heavy ion collisions, in particular viscous corrections at the transition between hy-drodynamics and transport.