RMF models with $\sigma$-scaled hadron masses and couplings for description of heavy-ion collisions below 2A GeV

TL;DR

This paper extends relativistic mean-field models with density-dependent hadron masses and couplings to finite temperatures, analyzing nuclear matter properties relevant for heavy-ion collisions below 2A GeV, including phase transitions and particle ratios.

Contribution

It introduces finite-temperature extensions of existing models with hadron mass and coupling dependencies, and studies their implications for heavy-ion collision observables and phase transitions.

Findings

Models reproduce thermodynamic characteristics of the liquid-gas phase transition.

Effects of $elta$ resonances and pion interactions on particle ratios are quantified.

Thermodynamic properties align with experimental data.

Abstract

Within the relativistic mean-field framework with hadron masses and coupling constants dependent on the mean scalar field we study properties of nuclear matter at finite temperatures, baryon densities and isospin asymmetries relevant for heavy-ion collisions at laboratory energies below 2$A$ GeV. Previously constructed (KVORcut-based and MKVOR-based) models for the description of the cold hadron matter, which differ mainly by the density dependence of the nucleon effective mass and symmetry energy, are extended for finite temperatures. The baryon equation of state, which includes nucleons and $\Delta$ resonances is supplemented by the contribution of the pion gas described either by the vacuum dispersion relation or with taking into account the $s$-wave pion-baryon interaction. Distribution of the charge between components is found. Thermodynamical characteristics on $T-n$ plane are considered. The energy-density and entropy-density isotherms are constructed and a dynamical trajectory of the hadron system formed in heavy-ion collisions is described. The effects of taking into account the $\Delta$ isobars and the $s$-wave pion-nucleon interaction on pion differential cross sections, pion to proton and $\pi^-/\pi^+$ ratios are studied. The liquid-gas first-order phase transition is studied within the same models in isospin-symmetric and asymmetric systems. We demonstrate that our models yield thermodynamic characteristics of the phase transition compatible with available experimental results. In addition, we discuss the scaled variance of baryon and electric charge in the phase transition region. Effect of the non-zero surface tension on spatial redistribution of the electric charge is considered for a possible application to heavy-ion collisions at low energies.

Tables

Table 1. Table 1: Coefficients of the energy expansion ( 8 ) for KVORcut03 and MKVOR models.

EoS	${ℰ}_0$	$n_0$	$K$	$m_N^{\ast }\left(n_0\right)$	$J$	$L$	$K^{\prime }$	$K_{\mathrm{sym}}$
	[MeV]	[fm-3]	[MeV]	$\left[m_N\right]$	[MeV]	[MeV]	[MeV]	[MeV]
KVORcut03	$-16$	0.16	275	0.805	32	71	422	-86
MKVOR	$-16$	0.16	240	0.730	30	41	557	-158

Equations

$$g_{Mb}^{\ast }=g_{Mb}{\chi }_{Mb}(\sigma ),m_i^{\ast }=m_i{\Phi }_i(\sigma )$$

$$\displaystyle E[{\{n_{b}\}},f,T]=\frac{m_{N}^{4}f^{2}}{2C_{\sigma}^{2}}\eta_{\sigma}(f)+\frac{1}{2m_{N}^{2}}\Big{[}\frac{C_{\omega}^{2}n_{V}^{2}}{\eta_{\omega}(f)}+\frac{C_{\rho}^{2}n_{I}^{2}}{\eta_{\rho}(f)}\Big{]}$$

$$+\sum _b{E}_{\mathrm{kin}}^b+E_{\mathrm{pion}},$$

$$n_V=\sum _b{x}_{\omega b}{n}_b,n_I=\sum _b{x}_{\rho b}{t}_{3b}{n}_b,$$

$$E_{\mathrm{kin}}^b={\gamma }_b\underset{0}{\overset{\infty }{\int }}\frac{p^{2}dp}{2{\pi }^2}{f}_b(p)\sqrt{p^2+m_b^{\ast 2}},{\gamma }_b=(2s_b+1),$$

$$f_b(p)=\frac{1}{e^{(\sqrt{p^2+m_b^{\ast 2}}-{\mu }_b^{\ast })/T)}+1},n_b={\gamma }_b\underset{0}{\overset{\infty }{\int }}\frac{p^{2}dp}{2{\pi }^2}{f}_b(p),$$

$${\mu }_b^{\ast }={\mu }_b-\frac{C_{\omega }^2{n}_V{x}_{\omega b}}{m_N^2{\eta }_{\omega }(f)}-\frac{C_{\rho }^2{n}_I{x}_{\rho b}{t}_{3b}}{m_N^2{\eta }_{\rho }(f)},$$

$$n_{\pi }=\underset{0}{\overset{\infty }{\int }}\frac{p^{2}dp}{2{\pi }^2}{f}_{\pi }(p),f_{\pi }(p)=\frac{1}{e^{({\omega }_{\pi }(p)+{\mu }_{\pi })/T}-1},$$

$$E_{\mathrm{pion}}=\sum _{\pi }\underset{0}{\overset{\infty }{\int }}\frac{p^{2}dp}{2{\pi }^2}{\omega }_{\pi }(p)f_{\pi }(p)$$

$${\eta }_{\omega }(f)={\Phi }_{\omega }^2(f)/{\chi }_{\omega N}^2(f),{\eta }_{\rho }(f)={\Phi }_{\rho }^2(f)/{\chi }_{\rho N}^2(f),$$

$${\eta }_{\sigma }(f)=\frac{{\Phi }_{\sigma }^2[\sigma (f)]}{{\chi }_{\sigma N}^2[\sigma (f)]}+\frac{2C_{\sigma }^2}{m_N^4{f}^2}U[\sigma (f)],$$

$${\Phi }_M=1-f,{\Phi }_b(f)=1-x_{\sigma b}{\xi }_{\sigma b}{m}_{N}f/m_b,$$

$${\xi }_{\sigma b}={\chi }_{\sigma b}/{\chi }_{\sigma N},f=g_{\sigma N}{\chi }_{\sigma N}\sigma /m_N,$$

$$x_{\omega \Delta }=x_{\rho \Delta }=1,x_{\varphi \Delta }=0.$$

$$U_b(n_0)=\frac{C_{\omega }^2{x}_{\omega b}{n}_0}{{\eta }_{\omega }(f(n_0))m_N^2}-(m_N-m_N^{\ast }(n_0))x_{\sigma b}.$$

$$\mathcal{E}={\mathcal{E}}_0+\frac{1}{2}K{ϵ}^2-\frac{1}{6}{K}^{\prime }{ϵ}^3+{\beta }^2{\tilde{\mathcal{E}}}_{\mathrm{sym}}(n)+...,$$

$${\tilde{\mathcal{E}}}_{\mathrm{sym}}(n)=J+Lϵ+\frac{K_{\mathrm{sym}}}{2}{ϵ}^2+\dots ,$$

$${\mu }_n-{\mu }_p=\frac{\partial E[n_p,n_n]}{\partial n_n}-\frac{\partial E[n_p,n_n]}{\partial n_p}\equiv {\mu }_Q.$$

$${\omega }^2=m_{\pi }^2+k^2+\text{mbox}Re\Pi (\omega ,k,n_b,T),$$

$${\Pi }_s^{{\pi }^{-}}=(n_n-n_p)\omega /(2f_{\pi }^2),{\Pi }_s^{{\pi }^0}=0,$$

$${\omega }_{{\pi }^{\mp }}(k)=\pm \frac{n_n-n_p}{4f_{\pi }^2}+\sqrt{m_{\pi }^2+k^2+{\left(\frac{n_n-n_p}{4f_{\pi }^2}\right)}^2},$$

$${\omega }_{{\pi }^0}(k)\equiv {\omega }_k=\sqrt{m_{\pi }^2+k^2}.$$

$$Y_Z=\frac{n_B^Q}{n}+\frac{n_{\pi }^Q}{n},\frac{n_B^Q}{n}=Y_p+\frac{2n_{{\Delta }^{++}}+n_{{\Delta }^{+}}-n_{{\Delta }^{-}}}{n},$$

$$\frac{n_{\pi }^Q}{n}=\frac{n_{{\pi }^{+}}-n_{{\pi }^{-}}}{n},$$

$${\mathcal{E}}_{c.m.}=m_N\sqrt{1+\frac{2A_{\mathrm{p}}{A}_{\mathrm{t}}{\mathcal{E}}_{\mathrm{lab}}}{(A_{\mathrm{p}}+A_{\mathrm{t}})m_N}}-m_N,$$

$$E(n_m,T_m)/n_m={\mathcal{E}}_{c.m.}+m_N+{\mathcal{E}}_{\mathrm{bind}}.$$

$$s=S(n(t),T(n(t)))/n\simeq S(n_m,T_m)/n_m.$$

$$\bar{v}=\frac{\sqrt{\Delta \mathcal{E}+2m\Delta \mathcal{E}}}{m+\Delta \mathcal{E}},\Delta \mathcal{E}={\mathcal{E}}_m-{\mathcal{E}}_{b.up},$$

$$\displaystyle\omega_{k}\frac{d\sigma}{d^{3}k}=15.64\,\frac{(2\pi)^{3}}{V_{{\rm b.up}}}\frac{A^{5/3}}{n_{\rm b.up}}\,\frac{{\rm{mb}}\cdot{\rm{GeV}}}{{\rm{sr}}\cdot({\rm{GeV}}/c)^{3}}\Big{\{}\omega_{k}\frac{d^{3}N_{\pi^{-}}}{d^{3}k}$$

$$\displaystyle+\omega_{k}\frac{d^{3}N_{\Delta}}{d^{3}k}\Big{\}}\,,$$

$$\frac{d^3{N}_{{\pi }^{-}}}{d^{3}k}=\frac{V_{b.up}}{(2\pi {)}^3}\frac{{\Gamma }_{{\pi }^{-}}(k)}{e^{[({\omega }_{{\pi }^{-}}(k)-{\mu }_Q)/T]}-1}.$$

$${\Gamma }_{{\pi }^{-}}(k)=\frac{2{\omega }_k}{\left(2\omega (k)-\frac{\partial \text{mbox}Re{\Pi }_{{\pi }^{-}}}{\partial \omega }{|}_{\omega ={\omega }_{{\pi }^{-}}(k)}\right)}.$$

$${\omega }_k\frac{d^3{N}_{\Delta }}{d^{3}k}=\frac{V_{b.up}}{(2\pi {)}^3}\sum _r\sum _{n=1}^{\infty }\frac{(-1{)}^{n-1}{b}_r{g}_r{m}_r^{\ast 2}}{k{\Lambda }^{1/2}}$$

$$\displaystyle\times\Big{[}\frac{nET+T^{2}}{n^{2}}e^{-[n(E-\mu_{r}^{*})/T]}\Big{]}^{E^{-}}_{E^{+}},$$

$$\lambda (x,y,z)=(x-y-z{)}^2-4yz,\Lambda \equiv \lambda (m_r^{\ast 2},m_{\pi }^2,m_N^{\ast 2}),$$

$$\displaystyle E^{\pm}(k)=\frac{1}{2m_{\pi}^{2}}\Big{[}(m_{r}^{*2}+m_{\pi}^{2}-m_{N}^{*2})E_{\pi}(k)\pm\Lambda^{1/2}k\Big{]}.$$

$$\frac{n_{\mathrm{p}-\mathrm{like}}}{n}=Y_p+\frac{1}{3n}(n_{{\Delta }^0}+2n_{{\Delta }^{+}}+3n_{{\Delta }^{++}}),$$

Full text

1 Introduction

Knowledge of the equation of state (EoS) of cold dense hadronic matter is required for the description of atomic nuclei and neutron stars after minutes-hours since their formation, and EoS of warm and hot hadron matter is required for the description of supernovae and heavy-ion collisions. Nowadays there exists a vast number of EoSs and a large set of experimental and observational constraints, which an appropriate EoS should fulfill Klahn:2006ir . No one of existing EoSs satisfies all the known constraints. Flexible phenomenological approaches to EoSs are introduced within relativistic mean-field (RMF) models with density dependent couplings, see Typel ; Typel2005 ; Voskresenskaya ; Oertel:2016bki and refs. therein, and with hadron masses and coupling constants dependent on the mean scalar field Kolomeitsev:2004ff . The latter model has been generalized to the description of the isospin-symmetric hot hadronic matter including various baryon resonances and bosonic excitations Khvorostukhin:2006ih ; Khvorostukhin:2008xn ; Khvorostukhin:2010aj and was applied to the description of heavy-ion collisions in a broad range of collision energies. Isospin-asymmetric matter (IAM) was not considered in mentioned works. Bosonic excitations were considered in the ideal gas model.

Recent measurements of masses of the most massive binary pulsars demonstrated that the maximum compact star mass, predicted by an EoS, should exceed $2M_{\odot}$. It was found that PSR J1614-2230 has the mass $M=1.928\pm 0.017\,M_{\odot}$ Demorest:2010bx ; Fonseca:2016tux and PSR J0348+0432 has the mass $M=2.01\pm 0.04\,M_{\odot}$ Antoniadis:2013pzd . These measurements rule out many soft EoSs of the purely nucleon matter. However additional degrees of freedom may appear in dense neutron-star interiors, such as hyperons, $\Delta$ isobars and meson condensates. This leads to an additional softening of the EoS of the beta-equilibrium matter (BEM) resulting in a decrease of the maximum neutron-star mass. On the other hand, description of the particle flow in heavy-ion collisions requires a rather soft EoS of the isospin-symmetric matter (ISM) Danielewicz:2002pu . Thereby it is challenging to construct an EoS, which would simultaneously satisfy the maximum neutron-star mass constraint together with the flow constraint.

The working model with hadron masses and coupling constants dependent on the mean scalar field was constructed in Kolomeitsev:2004ff and labeled in Klahn:2006ir as KVOR model. It satisfies the flow constraint for ISM and yields the maximum neutron-star mass $M\simeq 2.01M_{\odot}$ for BEM, however only if no baryons other than nucleons are included into consideration. In our subsequent works Maslov:2015msa ; Maslov:2015wba ; Kolomeitsev:2016ptu ; Kolomeitsev:2017gli we constructed RMF models of the cold hadronic matter of arbitrary isospin composition with effective hadron masses and coupling constants dependent on the scalar field with hyperons and $\Delta$ resonances taken into account Kolomeitsev:2016ptu , as well as with the charged $\rho$ condensate Kolomeitsev:2017gli , which successfully pass the maximum neutron-star mass constraint and the flow constraint simultaneously with many other constraints.

