Dynamical mechanisms for deuteron production at mid-rapidity in relativistic heavy-ion collisions from SIS to RHIC energies

The Parton-Hadron-Quantum Molecular Dynamics (PHQMD) has been recently conceived as a new type of microscopic transport approach which combines the characteristics of baryon propagation from the Quantum Molecular Dynamics (QMD) model Aichelin (1991); Aichelin et al. (1987, 1988); Hartnack et al. (1998) and the dynamical properties and interactions in and out of equilibrium of hadronic and partonic degrees of freedom of the Parton-Hadron-String-Dynamics (PHSD) approach Cassing and Bratkovskaya (2008); Cassing (2009); Cassing and Bratkovskaya (2009); Bratkovskaya et al. (2011); Linnyk et al. (2016).

In this section we provide a short summary of these two building blocks. For more details of the PHQMD model we refer to Ref. Aichelin et al. (2020).

I. In QMD the baryons are described by single-particle wave functions of Gaussian form with a time independent width. The Wigner density of each particle is obtained by a Fourier transformation of the density matrix. Then, the n-body Wigner density is constructed by the direct product of the single-particle Wigner densities and its propagation is determined by a generalized Ritz variational principle Feldmeier (1990). Contrary to mean-field approaches, where the n-body phase-space correlations are integrated out and the dynamics is reduced to a single-particle propagation in a mean-field potential, in QMD these correlations are preserved and the fluctuations not suppressed. This allows to investigate the dynamical formation of clusters, which are correlations between nucleons.

In PHQMD a baryon $i$ is represented by the single-particle Wigner density, which is given by

the Gaussian width $L$ is taken as $L=8.66$ fm².

The QMD equations of motion (EoMs) for the centroids $(\mathbf{r}_{i0},\,\mathbf{p}_{i0})$ are similar to those of the Hamilton equations for a classical particle Aichelin (1991)

where the difference lies in the fact that on the right hand side of both equations the expectation value of the quantum Hamiltonian of the many-body system appears. We note in passing that for a non-Gaussian form of the wave function the time evolution equations are different. The Hamiltonian is the sum of the kinetic energies of the particles and of their (density dependent) two-body interaction. The expectation value can be written as

The potential interaction consists of two parts: the Coulomb interaction and a local density dependent Skyrme potential $V_{i,j}=V_{Coul}+V_{Skyrme}$. The expectation value of the Coulomb interaction can be calculated analytically. The expectation value of the Skyrme part contains terms $\propto\rho^{2}$ and $\propto\rho^{\gamma}$, where $\rho$ is the local density. Their weights, as well as the exponent $\gamma$, are tuned to the Equation of State (EoS) of infinite nuclear matter $E(T=0,\rho/\rho_{0}=1)=-16$ MeV, where $\rho_{0}=0.16$ fm^-3 is the saturation density at zero temperature. This fixes two of the three parameters. In PHQMD two parameter sets have been introduced, which yield a “soft” and a “hard” EoS, respectively. For details on the realization of the QMD dynamics and the impact of different EoS on bulk and cluster observables we refer to Ref. Aichelin et al. (2020). For bulk and strangeness particle production in PHQMD with a “hard” and a “soft” EoS at low energy HICs and the comparison with other transport models see also Ref. Reichert et al. (2022). In our present study of the deuteron production mechanisms, we employ the PHQMD in its “hard” EoS setup like in Ref. Gläßel et al. (2022).

Finally, we want to stress two aspects: i) in this study, as in the previous PHQMD publications, we employ a “static”, i.e. momentum independent, Skyrme potential. A form of the Skyrme interaction which contains also a momentum dependent part will be reserved for future work. ii) The PHQMD approach aims to provide a dynamical description of cluster formation in HICs at low-energy, as well as at relativistic energies. The discussed QMD model uses non-relativistic quantum wave functions. The relativistic formulation of QMD dynamics as a n-body theory has been developed in Refs. Sorge et al. (1989); Marty and Aichelin (2013) by introducing extra constraints in order to reduce the 8N-dimensional phase-space to the (6N + 1)-dimensional phase-space in which particle trajectories can be defined. However, this method is computationally extremely expensive, requiring the inversion of a matrix of size $N\times N$ in each time step and thus, it is presently not applicable for high statistics simulations.
Therefore, in PHQMD the original framework of QMD is extended to relativistic energies by introducing a modified single-particle Wigner density for each nucleon $i$:

which accounts for the Lorentz contraction of the nucleus in the beam $z$-direction in coordinate and momentum space by including $\gamma_{cm}=1/\sqrt{1-v_{cm}^{2}}$, where $v_{cm}$ is the velocity of projectile and target in the computational frame, which is the center-of-mass system of the heavy-ion collision. Accordingly, the interaction density modifies as

We refer again to Ref. Aichelin et al. (2020) and Ref. Aichelin (1991) for a more detailed discussion and the explicit formulas. We note that in PHQMD the nuclei are initialized in their rest frame with the Gaussian distributions Eq. (1). The Lorentz squeezing of nuclei by the gamma factor $\gamma_{cm}$ is done after the initialization of the nuclei in their rest frame, so it does not affect the initialization.

II. As in the PHSD (Parton-Hadron-String Dynamics) approach Cassing and Bratkovskaya (2008); Cassing (2009); Cassing and Bratkovskaya (2009); Bratkovskaya et al. (2011); Linnyk et al. (2016), in PHQMD the strongly interacting medium is described by off-shell hadrons and off-shell massive quasi-particles representing the deconfined quarks and gluons of the QGP phase, which is created if the local energy density is larger than a critical value of $\epsilon_{c}\approx 0.5\!$ GeV/fm³. The propagation of these off-shell degrees of freedom, including their spectral functions, is based on the Kadanoff-Baym transport theory Kadanoff and Baym (1962) in first-order gradient expansion from which the Cassing-Juchem generalized off-shell transport equations Cassing and Juchem (2000a, b); Cassing (2009) in test-particle representation are derived (see Ref. Cassing (2021)). The hadronic part is taken from the early development of the Hadron-String-Dynamics (HSD) approach (see Ref. Cassing and Bratkovskaya (1999) for a detailed description of the baryon, meson and resonance species implemented).

The elementary baryon-baryon ($BB$), meson-baryon ($mb$) and meson-meson ($mm$) reactions for multi-particle production are realized according to the Lund string model Nilsson-Almqvist and Stenlund (1987) via the event generators FRITIOF 7.02 Nilsson-Almqvist and Stenlund (1987); Andersson et al. (1993) and PYTHIA 7.4 Sjostrand et al. (2006). Both generators are “tuned” for a better description of the experimental data for elementary $pp$ collisions at intermediate energies Kireyeu et al. (2020).

The partonic part, which describes the QGP phase, follows the description of the so-called Dynamical Quasi-Particle Model (DQPM) Cassing (2007a, b); Peshier and Cassing (2005). In the DQPM quarks and gluons are represented by massive, strongly interacting quasi-particles. They have spectral functions whose pole positions and widths are defined by the real and imaginary terms of parton self-energies. The parton masses and widths are functions of the temperature $T$ (and in the most recent extension Soloveva et al. (2020) also of the baryon-chemical potential $\mu_{B}$) through an effective coupling constant, which is fixed by fitting lQCD results from Refs. Aoki et al. (2009); Borsanyi et al. (2014); Bazavov et al. (2014); Cheng et al. (2008); Borsanyi et al. (2015). These DQPM partons are evolved with their self-energies according to the same off-shell transport equations and scatter by microscopically computed cross sections.

We recall that in PHQMD only the propagation of mesons and partons relies on the PHSD approach, while the baryons evolve according to the QMD dynamics. However, it is always possible to run PHQMD in the “(P)HSD-mode” by switching the baryon propagation back to the mean-field dynamics of HSD. Again we refer to Ref. Aichelin et al. (2020) for a description and detailed studies.

