XX Month, XXXX XX Month, XXXX XX Month, XXXX XX Month, XXXX XX Month, XXXX
OJSP.2023.1234567

CORRESPONDING AUTHOR: Paul Didier (e-mail: paul.didier@esat.kuleuven.be). This
research work was carried out in the frame of the European Union’s Horizon
2020 research and innovation programme under the Marie Skłodowska-Curie Grant
Agreement No. 956369: “Service-Oriented Ubiquitous Network-Driven Sound —
SOUNDS”. The research leading to these results has received funding from the
European Research Council under the European Union’s Horizon 2020 research and
innovation program / ERC Consolidator Grant: SONORA (no. 773268). This paper
reflects only the authors’ views and the Union is not liable for any use that
may be made of the contained information. The scientific responsibility is
assumed by the authors.

# Sampling Rate Offset Estimation and Compensation for Distributed Adaptive
Node-Specific Signal Estimation in Wireless Acoustic Sensor Networks

PAUL DIDIER1  TOON VAN WATERSCHOOT1  SIMON DOCLO2  AND MARC MOONEN1  STADIUS
Center for Dynamical Systems, Department of Electrical Engineering (ESAT),
KU Leuven, 3001 Leuven, Belgium  Signal Processing Group, Department of
Medical Physics and Acoustics and Cluster of Excellence Hearing4all,
University of Oldenburg, Oldenburg, Germany

###### Abstract

Sampling rate offsets (SROs) between devices in a heterogeneous wireless
acoustic sensor network (WASN) can hinder the ability of distributed adaptive
algorithms to perform as intended when they rely on coherent signal
processing. In this paper, we present an SRO estimation and compensation
method to allow the deployment of the distributed adaptive node-specific
signal estimation (DANSE) algorithm in WASNs composed of asynchronous devices.
The signals available at each node are first utilised in a coherence-drift-
based method to blindly estimate SROs which are then compensated for via phase
shifts in the frequency domain. A modification of the weighted overlap-add
(WOLA) implementation of DANSE is introduced to account for SRO-induced full-
sample drifts, permitting per-sample signal transmission via an approximation
of the WOLA process as a time-domain convolution. The performance of the
proposed algorithm is evaluated in the context of distributed noise reduction
for the estimation of a target speech signal in an asynchronous WASN.

Sampling rate offsets, coherence drift, signal enhancement, weighted overlap-
add, wireless acoustic sensor networks

##  1 Introduction

Wireless acoustic sensor networks (WASNs) have been a subject of great
interest in recent years as they provide a number of advantages over
centralised systems performing audio signal processing tasks [1]. Novel
algorithmic solutions aim to utilise the increased flexibility and scalability
of WASNs in order to tackle various audio signal processing challenges in a
distributed fashion, bypassing the need for a data fusion centre with which
all nodes communicate.

This paper focuses on distributed signal estimation, where each node in the
WASN aims to estimate a node-specific desired signal. The distributed adaptive
node-specific signal estimation (DANSE) algorithm was originally formulated in
[2, 3] to tackle this problem. This algorithm is designed to allow each node
of a fully connected WASN to achieve centralised performance upon convergence
by iteratively computing its own multichannel Wiener filter (MWF) while only
exchanging single-channel signals with other nodes. DANSE can significantly
reduce the required number of signals communicated between nodes in a WASN
with many sensors per node, compared to a centralised MWF where nodes
communicate with a single fusion centre. Although the DANSE algorithm has been
tested under various conditions and for different tasks [4], a key aspect
allowing its robust deployment in real-world scenarios has yet to be
addressed, namely, signals asynchronicity.

In many practical applications such as teleconferencing systems or smart
domotics, the WASN is heterogeneous, i.e., composed of various interconnected
devices such as laptops, tablets, or hearing aids. Each device samples the
incoming acoustic information at a specific rate via its own analog-to-digital
converter based on an internal clock, the sampling rate of which may differ
from the nominal value provided by the manufacturer [5]. The sampling rate
mismatch between two devices can be quantified as the sampling rate offset
(SRO), generally expressed in parts-per-million (PPM). SROs in the range of
$\pm$500 PPM have been measured between commonly used devices and reported in
[5]. The same study showed that SROs can slowly vary through time, e.g., when
the devices undergo significant temperature changes or fluctuations in supply
voltage.

SROs lead to an increasing time-drift between signals sampled by different
clocks, which inhibits their use in algorithms that rely on coherent signal
processing [6]. Notably, the performance of signal enhancement algorithms
based on the MWF such as the DANSE algorithm depends on the computation of
accurate spatial covariance matrices. DANSE can thus be expected to be
sensitive to a lack of synchronicity between locally recorded microphone
signals and signals received from other nodes. In fact, literature around
DANSE has so far assumed that all nodes involved in the algorithm have exactly
the same sampling rate [2, 3, 7, 8, 9]. The asynchronicity problem in WASNs
has recently been investigated in the context of algorithms other than DANSE
[10, 11].

In this paper, we propose a methodology to relax the synchronicity assumption
in DANSE, bringing this algorithm closer to robust deployment in real-life
scenarios. The presence of SROs is addressed in a fully connected WASN where
node and source positions are fixed. Time-invariant SROs are considered,
assuming that no temperature or supply voltage changes occur during the
convergence phase of the algorithm. Per-node-pair SRO estimation is performed
blindly based on a coherence-drift method [12, 13]. The weighted overlap-add
(WOLA) implementation of the generalised eigenvalue decomposition (GEVD-)DANSE
algorithm [14] is modified to permit detection of full-sample drifts (FSDs)
via per-sample signal broadcasting. This is achieved by approximating the WOLA
process used for local signal fusion (analysis, filtering in the short-time
Fourier transform (STFT) domain, and synthesis) as a single time-domain
convolution operation. This method allows to retain the low complexity of WOLA
processing for the more costly steps of GEVD-based filter update and desired
signal estimation. The estimated SROs and the detected FSDs are then
compensated for via phase shifts in the STFT-domain. The performance of the
proposed algorithm is evaluated in the context of distributed noise reduction
for the estimation of a target speech signal.

The paper is organised as follows. In Section 2, the centralised GEVD-MWF is
reviewed. The key aspects of the theory and implementation of the DANSE
algorithm are summarised in Section 3. The proposed method for SRO estimation
and compensation within the DANSE framework is presented in detail in Section
4. The performance of the proposed method is then analysed by means of
simulations in asynchronous WASNs in Section 5. Finally, conclusive remarks
are formulated in Section 6.

##  2 GEVD-MWF-based signal estimation

A WASN composed of $K$ nodes is considered, where each node
$k\in\mathcal{K}=\\{1,\ldots,K\\}$ has $M_{k}\geq 1$ microphones. The total
number of microphones in the network is denoted by
$M=\sum_{k\in\mathcal{K}}M_{k}$. In the acoustic scene, one localised static
desired signal source (e.g., a talker) and $J\geq 1$ localised static noise
sources are present. The signals recorded by node $k$ can be represented in
the STFT domain at frame $i$ and frequency bin $\nu$ via an additive-noise
signal model:

$\mathbf{y}_{k}[\nu,i]=\mathbf{s}_{k}[\nu,i]+\mathbf{n}_{k}[\nu,i],$ (1)

where $\mathbf{y}_{k}[\nu,i]$, $\mathbf{s}_{k}[\nu,i]$, and
$\mathbf{n}_{k}[\nu,i]$ are $M_{k}$-dimensional vectors corresponding to the
microphone signals, the desired signal components of these signals, and the
noise components, respectively. The additive noise is assumed to be
uncorrelated with the desired signal.

In centralised processing, the signal vector available at the fusion centre is
defined as an $M$-dimensional stacked version
$\mathbf{y}[\nu,i]=\begin{bmatrix}\mathbf{y}_{1}^{T}[\nu,i],\dots,\mathbf{y}_{K}^{T}[\nu,i]\end{bmatrix}^{T}$
of the node-specific microphone signals where $(\cdot)^{T}$ denotes the
transpose operation. Similarly to (1), this vector can be expressed as
$\mathbf{y}[\nu,i]=\mathbf{s}[\nu,i]+\mathbf{n}[\nu,i]$ with
$\mathbf{n}[\nu,i]=\begin{bmatrix}\mathbf{n}_{1}^{T}[\nu,i],\dots,\mathbf{n}_{K}^{T}[\nu,i]\end{bmatrix}^{T}$
and
$\mathbf{s}=\begin{bmatrix}\mathbf{s}_{1}^{T}[\nu,i],\dots,\mathbf{s}_{K}^{T}[\nu,i]\end{bmatrix}^{T}$.

The objective of node $k$ is then to estimate a local desired signal
$d_{k}[\nu,i]$ based on $\mathbf{y}[\nu,i]$. Define without loss of generality
(w.l.o.g.) the desired signal $d_{k}[\nu,i]$ at node $k$ to be the desired
signal component of the first local microphone signal, i.e.,
$d_{k}[\nu,i]=\mathbf{e}_{d_{k}}^{T}\mathbf{s}[\nu,i]$, where
$\mathbf{e}_{d_{k}}$ selects the appropriate channel of $\mathbf{s}[\nu,i]$.
An optimal filter $\bar{\mathbf{w}}_{k}[\nu,i]$ can be obtained by minimising
the mean squared error (MSE) between the desired signal and the filtered
microphone signals:

$\bar{\mathbf{w}}_{k}[\nu,i+1]=\underset{\mathbf{w}_{k}[\nu]}{\operatorname*{argmin}}\>E\left\\{\left|d_{k}[\nu,i]-\mathbf{w}_{k}^{H}[\nu]\mathbf{y}[\nu,i]\right|^{2}\right\\},$
(2)

where $\left(\cdot\right)^{H}$ denotes complex conjugation and $E\\{\cdot\\}$
the expected value operation. The closed-form solution of (2) is the MWF:

$\bar{\mathbf{w}}_{k}[\nu,i+1]=\left(\bar{\mathbf{R}}_{\mathbf{y}\mathbf{y}}[\nu,i]\right)^{-1}\bar{\mathbf{R}}_{\mathbf{s}\mathbf{s}}[\nu,i]\mathbf{e}_{d_{k}},$
(3)

where
$\bar{\mathbf{R}}_{\mathbf{y}\mathbf{y}}[\nu,i]=E\\{\mathbf{y}[\nu,i]\mathbf{y}^{H}[\nu,i]\\}$
is the network-wide microphone signal covariance matrix and
$\bar{\mathbf{R}}_{\mathbf{s}\mathbf{s}}[\nu,i]=\bar{\mathbf{R}}_{\mathbf{y}\mathbf{y}}[\nu,i]-\bar{\mathbf{R}}_{\mathbf{n}\mathbf{n}}[\nu,i]$,
where
$\bar{\mathbf{R}}_{\mathbf{n}\mathbf{n}}[\nu,i]=E\\{\mathbf{n}[\nu,i]\mathbf{n}^{H}[\nu,i]\\}$
is the network-wide noise-only covariance matrix. Assuming short-term
stationarity of the signals, the covariance matrices can be estimated by
averaging over observations of $\mathbf{y}[\nu,i]\mathbf{y}^{H}[\nu,i]$. In a
speech enhancement scenario with stationary noise, the on-off behaviour of the
desired signal can be exploited via a voice activity detector (VAD) [15, 16]
to isolate noise-only observations of $\mathbf{y}[\nu,i]$. The estimation of
$\bar{\mathbf{R}}_{\mathbf{y}\mathbf{y}}[\nu,i]$ and
$\bar{\mathbf{R}}_{\mathbf{n}\mathbf{n}}[\nu,i]$ can then be performed via
exponential averaging:

$\begin{split}\text{VAD = 1 :}\\\
\hat{\mathbf{R}}_{\mathbf{y}\mathbf{y}}[\nu,i]&=\beta\hat{\mathbf{R}}_{\mathbf{y}\mathbf{y}}[\nu,i-1]+(1-\beta)\mathbf{y}[\nu,i]\mathbf{y}^{H}[\nu,i],\\\
\text{VAD = 0 :}\\\
\hat{\mathbf{R}}_{\mathbf{n}\mathbf{n}}[\nu,i]&=\beta\hat{\mathbf{R}}_{\mathbf{n}\mathbf{n}}[\nu,i-1]+(1-\beta)\mathbf{y}[\nu,i]\mathbf{y}^{H}[\nu,i],\end{split}$
(4)

where the real-valued number $\beta$ acts as a forgetting factor,
$0\ll\beta\leq 1$, typically chosen close to 1 to preserve spatial coherence
between microphone signals [1].

In the presence of a single desired signal source, the signal model implies
that $\bar{\mathbf{R}}_{\mathbf{s}\mathbf{s}}[\nu,i]$ should be a rank-1
matrix [17]. However, the estimated
$\hat{\mathbf{R}}_{\mathbf{s}\mathbf{s}}[\nu,i]=\hat{\mathbf{R}}_{\mathbf{y}\mathbf{y}}[\nu,i]-\hat{\mathbf{R}}_{\mathbf{n}\mathbf{n}}[\nu,i]$
generally has a rank larger than 1. A GEVD-based approach was proposed in [17]
to obtain a rank-1 approximation of
$\hat{\mathbf{R}}_{\mathbf{s}\mathbf{s}}[\nu,i]$. The GEVD of the matrix
pencil
$\\{\hat{\mathbf{R}}_{\mathbf{y}\mathbf{y}}[\nu,i],\hat{\mathbf{R}}_{\mathbf{n}\mathbf{n}}[\nu,i]\\}$
yields:

$\begin{split}\hat{\mathbf{R}}_{\mathbf{y}\mathbf{y}}[\nu,i]&=\hat{\mathbf{Q}}[\nu,i]\hat{\mathbf{\Sigma}}[\nu,i]\hat{\mathbf{Q}}^{H}[\nu,i],\\\
\hat{\mathbf{R}}_{\mathbf{n}\mathbf{n}}[\nu,i]&=\hat{\mathbf{Q}}[\nu,i]\hat{\mathbf{Q}}^{H}[\nu,i],\end{split}$
(5)

