text:
# Strange quark production in a statistical effective model

## I INTRODUCTION

The recent observation that hadron multiplicities in e+e$`^{}`$and hadronic high energy collisions agree very well with a statistical- thermodynamical ansatz beca1 ; beca2  has revived the interest on the statistical models in high energy collisions, an idea dating back to the 50’s fermi ; hagedorn . As the apparent chemical equilibrium of all hadron species in these elementary collisions (EC) cannot be driven by inelastic collisions amongst hadrons after their formation, this finding has led to the idea of a pure statistical filling of multi-hadronic phase space of homogeneous hadronizing regions (clusters or fireballs) as an intrinsic feature of hadronization process, occurring at a critical value of energy density beca2 ; heinz ; stock . Otherwise stated, hadrons are born in an equilibrium state, as envisaged by Hagedorn hagedorn2 . It must be emphasized that, within this framework, temperature and other thermodynamical quantities have an essential statistical meaning which does not imply the occurrence of a thermalization process via inelastic hadronic collisions.

The same model, in different versions, has been successfully applied to a large set of heavy ion collisions (HIC) data hions ; becahions2 ; becahions  and this has been interpreted as a clue of the minor effect of post-hadronization inelastic rescattering heinz , an indication which is also supported by kinetic models calculations kinetic .

One of the nice features of the statistical hadronization model (SHM) is the very low number of free parameters required to reproduce a large number of hadronic multiplicities. Provided that the masses and charges of the assumed hadronizing clusters fluctuate according to a particular shape function beca2 ; becapt , all Lorentz-invariant quantities (like e.g. hadron average multiplicities) depend on two parameters: the sum $`V`$ of the volumes of the clusters and the temperature $`T`$. However, in order to get a satisfactory agreement with the data, the model has to be supplemented, both in EC and HIC, with one more phenomenological parameter, $`\gamma _S`$, suppressing the production of particles containing $`n`$ strange valence quarks by $`\gamma _S^n`$ with respect to the expected production in a fully equilibrated hadron gas 111Recently, a different parameterization of the extra strangeness suppression has been introduced for elementary collisions becapt ..

The behavior of $`\gamma _S`$ as a function of collision and center-of-mass energy was found to be rather odd as it turns out to be significantly higher in e+e$`^{}`$collisions than in hadronic collisions over a large energy range beca1 ; beca2  whereas it has a fairly stable value in heavy ion collisions becahions2 , which is approximately the same as in e+e$`^{}`$collisions. A clearer insight in the mechanism of strangeness production can be achieved by estimating the ratio between newly produced valence strange and up, down quark pairs, the so-called Wroblewski factor $`\lambda _S`$:

$$\lambda _S=\frac{2\mathrm{s}\overline{\mathrm{s}}}{\mathrm{u}\overline{\mathrm{u}}+\mathrm{d}\overline{\mathrm{d}}}$$
(1)
which shows a striking regularity, being in fact fairly constant in all kinds of elementary collisions becahions  over two orders of magnitude of center-of-mass energy and as twice as large in high energy heavy ion collisions (see Fig. 1). It should be stressed that this ratio is calculated by counting the primarily produced quark pairs, i.e. those belonging to directly emitted hadrons. As primary hadrons are not measurable, they must be calculated by using a model and $`\lambda _S`$ turns out to be a model dependent quantity especially with regard to the number of u and d quarks which significantly increases during the post-hadronization hadronic decay chain. However, if the number of measured species is large and the model accurately reproduces them (which is the case for the SHM), the thereby estimated $`\lambda _S`$ is expected to be reasonably close to the actual value. The need of an extra suppression of strangeness with respect to the full statistical equilibrium for the hadronic system urges the search for an explanation based on a more microscopic approach. This is indeed the main subject of the present work, where we try to calculate $`\lambda _S`$ within one statistical model scheme in both EC and HIC, employing a constituent quark model as the basic underlying structure, with the purpose of justifying the lack of strangeness chemical equilibrium at hadron level observed in the data. In fact, constituent quark models have already been used to calculate strangeness production in HIC on the basis of transport equation aichreb . Amongst other microscopic approaches to this problem and, more generally, to account for hadronization equilibrium features, a recent study has been performed in ref. dumitru  focussed on the role of massive Polyakov loops.

The paper is organized as follows: in Sect. II we introduce the physical picture and in Sect. III the full model is described. In Sect. IV the Nambu-Jona-Lasinio model with exact conservation of quantum charges is used to calculate $`\lambda _S`$ and compare it with the corresponding value obtained from fits of the statistical hadronization model to the data in both HIC and EC. Also, the stability of the results is addressed. Sect. V is devoted to conclusions.

## II THE PHYSICAL PICTURE

The physical picture keeps the same scheme of the SHM: formation of a set of clusters endowed with charge, momentum, mass and volume, in local statistical equilibrium. However, statistical equilibrium is now assumed to apply to a system of constituent quarks whilst hadrons are assumed to be produced via their coalescence, still in a purely statistical fashion (statistical coalescence) so as to give rise to a partially chemically equilibrated hadron gas. Furthermore, each cluster is required to be a color singlet. This requirement has of course no effect at hadron level, but it is crucial in a quark model.

Two more crucial assumptions are introduced. First, the thermodynamical parameters fitted in the SHM, in high energy EC and HIC, are interpreted as the critical values for color deconfinement and, consequently, for (approximate) chiral symmetry restoration heuristic ; karsch ; secondly, it is assumed that the produced s and light quarks, or at least the ratio s/u, i.e. $`\lambda _S`$, survive in the final hadrons after hadronization has taken place.

The first assumption is supported by the constancy of the fitted temperatures in the SHM for various processes in EC becabiele  and by its agreement with other estimates of the critical temperature; in view of this fact, the use of effective models such as the Nambu-Jona-Lasinio model njlvero ; kunihiro  and others bcd ; tricritical ; conmassa ; giap  embodying chiral symmetry breaking and restoration, looks well suited. It should be pointed out that the values of the thermodynamical parameters (i.e. temperature and baryon-chemical potentials) extracted within the SHM in HIC (see Table 1) lie in a range where all effective models and recent lattice simulations fodor  predict smooth cross-over transitions, far from a possible critical ending point. For the same reason, extensions of the NJL model klevalast ; aichelin  at high baryon chemical potential and low temperature, where color superconductivity takes placecfl , can be disregarded.

The second assumption allows us to calculate $`\lambda _S`$ within a quark model and compare it with that obtained from valence quark counting in the produced hadrons. Also, it implies that the lack of full chemical equilibrium at hadron level is indeed the consequence of full chemical equilibrium at constituent quark level.

In the grand-canonical framework, which applies to a large system, $`\lambda _S`$ would depend only on intensive quantities such as chemical potentials and temperature and its constancy would be tightly related to theirs. However, if the system is not large, the canonical (and possibly micro-canonical) ensemble, in which the exact conservation of cluster quantum numbers in the hadronic or quark system is enforced, must be used and, as a consequence, $`\lambda _S`$ gets a dependence also on the volume (the so-called canonical suppression). While the volume within which quantum charges (baryon number, electric charge and strangeness) are fixed can be taken as the sum of the proper volumes of all clusters, though under appropriate non-trivial assumptions beca2 ; becapt , the volume over which color must be neutralized is, as has been mentioned, that of a single hadronizing cluster. This very fact introduces a further parameter in the model, i.e. the average single cluster volume $`V_c`$ which can indeed be much smaller than the volume over which flavor is conserved. As quarks do carry color charge, $`\lambda _S`$ might be significantly affected by variations of $`V_c`$.

Therefore, the underlying idea is that the presence of a characteristic constant value of $`V_c`$ in EC, related to the typical distance over which color is neutralized, may sizeably reduce the value of $`\lambda _S`$ (together with flavor constraint) whereas color deconfinement in HIC may lead to an enhancement of s/u ratio. We thus argue that the main difference between hadronization in EC and HIC is to be found in the typical size of the pre-hadronization color-neutral region: it has to be something of the order of a hadron radius in EC (a sort of mini-QGP droplet) and much larger in HIC if a macroscopic (i.e. extending over many hadron volumes) QGP has been formed. Of course, it must be emphasized that the process stage at which the statistical equilibrium is achieved is expected to be deeply different in EC as compared with HIC: an early thermalization of partons is indeed envisaged in HIC, which is maintained until hadronization, whereas a late, pre-hadronization local equilibrium scenario, possibly driven by the strongly coupled non-linear evolution of QCD fields in the late soft regime is envisaged for EC.

Hence, we argue that relative strange quark production might probe these two scenarios because of a possible stronger color and flavor-canonical suppression of s quarks with respect to u and d quarks in EC which is absent in HIC due to much larger single-cluster volume and overall volume. The ultimate reason of the stronger canonical suppression of s quarks in comparison with u, d, resides in their different constituent mass values. The hereby proposed mechanism of strangeness enhancement in HIC is rather different from that suggested by Muller and Rafelski rafe  insofar as no difference in the time scale is invoked but an equilibrium situation is assumed in both EC and HIC, though over differently sized regions.

## III THE MODEL

According to the basic idea of the SHM and what has been said in the previous section, the formation of a set of pre-hadronic, color-singlet clusters at statistical equilibrium, with definite values of flavor quantum numbers is envisaged as the result of a complex dynamical evolution. In an equilibrium scheme, all physical observables relevant to a particular cluster can be calculated by means of suitable operations on its associated partition function. Provided that cluster masses are large enough so that micro-canonical effects can be neglected, the requirement of definite quantum numbers implies that the $`i^{\mathrm{th}}`$ cluster’s partition function has to be the canonical one rather than the more familiar grand-canonical:

$$Z_i=\underset{\mathrm{states}}{}e^{E/T}$$
(2)
where $`T`$ is the cluster’s temperature, $`E`$ its energy and the sum is meant to run over all allowed states, namely the states with flavor and color quantum numbers matching those of the cluster. The calculation of $`Z_i`$ can be done by using the symmetry group $`G`$ associated to the involved quantum numbers. If the continuous group $`G`$ is an exact internal symmetry for the Hamiltonian $`H`$ and the cluster state belongs to its irreducible representation $`\nu `$, then exact :

$$Z_i=𝑑\mu (\gamma _1,\mathrm{},\gamma _r)\chi _\nu ^{}(\gamma _1,\mathrm{},\gamma _r)\widehat{Z}_i(\gamma _1,\mathrm{},\gamma _r)=𝑑\mu (\gamma _1,\mathrm{},\gamma _r)\chi _\nu ^{}(\gamma _1,\mathrm{},\gamma _r)\mathrm{Tr}\left(e^{\beta H\mathrm{i}{\displaystyle \underset{l=1}{\overset{r}{}}}\gamma _lQ_l}\right)$$
(3)
where $`d\mu (\gamma _1,\mathrm{},\gamma _r)`$ is the invariant normalized measure of $`G`$, $`\chi _\nu `$ is the character of the representation $`\nu `$, $`r`$ is the rank of $`G`$, $`Q_l`$ are the mutually commutating generators of the Cartan subalgebra and $`\gamma _l`$ are the group parameters. The actual group of interest is SU(3)$`\times `$U(1)3, where each U(1) corresponds to conserved net flavors $`U`$, $`D`$ and $`S`$ and SU(3) to the color. For instance, the canonical partition function of the $`i^{\mathrm{th}}`$ cluster of a gas of free quarks (no antiquarks included) with flavor $`j`$ involves the following $`\widehat{Z}_i`$ in Eq. (3):

