Quantum liquids resulting from quark systems with four-quark interaction

Quark ensembles inﬂuenced by strong stochastic vacuum gluon ﬁelds are investigated within the four-fermion interaction approximation. The comparative analysis of several quantum liquid models is performed and this analysis leads to the conclusion that the presence of a gas–liquid phase transition is their characteristic feature. The problem of the instability of small quark number droplets is discussed and it is argued that it is rooted in the chiral soliton formation. The existence of a mixed phase of the vacuum and baryon matter is proposed as a possible explanation of the latter stability.


Introduction
The present paper is motivated by world-wide ambitious experimental programs to study the collisions of heavy ions at ultra-relativistic energies to create the quark-gluon plasma (QGP) and to investigate its properties. An obvious great challenge for the theoretical (and phenomenological) endeavors to perform this research is to be able to follow the evolution of a system created in such collisions from its initial state to a near-equilibrium plasma state and eventually to the final hadrons which are measured by detectors. Solving the problem of equilibrating and thermalizing in the framework of microscopic quark and gluon dynamics one needs to deal with the time-dependent strongly interacting systems not only at the asymptotically weak or strong coupling. Nowadays this task is obviously unrealizable and the current response on the requirements of progressing experiment is mainly based on the phenomenological considera-This article is dedicated to the memory of my young collaborator and friend Sergey Molodtsov. a Deceased b e-mail: Gennady.Zinovjev@cern.ch tion of relativistic heavy ion collisions as evolving through the certain stages.
As is believed, the initial heavy ion collisions at relativistic energies release the hadron constituents (partons) resulting in the highly non-equilibrium strongly interacting matter. The success of applying near-equilibrium hydrodynamics to describe the subsequent evolution of such a matter to describe the existing experimental data [1][2][3][4][5] suggests the system created is rather a quark-gluon liquid. Theoretically the transition between the nonequilibrium and equilibrium stages is still far from clear. But the data analysis foresees a very fast thermalization, i.e., sufficient degree of local equilibrium or, rather, isotropization, since the equations of hydrodynamics do not include the temperature of the produced matter with explicit collective properties. In recent years, there have appeared several scenarios of what could be the dynamics of a system transiting from the initial collision state to that when it becomes (almost) equilibrium [6]. However, this problem is not discussed in the present paper. Instead, we focus on another aspect of the problem, namely, the presence of the strong interaction in a system produced or, in other words, the small mean free path of its constituents, and on trying to understand the very nature of such interactions in a system whose dynamics is governed by the coupling constant, which is likely not too large (at the LHC energies the running coupling constant in QCD is α s ∼ 0.3-0.4) avoiding AdS/CFT duality (holographic QCD) arguments very popular at the moment [7,8].
The four-fermion (QCD-like) field theories still remain a most reliable source of quantitative information in the studies of the transport properties of strongly correlated systems and their thermodynamics, in particular, a chiral phase transition between massive hadrons and almost massless quarks. It is a thermodynamics that provides us with some general framework which lets one understand how the properties of macroscopic matter and, in particular, its collective behavior, emerge from the laws that govern microscopic dynamics. The results of this work allow us to suppose with a level of argumentation that is sufficiently high, in our view, that the picture based on the complex collective behavior of quarks (antiquarks, gluons), which is expressed in the presence of vacuum condensates even under normal conditions, can be set by the nontrivial thermodynamic properties of vacuum, which eventually determine the observable properties of strongly interacting matter. In our opinion, this possibility was not sufficiently widely discussed and, even more so, used already at the initial stage of studies of the quark-gluon matter due to purely accidental circumstances.