As already stated, in PHQMD the collision integral, which encodes the main scattering/dissipative processes of hadrons and partons, is adopted from the PHSD model. The main hadronic reactions have been implemented for many observables, like strangeness, dileptons, photons, heavy quarks, etc. (cf. examples in the reviews Linnyk et al. (2016); Bleicher and Bratkovskaya (2022)). It contains also in-medium effects, such as a dynamical modification of vector meson spectral functions by collisional broadening Bratkovskaya and Cassing (2008), and the modification of strange degrees of freedom in line with many-body G-matrix calculations Cassing et al. (2003); Song et al. (2021), as well as chiral symmetry restoration via the Schwinger mechanism for the string decay Cassing et al. (2016); Palmese et al. (2016) in a dense medium. The important and pioneering development in the PHSD is related to the formulation and development of the theoretical formalism in order to realize the detailed balance for $m\leftrightarrow n$ reactions based on covariant rates Cassing (2002). This formalism has been implemented in PHSD in Ref. Cassing (2002) for the description of baryon-antibaryon $B\bar{B}$ annihilation of $B=p,\Lambda$ and the inverse reaction of multi-meson fusion to $B+\bar{B}$ pairs; an extension of this first study accounting for all baryon-antibaryon combinations in PHSD has been presented in Ref. Seifert and Cassing (2018). We also mention that the implementation of detailed balance on the level of $2\leftrightarrow 3$ reactions is realized for the main channels of strangeness production/absorption by baryons ($B=N,\Delta,Y$) and pions Song et al. (2021).

One of the main goals of this work is to extend this formalism to the $2\leftrightarrow 2$ and $3\leftrightarrow 2$ processes, which are relevant for the production and disintegration of deuterons. Therefore, we dedicate a separate section to the detailed description of this formalism and its application to deuteron reactions.

As described, the collision processes involving the formation and the breakup of a deuteron are implemented in PHQMD by means of a covariant rate formalism which has been firstly developed in Ref. Cassing (2002) and applied within the PHSD microscopic approach in order to study the baryon-antibaryon production by multi-meson fusion Seifert and Cassing (2018); Seifert (2018). Following the steps of Ref. Cassing (2002), we start by writing the covariant Boltzmann transport equation for the single phase-space distribution function of an on-shell hadron $f_{i}(x,p)$

using the notation $x=(t,\vec{x})$ and $p=(E,\vec{p})$ with the on-shell condition $p^{2}=m_{i}^{2}$ ($m_{i}$ is the rest mass). The left hand side of Eq. (7) contains only the free streaming term because we neglect for simplicity any mean-field interaction. On the right hand side the collision integral $I_{coll}^{i}$ encodes all the multi-particle transitions which involve the hadron i either in the initial or in the final state. Hence, $I_{coll}^{i}$ can be written as a sum over all scattering processes with increasing number of participant particles

Each collision term $I^{i}_{coll}[n\leftrightarrow m]$ accounts for a particular forward scattering process ($\rightarrow$) with $n$-particles in the initial state and $m$-particles in the final state, as well as for the corresponding backward reaction ($\leftarrow$). The forward and backward reactions can be grouped together and in collision theory one usually distinguishes a gain and loss term. Therefore, the on-shell collision term $I^{i}_{coll}[n\leftrightarrow m]$ becomes

In Eq. (III) there are $n+m-1$ integrals over the initial $\vec{p}_{j}$ and final $\vec{p}_{k}$ momenta of all particles, excluding the tagged hadron i (the deuteron) whose momentum is denoted as $p_{1}$.

The quantity $W_{n,m}(p_{1},p_{j};i,\nu\mid p_{k};\lambda)$ is called transition amplitude and is related to the square of the scattering matrix for the transition $p_{1}+\sum_{j=2}^{n}p_{j}\rightarrow\sum_{k=n+1}^{n+m}p_{k}$ where $\nu$ and $\lambda$ denote a particular set of allowed discrete quantum numbers for the particles (except the hadron i) in the initial and final states. The $\delta$ function guarantees the energy and momentum conservation in each individual reaction. Finally, the single-particle distribution functions of the hadrons appear, in particular the functions $f_{1}$ and $f_{j}$ for the forward/loss term $n\rightarrow m$ and the function $f_{k}$ for the backward/gain term $n\leftarrow m$. We assign the arbitrary $\pm$ sign to distinguish between the gain and the loss term such that in our study the reaction which leads to the production of a deuteron (playing the role of the tagged hadron i) is associated to the backward/gain term, while the inverse reaction where the deuteron is destroyed, is associated to the forward/loss term. This choice agrees with the original formulation in Ref. Cassing (2002). In the collision integral of Eq. (III) we have also neglected the quantum statistical corrections, i.e. the Pauli-blocking or Bose-enhancement factors, which multiply the product of the phase-space distribution functions, as well as any anti-symmetrization procedure in the transition amplitude for the fermions involved in the reactions. Only the statistical factor $1/N_{id}!$, which counts the number of identical particles (either bosons or fermions) in the initial and final states, survives. The expression Eq. (III) can be straightforwardly generalized to the off-shell case by including an additional integration over the energy of each single hadron weighted by the associated spectral function (cf. Ref. Cassing (2002)). In PHQMD/PHSD such an off-shell version of the collision integral is adopted for the dynamical interaction of other baryons and mesons which are propagated with self-consistent off-shell transport equations.

