# Multicomponent holographic dark energy model with generalized entropy

Facundo Abaca1, Dario Javier Zamora2,3,
1 Departamento de Física, Facultad de Ciencias Exactas y Tecnología,
Universidad Nacional de Tucumán
Av. Independencia 1800, Tucuman, CP 4000, Argentina
2 Dipartimento di Scienza e Alta Tecnologia, Università degli Studi
dell’Insubria, Via Valleggio 11, 22100 Como, Italy
3 Instituto de Fisica del Noroeste Argentino, CONICET and Universidad Nacional
de Tucumán  E-mail<EMAIL_ADDRESS>
###### Abstract

Holographic dark energy cosmology, also known as entropic cosmology, provides
a concrete physical understanding of the late accelerated expansion of the
universe. The acceleration appears to be a consequence of entropy associated
with information storage in the universe. Therefore, the assumption of an ad-
hoc dark energy is not necessary. In this paper, we investigate the
implications of a multicomponent model (radiation and non-relativistic matter)
that includes a subdominant power-law term within a thermodynamically
admissible model. We use a generic power-law entropy and the temperature of
the universe horizon results from the requirement that the Legendre structure
of thermodynamics is preserved. We analyse the behaviour for different
combinations of the parameters and compare them with other cosmological
models, the observed redshift dependencies of the Hubble parameter $H$ and the
luminosity distance data obtained from supernovae. This is an early attempt to
analyse a multicomponent holographic dark energy model. Furthermore, the
analysis is based on a entropy scaling with an arbitrary power of the Hubble
radius instead of a specific entropy. This allows us to simultaneously infer
different models, compare them and conserve the scaling exponent as a
parameter that can be fitted with the observational data, thus providing
information about the form of the actual cosmological entropy and temperature.
We show that the introduced correction term is able to explain different
acceleration and deceleration periods in the late-time universe by solving the
model numerically. We discuss the advantages and disadvantages of holographic
dark energy models compared to mainstream cosmology.

##  1 Introduction

In mainstream cosmology, matter and spacetime emerged from a singularity and
evolved in four distinct periods, namely early inflation, radiation, dark
matter and late-time expansion. During the phases dominated by radiation and
dark matter, the universe is decelerating while the early and late-time
expansion are accelerating stages. The Lambda Cold Dark Matter ($\Lambda$CDM)
model assumes a cosmological constant $\Lambda$ and the existence of a
mysterious component called dark energy. This model is the simplest that can
explain an accelerated expansion of the late universe, which according to the
$\Lambda$CDM model is driven by dark energy. However, it implies some
theoretical peculiarities, such as the cosmic coincidence and the problem of
the cosmological constant [1, 2]. A possible connection between the
acceleration periods is still unknown, and interestingly, the most popular
candidate for dark energy driving the current acceleration phase ($\Lambda$\-
vacuum) relies on the cosmological constant and coincidence puzzle. To cope
with these difficulties, several alternative models have been proposed, see
for example [3, 4, 5, 6].

In this context, the holographic dark energy (HDE), also called entropic
cosmology, is an interesting attempt that can address this problem using the
holographic hypothesis. Horizon entropy is the backbone of HDE models, and
therefore any change in the horizon entropy affects the HDE model. Various
scenarios can be found, some of which are consistent with various astronomical
observations [7, 8, 9, 10, 11]. From this point of view, the controversial
dark energy component is not necessary. Here, the accelerated expansion in the
late period is based on the concept of entropic- force. Instead of dark
energy, we have the holographic principle and entropy as the source of the
late acceleration phase of the universe. An entropic force is an emergent
phenomenon that results from the natural tendency of a thermodynamic system to
increase its entropy, rather than from any particular underlying fundamental
force. There is no field associated with an entropic force. The force equation
is expressed in terms of a spatial dependence of the entropy $S$. The
cosmological entropic force $F$, is then given by

$F=-T\frac{dS}{dr_{H}},$ (1)

where $r_{H}$ is the Hubble radius.

Thermodynamic properties of the universe have always attracted attention [12,
13, 14] and in recent years especially entropic cosmology, see [15, 16, 17,
18] and citations therein. Here we follow the attempts of Easson, Frampton and
Smoot (EFS) [19, 20]. The first EFS model [19] assumes that the entropy and
temperature associated with the horizon of the universe are the Bekenstein-
Hawking entropy [21] and the Hawking temperature [22], respectively. Since
gravity is a long-range interaction, one can also use generalised statistical
mechanics to study gravitational systems [5, 23, 24, 25, 26, 27, 28, 29, 30,
31, 32, 33, 34, 35, 36, 37, 38]. This entropy was proposed in [38] in the
context of black holes. We would like to remind the reader that for such
systems the additive Bekenstein-Hawking entropy is proportional to the area of
the horizon, i.e. it is subextensive, while the non-additive
$S_{\delta=3/2}$\- entropy is proportional to the volume (at least in the case
of equal probabilities), i.e. it is extensive, as required by thermodynamics.
In the models [19, 20, 35], the expression for the temperature of the Hubble
horizon is not derived from a clean physical principle. It is simply assumed
to be the Hawking temperature expressed by the parameters of the universe,
namely

$T_{BH}(t)=\frac{\hbar c}{2\pi k_{B}r_{H}(t)}=\frac{\hbar H(t)}{2\pi k_{B}},$
(2)

where $c$ is the speed of light, $\hbar$ is the reduced Planck constant,
$k_{B}$ the Boltzmann constant, and $H(t)$ is the Hubble parameter. $H$ is
defined as

$H\equiv\frac{c}{r_{H}}=\frac{\dot{a}}{a},$ (3)

$a=a(t)$ being the scale factor.

Arbitrary combinations of entropy and temperature could violate the Legendre
structure of thermodynamics. This is not the case with the Bekenstein-Hawking
entropy and the Hawking temperature, as proposed by EFS [19]. This problem was
recently discussed in [39], where a physical principle for deriving the
thermodynamically consistent temperature associated with each class of entropy
was proposed. More precisely, the temperature must be derived from the
Legendre structure of thermodynamics

$\begin{split}G(V,T,p,\mu,...)&=U(V,T,p,\mu,...)-TS(V,T,p,\mu,...)\\\ &+pV-\mu
N(V,T,p,\mu,...)-...,\end{split}$ (4)

where $T,p,\mu$ are the temperature, pressure, and chemical potential, and
$U,S,V,N$ are the internal energy, entropy, volume, and the number of
particles of the system, respectively. This implies that, as detailed in [38],
in a Schwarzschild (3+1)-dimensional black hole, the relation

$\theta=1-d$ (5)

must hold, where $d$ is the dimension ($S\propto L^{d}$), and $\theta$ is the
corresponding exponent for the scaling of the temperature ($T\propto
L^{\theta}$), $L$ being a characteristic linear dimension of the d-dimensional
system. For the Bekenstein-Hawking entropy ($d=2$), the temperature depends on
$L^{-1}$ as the Hawking temperature.

The entropic force term, Eq. (1), affects the background evolution of the late
universe; in the present paper we do not focus on the inflation of the early
universe. It has been shown that models with entropic forces containing only
$H^{2}$ terms are not able to describe both decelerating and accelerating
phases on a single basis [39, 40, 41, 42]. Nevertheless, changes in the
deceleration parameter can be smoothly introduced by including subdominant
terms in the entropy of the horizon [20, 43]. Subdominant entropic terms are
not rare in nature, even if they are neglected in first approximations.
Moreover, they have a physical meaning and their inclusion is physically
motivated. With this additional term, the model predicts different periods of
acceleration and deceleration. With this in mind, we investigate here the
effects of adding a subdominant power-law term in the entropy and
simultaneously considering an multicomponent model with radiation and non-
relativistic matter.

Most of the models mentioned so far are single-component, and to be considered
a worthy alternative, HDE models must include different types of energy
densities in order to compare fairly with the mainstream model of cosmology,
i.e. the $\Lambda$CDM. Only a few papers have addressed the study of the
interaction of holographic dark energy with both radiation and non-
relativistic matter. The vast majority consider only matter, see for example
[44, 45, 46, 47, 48]. The present work is an early attempt to investigate a
multicomponent holographic dark energy model.

As an alternative to the hypothesis of the existence of dark energy, it is
important to compare HDE models with mainstream cosmology and observational
data. The aim of this work is to compare a multicomponent HDE model with the
$\Lambda$CDM. Although we recognize the need for observational validation,
this is beyond the scope of this work and is the subject of current research.

In this paper, we study how HDE models accommodate a viable cosmology for late
times without the consideration of dark energy by using an entropic force term
derived from a generic power-law entropy. Furthermore, we consider a
correction term necessary to obtain stages of acceleration and deceleration,
and derive the temperature of the horizon from the Legendre structure, making
our model compatible with thermodynamics. More details on all these
considerations can be found in the publications [39, 43]. We then go one step
further by considering a model in which both matter and radiation are present,
and we derive the Friedmann, acceleration and continuity equations. We also
present a different approach to the deduction of the Friedmann, continuity,
and acceleration equations. In previous works [39, 41, 43] the entropic
pressure is introduced in the continuity and acceleration equations, and from
there the Friedmann is deduced. Here, we deduce the acceleration equations
from the continuity and Friedmann ones. We show that both process takes to the
same result. Finally, we present numerical solutions to the equations and
discuss the results.

##  2 Multicomponent model

We consider a model under the assumption of a homogeneous, isotropic, and
spatially flat universe. Let us consider an entropy that scales with an
arbitrary positive power $d\in\mathbb{R^{+}}$ plus an additional term that
depends on a smaller power $0<\Delta<d$.

$\frac{S}{k_{B}}=A\left(\frac{r_{H}}{L_{P}}\right)^{d}+E\frac{\left(r_{H}/L_{P}\right)^{\Delta}-1}{\Delta},$
(6)

where $A$ and $E$ are dimensionless constants and $L_{P}=\sqrt{\hbar G/c^{3}}$
is the Planck length. By expressing the subdominant term in this way, we
obtain the logarithmic correction presented in [20] if we take
$\Delta\rightarrow 0$, since
$\lim_{\Delta\rightarrow 0}\frac{x^{\Delta}-1}{\Delta}=\ln{x}$. If we take
$E=0$ we also obtain the Bekenstein-Hawking entropy (d=2) and the $\delta=3/2$
entropy [38] (d=3) as special cases.

In order to make our model compatible with thermodynamics, according to Eq.
(5), the temperature must scale like $T\propto r_{H}^{1-d}$. Consequently, we
use

$T=\frac{T_{P}}{B}\left(\frac{r_{H}}{L_{P}}\right)^{1-d},$ (7)

as the temperature of the Hubble horizon, where $B$ is a dimensionless factor,
and $T_{P}=\sqrt{\hbar c^{5}/Gk_{B}^{2}}$ is the Planck temperature. The
entropic force is then given by

$\begin{split}F&\equiv-T\frac{dS}{dr_{H}}=-k_{B}\frac{d\,A}{B}.\frac{T_{P}}{L_{P}}\left[1+\frac{E}{d\,A}\left(\frac{r_{H}}{L_{P}}\right)^{\Delta-d}\right]\\\
&\equiv-C_{d}F_{P}(1+D_{d,\Delta}H^{d-\Delta}),\end{split}$ (8)