$$\widehat{Z}_i(\theta _1,\theta _2,\varphi _j)=\mathrm{exp}\left(\underset{n=1}{\overset{\mathrm{}}{}}\frac{(1)^{n+1}}{n}\chi _{1,0}(n\theta _1,n\theta _2)e^{\mathrm{i}n\varphi _j}\frac{gV}{(2\pi )^3}d^3pe^{{\displaystyle \frac{nϵ_j}{T}}}\right)$$
(4)
where:

$$\chi _{1,0}(\theta _1,\theta _2)=e^{\mathrm{i}\theta _1}+e^{\mathrm{i}\theta _2}+e^{\mathrm{i}\left(\theta _1+\theta _2\right)}$$
(5)
is the character of the fundamental representation of SU(3), with $`\theta _{1,2}[\pi ,\pi ]`$; $`\varphi _j`$ is the parameter of U(1)j $`V`$ is the volume, $`g`$ is the spin degeneracy factor, $`ϵ_j=\sqrt{p^2+m_j^2}`$ and the sum over discrete states $`k`$ has been approximated with its continuum limit. When including antiquarks and considering three flavors (i.e. the full group SU(3)$`\times `$Uu(1)$`\times `$Ud(1)$`\times `$Us(1)), $`\widehat{Z}^i`$ can be written as the product of three functions like those in Eq. (4):

$$\widehat{Z}_i(\theta _1,\theta _2,\varphi _u,\varphi _d,\varphi _s)=\underset{j=u,d,s}{}\widehat{Z}_{ij}(\theta _1,\theta _2,\varphi _j)\widehat{Z}_{ij}^{}(\theta _1,\theta _2,\varphi _j)$$
(6)
The irreducible representation $`\nu `$ can be labelled by three integer numbers, the net flavor $`U_i,D_i,S_i`$, and, as far as the SU(3) color group is concerned, $`\nu `$ is simply the singlet representation. Therefore, by introducing the vectors $`𝐐_i=(U_i,D_i,S_i)`$ and $`\mathit{\varphi }_i=(\varphi _{iU},\varphi _{iD},\varphi _{iS})`$, the character can be written as $`\chi _\nu =\mathrm{exp}[\mathrm{i}𝐐_i\mathit{\varphi }_i]`$ and:

$$Z_{ij}=\left[\underset{j=1}{\overset{3}{}}_\pi ^\pi \frac{d\varphi _j}{2\pi }\right]𝑑\mu (\theta _1,\theta _2)e^{\mathrm{i}𝐐_i\mathit{\varphi }_i}\widehat{Z}_i(\theta _1,\theta _2,\mathit{\varphi }_i)$$
(7)
where:

$$𝑑\mu (\theta _1,\theta _2)=\frac{1}{3!}_\pi ^\pi \frac{d\theta _1}{2\pi }_\pi ^\pi \frac{d\theta _2}{2\pi }\underset{j<k}{\overset{3}{}}4\mathrm{sin}^2\frac{\theta _j\theta _k}{2}$$
(8)
with the constraint $`_{k=1}^3\theta _k=0\mathrm{mod}2\pi `$, is the SU(3) invariant integration. The partition function hitherto considered is relevant to one cluster. Yet, there are several clusters in a single collision event and, as long as global observables are concerned, theoretical predictions require summation over all clusters. Moreover, clusters may well be produced with different configurations of flavor numbers $`(U_i,D_i,S_i)`$ in different events (provided that their sum fulfills the conservation law, i.e. it must be equal to the flavor numbers of the colliding particles), hence an integration over all possible configurations weighted by its probability must be carried out in order to calculate averages of physical quantities. However, the probability distribution of clusters flavor configuration is an unknown function, which cannot be predicted within the statistical model and which is most likely governed by the preceding dynamical process. We will therefore assume a configuration probability distribution $`w`$ which is the “maximum disorder” one and which allows a remarkable simplification of the expressions of averages quark multiplicities as well as any other Lorentz-invariant observables. This function $`w`$, for an event with $`N`$ clusters can be written as:

$$w(𝐐_1,\mathrm{},𝐐_N)=\frac{_{i=1}^NZ_i(𝐐_i)\delta _{𝐐,\mathrm{\Sigma }_i𝐐_i}}{_{𝐐_1,\mathrm{},𝐐_N}_{i=1}^NZ_i(𝐐_i)\delta _{𝐐,\mathrm{\Sigma }_i𝐐_i}}$$
(9)
where $`Z_i`$ as in Eq. (7), the $`\delta `$ ensures that global conservation of initial quantum numbers $`𝐐`$ is fulfilled. If one has to calculate the average multiplicity of a given quark species $`j`$ in a $`N`$ cluster event, it is advantageous to introduce a fictitious fugacity $`\lambda _j`$ and taking the derivative of $`\mathrm{log}Z_i`$ with respect to it, so that the overall multiplicity reads:

$`n_j`$ $`=`$ $`{\displaystyle \underset{𝐐_1,\mathrm{},𝐐_N}{}}w(𝐐_1,\mathrm{},𝐐_N){\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \frac{}{\lambda _j}}\mathrm{log}Z_i(𝐐_i)|_{\lambda _j=1}=`$ (10)
$`=`$ $`{\displaystyle \underset{𝐐_1,\mathrm{},𝐐_N}{}}w(𝐐_1,\mathrm{},𝐐_N){\displaystyle \frac{}{\lambda _j}}{\displaystyle \underset{i=1}{\overset{N}{}}}\mathrm{log}Z_i(𝐐_i)|_{\lambda _j=1}`$
and, by using Eq. (9):

$$n_j=\frac{}{\lambda _j}\mathrm{log}\underset{𝐐_1,\mathrm{},𝐐_N}{}\underset{i=1}{\overset{N}{}}Z_i(𝐐_i)\delta _{𝐐,\mathrm{\Sigma }_i𝐐_i}|_{\lambda _j=1}$$
(11)
Were not for the color-singlet constraint stuck to each $`Z_i`$, the sum over all configurations on the right hand side of Eq. (11) would be, for $`\lambda _j=1`$, the canonical partition function of a single cluster having as volume the sum of volumes of the individual clusters and with quantum number vector $`𝐐`$beca2 , provided that all clusters have the same temperature. Thereby, also the explicit dependence on $`N`$ would vanish. However, this is no longer true once the color singlet constraint is introduced. For a similar strong equivalence with one cluster to apply, single clusters should be in suitable color quantum states superpositions, which is definitely against the common belief of pre-confinement and local color neutralization. Nevertheless, it can be proved that a considerable simplification in Eq. (11) occurs if clusters in one event are assumed to have the same volume, i.e. $`V_iV_c`$. In that case, it can be proved (see Appendix A) that the right hand side of Eq. (11), which can be defined as the global partition function, reads:

$$Z=\left[\underset{j=1}{\overset{3}{}}_\pi ^\pi \frac{d\varphi _j}{2\pi }\right]e^{\mathrm{i}𝐐\mathit{\varphi }}\left[𝑑\mu (\theta _1,\theta _2)\widehat{Z}_i(\theta _1,\theta _2,\mathit{\varphi })\right]^{V/V_c}$$
(12)
where $`V=_iV_i=NV_c`$ (for the free quarks gas $`\widehat{Z}_i`$ is given by Eq. (6) with the single flavor term $`\widehat{Z}_{ij}(\theta _1,\theta _2,\varphi _j)`$ calculated from Eq. (4) with (5)). The introduction of the local color-singlet constraint gets the value of the global partition function severely reduced with respect to an only globally color-singlet constrained case where clusters are allowed to be in a quantum superposition of colored states. This can be quite easily understood for any extra constraint essentially involves a decrease of the number of available states in the system. If $`V_cV`$ and $`V`$ is very large, i.e. in the thermodynamical limit, it can be shown that the usual grand-canonical expressions of the partition function, average multiplicities etc. are recovered.

Besides the basic assumption of local statistical equilibrium (for a single cluster), the expression (12) ultimately relies on a particular choice of flavor distribution (i.e. (9)) among clusters, which is indeed a non-trivial assumption. This particular distribution allows the reduction to one equivalent global cluster just because it is the “maximum disorder” distribution, namely the probability of obtaining a configuration $`(𝐐_1,\mathrm{},𝐐_N)`$ by randomly splitting a statistically equilibrated global cluster into $`N`$ sub-cluster each with volume $`V_c`$ and in a color singlet state ; in fact, this can be proved by showing that the $`w`$’s in Eq. (9) minimize the free energy of the system (see Appendix B).

It should pointed out that hadronizing clusters might be too small for a local (i.e. in each of them) temperature to be defined, as we have tacitly assumed. Should this be the case, the appropriate treatment would be micro-canonical and not canonical, with clusters described by mass and volume instead of temperature and volume. Notwithstanding, as far as the calculation of Lorentz-invariants is concerned, a global temperature may still be recovered by resorting to a similar equivalent global cluster reduction procedure at micro-canonical level, which is described in detail in ref. becapt  222Therein, the proof has been given for a hadron gas but it can be easily extended to the present case., with a micro-canonical equivalent of the distribution (9). The proof is quite lengthy and the mathematical framework more involved, but it is no longer necessary to assume that each cluster ought to have the same temperature. In fact, only the equivalent global cluster or, equivalently, the global canonical partition function (12), must be large enough to allow a canonical treatment of the system as a whole.

We now have to face the problem of a suitable choice of a microscopic Hamiltonian to be used in Eq. (3) and, thereby, make a theoretical calculation of $`\lambda _S`$. As full QCD is out of question, we resort to QCD-inspired low energy models. The simplest model to start with is the Nambu-Jona-Lasinio (NJL) model, which has quarks as fundamental degrees of freedom and embodies essential features of dynamical chiral symmetry breaking ($`\chi `$SB) in the limit of vanishing current quark masses. Also, we confine ourselves to mean-field approximation, as presented in ref. kunihiro , using the effective Hamiltonian:

$$H_{MFA}^{^{NJL}}=V\left[g__S(\alpha ^2+\beta ^2+\gamma ^2)+4g__D\alpha \beta \gamma \right]+d^3x\overline{q}\left(\mathrm{i}\gamma +𝐌\right)q$$
(13)
where $`V`$ is the spatial volume, $`g__S`$ and $`g__D`$ are the four-fermion and the six-fermion $`U(1)_A`$-breaking couplings respectively kunihiro ; uauno , $`\alpha =\overline{u}u,\beta =\overline{d}d,\gamma =\overline{s}s`$ are the quark condensates and:

$`M_u`$ $`=`$ $`m_u2g__S\alpha 2g__D\beta \gamma `$
$`M_d`$ $`=`$ $`m_d2g__S\beta 2g__D\alpha \gamma `$
$`M_s`$ $`=`$ $`m_s2g__S\gamma 2g__D\alpha \beta `$ (14)
are the constituent quark masses which, owing to the contact four-fermion interaction, are linearly related to the current quark masses $`m_u,m_d,m_s`$. As far as the grand-canonical calculation is concerned, the variational parameters $`\alpha ,\beta ,\gamma `$ are determined by minimizing the effective potential kunihiro  obtained within standard methods cjt ; dolan . Of course, the same result can be readily obtained by calculating:

$`Z`$ $`=`$ $`\mathrm{Tr}e^{\beta \left(H{\displaystyle \underset{i}{}}\mu _j\widehat{N}_j\right)}`$ (15)
$`=`$ $`e^{\mathrm{\Gamma }|_{\mathrm{min}}}`$
with $`H=H_{MFA}^{^{NJL}}`$. In Eq. (15), $`\widehat{N}_j`$ are the conserved flavor operators and $`\mathrm{\Gamma }|_{\mathrm{min}}`$ is the effective action at the physical point. For constant fields $`\mathrm{\Gamma }|_{\mathrm{min}}=\beta V𝒱|_{\mathrm{min}}=\mathrm{log}Z`$, where $`𝒱`$ is the effective potential.
The quark over antiquark excess is obtained by taking the derivative of $`𝒱`$ with respect to the chemical potential for a given flavor:

$$n_j\overline{n}_j=\frac{𝒱|_{\mathrm{min}}}{\mu _j}$$
(16)
whereas the number of particles, in the grand-canonical ensemble, reads:

$$n_j=\frac{N_cV}{\pi ^2}_0^\mathrm{\Lambda }𝑑p\frac{p^2}{e^{\beta \left(\sqrt{p^2+M_j^2}\mu _j\right)}+1}$$
(17)
where $`\mathrm{\Lambda }`$ is an UV cutoff and $`N_c`$ the number of colors.

One of the remarkable features of the Hamiltonian $`H_{MFA}^{^{NJL}}`$ in Eq. (13) is that the previously obtained expressions in Eqs. (7), (12) for the free quark gas partition functions with exact flavor and color conservation, can be taken up by simply replacing the current masses with the constituent masses $`m_jM_j`$ and setting the cutoff $`\mathrm{\Lambda }`$ as upper bound for the momentum integration in Eq. (4).

The NJL model is the simplest choice in order to study strange quark production, yet it is not a compelling one. Actually, its validity is limited within a temperature range well below the UV cutoff $`\mathrm{\Lambda }`$ aichelin . However, we are reasonably confident that the main results found within this model would not essentially change when employing other effective models whose validity, in principle, extend to arbitrarily high temperatures. As an example, it is worth mentioning the model developed in refs. bcd ; tricritical ; conmassa , named in ref. kunihiro  as ladder-QCD. In this model the exchange of a gauge particle between quarks is considered, implying a momentum dependence of the self-energy and, consequently, an irreducibility to a simple model with free constituent quarks as in NJL. However, ladder-QCD is considerably different from NJL only in the UV regime, where a $`1/p^2`$ tail is added to the self-energy, which is otherwise constant in the small momentum regionbcd ; russi . The effect of the different behaviour in the UV regime on the integrated number of particles is expected to be rather small bcd .

## IV ANALYSIS AND RESULTS

The goal of the numerical analysis is the calculation of $`\lambda _S`$ in both HIC and EC, by performing the ratio between newly produced strange quark and u, d quark pairs in the NJL model. The main ingredient needed to calculate $`\lambda _S`$ are the constituent quark masses which, for a given temperature and baryon-chemical potential, have to be determined by the minimization of the free energy $`F`$ and thus depend on the parameters contained in the effective Hamiltonian (13), i.e. $`g_S`$, $`\mathrm{\Lambda }`$, $`g_D`$ and the current quark masses. Their values have been fitted at zero temperature and baryon density in ref. kunihiro  to many static and dynamical meson properties:

$`\widehat{m}`$ $`=`$ $`{\displaystyle \frac{m_u+m_d}{2}}=5.5\mathrm{MeV};m_s=135.7\mathrm{MeV}`$
$`\mathrm{\Lambda }`$ $`=`$ $`631.4\mathrm{MeV};g__S\mathrm{\Lambda }^2=3.67;g__D\mathrm{\Lambda }^5=9.29`$ (18)
Particularly for $`g_D`$, in ref. kunihiro  a temperature dependence is introduced through the phenomenological law:

$$g_D(T)=g_D(T=0)e^{\left(T/T_0\right)^2}$$
(19)
where $`T_0`$ is a further free parameter and $`g_D(T=0)`$ is quoted in Eq. (IV). Different values of $`T_0`$ give rise to different behaviors of the constituent u, d mass (see Fig. 2 later on) as a function of the temperature as well as different positions of the cross-over curve in the $`(\mu _B,T)`$ plane. According to the statement in Sect. II, we take the thermodynamical parameters $`T`$ and $`\mu _B`$ in Table 1 as an input, and therefore $`T_0`$ can be fixed by enforcing that the quark susceptibility for a light quark with mass $`\widehat{m}`$ in Eq. (IV):

$$\chi _m=\frac{\overline{u}u}{m}|_{m=\widehat{m}}$$
(20)
has a peak at the desired $`(\mu _B,T)`$ point. This has been done for the most accurately determined thermodynamical parameters, which are those in Pb-Pb collisions at $`\sqrt{s}_{NN}=17.3`$ GeV, yielding $`T_0=170`$ MeV (see Fig. 3). Thereby, all relevant parameters are now fixed. The only free parameter left is the single cluster proper volume $`V_c`$, which may vary as a function of center-of-mass energy and colliding system. The dependence of $`\lambda _S`$ on this parameter is a crucial issue to be studied because it can possibly establish a relationship between color deconfinement and strangeness enhancement.

### IV.1 Heavy Ions

In HIC, calculations are much simpler as the total volumes (namely the sum of proper volumes of all the clusters) turn out to be so large that the effect of the canonical suppression related to exact conservation of flavor quantum numbers can be disregarded keranen , provided that flavor quantum numbers are distributed among clusters according to the function (9). Indeed, it should be pointed out that different scenarios have been devised in which the volume within which strangeness exactly vanishes (strangeness correlation volume) is less than the total volume, but this approach cannot account for the low yield of $`\varphi `$ meson redlich , so we will stick to the identification between the strangeness correlation volume and the total volume. If $`V_c`$ is large enough, quark multiplicities can thus be estimated by means of Eq. (17), the grand-canonical formula, and $`\lambda _S`$ turns out to be independent of $`V`$. In fact, it must be emphasized that the total volumes fitted in the SHM for HIC, though affected by large errors, lie in a range (see Table 1) where canonical suppression is negligible and $`\lambda _S`$ is sensitive to $`V_c`$ if this is roughly below 10 fm3 (see Fig. 4). This little sensitivity to $`V_c`$ makes $`\lambda _S`$ certainly not a clear-cut probe for deconfinement, though hadronizing clusters as small as some fm3 can indeed be excluded.

Hence, by taking a full grand-canonical approach, the cross-over line for chiral symmetry breaking has been calculated by minimizing the light quark mass susceptibility $`\chi _m`$. As has been mentioned, by forcing the location of the minimum to coincide with the fitted $`T`$ and $`\mu _B`$ in Pb–Pb collisions at $`\sqrt{s}_{NN}=17.3`$ GeV, we have been able to set the parameter $`T_0`$ to 170 MeV (see Fig. 3). For this calculation, the u and d chemical potentials have been assumed to be equal, so that $`\chi _u`$ and $`\chi _d`$ coincide with $`\chi _m`$, whereas the strange quark chemical potential $`\mu _s`$ has been set to zero due to $`S=0`$ constraint. It must be noted that, in principle, $`\mu _u`$ and $`\mu _d`$ differ because:

$`\mu _u`$ $`=`$ $`{\displaystyle \frac{\mu _B}{3}}+{\displaystyle \frac{2}{3}}\mu _Q`$
$`\mu _d`$ $`=`$ $`{\displaystyle \frac{\mu _B}{3}}{\displaystyle \frac{1}{3}}\mu _Q`$
$`\mu _s`$ $`=`$ $`{\displaystyle \frac{\mu _B}{3}}{\displaystyle \frac{1}{3}}\mu _Q\mu _S`$ (21)
where $`\mu _Q`$ is the electrical chemical potential, which vanishes only in isospin symmetric systems. However, this term turns out to be negligible for all examined collisions (see Table 1) and we have simply taken $`\mu _u=\mu _d=\mu _B/3`$. Having set $`T_0`$, the cross-over curve can be predicted from the model and this is shown in Fig. 5 along with fitted SHM points which are in satisfactory agreement with the calculation.

It must be pointed out that in a grand-canonical finite system with fixed volume, the conservation of a definite initial values of baryon number, electric charge and strangeness necessarily leads to different values of the chemical potentials for the two different phases, such as the considered constituent quark phase and the hadron gas. Certainly, the volume is not supposed to be fixed in the transition, thus the assumption of equal $`\mu _B`$’s may be correct. On the other hand, once $`\mu _B`$ has been set, the constraints $`S=0`$ and $`Q/B=Z/A`$ lead to definite $`\mu _Q`$ and $`\mu _S`$ (beware the difference with $`\mu _s`$) values which, in principle, are different in the two phases. While retaining $`\mu _s=0`$, we have checked the accuracy of the assumption $`\mu _Q=0`$ by studying the effect on $`\lambda _S`$ of non-vanishing electrical chemical potentials in the quark phase. The procedure is as follows: first, we have calculated $`\lambda _S`$ and quark masses with the main assumption, i.e. $`\mu _Q=0`$; then, by using those masses, we have calculated $`\mu _Q`$ in the quark phase by enforcing $`Q/B=Z/A`$; finally, with the obtained $`\mu _Q`$, we have recalculated quark masses and the new $`\lambda _S`$. The difference between the $`\lambda _S`$’s calculated in the two ways ranges from 1% in Pb–Pb collisions at SPS energy up to 2.7% in Au–Au collisions at AGS energy and is therefore negligible throughout.

The parameter $`\lambda _S`$ has been calculated along the cross-over curve by using Eq. (17) and the comparison with the SHM fitted values is shown in Fig. 6. The main prediction of the NJL-based model is an increase of $`\lambda _S`$ for decreasing $`\sqrt{s}`$ which is driven by the increase of $`\mu _B`$ (see Table 1). This involves an enhancement of relative strange quark production with respect to u, d, the so-called Pauli blocking effect. On the other hand, the observed $`\lambda _S`$ does not keep growing but undergoes a dramatic drop at very small energies (see Fig. 1), a fact which can be possibly explained by the onset of a purely hadronic production mechanism at low (yet not easy to locate) center-of-mass energies with possible local strangeness conservation redlich . Looking at the deviations of the SHM-fitted central values of $`\lambda _S`$ in Au–Au and Si–Au from the theoretical curve in Fig. (6), it can be argued that this hadronic production mechanism takes over at energies as low as AGS’s. Certainly, error bars are large and a definite conclusion needs more precise data. It should also be pointed out that a mild strange canonical suppression sets in already in AGS Si–Au collisionsbecahions  which is not taken into account here. Another sizeable discrepancy between data and model shows up in Au–Au collisions at RHIC ($`\sqrt{s}=130`$ GeV). However, in this case, $`\lambda _S`$ has been determined by fitting ratios of hadronic yields measured at mid-rapidity, unlike all other quoted collisions for which full phase space multiplicities have been used. This method may have led to an overestimation of $`\lambda _S`$ if the midrapidity region is enriched in strangeness as observed in Pb–Pb collisions at SPSbecasqm .