with $\hat{\mathbf{Q}}[\nu,i]$ an $M\times M$ matrix of which the columns are
the generalised eigenvectors (GEVCs) and $\hat{\mathbf{\Sigma}}[\nu,i]$ is a
diagonal matrix of which the diagonal elements are the corresponding
generalised eigenvalues (GEVLs). The GEVLs in $\hat{\mathbf{\Sigma}}[\nu,i]$
are assumed to be ordered by decreasing magnitude. A rank-1 approximation of
$\bar{\mathbf{R}}_{\mathbf{s}\mathbf{s}}[\nu,i]$ can then be obtained by using
(5) and nullifying the $M-1$ smallest GEVLs. Substituting into (3) then leads
to the GEVD-MWF:

$\hat{\mathbf{w}}_{k}[\nu,i+1]=\hat{\mathbf{Q}}^{-H}[\nu,i]\mathbf{\Lambda}[\nu,i]\hat{\mathbf{Q}}^{H}[\nu,i]\mathbf{e}_{d_{k}},$
(6)

with
$\mathbf{\Lambda}[\nu,i]=\mathrm{diag}\\{1-1/\hat{\sigma}_{1}[\nu,i],0,\dots,0\\}$,
where $\mathrm{diag}\\{\cdot\\}$ transforms a vector into a diagonal matrix
and $\hat{\sigma}_{1}[\nu,i]$ is the largest GEVL. Finally, the desired signal
at frequency bin $\nu$ and frame $i$ is estimated as
$\hat{d}_{k}[\nu,i]=\hat{\mathbf{w}}_{k}^{H}[\nu,i+1]\mathbf{y}[\nu,i]$.

##  3 The DANSE algorithm

The DANSE algorithm [2] provides a distributed implementation of the MWF
described in Section 2 as an adaptive algorithm where nodes iteratively update
their local filter estimates. Different node-updating schemes exist: (i)
sequential updating [2], (ii) simultaneous updating [3], and (iii)
asynchronous updating [3]. Strategies (i) and (ii) rely on a network-wide
update protocol that coordinates the updates, unlike strategy (iii) [3]. Since
the presence of unknown SROs between nodes challenges the deployment of a
coordination protocol, asynchronous updating is assumed in the following. To
avoid limit cycles due to asynchronous updating [3], relaxed filter updates
can be performed [18]. Computational delays due to data transmission,
reception, and processing are assumed to be negligible in this paper.

As described in [14], the DANSE algorithm can be implemented using weighted
overlap-add (WOLA) processing to efficiently perform short-time Fourier
analysis and synthesis [19]. Using WOLA, time-domain microphone signals are
processed on a frame-by-frame basis. WOLA analysis consists of applying an
$N$-point DFT to a windowed frame of a time-domain signal, with the frame size
equal to the DFT size, effectively transforming the time-domain signal frame
into an STFT-domain signal frame. In DANSE, all filtering can be conducted in
the STFT domain, resulting in a lower computational complexity as compared to
a time-domain implementation [14]. Each new WOLA frame then corresponds to a
new DANSE iteration where the nodes update their filter estimate. As a WOLA
implementation of DANSE is assumed in the following, the variable $i$
simultaneously denotes the STFT frame index as well as the DANSE iteration
index, i.e., nodes update their filter estimates at each new $i$.

Although a variety of network topologies can exist, a fully connected WASN is
assumed in this paper. All the nodes that can communicate with node $k$ are
grouped in the set $\mathcal{K}_{k}=\mathcal{K}\backslash\\{k\\}$. The DANSE
algorithm in a fully connected WASN operates in two main stages: signals
fusion and broadcasting on the one hand, and filters updates on the other hand
[2]. At each frame $i$, each node $k$ fuses its $M_{k}$ local microphone
signals into a single-channel signal $z_{k}[\nu,i]$ $\forall\>\nu\in\\{1\dots
N\\}$ using a fusion vector $\mathbf{p}_{k}[\nu,i]$ before broadcasting it to
the other nodes, which reduces the per-node communication cost by a factor
$M_{k}$ (i.e., a factor $M/K$ over the entire network) compared to the
centralised MWF of Section 2. An appropriate definition of
$\mathbf{p}_{k}[\nu,i]$ guarantees convergence of the DANSE algorithm to the
centralised MWF solution [2]. In the WOLA implementation, a time-domain fused
signal denoted by $\dot{z}_{k}[n]$ is obtained via WOLA synthesis (inverse DFT
followed by windowing) and overlap-add of the fused signal frames
$z_{k}[\nu,i]$, where $n$ denotes the sample index. The time-domain signal
$\dot{z}_{k}[n]$ is then broadcast to other nodes, as summarised in Algorithm
1.

The STFT-domain signals available at node $k$ at iteration $i$ are grouped
into the vector:

$\tilde{\mathbf{y}}_{k}[\nu,i]=\begin{bmatrix}\mathbf{y}_{k}^{T}[\nu,i]\>|\>\mathbf{z}_{-k}^{T}[\nu,i]\end{bmatrix}^{T}$
(7)

where $\mathbf{y}_{k}[\nu,i]$ contains the local microphone signals and
$\mathbf{z}_{-k}[\nu,i]$ is a stacked version of all the microphone signals
received from other nodes. As in the centralised case (cfr. (1)),
$\tilde{\mathbf{y}}_{k}[\nu,i]$ can be written as a sum of a desired signal
component $\tilde{\mathbf{s}}_{k}[\nu,i]$ and a noise component
$\tilde{\mathbf{n}}_{k}[\nu,i]$. Node $k$ aims to compute the $i$-th STFT
frame of its desired signal estimate $\hat{d}_{k}[\nu,i]$ via multichannel
filtering of $\tilde{\mathbf{y}}_{k}[\nu,i]$. The filter at node $k$ is
denoted by
$\tilde{\mathbf{w}}_{k}[\nu,i]=[\mathbf{w}_{kk}^{T}[\nu,i]\>|\>\mathbf{g}_{k-k}^{T}[\nu,i]]^{T}$,
where $\mathbf{w}_{kk}[\nu,i]$ is applied to the local microphone signals
$\mathbf{y}_{k}[\nu,i]$ and $\mathbf{g}_{k-k}[\nu,i]$ is applied to the fused
microphone signals $\mathbf{z}_{-k}[\nu,i]$. The filter at node $k$ at
iteration $i+1$ is obtained by minimising the MSE between the desired signal
and its estimate:

$\tilde{\mathbf{w}}_{k}[\nu,i+1]=\underset{\mathbf{w}_{k}[\nu]}{\operatorname*{argmin}}\>E\left\\{\left|d_{k}[\nu,i]-\mathbf{w}_{k}^{H}[\nu]\tilde{\mathbf{y}}_{k}[\nu,i]\right|^{2}\right\\},$
(8)

and the $i$-th STFT frame of desired signal estimate is then
$\hat{d}_{k}[\nu,i]=\tilde{\mathbf{w}}_{k}^{H}[\nu,i+1]\tilde{\mathbf{y}}_{k}[\nu,i]$.
Equation (8) has the same structure as (2), be it with a different definition
of the filter $\mathbf{w}_{k}[\nu]$ and input vector
$\tilde{\mathbf{y}}_{k}[\nu,i]$, hence its solution again corresponds to an
MWF. With the covariance matrices
$\tilde{\mathbf{R}}_{\mathbf{y}_{k}\mathbf{y}_{k}}[\nu,i]$ and
$\tilde{\mathbf{R}}_{\mathbf{n}_{k}\mathbf{n}_{k}}[\nu,i]$ defined and
estimated per node instead of centrally as in (4), a GEVD is applied to the
matrix pencil
$\\{\tilde{\mathbf{R}}_{\mathbf{y}_{k}\mathbf{y}_{k}}[\nu,i],\tilde{\mathbf{R}}_{\mathbf{n}_{k}\mathbf{n}_{k}}[\nu,i]\\}$
and, similarly to (6), the filter is computed as:

$\tilde{\mathbf{w}}_{k}[\nu,i+1]=\tilde{\mathbf{Q}}_{k}^{-H}[\nu,i]\tilde{\mathbf{\Lambda}}_{k}[\nu,i]\tilde{\mathbf{Q}}_{k}^{H}[\nu,i]\mathbf{e}_{d_{k}},$
(9)

where, at frequency $\nu$ and iteration $i$, $\tilde{\mathbf{Q}}_{k}[\nu,i]$
is an $(M_{k}+K-1)\times(M_{k}+K-1)$ matrix of which the columns are the
GEVCs,
$\tilde{\mathbf{\Lambda}}_{k}[\nu,i]=\mathrm{diag}\left\\{1-1/\tilde{\sigma}_{k1}[\nu,i],0,\dots,0\right\\}$
and $\tilde{\sigma}_{k1}[\nu,i]$ is the largest GEVL. Finally, to ensure
convergence of $\mathbf{w}_{kk}[\nu,i]$ towards the corresponding elements of
the centralised MWF, the fusion vector at iteration $i$ is defined as
$\mathbf{p}_{k}[\nu,i]=\mathbf{w}_{kk}[\nu,i]$, such that:

$z_{k}[\nu,i]=\mathbf{w}_{kk}^{H}[\nu,i]\mathbf{y}_{k}[\nu,i].$ (10)

The WOLA implementation of the DANSE algorithm is summarised in Algorithm 1,
where the $M_{k}$ local time-domain microphone signals at node $k$ are denoted
by $\dot{\mathbf{y}}_{k}[n]$. The time-domain signal obtained after WOLA
synthesis and overlap-add of consecutive frames $\hat{d}_{k}[\nu,i]$ is
denoted by $\dot{\hat{d}}_{k}[n]$. The time-domain fused signals are grouped
in the vector
$\dot{\mathbf{z}}_{-k}[n]=\begin{bmatrix}\dot{z}_{1}[n]\>\dots\>\dot{z}_{k-1}[n]\>\dot{z}_{k+1}[n]\>\dots\>\dot{z}_{K}[n]\end{bmatrix}^{T}$.
The WOLA window shift, corresponding to the number of new samples recorded
between two consecutive DANSE iterations, is denoted by $N_{\mathrm{s}}$.

Algorithm 1 WOLA-based DANSE in a synchronised, fully connected WASN (50%
window overlap).

1: Initialise $\tilde{\mathbf{w}}_{k}[\nu,0]$ $\forall\>k\in\mathcal{K}$;

2: Each node $k\in\mathcal{K}$ performs, starting simultaneously:

3: for $i=1,2,3,\dots$ do

4: Record $N_{\mathrm{s}}$ new samples of $\dot{\mathbf{y}}_{k}[n]$ since
$i-1$;

5: WOLA analysis on $N$ most recent $\dot{\mathbf{y}}_{k}[n]$ samples to
obtain $\mathbf{y}_{k}[\nu,i]$;

6: Perform signal fusion via (10) to obtain $z_{k}[\nu,i]$;

7: WOLA synthesis on $z_{k}[\nu,i]$ and overlap-add with previous frame to
obtain $N_{\mathrm{s}}$ new $\dot{z}_{k}[n]$ samples;

8: Transmit $N_{\mathrm{s}}$ most recent $\dot{z}_{k}[n]$ samples to
$\mathcal{K}_{k}$;

9: Build $\dot{\mathbf{z}}_{-k}[n]$ from samples received from
$\mathcal{K}_{k}$;

10: WOLA analysis on $N$ most recent $\dot{\mathbf{z}}_{-k}[n]$ samples to
obtain $\mathbf{z}_{-k}[\nu,i]$;

11: Build
$\tilde{\mathbf{y}}_{k}[\nu,i]=[\mathbf{y}_{k}^{T}[\nu,i]\>|\>\mathbf{z}_{-k}^{T}[\nu,i]]^{T}$;

12: Compute $\tilde{\mathbf{R}}_{\mathbf{y}_{k}\mathbf{y}_{k}}[\nu,i]$ and
$\tilde{\mathbf{R}}_{\mathbf{n}_{k}\mathbf{n}_{k}}[\nu,i]$ via (4);

13: Compute the filter estimates $\tilde{\mathbf{w}}_{k}[\nu,i+1]$ via (9);

14: Compute new frame $\hat{d}_{k}[\nu,i]$;

15: WOLA synthesis on $\hat{d}_{k}[\nu,i]$ and overlap-add with the previous
frame to build $\dot{\hat{d}}_{k}[n]$.

16: end for

In the presence of SROs, the time misalignments between the local microphone
signals and the fused microphone signals from other nodes lead to incorrect
covariance matrix updates which, in turn, inhibit the computation of useful
filter estimates via (9). In the following section, we propose a method for
SRO estimation and compensation applicable to the WOLA implementation of the
DANSE algorithm.

##  4 SRO estimation and compensation

The SRO between node $k$ and $q$ is denoted by $\varepsilon_{kq}$ such that
$f_{\mathrm{s},q}=f_{\mathrm{s},k}(1+\varepsilon_{kq})$, where
$f_{\mathrm{s},k}$ and $f_{\mathrm{s},q}$ are the sampling rate of node $k$
and $q$, respectively. In the following, it is assumed that the SROs are time-
invariant and that signals recorded by the same node are synchronised.
Although a fully connected WASN is assumed, this SRO estimation and
compensation method can be also adopted in other network topologies.

###  4.1 Coherence-drift-based SRO estimation

In order to allow any node $k\in\mathcal{K}$ to blindly estimate the SROs
$\\{\varepsilon_{kq}\\}_{q\in\mathcal{K}_{k}}$ using the signals it can access
in the DANSE algorithm, we use a coherence-drift method based on principles
introduced in [12] and [13]. At frame $i$ and at node $k$, considering one
other node $q\in\mathcal{K}_{k}$, the available STFT-domain signals are (i)
the local microphone signals $\mathbf{y}_{k}[\nu,i]$ and (ii) the received
fused signal $z_{q}[\nu,i]$. The first local microphone signal
$y_{k,1}[\nu,i]$ is used in the following (w.l.o.g.).

