\section{Introduction}
In a recent article, Qin et al. \cite{QinPRL108} have investigated clustering
in low density nuclear matter. Equilibrium constants for $^{2}$H, $^{3}$H,
$^{3}$He and $^{4}$He have been measured following Guldberg and Waage law \cite{PaponLeblond} by selecting experimentally a sub-set of events
corresponding to a gas of neutrons and protons
in equilibrium with clusters.\par
Besides the fact that these results are important for the knowledge of the density dependence of the nuclear equation of state \cite{BaldoProgPartNuclPhys}, the measurements are related to density and temperature values that are of major importance for astrophysics since
nuclear equation of state plays a fundamental role in the understanding of core-collapse supernovae, mergers of compact stars and cooling proto-neutron stars for example.
In particular the chemical composition of proto-neutron stars influences the neutrino opacities and then their cooling \cite{HempelPRC91} \cite{PaisPRC97}.\par
This result therefore deserves to be confirmed or disproved using another experiment and equipment.\par
This paper presents the analysis framework of Qin et al. \cite{QinPRL108} adapted to our data concerning H and He isotopes. Then we raise the problem of cluster binding energy values which are assumed to be the vacuum values for the determination of the density which seems inconsistent to us.\par
In our experiment the use of four entrance channel systems with different neutron to proton ratios (N/Z) at the same bombarding energy, $^{136,124}$Xe+$^{124,112}$Sn, will allow to test different assumptions made during the analysis.
\section{Experimental details}
The $4\pi$ multi-detector INDRA
\cite{Pou95} was used to study four reactions with beams of $^{136}$Xe and $^{124}$Xe,
accelerated by the GANIL cyclotrons to 32 MeV/nucleon, and thin (530 $\mu$g/cm$^{2}$) targets of $^{124}$Sn and $^{112}$Sn.
INDRA is a charged product multidetector,
composed of 336 detection cells arranged in 17 rings centered
on the beam axis and covering 90\% of the solid angle.
The first ring (2$^{o}$ to 3$^{o}$) is made of 12 telescopes composed of
300 $\mu$m silicon wafer (Si) and CsI(Tl) scintillator (14 cm
thick). Rings 2 to 9 (3$^{o}$ to 45$^{o}$) are composed of 12 or 24
three-member detection telescopes : a 5 cm thick ionization
chamber (IC); a 300 $\mu$m or 150 $\mu$m silicon
wafer; and a CsI(Tl) scintillator (14 to 10 cm thick) coupled
to a photomultiplier tube.
Rings 10 to 17 (45$^{o}$ to 176$^{o}$) are composed of 24, 16 or 8 two-member telescopes:
an ionization chamber and a CsI(Tl) scintillator of 8, 6
or 5 cm thickness.
INDRA can identify in charge fragments from Hydrogen to Uranium and in mass light fragments with low thresholds.
Recorded event functionality was activated under a triggering factor based on a minimum number of fired
telescopes (M$_{trigger}^{min}$) over the detector acceptance (90\% of $4\pi$).
During the experiment
minimum bias (M$_{trigger}^{min}$=1) and exclusive (M$_{trigger}^{min}$=4) data were recorded.\par
This study is limited to the forward part of the centre of mass (hereinafter called c.m.)
for which excellent mass and charge identification performances are achieved for Hydrogen and Helium isotopes (hereinafter called lcp). Since only most violent collisions will be used in the following, only exclusive data are analysed.
\section{Event and sub-event selections}
In this paragraph, we present the event (central collisions) and the sub-event (c.m. angular cut) selections used in order to study a gas of neutrons and protons in equilibrium with clusters.\par
In a previous publication \cite{BougaultPRC97} concerning the same experiment, we have demonstrated that chemical equilibrium is achieved in central collisions. Therefore the most violent events, as in \cite{QinPRL108}, are employed in order to look for sub-events that correspond to a gas of neutrons and protons in equilibrium with clusters. Events with reduced impact parameters lower than 0.15 have been selected, using the impact parameter evaluator described in \cite{BougaultPRC97} (forward c.m. lcp total transverse energy greater than 200 MeV). As will be explained later, the gas of particles is assumed to have an isotropic momentum distribution. Since very central collisions have the advantage of minimizing the sideflow \cite{ReisdorfAnnRevNucl1997}, the centrality selection is thus appropriate.\par
In \cite{BougaultPRC97} we have also indicated for the studied systems the presence of two dominating lcp sources when analysing the forward part of the c.m. as a function of impact parameter : intermediate velocity (IV) and projectile-like (PL) sources. The PL source velocity is evolving with centrality while the IV source is located at the c.m. velocity The gas of neutrons and protons in equilibrium with clusters has to be looked for in the IV source \cite{QinPRL108} since the other source is mostly producing lcp by secondary decays.
A simple and efficient way to minimize contribution from PL source is to apply a 60$^{0}$-90$^{0}$ c.m. polar angular selection \cite{BougaultPRC97}.\par
In the following the data concerning lcp are selected via: (i) central events, (ii) emission almost perpendicular to the beam direction in the c.m.(60$^{0}$-90$^{0}$ polar angular range).
\section{Expanding gas source and surface velocity}
\begin{figure}[ht]
\begin{center}
\resizebox{0.48\textwidth}{!}
   \includegraphics{figure1a.eps}

   }
\resizebox{0.48\textwidth}{!}
   \includegraphics{figure1b.eps}

   }
\caption{Multiplicity of Z=1 (a) and Z=2 (b) isotopes as a function of their centre of mass energy for the four studied systems (60$^{0}$-90$^{0}$ polar angular range).}
\label{FigEcm}
\end{center}
\end{figure}
The framework of the analysis assumes that the characteristics of the selected lcp correspond to the outcome of a whole chain of evolution of the IV gas source. Clusters are formed during the space-time evolution of the expanding gas of nucleons.
It is assumed that, for each time step, neutrons and protons are in equilibrium with clusters.
Equilibrium means that
although each collision is a continuously evolving dynamical process with no notion of temporal equilibrium, after a certain time in each collision the composition of the gas of nucleons and light clusters is frozen and corresponds to that which is detected, this is called (chemical) freeze-out. Because the measured properties
\cite{BougaultPRC97} of the set of (sub-)events such as species multiplicity can be reproduced by assuming an unbiased population of the phase space available for the ensemble of freeze-out configurations produced by the collisions, then we may say that the observed (sub-)events are compatible with statistical equilibrium \cite{ColeISN1995} \cite{ChomazEPJA30} \cite{BorderiePPNP105}.\par
The key observable is the velocity (v$_{surf}$) of the particles in the IV frame prior to acceleration
by the Coulomb field of the remaining charged material
\cite{AwesPRC24}. Calculations \cite{WangPRC72} indicate that the surface velocity decreases with increasing average emission time and therefore v$_{surf}$ values may be used to select different time steps of the gas part evolution, i.e. of the clusterisation process, fastest particles corresponding to earliest emission times.\par
A fitting procedure of the observed energy spectra is performed in \cite{QinPRL108} in order to disentangle the different lcp productions and to deduce the Coulomb repulsion on lcp caused by the remaining charged matter. We have adopted another method since the barrier parameters are generaly poorly defined in a multi-fit procedure \cite{HagelPRC62}.\par
The characteristics of selected H and He elements are shown in figure \ref{FigEcm} through c.m. energy spectra for the four studied systems.
The spectra are double differential multiplicities. They are c.m. energy spectra for particles detected in the 60$^{0}$-90$^{0}$ c.m. polar angular range normalized to the number of central events for each studied system, normalized to the c.m. energy bin size and normalized to the 60$^{0}$-90$^{0}$ c.m. solid angle. This is then the number of $^{1}$H, $^{2}$H, $^{3}$H, $^{3}$He, $^{4}$He, $^{6}$He per event, per MeV and per steradian.
The chemistry of lcp production is system dependent since neutron-rich/poor isotope productions are following the neutron-richness/poorness of the entrance channel systems.\par
In table \ref{TableEc} we present a calculation using the spectra of figure \ref{FigEcm}.
Coulomb energy value per unit charge (E$_{C}$) is varied and
we calculate the particle emission Coulomb boost caused by the remaining charge after complete emission of the particles of the expanding gas. The average charge of the expanding gas in the 60$^{0}$-90$^{0}$ c.m. angular range is the charge integrated over the multiplicity spectra for E$_{c.m.}$ greater than Z times E$_{C}$ for each particle and then summed together. This average charge approximated to the nearest integer is presented in the second column.
The numbers quoted in the third column are the remaining charge over 4$\pi$ after the complete emission of the particles of the expanding gas assuming an isotropic emission in the c.m.; this charge is responsible for the Coulomb boost. Since we are considering an expanding gas of nucleons and clusters at very low densities therefore the boost from the gas itself is negligible and we only consider the boost from the remaining part which is considered to be at much higher density.
The amount of the boost per atomic charge (fourth column) is then calculated by using standard Coulomb barrier formula for surface emission of H or He from the remaining charge. Comparing the values between first and last columns which should be identical, only the third hypothesis is consistent (E$_{C}$=10 MeV).
The deduced Coulomb boost may be considered as an average value since, event by event, the size of the remaining charge varies.
\par
\begin{table}[h]
\caption{\label{TableEc}$\Sigma$Z$_{lcp}$: sum of atomic number of particles having c.m. energy greater than Z times E$_{C}$ (hypothesis) in the c.m. considered angular range. 104-4$\Sigma$Z$_{lcp}$: 4$\pi$ remaining atomic number in relation to the total system charge minus total atomic number of the gas. Last column: calculated Coulomb barrier per Z between lcp and 104-4$\Sigma$Z$_{lcp}$.}
\begin{indented}
\item[]\begin{tabular}{@{}llll}
\br
E$_{C}$&$\Sigma$Z$_{lcp}$&104-4$\Sigma$Z$_{lcp}$&E$_{C}$\\
hypothesis&[60$^{o}$,90$^{o}$]&4$\pi$&calculated\\
\mr
0 MeV&8&72&10 MeV\\
5 MeV&7&76&10 MeV\\
10 MeV&6&80&11 MeV\\
15 MeV&5&84&12 MeV\\
20 MeV&3&92&13 MeV\\
\br
\end{tabular}
\end{indented}
\end{table}
The key observable (v$_{surf}$) of the gas part is then calculated taking into account a Coulomb Barrier of 10 MeV per atomic number with an error estimated to be $\pm$ 2 MeV.\par
The sub-event selection is thus: (i) emission perpendicular to the beam direction in the c.m.(60$^{0}$-90$^{0}$ polar angular range), (ii) particles having c.m. energy greater than Z times 10 MeV.
The corresponding surface velocity spectra are presented in figure \ref{FigVsurf}, v$_{surf}$ being related to (E$_{c.m.}$-ZE$_{C}$)/A for a cluster (A, Z), with E$_{C}$=10 MeV.
\begin{figure}[ht]
\centering
\resizebox{0.48\textwidth}{!}
   \includegraphics{figure2a.eps}

   }