It is now worth discussing the robustness of these results in some detail. The calculation of $`\lambda _S`$ based on Eq. (17) depends on constituent quark masses and the UV cutoff parameter $`\mathrm{\Lambda }`$ besides the thermodynamical parameters $`T`$ and $`\mu _B`$ which have been set by using hadronic fits. This expression of $`\lambda _S`$ is indeed a general one for effective models with four fermion interactions in the mean-field approximation. It is then worth studying $`\lambda _S`$ by taking the constituent masses as free parameters instead of fixing them by means of a particular effective model as we have done so far. Whilst $`M_{u,d}(\mu ,T)`$ is strongly dependent on $`T_0`$, $`M_s(\mu ,T)`$ is fairly independent of it (see Fig. 2). Therefore, as far as the finite temperature sector is concerned, one can fairly conclude that this NJL model has actually one free parameter, either $`T_0`$ or $`M_{u,d}(\mu _{Bc},T_c)`$. In other words, fixing $`T_0`$ amounts to set a definite value of the light quark mass at the critical point.
However, the resulting $`\lambda _S`$ value does not strongly depend on this parameter either, as shown in Fig. 7, where the region in the $`(M_u,M_s)`$ plane (with both $`\mathrm{\Lambda }=631.4`$ MeV and $`\mathrm{\Lambda }=+\mathrm{}`$) allowed by the fitted $`\lambda _S`$ values in Pb–Pb collisions at $`\sqrt{s}_{NN}=17.3`$ GeV is shown: the bands pattern clearly indicates that $`\lambda _S`$ has a strong dependence on $`M_s`$ and much milder on $`M_{u,d}`$. Furthermore, the $`M_s`$ value at the critical temperature is quite constrained by its corresponding value at zero temperature and does not undergo strong variations as a function of the temperature and $`T_0`$, as shown in Fig. (2), which is indeed a quite general feature of effective models conmassa . Summarizing, we can state that this NJL model yields a $`\lambda _S`$ value in good agreement with the data (as determined through the SHM) essentially because it has a fit value of the constituent strange quark mass at zero temperature and density (about 500 MeV), whilst the particular value of $`T_0`$ has much less impact on it (see also Fig. (7)). On the other hand, for fixed constituent quark masses, $`\lambda _S`$ is sizeably affected by the cutoff $`\mathrm{\Lambda }`$. The calculation with $`\mathrm{\Lambda }=+\mathrm{}`$ is meant to give some upper bound on its variation due to neglected UV contributions, though a re-analysis within a renormalizable model would be desirable to have a more accurate estimation.

### IV.2 Elementary collisions

In EC the fitted total volumes are small and the calculations have to be carried out within the canonical formalism. As has been mentioned in the previous section, the partition function to be used is that in Eq. (12) with constituent masses replacing current masses and a $`\mathrm{\Lambda }`$ cut-off as upper bound in the momentum integration in Eq. (6). Unlike in the grand-canonical ensemble, $`\lambda _S`$ depends on the volume and a thorough comparison with SHM fitted parameters is much more involved. This is true in many respects: first, the fitted total volumes in the SHM are subject to large errors and obtained for point-like hadrons, thus reasonably underestimated. Secondly, for finite volumes one should in principle refit $`T_0`$ so as to obtain a cross-over point for a definite total volume $`V`$ and cluster volume $`V_c`$, for each process, at the desired temperature. Finally, one should minimize the full expression of free energy $`F=T\mathrm{log}Z`$ (with Z as in (12)) to determine the constituent masses at finite volume. For the present, this calculation is not affordable as the minimization of functions involving five-dimensional numerical integration as in Eq. (12) at each step with the due accuracy, implies exceedingly high computing times with presently available computers. For this reason, we have carried out a simpler calculation, with $`T_0`$ kept to 170 MeV and constituent quark masses calculated in the model at the thermodynamic limit.

Yet, even in this approximated calculation, the effect of flavor and color conservation over finite volumes can be studied and, hopefully, it can be verified if this mechanism of canonical suppression actually implies a reduction of $`\lambda _S`$ with respect to the thermodynamic limit consistently with the data. In fact, any conservation law enforced on finite volumes imply a reduction of heavy charged particle multiplicity stronger than lighter particles’ so that a decrease of $`\lambda _S`$ for decreasing volumes is expected. The ultimate reason of this effect is the reduced energy expense needed, in a finite system, to compensate any charge unbalance with light particles in comparison with heavy particles.

We can see this mechanism at play in Fig. 8 where $`\lambda _S`$ is plotted as a function of the total volume $`V`$ (over which flavor is exactly conserved), for different (color-singlet) cluster volumes $`V_c`$ for an initially completely neutral system, such as e+e$`^{}`$. The temperature has been set to $`T=160`$ MeV, which is a fair average of the SHM fitted values beca2 ; becapt ; becabiele  in e+e$`^{}`$and the constituent quark masses have been set to $`M_s=452`$ MeV and $`M_{u,d}=112`$ MeV, which are the thermodynamic limit values. Interestingly, $`\lambda _S`$ features a fair stability over a large range of total volumes greater than $`20`$ fm3 for $`V_c<10`$ fm3 within the data band, with a mild increase up to an asymptotic value which turns out to be closer to the thermodynamic limit (yet not equal to it, see also Fig. (4)) for larger $`V_c`$. As the volumes fitted in a point-like hadron gas model in e+e$`^{}`$and $`\mathrm{p}\overline{\mathrm{p}}`$collisions turn out to be larger than 15 fm3, but not larger than $`80`$ fm3, we can argue that a suitable $`V_c`$ can account for the observed values and the stability of $`\lambda _S`$ over a reasonably large volume range, taking into account that the actual volumes are certainly larger than their point-like estimates. Moreover, the mechanism of local color neutrality induces the reduction of $`\lambda _S`$ value with respect to the thermodynamic limit of $`0.31`$ which is needed to quantitative reproduce the data. Within this approach, the difference between $`\lambda _S`$ in high energy HIC and EC would be mainly the result of the non vanishing initial baryon density (from $`0.31`$ to $`0.45`$) and of local color neutrality in EC over small regions ($`510`$ fm3).

The situation is somewhat different for pp collisions, due to non vanishing initial baryon number and electric charge. In this system, one can observe two compensating effects, namely Pauli exclusion principle favoring relative s$`\overline{\mathrm{s}}`$ production and canonical suppression, favoring relative u$`\overline{\mathrm{u}}`$ and d$`\overline{\mathrm{d}}`$ production. As a result, $`\lambda _S`$ is nearly constant around 0.3, in apparent disagreement with the observed value which is about 0.2 (see Fig. 9). The reduction of $`\lambda _S`$ brought about by the decrease of $`V_c`$ is not sufficient to restore the agreement even for very low $`V_c`$. The disagreement can be possibly explained by the inadequacy of the taken assumptions, such as the statistical distribution of charges among clusters according to Eq. (9).

## V CONCLUSIONS AND OUTLOOK

We have presented a statistical model to explain the observed pattern of strangeness production in both elementary and heavy ion collisions. The basic idea is that full chemical equilibrium is locally achieved at the level of constituent quark degrees of freedom (within the framework of a simple effective model), at a different stage of the evolution process in elementary collisions (late) with respect to heavy ion collisions (early). In this approach, hadron formation takes place through the coalescence of constituent quarks and this accounts for the known observation of an incomplete strangeness equilibrium at hadron level in the statistical hadronization models. The underlying assumption is that at least the ratio s/u survives in the final hadrons.

Besides the effect of different initial excess of matter over antimatter, the smaller relative strangeness production in EC with respect to HIC has been related to the smaller overall system size and to color confinement over much smaller regions. Particularly, the existence of a characteristic cluster volume, roughly between 5 and 10 fm3 and independent of center-of-mass energy, can possibly account for the observed stability of $`\lambda _S`$, along with the constancy of temperature and the weak dependence of $`\lambda _S`$ itself on the total volume $`V`$ (if not below $``$60 fm3), at least in e+e$`^{}`$and $`\mathrm{p}\overline{\mathrm{p}}`$collisions. However, the dependence of $`\lambda _S`$ on cluster volume in HIC is so mild that it is not possible to use strangeness production to prove the formation of large color-neutral regions, i.e. color deconfinement.

The numerical analysis, carried out within an effective Nambu-Jona-Lasinio model, involved only one parameter to be adjusted which has been used to get the critical point agreeing with that fitted in Pb–Pb collisions. A satisfactory agreement with the data has been found in several heavy ion collisions and e+e$`^{}`$, $`\mathrm{p}\overline{\mathrm{p}}`$collisions. Deviations have been found at low energy HIC, such as Au–Au and Si–Au at AGS, where the strangeness production mechanism could be predominantly driven by hadronic inelastic collisions. Furthermore, a significant discrepancy has been found in pp, which might be perhaps cured by taking a different scheme of quantum numbers distribution among the hadronizing clusters. It should be pointed out that the analysis in EC has been performed with some approximations as a full consistent treatment of the Nambu-Jona-Lasinio model for finite volumes is presently beyond our possibilities.

The model-dependence of our results has been discussed. We are reasonably confident that the obtained results are quite general and stable within non-renormalizable models with four-fermion interactions treated in the mean field approximation. The exploration of renormalizable effective model is under investigation.

## APPENDIX A

We want to prove that:

$$\underset{𝐐_1,\mathrm{},𝐐_N}{}\underset{i=1}{\overset{N}{}}Z_i(𝐐_i)\delta _{𝐐,\mathrm{\Sigma }_i𝐐_i}$$
(22)
is the global partition function in Eq. (12), with $`Z_i`$ as in Eq. (7), provided that all clusters have the same volume $`V_c`$. We will first prove it in the canonical framework, with each cluster having also the same temperature $`T`$. One can rewrite Eq. (22) by using Eq. (7) and the integral representation of the $`\delta _{𝐐,\mathrm{\Sigma }_i𝐐_i}`$ as:

$$\underset{𝐐_1,\mathrm{},𝐐_N}{}\left[\underset{j=1}{\overset{3}{}}_\pi ^\pi \frac{d\varphi _j}{2\pi }\right]e^{\mathrm{i}𝐐\mathit{\varphi }\mathrm{i}\mathrm{\Sigma }_i𝐐_i\mathit{\varphi }}\underset{i=1}{\overset{N}{}}\left[\underset{j=1}{\overset{3}{}}_\pi ^\pi \frac{d\varphi _{ji}}{2\pi }\right]e^{\mathrm{i}𝐐_i\mathit{\varphi }_i}𝑑\mu (\theta _1,\theta _2)e^{V_cf(\theta _1,\theta _2,\mathit{\varphi }_i)}$$
(23)
where $`\widehat{Z}_i`$ has been written in the simple form $`\mathrm{exp}[V_cf(\theta _1,\theta _2,\mathit{\varphi }_i)]`$ according to Eq. (6). It must be stressed that this expression of $`\widehat{Z}`$ is very general and does not depend on the Hamiltonian of the considered model. Moreover, the function $`f`$ is the same for all clusters for it depends only on the temperature which has been assumed to be the same. We can now perform the summation over all $`𝐐_i`$ in Eq. (22) and get:

$$\underset{𝐐_1,\mathrm{},𝐐_N}{}e^{\mathrm{i}\mathrm{\Sigma }_i𝐐_i\left(\mathit{\varphi }_i\mathit{\varphi }\right)}=\underset{i=1}{\overset{N}{}}(2\pi )^3\delta ^3(\mathit{\varphi }\mathit{\varphi }_i)$$
(24)
so that the integration over $`\varphi _{ji}`$ in Eq. (22) can be easily performed and one is left with:

$$\left[\underset{j=1}{\overset{3}{}}_\pi ^\pi \frac{d\varphi _j}{2\pi }\right]e^{\mathrm{i}𝐐\mathit{\varphi }}\underset{i=1}{\overset{N}{}}𝑑\mu (\theta _1,\theta _2)e^{V_cf(\theta _1,\theta _2,\mathit{\varphi })}$$
(25)
The integral over SU(3) group is the same for all clusters and Eq. (25) can then be written as:

$`\left[{\displaystyle \underset{j=1}{\overset{3}{}}}{\displaystyle _\pi ^\pi }{\displaystyle \frac{d\varphi _j}{2\pi }}\right]e^{\mathrm{i}𝐐\mathit{\varphi }}\left[{\displaystyle 𝑑\mu (\theta _1,\theta _2)e^{V_cf(\theta _1,\theta _2,\mathit{\varphi })}}\right]^N`$
$`=`$ $`\left[{\displaystyle \underset{j=1}{\overset{3}{}}}{\displaystyle _\pi ^\pi }{\displaystyle \frac{d\varphi _j}{2\pi }}\right]e^{\mathrm{i}𝐐\mathit{\varphi }}\left[{\displaystyle 𝑑\mu (\theta _1,\theta _2)\widehat{Z}_i(\theta _1,\theta _2,\mathit{\varphi })}\right]^{V/V_c}`$ (26)
where $`V=\mathrm{\Sigma }_iV_i=NV_c`$, which is equal to the right hand side of Eq. (12).