The sampling rate mismatch can simply be approximated in the STFT-domain via
the linear phase drift (LPD) model [20, 21] at any frequency bin
$\nu\in\\{1,\dots,N\\}$, i.e.:

$\check{z}_{q}[\nu,i]\approx
z_{q}[\nu,i]\cdot\mathrm{exp}\left(\mathrm{j}\frac{2\pi}{N}\nu\varepsilon_{kq}[i]N_{\mathrm{c}}[i]\right),$
(11)

where $\check{z}_{q}[\nu,i]$ is the $\varepsilon_{kq}[i]$-compensated version
of $z_{q}[\nu,i]$ (synchronised with $y_{k,1}[\nu,i]$) and $N_{\mathrm{c}}[i]$
is the central sample index of frame $i$. The product
$\tau_{kq}[i]=\varepsilon_{kq}[i]N_{\mathrm{c}}[i]$ is the average accumulated
time-drift between $z_{q}[\nu,i]$ and $y_{k,1}[\nu,i]$. The LPD model relies
on the assumption that the SRO-induced time drift is constant within one
frame, implying that the model best approximates the effect of SROs for small
$\varepsilon_{kq}[i]$.

Based on (11), $\varepsilon_{kq}[i]$ can be estimated by node $k$ as follows.
First, we define the instantaneous estimate of the cross-power spectral
density (PSD) $\Psi_{kq}[\nu,i]$ as $\Psi_{kq}[\nu,i]=y_{k,1}[\nu,i]\cdot
z_{q}^{\ast}[\nu,i]$, where $\cdot^{\ast}$ denotes complex conjugation.
Similarly, the instantaneous auto-PSD estimates are defined as
$\Psi_{kk}[\nu,i]=|y_{k,1}[\nu,i]|^{2}$ and
$\tilde{\Psi}_{qq}[\nu,i]=|z_{q}[\nu,i]|^{2}$. An instantaneous estimate of
the coherence between $y_{k1}[\nu,i]$ and $z_{q}[\nu,i]$ can then be obtained
as:

$\Gamma_{kq}[\nu,i]=\frac{\Psi_{kq}[\nu,i]}{\sqrt{\Psi_{kk}[\nu,i]\cdot\tilde{\Psi}_{qq}[\nu,i]}}.$
(12)

The SRO can now be estimated by defining the product $P_{\Gamma,kq}[\nu,i]$
between the instantaneous coherence estimate at frame $i$ and at frame
$i-l_{\mathrm{d}}$ as:

$P_{\Gamma,kq}[\nu,i]=\Gamma_{kq}[\nu,i]\cdot\Gamma_{kq}^{\ast}[\nu,i-l_{\mathrm{d}}].$
(13)

Based on the LPD model and assuming static sources, it can be shown that an
SRO estimate $\hat{\varepsilon}_{kq}[i]$ proportional to the phase of
$P_{\Gamma,kq}[\nu,i]$ [12, 13] is obtained as:

$\angle\left\\{P_{\Gamma,kq}[\nu,i]\right\\}=\frac{2\pi}{N}\nu
l_{\mathrm{d}}N_{\mathrm{s}}\hat{\varepsilon}_{kq}[i],$ (14)

where $\angle\\{\cdot\\}$ denotes the phase. Increasing the value of
$l_{\mathrm{d}}$ is equivalent to estimating the average SRO over a longer
period of time, setting a trade-off between robust estimation and the ability
to track time-varying SROs. Since fixed SROs are considered here,
$l_{\mathrm{d}}$ may be safely set to a relatively large value, bearing in
mind that SRO estimation can only begin after $l_{\mathrm{d}}$ frames.
Temporal averaging can be applied before computing the phase to smoothen the
estimation:

$\bar{P}_{\Gamma,kq}[\nu,i]=\alpha\bar{P}_{\Gamma,kq}[\nu,i-1]+(1-\alpha)P_{\Gamma,kq}[\nu,i],$
(15)

where $\alpha$ is a scalar, $0\ll\alpha\leq 1$, set close to 1.

Since (14) and (15) are defined for all frequency bins $\nu$, the SRO can be
estimated, for example, as the least squares (LS) solution over all relevant
frequency bins [22]. Since this LS solution is, however, prone to inaccuracies
due to the periodicity of the phase, it has been proposed in [13] to interpret
$\bar{P}_{\Gamma,kq}[\nu,i]$ as a generalised cross-PSD. The integer time lag
$\lambda_{\mathrm{max}}[i]$ that maximises the absolute value of the
generalised cross-correlation
$\bar{p}_{\Gamma,kq}^{i}[\lambda]=\mathcal{F}^{-1}\\{\bar{P}_{\Gamma,kq}[\nu,i]\\}$,
with $\mathcal{F}^{-1}\\{\cdot\\}$ denoting the inverse DFT, can then be used
to estimate the SRO as:

$\hat{\varepsilon}_{kq}[i]=-\frac{\lambda_{\mathrm{max}}[i]}{l_{\mathrm{d}}N_{\mathrm{s}}}=-\frac{1}{l_{\mathrm{d}}N_{\mathrm{s}}}\cdot\underset{\lambda}{\operatorname*{argmax}}|\bar{p}^{i}_{\Gamma,kq}[\lambda]|.$
(16)

Higher SRO estimation accuracy can be obtained by determining the non-integer
value $\lambda[i]$ that maximises $|p_{\Gamma,kq}^{i}[\lambda]|$, via an
interpolation method such as a golden section search in the interval
$[\lambda_{\mathrm{max}}[i]-0.5,\lambda_{\mathrm{max}}[i]+0.5]$, as proposed
in [13], and substituting $\lambda_{\mathrm{max}}[i]$ by $\lambda[i]$ in (16).

###  4.2 SRO compensation and full-sample drifts

The SRO estimates obtained via the method described in Section 4.1 are now
used to perform SRO compensation on the elements of $\mathbf{z}_{-k}[\nu,i]$
before updating $\tilde{\mathbf{w}}_{k}[\nu,i]$ as described in Section 3.
Using the LPD model, SRO compensation can be performed at any node $k$ based
on $\\{\hat{\varepsilon}_{kq}[i]\\}_{q\in\mathcal{K}_{k}}$ by applying the
appropriate phase shift to each element of $\mathbf{z}_{-k}[\nu,i]$ as:

$\check{z}_{q}[\nu,i]=z_{q}[\nu,i]\cdot\mathrm{exp}\left(\mathrm{j}\frac{2\pi}{N}\nu\hat{\tau}_{kq}[i]\right)\>\forall\>q\in\mathcal{K}_{k},$
(17)

where
$\hat{\tau}_{kq}[i]=N_{\mathrm{s}}\sum_{\iota=0}^{i}\hat{\varepsilon}_{kq}[\iota]$
is the estimated average accumulated time-drift between $z_{q}[\nu,i]$ and
$y_{k1}[\nu,i]$ (cfr. (11)).

An important aspect comes into play
$t^{\mathrm{FSD}}_{kq}=1/(f_{\mathrm{s},k}|\varepsilon_{kq}|)$ seconds after
the simultaneous initialisation of the WASN, namely when the accumulated SRO-
induced time drift $\tau_{kq}[i]$ between node $k$ and node $q$ becomes
greater than one sample. Such event is referred to in the following as a full-
sample drift (FSD). At that time, if $\varepsilon_{kq}>0$, the growing time
drift between node $k$ and node $q$ leads to a situation where node $q$ has
recorded one more sample than node $k$, as depicted in Figure 1. Conversely,
if $\varepsilon_{kq}<0$, node $q$ has recorded one less sample than node $k$.

Figure 1: Schematic representation of a full-sample drift (indicated by the
circle) generated by an SRO $\varepsilon_{kq}>0$ between node $k$ and
$q\in\mathcal{K}_{k}$.

When correctly detected, an FSD can be compensated for by applying a
corrective phase shift $\phi_{kq}^{\mathrm{FSD}}[\nu,i]$ to $z_{q}[\nu,i]$ as:

$\phi_{kq}^{\mathrm{FSD}}[\nu,i]=\begin{cases}\mathrm{exp}\left(-\mathrm{j}\frac{2\pi}{N}\nu\right)&\text{if
one more sample at $q$,}\\\
\mathrm{exp}\left(\mathrm{j}\frac{2\pi}{N}\nu\right)&\text{if one less sample
at $q$,}\\\ 1&\text{otherwise.}\end{cases}$ (18)

The SRO estimation itself can be biased by the presence of one or more FSDs
between frame $i-l_{\mathrm{d}}$ and frame $i$. These can be accounted for by
multiplying $P_{\Gamma,kq}[\nu,i]$ by the accumulated FSD phase shift:

$\phi_{kq}^{\mathrm{ac}}[\nu,i]=\prod_{\iota=i-l_{\mathrm{d}}}^{i}\phi^{\mathrm{FSD}}_{kq}[\nu,\iota].$
(19)

The rest of the SRO estimation process remains unchanged, following (15) and
(16).

However, the accumulated effect of FSDs becomes particularly problematic when
considering the WOLA implementation of DANSE [14], where a fused time-domain
signal $\dot{z}_{q}[n]$ is transmitted in frames of $N_{\mathrm{s}}$ samples
from node $q$ to node $k$ (cfr. Algorithm 1). For clarity of exposition, we
assume an even DFT size $N$ and a 50% WOLA window shift such that
$N_{\mathrm{s}}=N/2$. A problematic phenomenon referred to as full-frame drift
(FFD) occurs when $N_{\mathrm{s}}$ uncompensated FSDs accumulate. If
$\varepsilon_{kq}>0$ (resp. $\varepsilon_{kq}<0$) and after
$t^{\mathrm{FFD}}_{kq}=N_{\mathrm{s}}t^{\mathrm{FSD}}_{kq}$ seconds, node $q$
has recorded $N_{\mathrm{s}}$ more (resp. less) samples than node $k$ since
the synchronous initialisation of both nodes. At that time, node $q$ has thus
transmitted two (resp. no) new $\dot{z}_{q}[n]$ frames since the last update
of node $k$ (see circles on Figure 2). Consequently, to perform its next
update, node $k$ skips (resp. duplicates) one $z_{q}[\nu,i]$ frame.

Figure 2: Schematic representation of a full-frame drift (highlighted by
circles), with $\varepsilon_{kq}>0$ and 50% WOLA window shift.

An FFD cannot be compensated for via a phase shift based on (18) if
$N_{\mathrm{s}}$ is close to $N$. For instance, with 50% WOLA window shift,
the corrective phase shift of (18) needed to compensate for $N_{\mathrm{s}}$
FSDs at once simplifies to $\mathrm{exp}(\pm\mathrm{j}\frac{2\pi}{N}\nu
N_{\mathrm{s}})=\pm 1$ $\forall\>\nu$. Even if FFDs compensation were
possible, before an FFD occurs node $k$ receives a single $N_{\mathrm{s}}$
samples-long frame of $\dot{z}_{q}[n]$ between two consecutive filter updates,
as in Figure 2. Node $k$ is, therefore, unable to detect FSDs by comparing the
number of local $\dot{\mathbf{y}}_{k}[n]$ samples with the number of received
$\dot{z}_{q}[n]$ samples since its previous update (as both are equal to
$N_{\mathrm{s}}$). The uncompensated growing drift between elements of
$\tilde{\mathbf{y}}_{k}[\nu,i]$ then leads to increasingly erroneous updates
of the covariance matrices.

If disregarded, FFDs can significantly perturb the convergence of DANSE as
well as the SRO estimation process. The detection of FSDs within the WOLA
implementation of DANSE is discussed in the following section.

###  4.3 Full-sample drift detection

In order to enable detection and compensation of FSDs within the WOLA
implementation of DANSE, we introduce a modification of the DANSE fusion and
broadcasting mechanism to allow per-sample transmission of fused signals
between nodes, while retaining WOLA frame-by-frame processing for the
computationally costly steps of GEVD-based filter update and desired signal
estimate computation. In principle, using this per-sample transmission, node
$k$ can easily detect FSDs for the $i$-th filter update by comparing the
number of local $\dot{\mathbf{y}}_{k}[n]$ samples with the number of received
$\dot{z}_{q}[n]$ samples from node $q$ since its previous update, then
compensate for them via the corrective phase shifts of (18).

We propose to approximate the WOLA filtering process (analysis, STFT-domain
filtering, and synthesis) by its so-called distortion function $T(\zeta)$
[23], where $\zeta$ is the $\mathcal{Z}$-transform variable. This function
relates the output of the WOLA filterbank to its input when no decimation and
expansion is performed, i.e., using maximal window overlap. At frame $i$, the
distortion function $T^{i}_{q,m}(\zeta)$ corresponding to the $m$-th
microphone of node $q$ can be obtained as:

$T_{q,m}^{i}(\zeta)=\frac{1}{N_{\mathrm{s}}}\begin{bmatrix}\zeta^{1-N}\,\dots\,1\end{bmatrix}\mathbf{D}_{q,m}^{i}\begin{bmatrix}1\,\dots\,\zeta^{1-N}\end{bmatrix}^{T},$
(20)

with
$\mathbf{D}_{q,m}^{i}=\mathbf{H}_{\mathrm{s}}\cdot\mathbf{F}^{-1}\cdot\mathrm{diag}\\{\mathbf{w}_{qq,m}[i]\\}\cdot\mathbf{F}\cdot\mathbf{H}_{\mathrm{a}}$,
where $\mathbf{F}^{-1}$ and $\mathbf{F}$ are the inverse DFT and DFT matrix,
respectively,
$\mathbf{w}_{qq,m}[i]=\begin{bmatrix}w_{qq,m}[1,i],\dots,w_{qq,m}[N,i]\end{bmatrix}^{T}$
denotes the local filter coefficients at frame $i$ for the $m$-th microphone
of node $q$ with all frequency bins stacked into one vector,
$\mathbf{H}_{\mathrm{s}}=\mathrm{diag}\\{\mathrm{flip}\\{\mathbf{h}_{\mathrm{s}}\\}\\}$,
and $\mathbf{H}_{\mathrm{a}}=\mathrm{diag}\\{\mathbf{h}_{\mathrm{a}}\\}$,
respectively, where $\mathbf{h}_{\mathrm{s}}$ and $\mathbf{h}_{\mathrm{a}}$
denote the WOLA synthesis and analysis time-domain windows, respectively, and
$\mathrm{flip}\\{\cdot\\}$ reverses the order of the elements of a vector.