\resizebox{0.48\textwidth}{!}
   \includegraphics{figure2b.eps}

   }
   \caption{Multiplicity of Z=1 (a) and Z=2 (b) isotopes as a function of their reconstructed surface velocity for the four studied systems (60$^{0}$-90$^{0}$ polar angular range).}
\label{FigVsurf}
\end{figure}
\section{Volume in momentum space}
The final-state coalescence model \cite{ButlerPRL7} proposed to explain deuteron formation was phenomenologically extended to heavier fragments \cite{ScPR129}. It is assumed that nucleons coalesce when located together inside a volume in momentum space of size $4\pi P_{0}^{3}/3$, with $P_{0}$ a parameter. To take into account for the Coulomb field generated by the remaining charge we use the Coulomb corrected coalescence model formalism of \cite{AwesPRC24} to determine coalescence parameter $P_{0}$ for each isotope as in \cite{QinPRL108}.\par
In the c.m. frame the relationship between the observed cluster and proton differential cross sections is
\vspace{12pt}
\begin{equation}
\label{EqCoalescence}
\fl
\frac{d^{2}M(A,Z)}{dE_{c.m.}d\Omega_{c.m.}} =R_{np}^{N}\frac{A^{-1}}{N!Z!}
\left(\frac{4\pi P_{0}^{3}}{3[2m^{3}(E^{proton}_{c.m.}-E_{C})]^{1/2}}\right)^{A-1}
\left(\frac{d^{2}M(1,1)}{dE^{proton}_{c.m.}d\Omega_{c.m.}}\right)^{A}
\end{equation}
\vspace{12pt}
where the cluster differential multiplicity, M(A,Z), having a Coulomb corrected energy $E_{c.m.}-ZE_{C}$ is related to the proton
differential multiplicity, M(1,1), at the same Coulomb corrected energy per nucleon, E$^{proton}_{c.m.}$-E$_{C}$.
E$_{C}$ is the Coulomb repulsion per unit charge, m is the nucleon mass,
$(E_{c.m.}-ZE_{C})/A=E^{proton}_{c.m.}-E_{C}$ thus $E_{c.m.}=AE^{proton}_{c.m.}-NE_{C}$ for a nucleus composed of Z protons and N=A-Z neutrons and R$_{np}$ is the neutron to proton ratio. The neutron spectra is assumed to be identical to the proton spectra once the Coulomb correction has been applied.\par
Contrary to the original coalescence analysis \cite{NagamiyaPRC24} which determined an average value of P$_{0}$, here the goal is to characterize the evolution of the volume of the gas as a function of time, therefore of v$_{surf}$ or (E$_{c.m.}$-ZE$_{C}$)/A.\par
Neutrons are not measured in our experiment but hypotheses concerning the neutron to proton ratio may be tested by comparing results given by the four entrance channels ($^{136,124}$Xe+$^{124,112}$Sn) since P$_{0}$ should not depend on them. In a first step we used equation \ref{EqCoalescence} regardless of the neutron to proton ratio (R$_{np}$=1). The parameter P$_{0}$ as a function of (E$_{c.m.}$-ZE$_{C}$)/A is presented in figure \ref{FigP0}(a). We see that P$_{0}$ depends on cluster size and we observe large differences for the four different entrance channels.\par
\begin{figure}[ht]
\centering
\resizebox{0.48\textwidth}{!}
   \includegraphics{figure3a.eps}

   }
\resizebox{0.48\textwidth}{!}
   \includegraphics{figure3b.eps}

   }
   \caption{P$_{0}$ as a function of cluster energy per nucleon in the source frame before Coulomb acceleration. Results for the four systems and the five isotopes are shown. R$_{np}$=1 (a). R$_{np}$ from free neutrons and protons (b).}
\label{FigP0}
\end{figure}
Following the coalescence model \cite{CsernaiPR131} and as in \cite{QinPRL108}, the unmeasured neutron multiplicity at a given time, determined by v$_{surf}$ or (E$_{c.m.}$-ZE$_{C}$)/A, is given by the proton multiplicity multiplied by $^{3}$H/$^{3}$He multiplicity ratio \cite{FamianoPRL97} for that time. The question asked here is which proton to neutron ratio should be applied for each time step: the global ratio taking into account all species or the ratio calculated with only free nucleons. It turns out that the only ratio which preserves entrance channel independence is the ratio calculated with free neutrons and protons as shown in figure \ref{FigP0}(b).
In the coalescence model the proton cross section used is the original production, before any composites are formed.
The observed proton cross section, figure \ref{FigEcm}, is the original one minus the protons bound in the clusters.
This is also true for the neutron to proton ratio. The observed proton cross section is not the original one and figure \ref{FigP0} demonstrates the validity of using observed nucleon spectra. There is therefore a contradiction. It may be that by using the original nucleon spectra the appropriate neutron to proton ratio is also the original ratio, but this is impossible to prove experimentally since original proton energy spectra cannot be reconstructed.\par
The chemical equilibrium model \cite{MekjianPRC17} \cite{MekjianPLB89} also predicts a power law for the momentum space density of a composite nucleus relative to protons but is not affected by the previous arguments, since the law of mass action involves the free proton and neutron concentrations already reduced by composite formation.\par
The results presented in this paragraph thus confirm the validity of the thermal model. A volume, V$_{0}$, which represents the spatial region (assumed spherical) over which chemical equilibrium is established will be calculated for each isotope using differential cross section power law\cite{MekjianPLB89} \cite{DasGuptaPR72}:
\vspace{12pt}
\begin{equation}\label{EqModeleEquilibre}
\fl
\frac{d^{3}M(A,Z)}{d^{3}p_{A}} =R_{np}^{N}\frac{(2s+1)~e^{B(A,Z)/T}}{2^{A}}\left(\frac{h^{3}}{V_{0}}\right)^{A-1}
\left(\frac{d^{3}M(1,1)}{d^{3}p}\right)^{A}
\end{equation}
\vspace{12pt}
where h is Planck's constant, s the spin of the cluster (Z, A, N=A-Z), B(A,Z) its ground state binding energy, T the temperature of the statistical ensemble describing the gas, $p$ the free proton momentum and $p_{A}=Ap$ the cluster momentum, both in the source frame before Coulomb acceleration ($p=mv_{surf}$).\par
The volume V$_{0}$ for each isotope (A,Z) is related to the volume in momentum space $4\pi P_{0}^{3}/3$ (P$_{0}$ of figure \ref{FigP0}(b)) by the relationship
\vspace{12pt}
\begin{equation}\label{EqVolume}
\fl
V_{0}=\frac{3h^{3}}{4\pi P_{o}^{3}}\left[(2s+1)\left(\frac{Z!N!A^{3}}{2^{A}}\right)e^{B(A,Z)/T}\right]
^{1/(A-1)}
\end{equation}
\vspace{12pt}
\section{Temperature and surface velocity confidence interval}
In this paragraph we will study the temperature evolution of the expanding gas source as a function of v$_{surf}$.
\par
The determination of the volume with equation \ref{EqModeleEquilibre} requires the measurement of the temperature. Following the statistical approach of \cite{AlbergoNCA89}, as in \cite{QinPRL108}, the temperature for a given v$_{surf}$ value is deduced from the multiplicities, M(A,Z), of $^{2}$H, $^{3}$H, $^{3}$He, $^{4}$He at the same v$_{surf}$ value with B(A,Z) the ground state binding energy of isotope (A,Z)
\vspace{12pt}
\begin{equation}\label{EqT}
\fl
T =\frac{B(4,2)+B(2,1)-B(3,2)-B(3,1)}{\ln(\sqrt{9/8}(1.59~R_{v_{surf}}))} MeV
~{\rm with}~R_{v_{surf}}=\frac{M(2,1)M(4,2)}{M(3,1)M(3,2)}
\end{equation}
\vspace{12pt}
This relationship, valid for particles emitted from a single source at temperature T and having a Maxwellian spectrum, is applied for each bin of v$_{surf}$.\par
The deduced temperature values as a function of v$_{surf}$ are shown in figure \ref{FigTemperature}. Temperatures are independent of the entrance channel, $^{136,124}$Xe+$^{124,112}$Sn, as expected since the bombarding energy is identical. This result was not self-evident because the multiplicities, especially $^{3}$H and $^{3}$He, are very system dependent (see figure \ref{FigVsurf}).
For V$_{surf}$ below 6.5 cm/ns, $^{2}$H and $^{4}$He productions are entrance channel system independent, the spectra are the same for the four studied systems. It is not the case for $^{3}$H and $^{3}$He productions; as expected $^{3}$H production is growing with neutron richness of the system while $^{3}$He is decreasing.
The temperature depends on $^{4}$He, $^{2}$H, $^{3}$He and $^{3}$H productions. Since only $^{3}$H and $^{3}$He productions depend on the entrance channel N/Z, therefore in order for the temperature to be identical for the four studied systems $^{3}$H and $^{3}$He multiplicities must counterbalance each other almost exactly.
This result further strengthens the validity of the equilibrium model for describing the data. For higher v$_{surf}$ values the error bars prevent us from drawing definite conclusions.\par
\begin{figure}[ht]
\centering
\resizebox{0.48\textwidth}{!}
   \includegraphics{figure4.eps}

   }
   \caption{Temperature versus surface velocity for the four systems.}
   \label{FigTemperature}
\end{figure}
The evolution of T as a function of v$_{surf}$ is similar to that presented in Figure 3 of \cite{WangPRC72} and figure 1 in \cite{KowalskiPRC75}. In our case the temperature is constant for v$_{surf}$ below 3 cm/ns, then it increases to reach a maximum around 9 MeV for v$_{surf}$ about 6.5 cm/ns. Then it seems to decrease with increasing v$_{surf}$. Figure \ref{FigTemperature} indicates that the scenario of the IV source which cools with time is valid for v$_{surf}$ greater than 3 cm/ns and lower than 6.5 cm/ns. The high limit is adopted for lack of cluster population while below 3 cm/ns several lcp production mechanisms are probably present.\par
In view of the results mentioned in this paragraph, the surface velocity interval between 3 and 6.5 cm/ns represents the confidence interval for the subsequent analysis.\par
\section{Volume from thermal model}
In this paragraph we will study the volume evolution of the expanding gas source as a function of v$_{surf}$ using equation \ref{EqModeleEquilibre} and \ref{EqT}.\par
\begin{figure}[ht]
\centering
\resizebox{0.48\textwidth}{!}
   \includegraphics{figure5.eps}

   }
   \caption{Volume from thermal model as a function of surface velocity for the five isotopes ($^{124}$Xe+$^{124}$Sn). The lines represent two fits of the average volumes of A$\ge$3. The full line is a fit of data using the ratio of free neutrons to free protons with equation \ref{EqNsurZfree} whereas the dashed line is a fit of data without the temperature dependence of free neutrons to free protons ratio.}
   \label{FigVolume}