Thermodynamics of the ensemble
In the present work we consider some aspects of the thermodynamical description of the quark ensemble with a fourfermion interaction Hamiltonian density (generated, as it is believed, by a strong stochastic gluon field) where j a μ =qt a γ μ q is the quark current, with corresponding quark operators q,q, taken in the spatial point x (the variables with prime corresponds to the y point), m is the current quark mass, t a = λ a /2 is the color gauge group SU (N c ) generators, μ, ν = 0, 1, 2, 3. We take the gluon field correlator A a μ A b ν in the simple form of a color singlet, with contact (in time) interaction (without retarding) 1 (we do not include the corresponding delta-function on time in this formula). This simple correlation function is a fragment of the corresponding ordered exponent and besides the four-fermion interaction accompanied by an infinite number of multi-fermion vertices arises. But for our purposes here it would be quite enough restrict ourselves with this simple form. The effective interactions mentioned above appear in a natural way by the coarse-grained description of the system with handling the corresponding averaging procedure, and having in mind that the vacuum gluon field changed stochastically (for example, in the form of instanton liquid; see [9][10][11]). But this elaboration of the effective Hamiltonian resulting from the first principles of quantum chromodynamics (QCD) will be unessential for us, as will be demonstrated below. The choice of correlation function in the simplest form with instantaneous interaction does not generate any problem at transforming the final results from the Minkovski space to the Euclidean one and the form factor F(x) is interpreted in a simple way as an interaction 'potential' of point-like particles. The correlation function itself looks, formally, like a gauge non-invariant object. 2 Nevertheless, there exists an effective way to significantly compensate for this shortcoming, if all similar 'potentials' are looked through, in some sense, (to be elucidated below). For example, this set would be quite perceptible, if it becomes possible to confront two limits opposite in physics, for example, starting from the form factor with a delta-like function in the coordinate space (the Nambu-Jona-Lasinio (NJL) model [12], the correlation length is finite in this case) and extended to the form factor of a delta-like function in the momentum space (clearly, the correlation length tends to infinity in this case), analogous to what is well known in condensed matter physics as the Keldysh model (KKB) [13][14][15][16].
It is worth to remark here that we will need only one of its properties, although one of exceptional importance, which is related to the fact that, due to the special form factor behavior, all the momentum integrations in the problem get factorized and effectively the problem becomes one-dimensional (then only integrations over energy are in play).
In the KKB model the fermion behavior is considered in the stochastic random field with an infinite correlation length (the NJL model corresponds to the 'white noise' with zero correlation length). In this case one is lucky enough to be able to 'sum up' an entire diverging series and, therefore, to demonstrate that fermions are in general not on mass shell. Actually, it looks like one of the possible scenarios of quark confinement although there exist the standard statements that the four-fermion interaction models do not lead to the confinement of quarks and are non-renormalizable. However, recently it was shown that the KKB model (which is akin to the NJL one) possesses the remarkable property of exact integrability (in the sense by Thirring or Luttinger) [17] and in spite of the presence of divergences (as in the well-known practical field theories) in the intermediate calculations the final results turn out to be finite. Another peculiar fact of such a model is that the Dirac sea displays a finite relative depth with momentum cutoff increasing. Another interesting quark confinement mechanism is discussed in the model in which the quarks are 'living' off mass shell, similarly to the KKB model, given in Ref. [18][19][20].
From this point of view, other models with an arbitrary form factor (including the NJL model) could be represented as a superposition of elementary blocks obtained by using the KKB model. The utmost distributions mentioned above can be considered as a limiting case for the corresponding Gaussian correlators in the coordinate and momentum spaces, 2 It is obvious that we are dealing with an approximate calculation of the corresponding generating functional for some specific conditions with the restricted area of applicability that does not imply the calculation of functional derivatives of arbitrary order. which, of course, look more realistic. The coupling constant scale G, which will turn out to be interesting for applications, can be tuned by using the corresponding PDG meson observables. Comparing the results obtained (by continuity arguments) one can draw some conclusions about the behavior of the system with practically any interaction potential.
We consider it necessary to comment briefly on the case with a linear potential, which was always giving hope to discover an unusual feature in quark behavior, thereby shedding some light on the nature of confinement. Meanwhile, at present, however, it appears that such a singular 'potential' is even superfluous for our purposes, since the properties, we are interested in, are already revealed in the KKB model, which, in a sense, is like 'half way' from the NJL model to that with a linear potential. Secondly, the quasiparticles in the model with a linearly increasing potential can not basically be distinguished from those in, for example, the NJL model, provided an integrable infrared singularity in the former is eliminated. As a result the same massive objects appear without the anomalies in the energy spectrum. Additionally, the analysis shows that the multi-fermion contributions present in the problem in the general case can be reduced to the fourfermion interaction in an acceptable way by inserting the respective vacuum expectation values. In other words, even the Hamiltonian of the form (1) seems to capture the essential features of quark interactions. Besides, already in the KKB model the imaginary parts of the appropriate effective scattering matrix in scalar, pseudoscalar, vector, and axial-vector channels are equal to zero, i.e. the free quarks do not appear in this scenario of dominating correlations. We give more comments on this point below in the chapter devoted to the polarization operator.
It is believed that at sufficiently large interaction the ground state of the system transforms from a trivial vacuum |0 (the vacuum of free Hamiltonian) to the mixed state (with quark-antiquark pairs with the opposite momenta and vacuum quantum numbers) which is presented as the Bogolyubov trial function (in that way some separate reference frame is introduced and the chiral phase becomes fixed) Here a + , a, and b + , b are the quark creation and annihilation operators, a|0 = 0, b|0 = 0. The dressing transformation T transmutes the quark operators to the creation and annihilation operators of quasiparticles A = T a T † , The thermodynamic properties of a quark ensemble are known to be determined by solving the following problem. It is required to find a statistical operator that at fixed mean charge d p = dp/(2π) 3 , (Q 0 =qγ 0 q), and fixed mean entropy S = − ln ξ , provides a minimal value of mean energy of the quark ensemble In other words, we are interested in the minimum of the following functional: where μ and T denote the Lagrangian multipliers for the chemical potential of the quark/baryon charge (which is usually taken to be three times larger than the baryon one in phenomenological considerations) and the temperature (β = T −1 ), respectively. V is the volume the system is enclosed in, γ = 2N c (in the case of several quark flavors γ = 2N c N f , where N f is the flavor number), n = Tr{ξ A + A}, n = Tr{ξ B + B} are the components of the corresponding density matrix. We restrict ourselves by considering the Bogolyubov-Hartree-Fock approximation in which the statistical operator is constructed on the basis of approximating the effective Hamiltonian H app , quadratic in creation and annihilation operators for quasiparticles acting in the corresponding Fock space with a vacuum state |σ . The average specific energy per quark w = E/(V γ ) results in [21][22][23] where , the primed variables, correspond to the integration over momentum q. The auxiliary angle θ m is determined from the relation sin θ m = m/ p 0 . The first term in Eq. (7) is introduced in view of normalizing in such a way as to have the zero energy of ground state when the interaction is switched off. This constant is unessential for the following consideration and may be omitted, however, it should be kept in mind that it will appear as a regularizer in singular expressions further in the text. The most stable extremals of the functional (7) are presented for comparison with the solid line for the NJL model and the dashed one for the KKB model under normal conditions (T = 0, μ = 0) in Fig. 1. For the delta-like potential in coordinate space (the NJL model) the expression (7) diverges and to obtain the reasonable results the upper limit cutoff in the momentum integration is introduced being one of the tuning model parameters along with the coupling constant G and current quark mass m. Below, we use one of the standard sets of the parameters for the NJL model [24]: = 631 MeV, G 2 /(2π 2 ) ≈ 1.3, m = 5.5 MeV, whereas the KKB model parameters are chosen in such a way that for the same quark current masses the dynamical quark ones in both NJL and KKB models coincide at vanishing quark momentum. The momentum p ϑ (parameter) corresponds to the maximal attraction between quark and antiquark. The value of this parameter inversed determines a characteristic size of the quasiparticle.
It is of order of p ϑ ∼ (m M q ) 1/2 , where M q is a characteristic quark dynamical mass for the models considered, i.e. the quasiparticle size is comparable with the size the of π meson (Goldstone particle). It is a remarkable fact that the quasiparticle, as is seen from Fig. 1, does not depend noticeably on the form factor profile or, in other words, on the scale, but rather depends on the coupling constant. Using the properties of extremals the functional expression (7) can be transformed to the form (see [21][22][23]) where P 0 = [ p 2 + M 2 q ( p)] 1/2 is the energy of the quark quasiparticle with the dynamical quark mass Below we omit often the arguments of the corresponding functions for the mass and quasiparticle energy. Varying the functional (8) with respect to the density of the induced quasiparticle mass M (in such a form it is convenient to take variational derivatives 3 ) we obtain the equation for the dynamical quark mass as (10) which corresponds exactly to the mean field approximation. In particular, under normal conditions (T = 0, μ = 0) the dynamical quark mass in the NJL model is M q ∼ 340 MeV, whereas the dynamical quark mass of the KKB model is determined by the following equation: In practice, it is convenient to use an inverse function p(M q ). Then in the chiral limit M q = (4G 2 − p 2 ) 1/2 , at | p| < 2G, and M q = 0 when | p| > 2G. In this case the quark states with momenta | p| < 2G are degenerate in energy P 0 = 2G. Figure 2 demonstrates three branches of solutions of Eq. (11) for the dynamical quark mass. The dots show the imaginary part of the solutions which are generated at the point where two real solution branches merge.