The covariant collision number for the $n(i)\rightarrow m$ scattering process is the number of forward reaction events in the covariant 4-volume $d^{4}x$ = $dVdt$ and, therefore, the covariant collision rate is obtained by dividing the loss term in Eq.(III) by the energy $E_{1}$, followed by the integration over the momentum $d^{3}\vec{p}_{1}/(2\pi)^{3}$ and and a summation over the quantum numbers $\tau(i)$ of the tagged hadron i in the initial state of the transition,

Similarly, the covariant collision rate of the backward reaction $n(i)\leftarrow m$ is obtained from the gain term of Eq. (III)

The transition amplitude $W_{n,m}$ in Eq. (III) and Eq. (III) is the same because of the equivalence of the scattering matrix under the detailed balance condition for forward and backward processes. Therefore, both expressions can be analytically or numerically solved knowing the dependence of the transition amplitude $W_{n,m}$ on the kinematic variables. This has been suggested in Ref. Cassing (2002) where, in particular, it has been shown that the collision probabilities of forward and backward reactions can be written in terms of the corresponding many-body phase-space integrals (cf. Appendix C) if one assumes that the transition amplitude $W_{n,m}$ is a function only of the invariant energy of the collision $\sqrt{s}=(\sum_{j=1}^{n}p_{j})^{2}=(\sum_{k=n+1}^{n+m}p_{k})^{2}$. We apply this procedure for deuteron reactions.

The goal of this work is to implement in the PHQMD transport approach the following deuteron reactions: i) the elastic $d\pi\rightarrow d\pi$ and $dN\rightarrow dN$ reactions, as well as $2\leftrightarrow 2$ inelastic $d\pi\leftrightarrow NN$ reactions; ii) the $2\leftrightarrow 3$ inelastic reactions $d\pi\leftrightarrow NN\pi$ and $dN\leftrightarrow NNN$ with all pion species $\pi=\pi^{+},\pi^{0},\pi^{-}$ and $N=p,n$.

Employing the covariant expressions in Eq. (III) and (III) with $n=m=2$ and $i=d$, the collision rate for the elastic and inelastic $d\pi\leftrightarrow NN$ reactions can be written as follows,

where $\sigma^{2,2}_{tot}$ is the total cross section for a two-to-two scattering process which is defined from the $W_{2,2}$ transition amplitude by the well known definition

with the flux factor $I_{flux}$ related the (relativistic) relative velocity $v_{rel}$ of the incident on-shell particles of masses $m_{1}$ and $m_{2}$

In Eq. (III) the sum is performed over the quantum numbers involved in the reaction and it includes also the statistical factor $1/N_{id}!$ which we absorb in the cross section.

The PHQMD/PHSD collision integral for the deuteron reactions is solved numerically by dividing the coordinate space in a grid of cells of volume $\Delta V_{cell}$ and sampling the on-shell single-particle distribution function $f(x,p)$ at each time step $\Delta t$ by means of the test-particle ansatz Cassing (2009)

where $\vec{r}_{j}(t)$ and $\vec{p}_{j}(t)$ are, respectively, the position and the momentum of the particle $j$ at time $t$. By inserting Eq. (15) in Eq. (12) we obtain the collision probability for the $2\leftrightarrow 2$ reactions in the unit volume $\Delta V_{cell}$ and unit time $\Delta t$