\end{figure}
The results are presented in figure \ref{FigVolume}. The gas source volumes are increasing with decreasing v$_{surf}$. This observation is compatible with the picture of an evolving expanding source with time.
There exists a dependence on the composite particle which has been interpreted in terms of the finite size of the composite \cite{DasGuptaPR72}.\par
In the coalescence model the free neutron to free proton ratio is given by the $^{3}$H/$^{3}$He multiplicity ratio. In the equilibrium model \cite{AlbergoNCA89} this ratio (R$_{np}$ of equation \ref{EqModeleEquilibre}) is given by
\vspace{12pt}
\begin{equation}\label{EqNsurZfree}
\fl
R_{np}=\frac{M(3,1)}{M(3,2)}~e^{((B(3,2)-B(3,1))/T)}
\end{equation}
\vspace{12pt}
where T is the temperature and B(A,Z) the ground state binding energy of isotope (A,Z).\par
The volume evolution of the expanding gas source as a function of v$_{surf}$ is calculated using equation \ref{EqNsurZfree} for R$_{np}$. Without this temperature dependence contained in the exponential term of equation \ref{EqNsurZfree} the deduced volumes of $^{3}$H and $^{3}$He are different. The validity of introducing this temperature dependence is confirmed by Relativistic Mean-Field approach calculations \cite{PaisPRC97} where deviations of about 20\% concerning the free neutron proton ratio can occur for temperatures around 4 MeV if the exponential term of equation \ref{EqNsurZfree} is not introduced \cite{PaisPrivate}. This dependence on temperature is a difference as compared to \cite{QinPRL108}.\par
The volumes extracted for deuterons are larger as in \cite{QinPRL108}. This could be explained by the fragility of the deuteron and its survival probability once formed \cite{QinPRL108}. Therefore, as in \cite{QinPRL108}, the average volumes over which chemical equilibrium is established are derived from A$\ge$3 clusters. The average volumes are represented by the full line in figure \ref{FigVolume}. The full line is the result of a fitting procedure of the average volumes using a Landau function for each studied system. The dashed line is the Landau fit of average volumes without the temperature dependence of R$_{np}$. We see that the temperature dependence of equation \ref{EqNsurZfree} is small for high v$_{surf}$ and becomes important for low v$_{surf}$ where the temperature is the lowest.\par
\section{Density}
In this paragraph we will study the density evolution of the expanding gas source as a function of v$_{surf}$.\par
The average volumes of figure \ref{FigVolume} are free volumes since they were calculated in a framework of point-like clusters. In order to deduce the total average volumes it is necessary to add the excluded volumes represented by the volumes occupied by the clusters themselves. The density is related to total average volume.\par
The evolving source volume for each bin of v$_{surf}$, for a given time t, is
\vspace{12pt}
\begin{equation}\label{EqVolumeTotal}
\fl
V_{t} = V_{0}+\frac{4}{3}\pi r_{0}^{3}A_{t}
\end{equation}
\vspace{12pt}
where V$_{t}$ is the total volume of the gas at time t, V$_{0}$ is the volume deduced from the thermal model at time t, A$_{t}$ is the total mass of the gas at time t and r$_{0}$ is the average particle radius (1.3 fm).
The density at time t is $\rho_{t}=A_{t}/V_{t}$.\par
\begin{figure}[ht]
\centering
\resizebox{0.6\textwidth}{!}
   \includegraphics{figure6a_6b.eps}

   }
\resizebox{0.6\textwidth}{!}
   \includegraphics{figure6c_6d.eps}

   }
\resizebox{0.6\textwidth}{!}
   \includegraphics{figure6e_6f.eps}

   }
   \caption{Fits of E$_{cm}$ isotope spectra: $^{1}$H (a), $^{2}$H (b), $^{3}$H (c), $^{3}$He (d), $^{4}$He (e), $^{6}$He (f). Full lines are the fits performed independently for each isotope, the dashed lines correspond to v$_{surf}$=3 cm/ns. $^{124}$Xe+$^{124}$Sn system.}
   \label{FigFitsEcm}
\end{figure}
The knowledge of A$_{t}$ at time t is related to the knowledge of A$_{source}$, the mass of the gas source at the beginning of the expansion process, through $A_{t}=A_{source}-A_{exp}$, where A$_{exp}$ is the mass expelled before time t. In other words, for a given v$_{surf}$, A$_{t}$ is the total mass integrated from 0 to v$_{surf}$.\par
It has been shown that the v$_{surf}$ interval between 0 and 3 cm/ns does not match with the picture of a cooling while expanding gas (confidence interval). Therefore we cannot rely directly on the data in this velocity range to calculate A$_{t}$. A$_{source}$ has thus been evaluated using a fitting procedure of E$_{cm}$ spectra of figure \ref{FigEcm}.\par
The spectra were fitted independently by an expanding Boltzmann distribution \cite{BondorfNPA296}
\vspace{12pt}
\begin{equation}\label{EqFit}
\fl
\frac{d^{2}M(A,Z)}{dE_{c.m.}d\Omega_{c.m.}} =C~e^{-(E_{c.m.}-ZE_{C}+E_{Reff})/T_{eff}}
~sinh(2\sqrt{(E_{c.m.}-ZE_{C})E_{Reff}}/T_{eff})
\end{equation}
\vspace{12pt}
where E$_{C}$ is the Coulomb repulsion per unit charge, E$_{Reff}$ is related to the radial expansion energy of the source, T$_{eff}$ is the temperature and C a normalization factor.\par
The fits are performed for E$_{cm}$ values corresponding to v$_{surf}$ greater than 3 cm/ns and E$_{C}$ is set at 10 MeV, the other three parameters being left free. E$_{Reff}$ and T$_{eff}$ are effective parameters since they represent an integration over the whole expansion process.\par
\begin{figure}[ht]
\centering
\resizebox{0.48\textwidth}{!}
   \includegraphics{figure7.eps}

   }
   \caption{Density versus surface velocity for the four studied systems.}
   \label{FigDensity}
\end{figure}
The results of the fitting procedure are given in figure \ref{FigFitsEcm} for the $^{124}$Xe+$^{124}$Sn system.
We see that the spectra are quite well reproduced for the v$_{surf}$ confidence interval whose lower limit is represented by the dashed line.
It should be mentioned that the proton spectrum is not so well reproduced between v$_{surf}$=3 and 4 cm/ns, which would tend to decrease the confidence interval. However, the latter has not been modified.
Results for the three other studied systems are similar.\par
The gas source size, integrated over 4$\pi$, is calculated using fit integration. v$_{surf}>0$ neutron multiplicity is deduced from equation \ref{EqNsurZfree} and as far as their contribution to the source is concerned, they are assigned the percentage of protons. The gas source size (A$_{source}$) is $51\pm2.2$, $53\pm1.6$, $49\pm1.6$ and $51\pm1.5$ for $^{124}$Xe+$^{124}$Sn, $^{136}$Xe+$^{124}$Sn, $^{124}$Xe+$^{112}$Sn, $^{136}$Xe+$^{112}$Sn respectively.
With regard to the confidence interval defined above (v$_{surf}$ values between 3 cm/ns and 6.5 cm/ns), as the mass of the evolving source (A$_{t}$) for a given v$_{surf}$ is taken as the total integrated mass between 0 and v$_{surf}$, it appears that the integrated mass between 6.5 cm/ns and infinity will have no influence on the following results.
\par
The density values ($\rho_{t}=A_{t}/V_{t}$) are presented in figure \ref{FigDensity} as a function of surface velocity in the confidence interval.
Density is increasing with surface velocity, i.e. is decreasing with time. As in \cite{QinPRL108} very low values are achieved.\par
Concerning the 4$\pi$ isotropic distribution of the gas source, we should note that this requirement could be restricted to the 60$^{o}$-90$^{o}$ angular range ($\pi$) without changing density values, the source size being divided by 4 as well as the volume.
\section{Gas source characteristics}
In this paragraph we present the characteristics of the gas source during its evolution.\par
The gas source characteristics for each time step are given by its temperature, its composition and its density.
\begin{figure}[ht]
\centering
\resizebox{0.48\textwidth}{!}
   \includegraphics{figure8a.eps}

   }
\resizebox{0.48\textwidth}{!}
   \includegraphics{figure8b.eps}

   }
\caption{Gas source characteristics. (a): temperature versus density for the four studied systems, squares are results from \cite{QinPRL108} with uncertainties represented by the lines.
(b): proton fraction as a function of density for the four studied systems; for clarity errors on density are not drawn. Density errors are indicated in figures \ref{FigDensity} and \ref{FigTvsDensity}(a).}
   \label{FigTvsDensity}
\end{figure}
In fact, the value of the temperature and the composition of the evolving source for each time step is given by the temperature and the composition calculated with lcp multiplicities contained in each bin of v$_{surf}$. Only the calculation of density requires the knowledge of the total mass of the source for each time step (A$_{t}$). The temperature is calculated with equation \ref{EqT} and the neutron multiplicity from the proton multiplicity through equation \ref{EqNsurZfree}.\par
The density dependence of the temperature is given in figure \ref{FigTvsDensity}(a) for the four studied systems. The temperature, ranging from 5 to 9 MeV, does not depend on the entrance channel. The figure presents the thermodynamical path of the cooling of the source. The INDRA data are confronted to results from \cite{QinPRL108} obtained at higher bombarding energy (47A MeV). This difference in beam energy may explain the gap between the two thermodynamical paths. For very low densities, the temperatures are closer. We may also attribute these differences to different collision dynamical paths and therefore different characteristics for the expanding gas present at mid-rapidity since in our case the entrance channel system is quasi-symmetric (Xe+Sn) while in the case of Qin et al. it is asymmetric (Ar+Sn and Zn+Sn).\par
The composition of the evolving source is given by the proton fraction, Z/A, whose density dependence is shown in figure \ref{FigTvsDensity}(b). The proton fraction values depend on the projectile plus target proton fractions, we note that they are identical for $^{136}$Xe+$^{112}$Sn and $^{124}$Xe+$^{124}$Sn systems in agreement with chemical equilibrium hypothesis.
The entrance channel proton fraction values are reached for low densities while there appears to be a trend towards larger proton fraction values as the density increases. This trend is larger for proton rich entrance channel systems.
This density dependence is essentially due to the $^{4}$He contribution to the global proton fraction. The $^{4}$He contribution evolves from 0.1 to 0.5 as the density decreases for the four studied systems, $^{4}$He clusters being more abundant for low densities. Therefore the global proton fraction is decreasing with decreasing density.\par
\section{Equilibrium constants}
In this paragraph we will compare the INDRA equilibrium constants to the data of Qin et al. \cite{QinPRL108}.\par
Equilibrium constant for a cluster of mass A and atomic number Z is defined as
\vspace{12pt}
\begin{equation}\label{EqKc1}
\fl
K_{c}(A,Z)=\frac{\rho_{pa}(A,Z)}{\rho_{pa}(1,1)^{Z}~\rho_{pa}(1,0)^{N}}
\end{equation}
\vspace{12pt}
where $\rho_{pa}(A,Z)$ is the (A,Z) particle partial density ($\rho_{t}=\sum_{i} \rho_{pa}(A_{i},Z_{i})$). K$_{c}$ should depend only on density and temperature of the equilibrated gas source which makes it a universal parameter. The isotope equilibrium constants are shown in figure \ref{FigKcvsDensity}(a) for the four studied systems as a function of the gas source density ($\rho_{t}$).\par
\begin{figure}[ht]
\centering
\resizebox{0.48\textwidth}{!}
   \includegraphics{figure9a.eps}

   }