Mean energy as a functional of quantum liquid theory
The goal that we pursued while passing from the expression for specific energy (7) to Eq. (8) was to derive such a form as would easily be recognized as the Landau energy functional of the Fermi-liquid theory [25][26][27][28][29]. Some aspects of this theory are interesting and useful apply for comparing the results obtained in the NJL and KKB models. We will also discuss the first order phase transition, which is apparently typical for interacting fermions (relativistic Fermi liquid). Thus, the second term in (8) describes the contributions coming from quark and antiquark quasiparticles with occupation numbers n andn respectively. The unity in the expression 1 − n −n corresponds to the vacuum fluctuations of quarks and antiquarks. The last term in (8) is due to the interaction of the quasiparticles. The presence of contributions coming from antiparticles and the relativistic form of the dynamics are the features which distinguish quark ensembles we study from the Fermi liquids considered in condensed matter physics. The first variation of the functional (8) with respect to the particle (antiparticle) density leads (as it should) to the energy of the quasiparticle: We consider, first, the situation of zero temperature and discuss some aspects of filling up the Fermi sphere by quarks. Let us assume that the momentum distribution of quarks (antiquarks) is determined by the following expressions taken at the β → 0 limit: that is by the Fermi step function: n = 1, at P 0 ≤ μ and n = 0 when P 0 > μ. It is clear that for antiquarksn = 0. The quark density is determined by using the Fermi momentum: with the quark chemical potential that coincides with the quasiparticle energy on the Fermi surface, as follows from Eq. (12), i.e.
The group velocity of quasiparticles on the Fermi surface v f = ∂ P 0 /∂p| |p|=P F is shown in Fig. 3 as a function of baryon (quark) density (by definition, the baryon density is three times smaller than the quark one ρ B = ρ/3). A solid line describes the NJL model, while a dashed one corresponds to the KKB model. There are points for comparison that show a version for the KKB model when the parameters are tuned in such a way that the π -meson masses coincide in the NJL and KKB models (a similar notation is used below in Figs. 4,5,6,7,8). The group velocity tending to unity in the region of normal nuclear densities corresponds to the chiral symmetry restoration when the induced quark mass tends to zero. The group velocity turns to zero for quarks with momenta |p| < 2G in the chiral limit in the KKB model. The negative group velocities in the NJL model correspond to the regions of instability (see below). The points in which the where P 0 F = P 0 | |p|=P F , N F = dρ/dμ. For more details on how to determine the parameter F 0 , see below. The interaction term in the functional (8) vanishes in an ideal gas and causes the derivative of the quark dynamical mass in the Fermi momentum to turn to zero: dM q /d P F = 0. Let us define the density of states of an ideal gas as then the relation (16) can be written in the form Another important characteristic is the compression coefficient Figure 4 demonstrates the data for the NJL and KKB models. They are consistent with the specific values obtained for a nuclear medium. One can also conclude that, in principle, these models admit a wide variety of equations of state including sufficiently restrictive ones. The negative values of the compression coefficient are not allowed and signal the region of instability. The first sound velocity, which is determined by the relation is shown in Fig. 5. When the baryon densities are somewhat higher than the density of normal nuclear matter, the sound velocity tends to its asymptotic value C 1 = 1/ √ 3 which is a natural manifestation of the chiral symmetry restoration. If the sound velocity of an ideal Fermi gas C 2 1 = v 2 F /3 is introduced in a way similar to the N F definition, then the expressions (16), (18) can be endowed with the form whose physical meaning is an equality of flow coming through the Fermi sphere of quasiparticles of an (imaginary) ideal Fermi gas and interacting Fermi liquid (that is, there basically is a relativistic analog of the Luttinger theorem [30,31]), The thermal conductivity at a constant volume and a low temperature is given by the expression Figure 6 shows the slope (the factor 1 3 π 2 N F in Eq. (20), N F = dρ/dμ), as a function of baryon/quark density, which demonstrates how informative it could be to measure the slope of a curve corresponding to the thermal conductivity. Yet another important characteristic of the Fermi liquid is defined by the second variational derivative, which for the functional (8) develops only a scalar component For the Fermi liquid at zero temperature, in particular, we have For example, in the NJL model In particular, in the chiral limit (when m = 0) we have The collective oscillation modes of the Fermi liquid, corresponding to the so-called zero sound (the collisionless mode), are found by using the parameter which is shown in Fig. 7. In particular, in the KKB model The zero sound oscillations are known to be determined by the solutions to the dispersion equation with a frequency parameter s (for details concerning this notation, see the section devoted to the polarization operator) of the form When there is a repulsion in a system and the factor is positive F 0 > 0, the solutions to the dispersion equation s = λ + iη describe continuous oscillations (η = 0). In the case of weak attraction, when −1 < F 0 < 0, the damped oscillations of zero sound are possible with a purely imaginary frequency (λ = 0), which is given by the solutions to the following equation: When the strong attraction is available and F 0 < −1, the solutions reside on a second sheet of the complex plane s and describe the damped oscillations which are found from the solution to the equation It should, however, be recalled that these states of a Fermi liquid are unstable (as will be discussed below). It is hardly possible to apply directly the consideration of zero sound given above to the situation of interacting quarks and antiquarks under study, because here the contribution of vacuum fluctuations of antiquarks, which form along with quarks a chiral condensate, was completely ignored. On the other hand, zero sound oscillations are known to be interpreted as a bound state of a particle and hole in the vicinity of the Fermi sphere. Therefore, the excitations in a Fermi liquid should be described (in our case) by taking into account the interference between the bound states of a quark and antiquark, as well as of a quark and a hole of the Fermi sphere (the quantum numbers of that hole allows one to consider it as an anti-particle). We are doing that on calculating the polarization operator.
Turning now to the chemical potential of quasiparticles presented in Fig. 8 let us emphasize that it is seen from the data for the NJL model that there is a region of occupied states almost degenerate with respect to the chemical potential with the vacuum chemical potential of a quasiparticle that quite naturally corresponds to the vanishing Fermi momentum.
In the pioneering papers of the NJL model there was a discontinuity in the chemical potential values [32] unlike the smooth curve in this plot.
Similarly, the chemical potential of occupied states in the KKB model differs from that in vacuum by a small quantity proportional to the quark current mass All the states with momentum |p| < 2G are degenerate with respect to the chemical potential in the chiral limit. We have . Such a behavior of the chemical potential is a consequence of a rapid decrease of the dynamical quark mass with increasing Fermi momentum [see also (15)]. It follows from Eq. (8) that the Fermi sphere is being filled as though from within. Those quarks with momenta smaller than the Fermi one |p| < P F do not take part in forming a condensate. As a result, the quark dynamical mass can only decrease with the Fermi momentum increasing. This dynamical mass is independent of the quark momentum in the NJL model because of the approximation assumed. This dependence should be taken into account in a more realistic case as an analysis of the KKB model shows.
It turns out that the pressure of some occupied states degenerate in the chemical potential almost coincides with that of vacuum (the pressure of a dilute Fermi gas; T = 0) where E = E/V is the specific energy. Below we analyze respective data in more detail including the situation with nonzero temperature. The energy (and, hence, the pressure) of the ensemble is a discontinuous functional of the quark current mass (see [9][10][11]) in the KKB model. The integrands in (8) are then estimated as follows: and we find a linearly diverging integral for the specific energy of ensemble despite the fact that the delta-like form factor in the momentum space is the strongest regularizer. It is paradoxical that any small value of the current mass m leads to the negative infinite energy of ensemble, while the expression w| m=0 is well defined in the chiral limit. Even more so, a similar divergence occurs in the case of a delta-like form factor in the coordinate space. This fact is concealed by introducing the cutoff momentum in the NJL model. Now it looks quite sensible to consider the relative pressure of quark ensemble in comparison with a (formally infinite) vacuum value because of the singular character (mentioned above) of the ensemble pressure in the KKB model. The pressure derivative in the ensemble density has the form: d P/dρ = ρ dμ/dρ. Therefore, one can conclude by using an estimate given in (23) that the occupied states with momenta |p| < 2G are observed to be degenerate with respect to the pressure (E = 2Gρ, μ = 2G) in the chiral limit in the KKB model. The deviations are proportional to the quark current mass beyond the chiral limit. Now, we are able to analyze some thermodynamic properties of a system and to consider, first, the pressure of the quark ensemble in detail By definition, the volume derivative should be calculated at the constant mean entropy, dS/dV = 0. Implementing this constraint, one can, for example, extract the volume derivative of the chemical potential, dμ/dV . However, this approach cannot be implied because the mean charge conservation might be broken. In fact, there is only one possibility to satisfy both conditions by introducing two independent chemical potentials for quarks and antiquarks separately. We use a symbol μ introduced earlier for the quark chemical potential, whereas the antiquark chemical potential is taken with an opposite charge and is denoted byμ. Then, we have n = 1 e β (P 0 −μ) + 1 ,n = 1 e β (P 0 +μ) + 1 for the quark and antiquark densities, respectively. Some nonequilibrium states of quark ensemble could also be described on this way (formally with a loss of covariance, just similar to the electrodynamics as for the situation of electron-positron gas). However, we are here interested only in the special configuration whenμ = μ. The partial derivative of the specific energy dw/dV can be presented in the following form: Dealing with the definition of an induced quark mass (9) and presenting the trigonometric factors via the quark dynamical mass we find the ensemble pressure as The condition of mean charge conservation gives the first equation that interrelates the derivatives dμ/dV and dμ/dV . Here, a regularized expression for the mean charge of quarks and antiquarks is assumed modulo respective vacuum contribution. Implementing the condition of constant mean entropy dS/dV = 0 in a similar way one can obtain the second equation of the chemical potential derivatives system as follows: Substituting the expressions T ln n 1 − n = μ − P 0 and T lnn 1 −n = −μ − P 0 into this equation and collecting similar terms we attain the following equation if the conditionμ = μ and Eq. (25) are satisfied. Finally, we have for the pressure Then the thermodynamic potential should obey the following thermodynamic identity: as it should. At low temperatures the antiquark contribution is small and the thermodynamic description can be approximately developed by using the chemical potential μ only. If the antiquark contribution becomes significant, the thermodynamic description is more sophisticated and should obviously include the chemical potentialμ with the additional condition μ = μ. Figure 9 shows the ensemble pressure P in MeV/fm 3 as a function of the charge density Q 0 /3V for various temperatures. The lowest curve is obtained at zero temperature. The next curves and the ones following upwards correspond to temperatures T = 10 MeV, T = 50 MeV (an upper curve) with a step T = 10 MeV. Let us also remember that the pressure of the vacuum for the NJL model was estimated in [9][10][11] to be 40-50 MeV/fm 3 , which is quite consistent with that obtained in the bag model. It was also demonstrated that there is a region of instability within a certain interval of the Fermi momenta generated by the anomalous behavior  of pressure d P/dn < 0 (see also [33][34][35][36]). Figure 10 displays fragments of isotherms shown in Fig. 9 (but now in different coordinates) in the form of a chemical potential as a function of the ensemble pressure. A top curve is obtained at zero temperature. The isotherms below are shown in steps of 10 MeV. The lowest curve is obtained at a temperature 50 MeV. It is clearly seen from the figure that there are states on the isotherms which are in thermodynamic equilibrium. The pressure and chemical potential are the same for these  Fig. 9 by a dashed curve. The points at which a dashed curve intersects with an isotherm give a boundary for a gasliquid phase transition. The respective line P = const cuts off nonequilibrium and unstable fragments of an isotherm and describes a mixed phase. The critical temperature turns out to be equal to T c ≈ 46 MeV with the critical charge den-sityQ 0 ≈ 0.12 charge/fm 3 for the above mentioned tuning parameters. Figure 11 shows the isobars. The pressure next to each curve is given in MeV/fm 3 . The vacuum pressure corresponds to approximately ∼50 MeV/fm 3 . It is possible to extrapolate isobars into the region of small charge densities, however, it is not really necessary. The figure clearly demonstrates the presence of dilute (a gas) and dense (a liquid) phases in the vicinity of the vacuum isobar. Figure 12 shows the quark dynamical mass M q (in MeV) as a function of the chemical potential μ (in MeV) for temperatures T = 0 MeV, T = 100 MeV in steps of T = 10 MeV. The right-most curve corresponds to zero temperature. At low temperatures, below 50 MeV, the quark dynamical mass is a multi-ciphered function of the chemical potential (but the discontinuous solution has been presented in the pioneering papers). Figure 13 shows the quark dynamical mass as a function of temperature at small charge density Q 0 ∼ 0. This picture is easily recognizable in the context of the NJL model. It is the latter that is implied in a scenario of chiral invariance restoration under extreme temperatures higher than 100 MeV and with a highly diluted quark ensemble. We have already noted (see also [9][10][11]) that the momentum p θ , which corresponds to the strongest quark-antiquark attraction d sin θ/d p = 0,  can be determined. For example, for the NJL model this parameter is equal to Its inverse value is given by the characteristic effective size of a quasiparticle r θ = p −1 θ . From the behavior of the quark dynamical mass as a function of temperature at small charge densities (see Fig. 13) one can conclude that the quasiparticle size grows with the energy increasing.
In [21][22][23] it was shown that if the quark chemical potential is defined as the energy necessary to add (remove) one quasiparticle, μ = dE/dN , then the chemical potential in vacuum coincides with the quark dynamical mass (see also (12), (15)). Therefore, it seems to be reasonable to consider a QCD phase diagram by starting from this value of the chemical potential, though formally it can be taken smaller than the quark dynamical mass. In particular, we exactly reproduce the standard picture [24], [32] by taking the chemical potential equal zero. The results obtained allows one to conjecture that the phase transition of the (partial) restoration of the chiral invariance could already be realized in nature as a mixed phase of physical vacuum and baryonic matter. An indirect confirmation of this hypothesis can be seen in degenerate excited states of some baryons (see, for instance, [37]). It is, however, clear that the data presented (in particular, on the temperature and density of a critical point position) should be understood as just estimates. The critical temperature of the gas-liquid transition for nuclear matter extracted from experiment is estimated to be about 20 MeV.
(Let us mention here that the hadronization temperature as extracted from the RHIC and LHC data analysis in the framework of statistical hadronization models looks like 165 MeV.) In addition, here (at T = 0) a gas component possesses the nonzero density of order of 0.01 of the normal nuclear density, whereas the observed value should correspond to physical vacuum, i.e., to zero baryon density. It should be noted that although such an uncertainty is inherent in the other predictions of the chiral symmetry restoration phase transition which are widely discussed in many papers, they are somewhere around 2-6 normal nuclear matter densities.