The time-domain equivalent of the distortion function $T_{q,m}^{i}(\zeta)$ in
(20) is a $(2N-1)$-tap impulse response denoted by $\mathbf{t}_{q,m}^{i}$.
From (20), it can be seen that each element of $\mathbf{t}_{q,m}^{i}$ is
obtained by summing over the corresponding diagonal of the matrix
$\mathbf{D}_{q,m}^{i}$. The complete WOLA analysis and synthesis process can
then be approximated by a convolution with $\mathbf{t}_{q,m}^{i}$. This means
that the $n$-th sample of the time-domain fused signal $\dot{z}_{q}[n]$ can be
obtained as:

$\dot{z}_{q}[n]=\sum_{m=1}^{M_{q}}\left(\dot{\mathbf{y}}_{q,m}^{(n)}\ast\mathbf{t}_{q,m}^{i}\right)[n+2N-1]$
(21)

where the time-domain vector $\dot{\mathbf{y}}_{q,m}^{(n)}$ contains the most
recent $N$ samples recorded by the $m$-th microphone of node $q$ and
$(\mathbf{a}\ast\mathbf{b})[c]$ denotes the $c$-th sample of the convolution
between time-domain signals $\mathbf{a}$ and $\mathbf{b}$. Note that the
distortion function does not need be computed at every frame $i$, especially
once the filters have converged after several DANSE iterations. The iteration
indices at which the distortion function is updated with the most recent
filter $\mathbf{w}_{qq}[\nu,i]$ are grouped in the set
$\mathcal{I}_{\mathrm{T}}$.

Although the proposed $T(\zeta)$-approximation introduces the same $N-1$
samples input-output delay as the standard WOLA implementation of DANSE [14],
it has the advantage to circumvent the $N_{\mathrm{s}}$ samples delay
introduced by frame-by-frame processing [19] since no downsampling is
performed. Additionally, the use of per-sample broadcasting reduces the amount
of transmitted data as each compressed signal sample is transmitted only once.
This differs from the usual WOLA scheme where the overlap-add operation
necessitates the transmission of $N_{\mathrm{s}}$ additional data points per
$N$-samples block of compressed signal (as in Algorithm 1).

Using the $T(\zeta)$-approximation, any node is able to broadcast its fused
signal on a per-sample basis. This modification of the DANSE algorithm,
although coming at the expense of some additional computational complexity
with respect to a purely WOLA-based implementation, enables the deployment of
DANSE in asynchronous WASNs where FSDs can be detected as soon as they occur.
An overview of the DANSE algorithm with per-sample fused signal broadcasting
using the $T(\zeta)$-approximation is provided in Figure 3.

Figure 3: Proposed WOLA-based DANSE processing at node $k\in\mathcal{K}$ with
per-sample fused signal broadcasting. [i] WOLA analysis applied to local
microphone signals $\dot{\mathbf{y}}_{k}[n]$ and fused signals from other
nodes $\dot{\mathbf{z}}_{-k}[n]$. [ii] Covariance matrix update and [iii]
computation of filter $\tilde{\mathbf{w}}_{k}[\nu,i+1]$. [iv] Computation of
new desired signal estimate frame $\hat{d}_{k}[\nu,i+1]$, followed by WOLA
synthesis and overlap-add. [v] Computation of distortion functions
$\\{T_{k,m}^{i}(\zeta)\\}_{m=1}^{M_{k}}$ from filter $\mathbf{w}_{kk}[\nu,i]$.
[vi] Computation of new $\dot{z}_{k}[n]$ samples and per-sample broadcasting
to other nodes.

###  4.4 Complete system

Figure 4: SRO estimation and compensation block-scheme at node
$k\in\mathcal{K}$, including full-sample drift detection (“flag”) and
compensation.

As SRO estimation is necessary for SRO compensation, both should be performed
in parallel. An open-loop strategy is proposed, as depicted in Figure 4, which
consists of three parts: SRO estimation, FSD detection, and SRO compensation.
First, the SRO-uncompensated fused signal $z_{q}[\nu,i]$ is used to estimate
$\hat{\varepsilon}_{kq}[i]$. Every time an FSD is detected, a flag is raised
and the FSD phase shift of (18) is included when performing SRO estimation and
compensation, leading to the signal $\check{z}_{q}[\nu,i]$, which is used to
update the DANSE filter.

Algorithm 2 provides a complete description of WOLA-based DANSE with SRO
estimation and compensation, including the $T(\zeta)$-approximation for FSD
detection. There, the STFT-domain SRO-compensated fused signals vector is
denoted by $\check{\mathbf{z}}_{-k}[\nu,i]$ and the SRO-compensated version of
$\tilde{\mathbf{y}}_{k}[\nu,i]$ becomes
$\check{\mathbf{y}}_{k}[\nu,i]=[\mathbf{y}_{k}^{T}[\nu,i]\>|\>\check{\mathbf{z}}_{-k}^{T}[\nu,i]]^{T}$.
The estimates of
$E\\{\check{\mathbf{y}}_{k}[\nu,i]\check{\mathbf{y}}_{k}^{H}[\nu,i]\\}$ and
$E\\{\check{\mathbf{n}}_{k}[\nu,i]\check{\mathbf{n}}_{k}^{H}[\nu,i]\\}$ are
denoted by $\check{\mathbf{R}}_{\mathbf{y}_{k}\mathbf{y}_{k}}[\nu,i]$ and
$\check{\mathbf{R}}_{\mathbf{n}_{k}\mathbf{n}_{k}}[\nu,i]$, respectively. The
filter estimate after SRO compensation is finally obtained similarly to (9),
i.e., performing a GEVD on the matrix pencil
$\\{\check{\mathbf{R}}_{\mathbf{y}_{k}\mathbf{y}_{k}}[\nu,i],\check{\mathbf{R}}_{\mathbf{n}_{k}\mathbf{n}_{k}}[\nu,i]\\}$,
and is denoted by $\check{\mathbf{w}}_{k}[\nu,i+1]$.

Algorithm 2 WOLA-based DANSE with SRO compensation in a fully connected
heterogeneous WASN.

1: Initialise $\tilde{\mathbf{w}}_{k}[\nu,0]$
$\forall\>(k,\nu)\in\mathcal{K}\times\\{1,\dots,N\\}$;

2: Each node $k\in\mathcal{K}$ performs, starting simultaneously:

3: for every new locally recorded sample $\dot{\mathbf{y}}[n]$ do

4: Compute $\dot{z}_{k}[n]$ via (21) and transmit to nodes in
$\mathcal{K}_{k}$.

5: end for

6: for $i=1,2,3,\dots$ do

7: if $i\in\mathcal{I}_{\mathrm{T}}$ then

8: Update $\\{T_{k,m}^{i}(\zeta)\\}_{m=1}^{M_{k}}$ via (20) using
$\mathbf{w}_{kk}[\nu,i]$;

9: else

10:
$\\{T_{k,m}^{i}(\zeta)\\}_{m=1}^{M_{k}}=\\{T_{k,m}^{i-1}(\zeta)\\}_{m=1}^{M_{k}}$;

11: end if

12: Shift WOLA window ($N_{\mathrm{s}}$ new samples since $i-1$);

13: WOLA analysis on local signals to obtain $\mathbf{y}_{k}[\nu,i]$;

14: WOLA analysis on fused signals to obtain $\mathbf{z}_{-k}[\nu,i]$;

15: for $q\in\mathcal{K}_{k}$ do

16: Detect FSDs based on number of new $\dot{z}_{q}[n]$ samples and compute
$\phi_{kq}^{\mathrm{ac}}[\nu,i]$ $\forall\>\nu$ via (18) and (19);

17: Compute $\hat{\varepsilon}_{kq}[i]$ via (16);

18: Compute $\check{z}_{q}[\nu,i]$ (Figure 4) and build
$\check{\mathbf{z}}_{-k}[\nu,i]$;

19: end for

20: Build
$\check{\mathbf{y}}_{k}[\nu,i]=[\mathbf{y}_{k}^{T}[\nu,i]\>|\>\check{\mathbf{z}}_{-k}^{T}[\nu,i]]^{T}$;

21: Compute $\check{\mathbf{R}}_{\mathbf{y}_{k}\mathbf{y}_{k}}[\nu,i]$ and
$\check{\mathbf{R}}_{\mathbf{n}_{k}\mathbf{n}_{k}}[\nu,i]$;

22: Compute $\check{\mathbf{w}}_{k}[\nu,i+1]$ via (9), then
$\hat{d}_{k}[\nu,i]$.

23: WOLA synthesis on $\hat{d}_{k}[\nu,i]$ and overlap-add with the previous
frame to build $\dot{\hat{d}}_{k}[n]$.

24: end for

##  5 Numerical experiments

The performance of Algorithm 2 is demonstrated and compared to Algorithm 1 via
numerical experiments. The acoustic environment is depicted in Figure 5. A
WASN of $K=4$ nodes is considered, with $\\{M_{k}\\}_{k=1}^{4}=\\{1,3,2,5\\}$
microphones with a 20 cm inter-microphone spacing. A 5$\times$5$\times$5 m3
room with a uniform absorption coefficient of 0.9 is considered, resulting in
a $T_{60}=0.15$ s reverberation time. One localised speech source and two
localised uncorrelated stationary white noise sources are present (note that
the validity of Algorithm 2 can also be demonstrated in the presence of a non-
stationary noise source such as babble noise). The speech signal consists of 3
s long LibriSpeech [24] snippets, each separated by 2 s of silence and
starting with 0.25 s of silence. The power of each source is set to obtain a
-3 dB signal-to-noise ratio (SNR) at the reference microphone of node 1. All
signals last 15 s and are simulated by convolving the source signals with 4096
samples room impulse responses obtained using the randomised image method
[25]. The nominal sampling rate is set to 16 kHz.

The filters are initialised as selecting the local reference microphone
signal, i.e.,
$\tilde{\mathbf{w}}_{k}[\nu,0]=[1\>\mathbf{0}]^{T}\>\forall\>k\in\mathcal{K}$.
The covariances matrices are updated using $\beta=0.978$ (cfr. (4)). All WOLA
processing is performed using $N=1024$-samples square-root Hann windows with
50% overlap (note that the conclusions presented here in terms of speech
enhancement are also valid for other frame lengths, e.g., $N=512$ or $2048$
samples). FSDs are detected using the WOLA approximation described in Section
4.3, where the distortion function $T_{k,m}^{i}(\zeta)$ in (20) is updated
based on the filter $\mathbf{w}_{kk}[\nu,i]$ every 30 DANSE iterations. The
covariance matrices are estimated via (4) assuming an ideal VAD, which avoids
the influence of VAD errors on the results. In practice, the VAD obviously
needs to be estimated from the microphone signals [15, 16].

Figure 5: Layout of acoustic scenario used in the simulations. Microphones
($\circ$) are grouped in nodes ($\bigcirc$) numbered 1 to 4. Two noise sources
(+) and one desired source ($\diamond$) are present.

The clock of node 1 is set as the reference, with $f_{\mathrm{s},1}=16$ kHz.
The SRO for all other nodes $k\in\\{2,3,4\\}$ is defined with respect to this
reference. Three degrees of network asynchronicity are considered based on the
measured SROs values reported in [5]. First, small SROs are considered by
setting $\\{\varepsilon_{1k}\\}_{k=2}^{4}=\\{20,-20,40\\}$ PPM. Second, more
asynchronicity is applied by setting
$\\{\varepsilon_{1k}\\}_{k=2}^{4}=\\{50,-50,100\\}$ PPM. Finally, a strongly
asynchronous network is simulated by setting
$\\{\varepsilon_{1k}\\}_{k=2}^{4}=\\{200,-200,400\\}$ PPM. Fixed SROs are
simulated at any node by resampling the signals appropriately. The SRO
estimation method uses $l_{\mathrm{d}}=10$ in (13) and $\alpha=0.95$ in (15),
resulting in a $\pm$ 3 PPM accuracy.

The performance at each node is quantified using the extended short-term
objective intelligibility (eSTOI) [26] with the clean speech component of the
first local microphone as reference. This metric is particularly relevant as
opposed to, e.g., SNR, as intelligible speech is of central interest in most
speech enhancement applications. The eSTOI is computed on the signal segment
starting from WOLA frame $i=15$ to reduce the impact of initial filter
updates. For each degree of asynchronicity, Figure 6 shows the eSTOI at each
node for the local reference microphone signal (without any noise reduction),
the desired signal estimate from Algorithm 1 without SRO compensation, and the
desired signal estimate from Algorithm 2 with the proposed SRO compensation
with or without compensating for FSDs. The eSTOI obtained using the
synchronised and centralised GEVD-MWF (cfr. (6)) is provided for comparison.

Figure 6: eSTOI obtained at each node from Figure 5 with small SROs (top),
moderate SROs (middle), or large SROs (bottom). Local noisy reference
microphone signals without any processing (lightest), DANSE estimates in the
presence of SROs (light), DANSE estimates with SRO estimation and
compensation, without FSD compensation (dark) and with (darkest), and MWF SRO-
free centralised estimates (dashed).