\resizebox{0.48\textwidth}{!}
   \includegraphics{figure9b.eps}

   }
   \caption{Equilibrium constants as a function of density. (a): for the studied isotopes and for the four studied systems. (b): $^{4}$He Equilibrium constant as a function of density, three results are shown for the $^{124}$Xe+$^{124}$Sn data to take into account for the error on the Coulomb barrier used to calculate v$_{surf}$, the squares are the Qin et al. \cite{QinPRL108} data with the uncertainties represented by the lines.}
   \label{FigKcvsDensity}
\end{figure}
K$_{c}$ is easily reformulated as a function of mass fractions $w(A,Z)$
\vspace{12pt}
\begin{equation}\label{EqKc2}
\fl
K_{c}(A,Z) = \frac{w(A,Z)}{w(1,1)^{Z}~w(1,0)^{(A-Z)}}~\rho_{t}^{-(A-1)}
~{\rm with}~w(A,Z)=\frac{A~M(A,Z)}{\sum_{i}{A_{i}~M(A_{i},Z_{i})}}
\end{equation}
\vspace{12pt}
M(A,Z) being the multiplicity of particle (A,Z) calculated in a given v$_{surf}$ bin.\par
This to say that K$_{c}$ has an explicit dependence on $\rho_{t}^{-(A-1)}$ which explains its global density behaviour.
From figure \ref{FigKcvsDensity} we note that K$_{c}$ is independent of A/Z entrance channel within the error bars.\par
In figure \ref{FigKcvsDensity}(b), K$_{c}$ density dependence is presented for $^{4}$He. In order to take into account for the error on the Coulomb barrier per atomic number used to calculate v$_{surf}$ (E$_{C}=10\pm 2$ MeV) we present three results (E$_{C}=8, 10~$and$~12$ MeV) whose superposition indicates the final error on our data.\par
In figure \ref{FigKcvsDensity}(b) our data is also compared to results of Qin et al. \cite{QinPRL108}. Since the thermodynamical paths (figure \ref{FigTvsDensity}(a)) are not the same, we do not expect the two data sets to overlap. This is indeed the case. At very low densities, where the temperatures of the two experiments are close, we could have expected a better agreement, however, as we indicated earlier, the confidence interval below v$_{surf}=4~cm/ns$ (density lower than 0.008 nucleon/fm$^{3}$) is questionable for our data. Also we remind the reader that Qin et al. did not use equation \ref{EqNsurZfree} whose impact is important at low temperatures. Therefore there is no contradiction between the two results. However the full compatibility test has to be done with a model comparison.\par
\section{The problem of binding energies}
In this paragraph we will comment the method used to extract the equilibrium constants and the density values.\par
Model comparisons with \cite{QinPRL108} data have been made \cite{HagelEurPhys50} \cite{HempelPRC91} \cite{PaisPRC97} and all the studies conclude that the measured equilibrium constants are not reproduced by an ideal gas modelisation.
Medium effects are present and the consequence is the existence of a shift of the cluster binding energies \cite{RopkeNPA867} \cite{RopkePRC92} \cite{TypelPRC81} for $\rho > 0$.
The decrease in the binding energy arises from the indistinguishability between the nucleons inside the
clusters and the free nucleons (Pauli blocking).
The procedure ignores the reduction in the binding energy of the clusters as
the density of the surrounding medium increases because it assumes values of the cluster binding energy in vacuum
in equations \ref{EqModeleEquilibre}, \ref{EqVolume}, \ref{EqT} and \ref{EqNsurZfree}.
\par
In summary, chemical constants measured by Qin et al. were compared to statistical models employing different treatments of the nucleon-nucleon and nucleon-cluster interactions with the purpose of constraining the expected cluster binding energy shifts in dense matter. Formulae used to extract the measured quantities, temperature and volume and free neutron number, explicitly assume that the cluster abundancies are uniquely governed by their vacuum binding energies, which is in contradiction with the purpose of the analysis. Furthermore in reference \cite{HempelPRC91} it is shown that the Qin et al. experimental points are not reproduced by an ideal gas hypothesis although all the experimental analysis is based supposing an ideal gas of classical clusters in thermodynamic equilibrium at temperature T in the grand-canonical ensemble. Indeed if in medium corrections were negligible, the measured chemical constants would agree with the ideal gas prediction.
There is therefore a contradiction in the method. This point has to be clarified.\par
\section{Conclusions}
Equilibrium constants as a function of density has been presented in the framework of the study made by Qin et al. \cite{QinPRL108}.
A subset of events has been selected whose characteristics are consistent with a chemically- and thermally-equilibrated expanding gas of nucleons and clusters. The power law for the momentum space density of composite nuclei relative to protons has been used to extract the source volumes. It turned out that the chemical equilibrium model is more adapted to our data as compared to the coalescence model. The temperature values of the expanding source have been deduced using nuclear statistical equilibrium framework and entrance channel independence ($^{136,124}$Xe+$^{124,112}$Sn) of the result confirmed the validity of its use and allowed us to define a confidence interval for the subsequent analysis. The volume values of the expanding source have been extracted and
an improvement, confirmed by theorical calculations, has been applied concerning the relationship between neutron and proton multiplicities as compared to the original work of Qin et al. Those volume values have been transformed in density values by determining the total source size. The thermodynamical path covered by the expanding source in the temperature-density plane is different to that of Qin et al. This is not abnormal because the bombarding energies of the two experiments are not the same. Equilibrium constant values as a function of density were presented for $^{2}$H, $^{3}$H, $^{3}$He, $^{4}$He, $^{6}$He and the result is not in contradiction with the results of Qin et al.\par
Finally, we raised the problem of the values of cluster binding energies. Indeed, the values used in the analysis are those in vacuum (zero density) while all comparisons with the models indicate the existence of a shift towards lower values and existence of in-medium effects. This point has to be clarified.
\section*{References}


\section{Introduction}\label{sec:intro}
Real-world environments are always contain stationary and/or time-varying noises that are received together with speech signals by recording devices \cite{michelsanti2017conditional}. The received noises inevitably degrade the performance of multi-channel (MC)-based human--human and human--machine interfaces, and this issue has attracted significant attention over the years \cite{8683662,bitzer2001multi,8683732,qi2020tensor,bitzer1999multimicrophone,liu1995room,grais2018raw}. In recent decades, numerous MC speech-enhancement (SE) approaches have been proposed to alleviate the effect of noise and improve the quality and intelligibility \cite{de2017adaptive,bagheri2019exploiting,8683768,1255457,heymann2016neural,chen2018building} of received speech signals. In general, most of these approaches have been proposed for use in a microphone array architecture, wherein multiple microphones are compactly placed in a small space. For example, the beam-forming algorithm, one of the most popular methods that exploit the spatial diversity of received signals to design a linear filter in the frequency domain, aims to preserve the signal received from the target direction while attenuating noise and interference from other directions \cite{hoshuyama1999robust,emura2018distortionless,chang2019mimo}. Recently, several novel approaches have combined deep-learning-based algorithms, such as deep neural networks (DNNs) and convolutional neural networks (CNNs) \cite{wang2018supervised}, with the beam-forming process to further promote the enhanced capability of an MC SE system \cite{zhang2017deep,wang2018all,tan2019real,araki2015exploring} so as to produce speech signals with even higher quality. In addition to beam-forming-based approaches, in \cite{tawara2019multi}, the multiple recordings are directly enhanced in the time domain along the specified spatial direction through a denoising auto-encoder, and this method benefits automatic speech recognition (ASR) systems by reducing recognition errors.

In contrast, some researchers pay more attention on performing SE on the architecture of distributed-microphones \cite{liu2019multichannel,hassani2017real,Matheja2013,li2018distributed,TU201596}. This physical configuration, consisting of many individual self-powered microphones or microphone arrays, can be deployed in a large area \cite{7503169,6101302}. Therefore, one or more received signals with a higher signal-to-noise ratio (SNR) and a direct-to-reverberant ratio can be used for the distributed-microphone enhancement system to achieve better sound quality and intelligibility.

In general, a fusion center (FC) and a distributed signal processor (also called ad-hoc) are two alternative forms used in the distributed-microphone architecture \cite{markovich2015optimal}. For the FC, each recording device can transmit the recorded sounds to a powerful central processor that aims at reconstructing nearly clean speech relative to the selected target speaker. Some successful approaches associated with this architecture include robust principal component analysis \cite{li2018distributed}, generalized eigen value decomposition \cite{markovich2009multichannel}, MC Wiener filter \cite{serizel2014low}, and optimal MC frequency domain estimators \cite{TRAWICKI2012345}. In comparison, the ad-hoc-based enhancement system, however, enhances the noisy input locally in the individual device and then shares the result with its neighbors for further refinement. Some well-known techniques of this type include distributed linearly constrained minimum variance \cite{bertrand2013distributed}, linearly constrained distributed adaptive node-specific signal estimation \cite{bertrand2011distributed}, distributed generalized sidelobe canceler \cite{markovich2012distributed}, and distributed maximum signal to interference-plus-noise filtering \cite{tavakoli2017distributed}.

