Towards a Unified Quark-Hadron Matter Equation of State for Applications in Astrophysics and Heavy-Ion Collisions

We outline an approach to a unified equation of state for quark-hadron matter on the basis of a $\Phi-$derivable approach to the generalized Beth-Uhlenbeck equation of state for a cluster decomposition of thermodynamic quantities like the density. To this end we summarize the cluster virial expansion for nuclear matter and demonstrate the equivalence of the Green's function approach and the $\Phi-$derivable formulation. For an example, the formation and dissociation of deuterons in nuclear matter is discussed. We formulate the cluster $\Phi-$derivable approach to quark-hadron matter which allows to take into account the specifics of chiral symmetry restoration and deconfinement in triggering the Mott-dissociation of hadrons. This approach unifies the description of a strongly coupled quark-gluon plasma with that of a medium-modified hadron resonance gas description which are contained as limiting cases. The developed formalism shall replace the common two-phase approach to the description of the deconfinement and chiral phase transition that requires a phase transition construction between separately developed equations of state for hadronic and quark matter phases. Applications to the phenomenology of heavy-ion collisions and astrophysics are outlined.


Introduction
The development of a unified equation of state (EoS) for quark-hadron matter, where the hadrons are not elementary degrees of freedom but rather appear as composites of their quark constituents (bound and scattering states of effective quark interactions) is a formidable task because it requires the implementation of dynamical mechanisms of quark confinement as well as chiral symmetry breaking which characterise the hadronic phase of matter and mechanisms for deconfinement and chiral symmetry restoration which determine the transition to the quark(-gluon) phase of strongly interacting matter.
The aspect of bound state formation and dissociation has been well studied in warm dense nuclear matter, where a cluster virial expansion has been developed as a most effective means of description of a system that contains clusters of different sizes with internal quantum numbers like in the nuclear statistical equilibrium picture, but generalizes it by including residual binary interactions among them (second cluster virial coefficient) and their dissociation due to compression and heating within the Mott effect. The consistent description of bound and scattering states on the same footing is provided by the generalized Beth-Uhlenbeck equation of state that uses phase shifts in order to describe correlations and their modifications in a hot, dense environment. It turns out that for the short-ranged strong interactions it is not the screening of the interaction which drives the dissociation of the bound states but rather the Pauli blocking effect that inhibits cluster formation when the phase space is densely populated. Being based on pure symmetry arguments, this effect is sufficiently general to apply to any fermionic many-particle system where-as a function of the density-the formation of bound state gets first favored (principle of Le Chatelier and Brown) until the Pauli principle inhibits scattering processes that would lead to cluster formation and a homogeneous phase of fermionic quasiparticles emerges: The dense nuclear matter.
Since nucleons themselves can be considered as clusters of quarks, the more fundamental degrees of freedom of quantum chromodynamics (QCD) as the underlying gauge field theory of strong interactions, we want to develop in this work the basics of a corresponding approach that will describe the transition from hadronic to quark matter as a Mott transition driven by the Pauli blocking effect. Despite this striking analogy there are also specifics for the quark degrees of freedom which need to be taken into account. These are mainly the internal quantum numbers (color, flavor, spin) that cause the confinement of colored states (quarks, gluons, diquarks etc.), the dynamical chiral symmetry breaking and combined symmetry requirements as well as relativistic kinematics.
In Figure 1 we show the phase diagram of QCD according to the present state of knowledge. In particular it addresses the aspect of cluster formation and dissociation in dense quark-hadron matter, which is the main goal of the approach to be outlined in this work.
To this end in Section 2 we first review the cluster virial approach for nuclear matter in the form of a generalized Beth-Uhlenbeck EoS. Then in Section 3 we suggest that it might be obtained from the Φ−derivable approach in its field theoretic formulation when for the Φ functional the class of all two-loop Feynman diagramns is chosen that can be constructed from cluster Green's functions and cluster T-matrices and consequently take the form of generalized "sunset" type diagrams.
In Section 4 we develop the approach for quark matter with mesons and baryons as clusters dominating the hadronic phase of matter since quarks and diquarks are suppressed by confining interactions that give rise to diverging selfenergies for those states in the low-density (confinement) phase. This mechanism is realized here within the so-called string-flip model that is generalized in relativistic form. As a first step, a selfconsistent mean-field approximation is performed which results in a new quark matter equation of state with confining features, thus superior to previous NJL model approaches and their Polyakov-loop generalized versions. We demonstrate how the quark Pauli-blocking effect for hadrons is already inherent in the Φ−derivable approach and can be made apparent by a perturbation expansion with respect to cluster selfenergies on top of the quasiparticle approximation. The Pauli-blocking effect can be mimicked by a hadronic excluded volume which should therefore be taken into account when a hadronic EoS is extrapolated from a low-density limit (hadron resonance gas, nuclear statistical equilibrium) to the vicinity of the deconfinement transition. In the quark matter phase, we implement higher order repulsive interactions which may be justified by multi-pomeron exchange interactions, thus non-perturbative effects of the gluon sector which we are not treating dynamically at this stage. These effects cause a stiffening of the high-density quark matter EoS that are essential for hybrid compact star applications.
In Section 5 we discuss the role of nuclear clusters for the astrophysics of supernovae and the deconfinement transition in compact stars while in Section 6 we consider the potential of a unified approach to quark-hadron matter for applications in heavy-ion collisions. We point out that nuclear cluster formation may occur directly at the hadronisation transition and thus result in their (sudden) chemical freeze-out together with other hadronic species as indicated by the ALICE experiment at the CERN LHC.
In Section 7 we draw conclusions for the further development of the approach towards a unified description of quark-hadron matter and its phenomenology in heavy-ion collisions and astrophysics.  [1], the position of the critical endpoint of the deconfinement transition is not known. The dashed line indicates the pseudocritical temperature T c (µ B ) of the crossover transition from hadronic to quark matter which has been determined at vanishing baryochemical potential µ B = 0 in lattice QCD simulations to T c (0) = 154 ± 9 MeV [2]. The red and white colored circles stand for neutrons and protons which at low densities can form nuclear clusters that get dissociated in the liquid phase. Similar to that case, nucleons themselves can be viewed as clusters of quarks symbolized by colored dots that get dissociated in the quark matter phase. The green shaded area is the domain of the phase diagram accessible by supernova explosions. The red shaded area is the region accessible by ab-initio simulations of QCD at finite temperatures and chemical potentials on space-time lattices. The green arrowed lines symbolize trajectories of heavy-ion collisions at different center of mass energies ranging form LHC over RHIC to the planned NICA/FAIR experiments. The green line on the chemical potential axis depicts the range of values accessible in compact star interiors. In this case, the third dimension of the charge (or isospin) chemical potential relevant for isospin asymmetric systems in supernovae, compact stars and their mergers is suppressed.