which depends on $\sqrt{s}$ through $v_{rel}$ and $\sigma^{2,2}_{tot}$. Employing sufficiently small values of $\Delta V_{cell}$ and $\Delta t$ we solve numerically the $2\leftrightarrow 2$ collisions for the deuterons by the stochastic method, i.e. by calculating the invariant energy $\sqrt{s}$ of each possible reaction and then the associated probability $P_{2,2}$ which is confronted with a random number between 0 and 1. To calculate $P_{2,2}$ for the inelastic $d\pi\leftrightarrow NN$ process we use parametrized expressions for the total cross section $\sigma^{2,2}_{tot}$ which are reported in the Appendix A.

Now we describe the inelastic reactions $d\pi\leftrightarrow NN\pi$ and $dN\leftrightarrow NNN$ and, in particular, how the backward reaction $2\leftarrow 3$ can be fully implemented within the covariant rate formalism adopted in our PHQMD approach. On the one hand, this is physically motivated by the fact that these are the dominant reactions for the production of deuterons in HICs due to their large cross sections, $\sigma_{tot}\simeq 200\,mb$, compared to the sub-dominant channel $NN\rightarrow d\pi$ with $\sigma_{tot}\simeq 10\,mb$. On the other hand, it provides an effective method to describe reactions with more than 2 particles in the entrance channel, which cannot be formulated in terms of cross sections as in Eq. (III).

For the forward reaction, the breakup of deuterons by an incident $N$ or $\pi$, the definition of the covariant collision rate follows straightforwardly and is given by Eq. (III) with $n=2$, $m=3$ and $i=d$,

where $\sigma^{2,3}_{tot}$ is the total inelastic cross section for either the $d\pi\rightarrow NN\pi$ or the $dN\rightarrow NNN$ scattering process, which is defined similarly to Eq. (III) with an extra integration over the momentum of the third particle in the final state. The sum over the internal quantum numbers appearing in Eq. (17) is also absorbed in $\sigma^{2,3}_{tot}$. In Appendix A we provide the parametrized expressions of such inelastic cross sections as a function of $\sqrt{s}$ and we describe in detail how they are obtained from the experimental measurements of the total inclusive cross section for $d\pi$ and $dN$ collisions. Employing again the test-particle ansatz in Eq. (17), we derive the collision probability for the forward reaction

which is a function of $\sqrt{s}$, and we sample stochastically the collisions in the unit volume $\Delta V_{cell}$ and the unit time $\Delta t$ for each PHQMD/PHSD parallel ensemble event. When a collision occurs, we construct the final state of three particles in the center-of-mass of the incident pair by means of standard kinematic routines Kajantie and Karimaeki (1971); Byckling and Kajantie (1971).

The covariant rate for the backward $NN\pi\rightarrow d\pi$ and $NNN\rightarrow dN$ reactions follows from Eq. (III), but we cannot write it in terms of a cross section. However, what is important for us is to obtain a collision probability $P_{3,2}$ in order to apply the stochastic method. With the assumption $W_{3,2}=W_{2,3}=W(\sqrt{s})$ Cassing (2002) the transition amplitude can be moved outside the integration over the momenta of the two particles in the final state. As a result, these integrations can be combined with the $\delta$ function into the two-body phase-space $R_{2}(\sqrt{s},m_{1},m_{2})$ (cf. Appendix C), so that we can write as intermediate step,

with the sum running over the quantum numbers and taking into account also the statistical factor for identical particles. Next, we introduce the $\sigma^{2,3}_{tot}$ of the forward process by inverting its definition from the transition amplitude and use again the condition $W(\sqrt{s})$ to isolate the three-body phase-space $R_{3}(\sqrt{s},m_{3},m_{4},m_{5})$ of the initial particles. Hence,

where $F_{spin}$ and $F_{iso}$ denote the factors coming from the sum over the spin and isospin quantum numbers in the transition matrix. For the spin contribution we have

where in $NN\pi(N)\rightarrow d\pi(N)$ reactions the particles are ordered as follows,

For the isospin part a separate discussion is given at the end of the section.

Combining Eq. (III) and Eq. (20) with the test-particle ansatz, we finally obtain the collision probability for the backward $3\rightarrow 2$ process in the unit volume $\Delta V_{cell}$ and the unit time $\Delta t$,