where $F_{P}\equiv k_{B}T_{P}/L_{P}=c^{4}/G$ is the Planck force, $C_{d}\equiv
d\,A/B$, and $D_{d,\Delta}\equiv E(L_{p}/c)^{d-\Delta}/(d\,A)$. Therefore, the
entropic pressure in the Hubble surface is

$p_{F}\equiv\frac{F}{4\pi r_{H}^{2}}=-\frac{C_{d}c^{2}}{4\pi
G}H^{2}(1+D_{d,\Delta}H^{d-\Delta}).$ (9)

Note that, when $d=2$ with $\Delta\rightarrow 0$, we obtain the $H^{4}$
correction term in [20], and when $D_{d,\Delta}=0$ we obtain the $H^{2}$-type
models [19, 39].

In our model, we consider a multicomponent universe with non-relativistic
matter, for which the energy density has the following dependence:

$\rho_{m}=\frac{3H_{0}^{2}\Omega_{m}}{8\pi
G}\left(\frac{a}{a_{0}}\right)^{-3},$ (10)

and radiation, with density:

$\rho_{r}=\frac{3H_{0}^{2}\Omega_{r}}{8\pi
G}\left(\frac{a}{a_{0}}\right)^{-4}.$ (11)

Remember that the key to the HDE models is that there is no dark energy. Hence
the Friedmann equation

$\left(\frac{H}{H_{0}}\right)^{2}=\frac{\Omega_{m}}{a^{3}}+\frac{\Omega_{r}}{a^{4}},$
(12)

where the constants $\Omega_{m}$ and $\Omega_{r}$ represent the fraction of
each type of component and they hold $\Omega_{r}+\Omega_{m}=1$, since the
assumption of the existence of dark energy is not necessary. Assuming an ideal
gas, the total pressure is

$p=\sum_{\omega}p_{\omega}=c^{2}\sum_{\omega}\omega\rho_{\omega},$ (13)

with $\omega$ being the parameter of the equation of state. $\omega$ takes the
values $0$ and $1/3$ for non-relativistic matter and radiation, respectively.

To obtain the continuity equation modified by the entropic force, we replace
the effective pressure $p^{\prime}=p_{m}+p_{r}+p_{F}$ in the usual continuity
equation

$\dot{\rho}+3\frac{\dot{a}}{a}\left(\rho_{m}+\rho_{r}+\frac{p^{\prime}}{c^{2}}\right)=0,$
(14)

thus arriving to

$\dot{\rho_{r}}+\dot{\rho_{m}}+3\frac{\dot{a}}{a}\left(\rho_{r}+\rho_{m}+\frac{p_{r}}{c^{2}}+\frac{p_{m}}{c^{2}}\right)=\frac{2C_{d}}{4\pi
G}H^{3}\left(1+D_{d,\Delta}H^{d-\Delta}\right).$ (15)

From Eq. (12) and (15), we can deduce the acceleration equation. To do so, we
multiply Eq. (12) by $a^{2}$ and take the time derivative, leading to:

$2\ddot{a}\dot{a}=\frac{8\pi
G}{3}(\dot{\rho}_{m}a^{2}+2\rho_{m}a\dot{a}+\dot{\rho}_{r}a^{2}+2\rho_{r}a\dot{a}).$
(16)

Dividing the last equation by $2\dot{a}a$, and replacing Eq. (15) we arrive to
the modified acceleration equation:

$\frac{\ddot{a}}{a}=-\frac{4\pi
G}{3}\left(\rho_{m}+\rho_{r}+\frac{3p_{m}}{c^{2}}+\frac{3p_{r}}{c^{2}}\right)+C_{d}H^{2}(1+D_{d,\Delta}H^{d-\Delta})$
(17)

Notice that in the case of a single component universe, following this
deduction, we would obtain the equations presented in [43], obtained in that
moment by a different approach. This shows consistency in the model. From the
three main eqs. (17), (15), and (12), we obtain

$\dot{H}=C_{d}D_{d,\Delta}H^{d-\Delta+2}+\frac{H^{2}}{2}\left[2C_{d}-3(\omega_{m}+1)\right]-\frac{3}{2}\frac{H_{0}^{2}\Omega_{r}}{a^{4}}(\omega_{r}-\omega_{m}),$
(18)

and whit a simple change of variables, we can obtain the differential equation

$\frac{dH}{da}=C_{d}D_{d,\Delta}\frac{H^{d-\Delta+1}}{a}+\frac{H}{2a}\left[2C_{d}-3(\omega_{m}+1)\right]-\frac{3}{2}\frac{H_{0}^{2}\Omega_{r}}{Ha^{5}}(\omega_{r}-\omega_{m}).$
(19)

Eq. (19) has a known form. First, we rewrite the equation by multiplying $2H$
and changing the variable $H^{2}=Z$:

$\frac{dZ}{da}=\frac{2C_{d}D_{d,\Delta}}{a}Z^{\frac{d-\Delta+2}{2}}+\frac{\left[2C_{d}-3(\omega_{m}+1)\right]}{a}Z-3\frac{H_{0}^{2}\Omega_{r}}{a^{5}}(\omega_{r}-\omega_{m}).$
(20)

This is a special case of the so-called Chini equation

$\frac{dy}{dt}=f(t)y(t)^{n}-g(t)y(t)+h(t),$ (21)

which generalizes the Riccati and the Abel equations and is, in general, not
solvable but some of its special cases are, see e.g. the book [49]. It is
known that there is a straightforward recipe for solving the equation in the
case that the Chini-invariant is independent of t.

$C=n^{-n}f(t)^{-1-n}h(t)^{1-2n}\left(f(t)h^{\prime}(t)-f^{\prime}(t)h(t)-nf(t)g(t)h(t)\right)^{n}.$
(22)

In the case examined here, however, this invariant depends on t, so that the
method in question is not applicable. We are therefore dealing with an
unsolvable equation. There are several ways of dealing with this equation,
e.g. working with approximations or solving the equation for different stages
(matter-dominated, radiation-dominated and entropic stages). In this paper, we
limit ourselves to the analysis of numerical methods to study the viability of
the model. More in-depth analyzes are the subject of current research and are
currently under development.

##  3 Results

The first thing we want to check is whether the model can accommodate stages
of acceleration and deceleration. And whether, in this case, they are similar
to those predicted by the $\Lambda$CDM model. To this end, we will study the
dependence of the scale factor $a$, with time. We obtain the corresponding
differential equation by changing the variables in equation (19):

$\left(\frac{da}{dt}a\right)^{2}\left[1+\frac{1}{2}[2C_{d}-3(1+\omega_{m})]\right]=\frac{3}{2}\Omega_{r}(\omega_{r}-\omega_{m})-\frac{C_{d}D_{d,\Delta}H_{0}^{d-\Delta}}{a^{d-\Delta-2}}\left(\frac{da}{dt}\right)^{d-\Delta+2}.$
(23)

We solve this equation numerically with the function NDsolve from Wolfram
Mathematica for different values of the parameters $C_{d}$, $D_{d,\Delta}$,
and $d-\Delta$. We set the other constants to
$H_{0}=67.4\frac{\frac{km}{s}}{mpc}$ (based on the Planck 2018 results [50]),
$\omega_{m}=0$, and $\omega_{r}=\frac{1}{3}$. Note that the relationship to
the underlying expansion parameters measured (e.g. via the angular diameter
distance) depends on the assumed cosmology [50]. As stated in page 17 of [50],
after presenting the value of $\Omega_{m}$: ”It is important to emphasize that
the values given in Eq. (13) assume the base-$\Lambda$CDM cosmology with
minimal neutrino mass. These estimates are highly modeldependent…”. Therefore,
the obtained values $\Omega_{m}=0.3...$ and $\Omega_{r}=0.000084$ might not be
valid for other models like the present one. In this work, due to the non-
existence of dark energy, the relationship $\Omega_{r}+\Omega_{m}=1$ must
apply. We have decided to set $\Omega_{r}=2.8\times 10^{-4}$, to maintain the
same ratio $\Omega_{r}/\Omega_{m}$ as in the $\Lambda$CDM model. Of course,
other considerations can be made to obtain values for $\Omega_{r}$ and
$\Omega_{m}$, for example by calculating them from the observational data, but
this is beyond the scope of this work. Since $d-\Delta$ is a difference of
indices, we expect to find whole numbers. Therefore, we set the parameter
$d-\Delta=\pm 1,\pm 2$. We find the numerical solution for different values of
$C_{d}$ and $D_{d,\Delta}$. We found 12 numerical solutions for different
combinations of values of the parameters in total. With all these
considerations in mind, some results of our model (in blue) together with the
$\Lambda$CDM solution are plotted in Fig. (1) for comparison.

(a)

(b)

Figure 1: Scale factor, $a$, as a function of time. Comparison of the
numerical solution of the eq. (23) (in blue) with the solution of the
$\Lambda$CDM model (red dashed). (a) $d-\Delta=-2$, $C_{d}=7$ and
$D_{d,\Delta}=-0.7$. (b) $d-\Delta=2$, $C_{d}=-1.5$ and $D_{d,\Delta}=2$

Several combinations of values of the parameters often gives similar curves,
so we show here only two characteristic behaviours. The first thing to notice
is that the numerical solution of the model shows different periods of
acceleration and deceleration, as expected since we included the correction
term. The similarity of the solution of our model with the $\Lambda$CDM model
is remarkable, especially in the late time, i.e. in the phase where the ’dark
energy’ is dominant. This suggests that HDE models are a serious alternative
to dark energy models to explain the accelerated expansion of the universe in
late time. The behaviour of our model in early time does not match that of the
$\Lambda$CDM model, but we can obtain qualitatively similar behaviours
depending on the values of the parameters, see for example Fig. (1.a). When
the parameter $d-\Delta$ takes on positive values, as in Fig. (1.b), the
change in concavity is not so evident, leading to a much older universe or to
behaviour similar to that described in [43]. Negative values of $d-\Delta$
seem more consistent with the $\Lambda$CDM model. However, this is an odd
combination of values as $d$ is the dimension and $\Delta$ is the scaling
exponent of the correction term. This means that the model needs a correction
term that is ’larger” than the main term in order to work. We will discuss
this with more detail in the next section.

In addition to the previous results, here we present an early simple
comparison between available observational data and the model. The aim of the
present paper is not to make a observational validation of the model, but we
think is important to present at least a simplified analysis. To do this, we
make a change to the variables in equation (19) and arrive at :

$\frac{dH}{dz}=\frac{3}{2}\frac{H_{0}^{2}}{H}(1+z)^{3}\Omega_{r}(\omega_{r}-\omega_{m})-\frac{C_{d}D_{d,\Delta}}{1+z}H^{d-\Delta+1}-\frac{[2C_{d}-3(1+\omega_{m})]}{2}\frac{H}{1+z}.$
(24)