Various characteristics of heavy-ion collisions at energies below a few $A$ GeV have been extensively studied within the expanding fireball framework Siemens:1978pb ; Friedman:1981qm ; Mishustin:1983nv and then within the ideal hydrodymamics and in various transport models, cf. Arsene:2006vf ; Buss:2011mx and refs. therein. Densities reached at such collision energies are typically $\stackrel{{\scriptstyle\scriptstyle<}}{{\phantom{}{}_{\sim}}}3n_{0}$, where $n_{0}$ is the nuclear saturation density, and temperatures are below the pion mass $m_{\pi}\simeq 140$ MeV. Only EoS for ISM has been studied within the expanding fireball framework Siemens:1978pb ; Friedman:1981qm ; Mishustin:1983nv . Reference Friedman:1981qm used a variational theory of nuclear matter for a description of nucleons in ISM and Voskresensky:1997vq used the original Walecka RMF model Walecka:1974qa , whereas pions were considered with taking into account $p$-wave pion-baryon medium polarization effects. Then Voskresensky:1989sn ; Migdal:1990vm ; Voskresensky:1993ud exploited a modified Walecka RMF model of Cubero:1987pr for the nucleon ISM. A comparison with the data on pion and nucleon spectra available at that time demonstrated advantages of such description. At the freeze-out stage an influence of the effects of the nuclear polarization on the pion spectra was considered within the prompt breakup model Senatorov:1989cg . However, being extended to describe BEM, such EoSs do not satisfy modern data on the high masses of cold compact stars.

In this work we generalize the models with $\sigma$-scaled hadron masses and couplings developed for the description of the cold BEM in Maslov:2015msa ; Maslov:2015wba ; Kolomeitsev:2016ptu ; Kolomeitsev:2017gli , now for the case of the ISM and IAM formed in heavy-ion collision reactions for collision energies $\stackrel{{\scriptstyle\scriptstyle<}}{{\phantom{}{}_{\sim}}}2A$ GeV, so the reached baryon densities are $\stackrel{{\scriptstyle\scriptstyle<}}{{\phantom{}{}_{\sim}}}3n_{0}$ and temperatures are below $m_{\pi}$. In heavy-ion collisions the strangeness is approximately conserved. Thereby, the hyperon contribution to the thermodynamical values, $\propto e^{-2m_{H}^{*}/T}$, where $m_{H}^{*}$ is the hyperon effective mass, can be neglected. Effect of boson $\sigma$, $\omega$, $\rho$ excitations can be also disregarded for $n\stackrel{{\scriptstyle\scriptstyle<}}{{\phantom{}{}_{\sim}}}3n_{0}$ and $T<m_{\pi}$. The temperature dependence can be then included only in the nucleon and $\Delta$ kinetic energy terms and in pion quantities. We use a simplified expanding fireball framework. In a subsequent work we plan to check the validity of our EoSs in actual hydrodynamical calculations. Up to now simulations of heavy-ion collisions have been done within ideal hydrodynamics with various EoSs of isospin symmetric matter with pions treated as particles obeying the vacuum dispersion law, cf. Ivanov:2016xev . As the first step, in the present work pions will be treated either as ideal gas of the particles obeying the vacuum dispersion law or as the quasiparticle gas with the $s$-wave pion-baryon interactions included. The latter contribute only for IAM. More involved effects of the $p$-wave pion-baryon interaction will be disregarded.

In the heavy-ion collisions at very low collision energies (${\cal{E}}_{\rm lab}\stackrel{{\scriptstyle\scriptstyle<}}{{\phantom{}{}_{\sim}}}(200-300)A$ MeV) in the expansion stage of the matter at nucleon densities $n(t)<n_{0}$ and at low temperatures, $T(t)\stackrel{{\scriptstyle\scriptstyle<}}{{\phantom{}{}_{\sim}}}(15-20)$ MeV, there may occur a spinodal instability during the first-order liquid-gas (LG) phase transition Ropke:1982vzx . Possible effects of the supercooled gas and superheated liquid phases, as well as the effects of the spinodal instabilities, have been discussed in Ropke:1982vzx ; Schulz:1983pz ; BS ; Panagiotou:1984rb . The nuclear LG transition phenomenon remained an arena of intense research both on theoretical and experimental sides during subsequent years, cf. Muller:1995ji ; Li:1997px ; Ducoin:2005aa ; Alam:2017krb . Occurrence of a negative specific heat at constant pressure was reported, as the first experimental evidence of the LG phase transition in heavy-ion collision reactions DAgostino:1999dod ; Schmidt:2000zs . For a review of this interesting topic see Chomaz:2003dz and more recent works Skokov:2008zp ; Skokov:2009yu ; Skokov:2010dd ; Voskresensky:2010gf ; Borderie:2018fsi . Isospin dependence was studied in Colonna:2002ti ; Ducoin:2005aa . Effects of a finite surface tension were disregarded. Below we apply our models with effective hadron masses and coupling constants dependent on the scalar field also to describe the LG first-order phase transition occurring at a low temperatures and densities. First we assume the surface tension to be zero and then include effects of the non-zero surface tension, which may result in preparation of the pasta-like structures in heavy-ion collisions.

The manuscript is organized as follows. In section 2 we present the RMF model with $\sigma$-scaled hadron masses and couplings generalized for description of the hadron matter at arbitrary isospin composition for $T\neq 0$, including pion gas. For specificity we further use the KVORcut03-based and MKVOR*-based models of EoS Kolomeitsev:2016ptu in the region of the baryon densities $\stackrel{{\scriptstyle\scriptstyle<}}{{\phantom{}{}_{\sim}}}3n_{0}$ and temperatures $T<m_{\pi}$. In Sect. 3 we present results obtained in a simplified expanding fireball model. This simplified description allows us to demonstrate many qualitative and quantitative effects. In Sect. 4 we focus on the description of the region of the LG phase transition first considering ISM and then IAM. Consideration is first performed within the RMF framework and then effects of fluctuations are discussed. Importance of the surface tension effects will be then emphasised. Then in Sect. 5 we formulate our conclusions. For completeness in Appendix Appendix we indicate effects of the polarization of the medium due to the $p$-wave pion-baryon interaction, which were disregarded in the present study.

2 RMF models with $\sigma$-scaled hadron masses and couplings. EoS of hadron matter

in the region $T<m_{\pi}$, $n<3n_{0}$, $0.4\leq Z/A\leq 0.5$

We use the framework proposed in Kolomeitsev:2004ff and then developed further in Maslov:2015msa ; Maslov:2015wba ; Kolomeitsev:2016ptu for $T=0$ and an arbitrary isospin composition and generalized in Khvorostukhin:2006ih ; Khvorostukhin:2008xn ; Khvorostukhin:2010aj for the case of $T\neq 0$, but only for ISM. In the present work we focus on the description of matter produced in heavy-ion reactions at collision energies less than few $A$ GeV. Thereby we study the ISM, for which $Y_{Z}=Z_{\rm tot}/A_{\rm tot}\simeq 0.5$, and the IAM matter, when $0.4\stackrel{{\scriptstyle\scriptstyle<}}{{\phantom{}{}_{\sim}}}Y_{Z}<0.5$, where $Z_{\rm tot}$ is the total charge of the colliding nuclei and $A_{\rm tot}$ is the corresponding baryon number.

Our model is a generalization on the case $T\neq 0$ of the non-linear Walecka model with effective coupling constants and hadron masses

[TABLE]

dependent on the scalar field $\sigma$. Here $M=(\sigma,\omega,\rho)$ denotes mesons, for which we use the mean-field solutions of the equations of motion, $b=(N,\Delta)$ lists the included baryon species, nucleons $N=(p,n)$ and $\Delta$ isobars, $i=(M,b)$. We neglect a contribution of hyperons and anti-baryons, which are tiny for collision energies under consideration, $\propto e^{-2m^{*}_{H}/T}$ for hyperons due to the strangeness conservation Khvorostukhin:2006ih ; Khvorostukhin:2008xn and $\propto e^{-2m^{*}_{b}/T}$ for anti-baryons. In absence of the hyperon occupations there is no contribution of the $\phi$ meson mean field. Besides, we include pions $\pi=\{\pi^{-},\pi^{0},\pi^{+}\}$, as lightest among pseudo-Goldstone particles. Other pseudo-Goldstone particles $G=(\pi,K,\eta)$ and their heavier partners $K^{*}$ and $\eta^{\prime}$ are not included, since their contributions, $\propto e^{-m_{G}/T}$, remain tiny for $T<m_{\pi}$. The quantities $\chi_{Mb}(\sigma)$ and $\Phi_{i}(\sigma)$ are the dimensionless scaling functions have been fitted in Maslov:2015msa ; Maslov:2015wba ; Kolomeitsev:2016ptu for the best description of the cold baryon matter.

Using mean-field solutions for meson fields we present the energy density of the hadronic system as Maslov:2015msa ; Maslov:2015wba ; Kolomeitsev:2016ptu

[TABLE]

where $s_{b}$ is the baryon spin, $\mu_{b}=\mu_{B}-Q_{b}\mu_{Q}$ is the chemical potential for the given baryon species $b$, $\mu_{B}$ is the baryon-charge chemical potential, $\mu_{Q}$ is the chemical potential of a negative electric charge, $Q_{j}$ is the electric charge of a particle $j$, $t_{3b}=Q_{b}-1/2$ is the isospin projection of baryon $b$.

For given $\mu_{B}$ and $\mu_{Q}$ eqs. (3) can be solved to find the particle densities and their effective chemical potentials $\mu_{b}^{*}$. Then the definitions of the total baryon density $n=\sum_{b}n_{b}$ and charge density $n_{Q}=\sum_{b}Q_{b}n_{b}+\sum_{\pi}Q_{\pi}n_{\pi}$ are treated as equations for finding the chemical potentials $\mu_{B},\mu_{Q}$ for given $n,n_{Q}$. Here $n_{\pi}$ are the pion number densities

[TABLE]

where $\omega_{\pi}(p)$ is the dispersion relation of a pion and $\mu_{\pi}=Q_{\pi}\mu_{Q}$ is its chemical potential. For ISM $\mu_{\pi}=0$. The term

[TABLE]

is the contribution of $\pi$ species to the energy density.

For any set of baryon concentrations and the temperature equation $\partial E[{\{n_{b}\}},f,T]/\partial f=0$ is solved to find the equilibrium value of the scalar field.

The effective hadron masses and scaling functions of mesons enter the volume part of the thermodynamic quantities only in combinations $C_{M}={g_{MN}m_{N}}/{m_{M}}$,

[TABLE]

where the self-interaction potential ${U}(\sigma)$ employed in standard RMF models is included into the definition of the scaling function $\eta_{\sigma}(f)$. The scaling functions of the mass are

[TABLE]

where

[TABLE]

and thereby $\Phi_{N}=\Phi_{M}$. We suppose that $\chi_{\omega b}(f)=\chi_{\omega N}(f)$ , $\chi_{\rho b}(f)=\chi_{\rho N}(f)$. Explicit expressions for the scaling functions $\eta_{M}(f)$ are presented in Kolomeitsev:2016ptu .

The coupling constant ratios are introduced as $x_{Mb}=g_{Mb}/g_{MN}$. The vector-meson coupling constants to $\Delta$s are chosen following the quark SU(6) symmetry:

[TABLE]

The $\Delta$ coupling constants with the scalar field are deduced from the values of the optical potentials $U_{b}$ in ISM at the saturation density $n=n_{0}$ given by

[TABLE]

The value of the $\Delta$ potential $U_{\Delta}(n_{0})$ is poorly constrained by the data. As in Kolomeitsev:2016ptu , we use $U_{\Delta}(n_{0})=-50\,{\rm MeV}$ as the most realistic estimate. Models including $\Delta$s will be denoted by ”$\Delta$” suffix.

Also we assume that the size of the system under consideration is such that the volume part of the thermodynamic quantities of our interest is much larger than the surface part. Moreover, we first disregard finite-size Coulomb effects compared to the strong-interaction effects. The former effects will be discussed in Sect. 4. Focusing on the description of heavy-ion collisions we do not include lepton terms.

In ref. Maslov:cut we demonstrated that within an RMF model the EoS becomes stiffer for $n>n^{*}>n_{0}$, if a growth of the scalar field as a function of the density is quenched and the nucleon effective mass becomes weakly dependent on the density for $n>n^{*}$. In Maslov:cut such a quenching was achieved by adding to the scalar potential a rapidly rising function of $f$ at $f>f^{*}$, where $f^{*}$ is $f$ corresponding to $n=n^{*}$. We called it the cut-mechanism. In Ref. Maslov:cut the cut-mechanism is realized in the $\sigma$ sector. We focus now on two models based on KVORcut03 and MKVOR* models proposed in Maslov:2015msa ; Maslov:2015wba ; Kolomeitsev:2016ptu . These models proved to satisfy many constraints on the hadronic EoS. In neutron-star matter for large densities the hyperons and $\Delta$ baryons appear in our models. These two models utilize the cut-mechanisms in the $\omega$ and $\rho$ sectors, respectively. The cut mechanism in $\rho$ sector is implemented in MKVOR-based models in order simultaneously to keep the EoS not too stiff in ISM (to satisfy the flow constraint from heavy-ion collisions Danielewicz:2002pu ) and to do the EoS as stiff as possible in the BEM to safely fulfill the constraint on the maximum mass of neutron stars. The $\rho$ mean field is coupled to the isospin density that makes the $f$-saturation mechanism very sensitive to the composition of the BEM. Incorporating $\Delta$ baryons we use the MKVOR* extension of the MKVOR model Kolomeitsev:2016ptu , which prevents the effective nucleon mass from vanishing at high density.

Free parameters of the model are fitted to reproduce properties of the cold nuclear matter near the saturation point. These parameters are defined as the coefficients of the Taylor expansion of the energy per nucleon for $T=0$,

[TABLE]

in terms of small $\epsilon=(n-n_{0})/3n_{0}$ and $\beta=(n_{n}-n_{p})/n$ parameters. The parameters for the MKVOR* and MKVOR models are identical. The properties of the KVORcut03 and MKVOR* models, which we exploit in this work, at the nuclear saturation density $n_{0}$ are illustrated in Table 1, where we collect coefficients of the expansion of the nucleon binding energy per nucleon near $n_{0}$.

For the difference of the neutron and proton chemical potentials we get

[TABLE]

The pion quasiparticle spectrum is determined as a solution of the dispersion equation Migdal:1990vm ; Voskresensky:1993ud

[TABLE]

where $\Pi(\omega,k,n_{b},T)$ is the pion polarization operator in the baryon medium. In this work we will consider pions as an ideal gas of quasiparticles including for IAM only their $s$-wave interaction with baryons, being determined by the so called Weinberg-Tomozawa term. For ISM the $s$-wave $\pi N$ interaction is suppressed Baym:1975tm ; Migdal:1978az ; Friedman:2019zhc ; Migdal:1990vm ; Voskresensky:1993ud . Including only $s$-wave pion-nucleon interactions the retarded pion polarization operator is given by

[TABLE]

