PHQMD: Parton-Hadron Quantum Molecular Dynamics

• PHQMD is a microscopic n-body transport framework uniting PHSD and QMD to dynamically simulate hadron production and preserve baryonic correlations for nuclear clusters and hypernuclei.
• It integrates collision integrals with QMD-type baryonic dynamics, enabling precise modeling of heavy-ion reactions across energies from SIS to RHIC with various effective EoS implementations.
• PHQMD employs advanced cluster-recognition algorithms such as MST, SACA, and aMST, allowing for the dynamic emergence of clusters and hypernuclei during the collision evolution.

Searching arXiv for PHQMD papers and related recent work. Parton-Hadron-Quantum-Molecular Dynamics (PHQMD) is a microscopic $n$-body transport approach for relativistic heavy-ion collisions that was developed to describe, within a single framework, general hadron production together with the dynamical formation of nuclear clusters and hypernuclei (Aichelin et al., 2019, Kireyeu et al., 2019). Its defining construction is the combination of the collision integral and partonic/hadronic dynamics of Parton-Hadron-String Dynamics (PHSD) with Quantum Molecular Dynamics (QMD)-type propagation of the baryonic sector, so that many-body correlations and fluctuations needed for fragment formation are preserved during the reaction (Aichelin et al., 2019, Gläßel et al., 2021). In this formulation, cluster formation is not introduced solely as a freeze-out afterburner; rather, correlated baryonic aggregates emerge dynamically from interactions and can then be identified by algorithms such as Minimum Spanning Tree (MST), the Simulated Annealing Clusterization Algorithm (SACA), or later advanced MST variants used for deuterons (Kireyeu et al., 2019, Bratkovskaya et al., 2022, Coci et al., 2023).

1. Origin, scope, and conceptual motivation

PHQMD was introduced to address a specific deficiency of standard heavy-ion transport descriptions: bulk hadron spectra and collective observables can often be described, but light nuclei and hypernuclei require the retention of explicit $n$-body baryonic correlations during the reaction (Aichelin et al., 2019, Gläßel et al., 2021). The model is intended for energies from SIS and AGS through SPS, RHIC Beam Energy Scan, and RHIC energies, while retaining the PHSD treatment of hadronic reactions, strings, deconfinement, partonic transport, and hadronization (Aichelin et al., 2019, Gläßel et al., 2021).

In the authors’ formulation, PHQMD “unites the collision integrals of the Parton-Hadron-String Dynamics (PHSD) approach with 2-body potential interactions between baryons similar as in the Quantum Molecular Dynamics (QMD) approach where baryons are described by Gaussian wave functions” (Kireyeu et al., 2019). This architecture is central because PHSD already contains explicit partonic degrees of freedom, quarks and gluons, a lattice-QCD-based equation of state for the deconfined phase, dynamical hadronization, and hadronic elastic and inelastic rescattering, whereas PHQMD replaces baryon mean-field propagation by QMD-like $n$-body dynamics in the baryonic sector (Kireyeu et al., 2019, Gläßel et al., 2021).

The motivating physics domain is the baryon-rich regime relevant for NICA-MPD, BM@N, and CBM, framed in one overview as $\sqrt{s_{NN}} \lesssim 11$ GeV, where compression and decompression of baryonic matter coexist with substantial hadronic chemistry and threshold strangeness production (Kireyeu et al., 2019). This suggests that PHQMD was designed not only for bulk hadrochemistry but also for observables that depend on the space-time history of baryonic correlations, including cluster production, hyperon capture, hypernuclear channels, and the interpretation of weakly bound objects in a hot medium (Kireyeu et al., 2019, Bratkovskaya et al., 2022).

2. Microscopic structure and dynamical degrees of freedom

PHQMD is a hybrid construction. Baryons are propagated by QMD, while mesonic and partonic degrees of freedom continue to follow PHSD dynamics (Gläßel et al., 2021). Accordingly, the relevant propagated degrees of freedom depend on the reaction stage: strings and partons at early times, hadrons during and after hadronization, and correlated baryons in the hadronic stage from which clusters and hypernuclei can be identified (Kireyeu et al., 2019).

