\section*{Acknowledgements}
We thank Xiaofeng Luo and Arghya Chatterjee for fruitful discussions.
This work supported by the U.S. National Science Foundation grant 1945296 and the U.S. Department of Energy grant DE-SC0020651.

\section{\label{sec:intro}Introduction}

A central goal of relativistic heavy ion collision physics is to probe the nonperturbative regime
  of quantum chromodynamics (QCD) by exploring the phase diagram of partonic matter as a function
  of thermodynamic variables.
At very high temperatures achievable at the Large Hadron Collider or Relativistic Heavy Ion
  Collider (RHIC), the quark-gluon plasma (QGP) is formed, a phase in which colored partons
  are the dynamical degrees of freedom~\cite{Adams:2005dq,Arsene:2004fa,Back:2004je,Adcox:2004mh}.
Collisions at progressively lower energies form matter with lower temperature and higher
  baryochemical potential, eventually reaching conditions at which the system is in the hadronic
  gas phase, in which partons are confined.
Understanding the nature of the transition between QGP and hadronic matter has been a
  central goal of the Beam Energy Scan (BES) program at RHIC~\cite{Stephanov:2007fk,Bowman:2008kc}.
Of particular interest is a possible critical point (CP)~\cite{Aoki:2006we,Gupta:2011wh,Fodor:2004nz,deForcrand:2002hgr,Qin:2010nq,Xin:2014ela,Shi:2014zpa,Fischer:2014ata,Lu:2015naa,Zhang:2017icm,Bazavov:2017tot,Fu:2019hdw,Fischer:2018sdj}
  at the terminus of the line characterizing a first-order phase transition on the QCD phase diagram,
  sketched in figure~\ref{fig:phase}.

 \begin{figure}
    \centering
    \includegraphics[width=0.45\textwidth]{QCD_phases.png}
    \caption{The QCD phase diagram, showing various states of QCD matter as functions of baryon chemical potential ($\mu_B$) and temperature (T). The CP location (and existence) is at this time not confirmed. \cite{Caines:2017search}}
    \label{fig:phase}
\end{figure}

In an infinite system in equilibrium, measured moments of conserved quantities such as charge or baryon number
  are directly sensitive to the correlation length of the system.
A non-monotonic energy dependence of higher-order moments (e.g. the kurtosis) of the net baryon
  distribution may reveal the a QCD critical point~\cite{Stephanov:1998dy,Bzdak:2019pkr}.
The STAR Collaboration has measured~\cite{Adamczyk:2013dal} the first four moments of the net-proton
  (a stand-in for net-baryon) distribution
  in heavy ion collisions at collision energies $\sqrt{s_{NN}}=7.7-200$~GeV, recently reporting non-monotonic
  behavior at the $3\sigma$ level~\cite{Adam:2020unf}.
This may suggest that the discovery of one of the long-sought features of bulk QCD is within reach.

However, great care must be taken in this analysis of small effects, to avoid nontrivial autocorrelations
  associated with event selection.
In particular, the protons and antiprotons used to construct the net-proton distribution are measured
  in the rapidity range $|y|<0.5$.
A critical feature of the observation is that the interesting behavior only occurs for the most
  central collisions, where centrality is estimated by the  multiplicity of charged
  particles-- {\it except for protons and antiprotons}--
  with pseudorapidity $|\eta|<1$.
In other words, one looks at (possibly correlated) fluctuations in the multiplicities of protons and antiprotons
  as a function of charged particle (except proton and antiproton) multiplicity.
Exclusion of the protons and antiprotons from the centrality measure does not trivially remove autocorrelations
  (by which we mean measurement of an event-wise quantity while selecting events with a measure directly related to that quantity),
  since very few (anti)protons are ``direct." Almost all are decay products of higher-mass baryons (e.g. $\Delta^{++}$ or $N^*$)
  whose sibling decay products are pions or other charged particles which are counted in the centrality measure.

A recent study by Chatterjee, et al.~\cite{Chatterjee:2019fey} based on a hadronic transport model concludes that the effect of these
  autocorrelations on the higher-order moments analysis is small.
However, the model employed in the study does not itself include critical fluctuations, so it is unclear how
  general this theoretical conclusion is.
In any event, it would be useful to check the robustness of the experimental observation by using a centrality
  estimator based on particles not emitted at midrapidity.

The STAR experiment has recently commissioned the Event Plane Detector (EPD)~\cite{Adams:2019fpo}
  to upgrade its capabilities in the BES program.
The EPD provides highly-segmented coverage for charged particles emitted in a wide pseudorapidity
  range ($2.1<|\eta|<5.1$) and full azimuth.
An analysis that used the EPD for a centrality estimate would be free of potential autocorrelations
  discussed above.

Chatterjee and collaborators~\cite{Chatterjee:2019fey} report that the summed multiplicity in the EPD
  acceptance would be a poor estimator of event centrality, because participant and spectator
  contributions are anticorrelated.
In what follows, we extend the study of Chatterjee and demonstrate that a slightly more sophisticated
  treatment of the signal in the EPD region improves the resolution substantially.
We show that this improvement persists even when the model is filtered through a realistic
  simulation of the STAR EPD, and that the detector can provide a measure of centrality
  with only slightly lower resolution than that from the STAR Time Projection Chamber (TPC, with coverage in the pseudorapidity range $|\eta| < 1$), but with no autocorrelation.

The paper is structured as follows: In Section~\ref{sec:methods} we describe the construction
  of global centrality estimators from collision data or the output of transport calculations.
We then discuss the EPD and limitations imposed by finite spatial and particle number resolution.
We then present an EPD-based centrality estimator that uses simple weights for the different
  regions of the EPD.
In Section~\ref{sec:results}, we compare the resolution of various centrality estimators
  both from midrapidity and forward rapidity regions.
In Section~\ref{sec:summary} we summarize and discuss the relevance these findings have for future centrality selection in the forward region of STAR.

\section{\label{sec:methods}Methods}
The Ultra-Relativistic Quantum Molecular Dynamic (UrQMD) model is a microscopic transport model used to simulate relativistic heavy ion collisions in the same energy ranges as the BES data \cite{Bass:1998ca}.
 For this study, we used 55k events per center of mass energy for Au+Au collisions at $\sqrt{s_{NN}} = $ 7.7, 11.5, 14.5, and 100k events for 19.6  GeV.

\subsection{Calculation of global centrality estimators}
\label{sec:GlobalVariables}
\begin{figure}
    \centering
    \includegraphics[width=0.45\textwidth]{dNdEta.pdf}
    \caption{Charged particle pseudorapidity distributions for BES energies of $\sqrt{s_{NN}} = $ 19.6 (magenta triangle), 14.5 (blue triangle), 11.5 (red square), 7.7 (black circle) GeV.  There is no centrality selection for these distributions.  The EPD acceptance is shown as a green box.}
    \label{fig:Rapidity}
\end{figure}

Our aim is to explore the connection between experimental global observables, generically designated $X$, with the impact parameter of the collision. Naturally, this connection depends on the collision model. Glauber models~\cite{Miller:2007ri} have been successfully used to make associate particle multiplicities at midrapidity in high-energy collisions. However, at lower energies (relevant to the RHIC BES program), the spectator-participant paradigm becomes less justified. Furthermore, the longitudinal dynamics of baryon stopping, which drives the physics of the forward direction covered by upgrade detectors~\cite{Adams:2019fpo}, requires a dynamical transport model.

We use the UrQMD transport model~\cite{Bass:1998ca,Bleicher:1999xi}, as was used by Chatterjee~\cite{Chatterjee:2019fey}. For a given impact parameter $b$, the model begins with realistic non-smooth initial conditions and a Monte-Carlo approach; subsequent evolution of the collision is based on string dynamics and hadronic rescattering to produce final-state particles that would be measured in a detector.

In the mid-rapidity region, we compute two quantities calculated by Chatterjee for reference. The first is $X_{\rm RM1}$, defined as the charged-particle multiplicity in the pseudorapidity range $|\eta|<0.5$. The second is $X_{\rm RM3}$, the charged particle multiplicity with $|\eta|<1.0$, excluding protons and antiprotons. These centrality estimators have been used in several analyses at RHIC~\cite{Adamczyk:2015eqo,Adamczyk:2013hsi,Adamczyk:2013dal}.

\begin{figure}
    \centering
    \includegraphics[width=0.45\textwidth]{AllETrNmipvsb.pdf}
    \caption{The sum of all charged particle yields in the EPD acceptance ($X_{FWD}$) vs impact parameter $b$ for the four collision energies considered in this paper.}

    \label{fig:FwdAllParticle}
\end{figure}

Within the EPD acceptance, the natural analog to these measures is $X_{\rm FWD}$, the charged-particle multiplicity in the range $2.1<|\eta|<5.1$. For completeness and to compare to Chatterjee, we calculate this quantity as well, but two problems with $X_{\rm FWD}$ must be considered: it cannot be unambiguously measured in the EPD as the EPD is essentially a segmented calorimeter~\cite{Adams:2019fpo}; which we discuss this in more detail in section~\ref{sec:EPDsim}, and most importantly, the nontrivial dynamics of stopping drives yields in the forward direction. These dynamics depend on the collision energy as the spectator region ranges within the EPD acceptance (see Figure~\ref{fig:Rapidity}). They also depend on centrality, which varies the relative contributions of participants and spectators. This can be seen in Figure \ref{fig:FwdAllParticle}, which is the total charged particle yield in the EPD acceptance ($X_{FWD}$) versus the impact parameter for a variety of energies.  In this figure we see that the correlation between forward multiplicity and impact parameter decreases with decreasing collision energy due to the contribution from spectators and participants within the EPD acceptance. However, due to the segmentation of the EPD in $\eta$, it is possible to separate the different contributions, as can be seen in Figure \ref{fig:RingbyRing}.  In section~\ref{sec:RingWeights}, we present a novel approach to use this nontrivial behavior to craft an observable centrality estimator with strong correlation to impact parameter at forward rapidity.

\begin{figure}
    \centering
    \includegraphics[width=0.45\textwidth]{Rings7_7.png}
    \caption{The sum of charged particle yields for a given EPD ring using simulated UrQMD events at $\sqrt{s_{NN}}$ = 7.7 GeV. Ring 1 is the ring closest to the beam-pipe and Ring 16 is the closest to mid-rapdity.  The precise acceptance boundaries for each ring are listed in Table \ref{table:rings}.}
    \label{fig:RingbyRing}
\end{figure}

\subsection{\label{sec:EPDsim}EPD Simulation}

The STAR TPC, which covers mid-rapidity in the experiment, measures individual charged particles with 90\% efficiency~\cite{anderson2003star}, so a centrality measure based on multiplicity is appropriate. This is not the case for the STAR EPD, which consists of 744 tiles of 1.2-cm-thick scintillator~\cite{Adams:2019fpo}. When a relativistic charged particle passes through a tile, it deposits a small amount of energy, $dE_1$, probabilistically according to the Landau distribution $\rho(dE_1)$. This distribution is characterized by two parameters: the most probable value ($dE_{\rm MPV}$) for energy loss and a width ($dE_{\rm WID}$); both the energy scale and the relative width ($dE_{\rm WID}/dE_{\rm MPV}$) are determined by the thickness of the scintillator. For the EPD, $dE_{\rm WID}/dE_{\rm MPV}\approx0.2$~\cite{Adams:2019fpo}.