where in the second line we employ the collision probability for the forward $2\rightarrow 3$ reaction from Eq. (III). We notice that on the right hand side of Eq. (20) and  (III) the energy of the produced particles $E_{1}^{f}$ and $E_{2}^{f}$ appear. That means that in our numerical implementation we have to sample the possible kinematics of the final state before the collisions take place. If the collision occurs, we reconstruct the kinematics of the emitted particles in the center-of-mass of the three interacting initial particles according to our previous sampling. In this sense, we can implement the forward and the backward reactions consistently within the same stochastic model.

In Ref. Oliinychenko et al. (2019) the deuteron reactions $\pi pn\rightarrow\pi d$ and $Npn\rightarrow Nd$ have been implemented in the SMASH transport approach for relativistic HICs. To do this numerically a fictitious $d^{\prime}$ resonance was introduced and the $3\rightarrow 2$ process was divided into a two 2-to-2 steps.

Later on, in Ref. Staudenmaier et al. (2021) the same multi-particle production mechanisms have been described according to the covariant rate formalism of Ref. Cassing (2002). In particular, Eq. (6) of Ref. Staudenmaier et al. (2021) shows the same probability for the stochastic treatment of the 3-to-2 process as the one we have just derived in Eq. (III). Comparing both expression we can make some comments:
i) The two- and three-body phase-spaces $R_{2}$ and $R_{3}$ appear in both equations as a function of $\sqrt{s}$ and particle masses. For $R_{2}$ we employ the well known analytic expression, while for $R_{3}$ we adopt the parametrization of Ref. Seifert (2018). We collect all formula in Appendix D. In Ref. Staudenmaier et al. (2021) it is done similarly, so we do not expect any discrepancy due to this part.
ii) The probability is proportional to the total cross section for the inverse 2-to-3 process. For deuteron disintegration into 3 particles by $\pi$ and $N$ scattering we employ a parametrization of the cross section as a function of $\sqrt{s}$, as reported in Appendix A, which is fitted to the available experimental data in the peak region. At high $\sqrt{s}$ we let our cross section to tend to zero because the 2-to-3 phase-space closes and other inelastic processes with final particles $m>3$ open. In Ref. Staudenmaier et al. (2021) cross sections from Ref. Oliinychenko et al. (2019) are used which differ only for a constant behavior at high $\sqrt{s}$. We investigated the possible difference arising from the different asymptotic behavior of the cross sections and we did not find any impact on deuteron yields and the $p_{T}$-spectra.
iii) Our spin factor $F_{spin}$ is the same as the one in Ref. Staudenmaier et al. (2021), while the isospin coefficient $F_{iso}$ does not appear there. This represents the novelty of our work.

In Ref. Staudenmaier et al. (2021) as from Ref. Oliinychenko et al. (2019) the $\pi$-catalysis is considered only for the channel where there is no charge difference between the initial and the final pion, i.e. $\pi d\leftrightarrow\pi pn$. We extend the deuteron production to all possible $\pi NN$ channels which fulfill the conservation of total isospin. We list all implemented channels in the table 1.

Since the deuteron has isospin zero, the state $\pi d$ is a state with defined isospin 1 provided by the pion. In general, a three particle state $\pi NN$ has not a definite value of total isospin (i.e. it is not an eigenstate of this quantum number), rather it is formed by a superposition of eigenstates. Therefore, for each channel of the table the $F_{iso}$ represents the projection of the state $\pi NN$ on the isospin 1 state. We perform the calculation in detail in Appendix C. For the inverse reaction $\pi d\rightarrow\pi NN$ the initial state has total isospin 1. The total cross section $\sigma_{tot}^{2,3}$ describes the reaction of $d\pi$ to any of the possible $\pi NN$ channels. In order to correctly evaluate the disintegration reaction, we weight the transition to one specific channel with the corresponding isospin factor calculated in Appendix C. Similarly we calculate the isospin factors $F_{iso}$ for the $N$-catalysis, where in this case there are no multiple channels available, i.e. $NNN\leftrightarrow dN$ does not account for other channels with respect to $Npn\leftrightarrow dN$.