In the QMD sector, a baryon is represented by a Gaussian single-particle Wigner density. One PHQMD formulation writes

$$f ({\bf r}_i, {\bf p}_i,{\bf r}_{i0},{\bf p}_{i0},t) = \frac{1}{\pi^3 \hbar^3 } {\rm e}^{-\frac{2}{L} ({\bf r}_i - {\bf r}_{i0} (t) )^2 } {\rm e}^{-\frac{L}{2\hbar^2} ({\bf p}_i - {\bf p}_{i0} (t) )^2},$$

with $L=8.66\ {\rm fm}^2$ in the relativistic cluster-production study (Gläßel et al., 2021), whereas the original 2019 formal development used $L=2.16\ \mathrm{fm}^2$ (Aichelin et al., 2019). The centroid variables $({\bf r}_{i0}(t),{\bf p}_{i0}(t))$ are propagated by a variational QMD dynamics derived from the Dirac–Frenkel–McLachlan variational principle,

$$\delta \int_{t_1}^{t_2} dt \, \langle\psi(t)|i\frac{d}{dt}-H|\psi(t)\rangle = 0,$$

which yields Hamilton-like equations for the centroids (Aichelin et al., 2019, Gläßel et al., 2021).

The expectation-value Hamiltonian contains kinetic energy and effective interaction terms. In the Skyrme-based implementations, the baryon-baryon interaction is written as a sum of local Skyrme and Coulomb terms,

$$\begin{aligned} V_{i,j}&= V({\bf r}_i, {\bf r}_j,{\bf r}_{i0},{\bf r}_{j0},t) = V_{\rm Skyrme}+ V_{\rm Coul} \ &=\frac{1}{2} t_1 \delta ({\bf r}_i - {\bf r}_j) + \frac{1}{\gamma+1}t_2 \delta ({\bf r}_i - {\bf r}_j) \, \rho^{\gamma-1}(\cdots) +\frac{1}{2} \frac{Z_i Z_j e^2}{|{\bf r}_i-{\bf r}_j|}, \end{aligned}$$

with the associated density-functional form

$n$0

Here $n$1 is the interaction density generated from Gaussian overlaps of surrounding baryons, and $n$2 encode the nuclear equation of state (EoS) (Aichelin et al., 2019, Kireyeu et al., 2019, Gläßel et al., 2021).

For relativistic beam energies, PHQMD uses a pragmatic Lorentz-contracted modification of the baryon Wigner density along the beam axis rather than a fully covariant molecular dynamics (Aichelin et al., 2019, Gläßel et al., 2021). The PHSD inheritance remains substantial: string excitation and fragmentation are treated via FRITIOF 7.02 and PYTHIA 6.4, hadronic collision channels are those of PHSD/HSD, and when the local energy density exceeds $n$3 GeV/fm$n$4, the system can enter the partonic sector described by the Dynamical Quasi-Particle Model (DQPM) (Aichelin et al., 2019, Gläßel et al., 2021).

3. Equation of state and effective interactions

The EoS enters PHQMD through the density dependence of the effective baryonic interaction. In the original Skyrme implementations, the model employed static soft and hard EoS parameterizations. One detailed formulation gives the soft set as

$n$5

and the hard set as

$n$6

(Aichelin et al., 2019). A later synopsis also associates the soft set with $n$7 MeV and $n$8, and the hard set with $n$9 MeV and $n$0 (Kireyeu et al., 2019). In this context, a low compression modulus $n$1 corresponds to a soft EoS and a high $n$2 to a hard EoS (Kireyeu et al., 2019).

A major later extension introduced a momentum-dependent baryonic potential and a soft momentum-dependent EoS, calibrated to proton–nucleus elastic scattering through the Schrödinger-equivalent optical potential (Kireyeu et al., 2024). In that work, the optical-potential-based two-body interaction is parameterized as

$n$3

with $n$4, and is promoted to a local interaction by multiplication with $n$5 (Kireyeu et al., 2024). The three EoS realizations used there are soft (S), hard (H), and soft momentum-dependent (SM), with parameter sets:

EoS	Parameters
S	$n$6 MeV, $n$7 MeV, $n$8, $n$9 MeV
H	$\sqrt{s_{NN}} \lesssim 11$0 MeV, $\sqrt{s_{NN}} \lesssim 11$1 MeV, $\sqrt{s_{NN}} \lesssim 11$2, $\sqrt{s_{NN}} \lesssim 11$3 MeV
SM	$\sqrt{s_{NN}} \lesssim 11$4 MeV, $\sqrt{s_{NN}} \lesssim 11$5 MeV, $\sqrt{s_{NN}} \lesssim 11$6, $\sqrt{s_{NN}} \lesssim 11$7 MeV

The momentum-dependent parameters are $\sqrt{s_{NN}} \lesssim 11$8, $\sqrt{s_{NN}} \lesssim 11$9, $$f ({\bf r}_i, {\bf p}_i,{\bf r}_{i0},{\bf p}_{i0},t) = \frac{1}{\pi^3 \hbar^3 } {\rm e}^{-\frac{2}{L} ({\bf r}_i - {\bf r}_{i0} (t) )^2 } {\rm e}^{-\frac{L}{2\hbar^2} ({\bf p}_i - {\bf p}_{i0} (t) )^2},$$0 in the units stated in that work (Kireyeu et al., 2024).

The EoS sensitivity is phenomenologically significant across several PHQMD studies. In Au+Au at $$f ({\bf r}_i, {\bf p}_i,{\bf r}_{i0},{\bf p}_{i0},t) = \frac{1}{\pi^3 \hbar^3 } {\rm e}^{-\frac{2}{L} ({\bf r}_i - {\bf r}_{i0} (t) )^2 } {\rm e}^{-\frac{L}{2\hbar^2} ({\bf p}_i - {\bf p}_{i0} (t) )^2},$$1 GeV, a hard EoS increases the slope of hadron spectra at large $$f ({\bf r}_i, {\bf p}_i,{\bf r}_{i0},{\bf p}_{i0},t) = \frac{1}{\pi^3 \hbar^3 } {\rm e}^{-\frac{2}{L} ({\bf r}_i - {\bf r}_{i0} (t) )^2 } {\rm e}^{-\frac{L}{2\hbar^2} ({\bf p}_i - {\bf p}_{i0} (t) )^2},$$2 compared with a soft EoS, while the proton slope remains slightly underestimated at large $$f ({\bf r}_i, {\bf p}_i,{\bf r}_{i0},{\bf p}_{i0},t) = \frac{1}{\pi^3 \hbar^3 } {\rm e}^{-\frac{2}{L} ({\bf r}_i - {\bf r}_{i0} (t) )^2 } {\rm e}^{-\frac{L}{2\hbar^2} ({\bf p}_i - {\bf p}_{i0} (t) )^2},$$3 (Kireyeu et al., 2019). At SIS energies, flow observables show that a soft momentum-dependent EoS can behave more like a hard static EoS than like a static soft EoS in $$f ({\bf r}_i, {\bf p}_i,{\bf r}_{i0},{\bf p}_{i0},t) = \frac{1}{\pi^3 \hbar^3 } {\rm e}^{-\frac{2}{L} ({\bf r}_i - {\bf r}_{i0} (t) )^2 } {\rm e}^{-\frac{L}{2\hbar^2} ({\bf p}_i - {\bf p}_{i0} (t) )^2},$$4 and $$f ({\bf r}_i, {\bf p}_i,{\bf r}_{i0},{\bf p}_{i0},t) = \frac{1}{\pi^3 \hbar^3 } {\rm e}^{-\frac{2}{L} ({\bf r}_i - {\bf r}_{i0} (t) )^2 } {\rm e}^{-\frac{L}{2\hbar^2} ({\bf p}_i - {\bf p}_{i0} (t) )^2},$$5, despite sharing $$f ({\bf r}_i, {\bf p}_i,{\bf r}_{i0},{\bf p}_{i0},t) = \frac{1}{\pi^3 \hbar^3 } {\rm e}^{-\frac{2}{L} ({\bf r}_i - {\bf r}_{i0} (t) )^2 } {\rm e}^{-\frac{L}{2\hbar^2} ({\bf p}_i - {\bf p}_{i0} (t) )^2},$$6 MeV with the soft case (Kireyeu et al., 2024). This suggests that compressibility alone does not exhaust the EoS dependence of PHQMD predictions.