where $f_{\pi}\simeq 92.4$MeV is the pion weak decay constant, cf. Migdal:1990vm ; Kolomeitsev:2002gc , and the $\pi^{\mp}$ spectrum is thereby as follows

[TABLE]

Models with pion quasiparticles treated following eq. (11) will be labeled by “$\pi_{\rm WT}$” suffix and models with pions described by the vacuum dispersion law we label by “$\pi_{\rm free}$” suffix, respectively. Important role of the $p$-wave pion-baryon interactions has been intensively studied in Baym:1975tm ; Migdal:1978az ; Migdal:1990vm . This issue will be briefly reviewed in the Appendix.

One of the main differences between the KVORcut-based and MKVOR*-based models we consider here is the behavior of the nucleon effective mass with the density. Therefore, on left panel in fig. 1 we show the baryon density dependence of the scaling function $\Phi_{N}(f)=m_{N}^{*}/m_{N}=m^{*}_{M}/m_{M}=1-f$ in ISM calculated for various temperatures in the models KVORcut03$\Delta\pi_{\rm free}$ and MKVOR*$\Delta\pi_{\rm free}$. We see that the density dependence of this quantity is significant for $n\stackrel{{\scriptstyle\scriptstyle<}}{{\phantom{}{}_{\sim}}}3n_{0}$, whereas the temperature dependence is moderate for $T<m_{\pi}$. In the KVORcut03$\Delta\pi_{\rm free}$ model $m^{*}_{\sigma}\stackrel{{\scriptstyle\scriptstyle>}}{{\phantom{}{}_{\sim}}}400$ MeV in the density and temperature interval under consideration, whereas in the MKVOR*$\Delta\pi_{\rm free}$ model we have $m^{*}_{\sigma}\stackrel{{\scriptstyle\scriptstyle>}}{{\phantom{}{}_{\sim}}}200$ MeV.111We should note that $m_{\sigma}^{*}$ here is the mass-coefficient of the $\sigma$ mean field rather than the effective mass of $\sigma$ excitations, $m_{\sigma}^{\rm{part}*}$, cf. Khvorostukhin:2006ih ; Khvorostukhin:2008xn . Thus the effective masses of the $\sigma,\omega,\rho$ mesons remain significantly larger than $m_{\pi}$, and the thermal contribution of the $\sigma$, $\omega$, $\rho$, $\phi$ excitations, which is $\propto e^{-m_{M}^{*}(n,T)/T}$, can be neglected for temperatures and densities we consider here. The curves computed for $\pi_{\rm free}$ and $\pi_{\rm WT}$ models prove to be visually almost not distinguishable for $0.4<Y_{Z}<0.5$. Therefore below we mainly focus consideration on $\pi_{\rm free}$ models. Besides, we note that the $Y_{Z}$ dependence of $1-f$ proves to be very weak in the interval $0.4\leq Y_{Z}\leq 0.5$ of our interest. The curves computed for $Y_{Z}=0.4$ and for $Y_{Z}=0.5$ are visually almost not distinguishable.

The temperature dependence of the nucleon energy is $\propto[1+(T/\epsilon_{{\rm F},N})^{2}]$ and thereby it is essential already for $T<\epsilon_{{\rm F},N}$, where $\epsilon_{{\rm F},N}$ is the nucleon Fermi energy ($\epsilon_{{\rm F},N}\sim 40(n/n_{0})^{2/3}$ MeV for $n\sim n_{0}$, $Y_{Z}=0.5$ and for the Landau nucleon effective mass $\simeq m_{N}$). The energy of the pion ideal gas is $\propto e^{-m_{\pi}/T}$ and becomes significant for $T\stackrel{{\scriptstyle\scriptstyle>}}{{\phantom{}{}_{\sim}}}0.5\,m_{\pi}$. The contribution of the $\Delta$s to the energy is suppressed compared to the nucleon one as $4e^{-(m^{*}_{\Delta}-m_{N}^{*}+Q_{\Delta}\mu_{Q})/T}$. All these contributions are included in our work. The first-order phase transition to the $\Delta$ resonance matter considered in Kolomeitsev:2016ptu for $T=0$ does not occur for $T<m_{\pi}$, $n\stackrel{{\scriptstyle\scriptstyle<}}{{\phantom{}{}_{\sim}}}3n_{0}$, for the $\Delta$ optical potential $U_{\Delta}(n_{0})=-50$ MeV that we use in this work. In Khvorostukhin:2008xn within the model, where pions are treated with the vacuum dispersion law, the effect of the nonzero $\Delta$ width was estimated as not significant. Therefore in what follows within our RMF-based model the $\Delta$ resonances will be treated as quasiparticles with an effective mass. Concluding this discussion, for temperatures $T$ below $m_{\pi}$ and for densities $n\stackrel{{\scriptstyle\scriptstyle<}}{{\phantom{}{}_{\sim}}}3n_{0}$ of our interest here, the temperature dependence can be included only in the nucleon and $\Delta$ isobar quasiparticle contributions and the pion kinetic energy terms, which within the $\pi_{\rm free}$ model are described by the free dispersion law.

Our KVORcut03-based and MKVOR*-based models differ also by the density dependence of the symmetry energy coefficient. On the right panel in fig. 1 we show the baryon density dependence of the symmetry energy coefficient $\widetilde{\mathcal{E}}_{\rm sym}(n)$ derived in our models KVORcut03 and MKVOR* for ISM. The density dependence of the symmetry energy in the KVORcut03-based models, where the cut-mechanism is implemented in the $\omega$ sector Kolomeitsev:2016ptu for $n>n^{*}=3n_{0}$, is rather smooth (dashed line). In the MKVOR*-based models, where the cut-mechanism is used in the $\rho$ sector Kolomeitsev:2016ptu responsible for the symmetry energy, the dependence of $\widetilde{\mathcal{E}}_{\rm sym}$ on $n$ becomes sharp for $n>(2.5-4)n_{0}$ (solid line). In Ma:2018xjw the dramatic change in the density dependence of the nuclear symmetry energy for $n$ above some value $n_{1/2}$ varied in the interval $(2-4)n_{0}$ was associated with the change of the topology mimicking the baryon-quark continuity. By the dash-dotted line in figure is shown the symmetry energy coefficient obtained in model Ma:2018xjw for $n_{1/2}=4n_{0}$. The resulting density dependence is similar to that we obtain in the MKVOR-based models.

In the heated dense nuclear system formed in the collision of nuclei, initial charge per baryon $n_{Z}$ is redistributed among all the involved electrically charged hadrons appeared for $T\neq 0$, cf. Muller:1995ji ; Li:1997px , following minimum of the free energy expressed in $n,T$ variables. We use the following decomposition of the system charge:

[TABLE]

where $Y_{Z}$ is the ratio of the initially fixed number of protons to the fixed total baryon number, $Y_{p}=n_{p}/n$ is the ratio of the number of protons to that of the baryons inside the thermal system, ${n^{Q}_{B}}/{n}$ is the ratio of the total charge of the baryon subsystem to the total baryon density, that includes contribution of charged $\Delta$ isobars, and ${n^{Q}_{\pi}}/{n}$ is the ratio of the total charge of the pions to the total baryon density.

In fig. 2 we show the charge per baryon for various species for the model KVORcut03$\Delta\pi_{\rm free}$ on the left panel and for the model MKVOR*$\Delta\pi_{\rm free}$ on the right panel. We take $Y_{Z}=0.4$ as an example relevant to the matter formed in heavy-ion collisions. At zero temperature $Y_{Z}=Z/A=Y_{p}$, since the $\Delta$ isobars do not appear in both our models for $n\stackrel{{\scriptstyle\scriptstyle<}}{{\phantom{}{}_{\sim}}}3\,n_{0}$ Kolomeitsev:2016ptu , which we consider here, and there are no pions for $T=0$. With an increase of the temperature the abundance of $\Delta$ resonances and pions increases. For $Y_{Z}<0.5$ the charge chemical potential $\mu_{Q}$ is positive, and there appears an excess of $\pi^{-}$ respectively to $\pi^{+}$. The ratio $-n_{Q}^{\pi}/n\propto e^{-m_{\pi}/T}(e^{\mu_{Q}/T}-e^{-\mu_{Q}/T})/n$ (lines labeled $\pi$) increases with the temperature $T$. The pion charge per unit of baryon density proves to be higher for smaller density. Unlike in the neutron-star matter, the $\Delta^{++},\Delta^{+},\Delta^{0},\Delta^{-}$ subsystem in heavy-ion collisions at $0.4\stackrel{{\scriptstyle\scriptstyle<}}{{\phantom{}{}_{\sim}}}Y_{Z}\stackrel{{\scriptstyle\scriptstyle<}}{{\phantom{}{}_{\sim}}}0.5$ remains positively charged. Indeed, $n_{Q}^{\Delta}/n\propto e^{(-m^{*}_{\Delta}+m_{N}^{*})/T}[2e^{-2\mu_{Q}/T}+e^{-\mu_{Q}/T}-e^{\mu_{Q}/T}]$. Since the value of $\mu_{Q}$ remains rather small for $n\stackrel{{\scriptstyle\scriptstyle<}}{{\phantom{}{}_{\sim}}}2\,n_{0}$, the $\Delta^{++},\Delta^{+},\Delta^{0}$ and $\Delta^{-}$ abundances remain close to each other. Indeed, as we see in fig. 2 excess of the positively charged $\Delta$ resonances (lines labeled $\Delta$) increases with an increase of $T$ in agreement with above estimate. The value $Y_{p}$ (lines labeled by $p$) slightly increases with increase of $T$ for $T<\epsilon_{{\rm F},N}$ and then it sharply decreases for higher $T$. The reason of a slight increase of $Y_{p}$ for low $T$ is that an enhancement of the proton fraction should compensate a small increase with $T$ of the negative pion charge to fulfill the charge conservation condition. However, at $T\stackrel{{\scriptstyle\scriptstyle>}}{{\phantom{}{}_{\sim}}}\epsilon_{{\rm F},N}$ the $\Delta$ isobar concentration increases noticeably and the proton fraction decreases because of the baryon charge conservation. The total charge of the baryon subsystem per baryon ${n^{Q}_{B}}/{n}$ is shown in fig. 2 by the lines labeled “$\Delta+p$”.

For MKVOR*$\Delta\pi_{\rm free}$ model at $n\sim 3\,n_{0}$ for $30\,{\rm MeV}\stackrel{{\scriptstyle\scriptstyle<}}{{\phantom{}{}_{\sim}}}T\stackrel{{\scriptstyle\scriptstyle<}}{{\phantom{}{}_{\sim}}}50\,{\rm MeV}$ the $\Delta$ system is negatively charged. This happens because in the MKVOR* model the symmetry energy coefficient and, correspondingly, $\mu_{Q}$ are larger at such densities than in the KVORcut03 model, see fig. 1. Therefore in agreement with the estimate given above the density of $\Delta^{-}$ becomes greater than the sum of densities of $\Delta^{+}$ and $\Delta^{++}$. The negative pion excess contributes less to the charge conservation at $n\sim 3\,n_{0}$ than for lower $n$ because the number of pions per baryon at a fixed temperature decreases with an increase of the baryon density.

For $n\stackrel{{\scriptstyle\scriptstyle<}}{{\phantom{}{}_{\sim}}}2n_{0}$ all the ratios obtained in both of our models are very close to each other. Only for $n\sim 3\,n_{0}$ the contribution of $\Delta+p$ to the charge excess computed within the KVORcut03$\Delta\pi_{\rm free}$ model proves to be a bit higher than that in MKVOR*$\Delta\pi_{\rm free}$ model.

Note that in MKVOR*$\Delta\pi_{\rm free}$ model a transition to the $\Delta$ resonance-enriched matter may occur in the dense medium. For $T=0$ and for $U_{\Delta}(n_{0})=-50$ MeV that we use in this work the $\Delta$s appear by the crossover at $n\stackrel{{\scriptstyle\scriptstyle>}}{{\phantom{}{}_{\sim}}}n_{c,\Delta}\simeq 4.5n_{0}$, cf. fig. 11 in Kolomeitsev:2016ptu . With an increase of the temperature in the MKVOR*$\Delta\pi_{\rm free}$ model $\Delta$s arise for a lower density and the phase transition becomes the transition of the first order. We found that for $U_{\Delta}(n_{0})=-50\,{\rm MeV}$ at $T=55\,{\rm MeV}$ the first-order phase transition occurs for $n_{c,\Delta}\simeq 3.6\,n_{0}$. The critical density decreases very smoothly with an increase of the temperature and we get $n_{c,\Delta}\simeq 3.3\,n_{0}$ at $T\simeq 140\,{\rm MeV}$, thus $n_{c,\Delta}$ remains above $3n_{0}$, i.e, outside the density range we consider in the given work. However we should point out that some papers argue that $U_{\Delta}$ could be a more attractive. For $U_{\Delta}(n_{0})=-100$ MeV the first-order phase transition to the $\Delta$ resonance matter would occur for $T=0$ already at the density $n\stackrel{{\scriptstyle\scriptstyle>}}{{\phantom{}{}_{\sim}}}n_{c,\Delta}\simeq 2.5n_{0}$, cf. a discussion in Kolomeitsev:2016ptu . In this work we use $U_{\Delta}(n_{0})=-50\,{\rm MeV}$ and consider $n\stackrel{{\scriptstyle\scriptstyle<}}{{\phantom{}{}_{\sim}}}3\,n_{0}$, thereby the first-order phase transition to the $\Delta$ reach matter does not occur within our MKVOR*$\Delta\pi_{\rm free}$ model. In the KVORcut03$\Delta\pi_{\rm free}$ model the phase transition to the $\Delta$ resonance matter does not occur for all relevant values of $U_{\Delta}$ at densities and temperatures we are interested in this work.

3 Simplified model for heavy-ion collisions at ${\cal{E}}_{\rm{lab}}<(1-2)$ $A$ GeV

3.1 EoS and the system dynamics

In a heavy-ion collision, nucleons can be subdivided on “participants,” which intensively interact with each other during the collision, and “spectators”, which remain practically unperturbed. Baryons-participants form a nuclear fireball, which first is compressed (for $t<0$) and then (for $t>0$) is expanded under action of the internal pressure.

The kinetic energy of the projectile nucleus per nucleon (per $A_{\rm p}=A$) in the laboratory system ${\cal E}_{\rm lab}$ is related to the kinetic energy per nucleon in the center-of-mass frame as follows

[TABLE]