Polarization operator
Returning to the discussion of zero sound and excitations of a chiral condensate we would like also to recall that this knowledge is necessary for a more consistent analysis of the transition gas-liquid layer. To this end, we will need to know a polarization operator of the form where is a respective density of the polarization operator in the channels = 1, iγ 5 , γ μ , γ 5 γ μ with the Green function of quark with the dynamical mass M q μ = μγ 0 , where p, q are the incoming and outgoing external momenta of the quark quasiparticles. It will be enough for our purposes to consider the quasiparticles with momenta where ε = p 0 − q 0 is the transferred energy, M p = M q ( p), E p = [ p 2 + M 2 q ( p)] 1/2 , the quantity F(k) is a form factor, and for the kinematics chosen p = k + Q/2 (here n p is the occupation number for a quasiparticle with momentum p). In particular, at zero temperature we have the Fermi step: n p = n(E p − μ). A similar notation is also used for a quasiparticle with momentum q = k − Q/2.
The first term in Eq. (32) corresponds to the quark and antiquark contributions, whereas the second one comes from the quark and hole configuration. It is easy to see that at F(k) = δ(k) (in the KKB model) we have cubic dispersion relations to determine the bound states: 1−2G π,σ (ε, Q) = 0. As an example, Fig. 14 shows the energies (in units of MeV) calculated of π (dashed line) and σ (solid line) mesons as functions of the momentum Q/2 (in units of MeV) at zero temperature for the gas with small baryon density corresponding to the Fermi momentum P F ∼ 130 MeV. The region of degeneracy, seen in Fig. 14 at low quark momenta, is a consequence of the above discussed fact that the Fermi sphere is filled from the inside, and quarks with momenta smaller than P F do not participate in forming the quark dynamical mass. Such a behavior is not observed in the NJL model because of the approximation adopted (the quark mass is independent of the momentum). Figure 15 shows the energies (in units of MeV) of π (dashed line) and σ (solid line) mesons as the functions of the baryon density (T = 0). The dots indicate a branch corresponding to the quark-hole bound state which appears to be degenerate for π and σ mesons. Just these branches correspond to the third additional root of the dispersion equation mentioned above, albeit so that there are only two roots (see the discussion of the NJL model). To be specific, the quark momentum is assumed to be 50 MeV larger than the Fermi one but the hole momentum is 50 MeV smaller than the latter in this example. As follows from Eq. (32), the polarization operator in the NJL model is defined by integrating over the running quark momentum k and is represented as a superposition of branches of the KKB model, which has already been mentioned in the introduction. The most significant contributions (for the kinematics we chose) are those coming from the terms denoted as a and c in Eq. (32). Integrating over the angle (it is more convenient to express the final formula by going to a nonsymmetric integration point; the corrections become negligibly small) one can obtain (at T = 0) the following results: π,σ = A π,σ + B π,σ , (22) with the parameter s = E F ε/(k F Q). The first component A π,σ results from the contribution coming from a quark-antiquark pair, whereas in the case of the second one B π,σ arises due to the coupling of quark and hole residing in the vicinity of the Fermi sphere. It should be noted that for a quark ensemble we consider the medium properties which are mainly governed by the term A π,σ responsible for the quark-antiquark condensate, contrary to what we have in condensed matter physics, where the dominant contribution, as is well known, is given by B π,σ . Therefore, the results obtained exclusively by using an analogy with the condensed matter physics should be taken with a grain of salt. In particular, in the present paper we have analyzed in detail a situation with the zero sound description taken as an example illustrating just this point. The zero sound would represent in itself the highly damped oscillations described by the only scalar parameter F 0 , while with no an antiquark present. A more accurate analysis shows that, for example, there is a stable branch of quark and hole excitations in the Fermi sphere in addition to a paired quark-antiquark state in the KKB model. We observe a regular mass convergence for π and σ mesons when the baryon density increases by performing a numerical integration in the NJL model. This effect is clearly related to the restoration of chiral symmetry. The influence of a bound quark-hole state in the Fermi sphere turns out to be insignificant. For instance, for the densities of the order of normal nuclear matter the dispersion law changes by a few MeV when the quark and hole momentum differs more than 200 MeV, but there are no damped oscillations as in the KKB model.
One of the drawbacks of the models studied so far is the lack of quark confinement that is understood here simply as the impossibility to observe a single particle state with a regular (real) dispersion law. We see formally that one quasiparticle can freely propagate, indeed. But adding just another quasiparticle can dramatically change the picture due to the existence of a bound channel. For example, in the KKB model the bound states in scalar, pseudoscalar, vector, and axial-vector channels appear at any quasiparticle momenta (details can be found in [38][39][40]). In particular, the bound state energy, obtained by using the dispersion equation 1 − 2G = 0, has the form in the π and σ channels (an upper sign corresponds to the pseudoscalar channel). The first term in this expression is the energy of free particle motion. The second one is strictly positive at any momenta p and q and plays a role of binding energy in π and σ channels (only in the configuration of q = p the binding energy vanishes for a scalar channel). Similarly, one can show that a quark and antiquark are always coupled in vector and axial-vector channels, i.e. the scattering matrix is always singular except for a tensor channel where it is trivial because of the initial interaction Hamiltonian, which is taken as a product of two color currents. Similar bound states exist in the diquark channel. As a consequence, the states with any number of quark quasiparticles turn out to be the bound states in the channels we have just mentioned. The same behavior is observed in the NJL model where the bound states appear for the quarks with momenta somewhat lower than the cutoff momentum, i.e. the scattering matrix is also singular within this momentum interval as in the KKB model. It seems the bound states appear rather due to the fermion correlations than the physical influence of a field that is familiar in the quantum electrodynamics. Then, in order to understand what may take place beyond the cutoff momentum one apparently has to study the appropriate nonlocal models.