4. Cluster-recognition strategies and dynamical cluster formation

A central feature of PHQMD is the distinction between dynamical formation and algorithmic recognition. The interactions among baryons generate correlated aggregates during the transport evolution; MST, SACA, or later advanced MST procedures identify those aggregates in the event record rather than producing them as an external afterburner prescription (Kireyeu et al., 2019, Gläßel et al., 2021, Kireyeu et al., 2024).

MST is the simplest recognition method. In the midrapidity-cluster studies, clusters are identified in coordinate space such that a proton-neutron pair belongs to a cluster when their distance in their rest frame is smaller than

$$f ({\bf r}_i, {\bf p}_i,{\bf r}_{i0},{\bf p}_{i0},t) = \frac{1}{\pi^3 \hbar^3 } {\rm e}^{-\frac{2}{L} ({\bf r}_i - {\bf r}_{i0} (t) )^2 } {\rm e}^{-\frac{L}{2\hbar^2} ({\bf p}_i - {\bf p}_{i0} (t) )^2},$$7

and additional momentum-space cuts were reported not to change the cluster distributions significantly (Bratkovskaya et al., 2022, Gläßel et al., 2021). Earlier PHQMD work on light clusters quoted an MST coordinate-space condition

$$f ({\bf r}_i, {\bf p}_i,{\bf r}_{i0},{\bf p}_{i0},t) = \frac{1}{\pi^3 \hbar^3 } {\rm e}^{-\frac{2}{L} ({\bf r}_i - {\bf r}_{i0} (t) )^2 } {\rm e}^{-\frac{L}{2\hbar^2} ({\bf p}_i - {\bf p}_{i0} (t) )^2},$$8

with the interpretation that MST is most reliable late in the reaction when fragments are already spatially separated (Aichelin et al., 2019). A plausible implication is that the concrete MST implementation evolved across PHQMD applications.

SACA is the more elaborate recognition method emphasized in the original cluster and hypernucleus overviews. At a given time $$f ({\bf r}_i, {\bf p}_i,{\bf r}_{i0},{\bf p}_{i0},t) = \frac{1}{\pi^3 \hbar^3 } {\rm e}^{-\frac{2}{L} ({\bf r}_i - {\bf r}_{i0} (t) )^2 } {\rm e}^{-\frac{L}{2\hbar^2} ({\bf p}_i - {\bf p}_{i0} (t) )^2},$$9, an initial pre-clustering is first obtained with MST, after which clusters and single nucleons are recombined into candidate partitions. Among these, the most bound configuration is sought stochastically by a Metropolis simulated-annealing procedure, and repeated iteration automatically leads to the most bound configuration (Kireyeu et al., 2019). SACA is based on the idea of Dorso and Randrup that the most bound configuration at early times evolves toward the final fragment distribution, so it can identify clusters before geometric separation is complete (Kireyeu et al., 2019, Bratkovskaya et al., 2019, Aichelin et al., 2019).

Later deuteron studies refined MST into an advanced Minimum Spanning Tree (aMST) procedure. There the coordinate-space condition remains

$L=8.66\ {\rm fm}^2$0

and the cluster must have negative binding energy,

$L=8.66\ {\rm fm}^2$1

with the cluster binding energy in its rest frame written as

$L=8.66\ {\rm fm}^2$2

The aMST stabilization logic freezes internal cluster degrees of freedom once the cluster is no longer in contact, by collision or potential interaction, with hadrons outside the cluster, and thereby removes the need for an ad hoc cluster-recognition time in deuteron analyses (Coci et al., 2023). A 2024 conference contribution describes a related “advanced Minimum Spanning Tree (aMST)” algorithm as identifying bound clusters with negative binding energies by correlations of baryons in coordinate space (Coci et al., 2024).