In figures (2(a)), (2(b)) and (3) we have plotted the Hubble parameter $H$ as
a function of redshift $z$ using the data points from Table 1 in [51], which
in turn come from Supernova SN1a databases. In Fig. (2(a)) we assume
$d-\Delta=-2$ and determine the constants $C_{d}$ and $D_{d,\Delta}$ by
performing a sweep of the values. We show the plots of the combinations of the
values $d-\Delta=-2$, $C_{d}=-0.5$ and $D_{d,\Delta}=15800$ (red line) and
$D_{d,\Delta}=16800$ (blue line). In Fig. (2(b)) we show solutions with the
combination $d-\Delta=1$, $C_{d}=1$, and $D_{d,\Delta}=-0.004$ (red line),
$D_{d,\Delta}=-0.005$ (blue line).

(a)

(b)

Figure 2: Hubble parameter $H$ versus redshift $z$. The open circle with bars
are data points taken from Table 1 in [51]. In all cases, the value of $H_{0}$
is set to be $67.4km/s/Mpc$ based on the Planck 2018 results [50]. Red and
blue curves correspond to different combinations of values of the parameters
of the model. (a) $d-\Delta=-2$, $C_{d}=-0.5$ and $D_{d,\Delta}=15800$ (red
line), $D_{d,\Delta}=16800$ (blue line). (b) $d-\Delta=1$, $C_{d}=1$, and
$D_{d,\Delta}=-0.004$ (red line), $D_{d,\Delta}=-0.005$ (blue line).

Finally, we show the combination $d-\Delta=2$, $C_{d}=1$ and
$D_{d,\Delta}=-0.00003$ (red), $D_{d,\Delta}=-0.00004$ (blue) in Fig. (3). We
see that the model has a satisfactory fit to the data points in all cases.

Figure 3: Hubble parameter $H$ versus redshift $z$. The open circle with bars
are data points taken from Table 1 in [51]. In all cases, the value of $H_{0}$
is set to be $67.4km/s/Mpc$ based on the Planck 2018 results [50]. Red and
blue curves correspond to different combinations of values of the parameters
of the model: $d-\Delta=2$, $C_{d}=1$ and $D_{d,\Delta}=-0.00003$ (red line),
and $D_{d,\Delta}=-0.00004$ (blue line).

##  4 Discussion and conclusions

In this paper, we have explored a multicomponent holographic dark energy
model, also known as entropic cosmology, as an alternative explanation for the
late accelerated expansion of the universe. The key idea behind this model is
to attribute the acceleration to entropy associated with information storage
in the universe, eliminating the need for an ad-hoc dark energy component. The
model includes a subdominant power-law term, introducing a correction to the
main entropy term. A noteworthy observation is that for the model to predict
similar behaviour to $\Lambda$ CDM model, the parameter $d-\Delta$ has to be
negative. This implies that the correction term must be larger than the main
term, emphasizing the importance of the correction in driving the dynamics.
The necessity of negative values for $d-\Delta$ requires careful consideration
of the physical implications. Maybe the cosmological entropy of the Hubble
horizon is not a known entropy yet, such as Bekenstein-Hawking, Boltzmann-
Gibbs, or Tsallis entropy, but a new different one, with a similar form to the
correction term here presented. A second problem is that several parameters in
the model lack clear physical interpretation, raising concerns about the
model’s theoretical foundation. All the new parameters in this model come from
the entropy considered. The parameter $d-\Delta$ has a clear interpretation,
but $C_{d}$ and $D_{d,\Delta}$ no. They are a combination of the exponents
$d$, $\Delta$, and of the coefficients $A$, $E$, and $B$ escorting the entropy
and temperature. In the simplest case in which we take $A=B=E=1$, $C_{d}$ and
$D_{d,\Delta}$ depend only in physical constant, acquiring in this way clear
physical interpretation. However, in this case, the model loses all its
predictive power, since the parameters are necessary to fit the available
observational data. The model’s predictions need to undergo observational
verification, acknowledging that this is an essential aspect yet to be
explored in detail. Comparison with available observational data, particularly
Hubble parameter and luminosity distance from supernovae, is necessary to
establish the model’s credibility. Observational validation is crucial, and
the model’s compatibility with a broader range of observational data remains
to be investigated.

The holographic dark energy model aims to address the shortcomings of the
$\Lambda$CDM model, which relies on a cosmological constant and introduces the
mysterious dark energy. The model can potentially explain different stages of
acceleration and deceleration without invoking dark energy and produces
results comparable to the $\Lambda$CDM model in late-time expansion. This is
an advantage since the model provides a unique and concrete physical
interpretation of the late accelerated expansion through entropy and the
holographic principle. The correction term allows for flexibility in
describing different periods of acceleration and deceleration.

In conclusion, while the holographic dark energy model presents intriguing
possibilities and potential advantages over the traditional $\Lambda$CDM
model, careful scrutiny, theoretical refinement, and observational testing are
essential for establishing its validity and addressing the identified
concerns. Further research and detailed observational comparisons will be
crucial in determining the true merits and limitations of this cosmological
approach.

## Acknowledgment

This work received financial support from CONICET, and CIN (Argentinian
agencies).

## References

  * [1] S. Weinberg, The Cosmological Constant Problem, Rev. Mod. Phys. 61, 1 (1989).
  * [2] S. M. Carroll, The Cosmological Constant, Living Reviews in Relativity 4. (2001)
  * [3] S. Weinberg, Cosmology (Oxford University Press, 2008).
  * [4] R. S. Ellis, R. Maartens, and M. A. H. MacCallum, Relativistic Cosmology (2012).
  * [5] N. Komatsu and S. Kimura, Evolution of the Universe in Entropic Cosmologies via Different Formulations, Phys. Rev. D 89, 123501 (2014).
  * [6] J. Solà, Cosmological Constant and Vacuum Energy: Old and New Ideas, in Journal of Physics: Conference Series, 453 (Institute of Physics Publishing, 2013), p. 012015.
  * [7] X. Zhang, F.Q. Wu, The Interacting Generalized Ricci Dark Energy Model in Non-Flat Universe, Phys. Rev. D 72, 043524 (2005)
  * [8] X. Zhang, F.Q. Wu, Constraints on holographic dark energy from the latest supernovae, galaxy clustering, and cosmic microwave background anisotropy observations, Phys. Rev. D 76, 023502 (2007)
  * [9] Q.G. Huang, Y.G. Gong, Supernova Constraints on a holographic dark energy model, JCAP 0408, 006 (2004)
  * [10] K. Enqvist, S. Hannestad, M.S. Sloth, Searching for a holographic connection between dark energy and the low-l CMB multipoles, JCAP 0502, 004 (2005)
  * [11] J.Y. Shen, B. Wang, E. Abdalla, R.K. Su, Constraints on the dark energy from the holographic connection to the small l CMB suppression, Phys. Lett. B 609, 200 (2005)
  * [12] P. C. W. Davies, Cosmological horizons and the generalised second law of thermodynamics, Class. Quantum Grav. 4 L225 (1987).
  * [13] I. Prigogine, J. Geheniau, E. Gunzig, and P. Nardone, Thermodynamics of cosmological matter creation, PNAS 85 (20) 7428-7432 (1988).
  * [14] I. Prigogine, Thermodynamics and cosmology. Int J Theor Phys 28, 927–933 (1989).
  * [15] Y. Cai, J. Liu, and H. Li, Entropic cosmology: A unified model of inflation and late-time acceleration, Physics Letters B 690, pp. 213-219 (2010).
  * [16] K. Karami, A. Sheykhi, M. Jamil, F. Felegary and M. M. Soltanzadeh, Thermodynamical description of interacting entropy-corrected new agegraphic dark energy, EPL 93, 69001 (2011).
  * [17] S. Mitra, S. Saha, and S. Chakraborty, Modified Hawking temperature and entropic force: A prescription in FRW model. Modern Physics Letters A 30, 1550058 (2015).
  * [18] S. Viaggiu, First law of thermodynamics for dynamical apparent horizons and the entropy of Friedmann universes, General Relativity and Gravitation 47, 1-16 (2015).
  * [19] D. A. Easson, P. H. Frampton, and G. F. Smoot, Entropic Accelerating Universe, Phys. Lett. B 696, 273 (2011).
  * [20] D. A. Easson, P. H. Frampton, and G. F. Smoot, Entropic Inflation, Int. J. Mod. Phys. A 27,1250066 (2012).
  * [21] J. D. Bekenstein, Black Holes and Entropy, Phys. Rev. D 7, 2333 (1973).
  * [22] S. W. Hawking, Black Hole Explosions?, Nature 248, 30 (1974).
  * [23] A. Majhi, Non-extensive statistical mechanics and black hole entropy from quantum geometry, Phys. Lett. B 775, 32 (2017)
  * [24] S. Abe, General pseudoadditivity of composable entropy prescribed by the existence of equilibrium, Phys. Rev. E 63, 061105 (2001)
  * [25] H. Touchette, When is a quantity additive, and when is it extensive?, Phys. A 305, 84 (2002)
  * [26] A. Sayahian Jahromi et al., Generalized entropy formalism and a new holographic dark energy model, Phys. Lett. B 780, 21 (2018)
  * [27] H. Moradpour, S.A. Moosavi, I.P. Lobo, et al. Thermodynamic approach to holographic dark energy and the Rényi entropy. Eur. Phys. J. C 78, 829 (2018).
  * [28] Komatsu, Eur. Phys. J. C 77, 229 (2017)
  * [29] H. Moradpour, A. Bonilla, E.M.C. Abreu, J.A. Neto, Phys. Rev. D 96, 123504 (2017)
  * [30] H. Moradpour, Int. J. Theor. Phys. 55, 4176 (2016)
  * [31] E.M.C. Abreu, J. Ananias Neto, A.C.R. Mendes, W. Oliveira, Phys. A 392, 5154 (2013)
  * [32] E.M.C. Abreu, J. Ananias Neto. Phys. Lett. B 727, 524 (2013)
  * [33] E.M. Barboza Jr., R.C. Nunes, E.M.C. Abreu, J.A. Neto, Phys. A Stat. Mech. Appl. 436, 301 (2015)
  * [34] R.C. Nunes et al., JCAP 08, 051 (2016)
  * [35] N. Komatsu and S. Kimura, Entropic Cosmology for a Generalized Black-Hole Entropy, Phys. Rev. D 88, 083534 (2013).
  * [36] N. Komatsu, S. Kimura, Phys. Rev. D 90, 123516 (2014)
  * [37] N. Komatsu, S. Kimura, Phys. Rev. D 93, 043530 (2016)
  * [38] C. Tsallis and L. J. L. Cirto, Black Hole Thermodynamical Entropy, Eur. Phys. J. C 73, 1 (2013).
  * [39] D. J. Zamora and C. Tsallis, Thermodynamically consistent entropic-force cosmology, Phys. Lett. B 827, 136967 (2022).
  * [40] E. L. D. Perico, J. A. S. Lima, S. Basilakos, and J. Solà, Complete Cosmic History with a Dynamical $\Lambda=\Lambda(H)$ Term, Phys. Rev. D 88, 063531 (2013).
  * [41] N. Komatsu and S. Kimura, Non-Adiabatic-like Accelerated Expansion of the Late Universe in Entropic Cosmology, Phys. Rev. D 87, 043531 (2013).
  * [42] S. Basilakos, D. Polarski, and J. Solà, Generalizing the Running Vacuum Energy Model and Comparing with the Entropic-Force Models, Phys. Rev. D 86, 043010 (2012).
  * [43] D. J. Zamora and C. Tsallis, Thermodynamically consistent entropic late-time cosmological acceleration, Eur. Phys. J. C 82, 689 (2022).
  * [44] M.A. Zadeh, A. Sheykhi, H. Moradpour, et al. Note on Tsallis holographic dark energy. Eur. Phys. J. C 78, 940 (2018).
  * [45] S. Chakrabarti and S. Chattopadhyay. Viewing the Cosmological Consequences of Modified Holographic Dark Energy in Various Interaction Scenarios. Zeitschrift für Naturforschung A, 73, pp. 251-257, (2018).
  * [46] W. Zimdahl and D. Pavon, Interacting holographic dark energy, Classical and Quantum Gravity 24 5461, (2007).
  * [47] Z. Zhang et al, Revisit of the interaction between holographic dark energy and dark matter, Journal of Cosmology and Astroparticle Physics 06, 009, (2012).
  * [48] J. Zhang, X. Zhang, H. Liu, Statefinder diagnosis for the interacting model of holographic dark energy, Physics Letters B 659, pp.26-33, (2008).
  * [49] E. Kamke, Differentialgleichungen: Lösungen und Lösungsmethoden, Band I: Gewöhnliche Differentialgleichungen, Leipzig, (1951).
  * [50] Planck Collaboration, Planck 2018 Results. VI. Cosmological Parameters, Astron. Astrophys. 641, (2018).
  * [51] A. Pradhan and A. Dixit, Tsallis Holographic Dark Energy Model with Observational Constraints in the Higher Derivative Theory of Gravity, New Astron. 89, 101636 (2021).