Quantum Statistical Approach and the Cluster Virial Expansion
A systematic approach to the cluster expansion of thermodynamic properties is obtained from quantum statistics. We consider a system consisting of species i with the (conserved) particle number N i and the corresponding chemical potential µ i , described by the hamiltonian H at equilibrium with temperature T. The index i will be dropped in the following. The grand canonical thermodynamic potential is given by where P is the pressure and V the volume. It can be represented by diagrams within a perturbation expansion [3,4], see also [5]. We have or (3) where λ is a scaling factor substituting the two-particle interaction V 2 by λ V 2 . G (0) 1 is the free single-particle propagator that gives the ideal part of the pressure P 0 . The full single-particle Green function G λ and the self-energy Σ λ are taken with the coupling constant λ. Depending on the selected diagrams, different approximations can be found. In particular, the second virial coefficient for charged particle systems has been investigated, see Ref. [5]. 1 An alternative way to derive the equation of state is to start from the expression for the total nucleon density where V is the system volume, the variable 1 = {p 1 , σ 1 , τ 1 } denotes the single particle states (here nucleons in momentum representation), τ 1 = n, p, and summation over spin direction is collected in the factor 2. Both the Fermi distribution function and the spectral function depend on the temperature and the chemical potentials µ p , µ n , which are not given explicitly. The spectral function S 1 (1, ω) of the single-particle Green function G 1 (1, iz ν ) is related to the single-particle self-energy Σ(1, z) according to where the imaginary part has to be taken for a small negative imaginary part in the frequency; E (0) (1) = p 2 1 /(2m 1 ). Both approaches are equivalent. The perturbation expansion can be represented by Feynman diagrams, see Equation (3). However, a finite order perturbation theory will not produce bound states, and partial summations of an infinite number of e.g., ladder diagrams must be performed to get bound states. As shown by Baym and Kadanoff [3,4], self-consistent approximations to the one-particle Green function can be given based on a functional Φ so that Different approximations for the generating functional Φ are discussed in the following Sections 3 and 4. The self-consistent Φ−derivable approximations not only lead to a fully-conserving transport theory. In the equilibrium case they also have the property that different methods to obtain the grand partition function such as integrating the expectation value of the potential energy with respect to the coupling constant λ, Equation (3), or integrating the density n with respect to the chemical potential µ, lead to the same result [6,7]. In particular, with also n = − ∂Ω ∂µ (8) holds in the considered approximation. The latter approach (using Equation (4)) has been extensively used in many-particle systems [8][9][10][11][12], in particular in connection with the chemical picture. The main idea of the chemical picture is to treat bound states on the same footing as "elementary" single particles. This way one describes correctly 1 Note that the interaction in nuclear systems is strong. However, the perturbation expansion is performed with respect to the imaginary part of the self-energy that is assumed to be small. Most of the interaction is already taken into account in the self-consistent determination of the quasiparticle energies. With increasing density, the Fermi energy will dominate the potential energy so that the correlations are suppressed. A quasiparticle description can be used to calculate the nuclear structure. the low-temperature, low-density region of many-body systems where bound states dominate. Within a quantum statistical approach, the chemical picture is realized considering the A-particle propagator. In ladder approximation, see Figure 2, the Bethe-Salpeter equation (BSE) is obtained, where z A is the A-particle Matsubara frequency, G 0 A is the product of single-particle propagators, The BSE is equivalent to the A-particle wave equation. Neglecting all medium effects we have where ν indicates the internal quantum number of the A-particle eigen states, including the channel c describing spin and isospin state. Different excitations are possible in each channel c, in particular bound states and scattering states. The chemical picture uses the eigen-representation of the A-particle propagator. In ladder approximation, neglecting medium effects, In particular, it contains the contribution of A-particle bound states (nuclei) similar to the "elementary" single particle propagator Within the representation of the perturbative expansion by Feynman diagrams, the chemical picture implies that diagrams containing single-particle propagators should be completed by adding A-particle propagators, at least for the bound states ν. However, also the scattering parts ν have to be considered for a full description.
Coming back to Equation (5), different approximations for the self-energy Σ 1 (1, z) can be considered. The main issues are the implementation of bound state formation which is of relevance in the low-density region, and the account for density effect if considering higher densities, such as mean-field approximations. To give some systematics with respect to the different approximations calculating the self-energy and the corresponding approximation for the EoS, in Ref. [13] the following overview was presented.
As seen in Table 1, to consider density effects, the ideal gas approximation is improved by the mean-field approximation. Quasiparticle shifts are introduced in a semi-empirical way within the relativistic mean-field (RMF) approximation 2 or considering Skyrme parametrizations. 2 A Padé approximation of the nucleon quasiparticle shifts, applicable in a wide temperature range, can be found in Ref. [14]. This approximation is assumed to give an adequate description of dense matter near and above the saturation density but fails to describe correlations, in particular bound state formation, in the low density region. Table 1. Aspects of the quantum statistical description of a many-particle system with bound and scattering states in the low-density limit (left column) and its modification due to medium effects at high densities (right column).
low density limit high density modification (medium effects) For this, the chemical picture is necessary, which is realized by a cluster decomposition of the self-energy, see [10]. In particular, considering only the bound state part of the free A-particle propagator, the nuclear statistical equilibrium (NSE) is obtained. However, also the scattering part of the free A-particle propagator must be considered to obtain the correct second virial coefficient in the case A = 2, as given by the Beth-Uhlenbeck formula [15]. A cluster Beth-Uhlenbeck formula has been worked out [16] which considers the chemical picture, where scattering states between two arbitrary clusters A and B are taken into account.
The inclusion of density effects for these approximations, as collected in Table 1, is not simple. A generalized Beth-Uhlenbeck formula has been worked out in Ref. [12]. In this approach, the Pauli blocking and the dissolution of bound states with increasing density is shown.
A more detailed approach should also take into account that the medium is not considered as uncorrelated, but correlations in the medium have to be taken into account. The cluster-mean-field (CMF) approximation has been worked out [11] where also clustering of the medium is treated. This is of relevance in the region where correlations and bound state formation dominate.
These effects have been considered recently [17] in warm dense matter. Improving the NSE, medium modifications of the bound states are discussed. For a consistent description, in particular of the correct second virial coefficient, scattering phase shifts have been implemented. If clustering is relevant, the medium cannot be approximated by an uncorrelated, ideal Fermion gas of quasiparticles, but correlations must be taken into account.
Here we focus on the cluster virial expansion [16]. We consider density effects such as single-particle quasiparticle shifts and Pauli blocking, which are responsible for the dissolution of bound states and the appearance of a new state of matter. In the present section we treat nuclear clusters which are dissolved to form a nuclear Fermi liquid when density increases. Later on we go a step forward to investigate the dissolution of hadrons by medium effects to form the quark-gluon plasma.
As shown in [17], using the cluster decomposition of the self-energy which takes into account, in particular, cluster formation, we obtain where P denotes the center of mass (c.o.m.) momentum of the cluster (or, for A = 1, the momentum of the nucleon). The internal quantum state ν contains the proton number Z and neutron number N = A − Z of the cluster, is the Bose or Fermi distribution function for even or odd A, respectively, that is depending on {T, µ n , µ p }. The integral over ω is performed within the quasiparticle approach, the quasiparticle energies E A,ν (P; T, µ n , µ p ) are depending on the medium characterized by {T, µ n , µ p }. Expressions for the in-medium modifications are given in [17]. In Equation (13) the sum is to be taken over the mass number A of the cluster, the center-of-mass momentum P, and the intrinsic quantum number ν. The summation over ν concerns the bound states as far as they exist, as well as the continuum of scattering states. Solving the few-body problem what is behind the calculation of the A-nucleon T matrices in the Green function approach, we can introduce different channels (c) characterized, e.g., by spin and isospin quantum numbers. We assume that these channels decouple. In contrast to the angular momentum which is not conserved, e.g., for tensor forces, the contribution of different channels to Equation (13) is additive. The remaining intrinsic quantum numbers will be denoted by ν c , it concerns the bound states as far as they exist (ground states and excited states), as well as the continuum of scattering states.
We analyze the contributions of the clusters (A ≥ 2), suppressing the thermodynamic variables {T, µ n , µ p }. We have to perform the integral over the c.o.m. momentum P what, in general, must be done numerically since the dependence of the in-medium quasiparticle energies E A,ν (P; T, µ n , µ p ) on P is complicated. We have in the non-degenerate case (15) with g A,c = 2s A,c + 1 the degeneration factor in the channel c. The partial density of the channel c at P contains the intrinsic partition function It can be decomposed in the bound state contribution and the contribution of scattering states z cont A,c (P; T, µ n , µ p ). We emphasize that the subdivision of the intrinsic partition function into a bound and a scattering contribution is artificial and not of physical relevance. The summations over A, c and P of (16) remain to be done for the EOS (13), and Z may be included in c. The region in the parameter space, in particular P, where bound states exist, may be restricted what is expressed by the step function Θ(x) = 1, x ≥ 0; = 0 else. The continuum edge of scattering states is denoted by E cont A,c (P; T, µ n , µ p ). The intrinsic partition function z cont A,c (P) contains the scattering state contribution (non degenerate case) Going beyond the ordinary Beth-Uhlenbeck formula [15], for nuclear matter the generalized Beth-Uhlenbeck formula has been worked out in Ref. [12]. Here, the single-particle contribution is described by quasiparticles, and to avoid double counting, the corresponding mean-field term must be subtracted from the contribution of scattering states what leads to the appearance of the term sin 2 δ c (E).
The NSE follows as a simple approximation where the sum is performed only over the bound states, and medium effects are neglected. We obtain the model of a mixture of non-interacting bound clusters, which can react so that chemical equilibrium is established. This approximation may be applicable for the low-density, low-temperature region where the components are nearly freely moving, and intrinsic excitations are not of relevance. However, the continuum of scattering states (the intrinsic quantum number ν c may contain the relative momentum) are of relevance at increasing temperature, which is clearly seen in the exact Beth-Uhlenbeck expression [15] for the second virial coefficient. Another argument to take the continuum of scattering states into account is the dissolution of bound states at increasing density. Here, the thermodynamic properties behave smoothly because the contribution of the bound state to the intrinsic partition function is replaced by the corresponding contribution of the scattering states.
It turns out advantageous for analyses of the thermodynamics of the Mott transition, to avoid the separation into a bound and a scattering state part of the spectrum and to include the discrete part of the spectrum into the definition of the phase shifts, so that these are merely parameters of a polar representation of the complex A−particle cluster propagator. The partition function then takes the form As an example let us consider the case A = 2. A consistent description of the medium effects should contain not only the mean-field shift of the quasiparticle energies but also the Pauli blocking, as shown in this work for conserving approximations. Within the generalized Beth-Uhlenbeck approach [12], Pauli blocking modifies the binding energy of the bound state (deuteron) in the isospin-singlet channel as well as the scattering phase shifts. The integrand of the intrinsic partition function is shown in Figure 3.
The intrinsic partition function z cont A,c (P) cannot be divided unambiguously into a bound state contribution and a contribution of scattering states. As seen from Figure 3, the sum of both contributions behaves smoothly when with increasing density the bound state is dissolved. Therefore, we emphasize that one should consider only the total intrinsic partition function including both, bound and scattering contribution, as expressed by the generalized phase shifts as shown in Figure 3.
Whereas the two-body problem can be solved, e.g., using separable interaction potentials, and the account of medium effects has been investigated [12], the evaluation of the contribution of clusters with mass numbers A > 2 to the EoS is challenging. The medium modification of the cluster binding energies has been calculated, e.g., using a variational approach [17]. Problematic is the inclusion of scattering states, in particular the treatment of different channels describing the decomposition of the A-particle cluster. In this work, we discuss the concept of a cluster-virial expansion and corresponding generalized cluster-Beth Uhlenbeck approaches. This concept is based on the chemical picture where bound states are considered as new components and can be treated in the same way as "elementary" particles. It is generally accepted that a second virial coefficient can be introduced for systems consisting of atoms, but also with molecules as components. The problem is the introduction of an effective interaction between the components (including the quantum symmetry postulates for fermions or bosons), and intrinsic excitations of the bound components are described in some approximations. Within the NSE, we describe a nuclear system as a mixture of "free" nucleons in single-particle states as well as of clusters (nuclei with A > 1). Taking into account the interaction between these components of a nuclear system, from measured scattering phase shifts (for instance α − α scattering) a virial EoS can be derived. Results have been presented in Ref. [19].
On a more microscopic level, we consider here interacting quarks which can form bound states (hadrons), and the general approach should include both cases, the region of the quark-gluon plasma and, after the confinement transition, the region of well established hadrons. A main difficulty is the introduction of an effective interaction which can be made by fitting empirical data. However, a systematic quantum statistical approach is needed to derive such effective interaction from a fundamental Lagrangian, and to introduce the cluster states performing consistent approximations and avoiding double counting of the contributions to the EoS and other physical properties.
The Green function technique as well as the path-integral approach are such systematic quantum statistical approaches. Different contributions to the EoS are represented by Feynman diagrams, and double counting is clearly excluded. A selection of diagram classes can be performed to recover the chemical picture. As seen from Equation (11), after separation of the center-of-mass momentum P, the propagator for the A-particle bound states is where ν covers only the bound state part of the internal quantum state of the A-particle cluster. As a new element, the bound state propagator is introduced as indicated in Figure 4. This bound state propagator has the same analytical form like the single particle propagator (12), besides the appearance of the internal wave function that determines the vertex function.