The results show that the presence of SROs in the WASN significantly
deteriorates the performance of WOLA-based GEVD-DANSE. The single-microphone
node ($M_{1}=1$) is particularly sensitive to the presence of SROs as it
heavily relies on the information provided by other nodes to compute its
desired signal estimate. This occurs regardless of the considered SRO, which
shows the negative impact of even relatively small SROs. Conversely, nodes
including many microphones (e.g., $M_{4}=5$), show almost no sensitivity to
SROs, suggesting that these nodes are able to rely solely on their locally
recorded signals to perform noise reduction with a comparable performance as
in the centralised case. For all considered SRO magnitudes, each node using
the proposed method with FSD compensation is able to restore the centralised
performance that GEVD-DANSE would showcase in an SRO-free WASN.

##  6 Conclusion

In this contribution, the WOLA-based implementation of the GEVD-DANSE
algorithm has been rendered robust to the presence of SROs by combining a
coherence-based SRO estimation technique with an approximation of the WOLA
process to allow FSDs detection and compensation via per-sample broadcasting
of fused signals. The performance of the proposed method has been assessed
through numerical experiments in the context of speech enhancement. in terms
of intelligibility of the desired signal estimate at each node. The results
show that even relatively small SROs (if not estimated and compensated for)
can have a detrimental impact on the ability of DANSE to recover the desired
signal at nodes that significantly rely on collaboration with other nodes.
However, it is shown that, in an asynchronous WASN, the proposed SRO
estimation and compensation method practically restores the performance that
the GEVD-DANSE algorithm would showcase in a fully synchronised network.

## References

  * [1] A. Bertrand,  “Applications and trends in wireless acoustic sensor networks: A signal processing perspective,”  in Proc. IEEE Symp. Commun. Veh. Technol., 2011, pp. 1–6.
  * [2] A. Bertrand and M. Moonen,  “Distributed adaptive node-specific signal estimation in fully connected sensor networks—Part I: Sequential node updating,”  IEEE Trans. Signal Process., vol. 58, no. 10, pp. 5277–5291, 2010\.
  * [3] A. Bertrand and M. Moonen,  “Distributed adaptive node-specific signal estimation in fully connected sensor networks—Part II: Simultaneous and asynchronous node updating,”  IEEE Trans. Signal Process., vol. 58, no. 10, pp. 5292–5306, 2010\.
  * [4] S. Ruiz, T. van Waterschoot, and M. Moonen,  “Distributed combined acoustic echo cancellation and noise reduction in wireless acoustic sensor and actuator networks,”  IEEE/ACM Trans. Audio, Speech, Language Process., vol. 30, pp. 534–547, 2022.
  * [5] M. Guggenberger, M. Lux, and L. Böszörmenyi,  “An analysis of time drift in hand-held recording devices,”  in Proc. Int. Conf. MultiMedia Model, 2015, pp. 203–213.
  * [6] R. Lienhart, I. Kozintsev, S. Wehr, and M. Yeung,  “On the importance of exact synchronization for distributed audio signal processing,”  in Proc. IEEE Int. Conf. Acoust., Speech, Signal Process., 2003, pp. 840–843.
  * [7] A. Hassani, A. Bertrand, and M. Moonen,  “GEVD-based low-rank approximation for distributed adaptive node-specific signal estimation in wireless sensor networks,”  IEEE Trans. Signal Process., vol. 64, no. 10, pp. 2557–2572, 2016\.
  * [8] J. Szurley, A. Bertrand, and M. Moonen,  “Topology-independent distributed adaptive node-specific signal estimation in wireless sensor networks,”  IEEE Trans. Signal Inf. Process. Netw., vol. 3, no. 1, pp. 130–144, 2017.
  * [9] R. Van Rompaey and M. Moonen,  “Distributed adaptive signal estimation in wireless sensor networks with partial prior knowledge of the desired sources steering matrix,”  IEEE Trans. Signal Inf. Process. Netw., vol. 7, pp. 478–492, 2021\.
  * [10] J. Zhang and P. Wu,  “Joint sampling synchronization and source localization for wireless acoustic sensor networks,”  IEEE Commun. Lett., vol. 24, no. 5, pp. 1020–1023, 2020.
  * [11] D. Hu, H. Zhang, F. Bao, and R. Wang,  “Distributed sampling rate offset estimation over acoustic sensor networks based on asynchronous network newton optimization,”  IEEE/ACM Trans. Audio, Speech, Language Process., vol. 31, pp. 301–312, 2023.
  * [12] J. Schmalenstroeer, J. Heymann, L. Drude, C. Boeddecker, and R. Haeb-Umbach,  “Multi-stage coherence drift based sampling rate synchronization for acoustic beamforming,”  in Proc. Int. Workshop Multimedia Signal Process., 2017, pp. 1–6.
  * [13] T. Gburrek, J. Schmalenstroeer, and R. Haeb-Umbach,  “On synchronization of wireless acoustic sensor networks in the presence of time-varying sampling rate offsets and speaker changes,”  in Proc. IEEE Int. Conf. Acoust., Speech, Signal Process., 2022, pp. 916–920.
  * [14] A. Bertrand and M. Moonen,  “Robust distributed noise reduction in hearing aids with external acoustic sensor nodes,”  EURASIP J. Adv. Signal. Process., vol. 2009, no. 1, pp. 530435, 2009\.
  * [15] A. Bertrand and M. Moonen,  “Energy-based multi-speaker voice activity detection with an ad hoc microphone array,”  in Proc. IEEE Int. Conf. Acoust., Speech, Signal Process., 2010, pp. 85–88.
  * [16] Y. Zhao, J. K. Nielsen, J. Chen, and M. G. Christensen,  “Model-based distributed node clustering and multi-speaker speech presence probability estimation in wireless acoustic sensor networks,”  J. Acoust. Soc. Am., vol. 147, no. 6, pp. 4189–4201, 2020.
  * [17] R. Serizel, M. Moonen, B. Van Dijk, and J. Wouters,  “Low-rank approximation based multichannel Wiener filter algorithms for noise reduction with application in cochlear implants,”  IEEE/ACM Trans. Audio, Speech, Language Process., vol. 22, no. 4, pp. 785–799, 2014.
  * [18] J. Szurley, A. Bertrand, and M. Moonen,  “Improved tracking performance for distributed node-specific signal enhancement in wireless acoustic sensor networks,”  in Proc. IEEE Int. Conf. Acoust., Speech, Signal Process., 2013, pp. 336–340.
  * [19] R. Crochiere,  “A weighted overlap-add method of short-time Fourier analysis/Synthesis,”  IEEE/ACM Trans. Audio, Speech, Language Process., vol. 28, no. 1, pp. 99–102, 1980.
  * [20] S. Miyabe, N. Ono, and S. Makino,  “Blind compensation of inter-channel sampling frequency mismatch with maximum likelihood estimation in STFT domain,”  in Proc. IEEE Int. Conf. Acoust., Speech, Signal Process., 2013, pp. 674–678.
  * [21] L. Wang and S. Doclo,  “Correlation maximization-based sampling rate offset estimation for distributed microphone arrays,”  IEEE/ACM Trans. Audio, Speech, Language Process., vol. 24, no. 3, pp. 571–582, 2016.
  * [22] M. Bahari, A. Bertrand, and M. Moonen,  “Blind sampling rate offset estimation for wireless acoustic sensor networks through weighted least-squares coherence drift estimation,”  IEEE/ACM Trans. Audio, Speech, Language Process., vol. 25, no. 3, pp. 674–686, 2017.
  * [23] P. P. Vaidyanathan,  Multirate systems and filter banks,  Prentice Hall, Englewood Cliffs, NJ, 1993.
  * [24] V. Panayotov, G. Chen, D. Povey, and S. Khudanpur,  “Librispeech: an ASR corpus based on public domain audio books,”  in Proc. IEEE Int. Conf. Acoust., Speech, Signal Process., 2015, pp. 5206–5210.
  * [25] E. De Sena, N. Antonello, M. Moonen, and T. van Waterschoot,  “On the modeling of rectangular geometries in room acoustic simulations,”  IEEE/ACM Trans. Audio, Speech, Language Process., vol. 23, no. 4, pp. 774–786, 2015.
  * [26] J. Jensen and C. Taal,  “An algorithm for predicting the intelligibility of speech masked by modulated noise maskers,”  IEEE/ACM Trans. Audio, Speech, Language Process., vol. 24, no. 11, pp. 2009–2022, 2016.


# Entropy from entangled parton states
and high-energy scattering behavior

Hans Günter Dosch<EMAIL_ADDRESS>Institut für Theoretische Physik der
Universität, D-69120 Heidelberg, Germany  Guy F. de Téramond
<EMAIL_ADDRESS>Laboratorio de Física Teórica y Computacional,
Universidad de Costa Rica, 11501 San José, Costa Rica  Stanley J. Brodsky
<EMAIL_ADDRESS>SLAC National Accelerator Laboratory, Stanford
University, Stanford, CA 94309, USA

###### Abstract

The relation between the gluon density in a hadron and entanglement entropy
can shed a new light on the high energy scattering behavior of hadrons. Using
the holographic light-front QCD framework the growth above the classical
geometric cross section is directly related to the increase of the internal
quantum entropy from the entangled parton distribution in hadrons. A rather
consistent picture emerges from the scale dependence of the Pomeron from the
QCD evolution of the gluon distribution function $g(x,\mu)$, the rising of the
integrated cross section in photoproduction of vector mesons, the deep
inelastic scattering (DIS) experiments at HERA, hadron multiplicity and
quantum entropy. We also point out a possible analogy between parton
entanglement entropy and the black-hole entropy of Bekenstein and Hawking.

††preprint: SLAC-PUB-17727

##  I Introduction

In their analysis of entanglement at the subnuclear level, Kharzeev and Levin
[1, 2] have drawn the attention to entropy. The quantum entropy of any bound
state is zero, but Kharzeev and Levin considered the entropy of the partons
resolved in a deep inelastic scattering (DIS) experiment

$\displaystyle S_{DIS}=\ln N(x,Q^{2}),$ (1)

where $N(x,Q^{2})$ is the number of partons in a hadron with longitudinal
light-front momentum fraction $x$ of the struck parton in the target hadron
and $Q^{2}=-q^{2}>0$ is the momentum transfer. The DIS entropy $S_{DIS}$ is
the logarithm of the number of degrees of freedom in the DIS measurement. It
represents the entropy of entanglement between the proton components probed by
deep inelastic scattering and the rest of the proton: the number of produced
final-state spectator quark and gluon partons. Since the partons cannot be
isolated as asymptotic states, the number of partons is not a directly
observable quantity, but depends on the virtuality scale of the process
$Q^{2}$, and therefore the quantum DIS entropy is not directly observable.

The quantum (von-Neumann) entropy of a state described by the density operator
$\rho$ is given in analogy to the classical (Boltzmann) entropy by the
expectation value of the trace of the statistical operator

$\displaystyle S_{Q}=-tr[\rho\,\ln\rho]=-\sum_{i}p_{i}\ln p_{i},$ (2)

where the $p_{i}$ are the eigenvalues of $\rho$; they give the probability to
find the system in the state $|i\rangle$. A pure state $|\psi\rangle$, like an
elementary particle, has therefore the quantum entropy 0. Under the assumption
that a hadron state $|\psi\rangle$ consists of $N$ interacting constituents
and that these partons do not factorize into hadronic substates, the partons
in the hadron are entangled.

In deep inelastic scattering (DIS) a measurement projects only on a single
parton $|j\rangle$. The pure state before the measurement is thus, after a
measurement, a mixture described by the statistical operator
$\rho=\sum_{j=1}^{N}\,p_{j}\,|j\rangle\langle j|$, where $p_{j}$ is the
probability of hitting the state $|j\rangle$; that is, after measurement,
which traces out the unobserved components of the state, the entropy is given
by the expression (2). For very slow partons the gluons are dominant, and we
restrain ourselves to this region in order to treat all partons on equal
footing. Therefore, we have equal probability for each parton to be hit and
thus one expects $p_{j}=1/N(x)$, hence $S_{Q}=\ln N(x)$, where $N(x)$ is the
number of partons (gluons) with longitudinal momentum fraction $x$, and the
result is (1).

In this letter we discuss the consequences of the relation of the DIS entropy
introduced in Ref. [1] with the determination of the gluon distribution given
in Ref. [3] from the high energy scattering of hadrons and holographic light-
front QCD [4, 5, 6]. This specific connection allows us to relate the increase
above the proton-proton geometric cross section to the number of entangled
partons in the proton probed in the DIS process. Furthermore, using the light-
front holographic framework and the QCD evolution results from Ref. [7], we
show that the scale dependence of the gluon distribution produces the Pomeron-
dominated energy dependence of the DIS cross section. It also accounts for the
logarithmic dependence on the observational scale $\mu$, which, for our
purposes, can be identified with the photon virtuality $\mu^{2}=Q^{2}$ in the
DIS measurement $\gamma^{*}(q^{2})+p\to X$; this is an essential prediction of
the scale-dependent Pomeron intercept discussed in Ref. [7].

In Refs. [1, 2, 8, 9] it was further assumed that the entropy of the final
hadronic state can, under certain kinematical conditions, be equated to the
DIS entropy. This assumption created delicate experimental and theoretical
questions, since parton states are not asymptotic states and the final state,
after a DIS event, consists of hadrons which are supposed to be formed after
the hadronization of the partons. By examining the hadron multiplicity in the
final state, we will briefly discuss how this problem can be also addressed
using the concept of a unique scale-dependent Pomeron introduced in Ref. [10],
and recently discussed in Ref. [7] in the framework of holographic light-front
QCD (HLFQCD).

##  II High energy scattering and entanglement entropy

###  II.1 High energy behavior of hadron cross sections

One of the most important concepts and tools in high energy scattering is the
notion of Regge poles. Originally introduced in scattering in quantum
mechanical potential theory by T. Regge [11], it was applied to particle
physics by Chew and Frautschi [12]. The guiding principle of the latter
authors was the notion of duality. This principle states, roughly speaking,
that the same fields which create the interaction between hadrons are the
quantum fields which show up as particles or resonances. It follows from
rather general principles of quantum field theory and is therefore independent
of the specific dynamics.