\begin{table}
\centering
\begin{tabular}{|c|c|c||c|c|c|}
 \hline
 & & & & & \\
 \hspace{0.1cm}EPD Ring\hspace{0.1cm} & $|\eta_l|$ & $|\eta_h|$ & \hspace{0.1cm}EPD Ring\hspace{0.1cm} & $|\eta_l|$ & $|\eta_h|$ \\ [1.5ex]
 \hline\hline
  & & & & & \\
 1 & \hspace{0.1cm}5.09\hspace{0.1cm} & \hspace{0.1cm}4.42\hspace{0.1cm} & 9 & \hspace{0.1cm}2.81\hspace{0.1cm} & \hspace{0.1cm}2.69\hspace{0.1cm} \\
 2 & \hspace{0.1cm}4.42\hspace{0.1cm} & \hspace{0.1cm}4.03\hspace{0.1cm} & 10 & \hspace{0.1cm}2.69\hspace{0.1cm} & \hspace{0.1cm}2.59\hspace{0.1cm} \\
 3 & \hspace{0.1cm}4.03\hspace{0.1cm} & \hspace{0.1cm}3.74\hspace{0.1cm} & 11 & \hspace{0.1cm}2.59\hspace{0.1cm} & \hspace{0.1cm}2.50\hspace{0.1cm} \\
 4 & \hspace{0.1cm}3.74\hspace{0.1cm} & \hspace{0.1cm}3.47\hspace{0.1cm} & 12 & \hspace{0.1cm}2.50\hspace{0.1cm} & \hspace{0.1cm}2.41\hspace{0.1cm} \\
 5 & \hspace{0.1cm}3.47\hspace{0.1cm} & \hspace{0.1cm}3.26\hspace{0.1cm} & 13 & \hspace{0.1cm}2.41\hspace{0.1cm} & \hspace{0.1cm}2.34\hspace{0.1cm} \\
 6 & \hspace{0.1cm}3.26\hspace{0.1cm} & \hspace{0.1cm}3.08\hspace{0.1cm} & 14 & \hspace{0.1cm}2.34\hspace{0.1cm} & \hspace{0.1cm}2.27\hspace{0.1cm} \\
 7 & \hspace{0.1cm}3.08\hspace{0.1cm} & \hspace{0.1cm}2.94\hspace{0.1cm} & 15 & \hspace{0.1cm}2.27\hspace{0.1cm} & \hspace{0.1cm}2.20\hspace{0.1cm} \\
 8 & \hspace{0.1cm}2.94\hspace{0.1cm} & \hspace{0.1cm}2.81\hspace{0.1cm} & 16 & \hspace{0.1cm}2.20\hspace{0.1cm} & \hspace{0.1cm}2.14\hspace{0.1cm} \\ [1ex]
 \hline
\end{tabular}
\caption{List of the low ($|\eta_l|$) and high ($|\eta_h|$) values for the $|\eta|$ ranges of each EPD ring \cite{Adams:2019fpo}.}
\label{table:rings}
\end{table}

As is the practice experimentally, we will quantify the signal in a tile by normalizing to the peak in the single-particle distribution:
 \begin{equation}
     \zeta \equiv \frac{dE}{dE_{\rm MPV}}  .
 \end{equation}
The single-particle distribution $\frac{dP_1}{d\zeta}$ then peaks at unity and has Landau width $\sim0.2$. Figure~\ref{fig:zSpectraLowFlux} shows the simulated distribution for a tile that is dominated by one or two particles crossing per collision event.

Depending on the collision energy and centrality, $N>1$ charged particles may pass through a tile in a given collision event, each leaving energy independently. The total energy in such events follows a distribution:
\begin{equation}
    \frac{dP_{N}}{d\zeta}\left(\zeta\right) = \int^{\infty}_{0}d\zeta_1 \frac{dP_1}{d\zeta}\left(\zeta_1\right)\frac{dP_{N-1}}{d\zeta}\left(\zeta-\zeta_1\right) .
\end{equation}
The total energy deposited into each tile is recorded; for every collision event, the distribution is a weighted sum of $\frac{dP_{N}}{d\zeta}$, seen for $N\leq5$ as the histogram in Figure~\ref{fig:zSpectraHighFlux} (distributions for $N=1,2,3,4,5$ are shown as shaded curves). This distribution has been shown to reproduce the EPD experimental measurements closely~\cite{Adams:2019fpo}.

In general, the more particles that pass through a tile in a collision event, the larger the measured $\zeta$ will be.
It can seem natural to construct EPD-based centrality measures analogous to those used in the TPC by substituting $\zeta$ for multiplicity. However, the most probable value for $\frac{dP_{N}}{d\zeta}$ is not found at $\zeta=N$, and there is clearly considerable overlap between $\frac{dP_{N}}{d\zeta}$ distributions. Indeed, for a low-flux tile (c.f. Figure~\ref{fig:zSpectraLowFlux}), it is most accurate to assume $N=1$ in any collision event for which $\zeta>0$. Assuming $N=\zeta$ (which is not even correct ``on average'' as the Landau distribution does not have a well-defined mean) in this case only builds in unwanted noise for an event-wise centrality measure. Even in the high-flux tile of Figure~\ref{fig:zSpectraHighFlux}, a $\zeta=4$ event is most likely caused by three particles having passed through the tile.

\begin{figure}
\includegraphics[width=0.45\textwidth]{LowMultPeak.png}
\caption{
A low-flux tile which is dominated by either one or two particles passing per collision events. The coloured peaks show the 1 (grey) and 2 (red) particle Landau distributions, whereas the spectrum (black circles) is a sum of the convoluted distributions (green line).
}
\label{fig:zSpectraLowFlux}
\end{figure}

\begin{figure}
\includegraphics[width=0.45\textwidth]{HighMultPeak.png}
\caption{
A high-flux tile which experiences a range of particle multiplicities, resulting in a spectrum (black circles) that is a sum of convoluted Landau distributions (green line). The coloured peaks show the 1 (grey), 2 (red), 3 (blue), 4 (purple) and 5 (dark grey) particle Landau distributions.}
\label{fig:zSpectraHighFlux}
\end{figure}

How these effects influence the correlation between a centrality estimator and the true impact parameter depends on the collision model and the detector geometry. Using the UrQMD model, the collisions are assumed to take place at the center of the STAR experiment (for simplicity) and, upon emission, charged particles are assumed to propagate in a straight line until they strike the EPD. A precise geometric model of the active elements of the EPD~\cite{Adams:2019fpo} is used to register the passage of charged particles through each tile. Each charged particle deposits energy according to the Landau distribution, as discussed above, with the net signal from a tile being the sum of all deposited energy.\footnote{Our simulation does not include the effects of the magnetic field in the STAR experiment; however, at these rapidities, this is relatively unimportant for measuring multiplicity as a function of $\eta$. Secondary scattering and production in surrounding material (magnet iron, beam pipe, etc) is a larger effect that is not accounted for, in our simulation. Exploring corrections due to these effects is outside the scope of the present study.}

Figure~\ref{fig:NmipRaw} shows the (anti-)correlation between the impact parameter from the model and the sum of the signals ($X_{\zeta} \equiv \sum_i\zeta_i$) from the 744 EPD tiles for 19.6~GeV Au+Au collisions. The ``noise'' effect of Landau fluctuations is clear in the extended tail of the distribution at large $X_{\zeta}$, reducing the correlation.

The effect of these fluctuations may be reduced by ``truncating'' the signal from each tile, replacing a tile's signal with:
\begin{equation}
      \zeta^\prime \equiv \begin{cases}
        \zeta, & \text{if~}\zeta < {\rm Mx} \\
    \textrm{Mx},              & \text{otherwise}
    \end{cases}
    \label{eq:zprime}