$A_{\rm t}$ is the nucleon number of the target nucleus. As in Voskresensky:1989sn ; Migdal:1990vm ; Voskresensky:1993ud , we assume that at energies less than a few $A$ GeV in the laboratory frame the energy in the center-of-mass frame of the nucleus-nucleus collision, ${\cal E}_{\rm c.m.}A_{\rm part}$, which corresponds to the nucleons-participants, is spent on the creation of an initially quasi-equilibrium nuclear fireball resting in the center-of-mass frame at the end of the compression stage, for $t=0$. The energy per baryon, $E(n,T)/n$, as a function of the baryon density $n$ at fixed $T$ has the concave shape in our KVORcut03$\Delta\pi$ and MKVOR*$\Delta\pi$ models, it decreases with increase of $n$, gets a minimum at $n=n_{m}$, and then begins to increase, see fig. 3 below. Following Voskresensky:1989sn ; Migdal:1990vm ; Voskresensky:1993ud we assume that the initial fireball state is characterized by the temperature $T(0)=T_{m}$, and the baryon density $n(0)=n_{m}$ corresponding to the minimum of the energy per baryon, $E(n,T_{m})/n$, as a function of the baryon density for $T=T_{m}$,

[TABLE]

The quantity ${\cal E}_{\rm bind}$ is the binding energy per baryon in a cold nucleus of the nucleon number $A_{\rm part}$, $-16\,{\rm MeV}\leq{\cal E}_{\rm bind}<0$. Below we will use values ${\cal E}_{\rm bind}$ which follow from (2) at ignorance of surface and Coulomb effects, i.e., as would be for very heavy nuclei. Thus we shall take ${\cal E}_{\rm bind}(Y_{Z}=0.5)\simeq-16\,{\rm MeV}$, ${\cal E}_{\rm bind}(Y_{Z}=0.4)=-14.7$ MeV.

The specific entropy is a decreasing function of $n$, see fig. 3 below. Thereby the state of the minimum of the energy on the right branch of ${\cal E}(n)$ (where ${\cal E}=E/n$ increases with $n$) corresponds to the maximum of the entropy on the given isotherm for the states belonging this branch. Moreover, the state $T_{m},n_{m}$ also corresponds to the maximum temperature among all available solutions of the eq. (13) for all curves $E(T={\rm const},n)$ and respectively this state corresponds to the maximum of the stirring of the degrees of freedom possible at assumption of the full stopping of the matter in the center-of-mass frame for $t=0$.

Note that for a weakly non-equilibrium system moving with the velocity $\vec{v}(t,\vec{r})$ the local pressure can be presented as Ivanov:2013uga $P_{\rm n.eq.}=P-\zeta\nabla\vec{v}$, where $\zeta$ is the bulk viscosity and $P$ is the quasi-equilibrium pressure depending on the local temperature and density $T(t,\vec{r})$ and $n(t,\vec{r})$ following a given EoS. Since $\zeta$ is a positive-definite quantity, the non-equilibrium correction to the pressure is positive during the compression stage of the nuclear system and it is negative during the expansion stage, see also Voskresensky:2010qu . So, on the stage of expansion of the fireball in the heavy ion collision the non-equilibrium pressure $P_{\rm n.eq.}$ is in reality smaller than the equilibrium one, $P$, being used in the ideal hydrodynamics. This means that a realistic equilibrium EoS to be used in non-ideal hydrodynamical calculations should be stiffer than the one used to fit experimental data within ideal hydrodynamical simulations. Besides, the entropy in the viscous process increases, whereas it stays constant within the ideal hydrodynamics.

The viscosity effects prove to be rather small at energies less than (1-2)$A$ GeV we study in this work. This conclusion is supported by the analyses of heavy-ion collisions performed in an expanding fireball model Voskresensky:1989sn ; Migdal:1990vm ; Voskresensky:1993ud , by calculations used ideal hydrodynamics in a broad energy range Ivanov:1991te ; Mishustin:1991sp ; Ivanov:2005yw ; Arsene:2006vf ; Buss:2011mx , by simulations done within transport codes Bass:1998ca ; Arsene:2006vf and by estimates of the viscosity Khvorostukhin:2010aj . In the ideal hydrodynamics dynamical trajectories of the system in heavy-ion collisions are characterized by constant initial values of the entropy per $A_{\rm part}$ (i.e. by total entropy density $S$ per net baryon density of baryons-participants, ${s}=S/n$). In our RMF approach with the contribution of the ideal pion gas included we neglect the inelastic processes and thereby the entropy is assumed to be conserved. Thus we have

[TABLE]

From this relation we obtain the dependence $T(n(t))$. Note that actually all the results of this model hold locally, and therefore are applicable to the case of either quasi-homogeneous fireball expansion, or inhomogeneous hydrodynamical expansion with $n=n(t,\vec{r})$ and $T=T(t,\vec{r})$ depending on the space point. Just for the illustration purposes, as in Voskresensky:1989sn ; Migdal:1990vm ; Voskresensky:1993ud , we will further assume that the expansion is uniform.

We use the KVORcut03$\Delta\pi$ and MKVOR∗$\Delta\pi$ models. The ideal pion gas either is described by the vacuum dispersion law in the $\pi_{\rm free}$ model or by the law (11) in the quasiparticle model $\pi_{\rm WT}$ for IAM. The energy density is given by eq. (2). As we have mentioned, from the relation (13) we unambiguously determine quantities $n_{m}$ and $T_{m}$. The values of the energy per baryon and specific entropy ${s}$ are shown in fig. 3 (left and right) for KVORcut03$\Delta\pi_{\rm free}$ model and in fig. 4 (left and right) for MKVOR*$\Delta\pi_{\rm free}$ model as functions of the density at various temperatures. Being computed with the dispersion law (11) in $\pi_{\rm WT}$ model, thermodynamic quantities prove to be visually almost not distinguishable from those calculated in $\pi_{\rm free}$ model for $0.4\leq Y_{Z}\leq 0.5$. Therefore we do not show the curves for $\pi_{\rm WT}$ model in figs. 3, 4. Values of the temperatures are indicated on lines in MeV. The horizontal arrows on the left panels denote the initial energy per baryon in the laboratory system. Thin lines indicate minima of the energies per baryon, which exist in our models for all $T$ and $n$ corresponding to ${\cal E}_{\rm lab}\stackrel{{\scriptstyle\scriptstyle<}}{{\phantom{}{}_{\sim}}}2.1A$ GeV.

With a knowledge of these quantities, the dynamics of the expanding nuclear fireball is determined by the constant value of the entropy per baryon (thin horizontal dashed lines on right panel). Solid bold curves are presented for $Y_{Z}=0.48$ (as for Ar+KCl collisions). Two dashed bold curves on each figure illustrate the case of 0.8$A$ GeV and 0.246$A$ GeV collisions of La+La $(Y_{Z}=0.4)$. Comparing the bold dashed and solid curves for the collisions with ${\cal E}_{\rm lab}=0.8A$ GeV we see that the effect of the $Y_{Z}$ dependence is tiny for $0.4\leq Y_{Z}\leq 0.5$. This is because the symmetry energy of asymmetric matter is approximately $\propto\beta^{2}$ and within the interval $0.4\leq Y_{Z}\leq 0.5$ the quantity $\beta^{2}$ changes from $0.04$ to 0 and the contribution remains negligible. For the neutron-star matter $\beta^{2}\simeq 0.9$ and thereby the symmetry energy gives significant contribution to the total energy. Also in figs. 3 and 4 by dotted bold lines we show the results for KVORcut03$\pi_{\rm free}$ and MKVOR*$\pi_{\rm free}$ models, i.e. without the inclusion of $\Delta$ resonances. We see that the contribution of $\Delta$s becomes noticeable for all densities already at $T_{m}\stackrel{{\scriptstyle\scriptstyle>}}{{\phantom{}{}_{\sim}}}(40-50)$ MeV, which roughly corresponds to ${\cal E}_{\rm lab}\stackrel{{\scriptstyle\scriptstyle>}}{{\phantom{}{}_{\sim}}}(200-300)$ MeV. Compared to the case without $\Delta$s, with $\Delta$s included the values $T_{m}$ and $n_{m}$ are lower for all ${\cal E}_{\rm lab}$, see the thin dotted line connecting minima of the energy. Contrary to that, the initial value of the entropy is larger for the models with $\Delta$s for all ${\cal E}_{\rm lab}$.

We assume that in the fireball, expanding with the velocity $v(t)$ (in reality with $v(\vec{r},t)$), the thermodynamical quasi-equilibrium is sustained up to a certain rather short breakup stage at which the nucleon and pion mean free paths become compatible with the fireball size, or more precisely, the typical expansion time becomes comparable with the inverse collision frequency Mishustin:1983nv ). After that the nucleon and pion momentum distributions can be considered as frozen. We assume that the breakup stage is characterized by the baryon density $n_{\rm b.up}$ and temperature $T_{\rm b.up}$. First fireball models estimated values of the freeze-out densities in the interval $0.5\,n_{0}<n_{\rm b.up}\stackrel{{\scriptstyle\scriptstyle<}}{{\phantom{}{}_{\sim}}}n_{0}$ Gosset:1988na ; DasGupta:1981xx ; Nagamiya:1981sd ; Barz:1982ed ; Mishustin:1983nv . Resonance gas model, cf. Randrup:2006nr , yields $0.3n_{0}<n_{\rm b.up}\stackrel{{\scriptstyle\scriptstyle<}}{{\phantom{}{}_{\sim}}}n_{0}$ in the whole interval of available collision energies. References Voskresensky:1989sn ; Migdal:1990vm ; Voskresensky:1993ud argued that at lowest SIS energies $0.5n_{0}<n_{\rm b.up}\stackrel{{\scriptstyle\scriptstyle<}}{{\phantom{}{}_{\sim}}}0.8n_{0}$ since for higher densities there appears a significant contribution to the $NN$ scattering amplitude from the $NN$ exchange by soft pions with momenta $k\sim p_{{\rm F},N}$, region of larger $n$ is usually called the liquid phase of the pion condensate Voskresensky:1993ud . Besides, the freeze-out densities can be estimated from analysis of the HBT pion interferometry, cf. fig. 2 in Mishra:2007xg . In figs. 3 and 4 the breakup moments are indicated by small dots. The choice of these points is explained further in the text.

The fireball expansion velocity is supposed to be zero at the initial moment and grows with time, since a part of the energy $E(n_{m},T_{m})/n_{m}-E(n(t),T(t))/n$ which is found using condition (14) is transformed to the kinetic energy of particles. Neglecting an energy loss due to the particle radiation in direct reactions and surface radiation during the fireball expansion up to its breakup, from an approximate conservation of the energy we may evaluate the velocity of the collective flow. They are slightly overestimated because of ignoring the mentioned effects. Resulting values of the mean flow velocities at the freeze-out can be estimated as

[TABLE]

where ${\cal E}_{m},\,{\cal E}_{\rm b.up}$ are respectively the maximum energy per participant nucleon reachable for a given ${\cal E}_{\rm lab}$ and the energy per participant nucleon at the breakup. Within the KVORcut03$\Delta\pi_{\rm free}$ model we obtain $\bar{v}=0.19$, 0.29, 0.39, and 0.54 for collision energies ${\cal E}_{\rm lab}=0.18$, 0.4, 0.8, and 2.1$A$ GeV respectively and within the MKVOR*$\Delta\pi_{\rm free}$ model we get values 0.21, 0.299, 0.40, 0.53, which differ only slightly from those obtained in the KVORcut03$\Delta\pi_{\rm free}$ model. Note that the values $\bar{v}$ obtained for $\stackrel{{\scriptstyle\scriptstyle>}}{{\phantom{}{}_{\sim}}}1\,A$ GeV are probably too high and would be smaller, if we included effects of the $p$-wave pion-nucleon interaction, cf. Voskresensky:1993ud , see also results of a fit of the data for higher energies Adamczyk:2017iwn .

On the left panel in fig. 5 we show the initial baryon density and temperature, $n_{m},T_{m}$, as functions of the collision energy in the laboratory system within MKVOR*-based models (solid and dash-dotted lines) and KVORcut03-based models (dashed and dotted lines). We see that for all ${\cal E}_{\rm lab}$ the initial density and temperature for models without $\Delta$s (see dash-dotted and dotted lines) are larger than those with $\Delta$s (see solid and dashed lines). Also for all ${\cal E}_{\rm lab}$ the values of $n_{m}$ in the MKVOR*$\Delta\pi_{\rm free}$ model are higher than those in the KVORcut03$\Delta\pi_{\rm free}$ model. However, it is remarkable that the initial temperature dependence on ${\cal E}_{\rm lab}$ proves to be almost model independent within the same particle set. Thus we see that in our approach to choosing the initial state the dependence of the initial thermodynamic state on the model for the EoS resides in the value of the maximum reachable baryon density, while the maximum fireball temperature depends only weakly on the employed model. On the right panel in fig. 5 we show the mean velocity of the fireball expansion at the breakup, $\bar{v}$, as a function of the collision energy in the laboratory system evaluated within our models. We see that values $\bar{v}$ evaluated in KVORcut03$\Delta\pi_{\rm free}$ and MKVOR*$\Delta\pi_{\rm free}$ models prove to be approximately the same.