In this study, three novel SE systems based on DNNs are introduced and investigated for the distributed-microphone architecture. For the first system, we train the DNN model for each microphone channel to enhance the corresponding noisy recordings. In other words, a large MC SE system is divided into several single-channel noise reduction tasks, and this system is called ``DNN--S''. Next, the second system follows a process similar to the work in \cite{tawara2019multi}, wherein the multiple noisy utterances received are transmitted to an FC and then aggregated and used as input to a DNN model for producing the final enhanced signals. Because an FC is used here, this system is called ``DNN--F''. Finally, the third system comprises two operational stages. The first stage acts locally in each channel device by enhancing the noisy input with a DNN, and in the second stage, all the enhanced local-channel signals are combined and further processed with a fusion DNN. We call this system ``DNN--C'' because it nearly combines the two previous systems to facilitate the enhancement. Notably, in our opinion the main contribution of this study is two-fold: First, we show that the deep learning-based algorithm can be successfully applied to the distributed microphone architecture and reveals its promising capability in reducing noise effect. Second, three types of deep neural network (DNN) structures are compared with each other when they are adopted in the distributed-microphone SE task, and we investigate the underlying reasons for the corresponding performance difference.

The rest of this paper is organized as follows. Section \ref{sec:threemodel} introduces the aforementioned three systems, namely DNN--S, DNN--F, and DNN--C. Experiments and the respective analysis are given in Section \ref{sec:experi}. Section \ref{sec:conclu} presented conclusions and a future avenue.

\section{Three different multi-channel enhancement system}\label{sec:threemodel}
\begin{figure}[!t]
{\centering \includegraphics[width=\columnwidth]{./ ovallflow.eps}
\caption{Block diagram of a DNN-based speech-enhancement system with a distributed-microphone architecture.}\label{fig:ovall}}
\end{figure}
Figure \ref{fig:ovall} shows the general diagram of a distributed-channel SE architecture common to the three presented systems. The original clean signal shown in this figure is first corrupted by the background diffuse noise and/or reverberation and is received by the distributed microphones. Next, the short-time Fourier transform (STFT) and logarithmic operation are performed on the received signals to obtain the log-power magnitude spectra (LPS) in the speech signal processor block. The MC SE system based on DNN, DNN--S, DNN--F, or DNN--C, is then used to generate the enhanced version from the noisy LPS input. Finally, in the enhanced signal processor block, the inverse STFT (ISTFT) is applied to the enhanced LPS, together with the original phase component preserved from the specific channel to provide the final enhanced speech waveform. In the following three sub-sections, we will introduce the detail process of the aforementioned three DNN-based SE systems.

\subsection{DNN--S enhancement systems}
\begin{figure}[!t]
{\centering \includegraphics[scale=0.95]{./ dnn_s.eps}
\caption{Block diagram of the DNN--S SE system.}\label{fig:dnns}}
\end{figure}

For the DNN--S architecture as shown in Fig. \ref{fig:dnns}, a DNN, presented by DNN--S$_p$, was trained for each of the $P$--microphone channels using the channel-wise noisy speech feature set $\mathbf{Y}_p$ with the respective clean counterpart $\mathbf{T}$, where $P$ is the total number of microphone channels, and $p$ is the channel index. Consider an $L$--layer DNN--S$_p$, an arbitrary $l$th layer is formulated in Eq. \eqref{eq:dnns} in terms of the input-output relationship ($\mathbf{z}^{(l-1)}$, $\mathbf{z}^{(l)}$):
\begin{equation}\label{eq:dnns}
\mathbf{z}^{(l)}=\sigma^{(l)}\left(h^{(l)}(\mathbf{z}^{(l-1)})\right),\quad l=1,\cdots,L,
\end{equation}
where $\sigma^{(l)}$ and $h^{(l)}$ are the ReLU activation and linear transformation functions, respectively. Notably, the input and output layers correspond to the first and $L$-th layers, respectively. In addition, for DNN--S$_p$, we have $\mathbf{z}^{(0)}=\mathbf{Y}_p$ and $\mathbf{z}^{(L)}=\mathbf{E}_p^S$, where $\mathbf{E}_p^S$ is the ultimate output of this system. The DNN parameters are obtained by means of supervised training that minimizes the mean squared error (MSE) between $\mathbf{E}_p^S$ and the noise-free counterpart $\mathbf{T}$.