Figure 1: The total proton-proton scattering cross section as a function of
$s$, the CM energy squared. The straight grey line in the figure shows the
energy dependence $s^{0.086}$, as predicted by a Pomeron with intercept
$\alpha_{P}(0)=1.086$. The pale blue data points are results from cosmic ray
experiments and the red are recent LHC data. Figure adapted from [13], Fig.
52.6.

In lowest order Regge theory of particle scattering (one pole exchange) the
dependence of the total cross section on the squared center of mass energy $s$
of the two scattered particles is obtained from the forward scattering
amplitude using the optical theorem. Therefore it is determined at $t=0$ by
the intercept of the trajectory in a Chew-Frautschi plot [12] and the cross
section is proportional to $s^{\alpha_{P}(0)-1}$. In Fig. 1 a double
logarithmic plot of the proton-proton total cross section vs the squared
center of mass energy $s$ is displayed. One sees that the Pomeron contribution
$s^{0.086}$ gives an excellent description of the cross section over 6 orders
of magnitude of $s$. For values of $s\lesssim 10^{3}$ GeV2 another Regge pole
exchange, the $\rho$ trajectory with an intercept 1/2 [14] gives an additional
contribution which vanishes as $s^{-1/2}$ at large $s$. By adding both
components

$\displaystyle\sigma_{pp}=a(s/s_{0})^{0.086}+b(s/s_{0})^{-0.5},$ (3)

with $a=24\,\rm{mb}$, $b=28\,\rm{mb}$ and $s_{0}\equiv 4m_{p}^{2}$, we find a
very good description of the total proton-proton cross section over 9 orders
of magnitude in terms of the Pomeron and the $\rho-\omega$ trajectories only.

The Froissart-Martin bound [15, 16] states, however, that the total cross
section of two colliding hadrons cannot grow faster than
$\ln^{2}\left(s/s_{0}\right)$. This result is based on very general
assumptions of quantum field theory, especially on the unitarity of the
scattering matrix. Therefore, the very high energy behavior described by
single Pomeron exchange has to be modified by unitarity corrections (two and
more Pomeron exchanges), but as can be seen from Fig. 1 these corrections
would affect only a region well above $10^{10}$ GeV2.

For later discussion we add here a general remark on the high energy behavior:
From classical, even wave mechanical, considerations one expects a constant
cross section determined by the geometrical extension of the scattered
objects. However, in the quantum field theoretical shock-wave model put
forward by Heisenberg [17], the total cross section increases as $\ln^{2}s$
due to the possibility of particle creation. It thus is equal to the upper
limit later predicted by the Froissart-Martin bound based on principles of
quantum theory.

###  II.2 DIS entropy and the Pomeron intercept

In QCD the gluon distribution at a given virtuality scale $\mu$ is obtained
from DGLAP evolution [18, 19, 20] from the behavior of the DIS cross sections
$\gamma^{*}+p\to X$. The gluon distribution $g(x,\mu)$ is not an observable,
and consequently, the parton distribution function depends also on the scale
$\mu$

$\displaystyle S_{DIS}=\ln\big{(}x\,g(x,\mu)\big{)}.$ (4)

This is in agreement with the result mentioned in the introduction, namely,
that also the DIS entropy is not a directly observable quantity.

Based on AdS/CFT duality concepts, the correspondence between a gravity theory
in a five-dimensional anti-de Sitter (AdS) space and conformal field theories
(CFT) in physical space-time [21, 22, 23], a semiclassical model has been
developed [4, 5, 6] which reproduces not only the hadron spectra from hadronic
two point functions but also allows one to obtain analytical expressions for
classical three point expectation values in 5-dimensional AdS metric. From
these expressions follow three point functions of hadronic states with vector
or tensor currents in the full momentum transfer range. From those one can
obtain form factors, Regge trajectories and particularly the intercept of the
trajectory, $\alpha(0)$. For the vector-isovector current the intercept is
fixed by the spectra to $\alpha_{\rho}(0)=0.5$ in fair agreement with
phenomenology [14]. For the gravitational, i. e. a spin-two current, Pomeron
exchange is identified as the graviton of the dual AdS theory [3, 24].

In the light front holographic approach the form factors can be expressed by
universal functions of the generalized parton distributions at zero skewness
[14] which yields, for the gluon density at small $x$, the relation [3]

$\displaystyle
x\,g(x,\mu)\sim\left(\frac{1}{x}\right)^{\alpha_{P}(0)-1},\qquad x\ll 1,$ (5)

in the hadronic domain $\mu\approx 1$ GeV. This calculated density sets the
initial scale for the DGLAP evolution equations. From (4) and (5) we then
obtain

$\displaystyle
S_{DIS}=\ln\big{(}x\,g(x,\mu)\big{)}\sim(\alpha_{P}(0)-1)\ln\Big{(}\frac{1}{x}\Big{)},\qquad
x\ll 1.$ (6)

It shows that a Pomeron intercept larger than 1 (hypercritical Pomeron), which
leads to the rising of the total proton-proton cross section at large energies
(Fig. 1), is uniquely related to a positive von Neumann entropy since the
rapidity $Y=\ln(1/x)$ is always positive for $x<1$. It is, in turn, a
consequence of the separation of entangled parton states in a DIS experiment
which gives rise to the entanglement entropy $S_{DIS}$.

The longitudinal light-front momentum fraction $x$ can be identified with the
Bjorken variable $x_{bj}$ measured in the deep inelastic lepton-hadron
scattering experiment

$\displaystyle x_{bj}=\frac{Q^{2}}{2p\cdot
q}=\frac{Q^{2}}{W^{2}+Q^{2}-m_{p}^{2}},$ (7)

where $W^{2}=(q+p)^{2}$ represents the total photon-hadron energy squared and
$Q^{2}=-q^{2}>~{}0$ is the photon virtuality. In the high-energy kinematic
domain, $W^{2}\gg Q^{2}$, Eq. 7 reduces to $x=\frac{Q^{2}}{W^{2}}$ and, for
fixed $Q^{2}$, we obtain

$\displaystyle
S_{DIS}\sim(\alpha_{P}(0)-1)\ln\left(\frac{W^{2}}{Q^{2}}\right),\qquad
W^{2}\gg Q^{2}.$ (8)

It leads to the remarkable result that one can identity the value of
$\alpha_{P}(0)$, which determines the gluon distribution $xg(x,\mu)$ at small
$x$ in Eq. (5), with the value of $\alpha_{P}(0)$, obtained from the high
energy Pomeron contribution to the $W$ dependence of the DIS cross section.

The quark distributions evaluated in the light front holographic approach are
in agreement with global fits [14, 25] and have recently been extended to
predict the strange and charm distributions in the proton [26, 27]. Likewise,
the gluon distribution depends crucially on the relation (5) and the value of
the soft Pomeron intercept $\alpha_{P}(0)=1.086$. It has been inserted as the
starting distribution at the hadron scale for the DGLAP evolution and also
compares favorably with global fits without any additional ingredient [3].
This opens the way for the attractive possibility that there is only a single
Pomeron with a scale dependent intercept: It manifests itself as a soft
Pomeron at soft scales and hard Pomeron at harder ones, as was already
considered in Refs. [10, 7]. In fact, as we discuss in the following
subsection, the relation between the Pomeron intercept and the DIS
entanglement entropy makes this possibility even more compelling. Also in the
framework of the holographic approach we further discuss in Sec. II D the
relation of hadron multiplicity with entanglement entropy first discussed in
[1, 2, 8, 9]. In the final section III we point out also other relations
between a hypercritical Pomeron intercept and rising entanglement entropy in
high energy scattering in an asymptotically free theory.

###  II.3 The scale dependent Pomeron

Results from quasi-elastic photo and electroproduction of vector mesons show,
for all processes at center of mass energies $W$ larger than 10 GeV above
threshold, an energy dependence of the integrated cross sections that can be
fitted reasonably well by a power law $\sigma_{\gamma^{*}+p\to\,V+p}\approx
W^{\delta}$. The value of $\delta$ depends strongly on the process and the
virtuality of the photon as shown in Table 1. For integrated exclusive cross
sections the whole trajectory contributes. In order to relate the exclusive
cross sections, integrated over the whole range of momentum transfer $t$ to
the Pomeron intercept at $t=0$ one has to take into account the shrinkage
corrections due to the slope of the trajectory. The latter can be estimated,
and from this, one can reconstruct the energy dependence of the squared
forward amplitude which, in Regge theory, is given by
$\delta_{Regge}=4(\alpha_{P}(0)-1)$. The values for this quantity are shown in
Table 1.

Reaction | $\delta$ | $\alpha_{P}-1$ | $\chi^{2}$
---|---|---|---
$\gamma\,p\to\rho\,p$ | 0.190 | 0.090 | 0.64
$\gamma^{*}\,p\to\rho\,p\,(Q^{2}=6\,{\rm GeV}^{2})$ | 0.486 | 0.135 | 0.66
$\gamma\,p\to J/\psi\,p$ | 0.675 | 0.173 | 0.81
$\gamma\,p\to\Upsilon\,p$ | 0.8626 | 0.260 | 0.46
Table 1:  Pomeron intercept from integrated cross sections for photo and
electroproduction of vector mesons. The cross sections were fitted by the
curve $\sigma=(W/W_{0})^{\delta}$. The intercept $\alpha_{P}$ was obtained by
taking into account the shrinkage corrections from the Regge trajectory slope,
see [10].

The situation is similar for the inclusive cross section for the DIS process
$\gamma^{*}+p\to X$. The proton structure function $F_{2}(x,Q^{2})$ has been
studied over a large range of values for $x$ (see (7)) and $Q^{2}$ and the
results for the $x$-dependence of the structure function could be fitted by
the power behavior [28]

$\displaystyle F_{2}\sim xg\left(x,Q^{2}\right)\sim
x^{-\lambda(Q^{2})},~{}~{}\mbox{with}~{}~{}\lambda(Q^{2})=0.048\ln\left(\frac{Q^{2}}{\Lambda^{2}}\right),$
(9)

and $\Lambda\simeq 290$ MeV. The structure function is related to the total
DIS cross section for the process $\gamma^{*}+p\to X$ by

$\displaystyle\sigma_{\gamma^{*}+p\to
X}(W^{2})=\frac{4\pi^{2}\alpha_{em}}{Q^{2}}F_{2}\left(x,Q^{2}\right).$ (10)

If one applies Regge theory to the elastic forward scattering amplitude
$\gamma^{*}p\to\gamma^{*}p$ and relates it by the optical theorem to the DIS
cross section $\gamma^{*}+p\to X$, one expects a dependence
$\sigma_{\gamma^{*}+p\to X}\sim(W^{2})^{\lambda}$, with
$\lambda=\alpha_{P}(0)-1$, whereas the observed value of $\lambda$ varies
between $\lambda=0.16\pm 0.015$ at $Q^{2}=2$ GeV2 and $\lambda=0.36\pm 0.11$
at $Q^{2}=150$ GeV2 [29, 30, 28].

This breakdown of Regge theory in electromagnetic processes led not to a
crisis of this theory, since Regge theory was created for purely hadronic
processes. But on the other hand, in quantum field theory a photon couples
inevitably to intermediate hadronic matter and the lifetime of these
intermediate states is large as compared to the time scale of a hard
scattering process [31]; In fact, in the case of vector meson production the
lifetime of the hadronic state is even infinitely long. Therefore DIS
scattering of (virtual) photons and diffractive electroproduction should be
treated like hadron-hadron scattering and Regge theory should be applicable.

It was indeed shown by Donnachie and Landshoff that Regge theory can be
applied to these processes by introducing two Pomerons [32]. One Pomeron was
the usual “soft” Pomeron, well established in hadron scattering, the other one
was tentatively related to the so called BFKL [33, 34] “hard” Pomeron [31].
This emerges in perturbative QCD by the exchange of gluon ladders and leads to
the effective power $\lambda_{h}$

$\displaystyle\lambda_{h}\equiv\alpha_{h}-1=\frac{4N_{c}\alpha_{s}}{\pi}\ln
2,$ (11)

where $\alpha_{h}$ is the hard Pomeron intercept and $\alpha_{s}$ is the QCD
coupling constant. This value is rather unstable against higher order
contributions, but with resummation corrections an intercept between 1.3 and
1.5 is plausible. With the conventional soft Pomeron and a BFKL Pomeron, with
an intercept $\alpha^{h}\approx 1.42$, all electron and photoproduction data
in the available energy range could be fitted with reasonable accuracy [35].
It turned out that the soft Pomeron couples dominantly to extended objects and
the hard one to smaller ones [36].

In DIS production the effective size of the intermediate hadronic system
coupled to the virtual photon decreases with increasing virtuality. If one
identifies the power $\lambda(Q^{2})$ obtained from DIS experiments (9) with
an effective scale-dependent Pomeron intercept

$\displaystyle\lambda(Q^{2})=\alpha_{P}(0,Q^{2})-1,$ (12)

this effective intercept increases indeed by a factor 2 from $Q^{2}=2$ GeV2 to
$Q^{2}=150$ GeV2. For the photoproduction of heavy vector mesons the
transverse size of the vector meson is relevant, and is roughly proportional
to the inverse mass and for a $J/\Psi$ and $\Upsilon$ meson, smaller than a
normal hadronic size. As mentioned above, this effective intercept could for
all energies available before LHC be explained by the two Pomeron model [32,
35, 37].

Figure 2:  In photoproduction of hadrons the scale is set by the virtuality of
the mass of the produced $J/\psi$ meson, $Q^{2}\simeq M^{2}_{J/\Psi}$. The
production cross section for $J/\psi$ mesons can be described by a single
Pomeron exchange with intercept $\alpha_{P}(0)\simeq 1.15$ (solid line),
distinctly higher than the intercept of conventional hadrons. The figure is
adapted from [10], LHC data from [38, 39].