G G
A A A, sc = + Figure 4. Splitting of the A-particle cluster propagator into a bound and scattering contribution. Note that the internal quantum number has been dropped.
As discussed above, different approximations are obtained such as the nuclear statistical equilibrium (NSE) and the cluster mean-field (CMF) approximation using the chemical picture. These approximations are based on the bound state part of the A-particle propagator. They give leading contributions in the low-density, low temperature range where bound states dominate the composition of the many-particle system.
To improve the approximation, the remaining part of continuum correlations has also been taken into account. From the point of view of the physical picture, these contributions arise in higher orders of the virial expansion of the equation of state. As example, the formation of the A-particle bound state is seen in the A-th virial coefficient, the mean-field shift due to a cluster B in the (A + B)-th virial coefficient. The chemical picture indicates which high-order virial coefficients of the virial expansion are essential, if the many-particle system is strongly correlated so that bound states are formed.
In an improved approximation, the scattering part of the A-particle propagator has to be considered. It contributes also to the A-th virial coefficient. The scattering processes within the A-particle system can have different channels. As an example we discuss here binary elastic scattering processes between sub-clusters A 1 and A 2 of the system of A particles, A = A 1 + A 2 . Binary phase shifts δ A 1 ,A 2 (E) are introduced that describe the corresponding scattering experiments. They can also be calculated within few-body theory. Besides the effective interaction between the sub-clusters that are depending on the internal wave function of the sub-clusters, also virtual transitions to excited states may be taken into account. In general, the effective interaction is non-local in space and time, i.e., momentum and frequency dependent.
A generalized cluster Beth-Uhlenbeck formula, see Equation (16), is obtained when in particle loops not the free propagator, but quasiparticle Green's functions are used. If the quasiparticle shift is calculated in Hartree-Fock approximation, the first order term of the interaction must be excluded from the ladder T ladder 2 matrix to avoid double counting. The bound state part is not affected, it is determined by an infinite number of diagrams. The scattering part is reduced subtracting the Born contribution as shown in Equation (18) by the 2 [sin(δ c )] 2 term; for the derivation see Ref. [12].
The continuum correlations that are not considered in the NSE give a contribution to the second virial coefficient in the chemical picture. We can extract from the continuum part two contributions [17]: resonances that can be treated like new particles in the law of mass action, and the quasiparticle shift of the different components contributing to the law of mass action. Both processes are expected to represent significant contributions of the continuum. After projecting out these effects, the residual contribution of the two-nucleon continuum is assumed to be reduced. One can try to parametrize the residual part, using the ambiguity in defining the bound state contribution. Eventually the residual part of the continuum correlations can be neglected.
The correct formulation of the cluster-virial expansion [16] is considered to be a main ingredient towards a unified EoS describing quark matter as well as nuclear matter. One has to introduce the interaction between the constituents in a systematic way. The account of correlations in the continuum is essential near the confinement phase transition where the hadronic bound states disappear. Conserving approximations may lead to acceptable results for the EoS. An important issue is the account of correlation in the medium, in particular when considering the Pauli blocking in the phase space of the elementary constituents.