The contrast with coalescence models is explicit in the PHQMD literature. Coalescence assembles nuclei at freeze-out by cuts in relative coordinate and momentum space, whereas PHQMD aims to preserve the interactions and correlations that make clusters emerge dynamically during the evolution (Aichelin et al., 2019, Bratkovskaya et al., 2022). At the same time, PHQMD has also been used as a host environment for embedded coalescence, enabling same-event comparisons between dynamical clustering and freeze-out coalescence prescriptions (Bratkovskaya et al., 2022, Kireyeu et al., 2024).

5. Hypernuclei, deuteron mechanisms, and the “ice in the fire” problem

Hypernuclei in PHQMD arise by extending the clusterization logic into the strange sector. Hyperons are produced through PHSD string, hadronic, or hadronization channels, propagate in the hadronic medium, and can become bound to nucleonic clusters; in the model’s own wording, “the capture of the produced hyperons by clusters of nucleons leads to the hypernuclei formation” (Kireyeu et al., 2019). In one relativistic cluster study, the hyperon–nucleon interaction is approximated by

$L=8.66\ {\rm fm}^2$3

which is explicitly described as a rough approximation and a limitation of the present hypernuclear treatment (Gläßel et al., 2021).

The deuteron sector became a major area of PHQMD development. A 2023 study formulated two microscopic mechanisms for deuteron production at midrapidity: “kinetic” deuterons produced by hadronic reactions and “potential” deuterons produced by the attractive interaction between nucleons and recognized by aMST (Coci et al., 2023). The explicit kinetic channels include

$L=8.66\ {\rm fm}^2$4

with subdominant $L=8.66\ {\rm fm}^2$5, plus elastic $L=8.66\ {\rm fm}^2$6 and $L=8.66\ {\rm fm}^2$7 channels in the full transport implementation (Coci et al., 2023). A 2024 proceeding likewise states that PHQMD includes catalytic hadronic reactions

$L=8.66\ {\rm fm}^2$8

for deuterons, while anti-clusters are produced by potential interactions of anti-baryons and recognized by aMST in an analogous way (Coci et al., 2024).

A key refinement in the kinetic deuteron treatment is the inclusion of the quantum nature of the deuteron in both coordinate and momentum space. PHQMD imposes a finite-size excluded-volume condition around the candidate deuteron and projects the relative momentum of the interacting nucleon pair onto the deuteron wave function in momentum space (Coci et al., 2023). In that study, the excluded radius is

$L=8.66\ {\rm fm}^2$9

while the text also mentions a deuteron rms radius $L=2.16\ \mathrm{fm}^2$0 (Coci et al., 2023). The combined effect is strong suppression of kinetic deuteron production in a dense medium: each of the excluded-volume and momentum-projection effects suppresses the yield at midrapidity by about a factor of $L=2.16\ \mathrm{fm}^2$1, both together give roughly another factor of $L=2.16\ \mathrm{fm}^2$2, the total suppression is about $L=2.16\ \mathrm{fm}^2$3 at $L=2.16\ \mathrm{fm}^2$4 GeV, and at top RHIC the suppression can reach an order of magnitude (Coci et al., 2023).

These time-resolved and geometry-resolved analyses underpin the PHQMD explanation of the “ice in the fire” puzzle. PHQMD studies report that deuterons remain at smaller transverse radii than free nucleons and “follow behind the front of the expanding baryonic fireball,” becoming spatially separated from the hotter and more collision-active front (Bratkovskaya et al., 2022). More generally, the baryons that end up in clusters are found, on average, closer to the reaction center and effectively behind the front of the fast expanding hadrons (Gläßel et al., 2021). This suggests that the survival of weakly bound nuclei is attributed in PHQMD not to emission from a homogeneous thermal source but to late formation in special spatial regions where collisions have largely ceased.

6. Phenomenology, comparisons to experiment, and limitations