Solitons and Normal Random Matrices

I. M. Loutsenko1, V. P. Spiridonov2,3 and O. V. Yermolayeva1

1 Laboratoire de Physique Mathématique, Centre de Recherches Mathématiques,

Université de Montréal, Montréal

2 Laboratory of Theoretical Physics, Joint Institute for Nuclear Research,
Dubna

3 National Research University “Higher School of Echonomics”, Moscow

###### Abstract

We discuss a general relation between the solitons and statistical mechanics
and show that the partition function of the normal random matrix model can be
obtained from the multi-soliton solutions of the two-dimensional Toda lattice
hierarchy in a special limit.

##  1 Introduction

Present paper is devoted to intersections between the theory of solitons and
statistical mechanics of two-dimensional “Coulomb-Dyson gases”, or normal
random matrices. The literature on the random matrix models and related
Coulomb gases is quite extensive due to the links to various important
problems of mathematics and physics. For an introduction to matrix models see,
e.g., [6, 20]. The normal matrix model was introduced in [3] in connection to
the quantum Hall effect and investigated further in detail in [4]. The
connection between the two-dimensional Toda lattice (2DTL) hierarchy and
normal random matrices has been established in [4], in [21] it emerged in the
context of the theory of Hele-Shaw flows, or Laplacian growth (for a review
see, e.g., [19]).

It is a well-known fact that there is deep connection between the theory of
various matrix ensembles and that of tau-function of integrable hierarchies,
see, e.g., [13, 14, 24] and references therein. Our “solitonic” models of
statistical mechanics (besides the normal matrix model discussed here, one can
mention a variety of ensembles considered in [16]) contain partition functions
of these ensembles as particular or limiting cases and therefore generalize
this connection: they introduce lattices and use non-trivial boundary
conditions to the Coulomb-gas models related to the matrix ensembles.

The Coulomg gas formalism is a powerful tool of treating two-dimensional
problems. Its in depth character is reflected in connecting seemingly
different physical notions like two-dimensional statistical mechanics and
quantum field theory models, random matrices and integrable systems [6]. It
naturally emerges in the investigation of critical phenomena in two-dimensiona
spin lattice systems [23] with further output to the two-dimensional conformal
field theory. In this sense, the interpretation of solitonic tau-functions as
lattice gas (Ising chains) partition functions described below is not an
occasional coincidence, but yet another indication on the universality of this
formalism with many physical manifestations.

We start from a brief introduction to basic concepts of statistical mechanics
of lattice gases and its connection to the theory of solitons established
previously in [15, 16, 17, 18] (see also, [31]). Statistical mechanics studies
macroscopic properties of systems with large number of degrees of freedom.
Such systems can exist in a discrete (but possibly infinite) set of
microstates. A microstate defines the values of all possible microscopic
variables.

Consider a lattice consisting of $N$ points (sites) on the plane which can be
occupied by some particles. In the lattice gas model no more than one particle
can occupy each site, i.e. the number of particles at the $i$-th site can be
$\nu_{i}=0$ or $\nu_{i}=1$. These filling numbers
$\nu_{1},\nu_{2},\dots,\nu_{N}$ constitute a set of microscopic variables.
Then, a microstate of the lattice gas corresponds to a binary string (i.e., a
sequence of zeros and ones) of the length $N$. The total number of microstates
of the gas equals $2^{N}$. Statistical mechanics considers a macroscopic
system, that is the $N\to\infty$ “thermodynamic” limit, and investigates the
corresponding behaviour of the various macroscopic variables, for example,
that of the total number of particles

$n=\sum_{i=1}^{N}\nu_{i}.$ (1.1)

This number $n$ can vary from 0, when the lattice is empty, to $N$ when the
lattice is fully occupied by particles.

The probability to find system in the microscopic state $\nu$ fixed by $N$
filling numbers $\nu=\\{\nu_{1},\nu_{2},\dots,\nu_{N}\\}$ equals

$p(\nu)=\frac{1}{Z}e^{-\beta(E(\nu)-\mu n)},$ (1.2)

where $E(\nu)=E(\nu_{1},\dots,\nu_{N})$ is the system energy and $n$ is
defined in (1.1). The parameter $\beta>0$ is called the inverse temperature,
and the parameter $\mu$ is the chemical potential. Since the sum of
probabilities of all possible micro-states of the system equals to 1, the
value $Z$ in the normalization factor in (1.2) equals to the following sum
over all possible system states:

$Z=\sum_{\nu\in{\cal S}}e^{-\beta(E(\nu)-\mu n)}.$ (1.3)

This sum is called the partition function of the grand canonical ensemble.
Here ${\cal S}$ stands for the set of all possible microscopic configurations
of the lattice gas given by $2^{N}$ binary strings $\nu$ of the length $N$.

From definitions (1.2) and (1.3), it follows that

$\frac{\partial\log Z}{\partial\beta}=-\langle E\rangle+\mu\langle
n\rangle,\quad\frac{\partial\log Z}{\partial\mu}=\beta\langle n\rangle,$

where

$\langle E\rangle=\sum_{\nu\in{\cal S}}p(\nu)E(\nu),\quad\langle
n\rangle=\sum_{\nu\in{\cal S}}p(\nu)n(\nu)$

are the average energy of the system and the average number of particles in
the system respectively. The quantities $\langle E\rangle$ and $\langle
n\rangle$ describe the results of measurements of the corresponding
macroscopic variables.

One can split the set of microstates of our system ${\cal S}$ into $N+1$
disjoint sets: ${\cal S}={\cal S}_{0}\cup{\cal S}_{1}\cup{\cal
S}_{2}\cup\dots\cup{\cal S}_{N}$, where ${\cal S}_{i}\cap{\cal
S}_{j}=\emptyset$ if $i\not=j$. Here, the set ${\cal S}_{0}$ corresponds to
the states with $n=0$ particles on the lattice (empty lattice), and the set
${\cal S}_{k}$ corresponds to the states with $n=k$ with ${\cal S}_{N}$ being
the state of fully occupied lattice ($n=N$, i.e. all $\nu_{i}=1$). Then from
(1.3) we have $Z=\sum_{n=0}^{N}{\cal Z}_{n}e^{\beta\mu n},$ where

${\cal Z}_{n}=\sum_{\nu\in{\cal S}_{n}}e^{-\beta E(\nu)}$ (1.4)

is the partition function for the system with fixed number of particles $n$.
The quantity (1.4) is called the partition function of the canonical ensemble,
or the $n$-particle partition function. For the canonical ensemble of $n$
particles we have

$\langle E\rangle=-\frac{\partial\log{\cal
Z}_{n}}{\partial\beta}=\frac{\partial\beta F}{\partial\beta},\quad
F=-\frac{1}{\beta}\log{\cal Z}_{n},$

where $F$ is called the “free energy” of the $n$-particle system. In the
thermodynamic limit $n\to\infty$, one usually is interested in the asymptotics
of free energy per particle $F/n$.

Let $\zeta_{i}$ be a complex coordinate of the $i$-th site of the lattice. We
consider a system where particles interact pairwise through the two-particle
potential $V(z,z^{\prime})=V(z^{\prime},z)$ as well as they interact with
external fields through a one-particle potential $w(z)$. In other words, the
energy of interaction of particles occupying the $i$-th and $j$-th sites
equals $V_{ij}=V(\zeta_{i},\zeta_{j})$. We consider only the two-particle
interaction with $V(z,z)=+\infty$, so that the condition that the filling
factor $\nu_{i}$ cannot be greater than 1 will be fulfilled automatically (see
below). The energy of interaction of $i$-th site’s particle with the external
fields is $w_{i}=w(\zeta_{i})$. The total energy of the gas then equals to

$E=\sum_{1\leq i<j\leq N}V_{ij}\nu_{i}\nu_{j}+\sum_{i=1}^{N}w_{i}\nu_{i}.$
(1.5)

This function, relations (1.1) and (1.3) define the grand partition function
of the gas

$Z=\sum_{\nu_{1}=0,1}\dots\sum_{\nu_{N}=0,1}e^{-\beta\left(\sum_{1\leq i<j\leq
N}V_{ij}\nu_{i}\nu_{j}+\sum_{i=1}^{N}(w_{i}-\mu)\nu_{i}\right)}.$ (1.6)

Note that the number $N$ stands here for the number of the lattice sites and
not for the number of particles $n$, which is not fixed in the grand canonical
ensemble (one should not confuse $Z$ with the $N$-particle partition function
${\cal Z}_{N}$).

For the $n$-particle partition function, defined by (1.4) and (1.5), we obtain

${\cal
Z}_{n}=\frac{1}{n!}\sum_{z_{1}\in\zeta}\dots\sum_{z_{n}\in\zeta}e^{-\beta(\sum_{1\leq
i<j\leq n}V(z_{i},z_{j})+\sum_{i=1}^{n}w(z_{i}))},$ (1.7)

where $\zeta$ stands for the set of lattice points
$\zeta=\\{\zeta_{1},\zeta_{2},\dots,\zeta_{N}\\}$. Note that, since
$V(z,z)=+\infty$, we do not have to care about the configurations in (1.7)
with filling numbers exceeding 1, because all summands with $z_{i}=z_{j}$
vanish there. The factor $1/n!$ appears in (1.7) because, for a given
configuration $\nu$ with the number of particles $n$, these $n$ particles are
identical in the sum (1.5), while it is not so in (1.7).

A crucial observation on the connection of lattice statistical mechanics
models with the solitons was done in [15, 16, 17]. Namely, soliton
$\tau$-functions of integrable hierarchies and partition functions of
particular Ising models or lattice gas models have identical structure. More
precisely, the $N$-soliton $\tau$-function of many integrable hierarchies can
be written in the following Hirota form [1, 10]