Taking into account the non-zero particle velocities at the breakup leads to a modification of their spectra. Experimental slopes of the spectra, which are determined by effective temperatures at freeze out, for nucleons are a bit higher than those for pions, cf. Nagamiya:1981sd . The mentioned difference is attributed to the fact that differential cross sections of massive nucleons are more affected by the presence of non-zero mean expansion velocity than differential cross sections of lighter pions, cf. Siemens:1978pb . Indeed, in a frame moving with the 3-velocity $\vec{v}$ the particle distribution is expressed through that in the rest frame by a replacement $p_{0}\to p_{\nu}u^{\nu}$, where $u^{\nu}$ is the 4-velocity of the frame. For non-relativistic particles and for $v\ll 1$ the transition to the moving reference frame is reduced to the replacement $(m+p^{2}/2m)/T\to(m+(\vec{p}-m\vec{v})^{2}/2mT$ and nucleon distributions are more affected than pion ones since $m_{N}v\gg m_{\pi}v$. In relativistic case in presence of the expansion velocity $\bar{v}$ the particle distributions are characterized by the effective temperatures $T_{{\rm ef},N}=T_{\rm b.up}\sqrt{1-\bar{v}^{2}}$ and by shifted momenta. Below we focus on pion distributions and determine the values $T_{\rm b.up}$ and $n_{\rm b.up}$ fitting the pion distributions. For $p\stackrel{{\scriptstyle\scriptstyle>}}{{\phantom{}{}_{\sim}}}m_{\pi}$ we may neglect $m_{\pi}v\stackrel{{\scriptstyle\scriptstyle<}}{{\phantom{}{}_{\sim}}}0.6\,m_{\pi}$ compared to $p$. Owing to this circumstance and taking into account that a slight decrease of $T_{{\rm ef},\pi}$ in comparison with $T_{\rm b.up}$ is partially compensated by the fact that $T_{\rm b.up}$ would be a bit higher, if we took into account an increase of the entropy in a realistic viscous expansion of the fireball, as in Voskresensky:1989sn ; Migdal:1990vm ; Voskresensky:1993ud , we simplifying put $T_{{\rm ef},\pi}\simeq T_{\rm b.up}$.

The momentum-dependent pion free path length proves to be short for pions with momenta $k\stackrel{{\scriptstyle\scriptstyle>}}{{\phantom{}{}_{\sim}}}1.5\,m_{\pi}$ up to the breakup stage and thereby such pions radiate from the fireball breakup Voskresensky:1991uv . Oppositely, pions with momenta $k<(1-1.5)\,m_{\pi}$ have larger mean-free path Voskresensky:1991uv and radiate from an intermediate stage of the fireball expansion, cf. Voskresensky:1993mw ; Voskresensky:1995wn . Moreover, a contribution to the pion yield comes from the decay of thermal $\Delta$-resonances at the breakup, or may be a bit later, if typical time of the reaction $\tau_{\Delta\to N\pi}$ at $n\sim n_{\rm b.up.}$, $T\sim T_{\rm b.up.}$ is larger than $\tau_{\rm b.up.}$ for thermal pions. Further we determine the values $n_{\rm{b.up}}({\cal{E}}_{\rm lab})$, $T_{\rm{b.up}}({\cal{E}}_{\rm lab})$ from the best fit of the differential pion cross sections for the momenta $k\stackrel{{\scriptstyle\scriptstyle>}}{{\phantom{}{}_{\sim}}}1.5m_{\pi}$.

3.2 Description of pion differential cross sections

The $\pi^{-}$ differential cross section in inclusive processes reads Senatorov:1989cg ; Voskresensky:1989sn

[TABLE]

where $\omega_{k}=\sqrt{m_{\pi}^{2}+{k}^{2}}$. The first term in the curly brackets is the contribution of thermal pions at the breakup stage, whereas the second term corresponds to the decay $\Delta\to N+\pi$ occurring at the breakup stage with direct radiation of free pions. For simplicity we consider collisions of nuclei with equal atomic weights $A=A_{\rm t}=A_{\rm p}$, $V_{{\rm b.up}}$ is the volume of the fireball at the breakup, $15.64\,A^{5/3}/n_{\rm{b.up}}$ is the geometric factor for inclusive processes. The numeric coefficient corresponds to all the quantities being measured in units of $m_{\pi}$, in particular $n_{\rm{b.up}}$ is measured in $m^{3}_{\pi}$, $n_{0}\simeq 0.45\,m_{\pi}^{3}$.

We assume that for $t>t_{0}$ during the breakup stage, which lasts for $\tau_{\rm b.up}(n_{\rm b.up},T_{\rm b.up})$, nucleons and pions decouple and for $t>t_{0}+\tau_{\rm b.up}$ pions, which before breakup stage were described by thermal distributions, can be considered as freely moving particles. In the quasiparticle model the distribution of pions is described by Senatorov:1989cg ; Voskresensky:1989sn

[TABLE]

The contribution from $\mu_{Q}\equiv\mu_{n}-\mu_{p}$ appears for $Y_{Z}\neq 0.5$ and distinguishes the $\pi^{-}$ from other pion species. The quantity $\omega(k)$ coincides with $\omega_{k}$, if one uses the vacuum dispersion law, and $\omega(k)$ is given by eq. (11) in the $\pi_{\rm WT}$ model. In the model of the sudden breakup (if the typical time for the pion sub-system breakup is $\tau_{\rm b.up}\ll 1/|\omega_{k}-\omega(k)|$) the value of $\Gamma_{k}$ is given by Senatorov:1989cg

[TABLE]

Using the spectrum (11) $\Gamma_{k}\simeq 1+O((n_{n}-n_{p})^{2}/n_{0}^{2})$ we can put $\Gamma_{k}\simeq 1$ at the breakup in both $\pi_{\rm WT}$ and $\pi_{\rm free}$ models. In the limit of a slow breakup, $\tau_{\rm b.up}\gg 1/|\omega_{k}-\omega(k)|$, one should use eq. (17) with $\omega(k)$ replaced by $\omega_{k}$. In present work where we disregard effects of the $p$-wave pion-baryon interaction the limit $\tau_{\rm b.up}\ll 1/|\omega_{k}-\omega(k)|$ looks as more realistic since effects of the $s$-wave interaction remain weak for $0.4\leq Y_{Z}\leq 0.5$ and thereby $1/|\omega_{k}-\omega(k)|\gg 1/m_{\pi}$, whereas a typical value of $\tau_{\rm b.up}$ is estimated as few $1/m_{\pi}$. However note that for $\omega(k)$ close to $\omega_{k}$ at breakup conditions the difference between pion distributions calculated in the limit of the prompt breakup and the slow breakup is not as significant.

For the contribution of $\Delta$ quasiparticles described with effective mass to the $\pi^{-}$ momentum distribution one gets Ivanov:2005yw

[TABLE]

Here $r$ runs through ($\Delta^{-},\,\Delta^{0}$). The branching ratios $b_{r}=(1,\,1/3)$ arise since $|N\pi>=\frac{1}{\sqrt{3}}|p\pi^{-}\rangle+\frac{2}{\sqrt{3}}|n\pi^{0}\rangle$ and therefore the rate of probabilities of reactions $\Delta^{0}\to n+\pi^{0}$, $\Delta^{0}\to p+\pi^{-}$ is $P(n\pi^{0})/P(p\pi^{-})=2$. The factor $g_{r}=4$ is the resonance spin degeneracy.

The in-medium effective masses should be used in (18), if the breakup time for pions to become freely moving particles and the time step at $n=n_{\rm b.up}$ and $T=T_{\rm b.up}$ for $\Delta$s to decay to pions and nucleons are shorter than the time typical for the nucleons to become freely moving particles. Otherwise during the breakup stage the in-medium $\Delta$s first transit to their vacuum branch and only after that decay into the freely moving nucleons and pions. In the latter case $m_{r}^{*}$ and $m_{N}^{*}$ in (18) should be replaced to $m_{r}$ and $m_{N}$.

Pion differential cross sections and yields in heavy-ion collisions at Bevalac and lowest GSI SIS energies were measured in Nagamiya:1981sd ; Miller:1987zz ; Hayashi:1988en ; Pelte:1997rg ; Sandoval:1980bm ; Harris:1987md ; Reisdorf:2006ie ; Hong:2005at with Ne, Ar, Ru, Zr, La, Au and Ru beams. In fig. 6 we show invariant differential cross sections of $\pi^{-}$ production in reactions Ne+NaF, Ar+KCl ($Y_{Z}\simeq 0.48$) and for La$+$La ($Y_{Z}\simeq 0.4$), computed within KVORcut03-(left panel) and MKVOR*-based (right panel) models with $\Delta\pi_{\rm free}$ particle set. Solid lines show results computed following eqs. (16-18). To be specific we use eq. (18) for $\Delta$s described with effective masses. The values $n_{\rm b.up},T_{\rm b.up}$ lie on the $s={\rm const}$ lines indicated on the right panels in figs. 3 and 4 and correspond to the best fit of the pion differential cross sections performed for $k>1.5m_{\pi}$. Dashed lines show the result obtained for the same $n_{\rm b.up},T_{\rm b.up}$ provided the contribution of $\Delta$ resonances is turned off. We see that $\Delta$ decays give a contribution centered around $E_{\pi}^{\rm kin}=(100-200)$ MeV, which is significant even for low temperatures due to the in-medium mean fields acting on the baryons. As expected, the cross sections for low pion kinetic energies $E_{\pi}^{\rm kin}\stackrel{{\scriptstyle\scriptstyle<}}{{\phantom{}{}_{\sim}}}(1-1.5)m_{\pi}$ are underestimated because the direct emission of low-momentum pions and the $p$-wave polarization effects are not taken into account. The differential cross sections computed for La $+$ La reactions for KVORcut03$\Delta\pi_{\rm WT}$ and MKVOR*$\Delta\pi_{\rm WT}$ models (i.e. using the dispersion law (11) and that $\mu_{{\pi}^{-}}=\mu_{n}-\mu_{p}$) are not shown since the difference with the distributions computed using the corresponding $\pi_{\rm free}$ models is not seen visually. More specifically, for ${\cal{E}}_{\rm lab}=800A$ MeV the difference changes within $(2-7)\%$, for lower collision energies it increases, and for ${\cal{E}}_{\rm lab}=246A$ MeV the difference reaches $16\%$, although on the logarithmic plot such a difference is almost not seen. So we conclude that our models with taking into account the $\Delta$ decays allow to describe the pion differential cross sections rather appropriately, even within $\pi_{\rm free}$ models, except for soft pions at low collision energies. Here we once more point out that we did not include a contribution Voskresensky:1993mw ; Voskresensky:1995wn of soft pions produced in direct reactions before the breakup and the $p$-wave polarization effects Voskresensky:1989sn ; Migdal:1990vm ; Voskresensky:1993ud , which inclusion should improve the description.

3.3 Description of ratios of pion yields to those of positively charged baryons

The mean number of protons per baryon measured after an event is given by

[TABLE]

which includes the contribution from the $\Delta$ decays with corresponding branching factors. Here we assume that all $\Delta$s at the breakup stage are transformed to pions and nucleons. The ratio of $\pi^{-}$-like to $p$ -like particle concentrations at the breakup is as follows

[TABLE]

where the first term in the brackets shows the contribution of thermal pions computed with the help of eq. (17), which at breakup stage are transformed to freely moving pions, and the second term given by (18) takes into account the decay of the thermal $\Delta$s into the freely moving pions and the nucleons at the breakup stage. As in eq. (16), $r$ runs through $(\Delta^{-},\Delta^{0})$ and the respective values $b_{r}$ are $(1,1/3)$. In our model, where we disregard the surface effects, the $R_{\pi^{-}Z}$ ratio does not depend on the fireball volume. In a case of the ISM we get $R_{\pi^{-}Z}=2(n_{\pi^{-}}+\frac{1}{3}n_{\Delta})/n$ and in the Boltzmann limit we have approximately:

[TABLE]

where $\Delta m=m^{*}_{\Delta}(n_{\rm{b.up}},T_{\rm{b.up}})-m^{*}_{N}(n_{\rm{b.up}},T_{\rm{b.up}})$ when we use eq. (18) with in-medium baryon masses, and $\Delta m=m_{\Delta}-m_{N}$, if we explore (18) with vacuum baryon masses.

In fig. 7 we show the results for the breakup density $n_{\rm b.up}$ (upper panels) and $T_{\rm b.up}$ (lower panels) within our models. The symbols denote the best fits using (19) to the available data on $R_{\pi^{-}Z}$ in Ne+NaF and Ar+KCl collisions, and the shaded areas demonstrate the associated uncertainties within the data error bars. The triangles and a star show the corresponding results for the differential cross-section given by eq. (16). We see that at ${\cal E}_{\rm lab}\stackrel{{\scriptstyle\scriptstyle>}}{{\phantom{}{}_{\sim}}}800$ MeV for each model both fits for Ne+NaF and Ar+KCl generally agree with each other within the error bars. For ${\cal E}_{\rm lab}\stackrel{{\scriptstyle\scriptstyle<}}{{\phantom{}{}_{\sim}}}800$ MeV the breakup densities and corresponding temperatures required by $R_{\pi^{-}Z}$ ratio are larger than those following from the fit of the differential cross sections demonstrated in fig. 6. To parameterize a possible breakup density dependence on ${\cal E}_{\rm lab}$ for the energies $200\,{\rm MeV}\stackrel{{\scriptstyle\scriptstyle<}}{{\phantom{}{}_{\sim}}}{\cal E}_{\rm lab}\stackrel{{\scriptstyle\scriptstyle<}}{{\phantom{}{}_{\sim}}}2100\,{\rm MeV}$ of our interest here we use the interpolation formula

[TABLE]

where $\theta(x)$ is the step function. Lines given by (21) pass through both the regions for Ne+NaF and Ar+KCl, and are shown for each model by solid lines in fig. 7.

In fig. 8 we show the fitted ratios of $\pi^{-}$-like yields to those of $p$-like yields, $R_{\pi^{-}Z}$, for particles emitted at the breakup of the nuclear fireball within the KVORcut03$\Delta\pi_{\rm free}$ model (left panel) and MKVOR*$\Delta\pi_{\rm free}$ model (right panel). Solid curves are calculated according to eq. (19) assuming $Y_{Z}=0.48$ with the $\Delta$ contribution calculated using the prompt breakup assumption, so all the $\Delta$s with effective masses $m^{*}_{\Delta}(n_{\rm b.up},T_{\rm b.up})$ are supposed to be transformed at the fireball breakup to free pions and nucleons with masses $m^{*}_{N}(n_{\rm b.up},T_{\rm b.up})$, which after switching off the mean fields go on the mass shell. The breakup parameters follow the solid curves in fig. 7 given by eq. (21), and we see that the $R_{\pi^{-}Z}$ data are indeed well described for such breakup parameters. Thin dashed lines show the $R_{\pi^{-}Z}$ ratios calculated taking into account only the thermal pion contribution to $R_{\pi^{-}Z}$ (first term in (20)) for the same breakup parameters. We see that $\Delta$s contribute essentially in $R_{\pi^{-}Z}$ ratios in the whole interval of energies $200\,{\rm MeV}\stackrel{{\scriptstyle\scriptstyle<}}{{\phantom{}{}_{\sim}}}{\cal E}_{\rm lab}\stackrel{{\scriptstyle\scriptstyle<}}{{\phantom{}{}_{\sim}}}2100\,{\rm MeV}$ demonstrated in figures. As in fig. 7, the triangles show the results with the breakup parameters following the best fits of the spectra. Data points are extracted from refs. Nagamiya:1981sd ; Sandoval:1980bm ; Harris:1987md . These data show that the pion yields are approximately linear functions of the participant number and that $R_{\pi^{-}Z}$ are approximately the same for different colliding nuclei with the same $Y_{Z}$. We see that the best fit values $(n_{\rm b.up},T_{\rm b.up})$ extracted from the $\pi^{-}$ spectra noticeably underestimate the $R_{\pi^{-}Z}$ ratios compared with the experimental data.

Figures 7 and 8 demonstrate that there is substantial difference in the breakup parameters inferred from $R_{\pi^{-}Z}$ ratios and the $\pi^{-}$ spectra. There are two reasons of this discrepancy within the current implementation of the thermodynamical model. It is known Voskresensky:1989sn ; Migdal:1990vm ; Voskresensky:1993ud that, if the in-medium modifications of the pion spectra due to the $p$-wave pion-baryon interaction are taken into account, the pion yields would increase, especially at low temperatures (for low ${\cal E}_{\rm lab}$). Another reason is a lack of soft pions with long mean-free paths in our calculations, being emitted in direct reactions before the fireball breakup. Their inclusion would not change the results of the fit of the spectra, because the fit is done for pions with kinetic energies $\stackrel{{\scriptstyle\scriptstyle>}}{{\phantom{}{}_{\sim}}}1\,m_{\pi}$. However, the direct pions will contribute substantially to $R_{\pi^{-}Z}$, which are largely underestimated, if the breakup parameters following from the $\pi^{-}$ spectra are used. As we expect, if the pion distributions and yields were calculated with taking into account mentioned effects, the values $(n_{\rm b.up},T_{\rm b.up})$ obtained in those fits would coincide.

3.4 Ratios of $\pi^{-}/\pi^{+}$ yields

Experimental data demonstrate that ratio of the $(\pi^{-}/\pi^{+})$ yields is larger than unity and increases with a decrease of the collision energy. In fig. 9 experimental ratios are shown by the squares. Previously $(\pi^{-}/\pi^{+})$ ratios have been computed within different models, cf. Xu:2013aza ; Hong:2013yva and refs. therein. A role of effects of $\Delta$ and nucleon optical potentials, $s$-wave and $p$-wave interactions has been discussed.

Since the $\Delta$ decays to nucleons and pions contribute to the ratio of the $(\pi^{-}/\pi^{+})$ we should compare the ratios with the data using the following relation

[TABLE]

As before, we assume that the breakup is prompt and the $\Delta$ densities are calculated with the in-medium effective masses and take into account the finite $\mu_{Q}>0$ for $Y_{Z}\neq 0.5$.

The solid curves in fig. 9 show results of our calculations performed using eqs. (22) for $Y_{Z}=0.4$. We use the same $n_{\rm b.up}({\cal E}_{\rm lab})$ given by (21) which we used to describe $R_{\pi^{-}Z}$ ratios. Triangles denote the results for $(n_{\rm b.up},T_{\rm b.up})$ obtained by fitting the differential cross sections shown in fig. 6. We employ the $\pi_{\rm free}$ model. If we used the WT model with spectra following eq. (11), the $s$-wave pion-nucleon interaction would affect the $\left(\frac{\pi^{-}}{\pi^{+}}\right)_{\rm like}$ ratios only slightly, since $e^{-(n_{n}-n_{p})/2f_{\pi}^{2}T}\simeq 1$ in the collision energy interval we show in fig. 9. The dashed lines demonstrate the results obtained within the models with the $\Delta$ contribution artificially suppressed. We see that without $\Delta$s the $(\pi^{-}/\pi^{+})$ ratios prove to be significantly overestimated. Results of calculations in both of our models KVORcut03$\Delta\pi_{\rm free}$ and MKVOR*$\Delta\pi_{\rm free}$ reproduce the experimentally observed increase of the $(\pi^{-}/\pi^{+})_{\rm like}$ ratios with a decrease of ${\cal E}_{\rm lab}$. However, the results of our calculations are roughly by $30\%$ higher than the experimental data. This can be again attributed to a lack of directly emitted soft pions and to the lack of $p$-wave pion-baryon interaction effects in our calculations.

4 First-order liquid-gas phase transition in KVORcut03 and MKVOR*- based models

The first-order nuclear LG phase transition in the matter prepared in low energy heavy-ion collisions may occur for densities $n<n_{0}$ and temperatures $T<T_{cr}^{\rm LG}<(15-23)$ MeV Ropke:1982vzx , see Dutra:2018amp and refs. therein. The values of the critical density and temperature depend on the model of the EoS. The $\Delta$ isobars do not contribute at such low densities and temperatures, and contribution of all thermal excitations of mesons (pions, $\sigma$, $\omega$, $\rho$ mesons and heavier ones) is negligibly small. The temperature dependence remains essential only in nucleon terms. Although we continue to use KVORcut03$\Delta\pi_{\rm free}$ and MKVOR*$\Delta\pi_{\rm free}$ models the results for densities $n<n_{0}$ and for temperatures $T\stackrel{{\scriptstyle\scriptstyle<}}{{\phantom{}{}_{\sim}}}(20-30)$ MeV are the same as within KVORcut03 and MKVOR models.

A general physical picture of the dynamics of the nuclear first-order LG phase transition has been formulated long ago, cf. Schulz:1983pz , where the EoS of ISM was considered in the original Walecka RMF model. The non-ideal hydrodynamics of the nuclear LG phase transition was then studied in Skokov:2008zp ; Skokov:2009yu ; Skokov:2010dd ; Voskresensky:2010gf on an example of the van der Waals EoS. The case of the ISM has been studied first. Depending on the collision energy, the system may enter the region of the supercooled gas, overheated liquid or the region of the spinodal instability. Thus depending on the experimental conditions, the system may undergo the hydrodynamical instability, the first-order phase transition from the supercooled gas to the liquid, or the transition from the overheated liquid to the gas. For ISM one deals with one conserved baryon charge. In case of IAM one deals with two conserved charges, the baryon charge and the isospin, that adds specific features to the problem, cf. Muller:1995ji ; Ducoin:2005aa ; Alam:2017krb .

4.1 Liquid-gas phase transition in ISM

4.1.1 RMF consideration

First let us consider the LG phase transition in ISM within our models of EoS. The pressure isotherms as functions of the density in the region of the first-order LG phase transition are demonstrated in fig. 10 for the models KVORcut03 (upper panel) and for MKVOR* (lower panel). The Gibbs conditions determining equilibrium of two phases (I-liquid and II-gas) are as follows:

[TABLE]

where $\mu_{B}=\mu_{n}=\mu_{p}$, $n=n_{n}+n_{p}$, $Y_{p}=1/2$ and thereby $n=2n_{n}$. The phase separation boundary in fig. 10 is denoted by a solid line. This line connects the endpoints of the isothermal Maxwell constructions. The isothermal spinodal (ITS) line, which delimits the hydrodynamically unstable region, is denoted by a dashed line. On the Maxwell construction (shown for an example by a horizontal dashed line for $T=10$ MeV) the total baryon density is composed from the densities of the liquid and gas phases as

[TABLE]

where the liquid phase fraction $\chi$ changes along the Maxwell construction horizontal line from [math] at $n=n^{\rm II}$ up to 1 for $n=n^{\rm I}$. The location of the critical point, which satisfies the conditions $(\partial P/\partial n)_{T}=0$, $(\partial^{2}P/\partial n^{2})_{T}=0$, is described by the critical values of $T=T_{cr}$, $n=n_{cr}$, $s={s}_{cr}$. For KVORcut03 model $T_{cr}=17.4\,{\rm MeV}$, $n_{cr}=0.054\,{\rm fm}^{-3}$, ${s}_{cr}=2.74$, $P_{\rm cr}=0.30\,{\rm MeV}/{\rm fm}^{3}$, and for MKVOR* model $T_{cr}=16.05\,{\rm MeV}$, $n_{cr}=0.051\,{\rm fm}^{-3}$, ${s}_{cr}=2.68$, $P_{\rm cr}=0.25\,{\rm MeV}/{\rm fm}^{3}$. These values lie within the $T_{\rm cr}=(16.4\pm 2.3)\,{\rm MeV}$ band given by microscopic calculations based on chiral nucleon-nucleon potentials Carbone:2018kji . The shaded region denotes the experimental bound taken from Elliott:2013pna : $n_{\rm cr}=(0.06\pm 0.02)\,{\rm fm^{-3}}$, and $P_{\rm cr}=(0.3\pm 0.1)\,{\rm MeV}/{\rm fm}^{3}$, $T_{\rm cr}=(17.9\pm 0.4)\,{\rm MeV}$. Both the KVORcut03 and MKVOR* models pass the constraints for critical density and pressure and KVORcut03 model passes marginally the constraint for critical temperature from Elliott:2013pna . The critical temperature calculated in the MKVOR* model is lower than the experimental bound extracted in Elliott:2013pna . Note that the traditional RMF models do not fulfill the existing experimental constraints, cf. fig. 2 in Dutra:2018amp . Two dotted lines demonstrate trajectories of the system characterized by the values of the constant specific entropy and are labeled by the values of $s$. The boundary of the hydrodynamically unstable region, as it occurs at constant value of the specific entropy – adiabatic spinodal (AS), is shown by a dash-dotted line.

In fig. 11 we show the phase diagram on $T(n)$ plane for ISM within the KVORcut03 (upper panel) and MKVOR* (lower panel) models. By solid and dotted curves we show boundaries of the Maxwell construction and the ITS regions, respectively. The dash-dotted curves indicate boundaries of the AS region. Dashed lines labeled by values of the laboratory collision energy show adiabatic trajectories of the system evolution computed within the approach described in the previous section. Dash-dotted lines ending with full dots show values of the temperature $T_{\rm comp}(n_{\rm comp})$, at which the pressure reaches zero at the maximum $n$ among three roots on each isotherm. For $T<T_{\rm comp}$ a compound nucleus can be formed, which by radiation of nucleons may reach the ground state. Shaded rectangle denotes the constraint from the experimental analysis of reactions with formation of a compound thermal nucleus, which emits neutrons, protons and heavier charged fragments Elliott:2013pna , and the hatched band denotes the result of calculations performed with chiral forces Carbone:2018kji .

After the ${s}\simeq{\rm const}$ trajectory of the expanding system crosses the border of the Maxwell construction, the system comes into a metastable region. One may expect that the typical fireball expansion time is less than the typical time for which in the system, being in the supercooled gas (overcooled liquid) state, an overcritical liquid drop (gaseous bubble) appears. At such assumptions we can still apply the thermodynamical description until the system trajectory did not reach the ITS line. Beyond the ITS line the system enters the hydrodynamically unstable region. The $s\simeq{\rm const}$ system trajectory crosses the ITS line already for the collision energies $E_{\rm lab}\stackrel{{\scriptstyle\scriptstyle<}}{{\phantom{}{}_{\sim}}}250-300$ MeV within our models, however the system breakup occurs at higher temperatures than $T_{cr}$, as we have seen from figs. 3, 4. Thereby, within our models the system breakup occurs in the region of the spinodal instability only for $E_{\rm lab}\stackrel{{\scriptstyle\scriptstyle<}}{{\phantom{}{}_{\sim}}}150$ MeV. For $E_{\rm lab}\stackrel{{\scriptstyle\scriptstyle<}}{{\phantom{}{}_{\sim}}}100$ MeV the system trajectory passes the region of $T-n$, where the pressure reaches zero and a compound nucleus can be formed. Experimentally such a region can be attained not only in heavy-ion collisions but also in collisions of light nuclei and particles with a heavier nucleus, cf. Elliott:2013pna .

The isothermal and adiabatic sound velocities, $u_{T}$ (at ${T}={\rm const}$) and $u_{{s}}$ (at ${s}={\rm const}$), satisfy the thermodynamic relation, cf. Skokov:2010dd ,

[TABLE]

Here $u_{{s}}^{2}=(\partial P/\partial E)_{{s}}$ and $u_{T}^{2}=(\partial P/\partial E)_{T}$, $c_{V}$ is the specific heat at fixed volume $V$. The quantity $u_{{s}}$ characterizes propagation of sound waves in ideal hydrodynamics. In non-ideal hydrodynamics with finite values of the shear and bulk viscosities, $\eta$ and $\zeta$, and thermal conductivity, $\kappa$, the propagation of sound waves is defined by the interplay between $u_{T}$ and $u_{{s}}$, cf. Skokov:2010dd . Condition $u_{T}=0$ defines the ITS line on the $T-n$ plane and $u_{{s}}=0$ determines the AS line. The maximum temperature points on these lines are the critical temperature $T_{cr}$ (for the $P(n)$ isotherms $T={\rm const}$), and the adiabatic maximum temperature $T_{P,max}$ (for $P(n)$ at ${s}={\rm const}$).

If the evolution is studied within the ideal hydrodynamics (for ${s}={\rm const}$), then the instability occurs at the crossing of the AS line, where $u_{{s}}=0$. Note that $u_{T}^{2}<0$ at the point $u_{{s}}=0$, whereas $u_{T}^{2}$ reaches zero first, when the system trajectory crosses the ITS line. Note that BS ; Panagiotou:1984rb and a number of subsequent works exploiting ideal hydrodynamics considered the crossing of the AS line by the trajectories ${s}$=const as the starting point for the clustering. In contrast, authors  Schulz:1983pz and later  FRS ; FRS1 ; PR ; LL ; Skokov:2009yu ; Skokov:2010dd argued that clustering in realistic systems, described by non-ideal hydrodynamics, may appear already at higher temperatures after the system trajectory crosses the ITS line.

4.1.2 Effects of fluctuations in ISM

Throughout this work we use the RMF framework. However, an inclusion of fluctuations Skokov:2009yu ; Skokov:2010dd results in the divergence of the specific heat at $T=T_{cr}$. Thereby it can be concluded from eq. (24), that with taking fluctuations into consideration both temperatures $T_{cr}$ and $T_{P,max}$ should coincide, cf. FRS ; FRS1 .

The heat capacity $C_{V}$ is related to the variance of the energy as

[TABLE]

within the grand canonical description. Papers  Stephanov:1999zu suggested that, if at some incident energy the trajectory passes in the vicinity of the critical point, the system may linger longer in this region due to strong thermodynamical fluctuations resulting in the divergence of the susceptibilities. For example, one could expect that in heavy-ion reactions the vicinity of the critical point could manifest via abnormal fluctuations of the energy. However ref. Skokov:2009yu argued that fluctuation effects in the vicinity of the critical point can hardly be pronounced in actual heavy-ion collisions, since the typical time for the development of critical fluctuations diverges in the critical point. Thereby, all relevant processes prove to be frozen near the critical point, whereas the system passes this region in a rather short time. Concluding this discussion, the time which the system spends in a vicinity of the critical point in the course of the heavy ion collisions might be insufficient for the formation of critical fluctuations.

Now we argue that the ITS line can be manifested in abnormal fluctuations of the conserved charges, such as the baryon number. In equilibrium one-component systems having the mean density $n=\overline{\hat{n}}({\bf r})$, the static structure factor describing fluctuations of the conserved particle number can be written in terms of the density derivative of the chemical potential at fixed temperature LL1980 ; Roepke:2017bad ,

[TABLE]

Thereby the normalized variance of the nucleon number $N$ for the ISM (in absence of other baryons except nucleons) is given by,

[TABLE]

Since $\left(\frac{\partial P}{\partial n}\right)_{T}=n\left(\frac{\partial\mu_{b}}{\partial n}\right)_{T}=0$ on the whole ITS line $T(n)$, the static structure factor and the variance of the baryon number diverge on this line. Thus fluctuations of the baryon number grow when the system trajectory comes close to the ITS line, which can be seen in corresponding event-by-event observables.

In fig. 12 for KVORcut03 model (upper panel) and MKVOR* model (lower panel) we show the dependence of $w(T)$ in ISM that we computed numerically for several values of $n$ in the region of the LG first-order phase transition. As we see, the normalized variance diverges on the ITS line, cf. fig. 11. The value $n=0.8n_{0}$ is outside the ITS region, thereby $w$ for $n=0.8n_{0}$ shown by the dotted line remains a smooth function of $T$.

4.2 Liquid-gas phase transition in IAM

Recently two reactions 124Xe $+^{112}$Sn and 136Xe $+^{112}$Sn have been experimentally studied at 32$A$ MeV and 45$A$ MeV bombarding energies to produce quasi-fusion hot nuclei, which undergo multifragmentation Borderie:2018fsi . Using charge correlations the fossil signature of spinodal instabilities, i.e. the abnormal presence of equal-sized fragments was established at a confidence level of around 6-7 sigma for both reactions at 32$A$ MeV incident energy.

In low-energy heavy-ion collisions, the baryon and electric charges obey two independent conservation laws, implying that the proton fraction $Y_{p}\simeq Y_{Z}$ is conserved. Therefore for the description of the LG first-order phase transition in the IAM formed in heavy-ion collisions we introduce two chemical potentials of the baryon charge $\mu_{B}=\mu_{n}$ and the electric charge $\mu_{Q}=\mu_{n}-\mu_{p}$. First, following the analysis performed in the literature in different models, consider an occurrence of the first-order LG phase transition at assumption that there is no surface tension on the spatial border of phases, cf. Ducoin:2005aa ; Alam:2017krb ; Glendenning:2001pe ; Poberezhnyuk:2018mwt .

As in fig. 10, in fig. 13 we show the $P(n)$ isotherms, but now for $Y_{p}=0.4$ within KVORcut03 EoS (upper panel) and MKVOR* EoS (lower panel). $P(n)$ isotherms are shown by solid lines. Numbers near the curves indicate values of the temperature in MeV. Contrary to the $Y_{Z}=0.5$ case, the equilibrium states within mixed phase are not described by horizontal Maxwell lines, but by the curves shown in fig. 13 by dashed lines. This is a general feature of the description of the first-order phase transitions in multi-component systems Glendenning:2001pe . Transitions in systems with more than one conserved charge are offen called Gibbs phase transitions or non-congruent phase transitions Iosilevskiy:2010qr ; Hempel:2013tfa .

In fig. 14 by solid lines we show the phase coexistence borders in the proton-neutron density plane for KVORcut03 (upper panel) and MKVOR* (lower panel) models. Bold dots show critical points. Numbers near the curves are the corresponding values of $T$ in MeV. Dashed lines denote borders of the isothermal spinodal instability regions, defined as the line where the matrix

[TABLE]

becomes singular. This corresponds to the lowest eigenvalue becoming negative, which is the known condition for the presence of the instability Chomaz:2003dz . We see that the allowed region of the LG phase transition on the proton-neutron density plane is a bit broader in MKVOR* model than in KVORcut03 model, whereas the critical temperature is a bit higher in the latter model.

On the pressure $P(n)$ isotherm for $T<T_{cr}$ the mixed phase begins at the point $G_{\rm eq}$ and ends at the point $L_{\rm eq}$. The point $G_{\rm eq}$ corresponds to the gaseous phase being in equilibrium with an infinitesimal fraction of the liquid phase $L_{0}$, which has a higher density at the same pressure and both chemical potentials:

[TABLE]

The point $L_{\rm eq}$ corresponds to the liquid phase, being at equilibrium with the gaseous phase $G_{0}$ of infinitely small fraction and a lower density at the same pressure and chemical potentials:

[TABLE]

At the points ${\rm G_{\rm eq}}$ and ${\rm L_{\rm eq}}$ the values of the pressure, as well as of both chemical potentials, are different, cf. similar discussion in Ducoin:2005aa . The baryon and electric charge densities in the mixed phase are connected to those in each phase as

[TABLE]

[TABLE]

where $n^{\rm G}$, $n_{Q}^{\rm G}$ and $n^{\rm L}$, $n_{Q}^{\rm L}$ are the baryon and charge densities of the gaseous and the liquid fractions, respectively, the liquid phase fraction is $0<\chi<1$, $Y_{p}\simeq n_{Q}/n$. For a fixed proton fraction $Y_{p}={\rm const}$ the system in the phase coexistence region follows the critical line, that allows us to determine the fraction of the liquid phase in the coexistence region as

[TABLE]

In fig. 15 we demonstrate the description of the LG phase transition in IAM in terms of $\mu_{p},\mu_{n}$ (left) and $n_{p},n_{n}$ for comparison (right). As an example, we present the results of calculations performed within the MKVOR* model at $T=10$ MeV and $Y_{p}=0.3$ for easier comparison with the result of calculation Ducoin:2005aa presented there in fig. 8 for the same values of $T$ and $Y_{p}$. Inside the coexistence region, the system with a given proton fraction $Y_{p}=0.3$ is decomposed into two phases, located at the intersections of the coexistence curve with the corresponding isotherm (left panel) and with constant $\mu_{n}-m_{N}$ curve (right panel). The constant $Y_{p}$ transformation does not exhibit a transition from liquid to gas at a single value of $\mu_{n}$. The intensive parameters change smoothly as the system is driven along the coexistence line with an increasing density. We can see on the right panel in fig. 15 that the liquid phase of a neutron-rich matter is closer to ISM than the gaseous phase, so our models exhibit the known isospin-distillation phenomenon Muller:1995ji ; Colonna:2002ti ; Ducoin:2005aa . When $Y_{p}$ is kept constant, the system is forced to follow the first-order phase-transition line Ducoin:2005aa . When the system reaches the coexistence border (point $G_{\rm eq}$), a liquid phase of infinitesimal fraction appears in point $L_{0}$ at the same values of $\mu_{n}$, $\mu_{p}$ and $T$. The gaseous phase goes along coexistence line from the point $G_{\rm eq}$ to point $G_{0}$, as it is indicated on right panel in fig. 15, while the dense phase goes on the other side of the coexistence border from $L_{\rm eq}$ to $L_{0}$. When the system reaches the state $L_{\rm eq}$ the gas is entirely transformed into a liquid, the phase transition from gas to liquid is over and the $Y_{Z}={\rm const}$ transformation corresponds to a homogeneous system again.

4.2.1 Effects of fluctuations in IAM

For the IAM the static structure factor can be presented as Burrows:1998cg ; Roepke:2017bad

[TABLE]

provided that $S_{np}=S_{pn}$. Partial normalized variances are

[TABLE]

where $n_{\tau}=\frac{\partial P[\mu_{\tau},\mu_{\tau^{\prime}},T]}{\partial\mu_{\tau}}$ and the matrix $C_{\tau\tau^{\prime}}$ is defined by (27). These expressions for $w_{\tau\tau^{\prime}}$ become manifestly equivalent if all $\partial n_{\tau}/\partial\mu_{\tau^{\prime}}$ are evaluated with all other $\mu_{i},\,i\neq\tau^{\prime}$ being fixed, and in turn $\partial\mu_{\tau}/\partial n_{\tau^{\prime}}$ are evaluated with all other $n_{i},i\neq\tau^{\prime}$ being held constant. Similarly we may introduce normalized variances of the baryon number

[TABLE]

and the charge

[TABLE]

In fig. 16 we show our results for $w_{B}$ (bold lines) and $w_{Q}$ (thin lines) as functions of the temperature for several values of the baryon density $n$ within the PT region. We see that both $w_{B}$ and $w_{Q}$ are divergent at the border of the isothermal spinodal instability region, because the matrix $C_{ij}$ becomes singular there.

4.3 Features of pasta phase in systems formed in heavy-ion

collisions

The conservation of the electric charge should be considered not locally, but globally Glendenning:2001pe . In the previous subsection following the standard description of the LG phase transition in IAM in finite-size systems in absence of the charge neutrality, cf. Ducoin:2005aa , effects of a non-zero surface tension $\sigma$ on the LG phase boundary were disregarded. Such a description assumes that $\sigma=0$. In reality $\sigma(n\sim n_{0},T=0)\sim 1$ MeV$/$fm2, cf. Baym:1971ax . In presence of $\sigma\neq 0$ the uniform BEM can be structured Voskresensky:2001jq . We will show that taking into account the surface tension between liquid and gas structures significantly modifies the physical picture of the LG transition phenomenon in isospin-asymmetric systems of finite size (here, in absence of the charge neutrality of the system).

Typical radius of the fireball is

[TABLE]

for a spherical expansion, and thereby $R_{\rm f}(t)\stackrel{{\scriptstyle\scriptstyle<}}{{\phantom{}{}_{\sim}}}(15-30)$ fm for central collisions of the heaviest atomic nuclei, and for $n\sim(0.1-0.5)n_{0}$. We will for simplicity assume that before the first-order LG phase transition occurs the fireball is homogeneously charged with the proton charged density $n_{Q}$. In order for the pasta phase to be formed it should be at least $R_{\rm WZ}<R_{\rm f}$, where $R_{\rm WZ}$ is the radius of the Wigner-Seitz cell for the given $d$ geometry, $d=3$ for droplets/bubbles, $d=2$ for rods and $d=1$ for slabs. The droplet fills the sphere $r_{3}=\sqrt{x^{2}+y^{2}+z^{2}}<R$, the rod fills space of the cylinder $r_{2}=\sqrt{x^{2}+y^{2}}<R$, where $R$ is now the transversal radius of the rod, and the slab fills the layer $r_{1}=|x|<R$. We assume $R<R_{\rm WZ}<R_{\rm f}$. Moreover, we will imply that the fireball expansion is so slow that there is enough time to prepare the pasta structures, which we further on consider in the static approximation. Thus, simplifying consideration we assume that all mentioned conditions are fulfilled although in actual heavy-ion collisions it might be not the case.

After the pasta phase has appeared, the charge density is redistributed as

[TABLE]

$n_{Q}=nY_{p}$, $\delta n_{Q}=\delta n_{Q}^{\rm I}$ for $r_{d}<R$, and $\delta n_{Q}=\delta n_{Q}^{\rm II}$ for $R<r_{d}<R_{\rm WZ}$, where we for specificity assumed that the interior of a structure is in the liquid phase and the exterior is in the gaseous phase. Thereby we start with consideration of the denser phase as the minor phase (for $0<\chi\stackrel{{\scriptstyle\scriptstyle<}}{{\phantom{}{}_{\sim}}}0.5$).

From eq. (29) we find the relation

[TABLE]

that is similar to the global charge neutrality condition in the pasta phase realized in BEM of neutron stars.

Then we present the electric potential well as $V=V_{0}+\delta V$, with $V_{0}$ and $\delta V$ obeying the Poisson equations

[TABLE]

$e>0$ is the proton charge. Simplifying, we assume the quantities $\delta n_{Q}^{\rm I,II}$ to be spatially uniform in each region. With the solution of eq. (34) at the condition (33) we are able to find the electric field contribution to the free-energy density (for fixed $T$) due to appearance of the Coulomb pasta Pethick:1983eft ; Heiselberg:1992dx , cf. also Norsen:2000wb ; Voskresensky:2002hu ,

[TABLE]

where

[TABLE]

and thereby $\Phi_{3}(\chi)=\frac{2-3\chi^{1/3}+\chi}{5}$, $\Phi_{2}(\chi)=\frac{\ln(1/\chi)+\chi-1}{4}$, and $\Phi_{1}(\chi)=\frac{(1-\chi)^{2}}{3\chi}$.

The surface free-energy density contribution is

[TABLE]

For $T\neq 0$ the surface tension $\sigma$ is given by Ravenhall:1984ss ; Schulz:1983pz , $\sigma=\sigma_{0}\left(\frac{T_{c}^{2}-T^{2}}{T_{c}^{2}+T^{2}}\right)^{5/4}$, $\sigma_{0}\sim 1$ MeV$/$fm2.

The total (Coulomb plus surface) contribution to the free-energy density due to Coulomb pasta effects $\delta F_{\rm C.pasta}$ is given by

[TABLE]

Now we are able to minimize $\delta F_{\rm C.pasta}^{d}(R)$ in $R$. Thus using the condition (33) we find the optimal size of a structure with a given geometry parameter $d$ and $\delta F_{\rm C.pasta}^{d}(R_{m}^{d})$,

[TABLE]

A simple perturbative way to take into account finite size effects is as follows Glendenning:2001pe ; Voskresensky:2001jq ; Voskresensky:2002hu . First for the given mean values $n$, $n_{Q}=nY_{p}$ at fixed $T<T_{cr}$ following eqs. (28), (29) we find values $n^{\rm G}$, $n^{\rm G}_{Q}=n_{Q}+\delta n^{\rm II}_{Q}$, $n^{\rm L}$, $n^{\rm L}_{Q}=n_{Q}+\delta n^{\rm I}_{Q}$, $\chi$, without including finite size contribution. Then, implying that the correction to the free-energy density owing to the finite-size term $\delta F_{\rm C.pasta}^{d}(R_{m}^{d})$ is small we add this term not modifying $\chi$.

Comparing $\delta F_{\rm C.pasta}^{d}(R_{m}^{d})$ for $d=1,2,3$ at a given $\chi$ and $\delta n_{Q}^{\rm I}$ we determine the energetically favorable geometry of the structures. The following structures are energetically favorable:

[TABLE]

The result for $\chi>0.5$ is obtained by the replacement $\chi\to 1-\chi$ and describes a bubble phase. Note Voskresensky:2001jq ; Voskresensky:2002hu that the Coulomb limit is actually valid only for $R_{m}^{d}\ll\lambda_{\rm D}^{\rm I,II}$, where $1/\lambda_{\rm D,I}^{2}=4\pi e^{2}\frac{\partial n_{Q}^{\rm I}}{\partial\mu_{Q}}$, $1/\lambda_{\rm D,II}^{2}=4\pi e^{2}\frac{\partial n_{Q}^{\rm II}}{\partial\mu_{Q}}$ are taken for $V=0$. Otherwise a spatial inhomogeneity of the charge distribution in both liquid and gas phases should be taken into account. Complete numerical calculations can be performed similar to those Maruyama:2005vb ; Maruyama:2005tb done for BEM in neutron stars.

For $Y_{p}=1/2$ the condition (29) coincides with (28). In this case $Y_{p}^{\rm I}=Y_{p}^{\rm II}=1/2$ and eq. (33) is fulfilled for $\delta n_{Q}^{\rm I}=\delta n_{Q}^{\rm I}=0$. Solutions with $\delta n_{Q}^{\rm I,II}\neq 0$ are energetically not profitable due to the symmetry energy. For the given mean density $n^{\rm II}<n<n^{\rm I}$ one finds the fraction $\chi$, now $\chi=(R/R_{\rm f})^{d}$, $R_{\rm WZ}=R_{\rm f}$. Thus one defines $R(\chi(n))$. Therefore for $Y_{p}=1/2$ there is no Coulomb pasta. The situation is rather similar to that occurring in the description of BEM in neutron stars by the Maxwell construction. Then there is only one boundary (following the Maxwell construction) that separates the phase I from the phase II. The minor phase (I for $\chi<1/2$) occupies the drop, the rod or the slab provided $\delta F_{\rm Cpasta}^{d}(R(\chi(n)))$ is the smallest among the drop, rod or slab geometries, respectively.

For $Y_{p}\neq 1/2$ for $R<R_{\rm WZ}<R_{\rm f}$ the pasta phase is energetically preferable compared to the case of the only one boundary (for $R_{\rm WZ}=R_{\rm f})$ since $\delta F_{\rm Cpasta}^{d}(R_{m}^{d})$ is less than that for $R\neq R_{m}$ (at least provided $R\ll R_{\rm f}$). It would be interesting to look for possible consequences of the formation of the pasta structures in actual heavy-ion reactions.

5 Conclusion

In this paper we constructed the equation of state (EoS) of the nuclear matter within the relativistic mean-field (RMF) framework with hadron masses and coupling constants dependent on the mean scalar field. We considered a range of thermodynamic parameters relevant for description of heavy-ion collisions at the laboratory energies per baryon ${\cal E}_{\rm lab}\stackrel{{\scriptstyle\scriptstyle<}}{{\phantom{}{}_{\sim}}}(1-2)A$ GeV, namely temperatures $T$ below $m_{\pi}$, densities $n\stackrel{{\scriptstyle\scriptstyle<}}{{\phantom{}{}_{\sim}}}3n_{0}$, and the isospin asymmetry in the range $0.4\leq Y_{Z}\leq 0.5$. We included the $\Delta$ isobars as the most important baryon resonances for this range of temperatures and densities with their effective masses. In the meson sector we included the lightest triplet of pions as an ideal gas of particles with either the vacuum dispersion law or (for $Y_{Z}\neq 0.5$) with the $s$-wave pion-nucleon interaction taken into account.

We used the fact, cf.  Maslov:cut , that within an RMF model the EoS becomes stiffer for $n>n^{*}>n_{0}$, if a growth of the scalar field as a function of the density is quenched and the nucleon effective mass becomes weakly dependent on the density for $n>n^{*}$. In refs.  Maslov:cut ; Maslov:2015msa ; Maslov:2015wba ; Kolomeitsev:2016ptu it was demonstrated how this cut-mechanism can be realized in the $\sigma$, $\omega$, or $\rho$ sectors. In KVORcut-based models the cut-mechanism is realized in $\omega$ sector, in MKVOR-based models the cut-mechanism is realized in $\rho$ sector. Thereby the density dependence of the symmetry energy proves to be essentially different in KVORcut-based and MKVOR-based models, see fig. 1. These models allowed to pass multiple constraints from properties of cold nuclear matter, neutron star observations, even with an inclusion of hyperons and $\Delta$ isobars, and heavy-ion collision flow analysis. In the given paper we applied the KVORcut03-based and MKVOR*-based models Kolomeitsev:2016ptu for the description of the matter formed in heavy-ion collisions. The generalization to the finite temperature case was done in the standard way by introducing the temperature dependence in the fermion distribution functions and adding the pion thermal excitation term. The MKVOR* extension of the MKVOR model Kolomeitsev:2016ptu prevents vanishing of the effective nucleon mass at high density.

We found a redistribution of the charge initially concentrated in the proton subsystem of colliding nuclei between components in isospin-asymmetric systems. It was found that with an increase of the temperature the light pion subsystem becomes more isospin asymmetric, while the baryon subsystem becomes more isospin symmetric, see fig. 2.

Thermodynamical characteristics on the $T-n$ plane were considered for $0.4\leq Y_{Z}\leq 0.5$ and the energy–density and entropy–density isotherms were constructed, shown in the figs. 3, 4. We further applied our results to the description of heavy-ion collisions for the collision energies $E\stackrel{{\scriptstyle\scriptstyle<}}{{\phantom{}{}_{\sim}}}2A$ GeV. To understand the physical picture of the phenomena clearer, we assumed validity of the expanding fireball model although our main results remain valid locally and can be used in hydrodynamical calculations. As in Voskresensky:1989sn ; Migdal:1990vm ; Voskresensky:1993ud , we assumed that at energies less than a few $A$ GeV the energy in the center-of-mass system of the nucleus-nucleus collision, ${\cal E}_{\rm c.m.}A_{\rm part}$, which corresponds to the nucleons-participants, is spent on the creation of a quasi-equilibrium nuclear fireball at the end of the compression stage (labeled $t=0$), characterized by the temperature $T(0)=T_{m}$ and the baryon density $n(0)=n_{m}$ corresponding to the minimum of the energy per baryon, $E(n,T_{m})/n$, as a function of the baryon density for $T=T_{m}$. Then exploiting an assumption of approximately isoentropic expansion of the system we found the temperature as a function of the baryon density in an expanding fireball and performed a best fit to describe the $\pi^{-}$ production differential spectra, see fig. 6, the ratios of $\pi^{-}$ to proton multiplicities $R_{\pi^{-}Z}$, see fig. 8, and $\pi^{-}$ to $\pi^{+}$ ratios, see fig. 9. The effects of taking into account the $\Delta$ isobars and the $s$-wave pion-nucleon interaction (for $Y_{Z}\neq 0.5$) on pion differential cross sections, pion to proton and $\pi^{-}/\pi^{+}$ ratios were studied. At the assumption of a prompt breakup the contribution to the $\pi^{-}$ yields from the $\Delta$ decays was evaluated using the in-medium effective masses of $\Delta$s and nucleons. Compared to the works mentioned above, we extended our consideration to the case of the charge ratio $Y_{Z}\neq 0.5$, for instance taking $Y_{Z}\simeq 0.4$ in case of Au+Au and La+La collisions. This isospin asymmetry, despite being small compared to that allowed in neutron stars or supernovae, plays a role in description of the heavy-ion collisions. For instance, taking into account the difference in the chemical potentials of neutrons and protons for $Y_{Z}\neq 0.5$ leads to a noticeable increase of the negative pion yields at the conditions of the fireball breakup.

The maximum and breakup temperatures of the fireball proved to be almost model independent quantities, see figs. 5, 7, whereas the maximum reachable density $n_{m}({\cal E}_{\rm lab})$ and the breakup density $n_{\rm b.up}({\cal E}_{\rm lab})$ proved to be higher in the MKVOR*-based model than in the KVORcut03-based model. The resulting values of the breakup temperatures and densities deduced from the pion spectra for the collision energy ${\cal E}_{\rm lab}\stackrel{{\scriptstyle\scriptstyle>}}{{\phantom{}{}_{\sim}}}800A\,{\rm MeV}$ proved to be lower than the ones required to describe the $R_{\pi^{-}Z}$. The reason of this is twofold. First, as was shown in Voskresensky:1993mw , a large yield of low-momentum pions originates from the reactions occurring before the fireball breakup, because of larger mean-free paths for soft pions. Taking into account the direct reactions could essentially affect the $R_{\pi^{-}Z}$ rates making values $n_{\rm b.up}$ and $T_{\rm b.up}$ deduced from $R_{\pi^{-}Z}$ lower. Second, we ignored effects of the $p$-wave pion-baryon interaction. Inclusion of $p$-wave effects leads to a substantial change of the temperature dependence of the pion distribution function Voskresensky:1993ud , which should be especially noticeable at low $T$ corresponding to low ${\cal E}_{\rm lab}$. Our evaluations show 30$\%$ larger $\pi^{-}/\pi^{+}$ ratios compared with the data for Au$+$Au collisions, that can be again attributed to ignoring the direct pion emission and $p$-wave effects in our calculations.

We also investigated various characteristics of the liquid-gas first-order phase transition in isospin symmetric and asymmetric systems within the same KVORcut03-based and MKVOR*-based models. Such a transition occurs for $T\stackrel{{\scriptstyle\scriptstyle<}}{{\phantom{}{}_{\sim}}}20$ MeV and for $n<n_{0}$. At these conditions $\Delta$ isobars and pions do not contribute. Thermodynamical characteristics were presented on $P-n$, $T-n$, $n_{n}-n_{p}$, $\mu_{n}-\mu_{p}$ planes, where $P$ is the pressure, $n_{n}$, $n_{p}$ and $\mu_{n}$ and $\mu_{p}$ are neutron and proton densities and chemical potentials.

The following values characterize the critical point in isospin-symmetric matter for KVORcut03 (MKVOR*) models: $n_{\rm cr}=0.34(0.32)\,n_{0}$, $T_{\rm cr}=17.4\,(16.05)\,{\rm MeV}$, and $P_{\rm cr}=0.30\,(0.25)\,$MeV$/{\rm fm}^{3}$. These values lie within the $T_{\rm cr}=(16.4\pm 2.3)\,{\rm MeV}$ band given by microscopic calculations with use of the chiral nucleon-nucleon potentials Carbone:2018kji . Besides, the KVORcut03 model passes the constraints on the critical density and pressure (and marginally on critical temperature) extracted in experimental analysis of reactions going through the compound nuclear decay and multi-fragmentation Elliott:2013pna : $n_{\rm cr}=(0.06\pm 0.02)\,{\rm fm^{-3}}$, and $P_{\rm cr}=(0.3\pm 0.1)\,{\rm MeV}/{\rm fm}^{3}$, $T_{\rm cr}=(17.9\pm 0.4)\,{\rm MeV}$. The MKVOR* model passes the constraints for the critical density and pressure, whereas the predicted critical temperature is lower than the bound extracted in Elliott:2013pna . We demonstrated the system trajectories for various heavy-ion collision energies on $P-n$ and $T-n$ planes demonstrating at which conditions the system passes instability regions, see figs. 10, 11.

We studied specifics of the liquid-gas phase transition in isospin-asymmetric matter within our models of the EoS. First, disregarding possible effects of the surface tension, following Ducoin:2005aa we solved the Gibbs conditions and constructed the equilibrium pressure in the mixed phase, see fig. 13. Our results are demonstrated in the $n_{n}-n_{p}$ and $\mu_{n}-\mu_{p}$ planes, see fig. 14, 15.

We evaluated the scaled variances of the baryon and electric charges within the phase transition region and demonstrated that they diverge at the onset of spinodal instabilities, see figs.12, 16.

In addition, taking into account the non-zero surface tension we formulated the novel possibility of the formation of a structured pasta phase in the isospin-asymmetric finite nuclear systems in the region of the liquid-gas first-order phase transition.

Concluding, in the given work we demonstrated efficiency of the KVORcut03-based and MKVOR*-based models of EoS, which have passed the check for the description of cold isospin asymmetric nuclear matter Maslov:2015msa ; Maslov:2015wba ; Kolomeitsev:2016ptu ; Kolomeitsev:2017gli , now for the description of heavy-ion collisions at the collision energies below few $A$ GeV. In the subsequent work we expect to probe our EoSs in actual hydrodynamical calculations. Besides, we will take into account effects of the $p$-wave pion-baryon interactions. Moreover, we will generalize our consideration to higher temperatures, densities and collision energies.

Acknowledgements

We would like to thank D. Blaschke, Yu. B. Ivanov, E. E. Kolomeitsev and A. S. Khvorostukhin for valuable discussions and M. Borisov and P. Lukyanov for the interest to this work. This work (Sect. 1-3) has been supported by the Russian Science Foundation under grant No. 17-12-01427. Work of K.A.M. on Sect. 4 was supported by the Foundation for the Advancement of Theoretical Physics and Mathematics “BASIS”. The work of D.N.V. on Sect. 4 was supported by the Ministry of Science and High Education of the Russian Federation within the state assignment, project No 3.6062.2017/6.7.

Appendix

In the paper body we have discussed only effects of the $s$-wave pion-nucleon interaction which are weak for $0.4\leq Y_{Z}\leq 0.5$, which we study in this work. The important role of the $p$-wave pion-baryon and kaon-baryon interactions has been intensively studied in the literature, cf. Baym:1975tm ; Migdal:1978az ; Migdal:1990vm ; Voskresensky:2016oee ; Voskresensky:2018ozf ; Kolomeitsev:1995xz ; Kolomeitsev:2002pg and refs. therein.

With taking into account the nucleon-nucleon hole, the $\Delta$ isobar-nucleon hole, and a residual interaction, the Dyson equation for the retarded pion Green function describing propagation of particles with pion quantum numbers is as follows

[TABLE]

Hatched vertices take into account $NN$ and $\Delta N$ correlations.

On the pion energy-momentum, $\omega-k$, plane there exist regions, where pions with a good accuracy can be treated as quasiparticles, and regions, where the pion spectral function differs significantly from the $\delta$ function Voskresensky:1989sn ; Migdal:1990vm ; Voskresensky:1993ud . There are three quasiparticle branches in the pion spectrum in ISM: the pion branch $\omega=\omega_{\pi}(k),\omega_{\pi}\to m_{\pi}$ for $k\to 0$ and for $n\ll n_{0}$, the $\Delta$ branch $\omega=\omega_{\Delta}(k)$, $\omega_{\Delta}(k)\to m_{\Delta}-m_{N}$ for $k\to 0$ and for $n\ll n_{0}$, and the spin-sound branch $\omega=\omega_{s}(k)$, $\omega_{s}(k)\simeq vk$ for $k\ll p_{{\rm F},N}$. A region, where the spectral function has a large width ($\omega<kv_{{\rm F},N}$ and $k\sim p_{{\rm F},N}$), describes virtual pions. Thereby the spectrum of particles with pion quantum numbers differs significantly from the dispersion law of the free pions, $\omega_{k}=\sqrt{m^{2}_{\pi}+k^{2}}$. Thus to consider pions in the baryon matter as obeying the vacuum dispersion law is an oversimplification. With taking into account the polarization in the ISM for $n=n_{0}$, $T=0$, the pion spectrum is shown in fig. 17,* cf.* Migdal:1990vm ; Voskresensky:2016oee ; Voskresensky:2018ozf ; Kolomeitsev:2000ie .

The medium effects on the pion emission from the fireball depend on a relation between the breakup time-scale $\tau_{\rm b.up}$ and the characteristic time of the quasiparticle absorption. The pion quasiparticle has enough time to transit to the vacuum spectrum branch during the system breakup, if

[TABLE]

where $\omega_{i}(k;n_{\rm{b.up}},T_{\rm b.up})$ is the pion energy on the $i$-th quasiparticle branch. In the opposite limit

[TABLE]

the breakup is prompt for pions from those branches and pion polarization should be included as before as at the fireball breakup stage, cf. Senatorov:1989cg . We point out that before the fireball breakup (for $n>n_{\rm{b.up}}$) the pion polarization should be taken into account in any case. Despite that, in the given work we focus on other features of our model and include only $s$-wave pion-nucleon effects. Generalization will be considered in the future work.