Transition layer between gas and liquid
The concept of a mixed phase of physical vacuum and baryonic matter would receive a substantial confirmation if we are able to demonstrate the existence of the boundary (transition) layer where a transformation of the quark ensemble from one aggregate state to another takes place. As was argued above the indicative characteristic to explore a homogeneous phase (at finite temperature) is the mean charge (density) of the ensemble. All the other characteristics, for example, a chiral condensate, dynamical quark mass, etc. can be reconstructed if one knows the ensemble mean charge. So, here we analyze a specific case of the surface (transition) layer at zero temperature. We assume that the quark ensemble parameters in a gaseous phase are approximately the same as those at zero charge ρ g = 0, i.e. as in vacuum (minor differences in pressure, chemical potential, and quark condensate are neglected). The dynamical quark mass develops here the maximal value, and it is M = 335 MeV for the parameter choice standard for the NJL model. Then, as the Van der Waals diagram shows, a liquid phase, being in equilibrium with a gas phase, gains the density ρ l = 3×0.185 charge/fm 3 (by a reason which becomes clear below we correct it to take the value ρ l = 3 × 0.157 charge/fm 3 ). The detached factor 3 here links again the magnitudes of quark and baryon matter densities. The quark mass is approximately * M≈ 70 MeV in this phase. Hereafter we focus on describing two adjoining semi-infinite layers (i.e. assuming a plane symmetry of the corresponding one-dimensional problem).
The previous experience teaches that an adequate description of heterogeneous states can be reached with the mean field approximation [41,42]. In our particular case it means making use of the corresponding effective quark-meson Lagrangian [43][44][45] (functional of the Ginzburg-Landau type) where and σ is the scalar field, V μ is the field of vector mesons, m σ , m v are the masses of scalar and vector mesons and g σ , g v are the coupling constants of the quark-meson interaction.
The U (σ ) potential includes the nonlinear σ field interaction terms up to the fourth order, for example. For the sake of simplicity we do not include the contributions coming from the pseudoscalar and axial-vector mesons. The meson component of such a Lagrangian should be self-consistently treated by considering the corresponding quark loops. Here we do not see any reason to go beyond the well elaborated and reliable one loop approximation (33) [43][44][45], although recently considerable progress has been reached in scrutinizing the non-homogeneous quark condensates by applying the powerful methods of exact integration [46][47][48][49][50][51]. Here we believe it is more practical to adjust phenomenologically the effective Lagrangian parameters based on the transparent physical picture. It is easy to see that handling (33) in one loop approximation we arrive, in actual fact, at the Walecka model [52][53][54] but adopted for the quarks. In what follows we are working with the designations of that model and we hope that it does not lead to misunderstandings.
In the context of our paper we propose to interpret Eq. (33) in the following way. Each phase might be considered, in a sense, with regard to another phase as an excited state which requires the additional (apart from a charge density) set of parameters (for example, the meson fields) for its complete description, and those are characterizing the measure of deviation from the equilibrium state. Then the crucial question becomes whether it is possible to adjust the parameters of the effective Lagrangian (33) to obtain the solutions in which the quark field interpolates between the quasiparticles in the gas (vacuum) phase and the quasiparticles of the filled-up states. For all that, the density of the filled-up state ensemble should asymptotically approach the equilibrium value of ρ l and should turn to the zero value in the gas phase (vacuum).
The scale inherent in this problem may be assigned by one of the masses referred to in the Lagrangian (33). In particular, we bear in mind the dynamical quark mass in the vacuum M. Besides, there are another four independent parameters in the problem and in order to compare them with the results of a study of the nuclear matter we employ the form characteristic for the (nuclear) Walecka model, Parameterizing the potential U (σ ) as b σ = 1.5 m 2 σ (g σ /M), c σ = 0.5 m 2 σ (g σ /M) 2 and we arrive at the sigma model whereas the choice b = 0, c = 0 results in the Walecka model. As to standard nuclear matter, application of the parameters b and c demonstrates the vital model dependent character and they are quite different from the parameter values of the sigma model. Truly, in that case their values are also regulated by additional requirement of an accurate description of the saturation property. On the other hand, for the quark Lagrangian (33) we could intuitively anticipate some resemblance with the sigma model and, hence, could introduce two dimensionless parameters η and ζ in the form of b = η b σ , c = ζ 2 c σ , which characterizes some fluctuations of the effective potential. Then the scalar field potential is presented as follows: The meson and quark fields are determined by solving the following system of stationary equations: where * M= M + g σ σ is the running value of the dynamical quark mass, E stands for the quark energy and V = −i V 4 . The density matrix describing the quark ensemble at T = 0 has the form in which p is the quasiparticle momentum and the Fermi momentum P F is defined by the corresponding chemical potential. The densities ρ s and ρ at the right hand sides of Eq. (34) are by definition Here we confine ourselves to the Thomas-Fermi approximation while describing the quark ensemble. Then the densities which we are interested in are given with some local Fermi momentum P F (x) as where γ is a quark gamma-factor which for one flavor is ) 1/2 and λ = * M /P F . Under the assumption adopted the ensemble chemical potential is constant and, therefore, the local value of the Fermi momentum is defined by the running value of dynamical quark mass and vector field as Now we should tune the Lagrangian parameters in Eq. (33). For asymptotically large distances (in a homogeneous phase) we may neglect the gradients of scalar and vector fields, and the equation for the scalar field of the system (34) leads to the first equation that relates the parameters C s , The vector field asymptotically is given by the ensemble density V = C 2 v ρ/(g v M 2 ). The second equation derived from the relation (36) for the chemical potential looks like If we know the liquid density we obtain the Fermi momentum (P F = 346 MeV) from (35). Applying the identities (37), (38) we have for the particular case b = 0, c = 0 the results C 2 s = 25.3, C 2 v = −0.471, i.e. the vector component C 2 v is small (compared to C 2 s ) and acquires a negative value that is unacceptable. Apparently, it looks necessary to abandon the contribution coming from the vector field or to reduce the dynamical quark mass * M up to the value which retains the identity (38) valid with positive C 2 v or even zero value. In the gaseous phase the dynamical quark mass can also be corrected to a value larger than the vacuum value. It is clear that in the situation of the liquid with the density ρ l = 3 × 0.185 ch/fm 3 the dynamical quark mass should coincide (or exceed) M = 346 MeV in the gaseous phase. However, here we correct the liquid density (as was argued above) to decrease its value up to ρ l = 3 × 0.157 ch/fm 3 , which is quite acceptable in the nucleation capacity. In fact, this possibility can be simply justified by another choice of the NJL model parameters. Thus, we obtain at * M= 70 MeV and b = 0, c = 0 the result that C 2 s = 28.4, C 2 v = 0.015, i.e. we have a small but positive value for the vector field coefficient. At the same time, aiming here to estimate the surface tension effects only we do not strive for the precise fit of the parameters. In the Walecka model these coefficients are C 2 s = 266.9, C 2 v = 145.7 (b = 0, c = 0). Moreover, there is another parameter set with C 2 s = 64, C 2 v ≈ 0 [55][56][57][58][59] but it is rooted in the essential nonlinearity of the sigma-field due to the nontrivial values of the coefficients b and c. The option (formally unstable) with negative c (b) has also been discussed.
The coupling constant of scalar field is fixed by the standard (for the NJL model) relation between the quark mass and the π -meson decay constant g σ = M/ f π (we put f π = 100 MeV) although there is not any objection to treat this coupling constant as an independent parameter. As a result of the whole agreement we have for the σ -meson mass m σ = g σ M/C s . In principle, we could even fix the σ -meson mass and coupling constant g σ , but all relations above mentioned lead eventually to quite suitable values of the σ -meson mass as will be demonstrated below. The vector field plays, as we see, a secondary role because of the small magnitude of the constant C v . Then taking the vector meson mass as m v ≈ 740 MeV (a slightly smaller value than the mass of the ω meson because of simple technical reasons only) we calculate the coupling constant of the vector field from the relation similar to the scalar field m v = g v M/C v . Amazingly, its value turns out to be steadily small being compared to the value characteristic for the NJL model, g v = √ 6g σ . However, at the values of constant C v which we are interested in it is very difficult to maintain the reasonable balance, and to be specific in this paper we prefer to choose the massive vector field. Actually, it is unessential because we need this parameter (as we remember) to estimate the vector field strength only.
The key point of our interest here is the surface tension coefficient [55][56][57][58][59] which can be defined as The parameter r o will be discussed in the next section on considering the features of the quark liquid droplet, and for the present we would like to notice only that for the parameters considered its magnitude for N f = 1 is around r o = 0.79 fm. Recalling the factor 3 1/3 which connects the baryon and quark numbers, we find the magnitude ( r o = 3 1/3 0.79 ≈ 1.14 fm) in full agreement with the magnitude standard for the nuclear matter calculations (in the Walecka model) In order to proceed we calculate E(x) in the Thomas-Fermi approximation as To give some idea for the 'setup' prepared we present here the characteristic parameter values for some fixed b and c with ρ l = 3 × 0.157 ch/fm 3 . In the liquid phase they are * M= 70 MeV (P F = 327 MeV) and e l = 310.5 MeV (the index l stands for a liquid phase and e(x) = E(x)/ρ(x) defines the density of the specific energy). Both relations (37) and (38)  Two other parameters, η, ζ , are fixed by looking through all the configurations in which the solution of the equation system (34) with a stable kink of the scalar field does exist and describes the transition of quarks from the gaseous phase to the liquid one. First, it is reasonable to scan the η, c (ζ = c η)plane, in order to identify the domain in which the increase of specific energy E − E l ρ/ρ l ≤ 0 is revealed at running through all possible states which provide the necessary transition (without taking into account the field gradients). In practice one needs to follow a simple heuristic rule. The state with P F ∼ 1 MeV (i.e. e and the corresponding ρ) and the state of the characteristic liquid energy E l (together with ρ l ) should be compared while scanning the Lagrangian parameters η and c. Just the domain in which they are commensurable could provide us with the solutions which we are interested in, and Fig. 5 shows its boundary. The curve could be continued beyond the value η = 2.5, but the values of corresponding parameter η are unrealistic and not shown in the plot. We calculate the solution of equation system (34) numerically by the Runge-Kutta method with the initial conditions σ (L) ≈ 0, σ (L) ≈ 0 imposed at the large distance L t, where t is a characteristic thickness of transition layer (about 2 fm). Such a simple algorithm appears to be quite suitable if the vector field contribution is considered as a small correction (this just takes place in the situation under consideration) and is presented as where the charge (density) ρ is directly defined by the scalar field. We considered the solutions including the contribution of the vector field as well and the corresponding results confirm the estimates obtained. A rather simple analysis shows that interesting solutions are located along the boundary of discussed domain. Some of those are depicted in Fig. 16 as the dots. Figure 17 shows the stable kinks of the σ field with the parameter c = 1.1 for two solutions with η ≈ 0.977 (m σ ≈ 468 MeV) (solid line) and η ≈ 1.813 (m σ ≈ 690 MeV) (dashed line). For the sake of clarity we consider the gas (vacuum) phase is on the right. Then the asymptotic value of σ field on the left hand side (σ ≈ 80 MeV) corresponds to * M= 70 MeV. The thickness of the transition layer for the solution with η ≈ 0.977 is t ≈ 2 fm, whereas for the second solution we have t ≈ 1 fm.
Characterizing the whole spectrum of the solutions obtained we should mention that there exist other more rigid (chiral) kinks, which correspond to the transition into the state with the dynamical quark mass changing its sign, i.e. M → −M. In particular, for the kink with the canonical parameter values η = 1, c = 1 is clearly seen (marked by the star in Fig. 16) and its surface tension coefficient is about 2m π (m π is the πmeson mass). The most populated class of solutions consists of those having metastable character. The system comes back to the starting point (after an evolution) pretty rapidly, and usually the σ field does not evolve to such an extent as to reach the asymptotic value (which corresponds to the dynamical quark mass in the liquid phase * M= 70 MeV). Switching on the vector field changes the solutions insignificantly (for our situation with small C v it does not exceed 2 MeV in the maximum).
The surface tension coefficient u s in MeV for the curve of stable kinks with parameter η ≤ 1.2 as a function of the parameter c (ζ = c η) is depicted in Fig. 18. The σ -meson mass at c ≈ 0 is m σ ≈ 420 MeV and changes smoothly up to the value m σ ≈ 500 MeV at c ≈ 1.16 (the maximal value of the coefficient c beyond which the stable kink solutions are not observed). In particular, m σ ≈ 450 MeV at c = 1. Two kink solutions with c = 1.1 for η ≈ 0.977 and for η ≈ 1.813 (shown in Fig. 17, and the second one is not shown in Fig. 18) have the tension coefficient values u s ≈ 35 MeV and u s ≈ 65 MeV, respectively. The maximal value of the tension coefficient for the normal nuclear matter does not exceed u s = 50 MeV. The nuclear Walecka model claims the value u s ≈ 19 MeV [55][56][57][58][59] as acceptable and calculable. The reason to have this higher value of the surface tension coefficient for the quarks is rooted in the different values of the mass deficit. Indeed, for nuclear matter it does not exceed * M≈ 0.5M albeit so that more realistic values are considered around * M≈ 0.7M and for the quark ensemble the mass deficit amounts to * M≈ 0.3M. We are also able to estimate the compression coefficient of the quark matter K , which occurs significantly larger than the nuclear one. Actually, we see a quite smooth analogy between the results obtained and the results of the bag soliton model [60][61][62][63]. The thermodynamic treatment developed in the present paper allows us to formulate the adequate boundary conditions for the bag in physical vacuum and to diminish considerably the uncertainties in searching the true soliton Lagrangian. We believe, and it was also shown here, that to single out one soliton solution among others (including even those obtained by the exact integration method [46][47][48][49][50][51]), which describes the transitional layer between two media, is not an easy problem if the boundary conditions above formulated are not properly imposed.