PHQMD has been benchmarked against both elementary and nucleus–nucleus data. In elementary inelastic $L=2.16\ \mathrm{fm}^2$5 collisions, PHQMD with the PHSD-tuned string sector reproduces the energy dependence of $L=2.16\ \mathrm{fm}^2$6 and $L=2.16\ \mathrm{fm}^2$7 multiplicities reasonably well and provides predictions for $L=2.16\ \mathrm{fm}^2$8 and $L=2.16\ \mathrm{fm}^2$9 (Kireyeu et al., 2019). In Au+Au at $({\bf r}_{i0}(t),{\bf p}_{i0}(t))$0 GeV, PHQMD reproduces the transverse-momentum spectra of $({\bf r}_{i0}(t),{\bf p}_{i0}(t))$1, $({\bf r}_{i0}(t),{\bf p}_{i0}(t))$2, $({\bf r}_{i0}(t),{\bf p}_{i0}(t))$3, $({\bf r}_{i0}(t),{\bf p}_{i0}(t))$4, $({\bf r}_{i0}(t),{\bf p}_{i0}(t))$5, and $({\bf r}_{i0}(t),{\bf p}_{i0}(t))$6 across centrality classes reasonably well, while the proton slope remains slightly underestimated at high $({\bf r}_{i0}(t),{\bf p}_{i0}(t))$7 (Kireyeu et al., 2019).

For clusters and hypernuclei, the phenomenology spans AGS, SPS, RHIC-BES, RHIC fixed-target, and SIS energies. At AGS energy $({\bf r}_{i0}(t),{\bf p}_{i0}(t))$8GeV, PHQMD describes E864 deuteron, triton, and $({\bf r}_{i0}(t),{\bf p}_{i0}(t))$9He rapidity and $$\delta \int_{t_1}^{t_2} dt \, \langle\psi(t)|i\frac{d}{dt}-H|\psi(t)\rangle = 0,$$0 spectra rather well and reproduces coalescence-like observables $$\delta \int_{t_1}^{t_2} dt \, \langle\psi(t)|i\frac{d}{dt}-H|\psi(t)\rangle = 0,$$1, $$\delta \int_{t_1}^{t_2} dt \, \langle\psi(t)|i\frac{d}{dt}-H|\psi(t)\rangle = 0,$$2, and the penalty factor up to $$\delta \int_{t_1}^{t_2} dt \, \langle\psi(t)|i\frac{d}{dt}-H|\psi(t)\rangle = 0,$$3 (Gläßel et al., 2021). At SPS energy $$\delta \int_{t_1}^{t_2} dt \, \langle\psi(t)|i\frac{d}{dt}-H|\psi(t)\rangle = 0,$$4 GeV, it reproduces NA49 deuteron and $$\delta \int_{t_1}^{t_2} dt \, \langle\psi(t)|i\frac{d}{dt}-H|\psi(t)\rangle = 0,$$5He rapidity and transverse-momentum spectra and the increase of cluster-production probability with transverse momentum (Gläßel et al., 2021). Across RHIC-BES energies up to top RHIC, PHQMD reproduces the excitation functions of proton, antiproton, and deuteron midrapidity yields reasonably well, together with deuteron $$\delta \int_{t_1}^{t_2} dt \, \langle\psi(t)|i\frac{d}{dt}-H|\psi(t)\rangle = 0,$$6 spectra and $$\delta \int_{t_1}^{t_2} dt \, \langle\psi(t)|i\frac{d}{dt}-H|\psi(t)\rangle = 0,$$7 systematics, although it overpredicts deuterons at $$\delta \int_{t_1}^{t_2} dt \, \langle\psi(t)|i\frac{d}{dt}-H|\psi(t)\rangle = 0,$$8 GeV and overpredicts $$\delta \int_{t_1}^{t_2} dt \, \langle\psi(t)|i\frac{d}{dt}-H|\psi(t)\rangle = 0,$$9 by about a factor of two if antinuclei are read out at the same time as nuclei (Gläßel et al., 2021).