$\tau_{N}=\sum_{\nu_{1}=0,1}\dots\sum_{\nu_{N}=0,1}\exp\left(\sum_{1\leq
i<j\leq N}A_{ij}\nu_{i}\nu_{j}+\sum_{i=1}^{N}\phi_{i}\nu_{i}\right),$ (1.8)

where the sum is performed over $2^{N}$ configuration of $N$ discrete
variables $\nu_{i}$, each of $\nu_{i}$ taking values $0$ or $1$. In equation
(1.8), $A_{ij}$ stands for the phase shift acquired due to interaction of the
$i$-th and $j$-th solitons, while $\phi_{i}$ stands for the $i$-th soliton
phase. The phase shift $A_{ij}$ depends on the momenta of corresponding
solitons and the phase $\phi_{i}$ depends on the momentum of the $i$-th
soliton and the full set of the integrable hierarchy times.

In the expression for the $N$-soliton $\tau$-function (1.8) one can recognize
the grand partition function of a gas models (1.6) on the lattice consisting
of $N$ sites: the variable $\nu_{i}$ is the filling factor of the $i$-th
lattice site and $\phi_{i}$ is proportional to the sum of the chemical
potential and the external potential at the $i$-th site. In this picture, the
phase shift $A_{ij}$ is proportional to the two-body interaction potential
between particles occupying the $i$-th and the $j$-th sites:

$A_{ij}=-\beta V_{ij},\quad\phi_{i}=-\beta\left(w_{i}-\mu\right).$ (1.9)

An important point is that this correspondence between the $\tau$-functions
and the partition functions has a restriction — the parameter $\beta$ appears
to be fixed in the $\tau$-functions (see below).

Before proceeding to the main topic, we mention that some connections between
soliton solutions of the integrable hierarchies and certain matrix models have
been established in [16]. In particular, the soliton solutions of the KP
hierarchy corresponding to the circular lattice are related to the Gaudin
matrix models, which interpolate between the Dyson and uniform unitary
ensembles. Also, soliton solutions of the KP and BKP hierarchies corresponding
to one-dimensional exponential lattices are related to the one-dimensional
translationally-invariant Ising models with non-local interactions. The
simplest models of such type are described by the self-similar potentials of
the one-dimensional Schrödinger equation [28] and they are related to special
infinite-soliton solutions of the KdV hierarchy [15]. Despite of the fact that
the inverse temperature $\beta$ is fixed in these models, through playing with
the model parameters it is possible to go to the zero temperature limit
exhibiting a particular critical behavior of the corresponding Ising chains in
the thermodynamic limit [18].

In this paper we consider Coulomb gases on general two-dimensional lattices
and discuss their relation to normal random matrices. Namely, we derive the
partition function of a discretized $n\times n$ normal matrix model (or a
Coulomb-Dyson gas of $n$ charges on $N$ sites on the plane) from the
$N$-soliton $\tau$-function of the two-dimensional Toda lattice (2DTL)
hierarchy.

##  2 Coulomb-Dyson Gases

Coulomb gas is a statistical ensemble of charged particles on the plane
interacting through the Coulomb potential $V$. On the $(x,y)$-plane the
potential created by the particle, placed at the point
$(x^{\prime},y^{\prime})$, is proportional to the Green function of the 2D
Laplace operator

$\Delta
V(x,y;x^{\prime},y^{\prime})=-2\pi\delta(x-x^{\prime})\delta(y-y^{\prime}),\quad
V(x,y;x^{\prime},y^{\prime})=V(x^{\prime},y^{\prime};x,y),$ (2.1)

with certain boundary conditions imposed on $V$. For the plane without
boundaries

$V(z,z^{\prime})=-\log|z-z^{\prime}|=-\frac{1}{2}\left(\log(z-z^{\prime})+\log(\bar{z}-\bar{z}^{\prime})\right),$

where the complex variables notation $z=x+\text{i}y$, $\bar{z}=x-\text{i}y$ is
used. If the system has an ideal dielectric boundary $\Gamma$, the normal to
$\Gamma$ component of the gradient of $V$ vanishes on this boundary, i.e.
$(n_{x}\partial_{x}+n_{y}\partial_{y})V(x,y,x^{\prime},y^{\prime})=0$, where
$(n_{x},n_{y})$ is the normal to $\Gamma$ at the point $(x,y)\in\Gamma$. The
tangential component of the gradient vanishes at $\Gamma$, if it is an ideal
conductor boundary, i.e. at fixed $z^{\prime}$ the conductor boundary $\Gamma$
is an equipotential surface of $V$. A useful way of solving (2.1) with such
boundary conditions is provided by the method of images. In what follows, we
consider systems where every charge has a finite number of images created by
the boundaries.

The electrostatic energy of a system of $n$ particles with charges ${\cal
Q}_{i}$, $i=1,2,\dots,n$, is

$E_{n}=\sum_{1\leq i<j\leq n}{\cal Q}_{i}{\cal
Q}_{j}V(z_{i},z_{j})+\sum_{1\leq i\leq n}{\cal
Q}_{i}^{2}\tilde{V}(z_{i})+\sum_{1\leq i\leq n}{\cal Q}_{i}W(z_{i}),$ (2.2)

where $z_{i}=x_{i}+\text{i}y_{i}$ are the coordinates of the particles. The
first term in (2.2) is the energy of interaction between different charges.
The second term is the sum of self-energies: the self-energy of a charge is
the energy of interaction between the charge and its own images. The third
term describes an interaction of charges with external fields.

The Dyson gas is a one-component (i.e., all particles have equal charges)
Coulomb gas at certain fixed temperatures corresponding to different random
matrix models. Here, we consider the case $\beta{\cal Q}_{i}=2$. Then, without
loss of generality, we can set all ${\cal Q}_{i}$ in (2.2) equal to unity and
let the inverse temperature to be $\beta=2$. Next, we consider the lattice
version of the Dyson gas where particles occupy a set of points
$\zeta=\\{\zeta_{1},\dots,\zeta_{N}\\}$ on the $N$ sites lattice (i.e.,
$z_{i}\in\zeta$) and no more than one particle can occupy each site. Then, the
grand partition function of the system corresponding to energy (2.2) equals

$Z=\sum_{\nu_{1}=0,1}\dots\sum_{\nu_{N}=0,1}e^{-\beta\left(\frac{1}{2}\sum_{1\leq
i\not=j\leq
N}V(\zeta_{i},\zeta_{j})\nu_{i}\nu_{j}+\sum_{i=1}^{N}(w(\zeta_{i})-\mu)\nu_{i}\right)},$

where $w(z)=\tilde{V}(z)+W(z)$ and $\mu$ stands for the chemical potential.

##  3 Ensembles Related to KP, BKP and 2DTL Hierarchies

It turned out that the Coulomb interaction $V(\zeta_{i},\zeta_{j})$ can be
identified with the phase shifts of different integrable hierarchies [17]. In
this section we describe connections between the ensembles of the Coulomb-
Dyson gases and soliton solutions of the KP (including KdV hierarchy as a
subcase) and the BKP hierarchies established in [17] and of the 2DTL hierarchy
found in [31].

###  3.1 KP Hierarchy

The $N$-soliton $\tau$-function of the Kadomtsev-Petviashvili (KP) hierarchy
can be written in the Hirota form (1.8) [10] with

$A_{ij}=\log\frac{(a_{i}-a_{j})(b_{i}-b_{j})}{(a_{i}+b_{j})(b_{i}+a_{j})},$
(3.1)
$\phi_{i}=\varphi_{i}+\sum_{n=1}^{\infty}\left(a_{i}^{n}-(-b_{i})^{n}\right)t_{n},$

where $t_{n}$, $n=1,2,3,\dots$, is an infinite set of independent variables
called the hierarchy “times”, $(a_{i},b_{i})$ is the two-dimensional momentum
of the $i$-th soliton and $\varphi_{i}$ is its initial phase. The first non-
trivial equation of the hierarchy, the celebrated KP-equation, involves three
independent variables $x=t_{1}$, $y=t_{2}$ and $t=t_{3}$ and has the form

$3\frac{\partial^{2}u}{\partial y^{2}}=\frac{\partial}{\partial
x}\left(4\frac{\partial u}{\partial t}+6u\frac{\partial u}{\partial
x}-\frac{\partial^{3}u}{\partial x^{3}}\right),\quad
u=-2\partial^{2}_{x}\log\tau.$ (3.2)

All equations of the hierarchy can be encoded in the single bi-linear Hirota
residue equation for the $\tau$-function

$\oint_{z=0}\frac{dz}{2\pi\text{i}}e^{\xi(\boldsymbol{t}^{\prime}-\boldsymbol{t},z)}\tau(\boldsymbol{t}^{\prime}-[z^{-1}])\tau(\boldsymbol{t}+[z^{-1}])=0,$
(3.3)

where

$\boldsymbol{t}=\\{t_{1},t_{2},t_{3}\dots\\},\quad\xi(\boldsymbol{t},z)=t_{1}z+t_{2}z^{2}+t_{3}z^{3}+\dots,\quad[z^{-1}]=\\{z^{-1},\tfrac{1}{2}z^{-2},\tfrac{1}{3}z^{-3},\dots\\},$
(3.4)

and $\oint_{z=0}dz/(2\pi\text{i})$ (also denoted as $\,\text{res}\,_{z=0}$)
means taking residue at $z=0$, i.e. a coefficient of $z^{-1}$ in the Laurent
series in $z$. A full (infinite) set of the bilinear differential Hirota
equations can be obtained from equation (3.3) as the conditions of vanishing
of various coefficients of the Taylor series in
$\varepsilon_{i}=t^{\prime}_{i}-t_{i}$.

Substituting

$a_{i}=\zeta_{i},\quad
b_{i}=-\bar{\zeta}_{i},\quad\text{Im}\,\zeta_{i}>0,\quad i=1,2,\dots,N,$ (3.5)

in (3.1), we obtain

$V(z,z^{\prime})=-\log|z-z^{\prime}|+\log|\bar{z}-z^{\prime}|.$

This potential satisfies equation (2.1). For a charge placed at a point
$z^{\prime}$, the real line $\text{Im}\,z=0$ is an equipotential surface of
the potential $V(z,z^{\prime})$ created by this charge. Therefore,
$V(z,z^{\prime})$ is the Coulomb potential in the upper half-plane with an
ideal conducting boundary along the real line. Equivalently, we may say that
this is a potential created by a positive charge at the point $z^{\prime}$ in
the upper half-plane and its reflection image of the opposite charge located
in the lower half-plane at $\bar{z}^{\prime}$. Thus, the $N$-soliton solution
of the KP hierarchy with the momenta defined by (3.5) and the phases

$\phi_{i}=-\beta\left(w(\zeta_{i})-\mu\right)=-\beta\left(\log|\bar{\zeta}_{i}-\zeta_{i}|+W(\zeta_{i})-\mu\right),$
(3.6)

describes a Coulomb gas at the temperature $\beta=2$ on a lattice in the upper
half-plane with the ideally conductive boundary along the $x$-axis (see the
first picture on Figure 1).