Droplet of quark liquid
The results of the previous sections have led us to put forward the challenging question about the creation and properties of finite quark systems or the droplets of quark liquid which are in equilibrium with the vacuum state. Thus, as a droplet we imply the spherically symmetric solution of the equation system (34) for σ (r ) and V (r ) with the obvious boundary conditions σ (0) = 0 and V (0) = 0 in the origin (the primed variables denote the first derivatives in r ) and rapidly decreasing at large distances, σ → 0, V → 0, when r → ∞.
A quantitative analysis of similar nuclear physics models which includes the detailed tuning of parameters is usually based on the comprehensive fitting of the available experimental data. This way is obviously irrelevant in studying the quark liquid droplets. This global difficulty dictates the specific tactics of analyzing. We propose to start, first of all, with selecting the parameters which could be worthwhile to play the role of physical observables. Naturally, the total baryon number which phenomenologically (via factor 3) related to the number of valence quarks in an ensemble, is a reasonable candidate for this role. Besides, the density of the quark ensemble ρ(r ), the mean size of the droplet R 0 , and the thickness of the surface layer t looks suitable for such an analysis.
It is argued above that the vector field contribution is negligible because of the small value of the coefficient C v compared to the C s magnitude, and we follow this conclusion (or assumption) albeit so that we understand it is scarcely justified in the context of a finite quark system. Thus, we put g v = 0, V = 0 in what follows and it simplifies all the calculations enormously. Figure 19 shows the number of solutions (σ field in MeV) to the system (34) at N f = 1 and Fig. 20 presents the corresponding distributions of ensemble density ρ (ch/fm 3 ). The parameters C s , C v , b, and c are derived by the same algorithm as in the previous section, i.e. the chemical potential of the quark ensemble M = 335 MeV (and σ → 0) is fixed at spatial infinity. The filled-up states (of a liquid) are characterized by the parameters * M= 70 MeV, ρ 0 = ρ l = 3×0.157 ch/fm 3 . The σ -meson mass and the coupling constant g σ are derived at fixed coefficients η and ζ , and they just define the behavior of the solutions σ (r ), ρ(r ), etc. The magnitudes of the functions σ (r ) and ρ(r ) at the origin are not strongly correlated with the values characteristic for the filled-up states and are practically determined by solving the boundary value problem for system (34). In particular, the solutions presented in Fig. 19 have been obtained with the running coefficient η at ζ = η. The most relevant parameter (instead of η) from the physical view point is the total number of quarks in the droplet N q (as discussed above) and it is depicted to the left of each curve. (The variation of * M, ρ 0 and f π could be con-  Table 1 The results of fitting the Fermi distribution with N f = 1, ρ 0 (ch/fm 3  sidered as well instead of the two mentioned parameters η and ζ .) Analyzing the full spectrum of solutions obtained by scanning, one can reveal a recurrent picture (at a certain scale) of kink droplets which are easily parameterized by the total number of quarks N q in a droplet and by the density ρ 0 . These characteristics are obviously fixed at completing the calculations. The property which allows us to single out these solutions is related to the value of droplet specific energy (see below). Table 1 exhibits the results of fitting the density ρ(r ) with the Fermi distribution where ρ 0 is the density in the origin, R 0 is the mean size of droplet, and the parameter b defines the thickness of the surface layer t = 4 ln(3)b. Besides, the coefficient r 0 , which is absorbed in the surface tension coefficient (39), the σmeson mass, R 0 = r 0 N 1/3 q , and the coefficient η at which all other values have been obtained are also presented in Table 1. The curves plotted in Fig. 19 and results of Table 1 justify us to conclude that the density distributions at N q ≥ 50 are in full agreement with the corresponding data typical for the nuclear matter. The thicknesses of the transition layers in both cases are also similar and the coefficient r 0 with the factor 3 1/3 included is in full correspondence with r 0 . The values of the σ -meson mass in Table 1 look quite reasonable  as well. Obviously, we are unable to tune our model, unlike the Walecka one, by drawing the experimental data on the quark drops. However, the subject is getting increasingly interesting because of possible astrophysical applications in studying the anomalous quark stars [55][56][57][58][59]. However, there is the possibility to trace some qualitative objective observations, in particular, the central density of drops is getting smaller with the quark quasiparticles number in the drop decreasing. Meanwhile, we know from the experiments that a nuclear matter develops an opposite tendency of a certain increase of nuclear density that becomes quite noticeable (factor 2.5) for helium and is much larger (factor 10) than the standard nuclear density for hydrogen. Figure 21 gives a summary of the approximate charge density distributions studied by R. Hofstadter in electron scattering off various nuclei.
Obviously, we understand that the Thomas-Fermi approximation which is used for evaluating becomes hardly justified at small number of quarks, and we should deal with the solutions of the complete equation system (34). However, one very encouraging hint comes from the chiral soliton model of nucleon [64], where it has been demonstrated that solving this system (34) the good description of nucleon and can be obtained.
In a sense we consider analyzing just a three quark system as a central result of our paper. Looking at Fig. 1 of Ref. [64] (we represent it here as Fig. 22) we see that the curve describing the behavior of the scalar field at large distances reaches its minimal value (according to the sign choice done in [64] it corresponds to the largest quark mass of order 300 MeV). It looks like it that by approaching the center of baryon, the chiral symmetry is partially restored and the scalar field in the region of ∼ 0.5 fm disappears. One of the possible scenarios for solving the system of equations (34) could be a variant in which the scalar field reaches maximal (zero) value (with a zero value of the derivative over coordinate) at this point (or the center of the baryon). Then a scalar field can, in principle, smoothly approach its minimal value coming to the center of the baryon. It allows us to conclude that we could deal with an "ordinary" quark with positive (zero) mass for the solutions of such a type. However, the baryon is getting a large width (size) in this scenario. There is another type of solutions, in which the "speed" of passing by the point 0.5 fm is not getting slower. In fact such a situation could be realized by doing a chiral rotation where a quark inside a baryon falls in the metastable region of negative quark masses. Such a solution develops an already quite suitable width of order ∼1 fm due to the presence of a massive (1 GeV) scalar field. Clearly, the problem of the existence of σ meson so heavy (strengthening the chiral effect) is crucial to collect the necessary information on the phase diagram of strongly interacting matter. Such solutions develop a surface tension coefficient which is larger by a factor 2 than the corresponding coefficient of a single kink and, as we believe, signals some instability of the single kink solution. Figure 23 displays the specific binding energy of the ensemble. It is defined by the expression similar to Eq. (39) in that the integration over the quark droplet volume is per- Fig. 23 The specific binding energy at N f = 1 and N f = 2 in MeV as a function of quark number N q formed. The specific energy is normalized (compared) to the ensemble energy at spatial infinity, i.e. in vacuum. Actually, Fig. 23 shows several curves in the upper part of the plot which corresponds to the calculations with N f = 1. The solid line is obtained by scanning over the parameter η and corresponds to the data presented in Table 1. The dashed curve is calculated at fixed η = 0.4 but by scanning over the parameter * M. It is clearly seen that if the specific energy data are presented as a function of quark number N q , then the solutions, which we are interested in, rally in the local vicinity of the curve where the maximal binding energy-|E b | is reached. A similar solution scanning can be performed over the central density parameter ρ 0 in origin. The corresponding data are dotted for a certain fixed * M and ρ 0 . It is interesting to notice that on scanning over any variable discussed a saturation property is observed and it looks like the minimum in e b at N q ∼ 200-250. The results for the specific binding energy as a function of particle number are in qualitative agreement with the corresponding experimental data. One may even observe quantitative agreement if the factor 3 (the energy necessary to remove one baryon) is taken into account. In fact, the equation system (34) represents an equation of balance for the current quarks circulating between liquid and gas phases.
Summarizing we would like to emphasize that in the present paper we have demonstrated how a phase transition of the liquid-gas kind (with reasonable values of parameters) emerges in the NJL-type models. The constructed quark ensemble displays some interesting features of the nuclear ground state (for example, the existence of the state degenerate with the vacuum one), and the results of our study are suggestive of speculating that the quark droplets could coexist in equilibrium with the vacuum under the normal conditions. It seems that if the quark drops really exist, then their features could be quite similar to the nuclear matter, and the scenario of partial chiral symmetry restoration has been realized as a chiral soliton. In the mean field (Ginsburg-Landau) approximation we are able to show the strengthening of chiral effects in the baryon that appears as the existence of a massive 1 GeV scalar meson, and a realistic phase diagram of strongly interacting matter may be obtained by studying the mechanism of chiral effects strengthening.

Conclusion
In the present paper we described quantum liquids (Landau Fermi liquids) resulting from the quark models with a four-fermion interaction. This consideration is based on the identity of results obtained in [21][22][23] by using a dressing Bogolyubov transformation and mean field approximation. We demonstrated that the mean energy of the ensemble serves as an energy functional of the Landau theory. It was shown that in a wide range of potentials interesting for applications one can expect the quantum liquids to behave in essentially the same way. For some of their properties a band of estimates was obtained. A comparison of NJL and KKB models, substantially different in many aspects, demonstrates that the properties of quantum liquids do not actually depend on the shape of the form factor (a natural interaction length); rather, they are mainly determined by the coupling constant of the interaction. It was shown that a common distinctive feature of ensembles is the presence of occupied states degenerate with respect to the vacuum in the chemical potential and pressure. Based on this observation, the inhomogeneous states, which allowed describing a transition layer, estimating the surface tension, as well as studying some properties of quark liquid droplets, were considered. It is noted that in the case of a small number of quarks in a droplet instability associated with lowering of the energy barrier, separating chiral phases, apparently manifests itself. This instability is seen in two kinks merging into one chiral soliton. The idea of a dynamical equilibrium of a mixed phase consisting of baryon matter and vacuum is discussed as a possible scenario for explaining the stability of nuclear matter.