\end{equation}
We chose the value of ${\rm Mx} = 3$ for all energies and centralities for this paper, though this could be tuned based on an analysis of the most probable value of the number of particles that will pass through a given tile.  The result of applying this methodology can be seen in Figure~\ref{fig:TnMip}, which plots $X_{\zeta'} \equiv \sum_i\zeta^\prime_i$ versus impact parameter. The Pearson coefficient is about 0.98, as compared to only 0.39 for the correlation in Figure~\ref{fig:NmipRaw}.

\begin{figure}
    \centering
    \includegraphics[width=0.45\textwidth]{NmipRawvsb.pdf}
    \caption{$X_{\zeta}$ value versus impact parameter for $\sqrt{S_{NN}}$ = 19.6 GeV UrQMD simulations. The extended tail of the distribution is due to Landau fluctuations.}
    \label{fig:NmipRaw}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=0.45\textwidth]{TrNmipvsb.pdf}
    \caption{$X_{\zeta'}$ value versus impact parameter for $\sqrt{s_{NN}}$ = 19.6 GeV UrQMD simulations. The truncation of $\zeta$ (from Equation \ref{eq:zprime}) reduces the noise from Landau fluctuations seen in the correlation between $X_{\zeta}$ and $b$ in Figure \ref{fig:NmipRaw}.}
    \label{fig:TnMip}
\end{figure}

\subsection{A new centrality estimator using ring weights}
\label{sec:RingWeights}
As we discuss in detail in section~\ref{sec:results}, the 16 rings of the EPD (corresponding to different $|\eta|$ ranges listed in Table \ref{table:rings}) can be affected quite differently as the impact parameter of the collision is varied. Indeed, depending on the collision energy, the signal in some rings may {\it increase} as $b$ is increased, while the signal in others may {\it decrease}. Thus their contributions to the simple sum $X_{\zeta'}$ discussed may partly cancel, reducing the sensitivity of this measure to collision centrality.

The differential response can be accounted for-- indeed, even be exploited to increase sensitivity-- by constructing a new simple measure that weights each ring's ``contribution'' differently.
Below, we consider two ways to quantify this contribution.
\subsubsection{Weighted sum of $\zeta'$ \label{sec:WeightsZeta}}
First, we define a ring's contribution to be the sum of the truncated signals in each tile in the ring:
\begin{equation}
\label{eq:RingContributionZetaPrime}
    C_r \equiv \sum_{\substack{\text{tile }j\text{ in}\\ \text{ring } r}}\zeta^\prime_j\qquad\text{(energy-loss based)}
\end{equation}
where we have indicated explicitly that the ring's contribution is based on energy loss (in the next section, we will discuss an analogous, particle based approach). Since we consider symmetric collisions occurring mid-way between the EPD wheels, ring $r$ on the wheel positioned at $z=-375$~cm is summed with ring $r$ on the wheel at $z=+375$~cm.

Based on these contributions, we define our centrality measure as a weighted sum:
\begin{equation}
\label{eq:Xzeta}
    X_{W,\zeta^\prime} \equiv \sum^{16}_{r=1} W_r C_r + W_{17}
\end{equation}
where $W_i$ are parameters determined below.

We wish to maximize the correlation between $X_{\rm \zeta'}$ and some global quantity $G$ (for the moment, $G$ is the impact parameter, $b$, but we generalize the discussion in section~\ref{sec:summary}). A figure of merit may be the squared residual, summed over events:
\begin{equation}\label{eqn:chi2}
    \chi^2 = \sum^{N_{events}}_{j=1} (X_{W,\zeta^\prime,j} - G_j)^2  ,
\end{equation}
where $X_{W,\zeta^\prime,j}$ and $G_j$ are respectively the values of the estimator and
  global quantity (e.g. impact parameter) for event $j$.

Maximizing $\chi^2$ yields 17 linear equations:
\begin{equation}
\sum_{q=1}^{17} A_{q,t} W_q = B_t
\end{equation}
where:
\begin{align}
    A_{q,t}    &= \sum_{j=1}^{N_{\rm events}} C_{q,j}C_{t,j} & {\rm for~}q,t=1...16 &~\\
    A_{17,t}   &= \sum_{j=1}^{N_{\rm events}}C_{t,j}         & {\rm for~}t=1...16 &~\\
    A_{17,17}  &= N_{\rm events} &~&~\\
    B_{t}    &= \sum_{j=1}^{N_{events}}G_j C_{t,j} & {\rm for~}t=1...16 &~ \label{eq:Bt}\\
    B_{17}   &= \sum_{j=1}^{N_{events}}G_j &~ &~ \label{eq:B17}
\end{align}
Since $A$ is a symmetric, real, $17\times17$ matrix it can be easily inverted to find the unique best parameters $W_t$.
\subsubsection{Weighted sum of particles}
As discussed above, an EPD tile measures in a given event the energy deposited in the tile, and not the number of particles that actually passed through the tile. However, it is natural to ask whether the sensitivity would be improved if a ring's contribution would be the number of charged particles passing through tiles in that ring, i.e. if $C_r$ in equation~\ref{eq:RingContributionZetaPrime} would be redefined:
\begin{equation}
\label{eq:RingContributionParticle}
    C_r \equiv \sum_{\substack{\text{tile }j\text{ in}\\ \text{ring } r}}N_j\qquad\text{(particle based)}
\end{equation}
Here, $N_j$ is the number of particles that passed through tile $j$ in the event.

Analogous to $X_{W,\zeta^\prime}$ a weighted-sum centrality measure $X_{W,FWD}$ may then be constructed from the contributions of equation~\ref{eq:RingContributionParticle}, with the weights determined by equations~\ref{eqn:chi2}-\ref{eq:B17}.

\section{\label{sec:results}Results}

\begin{figure}
    \centering
    \includegraphics[width=0.45\textwidth]{Weights.png}
    \caption{The weights for both forward $\eta$ particles and $\zeta'$ when using the linear weight scheme from Section \ref{sec:RingWeights}, by EPD ring. The sign change in ring weights is motivated by spectator proton intrusion into the EPD's acceptance window, as can be seen in Figure \ref{fig:Rapidity}.}
    \label{fig:Weights}
\end{figure}

With the mathematical formalism described in Section \ref{sec:RingWeights}, we can examine the performance of the EPD as it pertains to relating signals within the detector to the impact parameter using UrQMD.

\subsection{\label{sec:level2}Correlations between yields and $b$}

In order to apply the formalism described in the previous section, we will compare the two mid-rapidity observables ($X_{RM1}$ and $X_{RM3}$) with four forward rapidity observables ($X_{\zeta^\prime}, X_{W,\zeta^\prime}, X_{FWD}$, and $X_{W,FWD}$), where only the observables designated by $W$ use the linear weight method formalized above. The weights determined by this method versus EPD ring number can be seen in Figure \ref{fig:Weights}, with the acceptance of the EPD rings detailed in Table \ref{table:rings}.  The sign of the weights changes when moving from a distribution that is dominated by spectators versus participants.  In Figure \ref{fig:Weights}, we can see that the EPD ring this occurs in changes as the collision energy changes due to the change in beam rapidity (see Figure \ref{fig:Rapidity}).  The weights for EPD distributions with $\zeta^\prime$ are very similar to those required by the raw particle counts.  These weights were then used to determine the correlation between $\{X_{W,FWD}, X_{W,\zeta^\prime}\}$ and $b$, the results of which can be seen in Figure \ref{fig:Variations}.

Figure \ref{fig:FwdAllParticle} is in agreement with the conclusions of the Chatterjee analysis \cite{Chatterjee:2019fey}, and indicates that a summation of the yield over the EPD acceptance is a poor observable. However, when we apply the linear weighting technique described in Section \ref{sec:RingWeights}, we recover a much more usable correlation between forward $\eta$ particle yields and $b$. This is summarised in Figure \ref{fig:Variations}.  For all the energies under consideration, $X_{RM3}$ (middle row) is well correlated with the impact parameter, whereas  $X_{\zeta'}$ from the EPD (bottom row) shows a decreasing correlation with decreasing energy.  However, the linear weighted $X_{W,\zeta'}$ (top row) shows that the correlation is restored.

\begin{figure}
    \centering
    \includegraphics[width=0.45\textwidth]{Variations.pdf}
    \caption{Correlation between the impact parameter, $b$, and three observables: $X_{W,\zeta'}$ (top), $X_{RM3}$ (middle), $X_{\zeta'}$ (bottom), for four different collision energies.
        }
    \label{fig:Variations}
\end{figure}

Experimentally, the centrality is not determined by a model dependent relationship between a global observable and the impact parameter, but rather by considering the global observable's distribution quantiles (though it should be noted that real world analyses include a Glauber model as the efficiency of recording an event only approaches 100\% for the most central collisions~\cite{Miller:2007ri}).  This indicates that, as long as a global observable is reasonably correlated with the impact parameter, the centrality distribution based on this selection criteria will also be reasonable.

\begin{figure}
    \centering
    \includegraphics[width=0.45\textwidth]{ImpactParameterSlices.pdf}
    \caption{The impact parameter distributions for centrality selections 0 - 5\%, 20 - 30\% and 90 - 100\%.  The black histograms are the $b$ impact parameter distributions if the centrality selection is determined from the impact parameter directly.  The green circles are determined using $X_{RM1}$ (closed) and $X_{RM3}$ (open), the triangles are determined using $X_{\zeta'}$ (open blue) and $X_{W,\zeta'}$ (closed red), and the squares are determined using $X_{FWD}$ (open blue) and $X_{W,FWD}$ (closed red).
    }
    \label{fig:ImpSelection}
\end{figure}

In Figure \ref{fig:ImpSelection}, the impact parameter distributions that are determined by using the appropriate quantiles for the global variables are shown along with distributions that result from directly using the quantiles of the impact parameter distribution.  For $\sqrt{s_{NN}}$= 19.6 GeV, all methods performed similarly (quantified in Section \ref{sec:resRes}).  The mid-rapidity observables, $X_{RM1}$ and $X_{RM3}$, have $b$ distributions which peak within the $b$ distribution slices and do not change drastically as the collision energy decreases from 19.6 GeV to 7.7 GeV.  The forward observables without any weighting, $X_{FWD}$ and $X_{\zeta'}$, have distributions which no longer lie under the $b$ distribution slices at the lower energies; this suggests a poor centrality resolution for these observables in the lower energy ranges of the BES program. This potential loss in resolution, however, can be compensated for by applying the weighting scheme discussed in section \ref{sec:methods}. We see that the distributions for $X_{W,FWD}$ and $X_{W,\zeta'}$ are under the $b$ distribution slices for all collision energies under consideration.

\begin{figure}
    \centering
    \includegraphics[width=0.45\textwidth]{CentralityResolution.pdf}
    \caption{Centrality resolution for collision energies $\sqrt{s_{NN}}$ = 19.6, 14.5, 11.5, and 7.7 GeV.  The observables based on mid-rapidity multiplicity are $X_{RM1}$ (green closed circles) and $X_{RM3}$ (dark green open circles).  The observables based on unweighted distributions from the forward region are $X_{FWD}$ (blue open square) for particle yield and $X_{\zeta'}$ (blue open triangle) for truncated energy loss in the EPD.  The observables based on linear weights are $X_{W,FWD}$ (closed dark red square) and $X_{W,\zeta'}$ (closed red triangle).
    }
    \label{fig:CentralityResolution}
\end{figure}

\subsection{\label{sec:resRes}Centrality Resolution}
From the $b$ distributions in Figure \ref{fig:ImpSelection}, we determined the centrality resolution for all $X$. As in \cite{Chatterjee:2019fey}, we employ the centrality resolution metric $\Phi$:

\begin{equation}
    \Phi_i \equiv \frac{\sigma_{b_{X,i}}^2}{\sigma_{b_i}^2}
\end{equation}
where $\sigma_{b_i}^2$ is the variance from the impact parameter distribution for a given centrality range $i$ when centrality is determined by $b$, and $\sigma_{b_{X,i}}^2$ is the variance from the impact parameter distribution for the same centrality range $i$ where centrality is determined by $X$.

Results for $\Phi$ can be seen in Figure \ref{fig:CentralityResolution}.  At the top RHIC energies, all centrality methods perform similarly.  The resolution at lowest energies is poorest for unweighted distributions in the forward region ($X_{FWD}$ for particle yield, and $X_{\zeta'}$ for truncated energy loss in the EPD), but this resolution is clearly recovered when we apply the linear weights detailed in Section \ref{sec:RingWeights} ($X_{W,FWD}$ and $X_{W,\zeta'}$, respectively). Both $X_{RM1}$ and $X_{RM3}$ are based on mid-rapidity multiplicities. The anomalous, upward point at the 5-10\% range is due to the smaller centrality bins used for our two most central selections, which increases the standard deviation of the distributions.

The poor performance of $X_{FWD}$ and $X_{\zeta'}$, which agrees with the conclusions in \cite{Chatterjee:2019fey}, is due to spectator proton intrusion into the EPD's acceptance. If we do not weight the EPD rings which are dominated by the spectator protons yield, which have positive correlation with $b$, with a different sign than the EPD rings where we are dominate by participants, which have a negative correlation, they will cancel out. This is the entire purpose of the methodology in section \ref{sec:RingWeights}; thus the correlation weights found for rings with spectators will have inverted signs compared with those rings with only participants (Figures \ref{fig:Rapidity} and \ref{fig:Weights}). The linear weight method recovers the centrality resolution lost by the simple sum method of $X_{FWD}$ and $X_{\zeta'}$

\section{\label{sec:summary}Summary and Discussion}
Observables in heavy-ion experiments can suffer from auto-correlation effects if the particles used to construct the observable are from the same acceptance as the event centrality is defined (e.g. $\kappa$ and $S$ for net-proton multiplicity being studied as part of the STAR BES program)~\cite{Luo:2013bmi}. For observables analysed at mid-rapidity, such as those measured using the TPC in the STAR experiment, it would thus be preferable to select centrality using a forward detector. The EPD is the prime candidate for a forward centrality selection detector at STAR, but at the lower collision energies of the BES there is significant spectator proton intrusion into the $\eta$ acceptance window of the EPD. It had been suggested that this spectator intrusion could potentially degrade the centrality resolution of the EPD due to the inverse correlation effect on particle yields with $b$ in those $\eta$ regions where the EPD and spectator protons coincide \cite{Chatterjee:2019fey}.  However we showed that by treating the particle yield correlations from the EPD rings individually, instead of simply summing the yield from the total EPD acceptance, we can account for correlations with both participants and spectators in an event. This treatment leverages the spectator protons in the EPD acceptance range as a relevant marker for global quantities (such as $b$).

In this paper, we outlined a method of applying a linear weight to the rings (Equation \ref{eq:Xzeta}) as one example of a procedure that improves the centrality resolution by properly weighting the contribution from each EPD ring based on minimizing the residual between a global observable and the EPD rings. Results from this method were shown for correlations with $b$, but the method is sufficiently general that weighting may be found for any global quantity $G$; for instance, the method could also be employed to weigh the EPD ring contributions as they correlate with $X_{RM3}$ or V.

More sophisticated methods than the simplest one we described in section~\ref{sec:RingWeights} are certainly
  possible and are under development.

Using this method we recovered the centrality resolution in the forward $\eta$ region that is lost in the lower energy ranges of STAR BES when simply considering a sum of the yields over the entire EPD $\eta$ range. Further, EPD simulation shows no marked degradation of centrality resolution when comparing centrality using the linear weighted sum of particles in the EPD acceptance range ($X_{W,Fwd}$) versus simulated energy deposition in the EPD itself ({$X_{W,\zeta'}$}). This strongly suggests the EPD can be used as a reliable centrality detector in STAR BES energy ranges of $\sqrt{s_{NN}}$ = 7.7, 11.5, 14.5, 19.6, and 200 GeV, which would greatly reduce the possibility of autocorrelations in analyses of observables at mid-rapidity.

We conclude by briefly considering the implication of this study in an experimental analysis,
  in which the true impact parameter is unknown.
Any estimate of the impact parameter resolution of any measurable estimator is then completely model-dependent;
  $X_{\zeta^\prime}$ may in fact be the best estimator in reality, despite being the worst in the
  UrQMD calculations (c.f. figure~\ref{fig:CentralityResolution}).
In such a case-- especially if even a years-long program of high-quality data results in a subtle wiggle in fluctuations
  at only the $3\sigma$ level~\cite{Adam:2020unf}-- nontrivial effects of autocorrelations must be ruled
  out in a model-independent way, not relying solely on transport calculations~\cite{Chatterjee:2019fey}.
The experimental procedure must be to
  (1) perform the analysis, using a common estimator (e.g. $X_{RM3}$) to select on centrality;
  (2) in a kinematic region far from the fluctuation measurement, construct a weighted estimator by following the procedure of
      section~\ref{sec:WeightsZeta}, where the common estimator is used as the global quantity $G$ in equations~\ref{eq:Bt}-\ref{eq:B17};
  (3) repeat the analysis, using the new estimator (e.g. $X_{W,\zeta^\prime}$) to select on centrality.
Persistence of the signal when using the new estimator would lead to greater confidence that  autocorrelations
  are not influencing the signal itself.
Such confidence would be most welcome for the subtlest signals of fundamental physics of QCD.


\section{Introduction}
A ring $R$ is called \textit{left Archimedean}, if for each nonunit element $r\in R$ we have $\bigcap_{n\in \mathbb{N}}r^nR=\{0\}$. The right Archimedean rings are defined in a similar way. In \cite{Sheldon}, it is shown that if $R$ is a domain and $\bigcap_{n\in \mathbb{N}} a^nR=0$, then the quotient field $Q(R[[x/a]])$ has infinite transcendent degree over the quotient field $Q(R[[x]])$. A very important result in this direction is  that if $R$ is a domain (or is completely integrally closed) and satisfies \emph{ACCP}, then $R$ is Archimedean, see \cite{andersonanti,Sheldon}. However, each \emph{ACCP}-domain is not necessarily Archimedean, see for example, \cite[p. 1127]{Dumitrescu}. There are similar results  in \cite{ArnoldII,condo}.

Ribenboim in \cite{Ribenboimsemisimple} defined the ring of
generalized power series $R[[S]]$ consisting of all maps from $S$ to
$R$ whose support is Artinian and narrow  with  the pointwise
addition and the convolution multiplication. This construction
provided interesting examples of rings, for example see \cite{Ribenboim (1995a),Ribenboim (1995b)}. In \cite{Ribenboim}, Ribenboim gave a sufficient condition for the ring $R[[S]]$ (not necessarily commutative) to be Noetherian.  However, Frohn  gave an example to show that \emph{ACCP} does not rise to the power series ring in general (see \cite{counter}).

In \cite{zim}, R. Mazurek and M. Ziembowski, introduced a ``twisted" version of the Ribenboim
construction and studied when it produces a von Neumann regular ring.
They also proved  that if $R$ is a domain and $\omega$ an endomorphism of  $R$, then  $R$ satisfies \emph{ACCPL}, $S$ is an \emph{ACCPL}-monoid and $\omega_s$  is injective  for each $s\in S$ if and only if $R[[S,\omega]]$ is an \emph{ACCPL}-domain.

As a similar result, for a strictly ordered monoid $S$ and a monoid homomorphism $\omega : S \rightarrow \End (R)$, in \cite{paykan} it is proved that if the attached skew generalized power series ring $R[[S,\omega]]$ is a right Archimedean domain then $R$ is a right Archimedean domain, $S$ is a right Artinian and narrow monoid and $\omega_s$ is an injective which preserves nonunits of $R$ for any $s \in S$.
Next step in this direction,
we show that when $R$ is right (resp. left) Archimedean and satisfies the \emph{ACC} on  annihilators, then $R[[x;\alpha]]$ is a right (resp. left) Archimedean ring.
In section \ref{archh}, we give an example to justify that the condition given in Theorem \ref{rigidar} is necessary on the ring $R$, to be able to pass the property of being Archimedean between $R$ and $R[[x;\alpha]]$.

A ring $R$ is \textit{reduced}, if the condition $a^2=0$ results in $a=0$. According to Krempa \cite{krem}, an endomorphism $\alpha$ of a ring $R$ is said to be \textit{rigid} if $a\alpha(a)=0$ implies $a=0$, for $a \in R$. A ring $R$ is said to be $\alpha$-\textit{rigid} if there exists a rigid endomorphism $\alpha$ of $R$. Clearly, every domain $D$ with a monomorphism $\alpha$ is rigid. It is clear that, $\alpha$-rigid rings are reduced, (see e.g., \cite{ore}).
The following Lemma  by Hong in \cite{ore} is very useful  in the following sections.

\begin{lemma}\label{jav}
Let $R$ be a $\alpha$-rigid ring and $a,b \in R$. Then
\begin{itemize}
\item[(1)] $R$ is reduced.
\item[(2)] If $ab=0$, then  we have
$a\alpha^{n}(b)=\alpha^{n}(a)b=0$  for any positive integer $n$.

\end{itemize}
\end{lemma}
\begin{proof}
(1) Let $a^2=0$. Then $a\alpha(a)\alpha(a)\alpha^2(a)=0$. So $a\alpha(a)\alpha(a\alpha(a))=0$. Because $R$ is rigid $a\alpha(a)=0$. Another use of it gives us $a=0$.

(2) Let $ab=0$. So $ba=0$. This implies that $a\alpha^n(ba)\alpha^{2n}(b)=0$. So $a\alpha^n(b)\alpha^n(a\alpha^n(b))=0$. Since $R$ is $\alpha-$rigid, then $a\alpha^n(b)=0$. Similarly $\alpha^n(a)b=0$.
\end{proof}

\medskip
Throughout this paper, all monoids and rings are with identity element that is inherited by submonoids and subrings and preserved under homomorphisms, but neither monoids nor rings are assumed to be - necessarily -  commutative.\par

\section{Archimedean skew power series rings}\label{archh}
In this section, we give the conditions on the ring $R$ such that the property of being Archimedean can be transferred from $R$ to the skew power series ring $ R[[x;\alpha]]$ and vice versa. The following results in \cite[Theorems 4.1]{paykan} describe such a transfer with the assumption that $ R[[S,\omega]]$ is a domain.

Let  $\omega: S\longrightarrow \End(R)$ be a monoid homomorphism, where $R$ is a ring and $(S,\le)$ is a strictly ordered monoid. If $R[[S,\omega]]$ is a left (resp. right) Archimedean domain, then one can see that $R$ is a left (resp. right) Archimedean domain. It can be proved straightforward from the definition and the fact that $\bigcap _{n\in \mathbb{N}}(c_a)^nR[[S,\omega]]= 0$ for nonunit $a$. Also, $S$ is a left (resp. right) Artinian and  narrow monoid because of the fact that $\bigcap _{n\in \mathbb{N}}(\textbf{e}_s)^nR[[S,\omega]]=0$ for nonunit $s\in S$. Moreover, $\omega _s$ is injective (and preserve nonunits of $R$) for any $s\in S$.

As mentioned in previous section, if $(S,.,\leq)$ is a (in particular positive) strictly totally ordered monoid and $0\neq f \in R[[S,\omega]]$, then $\supp( f )$ is a nonempty well-ordered subset of $S$. The smallest element of $\supp(f)$ is denoted by $\pi( f )$ and we set $\theta(f):=f(\pi(f))$. In particular, for a power series $f(x)=\sum f_nx^n$, with the first nonzero coefficient is $f_m$, we have $\pi(f)=m$ and $\theta(f)=f_m$.
 Let $U(R)$ be the set of unit elements of $R$. Then it is shown  in \cite[Proposition 3.2]{zim} that if $\pi(f) \in U(S)$ and $\theta (f) \in U(R)$, then $f \in U(R[[S,\omega]])$, where $f\in R[[S,\omega]] \setminus \{0\}$.
\begin{theorem}\cite[Theorems 4.1]{paykan}\label{paykan}
Let $R$ be a ring, $(S,\le)$ a strictly ordered monoid and $\omega: S\longrightarrow \End(R)$ a monoid homomorphism. Then we have the following:
\begin{itemize}
\item[{\em (i)}] The ring $R[[S,\omega]]$ is a right Archimedean domain if and only if $R$ is a
right Archimedean domain, $\omega _s$ is injective and preserves nonunits of $R$ for any $s\in S$.

\item[{\em (ii)}] The ring $R[[S,\omega]]$ is a left Archimedean domain if and only if $R$ is a left Archimedean domain and $\omega _s$ is injective for any $s\in S$.
\end{itemize}
\end{theorem}

 Before stating the main result of the paper, we first recall a few definitions. An ideal $P$ of $R$ is \textit{completely prime} if $ab\in P$ implies $a\in P$ or $b\in P$ for $a,b\in R$.  Let $\alpha$ an endomorphism of $R$. An ideal $I$ of $R$ is called an \textit{$\alpha$-ideal} if $\alpha (I)\subseteq I$ and $I$ is called \textit{$\alpha$-invariant} if $\alpha ^{-1}(I)=I$. Note that if $I$ is an $(\alpha)$-ideal, then $\bar{\alpha}:R/I\longrightarrow R/I$ defined by $\bar{\alpha}(a + I)=\alpha (a)+I$ for $a\in R$ is an endomorphism of the factor ring $R/I$.

We prove a similar result in the realm of reduced rings. Reduced rings are natural generalization of domains, for example direct product of domains is reduced. Therefore, Archimedean rings can also keep certain properties similar to the Archimedean domains.

\begin{theorem}\label{rigidar}
Let  $R$ be a $\alpha$-rigid ring satisfying the ACC on annihilators and
$\alpha$  be a surjective endomorphism.
Then:
\begin{enumerate}
\item[\em{(1)}] $R$ is right Archimedean and $\alpha$ preserves nonunits of $R$ if and only if the reduced ring $R[[x;\alpha]]$ is right Archimedean .
\item[\em{(2)}] $R$ is left Archimedean if and only if the reduced ring $R[[x;\alpha]]$ is left Archimedean.
\end{enumerate}

\end{theorem}

\begin{proof}
We will give the proof for the right-sided version. The proof for the left-sided version follows on the same line if the order of the multiplication is changed from right to left throughout the proof. We recall that for a power series $g(x)=g_0+g_1x+\cdots $ we have $\pi(g):=m$ where $g_m\neq 0$  and $g_i=0$ for $0\leq i<m$. Also $\theta(g)=g_m$ is the leading term of the power series. We set $\pi(0)=-\infty$.

\medskip
 First, assume that $R$ is right Archimedean.  According to \cite{moussavi}, if $R$ is rigid, then $R[[x;\alpha]]$ is a reduced ring.  Let $A=R[[x;\alpha]]$ and
 \[
 \Gamma := \{g \in A\; | \;  \text{ $g$ is a nonunit in $A$ and}\; \bigcap_{n\in \mathbb{N}}Ag^n\neq \{0\}\}.
 \]
We need to show that $\Gamma=\emptyset$. Assume that $\Gamma \neq \emptyset$.
We define  the set of leading terms $I_f:=\{\theta (g)|g\in fA\}\cup \{0\}$. Then, clearly $I_f$ is a right ideal of $R$.
Let
\[
T:=\{\Ann\big(\bigcup_{i\in \mathbb{N}}\;I_{c_ig}\big)|
 \;  c_i \in R,\; g\in \Gamma \; \text{and}\;  \; c_ig\neq 0,\;c_ig^i=c_jg^j \; \forall \; i,j\in \mathbb{N}\}.
\]
 Then $T \neq \emptyset$ because we assumed that $\Gamma \neq \emptyset$. Note that $\Ann\big(I_{c_ig}\big)=\Ann\big(I_{c_ig^i}\big)$ because $R$ is rigid. Considering that $R$ satisfies \emph{ACC} on annihilators, and then applying Zorn's Lemma results in the fact that $T$ has a maximal element. Let $V:=\Ann\big(\bigcup_{i\in \mathbb{N}}I_{a_if}\big)$ be the maximal element of $T$ for some $f\in
\Gamma$ such that
\begin{align}\label{sharto}

 a_i \in R, \; a_if\neq 0,\; \text{and}\; a_if^i=a_jf^j \; \forall \; i,j\in \mathbb{N}.
\end{align}

\textbf{Step 1.} \textit{We claim that $V$ is a two-sided completely prime ideal of $A$. }

From Lemma \ref{jav}, we see that $R$ is reduced because it is rigid. It implies that annihilators in $R$ are two-sided ideals. Let $ab\in V$ for some $a,b\in R$ and $b\notin V$. We need to show that $a\in V$.
First, we show that $a\in \Ann(I_{ba_if})$ for all $i \in \mbox{$\mathbb{N}$}$. Observe that if $ba_ifh=0$ for all $h\in A$ and $i \in \mbox{$\mathbb{N}$}$, then $b\in V$ which contradicts our assumption that $b\notin V$.  Therefore, $ba_jfA \neq \{0\}$ for some $j \in \mbox{$\mathbb{N}$}$. Take $ba_jfh \in ba_jfA $ such that $ba_jfh \neq 0$. This means that $\bigcup_i I_{ba_if}\neq 0$. Also, for all $i\in \mathbb{N}$ and for all $h\in A$ we have  $ab\theta (a_ifhb)=0$ because $ab \in V$ and $a_ifhb\in a_ifA$. Assume that $\theta (a_ifhb):=tb$, then $\theta (ba_ifh)=bt$ because $R$ is reduced. We conclude that $ab\theta (ba_ifh)=0$, because $R$ is reduced  - since in every reduced ring $wyz=0$ implies $wzy=0$.

Take $w=a$ and $y=bt$ and $z=b$.
So $abbt=ab^2t=0$. Again by using that $R$ is reduced, we get $abt=0$. So $a \theta (ba_i fh)=0$ for all $i \in \mbox{$\mathbb{N}$}$ and $h\in A$; and we conclude that $a\in \Ann (\bigcup_{i\in \mathbb{N}} I_{ba_if})$.

Let $z\in V$. So $z\theta(a_ifl)=0$ for all $i>0$ and $l\in A$. In particular, $z\theta(a_ifhb)=0$ for all $h\in A$. Considering that $R$ is reduced, we get $z\theta(ba_ifh)=0$. This implies that $z\in \Ann(\bigcup_{i\in \mathbb{N}} I_{ba_if})$, which means that
\begin{align}\label{shar}
V\subseteq \Ann(\bigcup_{i\in \mathbb{N}} I_{ba_if}).
\end{align}
Also, since  $ba_ifA \neq 0$, for some $i \in \mbox{$\mathbb{N}$}$ and $ba_if^i=ba_jf^j \; \forall \; i,j\in \mathbb{N}$, we see that $\Ann(\bigcup_{i\in \mathbb{N}} I_{ba_if})\in T$. The equation
(\ref{shar}) and the maximal property of $V$ in $T$ yields that
$V=\Ann(\bigcup_{i\in \mathbb{N}} I_{ba_if})$. This implies that
$a\in V$ and that $V$ is a two sided completely prime ideal.

\textbf{Step 2.}
\textit{We show that $V$ is $\alpha$-invariant.}

Let $r\in \alpha^{-1}(V)$. Then there exists $v\in V$ such that
$\alpha^{-1}(v)=r$ and we have $\alpha(r)=v$. By definition
of $V$, we have $\alpha(r)u=0$ for all $u\in \bigcup_{i\in
\mathbb{N}}(I_{a_if})$. Then, from Lemma~\ref{jav}, we have  $ru=0$ and it shows that
$r\in V$ and $\alpha^{-1} (V)\subseteq V$.

To show the reverse inclusion, take $v\in V$. Then for all $u\in \bigcup_{i\in \mathbb{N}}(I_{a_if})$ we have $vu=0$. Again, by Lemma~\ref{jav}, we have $\alpha(v)u=0$ . It shows that $\alpha(v)\in V$. Since $\alpha$ is surjective and rigid, hence an
automorphism, we see that $v \in \alpha^{-1}(V)$ and therefore, $V\subseteq \alpha^{-1}(V)$. From this we conclude that $V$ is an $\alpha$-invariant.

\textbf{Step 3.} \textit{ We show that $W:=(\frac{R}{V})[[x;\overline{\alpha}]]$ is a well-defined Archimedean domain.}

We already know that $V$ is a two-sided completely prime ideal, so the
factor ring $U:=\frac{R}{V}$ is a well-defined domain. Then it is enough to show that $U$ is an Archimedean domain because then by using
Theorem \ref{paykan} we can conclude that $U[[x;\overline{\alpha}]]$ is an Archimedean
domain, where $\overline{\alpha}: R/V \longrightarrow R/V$ defined by $\overline{\alpha}(a+V)=\alpha(a)+V$ for all $a\in R$ is an endomorphism of factor ring $R/V$.

To show that $U$ is an Archimedean
domain, we take  $\overline{z}\in \bigcap _{i\in
\mathbb{N}}U\bar{y}^i $. Then $\overline{z}=\bar{z_1} \bar{y}=\bar{z_2} \bar{y}^2=\cdots$ for some $\overline{z_i}\in U$. It gives, \begin{align}\label{farzado}
z+l_0=z_1y+l_1=z_2y^2+l_2=\cdots.
\end{align}
for some $l_i \in V$. We choose $h\in A$ and $i\in
\mathbb{N}$. By multiplying  \ref{farzado}  with $\theta (a_ifh)$ and using that $l_i \in V$ for all $i$, we obtain
$\theta(a_ifh)z=\theta(a_ifh)z_1y=\cdots$. Then $\theta(a_ifh)z=0$, because $R$
is Archimedean and $y$ is not unit in $R$, otherwise $\overline{y}$ becomes a unit in $U$. As $h$ and $i$ are
chosen arbitrarily, we obtain that $z\in \Ann( I_{a_if})$ for all $i\in \mathbb{N}$ and therefore, $z\in V$. This shows, $\overline{z}=0$, as required. Therefore, $U$ is Archimedean.

\textbf{Step 4.} \textit{We claim that $\overline{f}$ is nonunit in $W=R/V[[x;\overline{\alpha}]]$.}

Assume that there exists $g\in A$
such that $\overline{f}\overline{g}-\overline{1}=\overline{0}$. It means that all coefficients of $fg-1$ belong to $V$. In particular, $\theta (fg-1) \in V$ and we have  $\theta (fg-1)\theta (a_if)=0$ for all $i\in\mathbb{N}$. Assume that $f(x)=\sum f_lx^l$ and $g(x)=\sum g_l x^l$. We show that $f_0g_0$ is not a unit.

\textbf{Case 1.}  Assume that $\pi(fg-1)>0$.
Then we have
 $f_0g_0=\theta(fg)=1$.
Thus,
$f_0$ is a unit.
This implies that $f$ is a unit, which contradicts to our assumption that $f$ is non-unit. Therefore, we conclude that $\pi(fg-1)> 0$ is not possible.

\textbf{Case 2.} Now, let $\pi(fg-1)=0$. Note that $(fg)_0-1\in V$ and  $(fg)_0-1\neq 0,-1$, because otherwise it shows that  $V=A$, which is impossible. Since $(fg)_0\neq 0$, then $(fg)_0$ is a nonunit, otherwise $fg$ is unit, \cite[Proposition 3.2]{zim}, which is not possible.
 Therefore, $f_0g_0$ is a nonunit.

We assumed that $\overline{f}$ is a unit, which implies that $\overline{f}^i$ is a unit for all $i \in \mbox{$\mathbb{N}$}$. For each $i \in \mbox{$\mathbb{N}$}$, let $\overline{b_i}$ be the inverse of $\overline{f}^i$. Then we have  $\overline{b_i}\overline{f}^i=\overline{b_j}\overline{f}^j=1$, for all $i,j \in \mbox{$\mathbb{N}$}$.
It means that for all $i,j \in \mbox{$\mathbb{N}$}$ we have $\overline{b_i}\overline{f}^i-\overline{b_j}\overline{f}^j=0$ and  all the coefficients of $b_if^i-b_jf^j$ belong to $V$. It leads us to have
$\theta(a_1f)(b_if^i-b_jf^j)=0$  - because $R$ is reduced.

Therefore, we obtain
\begin{align}\label{equ}
\theta (a_1f)=\theta (a_1f)b_1f=\theta (a_1f)b_2f^2=\cdots.
\end{align}

We know from equation (\ref{sharto}) that  $\theta (a_1f)\neq 0$, which gives $\theta (a_1f)(b_if^i)_0\neq
0$.
 Let $b_i(x)=\sum b_{i,j}x^j$ and  $t_i:=\theta (a_1f)b_{i,0}$ and $r:=f_0$.
One can see that $(b_if^i)_0
=b_{i,0}f^i_0=b_{i,0}r^i$.
So $0\neq \theta (a_1f)(b_if^i)_0=t_ir^i$.
Thus, $t_{i+1}r^{i+1}=t_ir^i\neq 0$ for all $i$. Note that $r\neq 0$ is a non-unit because $f$ is non-unit. It sums up to the result that $R$  is non-Archimedean which is contradiction to our hypothesis. Therefore, $\overline{f}$ must be a non-unit.

\textbf{Step 5.} \textit{Finally, we are ready to show that $A$ is an Archimedean reduced ring.}

From equation (\ref{sharto}), we have
$\overline{a_1}\overline{f}=\overline{a_2}\overline{f^2}=\cdots$.
In Step 3 and 4, we showed that $W$ is Archimedean domain and $\overline{f}$ is a non-unit in $W$, so we must have $\overline{a_1}\overline{f}=\overline{0}$ and consequently, $\overline{a_1}=\overline{0}$ or
$\overline{f}=\overline{0}$. If $\overline{a_1}=\overline{0}$, then all coefficients of $a_1=\sum a_{1,j}x^j$ belong to $V$ and $a_{1,j}\theta(a_1f)=0$ for all $j\in \mathbb{N}$.
Again, by using Lemma \ref{jav} we obtain
\begin{align*}
\big( \theta (a_1f)\big)^2 =\bigg(\theta(a_1f)\bigg)\bigg(\theta(a_1f)\bigg)
&=\theta(a_1f)\left(\sum_{x+y=\pi(a_1f)}a_{1,x}\alpha^x(f_y)\right)\nonumber\\
&=\sum_{x+y=\pi(a_1f)}\left(\theta(a_1f)a_{1,x}\alpha^x(f_y)\right)
\end{align*}
Because $\theta (a_1f)a_{1,x}=0$ for all  $x\in S$, we see that
$\theta (a_1f)^2=0$. It gives that $\theta (a_1f)=0$ because $R$ is reduced. This gives us a contradiction to our assumption that $a_if \neq 0$ for all $i \in \mbox{$\mathbb{N}$}$. Therefore, $\overline{a_1} \neq 0$.
Then we must have $\overline{f}=\overline{0}$, that is, all coefficient of $f$ must belong to $V$ and $f_j\theta (a_1f)=0$ for all $j\in \mathbb{N}$. Then, from Lemma~\ref{jav}, we see that $\theta(a_1f)\alpha^l\left(f_j\right)=0$. This gives,
\begin{align}
\bigg(\theta(a_1f)\bigg)^2=\sum_{x+y=\pi(a_1f)}a_{1,x}\alpha^x(f_y)\theta(a_1f)=0
\end{align}
which again implies that $\theta((a_1f))^2=0$ and contradicts our assumption that $a_if \neq 0$ for all $i \in \mbox{$\mathbb{N}$}$. Therefore, our primary assumption $\Gamma \neq \emptyset$ is false. It shows that $A$ is an Archimedean ring, which completes the proof of the "only if" part.

Conversely,  assume that $\bigcap _{n=1}^{\infty}Ra^n \neq 0$ for some $a\in R$. Then $0\neq \bigcap _{n=1}^{\infty}Ra^n \subseteq \bigcap _{n=1}^{\infty} Aa^n$. But $A$ is a right Archemidean ring, so $a$ should be a unit in $A$. Thus, $fa=1$ for some $f\in A$. This means that $(fa)_0=1$.
So we get $f_0a=1$. Hence $a$ is a unit in $R$, which leads us to the conclusion that $R$ is a right Archemidean ring.

We show that $\alpha$ preserves nonunits. If $b\in R$ is a nonunit element and $\alpha(b)$ is a unit, then $c\alpha(b)=1$ for some $c\in R$. Assume that $f(x)b=1$, then $f_0b=1$ which is not possible. So $b$ is not a unit in $A$. Also, $c^nxb^{n}=c^n\alpha(b^n)x=x=c^mxb^m$. So $x\in \bigcap_{n} Ab^n$, while $b$ is not a unit in $A$. This contradiction shows that $\alpha(b)$ should not be a unit in $R$.
\end{proof}

In the case $\alpha=id_R$ - but $R$ is still not necessarily commutative - we conclude the next corollary immediately.

\begin{corollary}\label{archcor}
Let $R$ be a right (resp. left) Archimedean  reduced ring which satisfies the ACC on annihilators.
Then the power series ring $R[[x]]$ is a right (resp. left)
Archimedean  reduced ring.
\end{corollary}

In Theorem~\ref{rigidar}, we proved that under some conditions on the ring $R$, it is possible to transfer the property of being Archimedean from $R$ to $R[[x;\alpha]]$.
It is natural to ask that what are the minimum conditions required on $R$ to allow such a transfer from $R$ to $R[[S,\omega]]$. With the help of the next example, we justify that the condition rigid required in Theorem~\ref{rigidar} cannot be changed to reduced.

\begin{example}\label{reducekhubu}
Let $\mathbb{F}$ be any field and put $S:=(\mathbb{N}\cup \{0\})\times (\mathbb{Z})$
with usual addition and usual order.

 Consider $\mathbb{F}[[S]]= \mathbb{F}[[(\mathbb{N}\cup \{0\})\times (\mathbb{Z})]]:=T$. If $\mathbb{F}[[\mathbb{Z}]]:=R$, then one can see that $T=R[[\mathbb{N}\cup \{0\}]]:=R[[S]]$. Also, we can see that $R$ is isomorphic to the topological closure of $\mathbb{F}((x))$ where the metric is induced from the p-adic valuation $\deg(\cdot)$ and $x:=e_{(1,0)}\in T$ - we use the standard notation $e_s$ for $s\in S$ in the  skew generalized  power series rings.

Define $x:=e_{(1,0)}\in T$ and $y:=e_{(0,1)}\in T$.
One can see that $T$ is an skew generalized power series rings of $R$ where $S$ is a positive stirctly totally ordered monoid.
Then, $T$ is a Noetherian integral domain (i.e. a reduced ring). Also, $T$ satisfies both \emph{ACCP} and \emph{ACC} on annihilators.
Take $\alpha :T\rightarrow T$ defined by
\begin{align}
\alpha (1)=1, \quad \alpha (x)=y^{-1}, \quad \alpha (y)=y.
\end{align}
Note that $\alpha$ is not a monomorphism, and therefore, not rigid because $\alpha (xy^2-y)=0$ and $xy^2-y \neq 0$.  So $T$ is not an $\alpha$-rigid ring.

Let $U:=T[[z;\alpha]]$ and $f_k(z)=x\sum _{n=k}^{\infty}y^n z \in T$. Then
\begin{align}
f_{k+1}(z).x=\big (\sum _{n=k+1}^{\infty} xy^nz \big )x=\sum _{n=k+1}^{\infty}x\alpha (x)y^n z=\sum _{n=k+1}^{\infty}xy^{n-1}z=\sum _{n=k}^{\infty}xy^nz=f_k(z).
\end{align}
It shows that $ f_k U \subseteq  f_{k+1} U$ and we obtain an ascending chain of principal ideals in $T$. In order to show that $f_k U \subsetneq  f_{k+1} U$, it is enough to show that $x$ is not unit in $U$. Suppose that $h(z)x=1$ for some $h(z) \in U$. Let $h(z)=\sum _{i=0}^{\infty}g_i(x,y,y^{-1})z^i$, then
$h(z)x=g_0(x,y,y^{-1})x+(\sum _{i=1}^{\infty}g_i(x,y,y^{-1})z^i)x=1$.
It implies, $g_0(x,y,y^{-1})x=1$ showing that $x$ is a unit in $T$, which is not true. Therefore, we obtain a non-stablizing chain in $U$ given by $f_1 T \subsetneq  f_{2} U\subsetneq \ldots$. Hence, $U$ does not satisfy \emph{ACCPR}.
\end{example}

In \cite{lantz}, the authors give another example which also shows that the $\alpha$-rigid condition in Theorem~\ref{rigidar} is not superfluous. A ring is $\alpha$-compatible if $a\alpha(b)=0$ results in $ab=0$ - see more information about $\alpha$-compatible rings in \cite{moussavi}. It is known that the set of $\alpha$-rigid rings are the intersection of the set of $\alpha$-compatible rings and reduced rings. We proved that Archimedean property can be shifted from the $\alpha$-rigid rings and propose a counterexample for the reduced rings. The remaining natural question which seeks an answer is:
\begin{question}
Let  $R$ be an $\alpha$-compatible ring satisfying the \emph{ACC} on annihilators. Furthermore, let  $\alpha$ be a surjective endomorphism in $R$. Then $R$ is right (resp. left) Archimedean if and only if the reduced ring $R[[x;\alpha]]$ is right (resp. left) Archimedean.
\end{question}


\section{Introduction}
\label{sec:intro}

In real acoustic environments, a speech source of interest is frequently corrupted by interfering sounds. Human auditory system excels at attending to a target speech source, and the cocktail party problem~\cite{cherry1953some} aims to develop such capabilities in man-made devices and systems. A critical aspect of the cocktail party problem is speaker separation which aims to separate multiple concurrent speech signals of interest from a sound mixture.

Conventionally, most of the speaker separation methods work in time-frequency (T-F) domain where T-F representations are typically computed using short-time Fourier transform (STFT). In recent years, the performance of T-F domain speaker separation has substantially improved due to the use of deep learning~\cite{hershey2016deep, isik2016single, yu2017permutation, kolbaek2017multitalker, chen2017deep, wang2018alternative, wang2018end, wang2019deep, liu2019divide}. Moreover, the advent of deep learning based speech separation has also ignited interest in time-domain approaches, which directly estimate the waveform of clean speech from that of the mixture without resorting to a T-F representation. A notable time-domain speaker separation approach is TasNet~\cite{luo2019conv}, which yields comparable scale-invariant signal-to-noise ratios (SI-SNRs) and signal-to-distortion ratios (SDRs) to the ideal ratio mask (IRM). Other related studies include~\cite{stoller2018wave, venkataramani2018end, luo2020dual, zhang2020furcanext, nachmani2020voice} and~\cite{zeghidour2020wavesplit}.

\begin{figure*}[t]
	\vspace{-10pt}
	\centering
	\includegraphics[width=14cm]{figure2.pdf}
	\vspace{-5pt}
	\caption{(Color Online). Diagram of SISO SAGRNN. The dotted lines represent the procedures that exist only during training.}
	\label{fig:siso_sagrnn}
	\vspace{-8pt}
\end{figure*}

While several time-domain monoaural speaker separation methods have been developed, very few works have focused on binaural separation. Moreover, most existing binaural separation systems have a multiple-input single-output (MISO) layout, which produce a mono estimate for each of the target speakers from a binaural mixture~\cite{zohourian2016binaural, liu2018iterative, dadvar2019robust}. Hence these systems do not preserve interaural cues such as interaural time differences (ITDs) and interaural level differences (ILDs), which are crucial for human listeners to perform sound localization and lateralization~\cite{domnitz1977lateral, hershkowitz1969interaural}.

On the T-F domain front, various techniques have been developed to preserve binaural cues in the estimated signals. One can apply a common real-valued T-F mask or spectral gain to both left and right channels~\cite{mandel2009model, zohourian2018gsc}. Alternatively, binaural cues can be preserved by applying adaptive beamformers with additional constraints that encourage interaural cue preservation~\cite{van2007binaural, marquardt2015theoretical, hadad2015theoretical}. However, these techniques sacrifice the separation performance. More recently, a multiple-input multiple-output (MIMO) TasNet~\cite{han2020real} was designed, which produces a binaural estimate for each speaker. MIMO TasNet yields significantly better speech quality than the single-channel TasNet while preserving both ITDs and ILDs.

In this paper, we propose a novel framework called \emph{multiple-input multiple-output self-attentive gated recurrent neural network} (MIMO SAGRNN) for binaural speaker separation. The proposed SAGRNN network architecture extends the gated RNN in~\cite{nachmani2020voice} by additionally incorporating self-attention mechanisms and dense connectivity (DC). We then derive MIMO SAGRNN from a single-input single-output (SISO) SAGRNN by first extending the SISO SAGRNN into a multiple-input single-output (MISO) layout by creating two encoders, one for the \emph{reference ear} input and the other for the \emph{non-reference ear} input. This MISO SAGRNN estimates the separated signals in the reference ear. The MIMO system is formulated by alternately treating each ear as the reference ear, yielding estimates for both ears in a symmetric manner.

The rest of this paper is organized as follows. Section~\ref{sec:algorithm} describes our proposed approach. The experimental results are presented in Section~\ref{sec:exp}, and Section~\ref{sec:conclusion} concludes this paper.

\vspace{-0.2cm}
\section{Algorithm Description}
\label{sec:algorithm}

We progressively develop a MIMO system for binaural speaker separation. Specifically, we start with a SISO SAGRNN architecture, and then present the MIMO setup.

\vspace{-10pt}
\subsection{SISO SAGRNN}
As in~\cite{nachmani2020voice}, the separation framework of a SISO SAGRNN comprises three stages: \emph{encoding and chunking}, \emph{block processing}, and \emph{decoding and overlap-add}. A time-domain input mixture is transformed into a set of overlapped chunks via encoding and chunking, which leads to a 3-D embedding. Subsequently, the 3-D embedding is passed into stacked RNN blocks to perform intra-chunk (local) and inter-chunk (global) modeling alternately and iteratively. The 3-D representation learned by the last RNN block is decoded and then transformed back to the time domain by an overlap-add operator.

Given a $T$-sample input waveform $y \in \mathbb{R}^T$, an encoder is used to segment and encode $y$ into $L$ overlapped time frames with a frame size of $P$ and a hop size of $P/2$, yielding a 2-D embedding $\mathbf{U} \in \mathbb{R}^{N \times L}$. Specifically, the encoder consists of a 1-D strided convolutional layer with $N$ kernels, followed by a rectified linear activation function. We divide the time frames into $S$ overlapped chunks with a chunk size of $R$ and a hop size of $R/2$. These chunks are then concatenated into a 3-D embedding $\widetilde{\mathbf{W}} = [\mathbf{W}_1, \dots, \mathbf{W}_S] \in \mathbb{R}^{N \times S \times R}$, where $\mathbf{W}_1, \dots, \mathbf{W}_S \in \mathbb{R}^{N \times R}$ are the 2-D chunks.

\begin{figure*}[t]
	\vspace{-11pt}
	\centering
	\includegraphics[width=16cm]{figure3.pdf}
	\vspace{-8pt}
	\caption{(Color Online). Diagram of the SA-MULCAT block. The dotted lines indicate that the procedures exist in all SA-MULCAT blocks except the first one. The symbol $\bigotimes$ represents the element-wise multiplication, and $\bigoplus$ the element-wise addition.}
	\label{fig:SA_MULCAT}
	\vspace{-15pt}
\end{figure*}

\begin{figure}[t]
	\centering
	\subfigure[]{%
		\centering
		\includegraphics[width=3.0cm]{figure4.pdf}
		\label{fig:subfig:self-attention}
	}\hspace{1mm}%
	\subfigure[]{%
		\centering
		\includegraphics[width=5.5cm]{figure5_t.pdf}
		\label{fig:subfig:miso_mimo}
	}%
	\vspace{-5pt}
	\caption{(Color Online). (a) Diagram of the self-attention block. (b) Overview of the MIMO separation system with an underlying MISO system.}
	\label{fig:sa_mimo}
	\vspace{-10pt}
\end{figure}

Subsequently, the 3-D embedding $\widetilde{\mathbf{W}}$ is fed into a series of $B$ RNN blocks for processing. To improve the information flow between blocks, we propose a dense connectivity pattern: each block receives the outputs of all preceding blocks, i.e. $\widetilde{\mathbf{W}}_b = \mathcal{H}_b ([\widetilde{\mathbf{W}}_0, \dots, \widetilde{\mathbf{W}}_{b-1}])$ for $b = 1, \dots, B$,
where $\mathcal{H}_b$ denotes the mapping function defined by the $b$-th block, and $\left[\cdot, \dots, \cdot \right]$ the concatenation operation. The output embedding of the $b$-th block is represented by $\widetilde{\mathbf{W}}_b$, where $\widetilde{\mathbf{W}}_0 = \widetilde{\mathbf{W}}$ and $\widetilde{\mathbf{W}}_b \in \mathbb{R}^{N \times S \times R}, \forall b$.

Similar to~\cite{nachmani2020voice}, we use a multi-scale loss for training, which necessitates producing a waveform estimate for each speaker
after each block. We decode the output embedding of each block with a decoder, which comprises a parametric rectified linear function~\cite{he2015delving} followed by a 2-D 1$\times$1 convolutional layer with $C\cdot N$ kernels. The decoded feature of size $CN \times S \times R$ is divided into $C$ 3-D representations of size $N \times S \times R$, corresponding to the $C$ speech sources. These 3-D representations are transformed back to waveforms by two successive overlap-add operations at the chunk level and the frame level, respectively. Note that the same decoder is applied to the output of each block. Fig.~\ref{fig:siso_sagrnn} depicts the SISO SAGRNN.

A series of multiply-and-concatenate (MULCAT) blocks are employed to model the intra-chunk and inter-chunk dependencies. In this study, we extend the MULCAT block by introducing self-attention~\cite{vaswani2017attention}, which amounts to a self-attention based MULCAT (SA-MULCAT) block illustrated in Fig.~\ref{fig:SA_MULCAT}. The concatenation of paths from the dense connections is fed into a linear projection layer for dimension reduction, yielding an embedding of size $N \times S \times R$. The resulting embedding is successively passed through two subblocks, one for intra-chunk modeling and the other for inter-chunk modeling. In each subblock, we employ a self-attention block followed by a gated RNN module, which consists of two bidirectional long short-term memory (BLSTM) layers coupled with each other. Each BLSTM contains $H$ units in each direction. The Hadamard product of their outputs is concatenated with the input to the gated RNN module, and then passed into a linear projection layer for dimension reduction. In addition, a skip connection is used to bypass the subblock. After the first subblock, the dimensions of the 3-D representation are re-permuted, so that sequential modeling can be performed across chunks in the second subblock. After the second subblock, the dimensions are re-permuted back.

The self-attention block is illustrated in Fig.~\ref{fig:subfig:self-attention}. We first divide a 3-D representation into a set of 2-D slices $\mathbf{Z} \in \mathbb{R}^{M \times N}$, where $M = R$ for intra-chunk modeling and $M = S$ for inter-chunk modeling. Each slice is linearly projected to a query matrix $\mathbf{Q}$, a key matrix $\mathbf{K}$ and a value matrix $\mathbf{V}$ by three different projection layers, where $\mathbf{Q}, \mathbf{K}, \mathbf{V} \in \mathbb{R}^{M \times D}$ and $D$ is set to 64. We apply a scaled dot-product attention function:
\vspace{-6pt}
\begin{equation}
Attention(\mathbf{Q}, \mathbf{K}, \mathbf{V}) = SoftMax(\frac{\mathbf{Q}\mathbf{K}^\top}{\sqrt{D}}) \mathbf{V},
\vspace{-6pt}
\end{equation}
where $SoftMax(\cdot)$ denotes the softmax function across columns. The output of the attention function is computed as a weighted sum of the values, where the weight assigned to each value is derived by measuring the similarities between the queries and the keys. Subsequently, all the attention output slices are merged and then linearly projected back to the size of the input 3-D representation. With a skip connection, this representation is concatenated with the input to the self-attention block, and then projected back to the original size. By leveraging the relevance among features at different time steps, the self-attention produces a dynamic representation in adapting to different acoustic conditions.

\vspace{-6pt}
\subsection{MIMO SAGRNN}
As shown in Fig.~\ref{fig:subfig:miso_mimo}, a \emph{reference encoder} and a \emph{non-reference encoder} are employed to process the binaural mixture waveforms. The resulting 2-D embeddings are concatenated and then linearly projected to the size of $N \times L$. Subsequently, we successively perform block processing, decoding and overlap-add, akin to the SISO system. In this MISO system, the separation outputs always correspond to the reference ear. We formulate the MIMO system by alternately treating each ear as the reference. Specifically, the separation outputs for the left ear are obtained by treating the left ear as the reference ear and the right ear as the non-reference. The separation outputs for the right ear are obtained by swapping the inputs of the two ears.

\vspace{-6pt}
\subsection{Training Objective}
We use the plain SNR rather than the widely-used SI-SNR~\cite{luo2019conv} as the training objective. The rationale is that SI-SNR training cannot preserve the ILD in the binaural estimates, as the power scale of the estimated signals is insusceptible to training due to the scale invariance.
The mean of the SNR losses from all SA-MULCAT blocks is used for training. The waveforms of the clean speech signals are used as the training target for calculating the losses from all blocks. In addition, we apply the permutation invariant training~\cite{kolbaek2017multitalker} criterion to the loss from each block individually, which allows the label permutation to change from one block to another.

\vspace{-0.2cm}
\section{Experiments}
\label{sec:exp}

\subsection{Experimental Setup}
We simulate a noise-free dataset and a noisy dataset from the WSJ0-2mix dataset~\cite{hershey2016deep}, which contains 20,000, 5,000 and 3,000 mixtures in the training, validation and testing sets, respectively. For both datasets, we convolve each pair of utterances in WSJ0-2mix with two randomly sampled head-related impulse responses (HRIRs) from the CIPIC HRTF Database~\cite{algazi2001cipic} respectively, which contains 45 subjects with 25~(azimuths)~$\times$~50~(elevations) directions for each subject. Specifically, we choose 35 subjects for training and cross validation, and use the 10 remaining subjects for testing. For the noisy dataset, we additionally simulate uncorrelated noise sources by randomly selecting HRIRs for them, where the number of noise sources is randomly sampled between 1 and 10. Note that all sound sources are placed in different directions. We use a set of roughly 65,000 noises from the DNS Challenge~\cite{reddy2020interspeech} for training and cross validation, and a different set~\cite{Robinhood76} of roughly 1,300 noises for testing. The SNR (w.r.t. the speech mixture in the left ear) is randomly chosen between -10~dB and 10~dB. All signals are sampled at 8~kHz.

We train the models on 4-second segments with the AMSGrad optimizer~\cite{reddi2018convergence} with a minibatch size of 4. The learning rate is initialized to 0.0002, which decays by 0.98 every 2 epochs. Gradient clipping with a maximum $\ell^2$ norm of 3 is applied during training. The network hyperparameters for MIMO SAGRNN are as follow: $P=8$, $N=128$, $R=126$, $H=128$, $D=64$ and $B=6$. Note that the value of $R$ is selected such that $R \approx S$ for the training segments.

We use several monaural and binaural separation models for comparison; the monaural models are trained and evaluated on each ear individually. Specifically, we use TasNet~\cite{luo2019conv}, dual-path RNN (DPRNN)~\cite{luo2020dual} and the gated RNN in~\cite{nachmani2020voice} as monaural baselines. We use MIMO TasNet~\cite{han2020real} as a binaural baseline. We slightly adjust the hyperparameter configurations of all baselines, so that they have comparable model sizes to our MIMO SAGRNN. For the noncausal temporal convolutional network (TCN) in TasNet~\cite{luo2019conv}, the number of repeated stacks is set to 4. For DPRNN, the number of kernels in the encoder and the decoder is set to 128, and the number of units in each direction for each BLSTM to 200. For MIMO TasNet, we replace the causal TCN by a noncausal TCN, with bottleneck size of 128. In addition, the number of kernels in the encoder and the decoder of MIMO TasNet is set to 512.

\vspace{-0.3cm}
\subsection{Experimental Results}
\subsubsection{Separation Results and Analysis}
\begin{table}[t]
	\vspace{-4pt}
	\caption{Comparison of different systems in the noise-free condition.}
	\centering
	\resizebox{0.48\textwidth}{!}{
		\begin{tabular}{|l|c|c|c|c|c|}
			\hline
			Metrics & \# param. & $\Delta$SDR (dB) & $\Delta$SNR (dB) & ESTOI (\%) & PESQ \\
			\hline
			Mixture & - & 0.00 & 0.00 & 56.10 & 1.99 \\
			\hline
			TasNet & 6.66~M & 15.78 & 15.80 & 90.25 & 3.30 \\
			DPRNN & 6.99~M & 18.25 & 18.28 & 92.64 & 3.49 \\
			Gated RNN & 7.56~M & 18.37 & 18.41 & 92.52 & 3.52 \\
			IBM & - & 13.69 & 13.61 & 88.49 & 3.36 \\
			IRM & - & 13.05 & 13.07 & 93.33 & 3.73 \\
			PSM & - & 16.77 & 16.58 & 95.55 & 3.91 \\
			\hline
			MIMO TasNet & 7.32~M & 21.14 & 20.69 & 95.53 & 3.73 \\
			Oracle MB-MVDR & - & 17.13 & 10.44 & 95.77 & 3.66 \\
			\hline
			\rowcolor{mygray}
			MIMO SAGRNN & 8.71~M & \textbf{27.19} & \textbf{26.88} & \textbf{98.08} & \textbf{4.06} \\
			\hspace{1mm}$-$ half multi-scale loss~(i) & 8.71~M & 24.35 & 23.99 & 96.70 & 3.90 \\
			\hspace{1mm}$-$ multi-scale loss~(ii) & 8.71~M & 22.93 & 22.45 & 95.69 & 3.79 \\
			\hspace{1mm}$-$ DC~(iii) & 8.38~M & 24.02 & 23.64 & 96.54 & 3.87 \\
			\hspace{1mm}$-$ SA~(iv) & 7.92~M & 23.13 & 22.75 & 96.12 & 3.82 \\
			\hspace{1mm}$-$ DC $-$ SA~(v) & 7.59~M & 21.97 & 21.64 & 95.66 & 3.73 \\
			\hline
		\end{tabular}
	}
	\label{tab:noise_free}
	\vspace{-10pt}
\end{table}

\begin{table}[t]
	\caption{Comparison of different systems in the noisy condition.}
	\centering
	\resizebox{0.4\textwidth}{!}{
		\begin{tabular}{|l|c|c|c|c|c|}
			\hline
			Metrics & $\Delta$SDR (dB) & $\Delta$SNR (dB) & ESTOI (\%) & PESQ \\
			\hline
			Mixture & 0.00 & 0.00 & 27.80 & 1.49 \\
			\hline
			TasNet & 11.71 & 13.33 & 54.65 & 2.12 \\
			DPRNN & 11.91 & 13.55 & 53.72 & 2.17 \\
			Gated RNN & 12.89 & 14.27 & 59.09 & 2.30 \\
			IBM & 14.11 & 14.64 & 71.19 & 2.63 \\
			IRM & 13.19 & 13.97 & 84.22 & 3.33 \\
			PSM & 17.41 & 17.62 & 89.40 & 3.55 \\
			\hline
			MIMO TasNet & 14.40 & 15.23 & 63.79 & 2.41 \\
			Oracle MB-MVDR & 4.98 & 4.90 & 42.71 & 1.79 \\
			\rowcolor{mygray}
			MIMO SAGRNN & \textbf{17.53} & \textbf{17.95} & \textbf{75.14} & \textbf{2.78} \\
			\hline
		\end{tabular}
	}
	\label{tab:noisy}
	\vspace{-5pt}
\end{table}

Tables~\ref{tab:noise_free} and~\ref{tab:noisy} show comparison among different approaches for the noise-free and noisy conditions respectively. The separation results are reported in terms of SDR improvement ($\Delta$SDR), SNR improvement ($\Delta$SNR), extended short-time objective intelligibility (ESTOI)~\cite{jensen2016algorithm}, and perceptual evaluation of speech quality (PESQ)~\cite{rix2001perceptual}. We can observe that MIMO TasNet produces consistently better results than the monaural baselines. Moreover, our proposed MIMO SAGRNN substantially outperforms MIMO TasNet in all the four metrics. For the noise-free condition, MIMO SAGRNN improves SDR by 6.05~dB and SNR by 6.19~dB over MIMO TasNet. Some demos can be found at https://jupiterethan.github.io/sagrnn.github.io/.

In addition, we compare MIMO SAGRNN with several oracle approaches, including ideal binary mask (IBM), ideal ratio mask (IRM), phase-sensitive mask (PSM)~\cite{erdogan2015phase} and an oracle masking-based minimum variance distortionless response (MB-MVDR) beamformer. We use an open-source implementation~\cite{beamformers} of the oracle MB-MVDR beamformer, with a frame length of 64~ms and a frame shift of 32~ms. The IRM is used to calculate the spatial covariance matrices. We alternately treat each channel as the reference channel to produce the binaural estimate. As shown in Table~\ref{tab:noise_free}, our approach consistently outperforms the ideal masks and the oracle beamformer in the noise-free condition. In the noisy condition (Table~\ref{tab:noisy}), our approach produces slightly higher SDR and SNR but lower ESTOI and PESQ than the PSM. Note that, ESTOI and PESQ improvements over the mixtures using the oracle MB-MVDR beamformer dramatically decrease in the noisy condition compared with the noise-free condition. This is likely because the directionality of the sound sources is smeared due to the presence of multiple noise sources. In contrast, our approach is more robust against the noise field.

\subsubsection{Ablation Study}
We conduct an ablation study to understand the contribution of each component in our approach. Several variants of MIMO SAGRNN are compared in Table~\ref{tab:noise_free}: (i) using multi-scale loss computed from only the last three RNN blocks; (ii) using loss computed from only the last RNN block; (iii) without dense connections; (iv) without self-attention blocks; (v) without dense connections and self-attention blocks. It is shown that self-attention and dense connectivity are crucial for MIMO SAGRNN. Without self-attention, for example, SDR decreases by 4.06~dB and SNR by 4.13~dB.

\begin{table}[t]
	\vspace{-14pt}
	\caption{Comparison between SISO and MIMO SAGRNNs on the better-ear channel for the noisy condition.}
	\centering
	\resizebox{0.4\textwidth}{!}{
		\begin{tabular}{|l|c|c|c|c|c|}
			\hline
			Metrics & $\Delta$SDR (dB) & $\Delta$SNR (dB) & ESTOI (\%) & PESQ \\
			\hline
			Mixture & 0.00 & 0.00 & 31.84 & 1.57 \\
			\hline
			SISO SAGRNN & 13.19 & 14.06 & 66.91 & 2.49 \\
			MIMO SAGRNN & 15.94 & 16.32 & 75.66 & 2.79 \\
			\hline
		\end{tabular}
	}
	\label{tab:better_ear}
	\vspace{-8pt}
\end{table}

\begin{table}[t]
	\caption{Evaluation of interaural cue preservation with binaural sound localization for the noise-free condition.}
	\centering
	\resizebox{0.48\textwidth}{!}{
		\begin{tabular}{|l|c|c|c|c|c|}
			\hline
			Metrics & $\Delta$Azimuth ($^\circ$) & $\Delta$ITD ($\mu s$) & \multicolumn{3}{c|}{$\Delta$ILD (dB)} \\
			\hline
			Frequency Channels & - & - & 2.07~kHz & 3.08~kHz & 3.75~kHz \\
			\hline
			Mixture & 26.03 & 255.85 & 4.31 & 4.85 & 5.12 \\
			\hline
			TasNet & 14.40 & 25.71 & 0.62 & 0.68 & 0.95 \\
			DPRNN & 13.64 & 25.71 & 0.62 & 0.68 & 0.93 \\
			Gated RNN & 13.82 & 23.42 & 0.91 & 0.73 & 1.22 \\
			IBM & 12.25 & 15.75 & 0.90 & 0.98 & 1.16 \\
			IRM & 9.43 & 55.00 & 0.29 & 0.28 & 0.39 \\
			PSM & 4.95 & 16.20 & 0.37 & 0.40 & 0.55 \\
			\hline
			MIMO TasNet & 6.45 & 20.10 & 0.88 & 0.77 & 1.13 \\
			Oracle MB-MVDR & 34.05 & 371.35 & 6.31 & 7.06 & 7.39 \\
			\rowcolor{mygray}
			MIMO SAGRNN & \textbf{5.88} & \textbf{14.95} & \textbf{0.53} & \textbf{0.45} & \textbf{0.70} \\
			\hline
		\end{tabular}
	}
	\label{tab:localization}
	\vspace{-10pt}
\end{table}

We also compare the SISO and MIMO SAGRNNs on the better-ear channel; {\it better ear} is defined as the ear that is closer to the target speech source. The azimuth position of the speech source is used to determine the better ear. Table~\ref{tab:better_ear} shows that MIMO SAGRNN yields significantly better results than SISO SAGRNN on all metrics, evidencing that the binaural inputs are effectively leveraged by our MIMO system.

\subsubsection{Evaluation of Interaural Cue Preservation}

\begin{figure}[t]
    \vspace{-14pt}
	\centering
	\subfigure[ITD Dist. (Clean)]{%
		\centering
		\includegraphics[width=4.25cm]{figure6a_t.pdf}
		\label{fig:subfig:itd_sphs}
	}%
	\subfigure[ITD Dist. (Estimated)]{%
		\centering
		\includegraphics[width=4.25cm]{figure6b_t.pdf}
		\label{fig:subfig:itd_sphs_est}
	}%

	\subfigure[ILD Dist. at 3.75~kHz (Clean)]{%
		\centering
		\includegraphics[width=4.25cm]{figure6c_t.pdf}
		\label{fig:subfig:ild_sphs}
	}%
	\subfigure[ILD Dist. at 3.75~kHz (Estimated)]{%
		\centering
		\includegraphics[width=4.25cm]{figure6d_t.pdf}
		\label{fig:subfig:ild_sphs_est}
	}%

	\caption{(Color Online). An example of ITD and ILD distributions.}
	\label{fig:hist_itd_ild}
	\vspace{-15pt}
\end{figure}

Lastly, we evaluate the preservation of interaural cues in the estimated binaural signals under noise-free condition. We apply a binaural sound localization algorithm~\cite{may2011probabilistic} to the binaural estimates, of which an open-source implementation is available. This implementation estimates the azimuth position of the sound source at the frame level, as well as the ITD and the ILD for each T-F unit of a 32-channel cochleagram based on a gammatone filterbank. The average frame-level azimuth errors are presented in Table~\ref{tab:localization}. Given the dominance of the ITD cue at low frequencies (below 1.5~kHz) in sound localization~\cite{wightman1992dominant}, we only take into account the frequency bands corresponding to gammatone filters with a maximum center frequency of approximately 1.5~kHz. Since ILD is highly frequency-dependent due to diffraction and attenuation of the sounds, we calculate the average ILD errors individually for three empirically selected frequency channels, corresponding to the gammatone filters with the center frequencies of roughly 2.07, 3.08 and 3.75~kHz. Given the fact that all sound sources are stationary in this study, we summarize only one ITD/ILD from an entire utterance in the following way. We plot a histogram of the T-F unit level ITDs/ILDs, and then estimate the ITD/ILD based on the center value of the highest bin. The number of bins is empirically set to 500 for ITD and 40 for ILD. An example of ITD and ILD histograms are presented in Fig.~\ref{fig:hist_itd_ild}. As shown in Table~\ref{tab:localization}, our approach reduces the azimuth error by 20.15$^\circ$ compared to the mixtures. Moreover, our approach yields consistently smaller azimuth, ITD and ILD errors than MIMO TasNet and the monaural baselines, showing that our approach preserves the interaural cues more effectively.

\vspace{-0.2cm}
\section{Conclusion}
\label{sec:conclusion}
We have proposed an end-to-end MIMO system for binaural speaker separation with interaural cue preservation. We developed a novel framework which relies on self-attention and dense connectivity for improved speaker separation. Our experimental results show that the proposed approach significantly outperforms a binaural separation approach (i.e. MIMO TasNet) in terms of $\Delta$SDR, $\Delta$SNR, ESTOI and PESQ. Moreover, our approach effectively preserves the auditory spatial cues of talkers. For future work, we would devote more efforts to the design of MIMO systems for real-time binaural speaker separation, as well as exploring separation of more concurrent speakers.

\newpage
\bibliographystyle{abbrv}