Figure 1: From the left to the right. KP hierarchy: a two-dimensional Coulomb
gas above an ideal conductor. BKP hierarchy: a two-dimensional Coulomb gas in
the corner between an ideal dielectric (the horizontal axis) and an ideal
conductor (the vertical axis). 2DTL hierarchy: a two-dimensional Coulomb gas
in the exterior of the unit disc with ideally conducting boundary. Positive
charges are shown as white squares while their negative images are shown as
black squares. The interactions between different charges are shown by dashed
lines, while the interactions between charges and their own images are shown
by solid lines.

In (3.6), the potential $w(z)$ is a sum of the self-interaction potential,
corresponding to the “charge-image” interaction,

$\tilde{V}(z)=\log|\bar{z}-z|,$

and of the external potential $W(z)$. The latter is a sum of a confining
potential $U(z)$, keeping particles in the compact plane domain and
determining the initial phases of solitons

$\varphi_{i}=-\beta(U(\zeta_{i})-\mu),$

and of a harmonic function, determined by the KP “times”, that corresponds to
the electric field created by some distant external charges, i.e.,

$W(z)=U(z)-\frac{1}{2}\sum_{p=1}^{\infty}(z^{p}-\bar{z}^{p})t_{p}.$

We draw attention to the fact that in order for the external potential to be
real all hierarchy “times” must be purely imaginary. The real axis is the
equipotential surface of the above harmonic function.

Putting all Coulomb particles on the vertical axis we come to the condition
that $a_{i}=b_{i}$ for the KP soliton momenta, which corresponds to the KdV-
hierarchy solitons. This system has a natural interpretation as a nonlocal
Ising chain with exactly computable free energy per spin in the
translationally invariant cases [15, 18].

###  3.2 BKP Hierarchy

Soliton solutions of the BKP hietrarchy in the Hirota form (1.8) are detrmined
by $N$ two-dimensional momenta $(a_{i},b_{i})$, an infinite number of odd BKP
times $t_{2p-1}$, $p=1,2,3,\dots,$ and $N$ initial phases of the solitons
$\varphi_{i}$, so that [11]

$A_{ij}=\log\frac{(a_{i}-a_{j})(b_{i}-b_{j})(a_{i}-b_{j})(b_{i}-a_{j})}{(a_{i}+a_{j})(b_{i}+b_{j})(a_{i}+b_{j})(b_{i}+a_{j})},$
(3.7)
$\phi_{i}=\varphi_{i}+\sum_{p=1}^{\infty}(a_{i}^{2p-1}+b_{i}^{2p-1})t_{2p-1}.$
(3.8)

The first non-trivial equation in the hierarchy (the BKP equation) involves
three independent variables $t_{1}$, $t_{3}$ and $t_{5}$. Similar to the KP
hierarchy (3.3) and (3.4), the BKP hierarchy can be encoded in a single Hirota
residue equation:

$\oint_{z=0}\frac{dz}{2\pi\text{i}z}e^{\xi(\boldsymbol{t}^{\prime}-\boldsymbol{t},z)}\tau(\boldsymbol{t}^{\prime}-2[z^{-1}])\tau(\boldsymbol{t}+2[z^{-1}])=\tau(\boldsymbol{t})\tau(\boldsymbol{t^{\prime}}),$

with only odd times involved, i.e.
$\boldsymbol{t}=\\{t_{1},t_{3},t_{5},\dots$} and
$[z^{-1}]=\\{z^{-1},\tfrac{1}{3}z^{-1},\tfrac{1}{5}z^{-5},\dots\\}.$

Choosing soliton momenta as follows

$a_{i}=\zeta_{i},\quad b_{i}=\bar{\zeta}_{i},\quad\text{Re}\,\zeta_{i}>0,$

the BKP phase shifts (3.7) yield the following two-particle interaction
potential

$V(z,z^{\prime})=-\log|z-z^{\prime}|-\log|\bar{z}-z^{\prime}|+\log|z+z^{\prime}|+\log|\bar{z}+z^{\prime}|,\quad\beta=2,$

which is a Coulomb potential in a quarter of the plane. For convenience, let
us choose the right upper corner of the plane, $\text{Re}\,z>0$,
$\text{Im}\,z>0$. Then the horizontal boundary of the corner $y=0$ is the
ideal dielectric, while the vertical boundary $x=0$ is the ideal conductor.
Each charge has now three images: one of the same sign and two of the opposite
signs, see picture 2 on Figure 1.

The phases of solitons have the form

$\phi_{i}=-\beta\left(\tilde{V}(\zeta_{i})+W(\zeta_{i})-\mu\right),$

where the self-interaction potential is

$\tilde{V}(z)=\log|\bar{z}-z|-\log|\bar{z}+z|-\log|2z|,$

and the external potential $W$ is a sum of the confining potential $U$ and a
harmonic function determined by the BKP times (see equations (3.8) and (3.5))

$W(z)=U(z)-\frac{1}{2}\sum_{p=1}^{\infty}(z^{2p-1}+\bar{z}^{2p-1})t_{2p-1}.$

This harmonic function satisfies all the necessary conditions on the
boundaries of the corner.

After setting $a_{i}=b_{i}$ in the phase shifts of the BKP hierarchy we obtain
the KdV solitons phase shifts multiplied by 2. This corresponds to the same
lattice gas or Ising chain as in the plain KdV case, but for the twice bigger
value of the inverse temperature $\beta=4$ [17]. Again, this value of the
temperature corresponds to the particular random matrix ensembles.

###  3.3 2TDL Hierarchy

The two-dimentional Toda lattice (2DTL) hierarchy consists of difference-
differential equations [29]. It involves an infinite number of continuous
independent variables (2DTL times) $t_{p}$, where $p$ runs over all positive
and negative integers, excluding zero, as well as a discrete variable
$m\in\mathbb{Z}$. The simplest equation in the hierarchy is the 2DTL equation:

$\frac{\partial^{2}u(m)}{\partial t_{1}\partial
t_{-1}}=e^{-u(m+1)}+e^{-u(m-1)}-2e^{-u(m)},$ (3.9)

where $u(m)$ depend on $t_{p}$, $p=\pm 1,\pm 2,\pm 3,\dots$. In terms of the
$\tau$-function the 2DTL equation reads as

$\frac{\partial\tau(m)}{\partial t_{1}}\frac{\partial\tau(m)}{\partial
t_{-1}}-\tau(m)\frac{\partial^{2}\tau(m)}{\partial t_{1}\partial
t_{-1}}=\tau(m-1)\tau(m+1),\quad
u(m)=\log\frac{\tau(m)^{2}}{\tau(m+1)\tau(m-1)}.$ (3.10)

The $N$-soliton $\tau$ function has the form (1.8) with the phase shifts

$A_{ij}=\log\frac{(a_{i}-a_{j})(b_{i}-b_{j})}{(a_{i}-b_{j})(b_{i}-a_{j})},$
(3.11)

which are essentially the KP phase shifts (one should replace $b_{i}$ to
$-b_{i}$ in (3.1)), and the soliton phases

$\phi_{i}=\varphi_{i}+m(\log a_{i}-\log
b_{i})+\sum_{p=-\infty,p\not=0}^{\infty}\left(a_{i}^{p}-b_{i}^{p}\right)t_{p}.$

Choosing soliton momenta as follows

$a_{i}=\zeta_{i}/R,\quad b_{i}=R/\bar{\zeta}_{i},\quad|\zeta_{i}|>R,$

where $R$ is a positive real number, we come to the interaction potential

$V(z,z^{\prime})=-\log|z-z^{\prime}|+\log|\bar{z}z^{\prime}-R^{2}|-\log
R,\quad\beta=2.$ (3.12)

We set now $R$ equal to 1, $R=1$, and obtain

$V(z,z^{\prime})=-\log|z-z^{\prime}|+\log|\bar{z}z^{\prime}-1|.$ (3.13)

Then $V(z,z^{\prime})$ satisfies equation (2.1) in the exterior of the ideally
conducting unit disc. In other words, for a charge placed at $z^{\prime}$, the
disc boundary (i.e., the unit circle $|z|=1$) is the equipotential surface of
(3.13). Equivalently, one can say that this is a potential created by a
positive charge at the point $z^{\prime}$, in the exterior of the disc, and
its reflection (where reflection is the inversion with respect to the unit
circle) image of the opposite charge located inside the disc at
$1/\bar{z}^{\prime}$, see the third picture on Figure 1.

The self-interaction potential in the 2DTL case is

$\tilde{V}(z)=\log\left|z-\bar{z}^{-1}\right|,$

while the external potential is a sum of the confining potential $U$ and the
harmonic function determined by the 2DTL times and the discrete variable $m$:

$W(z)=U(z)-m\log|z|-\frac{1}{2}\sum_{p=1}^{\infty}\left(\left(z^{p}-\bar{z}^{-p}\right)t_{p}+\left(\bar{z}^{p}-z^{-p}\right)\bar{t}_{p}\right).$

Note, that we have redefined the negative times as $t_{-p}=-\bar{t}_{p}$,
$p>0$ in order for the potential to be real. Obviously, the unit circle is an
equipotential surface of the harmonic field.

Concluding this section we note that the Coulomb interaction potential and the
boundary can be conformally transformed. Indeed, let $f(z)$ be a conformal
mapping from the exterior of a simple compact domain $\Omega$ to the exterior
of the unit circle $|z|=1$, then
$\tilde{V}(z,z^{\prime})=V(f(z),f(z^{\prime}))$ is also a Coulomb potential
whoose equipotential surface is $\partial\Omega$. Thus $N$-soliton solution
with momenta $a_{i}=f(\zeta_{i})$ and $b_{i}=\bar{f}(\bar{\zeta}_{i})$
corresponds to the grand partition function of the lattice Coulomb gas in an
exterior of conducting domain $\Omega$.

The simplest mapping $f(z)=z/R$ transforms (3.13) to (3.12). As an example of
less trivial mapping one can take $f(z)=z/2-\sqrt{z^{2}/4-1}$ (i.e. the
mapping inverse to $z\to z+1/z$) from the exterior of the real segment
$[-2,2]$ to that of the unit disc. In this case, the ideal conductor is placed
along that segment.

##  4 2DTL Solitons and Normal Random Matrices

Let us now remove the constraint $R=1$, and consider a Coulomb gas in the
exterior of the disc of an arbitrary radius $R$. Note that the 2DTL hierarchy
is invariant under the transformation

$m\to m-j,\quad t_{p}\to R^{p}t_{p},\quad\bar{t}_{p}\to
R^{p}\bar{t}_{p},\quad\tau(m)\to R^{m^{2}+cm+C}\tau(m),$

where $c,C$ and $j\in\mathbb{Z}$ are arbitrary constants. In terms of the
dependent variable $u$ this transformation produces the shift $u\to u-2\log
R$, see (3.10). In particular, the $\tau$-function

$\tilde{\tau}_{N}(m,t_{1},\bar{t}_{1},t_{2},\bar{t}_{2},\dots)=R^{m^{2}}\tau_{N}(m-1,Rt_{1},R\bar{t}_{1},R^{2}t_{2},R^{2}\bar{t}_{2},\dots)$