The two Pomeron picture, however, seems to be ruled out by the appearance of
data for quasi-elastic photoproduction of $J/\Psi$ [38, 39] and $\Upsilon$
[40] vector mesons at LHC energies. As can be seen from Fig. 2, the power law
for $J/\psi$ production holds over the range from 20 GeV to above 1 TeV. This
is impossibly to achieve with two Pomerons, since at higher energies the hard
one would dominate and the curve should show clear convexity. This leads to
the introduction of a single, but scale dependent Pomeron [10] on a purely
phenomenological basis. We thus write the $Q$-dependent effective power of the
proton structure function $F_{2}(x,Q^{2})\propto x^{\lambda(Q^{2})}$ with the
scale dependent Pomeron intercept $\alpha_{P}(0,Q^{2})$ given by (12), and the
experimentally determined function $\lambda(Q^{2})$ (9) from [28], as shown in
Figure 3.

In Ref. [7] we have examined the possibility that Eq. (5) is not only valid at
the hadronic scale, but at all scales, thus leading to a single scale
dependent Pomeron, $\alpha_{P}(0,\mu)$. Since the small-$x$ behavior is
determined uniquely by the scale-dependent intercept it can, in principle, be
extracted unambiguously.

Figure 3:  The effective power $\lambda(\mu)=\alpha_{P}(0,\mu)-1$ extracted
from the gluon distribution $g(x,\mu)$ for the pion (empty red circles) using
the procedure in [7], is compared with the values of $\lambda(Q^{2})$ obtained
from the measured proton structure function $F_{2}(x,Q^{2})$ under the
assumption $\mu^{2}=Q^{2}$. Experimental results from HERA [28] (filled black
circles).

For very small values of $x$ the structure function $F_{2}(x,Q^{2})$ is
strongly dominated by the gluon distribution and therefore the $x$ dependence
obtained in the HERA analysis [28], and parametrized in (9), yields a model
independent determination of the Pomeron intercept if one accepts the relation
between DIS-entropy, gluon density and Pomeron intercept expressed in (20).
This allows the following comparison: The gluon distribution function obtained
in [3] can be evolved using the DGLAP equations and the scale dependence of
the intercept can be obtained as described above. This is illustrated in Fig.
3, where we compare the values for our evolved results for the effective power
$\lambda(\mu)$, Eq. (12), in the range $0.0001\leq x\leq 0.00016$ following
the procedure described in [7], with the power $\lambda(Q^{2})$ of the proton
structure function $F_{2}$ from Ref. [28]. We show only the evolution results
for the pion, since the proton results present numerical instabilities at
lower evolution scales, currently under investigation. It should be noted,
however, that this is rather a consistency check than an independent
determination, since one has to relate the virtuality scale $\mu$ in the gluon
distribution function $g(x,\mu)$ to the photon virtuality $Q$ in the DIS
process [28, 7] and the measured structure function is an important input for
the evolution equations determining the $\mu$ dependence of the gluon
distribution $g(x,\mu)$.

###  II.4  Hadron entropy and multiplicity

It was suggested in Ref. [1] a relation between the hadron multiplicity in the
final state of the collision and the multiplicity of partons

$\displaystyle S_{partons}=\ln\big{(}xg(x,Q^{2})\big{)}\equiv S_{hadrons}.$
(13)

Such a conjectural relation would be an indication of the absence of a
significant increase of entropy in the hadronization process, as measured in
the hadron multiplicity distributions [1, 2, 41]. In the present approach,
this connection leads to the unexpected relevance of hadron multiplicity to
the scale dependence of the Pomeron intercept, and fixes an upper bound to the
growth’s rate of hadron multiplicity at higher energies.

The theoretical analysis of hadron multiplicity distributions
$P_{n}=\sigma_{n}/\sigma_{in}$, the ratio of the cross section to produce $n$
hadrons in a high-energy collision to the inelastic cross section, is a rather
complex undertaking [42, 43]. Most of the hadrons in the collision process are
produced at relatively low transverse momentum and one has to recur to
nonperturbative models or Monte Carlo event generator models of particle
production.

A simple model which incorporates some of the main features of the particle
multiplicity distributions is the parton cascade model of Levin and Lublinski
[44] for $P_{n}$:

$\displaystyle\frac{dP_{n}(Y)}{dY}=-\lambda_{h}nP_{n}(Y)+(n-1)\lambda_{h}P_{n-1}(Y),$
(14)

with the solution

$\displaystyle
P_{n}(Y)=e^{-\lambda_{h}Y}\left(1-e^{-\lambda_{h}Y}\right)^{n-1},$ (15)

where $Y=\ln(1/x)$ is the rapidity and $\lambda_{h}$ is the hard Pomeron
intercept (11). The mean multiplicity which follows from (14) and (15) is

$\displaystyle\langle
n\rangle=\sum_{n}nP_{n}(x)=e^{\lambda_{h}Y}=\left(\frac{1}{x}\right)^{\lambda_{h}},$
(16)

which reproduces the power increase of the cross section with energy [1,
44].111It was suggested in Ref. [45] that the quantum entropy carried by the
Pomeron to be at the origin of the entropy of the final state.

In the present approach the gluon distribution at small $x$ is given by Eq.
(5) and thus the multiplicity at the observational scale $\mu$ follows from Eq
(13) as discussed in [1]

$xg(x,\mu)\sim\langle
n\rangle=\sum_{n}nP_{n}=\left(\frac{1}{x}\right)^{\alpha_{P}(0,\mu)-1},$ (17)

but with the Pomeron intercept $\lambda(\mu)=\alpha_{P}(0,\mu)-1$ instead of
$\lambda_{h}$. At small values of $x$ we have [46, 47, 48]
$x=\frac{t_{0}}{s}$, where we take the virtuality scale $\mu$ as the average
transverse momentum $t_{0}$ of the produced hadrons or a particle mass in the
final state of the collision [46]. Eq. (17) suggests that

$\displaystyle\langle
n\rangle\sim\left(\frac{s}{t_{0}}\right)^{\alpha_{P}(0,t_{0})~{}-~{}1},$ (18)

where $t_{0}$ is independent of the collision energy $s$, but may differ for
different datasets [46]. We notice that the theoretical predictions for the
evolution of the effective power $\lambda(\mu)$, shown in Fig. 3, suggest that
this quantity saturates at higher scales, setting an asymptotic upper bound
for the rate of change of the hadron multiplicity (18) for $\lambda\simeq
0.4-0.5$, higher than the value $\lambda=1/3$ assumed in [1] from two-
dimensional conformal field theory.

Figure 4:  Left figure: CMS data with $p_{T}>100$ MeV and pseudorapidity
$|\eta|<0.5$ from [49]. Right figure: ALICE data with $p_{T}>150$ MeV and
$|\eta|<0.8$ from [50]. The grey lines in the figure correspond to Eq. (18)
with $\sqrt{t_{0}}\simeq 1.3-1.4$ GeV.

In a recent analysis of LHC data [46] it was found that the average
multiplicity is well described by a power law, but with different powers for
each dataset. As an illustration, we show in Fig. 4 the CMS results from Ref.
[49] which correspond to the Set I examined in Ref. [46] for $p_{T}>100$ MeV
and pseudorapidity $|\eta|<0.5$. The particles from this data set correspond
to gluon production. Set II from ALICE and ATLAS with larger rapidities
$|\eta|<2.4$ also includes the current fragmentation region, and Set III from
ATLAS and CMS, also with large rapidities, contains also hadrons produced in
the perturbative domain. We thus examine here only the results for Set I from
Ref. [46] which abides to the mechanisms and approximations discussed in this
article. We also show in Fig. 4 more recent results from ALICE [50] with
$p_{T}>150$ MeV and $|\eta|<0.8$ not included in [46]. Using the predictions
from Eq. (18) (grey line in Fig. 4), and comparing with the proton evolved
results shown in Fig. 3, we obtain for the average transverse momentum
$\sqrt{t_{0}}\simeq 1.3-1.4$ GeV.

Results for the charm structure function [51, 52] $F^{c\bar{c}}(x,Q^{2})$
indicate that in the presence of heavy quarks, the effective power
$\lambda(Q^{2})$ for the $x$ dependence (9) is larger than for the full
structure function $F_{2}(x,Q^{2})$ (see [36]). This leads to the suggestion
from Eq. (18), that also the multiplicity of hidden charm objects increases
faster than the multiplicity of hadrons with light quarks since, in this case,
the virtuality scale is determined by the charm mass.

##  III Summary and Conclusion

The light-front holographic framework, with essential constraints imposed by
the generalized Veneziano model [53, 54], successfully describes the gluon
distribution $g(x,\mu)$ in nucleons and mesons at the hadronic scale $\mu=1$
GeV. An essential input is its relation of the small-$x$ behavior of the gluon
distribution to the Pomeron intercept $\alpha_{P}(0)$ from high-energy
scattering [3].

The gluon distribution is a scale dependent quantity and the intrinsic
distribution evaluated at the hadronic scale can be extended to other scales
using QCD evolution equations. In Ref. [7] we have ventured that also the
Pomeron intercept follows this evolution and is therefore scale dependent.
Such a scale dependent Pomeron was proposed on purely phenomenological reasons
in Ref. [10], where it was shown that such a dependence is suggested by the
results for the photoproduction of $J/\psi$ mesons and not in contradiction
with the principles of Regge theory. A scale dependent Pomeron trajectory is
also a salient feature of the gauge-gravity correspondence [24]. At very small
values of $x$ the gluon dominance leads to universal behavior, well captured
by the present description, which allows us to relate different experimental
results in the high energy domain.

The relation (6) makes it clear that two features of high energy scattering,
the hypercritical Pomeron intercept $\alpha_{P}(0)>1$ and the single but scale
dependent Pomeron, here proposed, are closely related to typical features of a
non-Abelian theory like QCD. In fact, in analogy of the thermodynamical
definition of classical entropy for a system in equilibrium

$\displaystyle dS=\frac{dQ_{therm}}{T},$ (19)

the relation (6) can be written in the differential form

$\displaystyle dS_{DIS}=d\ln\left(x\,g(x,\mu)\right))\sim
2\,\left(\alpha_{P}(0)-1\right)\frac{dW}{W}.$ (20)

It shows that the hypercritical intercept $\alpha_{P}(0)$ is related to the
increase of the DIS-entropy and the gluon density $x\,g(x,\mu)$ with
increasing energy $W$, corresponding to decreasing $x$. The other salient
feature discussed here, the scale dependence of the intercept, is a
consequence of the scale dependence of the gluon distribution, which in turn
is the consequence of the renormalization group of an UV stable renormalizable
theory like QCD. It is therefore quite natural that in a nonperturbative
holographic theory corresponding to a non-Abelian gauge theory [6] these
features emerge in a transparent way.

If the quantum entropy $S_{DIS}$ of partons inside a hadron equals the
classical hadron entropy $S_{hadron}$ of the final state, as proposed in Ref.
[1], then the final state entropy is also scale dependent and, with it, also
the hadron multiplicity which depends on the average momentum of the hadrons
in the final state; here the transverse momentum could set the scale as
discussed in Sec. II.4. It describes the rate of change of the multiplicity in
terms of the scale-dependent Pomeron, and it sets an upper bound to its
increase from the saturation of the effective power $\lambda(\mu)$ at very
large virtualities predicted from the QCD evolution of the gluon distributions
in holographic QCD.

We conclude this letter with some speculative remarks on a possible analogy
between the black-hole (or Bekenstein-Hawking) entropy222For a short non-
technical review see Jacob D. Bekenstein, Bekenstein-Hawking entropy,
Scholarpedia, 3(10):7375 (2008). and $S_{DIS}$, the quantum entropy from deep
inelastic scattering.

The existence of black-hole entropy is demanded by the laws of thermodynamics
[55], but the observation of microstates inside the black hole seems to be
prevented by the horizon (no hair theorems). There is, however, a way out:
Using the temperature of the black hole, which can be determined from the
Hawking radiation [56] emitted from it, the Bekenstein-Hawkings entropy can be
related to the surface area $A$ of the black hole

$\displaystyle S_{BH}=\frac{A}{4L_{P}^{2}},$ (21)

where $L_{P}$ is the Planck length $L_{P}=\sqrt{G_{N}\hbar/c^{3}}$.

On the other hand, the DIS entropy $S_{DIS}$ is generated by the separation of
the hadron state into parton states by the deep inelastic process. It is,
however, not a directly observable quantity, since the partons cannot escape
the hadron due to confinement. It also has the consequence that the partons at
a fixed value of the Bjorken variable $x$ are not directly countable but the
number depends on the renormalization scale $\mu^{2}$.

In our approach to the gluon density [3] the parton number at small $x$ is
given by the Pomeron intercept, see (5). At first sight this seems like a
contradiction, since the gluon intercept, at least for energies in and below
the TeV region (i.e., the LHC range and below), can be observed
experimentally, e.g., from forward scattering experiments. The solution to
this problem consists in the introduction of a scale dependence of the Pomeron
intercept, $\alpha_{P}(t=0,\mu^{2})$, as discussed extensively in Sec. II.3
and proposed earlier on phenomenological grounds in Refs. [7, 10] .

In both cases, $S_{BH}$ and $S_{DIS}$, the determination of the entropy is
based on theoretical arguments. A direct counting of the microstates is
forbidden: By the event horizon in black holes and by confinement in DIS.
Additional theoretical arguments, however, allow a relation to observable
quantities, by Hawking radiation and high energy scattering. It may give
further support to the idea of a possible relation between confinement and
black holes as may be realized in a holographic context [57].

###### Acknowledgements.

We thank Tianbo Liu for providing us the DGLAP results for the light-front
holographic model shown in Fig. 3. We also thank Ulli Marquard for drawing our
attention to the possible parallel of black-hole entropy. SJB is supported in
part by the Department of Energy Contract No. DE-AC02-76SF00515.