\subsection{DNN--F enhancement system}
\begin{figure}[!t]
{\centering \includegraphics[scale=0.95]{./ dnn_f.eps}
\caption{Block diagram of the DNN--F speech-enhancement system.}\label{fig:dnnf}}
\end{figure}
Figure \ref{fig:dnnf} illustrates the DNN--F block diagram. According to this figure, we collect the channel-wise noisy features, $\mathbf{Y}_p$, $p=1, 2, ..., P$ from all of the $P$--microphone channels, and concatenate them to form a long feature $\mathbf{Y}'_p$ in a frame-wise manner, i.e., $\mathbf{Y}'=[\mathbf{Y}_1;\cdots,\mathbf{Y}_p;\cdots,\mathbf{Y}_P]$. Then, the long features in the training set are used to train a DNN model in FC, which can be formulated as follows
\begin{equation}\label{eq:dnnf}
\mathbf{E}^F=\mbox{\textit{DNN-F}}\{\mathbf{Y}'\},
\end{equation}
where $\mbox{\textit{DNN--F}}\{\cdot\}$ represents the operation of the used DNN model, and $\mathbf{E}^F$ is the corresponding enhanced output. DNN--F applies a fusion model to enhance the noisy features from all channels concurrently, whereas DNN--S exploits $P$-channel-dependent models.

\subsection{DNN--C enhancement system}
\begin{figure}[!t]
{\centering \includegraphics[width=0.9\columnwidth]{./ dnn_c.eps}
\caption{DNN--C block diagram depicted in (a) contains both DNN--FC and DNN--DP functions. The training target for DNN--FC is the clean LPS $\mathbf{T}$, while that for DNN--DP is the ground truth voice $\mathbf{G}_p^C$. In addition, the DNN--DP model is determined first in (b), and then fixed for performing DNN--FC in (a).}\label{fig:dnnc}}
\end{figure}

The third proposed system, DNN--C, consists of two stages, the distributed processing stage and the fusion stage. Both stages use DNN models and are therefore represented by ``DNN--DP'' and ``DNN--FC'', respectively. The general diagram of DNN--C is shown in Fig. \ref{fig:dnnc} (a), and the detailed configuration of the DNN--DP stage is shown in Fig. \ref{fig:dnnc} (b). According to Fig. \ref{fig:dnnc} (b), $P$-self-powered devices, each for every individual microphone channel, are employed and a channel-specific DNN model denoted by DNN--DP$_p$ is conducted on each of these devices for suppressing noise from the input $\mathbf{Y}_p$ to produce the enhanced features, denoted by $\mathbf{E}_p^C$. Like DNN--S and DNN--F models, each DNN--DP$_p$ is composed of $L$ layers with ReLU activation and linear transformation functions, and can be formulated as follows:
\begin{equation}\label{eq:dnndp}
\mathbf{E}_p^C=\mbox{\textit{DNN-DP}}_p\{\mathbf{Y}_p\}.
\end{equation}
Therefore, the $P$--DNN models---DNN--DP$_1$, DNN--DP$_2$,$\cdots$, DNN--DP$_P$---are first estimated in the DNN--DP stage. In particular, the target feature for each of the $P$--DNN models, denoted by $\mathbf{G}_p^C$, is created by the channel-dependent clean speech, which is recorded from the $p$th microphone in the noise-free environment. That is, the $\mathbf{Y}_p$--$\mathbf{G}_p^C$ speech pair is used to train the associated DNN--DP$_p$ model.

As for the second stage, ``DNN--FC'', the $P$--DNN outputs at the first stage are concatenated to form a new feature $\mathbf{E}^{C'}$, i.e., $\mathbf{E}^{C'}=[\mathbf{E}_1^{C};\cdots;\mathbf{E}_p^{C};\cdots;\mathbf{E}_P^{C}]$, that is used together with the clean target $\mathbf{T}$ to train the DNN--FC model, which can be formulated by:
\begin{equation}\label{eq:dnnfc}
\mathbf{E}^C=\mbox{\textit{DNN-FC}}\{\mathbf{E}^{C'}\}.
\end{equation}
Therefore, the overall operation of the DNN--C system can be represented by
\begin{equation}\label{eq:dnnc}
\mathbf{E}^C=\mbox{\textit{DNN-FC}}\{\textit{DNN-DP}_1\{\mathbf{Y}_1\},\cdots, \textit{DNN-DP}_P\{\mathbf{Y}_P\}\}.
\end{equation}

\subsection{Phase component}
Briefly, the MC SE systems presented here enhance multiple signal sources in the frequency domain. Therefore, an ISTFT is applied to the updated spectrogram to produce the enhanced time-domain signal. It is worth mentioning that the phase component in the enhanced spectrogram is generated from a noisy speech through STFT from a specific channel. For DNN--S, the enhanced signal processor shown in Fig. \ref{fig:ovall} is performed individually on each channel for reconstructing the speech waveform with the enhanced LPS and the preserved noisy phase. Conversely, the phase used for DNN--F and DNN--C is extracted from one of the noisy channels, where it is assumed that the recorded sound has the highest SNR and optimal speech quality among all channels.

\section{Experiments}\label{sec:experi}
In the following subsections, we first describe the experimental setup of the SE distributed-microphone tasks and then present the experimental results together with some discussions for the presented systems.

\subsection{Experimental setup}
\begin{figure}[!t]
{\centering \includegraphics[scale=0.9]{./ config.eps}
\caption{MC system consists of seven microphones (``m1'', ``m2'',$\cdots$, ``m7''); microphones m1, m2,$\cdots$, m6 were placed around the speaker T with a radius of 0.5 meter, and m7 was put behind m1 and oriented towards T at a distance of 1 meter.}\label{fig:conf}}
\end{figure}
The layout of the distributed MC system is shown in Fig. \ref{fig:conf}. Seven microphones ($P=7$) of the same brand and model (Sanlux HMT-11) were used and denoted by ``m1'', ``m2'',$\cdots$, ``m7'', respectively, while the target speaker was represented by ``$T$''. In this system, six microphones, m1, m2,$\cdots$, m6, were placed around the speaker and equally spaced by an angle $\pi/3$ in the median plane with a 0.5-meter radius. The m7 microphone was placed right behind the m1 microphone and oriented towards the $T$ speaker at 1 meter. Notably, recording speech signals in such a distributed MC configuration might introduce different amounts of time delays across microphone channels \cite{redondi2009geometric}. However, here we assume that the difference in time delay across channels is small and can be negligible for simplicity in analysis.

For the evaluation task, we used the Taiwan Mandarin Hearing in Noise Test (TMHINT) \cite{huang2005development} to prepare the speech dataset. According to the script provided by the TMHINT dataset, 300 phrases were selected as the training set, while the remaining 20 utterances were used for testing. Phrases in the training set were individually pronounced by a male and a female in a noise-free environment at a sampling rate of 16 kHz. These recordings were then corrupted by eight different types of noise (cockpit, machine gun, alarm, cough, PC fan, pink, babble, and engine) at eight different noise levels (ranging from $-5$ to $16$ dB SNRs with a 3 dB interval). Thus, 38,400 utterances ($300\times 2\times 8\times 8$) were reproduced and then received by each of the seven microphones. On the other hand, the test utterances were first recorded by another male and female speaker and then contaminated with siren and street noises at -5, 0 and 5 dB SNRs. Therefore, there were 240 noisy utterances ($20\times 2\times 2\times 3$) transmitted to each of the seven microphones.

For any of the seven microphone channels, each received utterance was first split into overlapped frames with a 32-ms frame duration and a 16-ms frame shift, and a series of 257-dim frame-wise LPS were constructed accordingly. The context feature for each frame was then created by concatenating the LPS of three neighboring frames. Therefore, the dimensions of the input feature were 771 ($257\times 3$) for each DNN--S$_p$, 771 ($257\times 3$) for DNN--DP$_p$, 5397 ($257\times 3\times 7$) for DNN--F, and 1799 ($257\times 7$) for DNN--FC, respectively. By contrast, the output dimensions of DNN--S$_p$, DNN--F, DNN--DP$_p$, and DNN--FC were 257. In addition, each DNN--S$_p$ and DNN--F model consisted of seven layers and 2,048 nodes per hidden layer. A DNN--DP$_p$ model was arranged to have five layers, whereas the DNN--FC model had four layers. The number of nodes for each hidden layer of DNN--DP and DNN--FC models was set to 2,048.

Three metrics were used to evaluate the enhanced utterances, including perceptual evaluation of speech quality (PESQ) \cite{rix2001perceptual}, short-time objective intelligibility (STOI) \cite{taal2011algorithm}, and segmental SNR improvement (SSNRI) \cite{mermelstein1979evaluation,plapous2006improved}. Higher scores for PESQ, STOI, and SSNRI indicate better enhanced performance.

\subsection{Experimental results}
The averaged STOI scores of the test utterances (input) and enhanced utterances (output) of individual microphone channels of the DNN--S system are listed in Tables \ref{tab:stoinoidnns}. It is clear that almost all DNN--S$_p$ models could improve the STOI score, except for channel m1 in DNN--S$_1$. One possible reason is the overfitting issue that may have degraded the generalization capability of the model in the testing environments. Meanwhile, the STOI scores of noisy testing utterances varied with different recording channels, which may be owing to the different locations of the microphones deployed in a space, as shown in Fig. \ref{fig:conf}, despite the fact that these microphones were at the same distance from the speaker. In addition, Table \ref{tab:sdinoidnns} lists the SSNRI scores achieved by individual DNN--S$_p$ models. Similar to the case in Table \ref{tab:stoinoidnns}, all the channel-wise DNN--S models bring significant SNR improvements, whereas DNN--S$_1$ (m1) was less effective than the other models (channels).

\begin{table}[!b]
\begin{center}
\caption{Average STOI scores of the noisy utterances tests (input) and enhanced utterances (output) of the individual DNN--S channels.}\label{tab:stoinoidnns}
\begin{tabularx}{\columnwidth}{>{\centering}m{0.9cm}|>{\centering}m{0.61cm}>{\centering}m{0.61cm}>{\centering}m{0.61cm}>{\centering}m{0.61cm}>{\centering}m{0.61cm}>{\centering}m{0.61cm}>{\centering\arraybackslash}X}
\toprule
\hline
\textbf{STOI}&\textbf{m1}&\textbf{m2}&\textbf{m3}&\textbf{m4}&\textbf{m5}&\textbf{m6}&\textbf{m7}\\
\hline
\textbf{Noisy}&\textbf{0.672}&0.581&0.631&0.668&0.671&0.666&0.663\\
{\scriptsize\textbf{DNN--S}}&0.670&\textbf{0.677}&\textbf{0.668}&\textbf{0.685}&\textbf{0.695}&\textbf{0.679}&\textbf{0.673}\\
\hline
\bottomrule
\end{tabularx}\vspace{-0.1cm}
\end{center}
\end{table}
\begin{table}[!b]
\begin{center}
\caption{Average SSNRI of enhanced utterances of individual channels of DNN--S.}\label{tab:sdinoidnns}
\begin{tabularx}{\columnwidth}{>{\centering}m{0.9cm}|>{\centering}m{0.61cm}>{\centering}m{0.61cm}>{\centering}m{0.61cm}>{\centering}m{0.61cm}>{\centering}m{0.61cm}>{\centering}m{0.61cm}>{\centering\arraybackslash}X}
\toprule
\hline
\textbf{SSNRI}&\textbf{m1}&\textbf{m2}&\textbf{m3}&\textbf{m4}&\textbf{m5}&\textbf{m6}&\textbf{m7}\\
\hline
\textbf{Noisy}&0.000&0.000&0.000&0.000&0.000&0.000&0.000\\
{\scriptsize\textbf{DNN--S}}&\textbf{8.722}&\textbf{13.179}&\textbf{12.859}&\textbf{12.132}&\textbf{8.699}&\textbf{12.530}&\textbf{12.127}\\
\hline
\bottomrule
\end{tabularx}\vspace{-0.1cm}
\end{center}
\end{table}
\begin{table}[!b]
\begin{center}
\caption{Average STOI and SSNRI scores of the enhanced utterances of DNN--S, DNN--F, and DNN--C under all noise conditions and microphone channels.}\label{tab:overall}
\begin{tabularx}{\columnwidth}{>{\centering}m{1cm}|>{\centering}m{1.9cm}>{\centering}m{1.9cm}>{\centering\arraybackslash}X}
\toprule
\hline
&\textbf{DNN--S}&\textbf{DNN--F}&\textbf{DNN--C}\\
\hline
\textbf{STOI}&0.678&0.764&\textbf{0.770}\\
\textbf{SSNRI}&11.464&12.760&\textbf{16.729}\\
\hline
\bottomrule
\end{tabularx}\vspace{-0.1cm}
\end{center}
\end{table}

In Table \ref{tab:overall}, the averaged STOI and SSNRI scores of the enhanced utterances of DNN--S, DNN--F, and DNN--C under all noise conditions and microphone channels are shown. For DNN--S, here we report the STOI and SSNRI scores averaged over all channels shown in Tables \ref{tab:stoinoidnns} and \ref{tab:sdinoidnns}. From these tables, it is clear that both DNN--F and DNN--C achieve superior evaluation scores than DNN--S, and these results further confirm that the two MC SE systems outperform DNN--S, a single-channel SE system, in improving the intelligibility and SNR of noisy utterances. Furthermore, both evaluation metrics indicate that DNN--C outperforms DNN--F, revealing the superiority of the two-stage SE architecture in promoting intelligibility and noise reduction for distorted signals.

Figure \ref{fig:detail} compares DNN--F and DNN--C under siren and street noise conditions under all SNR levels in terms of the averaged (a) STOI and (b) SSNRI evaluation metrics. From this figure, we further confirm that DNN--C provides higher intelligibility and more significant SNR improvements than DNN--F.

\begin{figure}[!t]
{\centering \includegraphics[width=0.9\columnwidth]{./ nostpe.eps}
\caption{The (a) STOI and (b) SSNRI scores of DNN--F and DNN--C enhanced noisy utterances in siren and street noise environmentss, with an average of more than three SNR levels.}\label{fig:detail}}
\end{figure}

Figures \ref{fig:specana}(a)-(d) show the spectrograms of an utterance under four conditions: (a) clean noise-free, (b) noise-corrupted (with a PESQ score of 1.642), (c) noise-corrupted and then enhanced by DNN--F (with a PESQ score of 1.834), and (d) noise-corrupted and then enhanced by DNN--C (with a PESQ score of 1.923). The utterance was corrupted with street noise at -5 dB SNR, and the noisy utterance was recorded by the m1 microphone. From these figures, we find that the spectrogram of the DNN--C-enhanced utterance in Fig. \ref{fig:specana}(d) is quite similar to that of the clean utterance in Fig. \ref{fig:specana}(a). In addition, by comparing Fig. \ref{fig:specana}(d) to Fig. \ref{fig:specana}(c), DNN--C reveals clearer spectral characteristics and sound structures than DNN--F, as indicated by the red blocks. These observations also explain why DNN--C achieved higher evaluation scores than DNN--F.

\begin{figure}[!t]
{\centering \includegraphics[width=0.9\columnwidth]{./ spec.eps}
\caption{Spectrograms of (a) a clean utterance, (b) the noisy signal recorded by m1, (c) the DNN--F enhanced speech, and (d) the DNN--C enhanced version.}\label{fig:specana}}
\end{figure}

\section{Conclusion}\label{sec:conclu}
In this study, we presented three DNN--based SE systems under a distributed-microphone scenario applied in the diffuse-noise field environment. These three systems (DNN--S, DNN--F, and DNN--C) were evaluated in the TMHINT dataset. Experimental results showed that all of these systems were able to reduce the noise effect to improve speech intelligibility. Meanwhile, the two-stage DNN--C system achieved the optimal objective intelligibility score among the three systems. In the future, we plan to perform the DNN-based distributed-microphone SE task by properly selecting the recording channels rather than using them all as the system input. Furthermore, in addition to diffuse-noise environments, the DNN-C algorithm will be applied to reverberant environments to examine its capability in reducing the respective distortion.

\vfill\pagebreak

\bibliographystyle{ieeetr}


\section{Background and Motivation}

The mathematical definitions of Kullback-Leibler (KL) Divergence and Cross Entropy do not imply an upper bound.
Consider a simple alphabet $\mathbb{Z}=\{\mathbf{z}_1, \mathbf{z}_2\}$, and two probability mass functions (PMFs) defined upon $\mathbb{Z}$: $P=\{0+\epsilon,1-\epsilon\}$ and $Q=\{1+\epsilon,0+\epsilon\}$, where $0 < \epsilon < 1$.
When $\epsilon \rightarrow 0$, we have $\mathcal{D}_{\text{KL}}(P||Q) \rightarrow \infty$ and $\mathcal{H}(P, Q) \rightarrow \infty$.

Chen and Golan proposed an information-theoretic measure for analyzing the cost-benefit of machine- and human-centrics in data intelligence workflows \cite{Chen:2016:TVCG}.
Given a process $\mathbf{P}_i$ with $\mathbb{Z}_i$ as its input alphabet and $\mathbb{Z}_{i+1}$ as its output alphabet, the cost benefit ratio of $\mathbf{P}_i$ is as follows:
\begin{equation}\label{eq:CBR}
\begin{split}
    \frac{\text{Benefit}}{\text{Cost}}
    &= \frac{\text{Alphabet Compression}-\text{Potential Distortion}}{\text{Cost}}\\
    &= \frac{\mathcal{H}(\mathbb{Z}_i) - \mathcal{H} (\mathbb{Z}_{i+1}) + \mathcal{D}_{\text{KL}}(\mathbb{Z}'_i||\mathbb{Z}_i)}{\text{Cost}}
\end{split}
\end{equation}
\noindent where $\mathbb{Z}'_i$ is an alphabet reconstructed based on $\mathbb{Z}_{i+1}$ by a reverse process $\mathbf{P}_i^{-1}$.
It has the same letters as $\mathbb{Z}_i$ but is likely to have a different PMF since in data intelligence, $\mathbf{P}_i$ will bring about information loss.

Because the alphabet compression part, $\mathcal{H}(\mathbb{Z}_i) - \mathcal{H} (\mathbb{Z}_{i+1})$, is bounded when $||\mathbb{Z}_i||$ and $||\mathbb{Z}_{i+1}||$ are finite but the potential distortion part, $\mathcal{D}_{\text{KL}}(\mathbb{Z}'_{i}||\mathbb{Z}_{i})$, is unbounded, it is difficult to interpret the meaning of the benefit when it approaches $-\infty$.

Chen and Golan discussed this issue when they were working on Equation \ref{eq:CBR}, and agreed that the issue could be addressed by using the conventional method employed in many practical applications of information theory, i.e., by fixing an $\epsilon$ as the lower bound of probability values in a PMF, $\forall p \in P$, $\epsilon \le p \le (1-\epsilon)$.

In the analysis of the results of an empirical study, Kijmongkolchai et al. converted the values of accuracy and response time to benefit and cost in Equation \ref{eq:CBR}.
They set the upper bound as $\max \bigl( \mathcal{D}_{\text{KL}}(P||Q) \bigr) = 2\mathcal{H}_\text{max}(P)$, such that the benefit is bounded by $[-\mathcal{H}_\text{max}(P),+\mathcal{H}_\text{max}(P)]$.
However, this raises a number of questions, such as:
\begin{enumerate}
    \item Should there be an upper bound? (See Section \ref{sec:Upper}.)
    \item If (1) is true, should $\mathcal{D}_{\text{KL}}(\mathbb{Z}'_i||\mathbb{Z}_i)$ be replaced with a bounded measure? (see Section \ref{sec:New}.)
    \item If (1) is true, what is the most suitable bounded divergence? (see Section \ref{sec:Analysis}.)
\end{enumerate}

\section{The Existence of an Upper Bound}
\label{sec:Upper}
During his research visit to Oxford in July 2017, Sbert suggested to examine the upper bound based on the semantics of coding associated with the cross entropy.
On 19 July 2017, they discovered that there is an upper bound of the value $||\mathbb{Z}||-1$ of $\mathcal{H}(P, Q)$.

Let $\mathbb{Z} = \{\mathbf{z}_1, \mathbf{z}_2, \ldots, \mathbf{z}_n\}$ be an alphabet, which is associated with a PMF, $Q$, such that:
\begin{align*}
q(\mathbf{z}_n) &= \epsilon, \quad\text{(where $0 < \epsilon < 2^{-(n-1)}$}),\\
q(\mathbf{z}_{n-1}) &= (1-\epsilon)2^{-(n-1)},\\
q(\mathbf{z}_{n-2}) &= (1-\epsilon)2^{-(n-2)},\\
&\cdots\\
q(\mathbf{z}_{2}) &= (1-\epsilon)2^{-2},\\
q(\mathbf{z}_{1}) &= (1-\epsilon)2^{-1} + (1-\epsilon)2^{-(n-1)}.
\end{align*}

When we encode this alphabet using an entropy binary coding scheme, we can be assured to achieve an optimal code with the lowest average length for codewords.
One example of such a code for the above probability is:

\begin{align*}
\mathbf{z}_1&: 0\\
\mathbf{z}_2&: 10\\
\mathbf{z}_3&: 110\\
&\cdots\\
\mathbf{z}_{n-1}&: 111\ldots10 \quad\text{(with $n-2$ ``1''s and one ``0'') }\\
\mathbf{z}_n&: 111\ldots11 \quad\text{(with $n-1$ ``1''s and no ``0'') }
\end{align*}

In this way, $\mathbf{z}_n$, which has the smallest probability, will always be associated with a codeword with the maximal length of $n-1$.
Regardless whatever a PMF is defined upon $\mathbb{Z}=\{\mathbf{z}_1, \mathbf{z}_2, \ldots, \mathbf{z}_n\}$,
there is no need to code any letter $z_i \in \mathbb{Z}$ with more than $n-1$ bits.
The entropy coding scheme is designed to minimize the number of bits to be transmitted over a communication channel for sending a ``very long'' sequence of letters in the alphabet.
The phrase ``very long'' implies that the string exhibits the PMF used in the coding.
However, for any other string that does not exhibit the PMF used in the coding, there will be inefficiency.
This inefficiency is usually measured using cross entropy.
Let $P$ be the PMF of such a string, and $P$ may be different from $Q$.
The inefficiency is measured by:
\[
    \mathcal{H}(P, Q) = -\sum_{i=1}^n p_i \log_2 q_i
    = \mathcal{H}(P) + \mathcal{D}_{\text{KL}}(P||Q)
\]

Clearly, the worst case is that the letter, which was encoded using the most number of bits, $n-1$, turns out to be the most frequently used letter.
It is so frequent that all letters in the string are of this letter.
The average codeword length of this string is thus of $n-1$ bits.
Since there is no informative variation in the PMF $P$ for this very long string, i.e., $\mathcal{H}(P) = 0$, in principle, the transmission of this string (of $n-1$ bits per letter) is unnecessary.
The situation cannot be worse.
Therefore $n-1=||\mathbb{Z}||-1$ is an upper bound $\top_{\text{CE}}$ for the cross entropy.

An upper bound, $\top_{\text{KL}}$, of the KL-Divergence between $P$ and $Q$, can be derived from:
\begin{equation}\label{eq:CEandKL}
	\mathcal{D}_{\text{KL}}(P||Q) = \mathcal{H}(P, Q) - \mathcal{H}(P) \le \top_{\text{CE}} - \min_{\forall P}\bigl( \mathcal{H}(P) \bigr)
\end{equation}
\noindent In the cases where all PMFs are possible, the minimal Shannon entropy is 0.
Hence we have $\top_{\text{KL}} = \top_{CE}$.

There is a special case worth mentioning.
In practice, it is common to assume that $Q$ is a uniform PMF, i.e., $q_i = 1/n, \forall q_i \in Q$, typically because $Q$ is unknown or varies frequently.
Hence the assumption leads to a code with an average length equaling the maximum entropy $\mathcal{H}_{\text{max}} = \log_2 (n)$.
In practice, the actual code length would be $\lceil \log_2 (n) \rceil$.

Under this special (but rather common) condition, all letters in a very long string have codewords of the same length.
The worst case is that all letters in the string turn out to the same letter.
Since there is no informative variation in the PMF $P$ for this very long string, i.e., $\mathcal{H}(P) = 0$, in principle, the transmission of this string is unnecessary.
The maximal amount of inefficiency is thus $\log_2 (n)$ or in practice $\lceil \log_2 (n) \rceil$.
This is indeed much lower than the upper bound $\top_{\text{CE}} = n-1$, justifying the assumption or use of a uniform $Q$ in many situations.

\section{A New Bounded Divergence}
\label{sec:New}
During his research visit to Madrid in July 2019, Chen discussed a possibly new measure with Sbert.
They conducted a literature study, and did not find anything similar in the literature.

Given two PMFs, $P$ and $Q$, associated with the same alphabet $\mathbb{Z}$, the new measure is:
\begin{equation} \label{eq:New}
    \mathcal{D}_{\text{new}}(P||Q) = \sum_{i=1}^n p_i \log_2 \bigl( |p_i - q_i| + 1 \bigr)
\end{equation}

\noindent
Like $\mathcal{D}_{\text{KL}}(P||Q)$, $\mathcal{D}_\text{new}(P||Q)$ is not commutative.
The cost-benefit ratio in Equation (\ref{eq:CBR}) does not require it to be commutative because it is about a reconstruction process that maps its output alphabet to its input alphabet but not vice versa.
Importantly, $\mathcal{D}_{\text{new}}(P||Q)$ is bounded.
Obviously, $\mathcal{D}_{\text{new}}(P||Q) \ge 0$, and it equals 0 if and only if $P \equiv Q$.
Because $|p_i - q_i| < 1$, we have
\[
    \mathcal{D}_{\text{new}}(P||Q) = \sum_{i=1}^n  p_i \log_2 \bigl( |p_i - q_i| + 1 \bigr)
    \le \sum_{i=1}^n p_i \log_2 (2) = 1
\]

If we replace $\mathcal{D}_{\text{KL}}(\mathbb{Z}'_i||\mathbb{Z}_i)$ with $\mathcal{D}_{\text{new}}(\mathbb{Z}'_i||\mathbb{Z}_i)$, it is necessary to scale its value with $\mathcal{H}_{\text{max}}(\mathbb{Z}_i)$.
In other words, Equation (\ref{eq:CBR}) can be rewritten as:
\begin{equation}\label{eq:NewCBR}
\begin{split}
    \frac{\text{Benefit}}{\text{Cost}}
    &= \frac{\text{Alphabet Compression}-\text{Potential Distortion}}{\text{Cost}}\\
    &= \frac{\mathcal{H}(\mathbb{Z}_i) - \mathcal{H}(\mathbb{Z}_{i+1}) + \mathcal{H}_{\text{max}} \mathcal{D}_{\text{new}}(\mathbb{Z}'_i||\mathbb{Z}_i)}{\text{Cost}}
\end{split}
\end{equation}

 This new version can be used to measure the benefit with values that can be interpreted intuitively.
 For example, consider a type of things (or situations) that have two states, good or bad.
 The actual PMF is $\{ 0.7, 0.3 \}$, that is, 70\% good and 30\% bad.
 A process $\mathbf{P}_i$ conveys a ``simplified'' binary message that these things are always good, i.e., with a PMF $\{ 1, 0 \}$.
 %
 \begin{itemize}
     \item If the audience is totally misled to think that the things are always good, we have:
     $\text{Benefit} =  0.88 - 0 - 0.38 = 0.5$ bits.
     \item If half of the audience have the knowledge about the ground truth, the reconstructed PMF becomes $(0.85, 0.15)$, the benefit of $\mathbf{P}_i$ improves from 0.5 to: $\text{Benefit} =  0.88 - 0 - 0.20 = 0.68$ bits.
     \item If all of the audience have the knowledge about the ground truth, the information loss due to the binary message does not cause any potential distortion. The benefit is thus: $\text{Benefit} =  0.88 - 0 - 0 = 0.88$ bits.
 \end{itemize}

 Consider another example, where a type of things (or situations) are always good, that is, for the input, the PMF $\{ 1, 0 \}$.
 A process $\mathbf{Q}_i$ conveys a misleading message that these things are always bad, i.e., for the output, the PMF is $\{0, 1\}$.
 %
 \begin{itemize}
     \item If the audience is totally misled to think that the things are always bad, we have:
     $\text{Benefit} =  0 - 0 - 1 = -1$ bits. Note that with the KL-divergence, the benefit would approach $-\infty$.
     \item If half of the audience notice a note of sarcasm in the message, and interpret that the things are actually good. The reconstructed PMF becomes $(0.5, 0.5)$ and the benefit is: $\text{Benefit} =  0 - 0 - 0.58 = -0.58$ bits.
     \item If all of the audience correctly interpret the message to be that the things are always good, the benefit is thus: $\text{Benefit} =  0 - 0 - 0 = 0$ bits.
 \end{itemize}

 A generalised version of $\mathcal{D}_{\text{new}}$ is
 %
\begin{equation} \label{eq:NewG}
    \mathcal{D}^k_{\text{newG}}(P||Q) = \sum_{i=1}^n p_i \log_2 \bigl( |p_i - q_i|^k + 1 \bigr)
\end{equation}

A commutative version of $\mathcal{D}^k_{\text{newG}}$ is
 %
\begin{equation} \label{eq:NewGC}
    \mathcal{D}^k_{\text{newGC}}(P||Q) = \frac{1}{2} \sum_{i=1}^n \bigl( p_i + q_i \bigr) \log_2 \bigl( |p_i - q_i|^k + 1 \bigr)
\end{equation}

\section{Comparing Several Bounded Measures}
\label{sec:Analysis}
During the meeting in Madrid in July 2019, Chen and Sbert also compared the measure with several other measures in the literature.
These include the Jensen-Shannon (JS) divergence and the Minkowaski distances.

\begin{figure}[t]
    \centering
    \includegraphics[width=.94\linewidth]{./Compare-NonZero.pdf}
    \caption{Comparing the measures of divergence in relation to an alphabet in the range [0.05, 0.7]. The curve segments of DKL and 0.3DKL in the range [0, 0.05] do not represent the actual shapes.}
    \label{fig:NonZero}
\end{figure}

Figure \ref{fig:NonZero} shows the measures of
\begin{itemize}
    \item the KL-divergence (as DKL),
    \item its scaled down measures (as 0.3DKL),
    \item the new divergence in Equation (\ref{eq:New}) (as New),
    \item the generalized and commutative version of the new divergence with $k=2$ (as NewGC),
    \item the JS-divergence (as DJS), and
    \item several Minkowaski distances with $k = 0.5, 1, 1.6, 2, 4, 256$.
\end{itemize}
The $x$-axis shows the value of $\epsilon$ in the range of [0.05, 0.7].
The two PMFs are set as $P=\{1-\epsilon, \epsilon\}$ and $Q=\{\epsilon, 1-\epsilon\}$.
Therefore, $P=Q$ when $\epsilon=0.5$, and their divergence increases when $\epsilon \rightarrow 0$ or $\epsilon \rightarrow 1$.

\begin{figure}[t]
    \centering
    \includegraphics[width=.94\linewidth]{./Compare-NearZero.pdf}
    \caption{Comparing the measures of divergence in relation to an alphabet in the range near zero. The curve segments of DKL and 0.3DKL in the range $[0, 0.1^{10}]$ do not represent the actual shapes. The ranges $[0, 0.1^{10}]$ and $[0.1, 0.5]$ are only for references as they do not use the same logarithmic scale as in the range $[0.1^{10}, 0.1]$.}
    \label{fig:NearZero}
\end{figure}

\clearpage

We can observe that the curve of DKL quickly moves about 1.0, which is the maximal entropy of $P$ and $Q$.
When we scale the values of DKL to the one-third of its values, we can observe that the curve of 0.3DKL is similar to those of DJS and NewGC.
Meanwhile, the Minkowaski distances do not seem to capture much features of DKL.
The curve of the basic version of the new divergence in Equation (\ref{eq:New}) seems to differ from those of 0.3DKL, DJS, and NewGC.

Figure \ref{fig:NearZero} shows the same set of measures in the range near zero, that is, $\epsilon$ varies from $0.1^{10}$ to $0.1$. The ranges $[0, 0.1^{10}]$ and $[0.1, 0.5]$ are there only for references as they do not have the same logarithmic scale as that in the range $[0.1^{10}, 0.1]$.
We can observe that in $[0.1^{10}, 0.1]$ the curve of 0.3DKL also rises quickly as DKL.
This confirms that simply scaling the KL-divergence is not an adequate solution.

The curves of New and NewGC converge earlier than that of DJS.
If the curve of 0.3DKL is used as a benchmark as in Figure \ref{fig:NonZero}, the curve of NewGC is closer to 0.3DKL than that of DJS.

From Figures \ref{fig:NonZero} and \ref{fig:NearZero}, we can see that DJS (the JS-divergence) and NewGC (the commutative version of the new divergence with $k=2$) are the better options as bounded measures to replace the KL-divergence in Equation (\ref{eq:CBR}).

\begin{figure}[t]
    \centering
    \includegraphics[width=0.95\linewidth]{./Compare-Element.pdf}
    \caption{Comparing the measures of divergence in relation to an individual letter $z_i \in \mathbb{Z}$.
    The $y$-axis shows only one term related to $p_i$ and $q_i$ but not the whole divergence measure.}
    \label{fig:Element}
\end{figure}

Figure \ref{fig:Element} compares the measures for an individual letter $z_i \in \mathbb{Z}$.
We set $p_i = 0.5$ while varying $q_i$ from 0 to 1.
The $x$-axis shows the difference $\delta = q_i - p_i = q_i - 0.5$.
The $y$-axis shows the measures returned by:
\begin{itemize}
    \item eDKL: $p_i \log_2 \frac{p_i}{q_i}$ --- an element of the KL-divergence;
    \item eNew: $p_i \log_2 (|p_i - q_i|+1)$ --- an element of the non-commutative version of the new divergence (i.e., $k=1$);
    \item eNewG: $p_i \log_2 (|p_i - q_i|^2+1)$ --- an element of the generalized version of the new divergence with $k=2$;
    \item eNewC: $\frac{1}{2} (p_i+q_i) \log_2 (|p_i - q_i|+1)$ --- an element of the commutative version of the new divergence;
    \item eNewGC: $\frac{1}{2} (p_i+q_i) \log_2 (|p_i - q_i|^2+1)$ --- an element of the generalized and commutative version of the new divergence with $k=2$;
    \item eDJS: $\frac{1}{2} \bigl( p_i \log_2 \frac{2 p_i}{p_i + q_i} + q_i \log_2 \frac{2 q_i}{p_i + q_i} \bigr)$ --- an element of the JS-divergence.
\end{itemize}

We can make the following observations:

\begin{enumerate}
    \item The trend towards infinity and the negative values of eDKL are not intuitive to interpret semantically.
    \item The asymmetric patterns shown by eNewC, eNewGC, and eDJS indicates that these measures are influenced not only by the difference or ratio between $p_i = 0.5$ and $q_i \in [0, 1]$, but also by the combined value of $p_i + q_i$. It is useful to note that the absolute value of eDKL, $|p_i \log_2 \frac{p_i}{q_i}|$, is asymmetric.
    \item The asymmetric patterns of eNewC and eNewGC are more intutive than that of eDJS since the right part of a curve (i.e., $\delta > 0$) is not expected to be lower than the left part (i.e., $\delta < 0$). For $q_a(z_a) = 0.5-\delta$ and $q_b(z_b) = 0.5+\delta$, the divergence at letter $z_b$ is expected to be no less than that at $z_a$ because the combined probability at $z_b$ is higher.
    \item the curve shapes of eDJS, eNewG, eNewGC correspond to that of eDKL better than eNew and eNewG.
\end{enumerate}

\section{Conclusions}
\label{sec:Conclusions}
In conclusion, given an alphabet $\mathbb{Z}$ with a finite number of letters and with a true PMF $P$ and an estimated PMF $Q$, the measure of their divergence should be bounded if one interprets the divergence based on the inefficiency caused by using $Q$ instead of $P$ in coding $\mathbb{Z}$.

This confirms that the unbounded term $\mathcal{D}_{\text{KL}}(\mathbb{Z}'_i || \mathbb{Z}_i)$ in Equation (\ref{eq:CBR}) should ideally be replaced with a bounded term. Semantically, the upper bound of the potential distortion should ideally be the maximum entropy of $\mathbb{Z}_i$, i.e., $\mathcal{H}_{\text{max}}(\mathbb{Z}_i)$. The lower bound should ideally be 0.
As the JS-divergence and different versions of the new divergence in Equations (\ref{eq:New}, \ref{eq:NewG}, and \ref{eq:NewGC}) are all bounded by [0, 1], they all meet the boundedness requirement with a scaling factor $\mathcal{H}_{\text{max}}(\mathbb{Z}_i)$.

If one prefers to preserve the curvature of the KL-divergence to some extent (e.g., based on the curve 0.3DKL in Figure \ref{fig:NonZero}), the JS-divergence and the generalized versions of the new divergence in Equations (\ref{eq:NewG},\ref{eq:NewGC}) (with $k=2$) are suitable candidates.

On the other hand, there is no fundamental reason to preserve the curvature of the KL-divergence.
The close-to-linear patterns shown by the basic version of the new divergence (i.e., Equation (\ref{eq:New})) in Figures \ref{fig:NonZero} and \ref{fig:NearZero} suggest that it may be easy to estimate or interpret such a measure mentally. Hence the new divergence in Equations (\ref{eq:NewG},\ref{eq:NewGC}) (with $k=1$) may have some advantages.

Meanwhile, if one prefers to preserve the non-commutative property of the KL-divergence to some extent, the non-commutative version of the new divergence in Equations (\ref{eq:NewG}) is a suitable candidate.
On the other hand, if one prefers a commutative divergence measure, the JS-divergence and the commutative version of the new divergence in Equation (\ref{eq:NewGC}) are suitable candidates.

In addition, for measuring the divergence of an individual letter $z_i \in \mathbb{Z}$, the KL-divergence, JS-divergence, and the new divergence all have the common informative term in the form of $x \log (y)$, where $x$ and $y$ are some quantities.
As discussed in conjunction with Figure \ref{fig:Element}, if it is sensible to assume that this term should correlate to $x$ positively when $y$ is fixed, the commutative version of the new divergences in Equation (\ref{eq:NewGC}) exhibits such a positive correlation, while the JS-divergence exhibits a negative correlation.

We appreciate that there will be different preferences in different applications, and it may take many years for different preferences to converge.

\bibliographystyle{abbrv}
\nocite{*}