is also a solution of the hierarchy. This follows from the bilinear equation
(notations are defined in (3.4))

$\oint_{z=0}\frac{dz}{2\pi\text{i}}z^{m^{\prime}-m}e^{\xi(\boldsymbol{t}^{\prime}-\boldsymbol{t},z)}\tau(m^{\prime},\boldsymbol{t}^{\prime}-[z^{-1}],\boldsymbol{\bar{t}}^{\prime})\tau(m,\boldsymbol{t}+[z^{-1}],\boldsymbol{\bar{t}})=$
$=\oint_{z=0}\frac{dz}{2\pi\text{i}}z^{m^{\prime}-m}e^{\xi(\boldsymbol{\bar{t}}^{\prime}-\boldsymbol{\bar{t}},z^{-1})}\tau(m^{\prime}+1,\boldsymbol{t}^{\prime},\boldsymbol{\bar{t}}^{\prime}-[z])\tau(m-1,\boldsymbol{t},\boldsymbol{\bar{t}}+[z])$

which encodes the whole hierarchy (for a review see, e.g., [29]).

Repeating the above computations for an arbitrary $R$, we get the transformed
$N$-soliton $\tau$-function

$\tilde{\tau}_{N}(m)=R^{m^{2}}\sum_{\nu_{1}=0,1}\dots\sum_{\nu_{N}=0,1}e^{-\beta
E(\nu_{1},\dots,\nu_{N})},\quad\beta=2,$ (4.1)

where

$E=-\sum_{1\leq i<j\leq
N}\left(\log|\zeta_{i}-\zeta_{j}|-\log|\zeta_{i}\bar{\zeta}_{j}-R^{2}|+\log
R\right)\nu_{i}\nu_{j}$
$-\sum_{i=1}^{N}\left(\frac{\varphi_{i}}{2}+(m-1)(\log|\zeta_{i}|-\log
R)+{\cal U}(\zeta_{i})\right)\nu_{i}$

and ${\cal U}$ is the harmonic function of the form

${\cal
U}(z)=-\frac{1}{2}\sum_{p=1}^{\infty}\left(\left(z^{p}-\frac{R^{2p}}{\bar{z}^{p}}\right)t_{p}+\left(\bar{z}^{p}-\frac{R^{2p}}{z^{p}}\right)\bar{t}_{p}\right).$

Fixing the initial phases of solitons as

$\varphi_{i}=-\beta U(\zeta_{i})-\log R,\quad\beta=2,$ (4.2)

where $U(z)$ is a confining potential, we then take the limit of small $R$. In
this limit, when the disc contracts to a point, the harmonic function becomes
analytic at $z=0$

${\cal
U}(z)=-\frac{1}{2}\sum_{p=1}^{\infty}\left(z^{p}t_{p}+\bar{z}^{p}\bar{t}_{p}\right),\quad
R\to 0,$ (4.3)

and the energy of the gas becomes

$E=-\sum_{1\leq i<j\leq
N}\nu_{i}\nu_{j}\log|\zeta_{i}-\zeta_{j}|-\sum_{i=1}^{N}\left(U(\zeta_{i})+{\cal
U}(\zeta_{i})+(m-n)\log|\zeta_{i}|\right)\nu_{i}-n\left(\frac{n}{2}-m\right)\log
R.$ (4.4)

Here

$n=\sum_{i=1}^{N}\nu_{i}$

is the number of particles. In the $R\to 0$ limit the last term in (4.4), i.e.
$-n\left(\frac{n}{2}-m\right)\log R$, is the dominant one. Combining it with
the factor $R^{m^{2}}$ in (4.1), we obtain the common factor of the
$n$-particle terms in the $\tilde{\tau}_{N}(m)$-function

$R^{m^{2}}e^{\beta n\left(\frac{n}{2}-m\right)\log R}=R^{(m-n)^{2}},$

since $\beta=2$. We see that only the $n$-particle terms with
$n=\sum_{i}\nu_{i}=m$ are finite for $R\to 0$, while other terms vanish. Also,
since the total number of gas particles $n$ cannot be negative or exeed the
number of lattice sites $N$, the lowest vanishing rate of
$\tilde{\tau}_{N}(m)$ for $m\leq 0$ is reached at $n=0$, while for $m>N$ it is
reached at $n=N$ (i.e., at one of the ends of the interval $0\leq n\leq N$),
so that all the terms in (4.1) for $m$ outside of this interval vanish. Thus,
in the limit $R\to 0$, we obtain

$\tilde{\tau}_{N}(m)={\cal Z}_{m},\quad 0\leq m\leq
N,\quad\tilde{\tau}_{N}(m<0)=\tilde{\tau}_{N}(m>N)=0,$

where ${\cal Z}_{m}$ is the partition function of the gas of $m$ Coulomb
particles on the lattice $\zeta=\\{\zeta_{1},\zeta_{2},\dots,\zeta_{N}\\}$ in
the complex plane without boundaries. No more than one particle can occupy a
lattice site and

${\cal
Z}_{m}=\frac{1}{m!}\sum_{z_{1}\in\zeta}\dots\sum_{z_{m}\in\zeta}e^{-\beta
E_{m}(z_{1},\dots,z_{m})},\quad\beta=2,$ (4.5)

with the gas energy being

$E_{m}(z_{1},\dots,z_{m})=-\sum_{1\leq i<j\leq m}\log|z_{i}-z_{j}|+\sum_{1\leq
i\leq m}\left(U(z_{i})+{\cal U}(z_{i})\right).$ (4.6)

Finally, we can write ${\cal Z}_{0}=1$ and

${\cal
Z}_{m}=\frac{1}{m!}\sum_{\ell_{1}=1}^{N}\dots\sum_{\ell_{m}=1}^{N}\prod_{1\leq
i<j\leq
m}\left|\zeta_{\ell_{i}}-\zeta_{\ell_{j}}\right|^{2}\prod_{i=1}^{m}\exp\left(-2U(\zeta_{\ell_{i}})-2{\cal
U}(\zeta_{\ell_{i}})\right),\quad m>0.$ (4.7)

Thus, we have shown that the $m$-particle partition function ${\cal Z}_{m}$,
determined in (4.5), (4.6) or (4.7), with the harmonic external potential
${\cal U}$ fixed in (4.3), is a $\tau$-function of the two-dimensional Toda
lattice of the length $N-1$. In other words, the function