Conference-overview studies of midrapidity cluster formation reported that proton and antiproton midrapidity yields and deuteron $$\begin{aligned} V_{i,j}&= V({\bf r}_i, {\bf r}_j,{\bf r}_{i0},{\bf r}_{j0},t) = V_{\rm Skyrme}+ V_{\rm Coul} \ &=\frac{1}{2} t_1 \delta ({\bf r}_i - {\bf r}_j) + \frac{1}{\gamma+1}t_2 \delta ({\bf r}_i - {\bf r}_j) \, \rho^{\gamma-1}(\cdots) +\frac{1}{2} \frac{Z_i Z_j e^2}{|{\bf r}_i-{\bf r}_j|}, \end{aligned}$$0 are well reproduced from the lowest SPS energies to the highest RHIC energies, that $$\begin{aligned} V_{i,j}&= V({\bf r}_i, {\bf r}_j,{\bf r}_{i0},{\bf r}_{j0},t) = V_{\rm Skyrme}+ V_{\rm Coul} \ &=\frac{1}{2} t_1 \delta ({\bf r}_i - {\bf r}_j) + \frac{1}{\gamma+1}t_2 \delta ({\bf r}_i - {\bf r}_j) \, \rho^{\gamma-1}(\cdots) +\frac{1}{2} \frac{Z_i Z_j e^2}{|{\bf r}_i-{\bf r}_j|}, \end{aligned}$$1 is reproduced rather well, and that $$\begin{aligned} V_{i,j}&= V({\bf r}_i, {\bf r}_j,{\bf r}_{i0},{\bf r}_{j0},t) = V_{\rm Skyrme}+ V_{\rm Coul} \ &=\frac{1}{2} t_1 \delta ({\bf r}_i - {\bf r}_j) + \frac{1}{\gamma+1}t_2 \delta ({\bf r}_i - {\bf r}_j) \, \rho^{\gamma-1}(\cdots) +\frac{1}{2} \frac{Z_i Z_j e^2}{|{\bf r}_i-{\bf r}_j|}, \end{aligned}$$2 has the correct shape as a function of $$\begin{aligned} V_{i,j}&= V({\bf r}_i, {\bf r}_j,{\bf r}_{i0},{\bf r}_{j0},t) = V_{\rm Skyrme}+ V_{\rm Coul} \ &=\frac{1}{2} t_1 \delta ({\bf r}_i - {\bf r}_j) + \frac{1}{\gamma+1}t_2 \delta ({\bf r}_i - {\bf r}_j) \, \rho^{\gamma-1}(\cdots) +\frac{1}{2} \frac{Z_i Z_j e^2}{|{\bf r}_i-{\bf r}_j|}, \end{aligned}$$3 but is overestimated by about a factor of two because of overpredicted $$\begin{aligned} V_{i,j}&= V({\bf r}_i, {\bf r}_j,{\bf r}_{i0},{\bf r}_{j0},t) = V_{\rm Skyrme}+ V_{\rm Coul} \ &=\frac{1}{2} t_1 \delta ({\bf r}_i - {\bf r}_j) + \frac{1}{\gamma+1}t_2 \delta ({\bf r}_i - {\bf r}_j) \, \rho^{\gamma-1}(\cdots) +\frac{1}{2} \frac{Z_i Z_j e^2}{|{\bf r}_i-{\bf r}_j|}, \end{aligned}$$4 (Bratkovskaya et al., 2022). At $$\begin{aligned} V_{i,j}&= V({\bf r}_i, {\bf r}_j,{\bf r}_{i0},{\bf r}_{j0},t) = V_{\rm Skyrme}+ V_{\rm Coul} \ &=\frac{1}{2} t_1 \delta ({\bf r}_i - {\bf r}_j) + \frac{1}{\gamma+1}t_2 \delta ({\bf r}_i - {\bf r}_j) \, \rho^{\gamma-1}(\cdots) +\frac{1}{2} \frac{Z_i Z_j e^2}{|{\bf r}_i-{\bf r}_j|}, \end{aligned}$$5 GeV, PHQMD describes the $$\begin{aligned} V_{i,j}&= V({\bf r}_i, {\bf r}_j,{\bf r}_{i0},{\bf r}_{j0},t) = V_{\rm Skyrme}+ V_{\rm Coul} \ &=\frac{1}{2} t_1 \delta ({\bf r}_i - {\bf r}_j) + \frac{1}{\gamma+1}t_2 \delta ({\bf r}_i - {\bf r}_j) \, \rho^{\gamma-1}(\cdots) +\frac{1}{2} \frac{Z_i Z_j e^2}{|{\bf r}_i-{\bf r}_j|}, \end{aligned}$$6 spectra of $$\begin{aligned} V_{i,j}&= V({\bf r}_i, {\bf r}_j,{\bf r}_{i0},{\bf r}_{j0},t) = V_{\rm Skyrme}+ V_{\rm Coul} \ &=\frac{1}{2} t_1 \delta ({\bf r}_i - {\bf r}_j) + \frac{1}{\gamma+1}t_2 \delta ({\bf r}_i - {\bf r}_j) \, \rho^{\gamma-1}(\cdots) +\frac{1}{2} \frac{Z_i Z_j e^2}{|{\bf r}_i-{\bf r}_j|}, \end{aligned}$$7 well, underestimates the rapidity distribution of light clusters, reproduces the shape of the $$\begin{aligned} V_{i,j}&= V({\bf r}_i, {\bf r}_j,{\bf r}_{i0},{\bf r}_{j0},t) = V_{\rm Skyrme}+ V_{\rm Coul} \ &=\frac{1}{2} t_1 \delta ({\bf r}_i - {\bf r}_j) + \frac{1}{\gamma+1}t_2 \delta ({\bf r}_i - {\bf r}_j) \, \rho^{\gamma-1}(\cdots) +\frac{1}{2} \frac{Z_i Z_j e^2}{|{\bf r}_i-{\bf r}_j|}, \end{aligned}$$8 and $$\begin{aligned} V_{i,j}&= V({\bf r}_i, {\bf r}_j,{\bf r}_{i0},{\bf r}_{j0},t) = V_{\rm Skyrme}+ V_{\rm Coul} \ &=\frac{1}{2} t_1 \delta ({\bf r}_i - {\bf r}_j) + \frac{1}{\gamma+1}t_2 \delta ({\bf r}_i - {\bf r}_j) \, \rho^{\gamma-1}(\cdots) +\frac{1}{2} \frac{Z_i Z_j e^2}{|{\bf r}_i-{\bf r}_j|}, \end{aligned}$$9 $n$00 spectra, overestimates the differential yield of $n$01, and reproduces the yield of $n$02 (Bratkovskaya et al., 2022).