Φ-Derivable Approach to the Cluster Virial Expansion for Nuclear Matter
Recently, it has been suggested [20] that the cluster virial expansion for many-particle systems [16] can be formulated within the Φ-derivable approach [3,4]. This approach is straightforwardly generalized to A-particle correlations in a many-fermion system where the full A−particle Green's function obeys the Dyson equation where G A is the free A−particle Green's function and the selfenergy is defined as a functional derivative of the two-cluster irreducible Φ functional This generalization of the Φ−derivable approach fulfills by its construction the conditions of stationarity of the thermodynamical potential with respect to variations of the cluster Green's functions The Φ functional for our purpose of defining a cluster decomposition of the system with inclusion of residual interactions among the clusters captured by a second virial coefficient is given by a sum of all two-loop diagrams that can be drawn with cluster Green's functions G i+j and their subcluster Green's functions G i and G j for a given bipartition with the appropriate vertex functions Γ i+j;ij . This generalization of the so-called "sunset" diagram case is depicted diagrammatically in Figure 5. Herewith we have generalized the notion of the Φ−derivable approach to that of a system where the hierarchy of higher order Green functions is built successively from the tower of all Greens functions starting with the fundamental one G 1 . The open question is how to define the vertex functions joining the cluster Greens functions. As a heuristic first step one may always introduce local coupling constants, like in the Lee model discussed in the context of a Φ−derivable approach by Weinhold et al. [6]. Introducing nonlocal formfactors at the vertices will correspond to a separable representation of the interaction. The definition of the interaction can be absorbed in the introduction of the cluster T-matrix in the corresponding channel which, in the ladder approximation, will reproduce the analytical properties encoded in the full cluster Green's function concerning bound state poles and scattering state continuum. This equivalence is depicted in Figure 6. The T A matrix fulfills the Bethe-Salpeter equation in ladder approximation which in the separable approximation for the interaction potential, leads to the closed expression for the T A matrix where the generalized polarization function has been introduced and the one-frequency free i−particle Green's function is defined by the (i − 1)-fold Matsubara sum Note that for these Green's functions holds the relationship ( Another set of useful relationships follows from the fact that in the ladder approximation both, the full two-cluster (i + j particle) T matrix (26) and the corresponding Greens' function have similar analytic properties determined by the i + j cluster polarization loop integral (27) and are related by the identity which is straightforwardly proven by multiplying Equation (26) with G i+j and using Equation (30). Since these two equivalent expressions in Equation (31) are at the same time equivalent to the two-cluster irreducible Φ functional introduced above in Equations (20) and Figure 5, the functional relations follow and may become useful in proving cancellations that are essential for the relationship to the Generalized Beth-Uhlenbeck approach as discussed below.

Generalized Beth-Uhlenbeck EoS from the Φ−Derivable Approach
Now we return to the question how the relation between the cluster Φ−derivable approach to the partition Function (20) and the generalized Beth-Uhlenbeck equation for the cluster density may be established. To this end we consider the partial density of the A−particle state defined as Taking into account that any analytic complex function F(ω) has the spectral representation we perform the Matsubara summation 3 in Equation (20) ∑ Using the relation ∂ f A (ω)/∂µ = −∂ f A (ω)/∂ω we get for Equation (34) now where a partial integration over ω has been performed and the degeneracy factor d A for cluster state has been introduced, stemming from the trace operation in the internal spaces. Now we use the fact that for two-loop diagrams of the sunset type a cancellation holds [21,22] which we generalize here for cluster states Using generalized optical theorems [8,12] we can show that where the phase shifts δ A have been introduced via the polar representation of the complex A−particle propagator G A = |G A | exp(iδ A ). With these ingredients follows from the cluster Φ−derivable approach the cluster virial expansion for the density in the form of a generalized Beth-Uhlenbeck EoS 3 For odd A, z n = (2n + 1)πT + µ are the fermionic Matsubara frequencies and for even A, z n = 2nπT + µ the bosonic ones.
In this way we have drawn the connection between the cluster virial expansion of Ref. [16] reviewed in the previous section with the Φ−derivable approach [3,4]. In the following subsection we consider the example of deuterons in nuclear matter in order to elucidate the application of the approach.

Deuterons in Nuclear Matter
Within the Φ−derivable approach [3,4] the grand canonical thermodynamic potential for a dense fermion system with two-particle correlations is given as where the full propagators obey the Dyson-Schwinger equations with selfenergies which are defined by the choice for the Φ functional, a two-particle irreducible set of diagrams such the ones in Figure 7.  Figure 5 with the two-cluster potential V i+j .
The functional for the thermodynamic potential (41) is constructed such that the requirement of its stationarity, in thermodynamic equilibrium is equivalent to Equation (4) where S 1 (1, ω) = 2 G 1 (1, ω + iη) is the fermion spectral function and Equation (45) expresses particle number conservation in a system with volume V.
Having introduced the notion of a cluster expansion of the Φ functional we want to suggest a definition which eliminates the unknown vertex functions in favour of the T A+B matrix which describes the nonperturbative binary collisions of A− and B− particle correlations in the channel A + B, see Figure 6. The application of this scheme to the simplest case of two-particle correlations in the deuteron channel in nuclear matter results in the selfenergy [12] where the correlation density contains besides the bound state a scattering state contribution as can be seen from examining the derivative of the phase shift shown in Figure 3. The one-particle density of free quasiparticle nucleons n q u is reduced in order to fulfil the baryon number conservation in the presence of deuteron correlations and contains a selfenergy contribution due to the deuteron correlations in the medium. This improvement of the quasiparticle picture due to the correlated medium accounted for by the consistent definition of the selfenergy as a derivative of the Φ Functional (20) is the reason the continuum correlations (48) are reduced by the factor 2 sin 2 δ as compared to the traditional Beth-Uhlenbeck formula [15,23]. For details, see [8,12]. With the definition of the Φ functional via the T matrix in Figure 6 we were able to show the correspondence between the generalized Beth-Uhlenbeck approach and the Φ−derivable approach for the nonrelativistic potential model approach to two-particle correlations in a warm, dense Fermion system [12,16]. Now we would like to discuss its application to a relativistic model for correlations in quark matter: mesons, diquarks and baryons.

Cluster Virial Expansion for Quark-Hadron Matter within the Φ Derivable Approach
Finally, we would like to sketch how the Φ derivable approach can be employed to define a cluster virial expansion for quark-hadron matter consisting of quarks (Q), mesons (M), diquarks (D) and baryons (B) that can represent a unified quark-hadron matter EoS. The thermodynamical potential for this system obtains the form very similar to the case of clustered nuclear matter, i.e.
where c i = 1/2 (c i = −1/2) for bosonic (fermionic) states and d i are the degeneracy factors that stem from the trace operation in the internal spaces of the quark, meson, diquark and baryon states. We suggest that in going from Equation (49) to (50) the same cancellations apply that were used above for the density formula and that are known to apply also for the entropy [21,22] would allow to derive this generalized Beth-Uhlenbeck equation of state for the thermodynamic potential, i.e., the negative pressure, once we restrict ourselves to the minimal set of two-particle irreducible diagrams in defining the Φ functional by the class of sunset type diagrams only, as given in Figure 8. From this Φ functional follow the selfenergies defining the full Greens functions of the system by functional derivation The resulting Feynman diagrams for the selfenergy contributions are given in Figure 9. Note that it is immediately plain from this formulation that in the situation of confinement, when the propagators belonging to colored excitations (quarks and diquarks) and thus to states that could not be populated would be cancelled, the system simplifies considerably. When all closed loop diagrams containing quarks and diquarks are neglected, this system reduces to a meson-baryon system. Out of the five closed-loop diagrams of Figure 8 remains then only the rightmost one from which the two selfenergy diagrams in Figure 9 emerge that contain only meson and baryon lines.

Relativistic Density Functional Approach to Nuclear Matter
In the limit of quark (and gluon) confinement, the meson-baryon system can be further reduced when the mesonic degrees of freedom are not considered as dynamic ones but just as their meanfield values coupled to the baryon degrees of freedom with effective, possibly density-and temperature-dependent couplings, as sketched in Figure 10. U = Figure 10. Effective density-functional for the interaction of species i and j given by a density-dependent coupling G i+j as the local limit of a Φ functional.
In this case, the Φ−derivable approach reduces to a selfconsistent relativistic meanfield theory of nuclear matter, where the functional U n i,j stands for the interaction in the system and is expressed in terms of the couplings of bilinears of the baryon spinors that define their scalar and vector densities n i,S and n i,V , resp. The thermodynamic potential (52), where the role of the Φ functional is now played by the U functional that depends not on the propagators but on the densities and is therefore called a density functional. The thermodynamic potential (52) shares with the Φ−derivable approach the fulfillment of the conditions of stationarity and selfconsistency, which is provided by the fact that the (now real) selfenergies in the Dyson equations are obtained as derivatives of the U functional The baryon quasiparticle propagators S i,qu fulfill the Dyson equations To this class of relativistic density functional models for baryonic matter belong the density-dependent relativistic meanfield models known as, e.g., DD2, NL3, KVOR, TM1. Their nonrelativistic relatives are the density functional models of the Skyrme type.
Recently, these models have been augmented with the inclusion of hadronic excluded volume effects. For an elaborate version of such corrections, see [24]. The origin of the excluded volume effects is the quark substructure of the baryons which entails the quark exchange effects between baryons that are a consequence of the Pauli principle on the quark level of description.

Quark Pauli Blocking in Hadronic Matter
The step to cancel all diagrams that contain propagators of colored excitations with the argument of confinement may be a too drastic step when we want to describe matter with a high density so that phase space occupation effects shall become important and the hadronic bound states "remember" their quark substructure. Formally these effects are included into the selfconsistent description of quark-hadron matter within the Φ−derivable approach outlined above in Equation (50) with the Φ functional given by the diagrams in Figure 8. In the limit of the confined phase, however, it is important to restore these quark substructure effects as they will drive the system towards deconfinement for sufficiently large densities. This can be accomplished by including selfenergy diagrams containing quark and diquark lines in Figure 9 in a perturbative manner by considering their propagators not fully selfconsistently, but in first order with respect to their selfenergies.
Here the quark and diquark quasiparticle propagator lines in Equations (56) and (57) are defined by a Dyson-Schwinger equation as given generically in Equation (55), but with a density functional for the effective interaction that is appropriate for the decription os quark matter and will be discussed in detail in Section 4.4.
In such a way two quark-diquark substructure contributions to the baryon selfenergy appear due to the corrections (56) and (57) beyond the quasiparticle approximation for quark and diquark propagators = Ω(Σ (0) ) + + + O(Σ (2) ) . (58) They contain one closed baryon line and are therefore of first order in the baryon density. By functional derivative w.r.t. the baryon propagator line (cutting) an effective quark and diquark exchange interaction can be obtained from those contributions to the baryon selfenergy shown in Equation (58). The diagram for the quark exchange interaction between baryons resulting from cutting the baryon line in the first of the two diagrams in Equation (58) is shown in (59) in two forms which are topologically equivalent, Analogously, the diquark exchange interaction is obtained by cutting the baryon line in the second of the two diagrams in Equation (58).
The quark Pauli blocking effects in nuclear matter have been evaluated in a nonrelativistic approximation with a confining potential model in Refs. [25,26] where it was found that the result for the repulsive density-dependent nucleon-nucleon interaction corresponds well to the repulsive part of the effective Skyrme interaction functional in Ref. [27]. Note that the resulting EoS has been successfully applied in predicting massive hybrid stars with quark matter cores [28]. A flaw of the nonrelativistic calculations of the quark Pauli blocking effect is that the quark mass is a fixed parameter so that partial chiral symmetry restoration in dense hadronic matter as a selfenergy effect on the quark propagator (consistent with the quark exchange) is not accounted for. This question has recently been taken up by Blaschke, Grigorian and Röpke who demonstrated that a chirally improved calculation of the quark Pauli blocking effect results in an EoS for nuclear matter that is similar to the DD2 model with excluded volume [24].
Strange quarks as well as strange hadrons belong to the system of our model, as it is formulated for any flavour. Interesting new aspects, which are expected from this approach concern, for instance, the Pauli blocking between baryons, including hyperons, in dense matter. This shall be of relevance for the discussion of the hyperon puzzle in compact stars [29,30].
The quark Pauli blocking effect applies also to mesons and corresponding expressions for selfenergy effects can be extracted from the Φ−derivable approach in a similar manner as for the baryons. This has been outlined in Ref. [31]. In a nonrelativistic potential model calculation, an effective quark exchange potential for the π-π interaction has been derived [32] which reproduces the scattering length of the pion interaction in the isospin = 2 channel, see also [33].
Let us now turn to the other limit of the Φ−derivable approach to a unified EoS for quark-hadron matter, the case of deconfined quarks. In this case, also chiral symmetry is restored, and due to the resulting lowering of the mass threshold the meson and baryon states become unbound (Mott effect). Their contribution to the thermodynamics as captured in the corresponding phase shift functions is gradually vanishing at high temperatures and chemical potentials with just chiral quark matter remaining asymptotically. As a paradigmatic example for the treatment of Mott dissociation of hadrons in hot, dense matter let us consider the case of pion dissociation in hot quark matter.

Mott Dissociation of Pions in Quark Matter
In order to describe the problem of mesons in quark matter within the Φ−derivable approach we define the Φ functional and the corresponding selfenergy in Figure 11. The meson polarization loop Π M (q, z) in the middle panel of Figure 11 enters the definition of the meson T matrix (often called propagator) which in the polar representation introduces a phase shift δ M (q, ω) = arctan( T M / T M ), that results in a generalized Beth-Uhlenbeck equation of state for the thermodynamics of the consistently coupled quark-meson system [34].
where the selfconsistent quark meanfield contribution is with the quasiparticle energy shift for quarks (antiquarks) due to mesonic correlations given by Σ ± = ∑ M=π,σ tr D [Σ M Λ ± γ 0 ] /2 and the positive (negative) energy projection operators Λ ± = (1 ± γ 0 )/2. The mesonic contribution to the thermodynamics is where similar to the case of deuterons in nuclear matter the factor 2 sin 2 δ M accounts for the fact that mesonic correlations in the continuum are partly already accounted for by the selfenergies Σ M defining the improved selfconsistent quasiparticle picture. In the previous works of Refs. [34][35][36][37][38] on this topic, however, the effect of the backreaction from mesonic correlations on the quark meanfield thermodynamics had been disregarded. In Figure 12 we show the phase shift of the pion as a quark-antiquark state for different temperatures, below and above the Mott dissociation temperature. The shape of these functions and their evolution with increasing temperature over the Mott dissociation resembles the similar behaviour of the deuteron phase shift in nuclear matter at increasing density, see Figure 3. Here we note from the Φ−derivable approach that for consistency the quark propagator in the quark meanfield thermodynamic potential shall contain effects from the selfenergy Σ M due to the coupling to the mesonic correlations as in the right panel of Figure 11. This total quark selfenergy is the given by Σ(p, p 0 ) = σ MF + Σ M (p, p 0 ), where for a local NJL model with scalar coupling constant G S the meanfield contribution is and the contribution due to scalar/pseudoscalar mesons (corresponding to the diagram shown in the rightmost panel of Figure 11) is given by [40] where M = (−1/π) T M (q, ω + iη) is the meson spectral density and E q = q 2 + m 2 is the quark dispersion law with the quark mass m = m 0 + σ MF . One can observe the similarity of this result (65) with that for a Dirac fermion coupled to a pointlike scalar meson, as given in [41]. Finally, let us consider the aspect of quark deconfinement in dense matter within a relativistic density functional approach which can be considered as a local limit of the Φ−derivable approach to quark-hadron matter (50) that obtains a form similar to the relativistic density functional theory for hadronic matter. The important difference lies in the density functional that in the case of quark matter shall account for confining effects, see [25,42]. In the final part of this section, we outline the recently developed relativistic version of the so-called string-flip model that proved rather successful for phenomenological applications to compact stars, supernovae and neutron star mergers.

Relativistic Density Functional Approach to Quark Matter
Hence, one obtains a contribution to the energy density functional of quark matter correspondingly [43]. In analogy to the Walecka model of nuclear matter [44], the relativistic density-functional approach to interacting quark matter can be obtained from the path integral approach based on the partition function [42], with q = (q u , q d ) T ,μ = diag(µ u , µ d ) and effective Lagrangian density L eff = L free − U(qq,qγ 0 q). The interaction is given by the potential U(qq,qγ 0 q), which is a nonlinear functional of the scalar and vector quark field-currents. In the mean-field approximation, this potential can be expanded around the expectation values of the field currents, n s = qq and n v = qγ 0 q respectively, with scalar and vector self-energies, Σ s and Σ v . By appropriate rearranging of the quantities and performing the path integrals of Equation (66) one gets the thermodynamic potential The quasi-particle term (for the case of isospin symmetry and degenerate flavors) can be calculated by using the ideal Fermi gas distribution for quarks with the quasiparticle energy E * = p 2 + M 2 , the effective mass M = m + Σ s and effective chemical potential µ * = µ − Σ v . The self energies are determined by the density derivations In this approach the stationarity of the thermodynamical potential is always fulfilled. For the case of isospin asymmetry, see Ref. [42].
In the mean-field approximation, the correlation energy can be obtained by folding the string-length distribution function for a given density with some interaction potential [25,45,46]. Moreover, the average string length between quarks in uniform matter is related to the scalar density, n s , being proportional to n −1/3 s . To capture this phenomenology, the following density functional of the interaction is adopted, The first term captures aspects of (quark) confinement through the density dependent scalar self-energy, defining the effective quark mass M = m + Σ s . We also employ higher-order quark interactions [47], by inclusion of the third term in Equation (73), for the description of hybrid stars (neutron stars with a quark matter core) in order to obey the observational constraint of 2 M . To this end, the denominator in the last term of Equation (73) guarantees that for the appropriate choice of the parameters b and c, causality is not violated (i.e., the speed of sound c s = √ ∂P/∂ε does not exceed the speed of light). All terms in Equation (73) that contain the vector density contribute to the shift defining the effective chemical potentials µ * = µ − Σ V , where The SFM modification takes into account the effective reduction of the in-medium string tension, D(n v ) = D 0 φ(n v ; α). It is understood as a consequence of the modification of the pressure on the color field lines by the dual Meissner effect, since the reduction of the available volume corresponds to the reduction of the non-perturbative dual superconductor QCD vacuum that determines the strength of the confining potential between the quarks. The reduction of the string tension is modeled via a Gaussian function of the baryon density n v , similar the available volume fraction in dense nuclear matter (see [24] and Section 5.3 below), as it shall be related to the available volume of nonperturbative QCD vacuum that according to the dual superconductor model is responsible for the formation of stringy color field configurations because of the dual Meissner effect. A detailed discussion of the role of the parameters a, b, c and α is given in Ref. [42].

Applications in Core-Collapse Supernovae and Neutron Stars
In the recent years the equation of state (EoS) has been significantly constraint, in particular at densities around and in excess of nuclear saturation density (ρ 0 ). At densities ρ ≤ ρ 0 , chiral effective field theory is the ab-initio approach to the nuclear many-body problem of dilute neutron matter [48][49][50][51][52][53][54]. Moreover, the high-precision observation of massive neutron stars of about 2 M [55][56][57] constraints the supersaturation-density EoS, i.e., sufficient stiffness is required. The latter aspect challenges the appearance of additional particle degrees of freedom, e.g., hyperons and quarks, which tend to soften the EoS at ρ > ρ 0 .
While neutron stars feature matter in β-equilibrium at zero temperature, the challenge lays in the development of EoS for core-collapse supernova applications, which cover a large domain of temperature, density and isospin asymmetry (cf. Figure 1a of Ref. [58]). At T ≤ 0.5 MeV, time-dependent processes determine the nuclear composition, where heavy nuclei dominate. With increasing temperature, towards T 0.5 MeV, complete chemical and thermal equilibrium known as NSE (nuclear statistical equilibrium) is achieved, where, the nuclear composition is determined from the three independent variables: T, ρ (or n B ), y C (baryonic charge fraction). Note that there is a strong density dependence, i.e., these heavy nuclei, originally belonging to the iron group with A 56, become increasingly heavier with increasing density. This phenomenon is well known also from the neutron star crust, where due to the low temperatures only a single nucleus exists for at given density, instead of a (broad) distribution as in the supernova case. At ρ = ρ 0 and at temperatures above T = 5 − 10 MeV, nuclei dissolve at the liquid-gas phase transition into homogeneous nuclear matter composed of (quasi-free) nucleons [16,59,60].
It becomes evident that first-principle calculations covering the entire domain are presently inexistent. Instead, model EOS are being developed for astrophysical applications. These combine several domains with different degrees of freedom, i.e., heavy nuclei at low temperatures, inhomogeneous nuclear matter with light and heavy nuclei together with the free nucleons (mean field), and homogeneous matter at high temperatures and densities. The latter has long been subject to investigations of a possible phase transition to the quark-gluon plasma.
The role of the EOS in simulations of core-collapse supernovae was explored within the failed scenario and consequently the formation of black holes, with focus on the dynamics and the neutrino signal [61][62][63][64]. In the multi-dimensional framework, neutrino-driven supernova explosions were the subjects of investigation [65][66][67], where it was found that such explosions are favored for soft EOS [68] with an earlier onset of shock revival and generally higher explosion energies, in comparison to stiff EOS [69]. Moreover, the role of the nuclear symmetry energy has been reviewed [70]. This an important nuclear matter parameter has a strong density dependence and becomes more tightly constrained by experiments, nuclear theory and observations [71,72].

Heavy Nuclear Clusters with A 56
At low temperatures (T < 0.5 MeV) time-dependent nuclear processes determine the evolution, which corresponds to the outer core of the stellar progenitor, with the nuclear composition of dominantly silicon, sulfur as well as carbon and oxygen. In some cases, even parts of the hydrogen-rich helium envelope are taken into account, e.g., in simulations of supernova explosions following the shock evolution for tens of seconds through parts of the stellar envelope. Therefore, small nuclear reaction networks are commonly employed in supernova studies [73,74]. They are sufficient for the nuclear energy generation.
At T 0.5 MeV, NSE is fulfilled and the relation µ (A,Z) = Zµ p + (A − Z)µ n between the chemical potential of nuclei µ (A,Z) , with atomic mass A and charge Z, and the chemical potentials of neutron µ n and proton µ p holds. The NSE conditions found in the collapsing stellar core feature a broad distribution of nuclei with a pronounced peak around the iron-group, at low densities (see Figure 2 of Ref. [58]). In simulations of supernovae, this nuclear distribution is classified by the NSE average, including nuclear shell effects as discussed [75], which extends beyond the commonly used single-nucleus approximation. This is important for the consideration of weak processes with heavy nuclei. In particular, rates for electron captures on protons bound in these heavy nuclei [76] are averaged over the NSE composition and provided to the community as a table. In addition, coherent neutrino-nucleus scattering is considered [77], where even inelastic contributions are take into account [78], as well as nuclear (de)excitations [79,80]. The profound understanding of weak processes associated with nuclear transitions is also important for the understanding of the physics of the neutron star crust (cf. Ref. [81] and references therein), in particular for accreting neutron stars leading to the phenomenon of deep-crustal heating [82][83][84].

Light Nuclear Clusters with A ≤ 4
In the domain corresponding to NSE, there is a rather narrow density domain where light nuclear clusters such as 2 H, 3 H, 3 He and 4 He can exist, at finite temperatures on the order of few MeV [59]. There are two aspects related to the presence of these light nuclear clusters, modification of the nuclear EOS due to these degrees of freedom, and the neutrino response due to the inclusion of a large variety of weak processes (cf. Table 1 in Ref. [85]). Common model supernova EOS with (light) clusters are based on the modified NSE [75]. However, the dissolving of the clusters towards high density is mimicked by the geometric excluded volume approach, as well as by hand with increasing temperature. This applies equally for light and heavy nuclear clusters. In comparison with the state-of-the-art quantum statistical approach [86,87] the deficits of NSE are revealed [58], while the cluster-virial EOS can provide the constraint at low densities [16]. It relates to the overestimation of the abundance of light clusters and in particular the too late dissolving of clusters into homogeneous matter.
Besides the NSE approach with 'all' (light) clusters included, in simulations of core-collapse supernovae the simplified nuclear composition (n, p, α, A, Z ), with only 4 He as representative light cluster, has long been employed [68,69]. It leads to an overestimate of the abundances of the unbound baryons and 4 He. This has important consequences for the supernova results, since such simplistic approach overestimates the neutrino response with the neutrons and protons, which in turn changes the supernova neutrino fluxes and spectra, however, with only a mild impact on the overall supernova dynamics [58].
In supernova simulations the temperatures are generally too high for any significant abundance of any light cluster. Only when the supernova explosion onset has been launched and the remnant central proto-neutron star deleptonizes and cools via the emission of neutrinos, light clusters can start to play a role. However, it has been demonstrated that neutral-current reactions are dominated by scattering on free neutrons, which is the most abundant nuclear species due to the generally large neutron excess (>95%) of supernova matter. Scattering reactions with light clusters play a negligible role. On the other hand, the neutrino response for light clusters is dominated by charged-current (break-up) reactions involving deuteron, triton and helium-3. However, the ν e charged-current opacity is dominated by absorption on free neutrons. The situation is different forν e , since it was shown that protons and light clusters have similar abundances [58,85]. Taking properly into medium modifications for the cross sections and the proper phase-space of the contributing particles, the final rates are generally small. They never reach values as for the standard rates with protons, and hence the impact from charged-current weak processes with light clusters was found to be negligible [58,85]. Note that this study was based on the NSE approach which generally overestimates the abundance of light clusters. Hence, any improved treatment of the nuclear composition will most likely even reduce the impact of light clusters and the associated weak processes.

Homogeneous Matter at Supersaturation Density and Phase Transition to Quark Matter
With increasing densities, around ρ 0 (depending on the temperature), the transition to homogeneous nuclear matter proceeds where the EOS becomes less and less constrained by nuclear physics. As discussed above, the quark substructure of baryons shall become apparent at supersaturation densities and manifest itself by a nucleonic hard core repulsion due to quark Pauli blocking. This effect is likely to be enhanced by chiral symmetry restoration at high densities. To explore its role for the nuclear EOS at supersaturation densities the geometric excluded volume mechanism can be employed [24], where the available volume of the nucleons, V N = V φ(ρ; v), is proportional to the total volume V of the system where as proportionality factor occurs the available volume fraction This is a density functional taken here according to [24] in a Gaussian form similar to (76). Depending on the sign of the excluded volume parameter, v, it allows to model both, stiffening and softening of the supersaturation-density EOS based on some reference model. This approach has been applied to confirm the independence of supernova simulations, e.g., the supernova shock dynamics as well as the evolution of neutrino luminosities and average energies, on the supersaturation-density EOS [88]. In studies of neutron stars, this approach results in significantly altered neutron star properties, e.g., maximum masses and radii [42,47]. The latter property is currently constraint only poorly, from the analysis of observations of low-mass x-ray binary systems [89][90][91].
Another uncertain aspect of the supersaturation-density EoS is the possibility of a phase transition from nuclear matter, with hadrons as degrees of freedom, to the deconfined quark gluon plasma with quarks and gluons as the new degrees of freedom. This has long been explored in the context of cold neutron stars. Unfortunately, the evaluation of the partition function of Quantum Chromodynamics (QCD)-the theory of strongly interacting matter-is possible only at vanishing baryon density by means of large-scale numerical simulations of this gauge field theory in a representation on space-time lattices [92,93]. These numerical ab-initio solutions predict a smooth cross-over transition at a temperature of T = 154 ± 9 MeV at µ B 0, see Refs. [94][95][96]. Consequently, to study the role of quark degrees of freedom at high baryon densities, effective models for low-energy QCD have to be employed. Generally, the nuclear and quark matter phases are modeled separately and a phase transition construction is employed. This so-called two-phase approach results in a first-order phase transition by design. Note further that perturbative QCD, which is valid in the limit of asymptotic freedom, where the smallness of the coupling between quarks allows the usage of perturbative methods and corrections to the behaviour of an ideal gas of ultra-relativistic particles [97] are small, become applicable only at extremely high temperatures and densities exceeding by far the values attainable in compact stars or ultrarelativistic heavy-ion collisions. Instead, for studies of neutron stars and supernovae, effective quark matter models have been commonly employed, such as the thermodynamic bag model [98], models based on the Nambu-Jona-Lasinio (NJL) type [99][100][101], the recently developed vector-interaction enhanced chiral bag model [102,103] and in particular the density-functional-based DD2-SFM hybrid EoS approach [42] that has been outlined above. With the finite-temperature extension of the DD2F-SFM EoS it was possible to show that the deconfinement phase transition may serve as an explosion mechanism for massive (∼50 M ) blue-supergiant stars [104] long sought-for. At the same time, it explains the occurrence of a population of neutron stars born with high masses.
It has been realized that repulsive interactions are the necessary ingredient to provide a sufficient stiffness for the EOS at high densities in order to yield massive neutron stars with quark-matter core (known as hybrid stars) in agreement with the current constraint of 2 M . Moreover, higher-order vector repulsion terms [42,47] can lead to the 'twin' phenomenon (cf. [105][106][107][108] and references therein). It relates to the existence of compact stellar objects with similar-to-equal masses but different radii, to a strong phase transition with a large latent heat. As a consequence, the hybrid stars for such an EoS appear in the mass-radius diagram on a disconnected "third family" branch, separated from the branch of neutron stars (second family) by a sequence of unstable configurations. As has been demonstrated impressively for the DD2-SFM hybrid EoS [42], by varying only the available volume parameter α describing the screening of the confining interaction in dense matter, the twin phenomenon can be obtained at high masses of ∼2 M as well as for typical pulsar masses of 1.3-1.4 M or even below. This feature allows to discuss such hybrid stars in the context of the first multi-messenger observation of the binary compact star merger GW170817 [109], where such a scenario appears as an alternative to the conservative one of a binary neutron star merger with a relatively soft EoS [110][111][112]. In this interesting situation the NICER (Neutron Star Interior Composition Explorer) 4 NASA mission has the potential to rule out the soft EoS scenario, when it would measure for its primary target, the nearest millisecond pulsar PSR J0437-4715 with a mass of 1.44 ± 0.07 M a large radius of, say, 14 km with the expected high precision of 0.5 km. Such a measurement would contradict the constraint on compactness of neutron stars extracted from GW170817 in the double neutron star merger scenario that constrains the radius in the mass range of 1.4 M to R < 13.4 km, see [113]. Such a measurement would imply the discovery of the third family of compact stars [112] with an onset mass in the mass range 1.16-1.60 M extracted from the gravitational wave signal of binary inspiral [109].

Light Cluster Formation and Symmetry Energy in Low-Energy Heavy-Ion Collisions
The only possibility to probe the properties of hot and dense nuclear matter in the laboratory are heavy ion collisions (HIC). The fragment distributions, their energy spectra and correlations measured in the detectors are used to infer the properties of the initial state ("fireball") produced in HIC. To reconstruct the initial state, one has to model the time evolution of the expanding hot and dense matter, including the formation of correlations and clusters.
A strict quantum statistical approach to this nonequilibrium process is not available at present. A possible method is the Zubarev nonequilibrium statistical operator which is able to describe the formation of correlations in the expanding hot and dense matter [13,114]. First steps to describe cluster formation (A < 4) in expanding matter have been performed [115], but have to be worked out further, in particular to include larger clusters (A ≥ 4). Kinetic codes based on a single-particle description, but also QMD and AMD codes which simulate cluster formation using the coalescence model and simpler concepts, have to be improved to obtain a microscopic description of cluster formation. Work in this direction is in progress [14,116,117]. Alternatively, the freeze-out concept has been used to model the expansion process of the "fireball". Within the chemical freeze-out approach, it is assumed that at a certain instant of the expanding and down-cooling fireball, the composition is frozen because the "chemical" reactions become slow so that the composition remains unchanged. This concept of a sudden freeze-out has been applied to HIC not only at moderate energies, but also at very high energies [14]. We give a short summary of some results obtained from HIC experiments at moderate energies (≈35 MeV/A), where the production of light clusters was measured and interpreted within a freeze-out model. The main issue of these investigation was to show the relevance of clustering in nuclear systems and the necessity to describe medium effects. In particular, the equation of state of nuclear matter is of interest, at moderate temperatures (T ≤ 20 MeV) and at subsaturation densities.
In contrast to the standard treatment of symmetry energy within mean-field approaches, see, for instance, [123], it does not vanish in the zero density limit, but is significantly determined by cluster formation in the low-density region. This is a very obvious result, and more or less trivial in the density region where medium effects can be neglected. However, going to higher densities, medium effects, in particular the dissolution of bound states, must be included so that a smooth transition to the near-saturation density region is expected, where mean-field approaches can be applied. It has been shown in [118][119][120][121][122] that general expressions for the symmetry energy can be obtained which reproduce the results of experiments at low densities, determined by cluster formation, but agree with mean-field approaches at high densities. A cluster-virial approach may improve the description in the intermediate region.
The direct observation of medium effects in the nuclear matter EoS from HIC experiments is more involved. The composition, in particular the yield of light clusters (A ≤ 4), is obtained in the low-density limit from a mass-action law using the binding energies of free nuclei. We discussed this approach above, Section 2, as nuclear statistical equilibrium (NSE). A special ratio of cluster yields Y i , the so-called chemical constant K A = Y A /(Y Z p Y N n ), can be considered which in this low-density limit is solely a function of the temperature and the volume, but not the chemical potentials. Because in chemical equilibrium the chemical potentials of the cluster A, consisting of Z protons and N = A − Z neutrons, are related as µ A = Zµ p + Nµ n , the chemical potentials cancel in the low-density limit. This simple approach, neglecting any density effects, has been disproved by experiments with HIC [124]. The reason is the modification of the binding energies if the density is increasing. In particular, self-energy shifts and Pauli blocking lead to the reduction of the binding energy and the dissolution at the so called Mott density. The measured chemical constants can be used for the experimental determination of in-medium cluster binding energies and Mott points in nuclear matter [125].
We can discus these experimental results as clear indication for the need to consider medium effects. Within a QS approach including quasiparticle shifts and correlations in the continuum [17], it was possible to reproduce the data for the chemical constants of the light elements d, t, h, α obtained from the cluster yields. More simple models for medium corrections such as the semiempirical excluded volume concept [75] can be adapted to reproduce the data [126].
Also for these experiments, a nuclear matter EoS is needed which describes cluster formation and the medium modification, as well as the treatment of continuum correlations. A cluster-virial approach may improve the calculation of the EoS in a wide region of the phase diagram.
In conclusion, medium effects, in particular self-energy shifts and Pauli blocking for light clusters, are verified by recent HIC experiments. An improved cluster viral approach should be worked out to describe adequately the contributions of correlations in the continuum, expressed by in-medium scattering phase shifts between the different constituents. The freeze-out model to describe the expanding fireball may be considered as an approximation to treat the nonequilibrium process. Kinetic approaches, for instance transport codes allowing for cluster production (such as QMD or AMD), may be developed further to include in-medium correlation effects. This would allow for a systematic and consistent treatment of HIC experiments. Nevertheless, the correct description of the thermodynamic equilibrium, in particular the cluster viral approach, is a benchmark for all nonequilibrium approaches, and should be advanced in the future.

Deconfinement Transition in Relativistic Heavy-Ion Collisions
The main goal of the theoretical developments towards a unified EoS for quark-hadron matter is to achieve a most reliable prediction for the behavior of warm, dense strongly interacting matter including the deconfinement transition, which is described here as a Mott dissociation of baryonic and mesonic bound states, triggered by the restoration of the dynamically broken chiral restoration. Dynamical chiral symmetry breaking is a rather robust phenomenon that can be described in a broad variety of chiral quark models of the Nambu-Jona-Lasinio type, i.e., with a four-fermion interaction of the current-current form with a sufficiently strong coupling to allow for a nontrivial solution of the gap equation for the quark mass. This is a nonperturbative effect that cannot be obtained in any finite order of perturbation theory and is nowadays an obligatory element of modern descriptions of quark matter. At finite temperatures, however, such models often fail to provide a sensible description of the EoS of QCD matter since the quark matter pressure dominates over the hadronic one already at unphysically low temperatures, due to the lack of a quark confinement mechanism in those models. A simple way out is provided by adopting a bag pressure for mimicking confinement. This is a too simple concept and spoils the beauty of a dynamical description. While the details of the confinement mechanism in QCD are still debated, a viable compromise is provided by the concept of a confining density functional that is based on string-type interactions between color charges and, together with the concept of saturation of color interactions within nearest neighbors (string-flip model) allows for the treatment of such confining interactions in quark matter within a relativistic, selfconsistent quasiparticle model.
The success of the DD2(DD2F)-SFM hybrid EoS for astrophysical applications has been summarized in the previous section. The applications for heavy-ion collisions, in the isospin-symmetric case, are still under way. A necessary requirement for a sensible description of the complex system of an ultrarelativistic heavy-ion collision calls for a numerical simulation code like THESEUS, the three-fluid hydrodynamics-based event simulator extended by UrQMD simulations of final-state interactions [127]. Such a description is most appropriate for the investigation of possible effects of a phase transition in the baryon stopping regime, i.e., when the projectile and target fluids collide and form a highly compressed baryonic matter system, i.e., in the collision energy range √ s NN = 2 − 20 which is covered by the THESEUS program and by the high-statistics collison experiments like NA61, RHIC BES or RHIC FXT and the upcoming NICA and FAIR experiments. Therefore, we have chosen THESEUS as the tool for identifying the QCD phase transition and for investigating the effects of a first-order phase transition on heavy-ion collision observables. Previous studies of this question have been performed with THESEUS [128] and with the three-fluid hydrodynamics code in it [129][130][131] using model EoS of three kinds: purely hadronic, crossover and with a strong first-order transition 5 . The upgrade of the EoS with the DD2(DD2F)-SFM hybrid EoS as described in this work is under way. In particular, one expects modifications from the lower density region of the phase transition and its temperature dependence (for a direct comparison, see the right panel of Figure 2 in [133]) and a better description of flow observables due to the increased stiffening of the high-density part of the new hybrid EoS when compared to the one in Ref. [127]. A main goal of the unified approach to the quark-hadron EoS as outlined in this work is the possibility to obtain an EoS with a second critical endpoint in the QCD phase diagram associated with the chiral/deconfinement transition. The first one, corresponding to the liquid-gas phase transition in nuclear matter is already an integral part of the RDF description of nuclear matter within the DD2(DD2F) part of the present approach. We want to point out that with this approach one may have achieved a systematic formulation of a theory for quark-hadron matter that allows to address also the presently puzzling questions of chemical freezeout of hadrons and nuclei like: 1. Can the success of the thermal statistical model in describing the production of nuclear clusters as measured by the ALICE experiment at LHC [134] be interpreted so that they freeze out directly when hadronizing the QGP so that they may be viewed as preformed multiquark systems already in the QGP? 2. What are the necessary ingredients to understand chemical freezeout of hadrons and clusters kinetically [135]?
With these prospects for the development of a unified quark-hadron matter EoS we want to conclude the present work.

Conclusions
We have outlined an approach to a unified equation of state for quark-hadron matter on the basis of a Φ−derivable approach to the generalized Beth-Uhlenbeck equation of state for a cluster decomposition of thermodynamic quantities like the density. To this end we have summarized the cluster virial expansion for nuclear matter and demonstrated the equivalence of the Green's function approach to the Φ−derivable formulation. As an example, the formation and dissociation of deuterons in nuclear matter was discussed. We have formulated the cluster Φ−derivable approach to quark-hadron matter which allows to take into account the specifics of chiral symmetry restoration and deconfinement in triggering the Mott-dissociation of hadrons. Applications to the phenomenology of nuclear clusters and quark deconfinement in the astrophysics of supernovae and compact stars as well as in heavy-ion collisions are outlined.
This approach unifies the description of a strongly coupled quark-gluon plasma with that of a medium-modified hadron resonance gas description which are shown to be its limiting cases. The developed formalism shall replace the common two-phase approach to the description of the deconfinement and chiral phase transition, where separately developed equations of state for hadronic and quark matter are matched with Gibbs conditions of phase equilibrium. Roughly speaking, one would develop a Ginzburg-Landau-type density functional which allows first, second and higher-order transitions, including crossovers. Examples are the van-der-Waals EoS which has a first-order transition with critical endpoint or the RMF models of nuclear matter for the liquid-gas transition. The cluster virial expansion shall allow a formulation of the quark-hadron transition in a similar way. 5 Recently, a thermodynamically consistent generalization of the excluded-volume improved RDF approach to the hadronic EoS has been suggested which employs a density-and temperature-dependent excluded volume parameter. Within this setting, a second first-order phase transition with a critical endpoint in the QCD phase diagram has been obtained [132]. Such a formulation may be most convenient, e.g., for Bayesian studies of the structure of the QCD phase diagram to be extracted from data of heavy-ion collision experiments.