## APPENDIX B

We show that the configuration probabilities $`w(𝐐_1,\mathrm{},𝐐_N)`$ in Eq. (9) minimize the free energy of a system with volume $`V`$ and temperature $`T`$ which is split into $`N`$ color-singlet clusters with volumes $`V_1,\mathrm{},V_N`$ such that $`_iV_i=V`$. Let $`p`$ be the full probability of a single state in this system and $`w`$ the probability of a configuration $`(𝐐_1,\mathrm{},𝐐_N)`$. Then:

$$p_{\mathrm{state}}=w(𝐐_1,\mathrm{},𝐐_N)\underset{i=1}{\overset{N}{}}\frac{\mathrm{exp}(E_i/T)}{Z_i(𝐐_i)}\delta _{𝐐_i,𝐐_{i;\mathrm{s}}}\delta _{\mathrm{singlet}_\mathrm{i}}$$
(27)
where $`Z_i`$ is given by Eq. (7), $`E_i`$ is the energy of a cluster, $`𝐐_{i;\mathrm{s}}`$ is the vector of flavor quantum numbers of the state and $`\delta _{\mathrm{singlet}_\mathrm{i}}`$ signifies the color singlet constraint on each cluster. The entropy reads:

$`S`$ $`=`$ $`{\displaystyle \underset{\mathrm{states}}{}}p\mathrm{log}p`$
$`=`$ $`{\displaystyle \underset{𝐐_1,\mathrm{},𝐐_N}{}}{\displaystyle \underset{\mathrm{states}_1}{}}\mathrm{}{\displaystyle \underset{\mathrm{states}_\mathrm{N}}{}}w(𝐐_1,\mathrm{},𝐐_N){\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \frac{e^{E_i/T}}{Z_i(𝐐_i)}}\mathrm{log}w(𝐐_1,\mathrm{},𝐐_N)\left[{\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \frac{\mathrm{exp}(E_i/T)}{Z_i(𝐐_i)}}\right]`$
where the flavor and color constraint on each cluster are now implied in the sum over the states.

Eq. (APPENDIX B) can be worked out:

$`S`$ $`=`$ $`{\displaystyle \underset{𝐐_1,\mathrm{},𝐐_N}{}}{\displaystyle \underset{\mathrm{states}_1}{}}\mathrm{}{\displaystyle \underset{\mathrm{states}_\mathrm{N}}{}}w(𝐐_1,\mathrm{},𝐐_N){\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \frac{e^{E_i/T}}{Z_i(𝐐_i)}}`$
$`\times `$ $`\left[\mathrm{log}w(𝐐_1,\mathrm{},𝐐_N){\displaystyle \frac{E_i}{T}}\mathrm{log}{\displaystyle \underset{i=1}{\overset{N}{}}}Z_i(𝐐_i)\right]`$
$`=`$ $`{\displaystyle \underset{𝐐_1,\mathrm{},𝐐_N}{}}w(𝐐_1,\mathrm{},𝐐_N)\left[\mathrm{log}w(𝐐_1,\mathrm{},𝐐_N)+{\displaystyle \frac{_{i=1}^NE_i_{(𝐐_1,\mathrm{},𝐐_N)}}{T}}+\mathrm{log}{\displaystyle \underset{i=1}{\overset{N}{}}}Z_i(𝐐_i)\right]`$
where $`_iE_i_{(𝐐_1,\mathrm{},𝐐_N)}`$ is the average energy for the $`i^{\mathrm{th}}`$ cluster and a fixed configuration $`(𝐐_1,\mathrm{},𝐐_N)`$. The total average energy of the system, to be identified with the internal energy $`U`$, evidently reads:

$$U=\underset{𝐐_1,\mathrm{},𝐐_N}{}w(𝐐_1,\mathrm{},𝐐_N)\underset{i=1}{\overset{N}{}}E_i_{(𝐐_1,\mathrm{},𝐐_N)}$$
(30)
so that entropy can be rewritten as:

$$S=\underset{𝐐_1,\mathrm{},𝐐_N}{}w(𝐐_1,\mathrm{},𝐐_N)\mathrm{log}w(𝐐_1,\mathrm{},𝐐_N)+\frac{U}{T}+\underset{𝐐_1,\mathrm{},𝐐_N}{}w(𝐐_1,\mathrm{},𝐐_N)\mathrm{log}\underset{i=1}{\overset{N}{}}Z_i(𝐐_i)$$
(31)
and the free energy $`F`$ as:

$$F=UTS=T\underset{𝐐_1,\mathrm{},𝐐_N}{}w(𝐐_1,\mathrm{},𝐐_N)\left[\mathrm{log}w(𝐐_1,\mathrm{},𝐐_N)\mathrm{log}\underset{i=1}{\overset{N}{}}Z_i(𝐐_i)\right]$$
(32)
Now we seek for the probabilities $`w`$ which minimize the free energy with the constraint $`w=1`$. Therefore we have to introduce a Lagrange multiplier $`\lambda `$ and look for the extremals of $`F+\lambda (w1)`$ with respect to each $`w(𝐐_1,\mathrm{},𝐐_N)`$ and $`\lambda `$. This leads to:

$$\mathrm{log}w(𝐐_1,\mathrm{},𝐐_N)=\mathrm{log}\underset{i=1}{\overset{N}{}}Z_i(𝐐_i)\frac{\lambda }{T}$$
(33)
and

$$w(𝐐_1,\mathrm{},𝐐_N)=\frac{_{i=1}^NZ_i(𝐐_i)\delta _{𝐐,\mathrm{\Sigma }_i𝐐_i}}{_{𝐐_1,\mathrm{},𝐐_N}_{i=1}^NZ_i(𝐐_i)\delta _{𝐐,\mathrm{\Sigma }_i𝐐_i}}$$
(34)
after proper normalization.

###### Acknowledgements.

We are grateful to J. Aichelin, A. Barducci, R. Casalbuoni, U. Heinz, J. Lenaghan, O. Scavenius for useful discussions.


format: text

text:
# An exact formula for general spectral correlation function of random Hermitian matrices

## 1 Introduction

It is a classical result of analysis known from the nineteenth century that orthogonal polynomials can be represented as multi-variable integrals. For example, let $`\pi _j(x)`$ denote a monic orthogonal polynomial of $`j^{\text{th}}`$ degree with respect to the measure $`d\mu (x)=\mathrm{exp}(NV(x))dx`$ on the real axis,

$$\pi _j(x)\pi _k(x)e^{NV(x)}𝑑x=c_jc_k\delta _{jk}$$
(1)
and $`\pi _j(x)=x^j+\text{lower}\text{degrees}`$. Then for the monic orthogonal polynomial $`\pi _N(x)`$ there exists an integral representation (see Szeg$`\ddot{\text{o}}`$ )

$$\pi _N(x)=\frac{1}{Z_N}\underset{j=1}{\overset{N}{}}(xx_j)e^{N\underset{j=1}{\overset{N}{}}V(x_j)}\mathrm{}^2(x_1,\mathrm{},x_N)dx_1\mathrm{}x_N$$
(2)
In the formula above $`\mathrm{}(x_1,\mathrm{},x_N)`$ stands for the Vandermonde determinant, and the constant $`Z_N`$ is given by the product $`Z_N=N!\underset{j=0}{\overset{N1}{}}c_j^2`$.

The classical integral representation of the monic orthogonal polynomial $`\pi _N(x)`$ given in (2) suggests a random matrix interpretation. Indeed if we consider an ensemble of random $`N\times N`$ Hermitian matrices $`H`$ with the joint probability density of eigenvalues (see Mehta )

$$𝒫^{(N)}(x_1,\mathrm{},x_N)=\frac{1}{Z_N}\mathrm{exp}\left\{N\underset{i=1}{\overset{N}{}}V(x_i)\right\}\mathrm{}^2(x_1,\mathrm{},x_N)$$
(3)
the monic orthogonal polynomial $`\pi _N(x)`$ can be understood as an average of the characteristic polynomial $`𝒵_N[x,H]=\text{det}(xH)`$ over the ensemble

$$\pi _N(x)=𝒵_N[x,H]_H$$
(4)
A natural question which arises at this point is the following: what do we obtain if instead of one characteristic polynomial we average a combination of characteristic polynomials (for example, their product). This question was recently addressed in Brezin and Hikami , Mehta and Normand  who derived a generalization of formula (4). The authors considered the ensemble average of products of characteristic polynomials. It was found that this average is expressed as a determinant:

$$\underset{j=1}{\overset{L}{}}𝒵_N[\lambda _j,H]_H=\frac{1}{\mathrm{}(\widehat{\lambda })}\text{det}\left|\begin{array}{cccc}\pi _N(\lambda _1)& \pi _{N+1}(\lambda _1)& \mathrm{}& \pi _{N+L1}(\lambda _1)\\ \pi _N(\lambda _2)& \pi _{N+1}(\lambda _2)& \mathrm{}& \pi _{N+L1}(\lambda _2)\\ \mathrm{}& & & \\ \pi _N(\lambda _L)& \pi _{N+1}(\lambda _L)& \mathrm{}& \pi _{N+L1}(\lambda _L)\end{array}\right|$$
(5)
where $`\mathrm{\Delta }(\widehat{\lambda })\mathrm{\Delta }(\lambda _1,\mathrm{},\lambda _L)`$. (See also a related paper by Forrester and Witte  which deals with a particular case of positive integer moments: $`\lambda _1=\mathrm{}=\lambda _L`$).

In the present paper we obtain a further generalization of formulas (4) and (5). Namely we compute the general spectral correlation function of characteristic polynomials defined as

$$𝒦_N(\widehat{ϵ},\widehat{\mu })=\frac{\underset{j=1}{\overset{L}{}}𝒵_N[\mu _j,H]}{\underset{j=1}{\overset{M}{}}𝒵_N[ϵ_j,H]}_H$$
(6)
Here the symbols $`\widehat{ϵ},\widehat{\mu }`$ denote the vectors with the components:

$`\widehat{ϵ}=(ϵ_1,\mathrm{},ϵ_M);\widehat{\mu }=(\mu _1,\mathrm{},\mu _L)`$ (7)
For the function $`𝒦_N(\widehat{ϵ},\widehat{\mu })`$ to be well defined we assume that $`\text{Im}ϵ_k0`$, and also assume $`NM,NL`$. Then the result we obtain for the correlation function $`𝒦_N(\widehat{ϵ},\widehat{\mu })`$ is

$$𝒦_N(\widehat{ϵ},\widehat{\mu })=\frac{_{j=NM}^{N1}\gamma _j}{\mathrm{}(\widehat{\mu })\mathrm{}(\widehat{ϵ})}\text{det}\left|\begin{array}{cccc}h_{NM}(ϵ_1)& h_{NM+1}(ϵ_1)& \mathrm{}& h_{N+L1}(ϵ_1)\\ \mathrm{}& & & \\ h_{NM}(ϵ_M)& h_{NM+1}(ϵ_M)& \mathrm{}& h_{N+L1}(ϵ_M)\\ \pi _{NM}(\mu _1)& \pi _{NM+1}(\mu _1)& \mathrm{}& \pi _{N+L1}(\mu _1)\\ \mathrm{}& & & \\ \pi _{NM}(\mu _L)& \pi _{NM+1}(\mu _L)& \mathrm{}& \pi _{N+L1}(\mu _L)\end{array}\right|$$
(8)
where $`h_k(ϵ)`$ denotes the Cauchy transform of the monic orthogonal polynomial $`\pi _k(x)`$,