## References

  * Kharzeev and Levin [2017] D. E. Kharzeev and E. M. Levin, Deep inelastic scattering as a probe of entanglement, Phys. Rev. D 95, 114008 (2017), arXiv:1702.03489 [hep-ph] .
  * Kharzeev and Levin [2021] D. E. Kharzeev and E. Levin, Deep inelastic scattering as a probe of entanglement: Confronting experimental data, Phys. Rev. D 104, L031503 (2021), arXiv:2102.09773 [hep-ph] .
  * de Téramond  _et al._ [2021] G. F. de Téramond, H. G. Dosch, T. Liu, R. S. Sufian, S. J. Brodsky, and A. Deur (HLFHS), Gluon matter distribution in the proton and pion from extended holographic light-front QCD, Phys. Rev. D 104, 114005 (2021), arXiv:2107.01231 [hep-ph] .
  * Brodsky and de Téramond [2006] S. J. Brodsky and G. F. de Téramond, Hadronic spectra and light-front wave functions in holographic QCD, Phys. Rev. Lett. 96, 201601 (2006), arXiv:hep-ph/0602252 .
  * de Téramond and Brodsky [2009] G. F. de Téramond and S. J. Brodsky, Light-front holography: A first approximation to QCD, Phys. Rev. Lett. 102, 081601 (2009), arXiv:0809.4899 [hep-ph] .
  * Brodsky  _et al._ [2015] S. J. Brodsky, G. F. de Téramond, H. G. Dosch, and J. Erlich, Light-front holographic QCD and emerging confinement, Phys. Rept. 584, 1 (2015), arXiv:1407.8131 [hep-ph] .
  * Dosch  _et al._ [2022] H. G. Dosch, G. F. de Téramond, T. Liu, R. S. Sufian, S. J. Brodsky, and A. Deur (HLFHS), Towards a single scale-dependent Pomeron in holographic light-front QCD, Phys. Rev. D 105, 034029 (2022), arXiv:2201.09813 [hep-ph] .
  * Tu  _et al._ [2020] Z. Tu, D. E. Kharzeev, and T. Ullrich, Einstein-Podolsky-Rosen paradox and quantum entanglement at subnucleonic scales, Phys. Rev. Lett. 124, 062001 (2020), arXiv:1904.11974 [hep-ph] .
  * Hentschinski  _et al._ [2023] M. Hentschinski, D. E. Kharzeev, K. Kutak, and Z. Tu, Probing the onset of maximal entanglement inside the proton in diffractive deep inelastic scattering, Phys. Rev. Lett. 131, 241901 (2023), arXiv:2305.03069 [hep-ph] .
  * Dosch and Ferreira [2015] H. G. Dosch and E. Ferreira, Diffractive electromagnetic processes from a Regge point of view, Phys. Rev. D 92, 034002 (2015), arXiv:1507.03442 [hep-ph] .
  * Regge [1959] T. Regge, Introduction to complex orbital momenta, Nuovo Cim. 14, 951 (1959).
  * Chew and Frautschi [1961] G. F. Chew and S. C. Frautschi, Principle of equivalence for all strongly interacting particles within the $S$-matrix framework, Phys. Rev. Lett. 7, 394 (1961).
  * Zyla  _et al._ [2020] P. A. Zyla _et al._ (Particle Data Group), Review of Particle Physics, PTEP 2020, 083C01 (2020).
  * de Téramond  _et al._ [2018] G. F. de Téramond, T. Liu, R. S. Sufian, H. G. Dosch, S. J. Brodsky, and A. Deur (HLFHS), Universality of generalized parton distributions in light-front holographic QCD, Phys. Rev. Lett. 120, 182001 (2018), arXiv:1801.09154 [hep-ph] .
  * Froissart [1961] M. Froissart, Asymptotic behavior and subtractions in the Mandelstam representation, Phys. Rev. 123, 1053 (1961).
  * Martin [1965] A. Martin, Extension of the axiomatic analyticity domain of scattering amplitudes by unitarity-I., Nuovo Cim. A 42, 930 (1965).
  * Heisenberg [1952] W. Heisenberg, Mesonenerzeugung als stoßwellenproblem, Z. Physik 133, 65 (1952).
  * Altarelli and Parisi [1977] G. Altarelli and G. Parisi, Asymptotic freedom in parton language, Nucl. Phys. B 126, 298 (1977).
  * Dokshitzer [1977] Y. L. Dokshitzer, Calculation of the structure functions for deep inelastic scattering and $e^{+}e^{-}$ annihilation by perturbation theory in quantum chromodynamics, Sov. Phys. JETP 46, 641 (1977).
  * Abramowicz  _et al._ [2015] H. Abramowicz _et al._ (H1, ZEUS), Combination of measurements of inclusive deep inelastic ${e^{\pm}p}$ scattering cross sections and QCD analysis of HERA data, Eur. Phys. J. C 75, 580 (2015), arXiv:1506.06042 [hep-ex] .
  * Maldacena [1998] J. M. Maldacena, The Large-N limit of superconformal field theories and supergravity, Adv. Theor. Math. Phys. 2, 231 (1998), arXiv:hep-th/9711200 .
  * Gubser  _et al._ [1998] S. S. Gubser, I. R. Klebanov, and A. M. Polyakov, Gauge theory correlators from noncritical string theory, Phys. Lett. B 428, 105 (1998), arXiv:hep-th/9802109 .
  * Witten [1998] E. Witten, Anti-de Sitter space and holography, Adv. Theor. Math. Phys. 2, 253 (1998), arXiv:hep-th/9802150 .
  * Brower  _et al._ [2007] R. C. Brower, J. Polchinski, M. J. Strassler, and C.-I. Tan, The Pomeron and gauge/string duality, J. High Energy Phys. 12, 005, arXiv:hep-th/0603115 .
  * Liu  _et al._ [2020] T. Liu, R. S. Sufian, G. F. de Téramond, H. G. Dosch, S. J. Brodsky, and A. Deur (HLFHS), Unified description of polarized and unpolarized quark distributions in the proton, Phys. Rev. Lett. 124, 082003 (2020), arXiv:1909.13818 [hep-ph] .
  * Sufian  _et al._ [2018] R. S. Sufian, T. Liu, G. F. de Téramond, H. G. Dosch, S. J. Brodsky, A. Deur, M. T. Islam, and B.-Q. Ma, Nonperturbative strange-quark sea from lattice QCD, light-front holography, and meson-baryon fluctuation models, Phys. Rev. D 98, 114004 (2018), arXiv:1809.04975 [hep-ph] .
  * Sufian  _et al._ [2020] R. S. Sufian, T. Liu, A. Alexandru, S. J. Brodsky, G. F. de Téramond, H. G. Dosch, T. Draper, K.-F. Liu, and Y.-B. Yang, Constraints on charm-anticharm asymmetry in the nucleon from lattice QCD, Phys. Lett. B 808, 135633 (2020), arXiv:2003.01078 [hep-lat] .
  * Adloff  _et al._ [2001] C. Adloff _et al._ (H1), On the rise of the proton structure function $F_{2}$ towards low $x$, Phys. Lett. B 520, 183 (2001), arXiv:hep-ex/0108035 .
  * Abt  _et al._ [1993] I. Abt _et al._ (H1), Measurement of the proton structure function $F_{2}(x,Q^{2})$ in the low-$x$ region at HERA, Nucl. Phys. B 407, 515 (1993).
  * Derrick  _et al._ [1993] M. Derrick _et al._ (ZEUS), Measurement of the proton structure function $F_{2}$ in $ep$ scattering at HERA, Phys. Lett. B 316, 412 (1993).
  * Ioffe [1969] B. L. Ioffe, Space-time picture of photon and neutrino scattering and electroproduction cross-section asymptotics, Phys. Lett. B 30, 123 (1969).
  * Donnachie and Landshoff [1998] A. Donnachie and P. V. Landshoff, Small $x$: Two pomerons!, Phys. Lett. B 437, 408 (1998), arXiv:hep-ph/9806344 .
  * Kuraev  _et al._ [1977] E. A. Kuraev, L. N. Lipatov, and V. S. Fadin, The Pomeranchuk singularity in nonabelian gauge theories, Sov. Phys. JETP 45, 199 (1977).
  * Balitsky and Lipatov [1978] I. I. Balitsky and L. N. Lipatov, The Pomeranchuk singularity in quantum chromodynamics, Sov. J. Nucl. Phys. 28, 822 (1978).
  * Donnachie  _et al._ [2004] S. Donnachie, H. G. Dosch, O. Nachtmann, and P. Landshoff,  _Pomeron physics and QCD_ , Vol. 19 (Cambridge University Press, 2004).
  * Donnachie and Dosch [2002] A. Donnachie and H. G. Dosch, Comprehensive approach to structure functions, Phys. Rev. D 65, 014019 (2002), arXiv:hep-ph/0106169 .
  * Chwastowski and Figiel [2004] J. Chwastowski and J. Figiel, Photoproduction at HERA, Phys. Part. Nucl. 35, 619 (2004), arXiv:hep-ex/0311044 .
  * Aaij  _et al._ [2014] R. Aaij _et al._ (LHCb), Updated measurements of exclusive $J/\Psi$ and $\Psi$(2S) production cross-sections in $pp$ collisions at $\sqrt{s}=7$ TeV, J. Phys. G 41, 055002 (2014), arXiv:1401.3288 [hep-ex] .
  * Abelev  _et al._ [2014] B. B. Abelev _et al._ (ALICE), Exclusive $\mathrm{J/}\psi$ photoproduction off protons in ultraperipheral $p$-Pb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV, Phys. Rev. Lett. 113, 232504 (2014), arXiv:1406.7819 [nucl-ex] .
  * Aaij  _et al._ [2015] R. Aaij _et al._ (LHCb), Measurement of the exclusive $\Upsilon$ production cross-section in pp collisions at $\sqrt{s}=7$ TeV and 8 TeV, J. High Energy Phys. 09, 084, arXiv:1505.08139 [hep-ex] .
  * Gotsman and Levin [2020] E. Gotsman and E. Levin, High energy QCD: multiplicity distribution and entanglement entropy, Phys. Rev. D 102, 074008 (2020), arXiv:2006.11793 [hep-ph] .
  * Koba  _et al._ [1972] Z. Koba, H. B. Nielsen, and P. Olesen, Scaling of multiplicity distributions in high-energy hadron collisions, Nucl. Phys. B 40, 317 (1972).
  * Grosse-Oetringhaus and Reygers [2010] J. F. Grosse-Oetringhaus and K. Reygers, Charged-particle multiplicity in proton-proton collisions, J. Phys. G 37, 083001 (2010), arXiv:0912.0023 [hep-ex] .
  * Levin and Lublinsky [2004] E. Levin and M. Lublinsky, A linear evolution for non-linear dynamics and correlations in realistic nuclei, Nucl. Phys. A 730, 191 (2004), arXiv:hep-ph/0308279 .
  * Stoffers and Zahed [2013] A. Stoffers and I. Zahed, Holographic Pomeron and entropy, Phys. Rev. D 88, 025038 (2013), arXiv:1211.3077 [nucl-th] .
  * Germano and Navarra [2022] G. R. Germano and F. S. Navarra, Energy dependence of the multiplicity moments at the LHC, Phys. Rev. D 105, 014005 (2022), arXiv:2110.12028 [hep-ph] .
  * Bartels  _et al._ [2003] J. Bartels, E. Gotsman, E. Levin, M. Lublinsky, and U. Maor, The dipole picture and saturation in soft processes, Phys. Lett. B 556, 114 (2003), arXiv:hep-ph/0212284 .
  * Carvalho  _et al._ [2008] F. Carvalho, F. O. Duraes, V. P. Goncalves, and F. S. Navarra, Gluon saturation and the Froissart bound: A Simple approach, Mod. Phys. Lett. A 23, 2847 (2008), arXiv:0705.1842 [hep-ph] .
  * Khachatryan  _et al._ [2011] V. Khachatryan _et al._ (CMS), Charged particle multiplicities in $pp$ interactions at $\sqrt{s}=0.9$, 2.36, and 7 TeV, J. High Energy Phys. 01, 079, arXiv:1011.5531 [hep-ex] .
  * Acharya  _et al._ [2023] S. Acharya _et al._ (ALICE), Multiplicity dependence of charged-particle production in pp, p-Pb, Xe-Xe and Pb-Pb collisions at the LHC, Phys. Lett. B 845, 138110 (2023), arXiv:2211.15326 [nucl-ex] .
  * Adloff  _et al._ [1996] C. Adloff _et al._ (H1), Inclusive $D^{0}$ and $D^{*\pm}$ production in deep inelastic $ep$ scattering at HERA, Z. Phys. C 72, 593 (1996), arXiv:hep-ex/9607012 .
  * Breitweg  _et al._ [2000] J. Breitweg _et al._ (ZEUS), Measurement of $D^{*\pm}$ production and the charm contribution to $F_{2}$ in deep inelastic scattering at HERA, Eur. Phys. J. C 12, 35 (2000), arXiv:hep-ex/9908012 .
  * Ademollo and Del Giudice [1969] M. Ademollo and E. Del Giudice, Nonstrong amplitudes in a Veneziano-type model, Nuovo Cim. A 63, 639 (1969).
  * Landshoff and Polkinghorne [1970] P. V. Landshoff and J. C. Polkinghorne, The scaling law for deep inelastic scattering in a new Veneziano-like amplitude, Nucl. Phys. B 19, 432 (1970).
  * Bekenstein [1972] J. D. Bekenstein, Black holes and the second law, Lett. Nuovo Cim. 4, 737 (1972).
  * Hawking [1975] S. W. Hawking, Particle creation by black holes, Commun. Math. Phys. 43, 199 (1975), [Erratum: Commun. Math. Phys. 46, 206 (1976)].
  * Witten [2001] E. Witten, Black holes and quark confinement, Available at https://www.ias.edu/sites/default/files/sns/CurrentScienceVol81(4).pdf (2001).