$\tau(m,t_{1},\bar{t}_{1},t_{2},\bar{t}_{2},\dots)=\left\\{\begin{array}[]{c}{\cal
Z}_{m}(t_{1},\bar{t}_{1},t_{2},\bar{t}_{2},\dots),\quad{\rm if}\quad
m=0,1,2,\dots,N,\\\ 0,\quad{\rm otherwise,}\end{array}\right.$ (4.8)

is a solution of the 2DTL hierarchy. In terms of $u(m)$, the 2DTL equation
(3.9) becomes a system of $N-1$ differential equations for
$u(1),\dots,u(N-1)$. For $m=2,3,\ldots,N-2$, the equations are fixed in (3.9),
while for $m=1$ and $m=N-1$ we have
$\partial_{t_{1}}\partial_{t_{-1}}u(1)=e^{-u(2)}-2e^{-u(1)}$ and
$\partial_{t_{1}}\partial_{t_{-1}}u(N-1)=e^{-u(N-2)}-2e^{-u(N-1)}$,
respectively.

The $m$-particle partition function (4.7) (or, equivalently, (4.5), (4.6)) can
be rewritten in the form

${\cal Z}_{m}=\frac{1}{m!}\int\prod_{1\leq i<j\leq
m}|z_{i}-z_{j}|^{2}\prod_{j=1}^{m}e^{-2{\cal
U}(z_{j})}\rho(z_{j})dx_{j}dy_{j},$ (4.9)

where

$\rho(z)=\sum_{i=1}^{N}e^{-2U(z)}\delta(x-X_{i})\delta(y-Y_{i}),\quad\zeta_{i}=X_{i}+\text{i}Y_{i}.$

In the continuous $N\to\infty$ limit, the spacing between sites tends to zero,
i.e. $\zeta_{i}=\epsilon\xi_{i}$, $\epsilon\to 0$, while the area of the
lattice, which is of the order of $\epsilon^{2}N$, remains finite. In this
limit the measure $\rho$ tends to the continuous measure

$\rho(z)=\epsilon^{-2}e^{-2U(z)}\varrho(z),$ (4.10)

where $\varrho(z)$ is a normalized finite density of the lattice sites. We
recall that the $\tau$-function is defined modulo a gauge factor, and, in
particular, it can be multiplied by any constant. Therefore, the diverging
factor $\epsilon^{-2}$ can be discarded in the above equality.

The $\tau$-function (4.9) is the partition function of the random normal
$m\times m$ matrix model Depending on the measure, we obtain either continuous
or discrete version of the model. Note that the standard hermitian matrix
model is a special case of the normal matrix model, when the measure $\rho$ is
concentrated on the real line, i.e. this model enters our considerations as
well.

By choosing $\varphi_{i}=-\beta U(\zeta_{i})+(2\ell-1)\log R$,
$\ell\in\mathbb{Z}$, instead of (4.2), we will get a partition function of the
Coulomb gas in the presence of a fixed point charge of the value $\ell$ placed
at the point $z=0$. Determinantal representations of the partition functions
of gases with several fixed charges were given, e.g., in [4, 9].

Note that the partition function (4.9) of the normal matrix model with an
arbitrary (continuous or/and discrete) measure can be obtained without the
lattice scaling procedure, by applying our limiting procedure directly to the
well-known solution of the hierarchy (see e.g. [12, 31])

$\tau=\sum_{n\geq 0}\frac{1}{n!}\int\dots\int\prod_{1\leq i<j\leq
n}\frac{(p_{i}-p_{j})(q_{i}-q_{j})}{(p_{i}-q_{j})(q_{i}-p_{j})}\prod_{i=1}^{n}e^{m(\log
p_{i}-\log
q_{i})+\sum_{k=-\infty\atop\makebox[-5.0pt]{}k\not=0}^{\infty}\left(p_{i}^{k}-q_{i}^{k}\right)t_{k}}r(p_{i},q_{i})dp_{i}dq_{i},$

for which the $N$-soliton solution is a particular case corresponding to
$r(p,q)=\sum_{i=1}^{N}r_{i}\delta(p-a_{i})(q-b_{i})$. The procedure is
essentially the same as in the case of the $N$-soliton solutions, except the
continuum limit at the final stage is not needed.

##  5 Conclusions

The grand partition function of the Coulomb-Dyson gas and the $m$-particle
partition function related to the normal matrix model are usually obtained by
using different “group-like” elements in the free-fermion approach to
integrable hierarchies [31]. In this article, we have shown that the latter
can be obtained from the former one by a special limiting procedure.

It is well known that the partition function of the standard normal random
matrix model having the continuous measure emerges in the dispersionless limit
of the 2DTL hierarchy corresponding to a special $m\to\infty$ scaling limit.
In the dispersionless limit one introduces the scaled (“slow”) times $T_{i}$,
such that $T_{\pm i}=ht_{\pm i}$, and $T_{0}=hm$. Here, one considers the
double scaling limit, such that $h\to 0$, $m\to\infty$ with $T_{0}$ remaining
finite. The scaling parameter $h$ plays the role of spacing on the Toda
lattice (i.e., the minimal step in “time” $T_{0}$, not to confuse with spacing
on the soliton momentum lattice). When $h\to 0$, the time $T_{0}$ becomes the
continuous variable for slowly varying solutions. For instance, in this limit
the Toda lattice equation (3.9) becomes

$\frac{\partial^{2}u}{\partial
T_{1}\partial{T_{-1}}}=\frac{\partial^{2}e^{-u}}{\partial T_{0}^{2}}.$

Commutators in the Lax representation of the hierarchy degenerate to the
Poisson brackets and the Lax operators become functions, i.e. the
dispersionless limit is a kind of a quasiclassical limit with $h$ playing a
role of the “Plank constant” (for more details see e.g. [30]).

From the point of view of the matrix models, introduction of the “slow” times
corresponds to the scaling of the harmonic potential ${\cal U}(z,t)={\cal
U}(z,T)/h$ (see eq. (4.3)). Also by making the measure (4.10) dependent on the
parameter $h$ through the $h$-dependence of the confining potential
$U(z,h)=v(z)/h$, we get the matrix model with the scaled total potential
$U(z,h)+{\cal U}(z,t)=\frac{1}{h}\left(v(z)+{\cal U}(z,T)\right)$. The
potential diverges as $h\to 0$, while $hm$ (where $m$ is the size of matrix)
remains finite. In this limit, for a wide class of potentials, the eigenvalues
of random normal matrices occupy a compact domain in the complex plane, called
a “droplet”. In the case of the confining potential created by the uniform
neutralizing background, the evolution of boundaries of this droplet with
respect to the matrix size, i.e. with respect to the time $T_{0}$, is a
solution to the Laplacian growth, or Hele-Shaw moving boundary problem of the
2D fluid dynamics (with the area of the droplet proportional to $T_{0}$, see,
e.g., [8, 19, 22, 21]).

The Laplacian growth is a model of evolution of a droplet whose boundary is
driven by a harmonic field being a potential for the growth velocity. The
field is a Green function of the exterior of the droplet which vanishes on the
droplet boundary. The normal velocity of the boundary is proportional to the
normal derivative of the field. The problem in which the droplet is expanding
is linearly unstable and ill-posed for almost any initial condition (see, e.g.
[7] and references therein). The finite $m$ normal matrix model provides a
sort of integrable regulaization of the above ill-posed problem. The
discretization of the matrix model, related to $N$-soliton solutions, could
provide another type of the regularization, a kind of the lattice
regularization. In this case, the continuous Laplacian growth is recovered in
a pair of double scaling limits $N\to\infty,\epsilon\to 0$, $m\to\infty,h\to
0$.

It is worth mentioning that apart from the standard normal matrix model, which
is a generalization of the hermitian random matrix model to normal matrices,
non-hermitian generalizations of the symmetric and quarternion-real models
were considered in the literature as well, see [24] and references therein.
These models, called “generalized Ginibre ensembles”, turn out to be related
to the so-called “large BKP” and “large 2-BKP” hierarchies. It would be
interesting to study solitonic $\tau$-functions of these hierarchies in the
context of the statistical mechanics of Coulomb gases with boundaries. Some
soliton solutions of the large (“fermionic”) BKP hierarchy were written down
in [26] and this may be useful for further studies.

It can be also interesting to try to apply our method to the Mehta-Pandey
interpolating ensembles, to circular $\beta=1,4$ ensembles [20] and also to
certain new solvable ensembles of random matrices such as recently considered
in [2] and [25].

Note that since the KP and 2DTL phase shifts are essentialy the same, the
Coulomb potential related to the 2DTL could be already obtained from the KP
phase shifts (3.1) by setting $a_{i}=\zeta_{i}/R$, $b_{i}=-R/\zeta_{i}$. One
might try to apply this choice to the BKP phase shifts (3.7). However, the
corresponding interaction potential has no clear physical meaning in this
case: here not only images which are reflections $\zeta_{i}\to
R^{2}/\bar{\zeta}_{i}$ with respect to the circle, but also the inversions
$\zeta_{i}\to-\zeta_{i}$ with respect to the origin are present (see Figure
2).

Figure 2: Positive charges/images are shown as white squares while negative
images are shown as black squares. The interactions between different charges
are shown by dashed lines, while the interactions between charges and their
own images are shown by solid lines.

Concluding this article we would like to remind that, while useful from the
point of view of random matrices and Dyson gases, the $\tau$-function approach
has an essential drawback in the framework of general statistical mechanics.
Namely, the partition functions derived from the $\tau$-functions of
integrable hierarchies correspond to the Coulomb gas (or Ising models [15,
16]) models at fixed (inverse) temperatures $\beta=2$. For the translationally
invariant Ising models with non-local interaction in one dimension, related to
the self-similar potentials [28] and certain soliton solutions of the KdV and
BKP hierarchies [15, 16], we have temperatures $\beta=2$ and $\beta=4$,
respectively. Restriction to fixed temperatures is a consequence of the fact
that the integrable hierarchies are nothing, but the Plucker relations on an
infinite dimensional Grassmanian [5, 12, 27] that can be obtained in the
framework of the free-fermion formalism.

Acknowledgments. This study has been partially funded within the framework of
the HSE University Basic Research Program. I. Loutsenko and O. Yermolaeva
would like to thank Centre de Recherches Mathématiques for a support. The
authors thank also the referees for their helpful remarks.

## References

  * [1] M. J. Ablowitz and H. Segur, Solitons and the Inverse Scattering Transform, SIAM, Philadelphia, 1981.
  * [2] N. Amburg, A. Yu. Orlov, D. Vasiliev, On products of random matrices, Entropy 22 (2020), no. 9, 972; doi:10.3390/e22090972
  * [3] L.-L. Chau and Y. Yu, Unitary polynomials in normal matrix model and wave functions for the fractional quantum Hall effect, Phys. Lett. A 167 (1992), 452–458.
  * [4] L.-L. Chau and O. Zaboronsky, On the structure of correlation functions in the normal matrix models, Commun. Math. Phys. 196 (1998), 203–247.
  * [5] E. Date, M. Jimbo, M. Kashiwara, and T. Miwa, Operator approach to the Kadomtsev-Petviashvili equation - transformation groups for soliton equations III, J. Phys. Soc. of Japan, 50 (11) (1981), 3806–3812.
  * [6] P. Forrester, Log-Gases and Random Matrices, London Math. Soc. Monographs, Princeton Univ. Press, 2010.
  * [7] B. Gustafsson, R. Teodorescu, and A. Vasil’ev, Classical and Stochastic Laplacian Growth, Advances in Mathematical Fluid Mechanics, Birkhauser, Cham, 2014.
  * [8] J. Harnad, I. Loutsenko, and O. Yermolayeva, Constrained reductions of two-dimensional dispersionless Toda hierarchy, Hamiltonian structure, and interface dynamics, J. Math. Phys. 46 (2005), 112701.
  * [9] J. Harnad and A. Yu. Orlov, Fermionic approach for evaluating integrals of rational symmetric functions, Theor. Math. Phys. 158 (2009), 17–39.
  * [10] R. Hirota, Direct methods of finding exact solutions of nonlinear evolution equations, in: Backlund Transformations, R. M. Miura, ed. Lecture Notes in Mathematics 515, Springer-Verlag, New York (1976), pp. 40–68.
  * [11] R. Hirota, Soliton solutions to the BKP equations. I. the Pfaffian technique, J. Phys. Soc. Japan 58 (1989), 2285–2296.
  * [12] M. Jimbo and T. Miwa, Solitons and infinite-dimensional Lie algebras, Publ. RIMS, European Math. Soc. Publ. House 19 (3) (1983), 943–1001.
  * [13] S. Kharchev, A. Marshakov, A. Mironov, A. Orlov, A. Zabrodin, Matrix models among integrable theories: Forced hierarchies and operator formalism, Nucl. Phys. 366 (1991),569–601.
  * [14] I. K. Kostov, I. Krichever, M. Mineev-Weinstein, P. B. Wiegmann, A. Zabrodin, $\tau$-function for analytic curves, in: Random Matrices and Their Applications (MSRI publications, vol. 40), Cambridge (2001), pp. 285–299; https://arxiv.org/abs/hep-th/0005259
  * [15] I. M. Loutsenko and V. P. Spiridonov, Self-similar potentials and Ising models, JETP Lett. 66 (1997), 747–753.
  * [16] I. M. Loutsenko and V. P. Spiridonov, Spectral self-similarity, one-dimensional Ising chains and random matrices, Nucl. Phys. B 538 (1999), 731–758.
  * [17] I. M. Loutsenko and V. P. Spiridonov, Soliton solutions of integrable hierarchies and Coulomb plasmas, J. Stat. Phys. 99 (2000), 751–767.
  * [18] I. M. Loutsenko and V. P. Spiridonov, A critical phenomenon in solitonic Ising chains, SIGMA 3 (2007), 059, 11 pp.
  * [19] I. Loutsenko and O. Yermolayeva, On integrability and exact solvability in deterministic and stochastic Laplacian growth, Math. Model. Nat. Phenom. 15 (2020), 3, 24 pp.
  * [20] M. L. Mehta, Random Matrices, Elsevier, 2004.
  * [21] M. Mineev-Weinstein, P. Wiegmann, and A. Zabrodin, Integrable structure of interface dynamics, Phys. Rev. Lett. 84 (2000), 5106–5109.
  * [22] M. Mineev-Weinstein, M. Putinar, and R. Teodorescu, Random matrices in 2D, Laplacian growth and operator theory, J. Phys. A: Math. Theor. 41 (2008), 263001.
  * [23] B. Nienhuis, Coulomb gas formulation of two-dimensional phase transitions, in: Phase transitions and critical phenomena 11 (1987), 1–53.
  * [24] A.Yu. Orlov, Deformed Ginibre ensembles and integrable systems, Phys. Lett. A 378 (2014) 319–328
  * [25] A. Yu. Orlov, New solvable matrix models III, https://arxiv.org/abs/2112.14741
  * [26] A. Yu. Orlov, T. Shiota, K. Takasaki, Pfaffian structures and certain solutions to BKP hierarchies I. Sums over partitions, https://arxiv.org/abs/1201.4518
  * [27] M. Sato, Soliton equations as dynamical systems on infinite dimensional Grassmann manifolds, Kokyuroku, RIMS, Kyoto Univ. 439 (1981), 30–46.
  * [28] V. Spiridonov, Universal superpositions of coherent states and self-similar potentials, Phys. Rev. A 52 (1995), 1909–1935, arXiv:quant-ph/9601030.
  * [29] K. Takasaki, Toda hierarchies and their applications, J. Phys. A: Math. Theor. 51 (2018), 203001.
  * [30] K. Takasaki and T. Takebe, Integrable hierarchies and dispersionless limit, Rev. Math. Phys. 7 (1995), no. 5, 743–808.
  * [31] A. Zabrodin, Canonical and grand canonical partition functions of Dyson gases as tau-functions of integrable hierarchies and their fermionic realization, Complex Analysis and Operator Theory 4 (2010), 497–514.