$$h_k(ϵ)=\frac{1}{2\pi i}\frac{\pi _k(x)e^{NV(x)}dx}{xϵ},\text{Im}ϵ0$$
(9)
and

$$\gamma _{n1}=\frac{2\pi i}{c_{n1}^2}$$
(10)
The formula above is of interest since it shows the place of the correlation function $`𝒦_N(\widehat{ϵ},\widehat{\mu })`$ in the theory of orthogonal polynomials. Furthermore it provides new insight on $`𝒦_N(\widehat{ϵ},\widehat{\mu })`$, for example it relates $`𝒦_N(\widehat{ϵ},\widehat{\mu })`$ with Riemann-Hilbert problem for orthogonal polynomials (see section 5)

Correlation functions $`𝒦_N(\widehat{ϵ},\widehat{\mu })`$ are interesting objects themselves since they contain a very detailed information about spectra of random matrices. In particular, knowledge of such correlations enables to extract the n-point correlation function of spectral densities. Moreover, distributions of some interesting physical quantities such as, e.g. level curvatures are expressed in terms of $`𝒦_N(\widehat{ϵ};\widehat{\mu })`$, (see e.g  and Appendix A of  for more examples).

The investigation of correlation functions of characteristic polynomials and of their moments is also motivated by hope to relate statistics of zeroes of the Riemann zeta function to that of eigenvalues of large random matrices . The spectral determinants are also relevant for interesting combinatorial problems, see . Other numerous applications of spectral determinants are in the theory of quantum chaotic and disordered systems, and in quantum chromodynamics, see - for an extensive list of references.

There are several analytical techniques for dealing with the integer moments (positive or negative) of characteristic polynomials and more general correlations functions. Their applicability varies with the nature of the underlying random matrix ensemble. The particular case of the Gaussian measure - Gaussian Unitary Ensemble (GUE), as well as its chiral counterpart, chGUE - can be studied very efficiently by a modification of the standard supersymmetry technique  suggested recently by the authors . One starts, as usual, by representing each of the characteristic polynomials as a Gaussian integral over anti-commuting (Grassmann) variables and inverse characteristic polynomial by Gaussian integral over complex variables. This allows to average the resulting expressions straightforwardly. At the next step one employs the so-called Hubbard-Stratonovich transformation for ”fermionic” degrees of freedom, and then exploits the Itzykson-Zuber-Harish-Chandra integrals together with their natural non-compact extensions. A detailed account of the method and results as well as related references can be found in the recent papers -.

The resulting expressions in  revealed a very attractive determinantal structure which was especially evident in the case of chiral GUE (Laguerre ensemble) , see also . When one deals with the positive moments only such structures naturally arises in the framework of the orthogonal polynomial method . This fact is suggestive of the idea that the determinantal structure in the general case (i.e. when both products and ratios of characteristic polynomials are involved) should be valid for an arbitrary unitary-invariant potential. In present paper we show that this is indeed the case.

The paper is organized as follows. In section 2 we derive an algebraic identity which represents $`\underset{j=1}{\overset{M}{}}𝒵_N^1[ϵ_j,H]`$ as a sum over permutations. This identity together with formula (5) will permit us to rewrite $`𝒦_N(\widehat{ϵ},\widehat{\mu })`$ as an integral over $`M`$ variables with subsequent derivation of formula (8). Section 3 illustrates our approach on the simplest case $`M=1`$ and $`L=0`$. Here we derive a multi- variable integral representation for Cauchy transforms of monic orthogonal polynomials. The computation of the general correlation function $`𝒦_N(\widehat{ϵ},\widehat{\mu })`$ for any integers $`0L,MN`$ is given in section 4. In section 5 we establish relation with the Riemann-Hilbert problem for orthogonal polynomials proposed by Fokas, Its and Kitaev  (see also ), and outline a way to investigate the large $`N`$ asymptotic of $`𝒦_N(\widehat{ϵ},\widehat{\mu })`$ by Deift-Zhou steepest descent/stationary phase method . A detailed asymptotic analysis of $`𝒦_N(\widehat{ϵ},\widehat{\mu })`$ is beyond the scope of this paper and will be given in the subsequent publication . Finally in section 6 we briefly discuss some ”duality relations” for matrix integrals emerging for the case of the Gaussian potential.

## 2 Inverse of products of characteristic polynomials as a sum over permutations

Let $`MN`$, and $`x_1,\mathrm{},x_N`$ denote the eigenvalues of a matrix $`H`$ ($`\text{dim}H=N\times N`$). Then the following algebraic identity holds:

$`{\displaystyle \underset{l=1}{\overset{M}{}}}{\displaystyle \frac{ϵ_l^{NM}}{𝒵_N[ϵ_l,H]}}=`$ (11)
$`{\displaystyle \underset{\sigma 𝖲_N/𝖲_{NM}\times 𝖲_M}{}}\left({\displaystyle \frac{\underset{i=1}{\overset{M}{}}x_{\sigma (i)}^{NM}}{\underset{i,j=1}{\overset{M}{}}ϵ_jx_{\sigma (i)}}}\right){\displaystyle \frac{\mathrm{}(x_{\sigma (1)},\mathrm{},x_{\sigma (M)})\mathrm{}(x_{\sigma (M+1)},\mathrm{},x_{\sigma (N)})}{\mathrm{}(x_{\sigma (1)},\mathrm{},x_{\sigma (N)})}}`$
where $`𝖲_N`$ is the permutation group of the indices $`1,\mathrm{},N`$; $`𝖲_M`$ is the permutation group of the first $`M`$ indices and $`𝖲_{NM}`$ is the permutation group of the remaining $`NM`$ indices. Identity (11) follows as a consequence of the Cauchy-Littlewood formula

$$\underset{j=1}{\overset{M}{}}\underset{i=1}{\overset{N}{}}(1x_iy_j)^1=\underset{\lambda }{}s_\lambda (x_1,\mathrm{},x_N)s_\lambda (y_1,\mathrm{},y_M)$$
(12)
and the Jacobi-Trudi identity :

$$s_\lambda (x_1,\mathrm{},x_N)=\frac{\text{det}\left(x_i^{\lambda _jj+N}\right)}{\mathrm{}(x_1,\mathrm{},x_N)}$$
(13)
where the Schur polynomial $`s_\lambda (x_1,\mathrm{},x_N)`$ corresponds to the partition $`\lambda `$ and the indices $`i,j`$ take the values from $`1`$ to $`N`$. In order to prove (11) we rewrite the Schur polynomial $`s_\lambda (x_1,\mathrm{},x_N)`$ in (12) as

$$s_\lambda (x_1,\mathrm{},x_N)=\frac{\underset{\pi \text{S}_N}{}()^{\nu _\pi }x_{\pi (1)}^{\lambda _11+N}\mathrm{}x_{\pi (M)}^{\lambda _MM+N}x_{\pi (M+1)}^{M1+N}\mathrm{}x_{\pi (N)}^0}{\mathrm{}(x_1,\mathrm{},x_N)}$$
(14)
(As $`s_\lambda (y_1,\mathrm{},y_M)=0`$ for any partition with the number of rows larger then $`M`$ only the partitions with $`\lambda _{M+1}=\mathrm{}=\lambda _N=\mathrm{}=0`$ contribute in equation (12)). In equation (14) $`\nu _\pi =1`$ (or $`\nu _\pi =0`$) if the permutation $`\pi \text{S}_N`$ is odd (even).

Any permutation $`\pi \text{S}_N`$ in the above sum can be decomposed as a product of three subsequent permutations, i.e. $`\pi =\sigma \pi _2\pi _1`$. The first one, $`\pi _1\text{S}_M`$ is a permutation of elements of the set $`[1,2,\mathrm{},M]`$; the second permutation, $`\pi _2\text{S}_{NM}`$ is that of elements of the set $`[M+1,M+2,\mathrm{},N]`$; and the third permutation, $`\sigma \text{S}_N/\text{S}_M\times \text{S}_{NM}`$, is an exchange of elements between these two sets. For example, let us take $`M=3`$, $`N=5`$, and permutation $`\pi =\left(\begin{array}{ccccc}1& 2& 3& 4& 5\\ 1& 4& 2& 5& 3\end{array}\right)`$. This permutation can be represented as follows:

$`\pi =\sigma \pi _2\pi _1`$
$`\pi _1=\left(\begin{array}{ccccc}1& 2& 3& 4& 5\\ 1& 3& 2& 4& 5\end{array}\right),\pi _2=\left(\begin{array}{ccccc}1& 2& 3& 4& 5\\ 1& 2& 3& 5& 4\end{array}\right),\sigma =\left(\begin{array}{ccccc}1& 2& 3& 4& 5\\ 1& 5& 3& 4& 2\end{array}\right)`$ (21)
We now rewrite the denominator of the right-hand side of equation (14) as

$`\mathrm{}(x_1,\mathrm{},x_N)=()^{\nu _\pi }\mathrm{}(x_{\pi (1)},\mathrm{},x_{\pi (N)})=`$ (22)
$`()^{\nu _\pi }\mathrm{}(x_{\pi (1)},\mathrm{},x_{\pi (M)})\mathrm{}(x_{\pi (M+1)},\mathrm{},x_{\pi (N)}){\displaystyle \underset{i[1,\mathrm{},M],j[M+1,\mathrm{},N]}{}}(x_{\pi (i)}x_{\pi (j)})`$
and observe that

$$\underset{\pi _1\text{S}_M}{}\frac{x_{\pi (1)}^{\lambda _11+N}\mathrm{}x_{\pi _1(M)}^{\lambda _MM+N}}{\mathrm{}(x_{\pi _1(1)},\mathrm{},x_{\pi _1(M)})}=(x_1\mathrm{}x_M)^{NM}s_\lambda (x_1,\mathrm{},x_M)$$
(23)
and

$$\underset{\pi _2\text{S}_{NM}}{}\frac{x_{\pi _2(M+1)}^{M1+N}\mathrm{}x_{\pi _2(N)}^0}{\mathrm{}(x_{\pi _2(M+1)},\mathrm{},x_{\pi _2(N)})}=1$$
(24)
After the permutations of type $`\pi _1`$ and $`\pi _2`$ have been performed in equation (14) we remain with the following sum over permutations $`\sigma \text{S}_N/\text{S}_M\times \text{S}_{NM}`$:

$`s_\lambda (x_1,\mathrm{},x_N)={\displaystyle \underset{\sigma \text{S}_N/\text{S}_M\times \text{S}_{NM}}{}}{\displaystyle \frac{\left(x_{\sigma (1)}^{NM}\mathrm{}x_{\sigma (M)}^{NM}\right)s_\lambda (x_{\sigma (1)}\mathrm{}x_{\sigma (M)})}{\underset{i[1,\mathrm{},M],j[M+1,\mathrm{},N]}{}(x_{\sigma (i)}x_{\sigma (j)})}}=`$ (25)
$`{\displaystyle \underset{\sigma \text{S}_N/\text{S}_M\times \text{S}_{NM}}{}}{\displaystyle \frac{\left(\underset{j=1}{\overset{M}{}}x_{\sigma (j)}^{NM}\right)s_\lambda (x_{\sigma (1)}\mathrm{}x_{\sigma (M)})\mathrm{}(x_{\sigma (1)}\mathrm{}x_{\sigma (M)})\mathrm{}(x_{\sigma (M+1)}\mathrm{}x_{\sigma (N)})}{\mathrm{}(x_{\sigma (1)}\mathrm{}x_{\sigma (N)})}}`$
We insert expression (25) to the sum in equation (12) and apply the Cauchy-Littlewood formula. It gives