At SIS energies, the focus shifted to flow. The momentum-dependent EoS study finds that the best description of HADES and FOPI data on the directed and elliptic flow coefficients of protons and light clusters is obtained with the momentum-dependent EoS, and it also reports a scaling behavior of $n$03 versus $n$04 with atomic number $n$05 (Kireyeu et al., 2024). The same work argues that flow observables can help identify cluster production mechanisms because deuterons from kinetic plus MST production differ from deuterons formed by coalescence when analyzed on the same underlying events (Kireyeu et al., 2024).

PHQMD’s limitations are stated repeatedly. The approach is semiclassical, so clusters are not exact quantum ground states and can “spontaneously” lose nucleons; cluster multiplicities therefore depend somewhat on observation time unless stabilization procedures are applied (Bratkovskaya et al., 2022, Coci et al., 2023). Hypernuclear yields remain sensitive to the poorly constrained $n$06 interaction, often represented only by $n$07 (Gläßel et al., 2021). Earlier PHQMD versions used static momentum-independent Skyrme interactions and approximate relativistic extensions of QMD (Aichelin et al., 2019, Coci et al., 2023). Computational cost is also substantial: one study notes that extending cluster calculations from $n$08 to $n$09 would require roughly three orders of magnitude more CPU time (Gläßel et al., 2021).

Taken together, the PHQMD literature presents the model as a transport framework in which heavy-ion collisions, partonic and hadronic reaction dynamics, baryonic mean-field and correlation effects, and the emergence of nuclear clusters and hypernuclei are treated within a common microscopic scheme (Aichelin et al., 2019, Kireyeu et al., 2019). Its enduring significance lies in the attempt to make cluster and hypernucleus formation a consequence of the transport dynamics itself rather than a purely external prescription, while retaining contact with bulk hadron phenomenology from SIS and AGS energies through SPS and RHIC (Gläßel et al., 2021, Bratkovskaya et al., 2022).