$`{\displaystyle \frac{1}{\underset{j=1}{\overset{M}{}}\underset{i=1}{\overset{N}{}}\left(1x_iy_j\right)}}=`$ (26)
$`{\displaystyle \underset{\sigma \text{S}_N/\text{S}_{NM}\times \text{S}_M}{}}\left({\displaystyle \frac{\underset{i}{\overset{M}{}}x_{\sigma (i)}^{NM}}{\underset{i,j=1}{\overset{M}{}}\left(1x_{\sigma (i)}y_j\right)}}\right){\displaystyle \frac{\mathrm{}(x_{\sigma (1)},\mathrm{},x_{\sigma (M)})\mathrm{}(x_{\sigma (M+1)},\mathrm{},x_{\sigma (N)})}{\mathrm{}(x_{\sigma (1)},\mathrm{},x_{\sigma (N)})}}`$
It is easy to see that the formula above is equivalent to the formula (11).

## 3 Multi-variable integral representation for Cauchy transforms of orthogonal polynomials

In this section we compute $`𝒵_N^1[ϵ,H]_H`$, $`\text{dim}H=N`$, and show that this average is equal to $`\gamma _{N1}h_{N1}(ϵ)`$. Here $`\gamma _{N1}`$ is given by equation (10) and $`h_{N1}(ϵ)`$ is the Cauchy transform of the monic polynomial $`\pi _{N1}(x)`$ (equation (9)). The aim is just to illustrate our approach to the computation of the general correlation function $`𝒦_N(\widehat{ϵ},\widehat{\mu })`$ of characteristic polynomials on the simple example $`L=0,M=1`$.

Once $`M=1`$ algebraic identity (11) takes the form:

$$\frac{ϵ^{NM}}{𝒵_N[ϵ,H]}=\underset{\sigma \text{S}_N/\text{S}_{N1}}{}\left(\frac{x_{\sigma (1)}^{N1}}{ϵx_{\sigma (1)}}\right)\frac{\mathrm{}(x_{\sigma (2)},\mathrm{},x_{\sigma (N)})}{\mathrm{}(x_{\sigma (1)},\mathrm{},x_{\sigma (N)})}$$
(27)
Averaging the expression above with the probability measure $`𝒫^{(N)}(x_1,\mathrm{},x_N)`$ (equation (3)) we notice that each term yields the same contribution ($`𝒫^{(N)}(x_1,\mathrm{},x_N)`$ is symmetric with respect to the permutations). Thus we have

$`𝒵_N^1[ϵ,H]_H=`$
$`{\displaystyle \frac{ϵ^{MN}N}{Z_N}}{\displaystyle \frac{x_1^{N1}}{ϵx_1}\mathrm{}(x_2,\mathrm{},x_N)\mathrm{}(x_1,\mathrm{},x_N)e^{N\underset{i=1}{\overset{N}{}}V(x_i)}𝑑x_1\mathrm{}𝑑x_N}`$ (28)
The integral above can be further rewritten as

$$Z_{N1}𝑑x_1\frac{x_1^{N1}e^{NV(x_1)}}{ϵx_1}𝒵_{N1}[x_1,\stackrel{~}{H}]_{\stackrel{~}{H}},\text{dim}\stackrel{~}{H}=N1$$

Using the formula (4) we then obtain

$$𝒵_N^1[ϵ,H]_H=Nϵ^{NM}\frac{Z_{N1}}{Z_N}\frac{x^{N1}\pi _{N1}(x)e^{NV(x)}dx}{ϵx}$$
(29)
Now we note that

$$\frac{x^{N1}}{ϵx}=\frac{x^{N1}ϵ^{N1}}{ϵx}+\frac{ϵ^{N1}}{ϵx}$$

The first term above is a polynomial of degree $`N2`$ which is orthogonal to the polynomial$`\pi _{N1}(x)`$. Thus only the second term contributes to the integral. It gives

$$𝒵_N^1[ϵ,H]_H=\gamma _{N1}h_{N1}(ϵ)$$
(30)
or

$$\gamma _{N1}h_{N1}(ϵ)=\frac{1}{Z_N}\underset{j=1}{\overset{N}{}}(ϵx_j)^1e^{N\underset{j=1}{\overset{N}{}}V(x_j)}\mathrm{}^2(x_1,\mathrm{},x_N)dx_1\mathrm{}x_N$$
(31)
where we have used formulas (9), (10) and have expressed $`Z_N`$ as $`N!\underset{j=0}{\overset{N1}{}}c_j^2`$.

We expect that the representation of Cauchy transforms of orthogonal polynomials as multi-variable integrals (31) might be known, but we failed to trace it in the standard monographs on the subject.

## 4 Derivation of the formula for $`𝒦_N(\widehat{ϵ},\widehat{\mu })`$

Let us now consider the general case. The correlation function $`𝒦_N(\widehat{ϵ},\widehat{\mu })`$ defined by equation (6) is an integral over $`N`$ variables $`x_1,\mathrm{},x_N`$ (which are eigenvalues of $`H`$ ) with an integrand symmetric under permutations. It then follows that each component of the sum in equation (11) gives the same contribution to the correlation function $`𝒦_N(\widehat{ϵ},\widehat{\mu })`$. Thus we have

$`𝒦_N(\widehat{ϵ},\widehat{\mu })={\displaystyle \frac{N!}{(NM)!M!}}{\displaystyle \underset{j=1}{\overset{M}{}}}ϵ_j^{MN}{\displaystyle 𝑑x_1\mathrm{}𝑑x_N𝒫^{(N)}(x_1,\mathrm{},x_N)}`$
$`\times \left({\displaystyle \underset{i,j=1}{\overset{M}{}}}{\displaystyle \frac{x_i^{NM}}{ϵ_jx_i}}\right){\displaystyle \frac{\underset{j=1}{\overset{L}{}}𝒵_N[\mu _j,H]}{\underset{i=1}{\overset{M}{}}\underset{j=M+1}{\overset{N}{}}(x_ix_j)}}`$ (32)
We decompose the eigenvalue probability density function as

$`𝒫^{(N)}(x_1,\mathrm{},x_N)={\displaystyle \frac{Z_{NM}Z_M}{Z_N}}{\displaystyle \underset{i=1}{\overset{M}{}}}{\displaystyle \underset{j=M+1}{\overset{N}{}}}(x_ix_j)^2`$ (33)
$`\times 𝒫^{(M)}(x_1,\mathrm{},x_M)𝒫^{(NM)}(x_{M+1},\mathrm{},x_N)`$
With the decomposition above the integral expression for the correlation function $`𝒦_N(\widehat{ϵ},\widehat{\mu })`$ can be rewritten as

$`𝒦_N(\widehat{ϵ},\widehat{\mu })={\displaystyle \frac{N!}{(NM)!M!}}\left({\displaystyle \underset{j=1}{\overset{M}{}}}ϵ_j^{MN}\right){\displaystyle \frac{Z_{NM}Z_M}{Z_N}}`$ (34)
$`{\displaystyle 𝑑x_1\mathrm{}𝑑x_M𝒫^{(M)}(x_1,\mathrm{},x_M)F_{\widehat{ϵ},\widehat{\mu }}(x_1,\mathrm{},x_M)I^{(NM)}(x_1,\mathrm{},x_M,\widehat{\mu })}`$
where we have introduced two functions of $`M`$ variables $`x_1,\mathrm{},x_M`$. The first function is

$$F_{\widehat{ϵ},\widehat{\mu }}(x_1,\mathrm{},x_M)=\underset{i,j=1}{\overset{M}{}}\underset{k=1}{\overset{L}{}}\frac{x_i^{NM}(\mu _kx_i)}{ϵ_jx_i}$$
(35)
The second function, $`I^{(NM)}(x_1,\mathrm{},x_M,\widehat{\mu })`$, can be understood as an averaged value of the products of characteristic polynomials over an ensemble of $`NM`$ dimensional Hermitian matrices $`\stackrel{~}{H}`$ with eigenvalues $`x_{M+1},\mathrm{},x_N`$,

$`I^{(NM)}(x_1,\mathrm{},x_M,\widehat{\mu })=`$ (36)
$`{\displaystyle \underset{i=1}{\overset{M}{}}}{\displaystyle \underset{k=1}{\overset{L}{}}}𝒵_{NM}[x_i,\stackrel{~}{H}]𝒵_{NM}[\mu _k,\stackrel{~}{H}]_{\stackrel{~}{H}}`$
It is convenient to introduce $`L+M`$-dimensional vector $`\widehat{r}`$ whose components are the integration variables $`x_1,\mathrm{},x_M`$ and the elements of the vector $`\widehat{\mu }`$, i.e.

$$\widehat{r}=(x_1,\mathrm{},x_M,\mu _1,\mathrm{},\mu _L)$$
(37)
We exploit formula (5) which for $`I^{(NM)}(x_1,\mathrm{},x_M,\widehat{\mu })I^{(NM)}(\widehat{r})`$ gives the following expression

$$I^{(NM)}(\widehat{r})=\frac{1}{\mathrm{}(\widehat{r})}\text{det}\left[\pi _{NM+j1}(r_i)\right]_{1i,jL+M}$$
(38)
The Vandermonde determinant $`\mathrm{}(\widehat{r})`$ can be factorized as

$$\mathrm{}(\widehat{r})=\left\{\underset{k=1}{\overset{L}{}}\underset{j=1}{\overset{M}{}}(x_j\mu _k)\right\}\mathrm{}(x_1,\mathrm{},x_M)\mathrm{}(\mu _1,\mathrm{},\mu _k)$$
(39)
Now we insert equations (38) and (39) to the formula for the correlation function $`𝒦_N(\widehat{ϵ},\widehat{\mu })`$ (equation (34)). It gives:

$$𝒦_N(\widehat{ϵ},\widehat{\mu })=\frac{N!}{(NM)!M!}\left(\underset{j=1}{\overset{M}{}}ϵ_j^{MN}\right)\frac{Z_{NM}}{Z_N}\frac{1}{\mathrm{}(\widehat{\lambda },\widehat{\mu })}J_M(\widehat{ϵ},\widehat{\mu })$$
(40)
where $`J_M(\widehat{ϵ},\widehat{\mu })`$ is the $`M`$\- fold integral

$`J_M(\widehat{ϵ},\widehat{\mu })=`$ (41)
$`{\displaystyle 𝑑x_1\mathrm{}𝑑x_m\text{det}\left[\mathrm{\Phi }_i(x_j)\right]_{1i,jM}\text{det}\left[\pi _{NM+j1}(r_i)\right]_{1i,jL+M}}`$
where we have introduced $`M`$ functions $`\mathrm{\Phi }_i(x)`$

$$\mathrm{\Phi }_i(x)=\frac{x^{i1}G(x)}{\underset{j=1}{\overset{M}{}}(ϵ_jx)},G(x)=e^{NV(x)}x^{NM}$$
(42)
Thus the computation of the correlation function $`𝒦_N(\widehat{ϵ},\widehat{\mu })`$ is reduced to that of the integral of $`M`$ variables. This integral can be further transformed to a determinant form. To proceed we first note that $`\text{det}\left[\mathrm{\Phi }_i(x_j)\right]`$ can be simplified<sup>3</sup><sup>3</sup>3 This follows from the algebraic relation
$$\frac{x_j^{i1}G(x_j)}{\underset{l=1}{\overset{M}{}}(ϵ_lx_j)}=()^M\underset{K=1}{\overset{M}{}}\frac{ϵ_k^{i1}}{\underset{lk}{}(ϵ_kϵ_l)}\frac{G(x_j)}{x_jϵ_k}$$

which can be derived from identity (11). With the equation above we have
$$\text{det}\left[\mathrm{\Phi }_i(x_j)\right]=()^{M^2}\text{det}\left(\frac{G(x_i)}{ϵ_jx_i}\right)\text{det}\left(ϵ_j^{i1}\underset{lj}{}\frac{1}{ϵ_jϵ_l}\right)$$
Noting that
$$()^{M^2}\text{det}\left(ϵ_j^{i1}\underset{lj}{}\frac{1}{ϵ_jϵ_l}\right)=\frac{1}{\mathrm{}(\widehat{ϵ})}$$
we obtain equation (43)

$$\text{det}\left[\mathrm{\Phi }_i(x_j)\right]=\frac{1}{\mathrm{}(\widehat{ϵ})}\text{det}\left(\frac{G(x_i)}{ϵ_jx_i}\right).$$
(43)
Insert the expression (43) to the integral (41), rewrite the determinants as sums over permutations and find in this way that

$$J_M(\widehat{ϵ},\widehat{\mu })=\frac{M!}{\mathrm{}(\widehat{ϵ})\mathrm{}(\widehat{\mu })}\text{det}\left|\begin{array}{ccc}\stackrel{~}{\pi }_{NM}(ϵ_1)& \mathrm{}& \stackrel{~}{\pi }_{N+L1}(ϵ_1)\\ \mathrm{}& & \\ \stackrel{~}{\pi }_{NM}(ϵ_M)& \mathrm{}& \stackrel{~}{\pi }_{N+L1}(ϵ_M)\\ \pi _{NM}(\mu _1)& \mathrm{}& \pi _{N+L1}(\mu _1)\\ \mathrm{}& & \\ \pi _{NM}(\mu _L)& \mathrm{}& \pi _{N+L1}(\mu _L)\end{array}\right|$$
(44)
where

$$\stackrel{~}{\pi }_k(ϵ)=\frac{e^{NV(x)}x^{NM}}{ϵx}\pi _k(x)𝑑x,k[NM,NM+L1].$$
(45)
Using the orthogonality of monic polynomials $`\pi _k(x)`$ with respect to the measure $`e^{NV(x)}dx`$ we observe the relation between $`\stackrel{~}{\pi }_k(ϵ)`$ and the Cauchy transforms $`h_k(ϵ)`$:

$$\stackrel{~}{\pi }_k(ϵ)=2\pi iϵ^{NM}h_k(ϵ)$$
(46)
We insert expressions (44), (46) to formula (40) for the correlation function $`𝒦_N(\widehat{ϵ},\widehat{\mu })`$. Using the relations between coefficients $`Z_j`$, $`c_j`$, $`\gamma _j`$ we finally prove (8).

## 5 Correlation function $`𝒦_N(\widehat{ϵ},\widehat{\mu })`$ and the Riemann-Hilbert problem for orthogonal polynomials

Our formula (8) enables us to express the correlation function $`𝒦_N(\widehat{ϵ},\widehat{\mu })`$ in terms of solutions of Riemann-Hilbert problems for orthogonal polynomials proposed by Fokas, Its and Kitaev . It follows from (8) that the correlation function $`𝒦_N(\widehat{ϵ},\widehat{\mu })`$ is determined by monic orthogonal polynomials and their Cauchy transforms. In turn, the monic orthogonal polynomials and their Cauchy transforms can be understood as elements of a (matrix valued) solution of the following Riemann-Hilbert problem. Let contour $`\mathrm{\Sigma }`$ be the real line oriented from left to right. The upper side of the complex plane with respect to the contour will be called the positive one and the down side- the negative one. Once an integer $`n0`$ is fixed the Riemann-Hilbert problem is to find a $`2\times 2`$ matrix valued function $`Y=Y^{(n)}(z)`$ satisfying the following conditions:

* $`Y^{(n)}(z)\text{analytic}\text{in}\text{C}\mathrm{\Sigma }`$
* $`Y_+^{(n)}(z)=Y_{}^{(n)}(z)\left(\begin{array}{cc}1& e^{nV(z)}\\ 0& 1\end{array}\right),z\mathrm{\Sigma }`$
* $`Y^{(n)}(z)\left(I+𝒪(z^1)\right)\left(\begin{array}{cc}z^n& 0\\ 0& z^n\end{array}\right)\text{as}z\mathrm{}`$

Here $`Y_\pm ^{(n)}(z)`$ denotes the limit of $`Y^{(n)}(z^{})`$ as $`z^{}z\mathrm{\Sigma }`$ from the positive/negative side. As is proved by Fokas, Its and Kitaev  the solution of the Riemann-Hilbert problem is unique and is given by

$$Y^{(n)}(z)=\left(\begin{array}{cc}\pi _n(z)& h_n(z)\\ \gamma _{n1}\pi _{n1}(z)& \gamma _{n1}h_{n1}(z)\end{array}\right),\text{Im}z0$$
(47)
On comparing the formulae (8) and (47) we observe that the correlation function $`𝒦_N(\widehat{ϵ},\widehat{\mu })`$ can be expressed in terms of the matrix elements of the solution of the above Riemann-Hilbert problem. The relation provides us with a possibility to investigate the asymptotics of $`𝒦_N(\widehat{ϵ},\widehat{\mu })`$ at large $`N`$ for essentially any potential function $`V(x)`$ entering the probability distribution (3). The details will be presented in a forthcoming publication . Here we just outline the main steps. First we will show that the formula (8) can be rewritten as a determinant whose entries are two kernel functions. Those kernel functions are expressible in terms of the solution of the Riemann-Hilbert problem for orthogonal polynomials. The large $`N`$ asymptotics in the Dyson scaling limit can be studied by the steepest descent/stationary phase method for Riemann-Hilbert problems introduced by Deift and Zhou  and developed further in . As a result we will prove the universality of various quantities related to $`𝒦_N(\widehat{ϵ},\widehat{\mu })`$. In particular we shall be able to prove the universality of the distributions of the level curvatures and the local density of states, as well as to derive the Poisson kernel distribution of S-matrix in a random matrix model of quantum chaotic scattering. Understanding of universality of objects arising in random matrix theory in various scaling limits was recently a subject of intensive work both in physical  as well as mathematical  communities.

## 6 ”Duality relations” for the Gaussian case

For the particular case of the Gaussian measure $`V(x)=x^2/2`$ equation (5) for the ensemble average of products of characteristic polynomials can be written in an equivalent form:

$$_{N\times N}d[\widehat{H}]e^{\frac{N}{2}\text{Tr}\widehat{H}^2}det\left(\widehat{\lambda }\mathrm{𝟏}_N\mathrm{𝟏}_L\widehat{H}\right)_{L\times L}[d\widehat{Q}]e^{\frac{N}{2}\text{Tr}\widehat{Q}^2}det\left(\widehat{\lambda }i\widehat{Q}\right)^N$$
(48)
Here the integration in the left-hand side goes over the manifold of Hermitian matrices of the size $`N\times N`$, whereas in the right-hand side it goes over $`L\times L`$ Hermitian matrices.

This integral identity can be viewed as a certain ”duality” relation for matrix integrals. It emerged in various physical contexts, most notably in the context of matrix models of the string theory (see e.g. , p. 27) where it played an important role in understanding equivalence between one-matrix models of the quantum gravity and the so-called Kontsevich-type models . To understand equation (48) one recalls that the orthogonal polynomials for the Gaussian case are Hermite polynomials possessing an integral representation:

$$\pi _k(\lambda )e^{N\lambda ^2/2}_{\mathrm{}}^{\mathrm{}}𝑑qq^ke^{N\left(\frac{q^2}{2}i\lambda q\right)}$$
(49)
Substituting (49) into (5) one easily brings the right-hand side of the latter to the form

$`{\displaystyle \underset{j=1}{\overset{L}{}}}𝒵_N[\lambda _j,H]_H{\displaystyle \frac{e^{\frac{N}{2}\text{Tr}\widehat{\lambda }^2}}{\mathrm{\Delta }\{\lambda \}}}{\displaystyle 𝑑\widehat{Q}_\lambda \mathrm{\Delta }\{\widehat{Q}_\lambda \}e^{N\left[\frac{1}{2}\text{Tr}\widehat{Q}_\lambda ^2i\text{Tr}\widehat{Q}_\lambda \widehat{\lambda }\right]}det\widehat{Q}_\lambda ^N}`$ (50)
where $`\widehat{Q}_\lambda =\text{diag}(q_1,\mathrm{},q_L)`$. From the other hand the expression above can be obtained after shifting $`\lambda i\widehat{Q}i\widehat{Q}`$ in equation (48), diagonalizing the $`L\times L`$ matrix $`\widehat{Q}`$ and integrating out the angular degrees of freedom with help of the Itzykson-Zuber-Harish-Chandra formula. A particular case of the formula (50) was mentioned in  in the context of symmetric Jack polynomials.

It is easy to check that the Cauchy transforms of the Hermite polynomials are given by

$$h_k(ϵ)𝑑qq^ke^{N\left(\frac{q^2}{2}i\text{sgn}(\text{Im}ϵ)ϵq\right)}$$
(51)
Here the integration domain is $`0<q<\mathrm{}`$ for $`\text{Im}ϵ>0`$ and $`\mathrm{}<q<0`$ for $`\text{Im}ϵ<0`$. This gives us a possibility to rewrite our main object - correlation function (8) - in the form of $`M+L`$-fold integral :

$`𝒦_N(\widehat{ϵ},\widehat{\mu }){\displaystyle \frac{e^{\frac{N}{2}\text{Tr}\widehat{\mu }^2}}{\mathrm{\Delta }\{\mu \}\mathrm{\Delta }\{ϵ\}}}{\displaystyle 𝑑\widehat{Q}\mathrm{\Delta }\{\widehat{Q}\}e^{N\left[\frac{1}{2}\text{Tr}\widehat{Q}^2i\text{Tr}\widehat{Q}\widehat{E}\right]}det\widehat{Q}^{NM}}`$ (52)
where $`\widehat{Q}=(q_1,\mathrm{},q_L,q_{L+1},\mathrm{},q_{M+L})`$ and $`\widehat{E}=(\widehat{\mu },\widehat{ϵ}_+,\widehat{ϵ}_{})`$, with $`ϵ_+`$ and $`ϵ_{}`$ denoting spectral parameters $`ϵ_k`$ with positive (negative) imaginary part, respectively. The integration domain is $`\mathrm{}<q_l<\mathrm{}`$ for $`1lL`$ but $`0<q_l<\mathrm{}`$ or $`\mathrm{}<q_l<0`$ for $`L<lM+L`$, depending on the sign of the imaginary part of the corresponding spectral parameter.

Expression (52) is a generalization of duality relation (50). It arises most naturally in the method based on Gaussian integral representations and Itzykson-Zuber type integrations . Similar identities hold for the chiral ensemble (Laguerre polynomials), see  for more details.

## 7 Conclusions

In this paper we have found an exact formula for the general correlation function $`𝒦_N(\widehat{ϵ},\widehat{\mu })`$ containing both products and ratios of characteristic polynomials. Our result is valid for an arbitrary ensemble of Hermitian matrices of $`\beta =2`$ class. The obtained formula establishes a correspondence between the correlation function $`𝒦_N(\widehat{ϵ},\widehat{\mu })`$ and the Riemann-Hilbert problem for orthogonal polynomials. It is remarkable as it enables us to study the large $`N`$ asymptotics of this correlation function via the Riemann-Hilbert approach. Among interesting prospects for future research we would like to mention a challenging problem of extending our calculations to other symmetry classes $`\beta =1,4`$ (cf. ) as well as to ensembles of non-Hermitian random matrices important for the problems of quantum chaotic scattering .

## 8 Acknowledgements

We would like to thank A. Its, J. Keating, B. A. Khoruzhenko for valuable comments and discussions. This research was supported by EPSRC grant GR/13838/01 ”Random matrices close to unitarity or Hermitian”.