\section{Introduction}

Reservoir Computing (RC) \cite{Lukosevicius2009,Verstraeten2007experimental} delineates a class of Recurrent Neural Network (RNN) models based on the idea of separating the non-linear dynamical component of the network, i.e. the recurrent hidden reservoir layer, from the feed-forward linear readout layer.
The reservoir is initialized randomly under stability constraints and then is left untrained, leaving the burden of training to fall only on the readout part of the architecture, hence resulting in a strikingly efficient approach to RNN design.
In this context, the Echo State Network (ESN) model \cite{Jaeger2004,Jaeger2001} is a popular realization of the RC paradigm based on implementing the reservoir in terms of a discrete-time non-linear dynamical system.
Being featured by untrained dynamics, ESNs represent an important tool to understand and characterize the operation and potentialities of recurrent neural models. Shaping the reservoir architecture in order to achieve desired properties and  optimized performance in applications, even in the absence of training of the recurrent connections,
is
one of the key goals of RC research \cite{Gallicchio2018randomized}.

In this paper we bring together two major trends in the area of ESN architectural studies.
The first one focuses on the pattern of connectivity among the recurrent units. In this case, the aim is to constrain the  random reservoir initialization process towards topologies that determine specific algebraic properties of the resulting recurrent weight matrices. A relevant class of reservoir variants in this regard is given by ESNs with orthogonal recurrent matrices \cite{White2004,Hajnal2006}, which were shown to lead to improved performance with respect to random reservoirs both in terms of memorization skills and in terms of predictive performance on non-linear tasks. In particular, reservoirs whose structure is based on permutation matrices represent particularly appealing instances of orthogonal ESNs \cite{Hajnal2006,Strauss2012}, entailing a simple and very sparse pattern of connectivity among the recurrent units.
Other relevant architectural variants are given by reservoirs structured according to a ring topology or to form a chain of units \cite{Rodan2011,Strauss2012}.
The second major  line of research that we consider regards the construction of hierarchically structured reservoir models. While initial studies in this context focused on composing multiple ESN modules
to form ad-hoc architectures
\cite{jaeger2007discovering,triefenbach2010phoneme}, recent works started analyzing the effects of stacking multiple untrained reservoir layers with the introduction of the DeepESN model \cite{Gallicchio2017DeepESN}.
On the one hand, the analysis of DeepESN dynamics contributes to uncover the intrinsic computational properties of deep neural networks in the temporal domain \cite{Gallicchio2017DeepESN,Gallicchio2018local}. On the other hand, a proper architectural design of deep reservoirs might have a huge impact in real-world applications \cite{Gallicchio2018design}, enabling effective multiple time-scales processing and at the same time preserving the training efficiency typical of RC models.

In this paper we analyze
the impact on the predictive performance given by a constrained
reservoir topology in DeepESNs. Specifically,
we consider deep architectures in which the individual reservoir layers are implemented based on permutation matrices, as well as on ring and on chain topologies. Our study is conducted in comparison to shallow ESN counterparts through numerical experiments on several benchmarks in the RC area.

The rest of this paper is structured as follows.
The DeepESN model is introduced in Section~\ref{sec.DeepESN}, while the investigated reservoir topologies
are described in Section~\ref{sec.topology}. The experimental analysis
is reported in Section~\ref{sec.experiments}. Finally, Section~\ref{sec.conclusions} draws conclusions and delineates future research directions.

\section{Deep Echo State Networks}
\label{sec.DeepESN}
A DeepESN is an RC model in which the reservoir part is organized into a stacked composition of multiple untrained recurrent hidden layers. The external input is propagated only to the first reservoir layer, while each successive level in the deep architecture is fed by the output of the previous one, as graphically illustrated in Figure~\ref{fig.deepreservoir}.

\begin{figure}[bth]
  \centering
	  \includegraphics[width=1\textwidth]{reservoir.eps}
	  \caption{Hierarchical reservoir architecture in a DeepESN.}
  \label{fig.deepreservoir}
\end{figure}

To fix our notation,  we use $L$ to indicate the number of layers in the deep reservoir, while $N_U$ and $N_Y$ respectively denote the sizes of input and output spaces.
For the sake of simplicity in the presentation of the DeepESN model, here we make the assumption that all the reservoir layers are featured by the same number of units, indicated by $N_R$.
The operation of each reservoir layer can be described in terms of a discrete-time non-linear dynamical system,
whose state update equation is given in the form of an iterated mapping. In particular, at time-step $t$, the state of the first layer, i.e. $\xl{1}(t) \in \mathbb{R}^{N_R}$, is computed as follows:
\begin{equation}
\label{eq.layer1}
\xl{1}(t) = \tanh(\W_{in} \u(t) + \hat{\W}^{(1)} \xl{1}(t-1)),
\end{equation}
while the state of each successive layer $l>1$, i.e. $\xl{l}(t) \in \mathbb{R}^{N_R}$, is given by:
\begin{equation}
\label{eq.layerl}
\xl{l}(t) = \tanh(\mathbf{W}^{(l)} \xl{l-1}(t) + \hat{\W}^{(l)} \xl{l}(t-1)).
\end{equation}
Here, $\tanh$ indicates the element-wise application of the hyperbolic tangent non-linearity, $\u(t) \in \mathbb{R}^{N_U}$ represents the external input at time-step $t$, while $\W_{in}$, $\mathbf{W}^{(l)}$ and $\hat{\W}^{(l)}$ respectively denote the input weight matrix (that modulates the external input stimulation to the first layer), the inter-layer weight matrix for layer $l$ (that modulates the strength of the connections from layer $l-1$ to layer $l$), and the recurrent reservoir weight matrix for layer $l$. In both the above equations~\ref{eq.layer1} and \ref{eq.layerl} we omitted the bias terms for the ease of notation. The interested reader can find in \cite{Gallicchio2017DeepESN} a more detailed description of the deep reservoir equations, framed in the more general context of leaky integrator reservoir units. In order to set up initial conditions for the state update
equations~\ref{eq.layer1} and \ref{eq.layerl}, at time-step $0$ all reservoir layers are set to a null state, i.e. $\mathbf{x}^{(l)}(0) = \mathbf{0}$ for all $l = 1,\ldots, L$.
Given this framework, it is worth noticing that a standard shallow ESN model can be seen as a special case of DeepESN in which a single reservoir layer is considered, i.e. $L = 1$.

As in standard RC approaches, the parameters of the entire reservoir component, i.e. the elements in all the weight matrices in equations~\ref{eq.layer1} and \ref{eq.layerl}, are left untrained after initialization subject to stability constraints. These are required to avoid the system dynamics to fall into unstable regimes, which would make them unsuitable for robust processing of time-series data. In the context of ESNs, the analysis of asymptotic stability is usually described in terms of the Echo State Property (ESP) \cite{Jaeger2001,Lukosevicius2009}, providing simple algebraic conditions for the initialization of reservoir weight matrices that have been recently extended to cope with the case of deep reservoirs in \cite{Gallicchio2017echo}.
Under a practical view-point, the analysis in \cite{Gallicchio2017echo} suggests to carefully control the spectral radius
of all the reservoir weight matrices in the deep reservoir. In this paper, we use $\rho^{(l)}$ to denote the spectral radius in layer $l$, i.e. the largest among the absolute values of the eigenvalues of $\hat{\W}^{(l)}$.
A simple initialization procedure for the reservoir of a DeepESN then consists in
choosing the elements in $\hat{\W}^{(l)}$ randomly from a uniform distribution on $[-1,1]$, subsequently re-scaling them to
achieve desired values of $\rho^{(l)}$, typically not above unity. Similarly, the elements in $\W_{in}$ and those in $\mathbf{W}^{(l)}$ (for $l>1$) are initialized randomly from a uniform distribution on $[-1,1]$, and then are re-scaled to control
the input scaling hyper-parameter $\omega_{in} = \|\W_{in}\|_2$, and the set of inter-layer scaling hyper-parameters
$\omega_{il}^{(l)} = \|\mathbf{W}^{(l)}\|_2$.

The output of the DeepESN is computed by a simple readout tool, which linearly combines the reservoir representations developed in all the layers of the deep architecture. In formulas, the output at time-step $t$, denoted as $\mathbf{y}(t) \in \mathbb{R}^{N_Y}$, is computed by the following equation:
\begin{equation}
\label{eq.readout}
\mathbf{y}(t) = \W_{out}  \, [\mathbf{x}^{(1)}(t); \ldots ; \mathbf{x}^{(L)}(t)],
\end{equation}
where $\W_{out}$ is the output weight matrix, and $[\mathbf{x}^{(1)}(t); \ldots ; \mathbf{x}^{(L)}(t)]$ represents the global deep reservoir state at time-step $t$, expressed as the concatenation of all the states in the architecture.
The elements in $\W_{out}$ represent the only learnable weights of the DeepESN, and are typically adjusted to fit a training set by exploiting non-iterative training algorithms as in the case of standard RC models \cite{Lukosevicius2009}.
Notice that, although different patterns of reservoir-to-readout connectivity are possible \cite{Pascanu2014}, the one employed here, where all reservoir layers are used to feed the readout, has the advantage to allow the training algorithms to modulate and exploit differently the variety of representations provided by the different levels in the deep architecture.

A more comprehensive description of the characteristics and advantages of the DeepESN approach can be found in \cite{Gallicchio2018Layering}, while a constantly updated overview on the advancements achieved in this research field is given in \cite{Gallicchio2017survey}.
To date, software implementations of the DeepESN model are made publicly available as libraries for
Python\footnote{\url{https://github.com/lucapedrelli/DeepESN}},
Matlab\footnote{\url{https://it.mathworks.com/matlabcentral/fileexchange/69402-deepesn}} and
Octave\footnote{\url{https://github.com/gallicch/DeepESN_octave}}.

\section{Reservoir Topology}
\label{sec.topology}
We consider DeepESN architectural variants where the recurrent weight matrix in each layer $l$, i.e. $\hat{\W}^{(l)}$, is characterized by a specific structure, according to the topologies described in the following.
The resulting patterns of reservoir connectivity are graphically exemplified in Figure~\ref{fig.topology}.

\begin{figure}[bth]
  \centering
	  \includegraphics[width=0.9\textwidth]{topology.eps}
	  \caption{Reservoir topologies of DeepESN layers.}
  \label{fig.topology}
\end{figure}

\begin{description}
\item[Sparse] Each reservoir unit is randomly connected to a subset of the others, determining a \emph{sparse} recurrent matrix $\hat{\W}^{(l)}$ (see Figure~\ref{fig.topology}(a)). This corresponds to a common setting used in RC practice and serves here as baseline for our analysis.
\item[Permutation] The structure of the recurrence matrix $\hat{\W}^{(l)}$ is given by a \emph{permutation} matrix $\mathbf{P}$, i.e. we have:
\begin{equation}
\label{eq.permutation}
\hat{\W}^{(l)} = \lambda \, \mathbf{P},
\end{equation}
where $\mathbf{P}$ is obtained by randomly permuting the columns of the identity matrix, and $\lambda$ is a multiplicative constant that specifies the value of the non-zero recurrent weights. In this case, the spectral radius of $\Wh^{(l)}$ is determined by the value of $\lambda$, i.e. $\rho^{(l)} = \lambda$.
The permutation  topology implies that each row and each column of the recurrence matrix have exactly one non-zero element, resulting into a reservoir architecture that presents a variable number of disjoint cyclic structures, as graphically exemplified in Figure~\ref{fig.topology}(b).
The levels in the deep reservoir architecture are allowed to employ different permutations, i.e. the number of cycles in each reservoir layer can be different.

In the context of shallow ESNs, this kind of topology has been empirically studied
in \cite{Boedecker2009}, where it was shown to achieve good memorization skills at the same time improving the performance of randomly initialized reservoirs in tasks involving non-linear mappings. Interestingly, the permutation topology has been investigated in \cite{Hajnal2006} as a way to implement orthogonal reservoir matrix structures, under the name of Critical ESNs.

\item[Ring] The reservoir units  are organized to form a single \emph{ring}, as shown in Figure~\ref{fig.topology}(c). Accordingly, the recurrent weight matrix $\Wh^{(l)}$ is expressed as:
\begin{equation}
\label{eq.ring}
\Wh^{(l)} = \lambda \,
\left[
\begin{array}{c c c c}
0 & 0 & \ldots & 1\\
1 & 0 & \ldots & 0\\
\vdots & \ddots & \ddots & \vdots\\
0 & \ldots & 1 & 0\\
\end{array}
\right],
\end{equation}
where $\lambda$ is the value of non-zero recurrent weights, and determines the spectral radius of $\Wh^{(l)}$, i.e. $\rho^{(l)} = \lambda$. The ring topology can be easily seen as a special case of the permutation topology, where the pattern of reservoir connectivity is ruled by the specific permutation matrix in equation~\ref{eq.ring}, and the reservoir units are all part of the same cyclic structure.

Reservoirs following this architectural organization have been subject of several studies in literature on shallow RC. Notable instances in this regard are given by the work in \cite{Strauss2012}, in which the ring topology is studied in the context of orthogonal reservoir structures, and by the work in \cite{Rodan2011},
where the study is carried out under the perspective of architectural design simplification for minimum complexity ESN construction. One interesting outcome of previous analysis on the ring topology is that, compared to randomly initialized reservoirs, it shows superior memory capacity that, at least in the linear case, approaches the optimal value \cite{Rodan2011}.
While this optimal memory characterization has been extensively analyzed in literature for the more general class of orthogonal recurrent weight matrices (see e.g. \cite{White2004,Jaeger2001short,Farkas2016}), the ring topology presents the advantage of a strikingly simple (and sparse) dynamical network construction.

\item[Chain] The recurrent units are arranged in a pipeline, where each unit - except for the first one - receives in input the activation of the previous one, forming a \emph{chain} as in the example in Figure~\ref{fig.topology}(d).
The only non-zero elements in $\Wh^{(l)}$ are located in the lower sub-diagonal, i.e. we have:
\begin{equation}
\label{eq.chain}
\Wh^{(l)} = \lambda \,
\left[
\begin{array}{c c c c}
0 & 0 & \ldots & 0\\
1 & 0 & \ldots & 0\\
\vdots & \ddots & \ddots & \vdots\\
0 & \ldots & 1 & 0\\
\end{array}
\right],
\end{equation}
where as in previous cases $\lambda$ identifies the value of non-zero weights.
Although in this case $\Wh^{(l)}$ is nilpotent and hence its spectral radius is always $0$, we still operate on
$\lambda$ to control the magnitude of recurrent weights. As such, with a slight abuse of notation, also in this case we set $\rho^{(l)} = \lambda$.
Overall, the chain topology results in a particularly simple design strategy that, from the architectural perspective, applies a further simplification to the ring topology by removing one of the connections between the internal units.

Literature works on shallow ESN models pointed out the merits of reservoir organizations based on a chain topology (also called delay-line reservoirs), as a very simple approach to the architectural design of the network, resulting in a model that is easier to analyze \cite{White2004} and that leads to comparable or even better performance than standard ESNs \cite{Rodan2011,Strauss2012}.
\end{description}

\section{Experiments}
\label{sec.experiments}
In this section we illustrate the experimental analysis conducted in this paper.
Specifically, in Section~\ref{sec.settings} we detail the datasets considered and the experimental settings
adopted in our work, whereas in Section~\ref{sec.results} we report and discuss the achieved numerical results.

\subsection{Datasets and Experimental Settings}
\label{sec.settings}

In our experiments, we considered benchmark datasets featured by univariate time-series (i.e., $N_U = N_Y = 1$).

The first dataset is obtained from a non-linear auto-regressive moving average system of the 10-th order (NARMA10).
At each time-step, the input $u(t)$ comes from a uniform distribution over $[0, 0.5]$, whereas the corresponding target output
$y_{tg}(t)$ is given by the following relation:
\begin{equation}
\label{eq.narma10}
y_{tg}(t) = 0.3 \, y_{tg}(t-1)+0.05 \, y_{tg}(t-1) \sum\limits_{i = 1}^{10} y_{tg}(t-i)
                 + 1.5 \, u(t-10) \, u(t-1)+0.1.
\end{equation}

The second dataset that we considered is the Santa Fe Laser time-series \cite{Weigend2018time}, where the input values $u(t)$
are sampled intensities from a far-infrared laser in chaotic regime, re-scaled by a factor of 0.01.
We used the Laser dataset to define a next-step prediction task, where $y_{tg}(t) = u(t+1)$ for each time-step $t$.

The last two datasets are instances of the Mackey-Glass (MG) \cite{Mackey1977oscillation,Farmer1982chaotic} time-series, obtained by discretizing the following non-linear differential equation:
\begin{equation}
\label{eq.mg}
\frac{\delta u(t)}{\delta t} = \frac{0.2 \, u(t-\tau)}{1+u(t-\tau)^{10}} - 0.1 \, u(t),
\end{equation}
where $\tau$ is a parameter of the system influencing its dynamical behavior.
We generated two MG time-series using $\tau = 17$ (MG17) and $\tau = 30$ (MG30), representing cases with increasingly complex chaotic behavior. In both cases, the elements of the time-series where shifted by -1 and then passed through the $\tanh$ squashing function as in \cite{Jaeger2004,Jaeger2001}. The two MG time-series allowed us to set up two next-step prediction tasks, where $y_{tg}(t) = u(t+1)$ for each time-step $t$.

For NARMA10, MG17 and MG30 we generated datasets with 10000 time-steps, while the Laser dataset contained a number of 10092 samples. In all the cases, the available data was split into a training set, comprising the first 5000 time-steps, and a test set, comprising the remaining samples. The first 100 time-steps were used as transient to wash out the initial conditions. The performance of the considered RC models was evaluated in terms of mean squared error (MSE) in all the tasks.

In our experiments, we considered DeepESNs with a total number of $500$ recurrent reservoir units,
distributed evenly across the layers of the deep architecture, varying the number of layers $L$ from 2 to 5\footnote{With the only exception of the case $L=3$, where the first two layers contained 167 units and the last one contained 166 units.}.
All the reservoir layers in the deep architecture shared the same values for the scaling hyper-parameters $\rho$ and $\omega_{il}$, i.e. $\rho = \rho^{(1)} = \ldots = \rho^{(L)}$ and $\omega_{il} = \omega_{il}^{(2)} = \ldots = \omega_{il}^{(L)}$.
To account for sparsity, each reservoir unit was randomly connected to 5 units in the previous layer and to 5 units in the same layer. Of course, when considering permutation, ring and chain reservoir topologies, the connectivity of the reservoir units in each layer followed the corresponding specific structure described in Section~\ref{sec.topology}. In all the cases, we used fully-connected input weight matrices.
For every task and choice of the reservoir topology, the DeepESN hyper-parametrization was chosen by model selection on a validation set comprising the last 1000 time-steps of the training split. To this end, we performed a random search with 50 networks configurations for each number of layers, sampling the value of $\rho$ from a uniform distribution in $(0.1, 1]$, and the values of $\sin$ and $\omega_{il}$ from uniform distributions in $(0.1, 2]$. The achieved results were averaged on 10 network guesses for each hyper-parametrization explored, and readout training was performed by using pseudo-inversion.
Finally, our experimental analysis was conducted in comparison with shallow ESN setups, considering the same reservoir topologies investigated in the DeepESN case, and using the same experimental setting described above, with the only crucial exception that all the available reservoir units were organized into a single layer (i.e. $L = 1$).

\subsection{Results}
\label{sec.results}
The test MSE values obtained by DeepESNs in correspondence of all the considered types of layer-wise reservoir topology are reported in Table~\ref{tab.results}. For the sake of comparison, the same table shows the results achieved by shallow ESNs under the same architectural conditions examined in the deep case. In all the cases, the sparse reservoir topology is considered as a baseline setup for our analysis.

\begin{table}[tb]
  \centering
  \begin{tabular}{l l  l l}
    \multicolumn{4}{c}{\textbf{NARMA10}}  \\
    \hline
    %
    \textbf{Topology} $\quad$ & \textbf{ESN}  & \textbf{DeepESN} & Layers\\
  	Sparse & $1.658 \; 10^{-4} \, (3.367 \; 10^{-5}) \quad$  & $1.647 \; 10^{-4} \,(3.415 \; 10^{-5})$ & 2\\
    Permutation & $1.354 \; 10^{-4} \, (1.589 \; 10^{-5})$  & $\underline{1.243 \; 10^{-4}} \,(1.464 \; 10^{-5})$ & 2\\
    Ring & $1.494 \; 10^{-4} \, (1.547 \; 10^{-5})$  & $1.482 \; 10^{-4} \,(1.713 \; 10^{-5})$ & 2\\
    Chain & $1.571 \; 10^{-4} \, (1.780 \; 10^{-5})$  & $1.569 \; 10^{-4} \,(2.594 \; 10^{-5})$ & 2\\
    \\
  \end{tabular}
  %
  \vspace{3mm}
  %
  \begin{tabular}{l l  l l}
    \multicolumn{4}{c}{\textbf{Laser}}  \\
    \hline
    %
    \textbf{Topology} $\quad$ & \textbf{ESN}  & \textbf{DeepESN} & Layers\\
  	Sparse & $1.226 \; 10^{-3} \, (1.037 \; 10^{-4}) \quad$  & $8.228 \; 10^{-4} \,(2.309 \; 10^{-4})$ & 5\\
    Permutation & $1.312 \; 10^{-3} \, (1.385 \; 10^{-4})$  & $\underline{6.633 \; 10^{-4}} \,(8.861 \; 10^{-5})$ & 5\\
    Ring & $1.161 \; 10^{-3} \, (7.541 \; 10^{-5})$  & $7.640 \; 10^{-4} \,(4.331 \; 10^{-5})$ & 4\\
    Chain & $9.496 \; 10^{-4} \, (1.183 \; 10^{-4})$  & $8.555 \; 10^{-4} \,(8.302 \; 10^{-5})$ & 4\\

  \end{tabular}

  \vspace{3mm}

  \begin{tabular}{l l  l l}
    \multicolumn{4}{c}{\textbf{MG17}}  \\
    \hline
    %
    \textbf{Topology} $\quad$ & \textbf{ESN}  & \textbf{DeepESN} & Layers\\
  	Sparse & $3.739 \; 10^{-9} \, (1.387 \; 10^{-9}) \quad$  & $2.328 \; 10^{-9} \,(8.299 \; 10^{-10})$ & 2\\
    Permutation & $3.093 \; 10^{-9} \, (3.241 \; 10^{-10})$  & $\underline{4.576\; 10^{-10}} \,(6.280 \; 10^{-10})$ & 5\\
    Ring & $1.585 \; 10^{-9} \, (2.989 \; 10^{-10})$  & $5.043 \; 10^{-10} \,(3.891 \; 10^{-10})$ & 5\\
    Chain & $1.950 \; 10^{-9} \, (3.745 \; 10^{-10})$  & $4.913 \; 10^{-10} \,(2.535 \; 10^{-10})$ & 3\\
    \\
  \end{tabular}
  %
  \vspace{3mm}

  \begin{tabular}{l l  l l}
    \multicolumn{4}{c}{\textbf{MG30}}  \\
    \hline

    \textbf{Topology} $\quad$ & \textbf{ESN}  & \textbf{DeepESN} & Layers\\
  	Sparse & $1.476 \; 10^{-8} \, (1.781 \; 10^{-9}) \quad$  & $1.172 \; 10^{-8} \,(1.406 \; 10^{-9})$ & 2\\
    Permutation & $1.027 \; 10^{-8} \, (5.412 \; 10^{-10})$  & $\underline{8.618\; 10^{-9}} \,(1.457 \; 10^{-9})$ & 3\\
    Ring & $1.197 \; 10^{-8} \, (1.549 \; 10^{-9})$  & $1.078 \; 10^{-8} \,(2.066 \; 10^{-9})$ & 5\\
    Chain & $1.086 \; 10^{-8} \, (9.519 \; 10^{-10})$  & $9.096 \; 10^{-9} \,(1.803 \; 10^{-9})$ & 3\\

    \\
  \end{tabular}
  \caption{Test MSE (and std) achieved by shallow ESN and DeepESN settings
  for different choices of the reservoir topology.
  The last column reports the number of layers selected for DeepESN.
  Best results for each task are underlined.
  }
\label{tab.results}
\end{table}

The performance values reported in Table~\ref{tab.results} allow us to draw several lines of observations.
First of all, our results confirm the goodness of the considered reservoir architectural variants already in the shallow setup,
showing improved performance (i.e., a smaller MSE) with respect to the sparse baseline in all the cases (with the sole exception of permutation shallow reservoirs on Laser).
Second,
we observe that the performance of DeepESN with constrained topology (i.e. permutation, ring and chain) enhances that one of sparse DeepESN in all the considered tasks (with the only exception of deep reservoirs with chain architecture on Laser). Moreover, we can see that DeepESN improves the results of shallow ESN in all the tasks and for all the choices of reservoir topology, both in the constrained architectural cases and for the base sparse reservoir setup.
Taken together, results in Table~\ref{tab.results} clearly indicate the performance advantage arising from the synergy between deep organization and constrained topology as factors of architectural design of reservoirs.
Giving a structure to the architecture of reservoirs both at a coarser level, i.e. organizing the recurrent units into layers, and at a finer level, i.e. organizing individual layers' units into cyclic or chain structures,
amplifies the benefits brought by the two factors individually.
Finally, we notice that the best performing architecture in our experiments is the DeepESN with permutation reservoir topology, which obtained the smallest errors on all the tasks, and is put forward here as a particularly effective (yet sparse and efficient) approach to the architectural design of reservoir models.

\section{Conclusions}
\label{sec.conclusions}
In this paper we have investigated the role of reservoir topology in the architectural design of DeepESNs.
Specifically, we focused on analyzing the effects of constraining the recurrent weight matrix of each layer according to permutation, ring and chain topologies.
Numerical results on several RC benchmarks pointed out a striking beneficial effect arising from the combination of a deep reservoir construction with a structured organization of the recurrent units in each layer. Our results indicate that DeepESN with reservoir units arranged to obey a permutation scheme (i.e., forming multiple rings) provides a particularly advantageous design strategy for reservoirs, leading to the best performance in all the explored tasks.

While already giving interesting empirical evidences on the potentialities of deep RC architectures, the study presented in this paper opens the way to several directions for further research.
First of all, the experimental analysis described here suggests that the use of simplified deep RC models has a great potential that can be exploited massively in real-world applications. Leveraging the parsimonious design approach resulting from structured sparsity of reservoir units, the class of deep neural models studied in this work seems an ideal candidate e.g. for embedding advanced learning capabilities on resource-constrained computing devices.
On the methodological side, a natural extension of the work in this paper is to analyze the effect of a broader pool of reservoir architectural variants, including e.g. small-world \cite{Kawai2019small}, cycles with regular jumps \cite{Rodan2012} and concentric \cite{bacciu2018concentric} reservoirs.
Moreover, future research could pursue even further the simplification of architectural construction in deep RC models, reducing the impact of randomness in the network initialization in the same vein as the works on minimum complexity ESNs \cite{Rodan2011,Rodan2012}. Simplifying the reservoir structure locally to each layer can also be exploited from a more theoretically-oriented perspective, easing the mathematical analysis of dynamical properties naturally emerging in deep RNNs. In this concern, it is certainly interesting to extend fundamental mathematical results,
e.g. pertaining to short-term memory capacity \cite{Rodan2011,White2004,Jaeger2001short}, or
to  approximation properties \cite{grigoryeva2018echo} of shallow reservoirs to the case of DeepESN.
In addition to this, we believe that the role of orthogonality in deep reservoirs, studied in this paper in relation to the individual layers of the architecture, is an intriguing concept that deserves to be investigated also at the level of global (instead of local) DeepESN dynamics. Finally, the advantages of constrained DeepESN architectures delineated in this paper can be extended to larger classes of models, including e.g. deep RC for complex data structures \cite{gallicchio2019deep}, as well as fully trained deep RNNs.

\bibliographystyle{splncs04}


\section{Introduction}
\label{intro}

Dense matter is customarily described using some variation of the mean field theory approach (MFT) based on the extensive Boltzmann-Gibbs statistics (BS) (see \cite{Symmetry} for a recent review and references). However, in most cases the systems considered are not extensive because there are phenomena like, for example, long-range correlations or intrinsic fluctuations, not accounted for by the MFT used; in such cases the conditions justifying the use of BG statistics are not fulfilled. The usual remedy in such cases is to enrich the original dynamics by adding to the initial MFT some new elements accounting for these factors. The other possibility is to keep the original MFT intact, but to replace the BG statistics by its nonextensive counterpart which, by definition has these factors built in; usually by the Tsallis statistics (TS) \cite{PSA,PSA1,CSRH,JRGW1,JRGW2,JRGW3}. It is characterized by a nonextensivity parameter $q\neq 1$ such that, for $q \to 1$ one returns to the extensive situation of BG statistics. In such an approach modifications caused by introduction of nonextensive statistics are supposed to sum up the actions of all factors violating extensivity, both dynamical and caused by the environment. It is therefore expected that when these factors are gradually identified and their impact is accounted for by a suitable modification of the dynamics of the original model, the $q$ needed to fit data gradually tends to unity and one recovers the usual extensive situation of the BS, albeit now used with a modified MFT \cite{MB}. Note that this means that there is no such thing as a non-extensive free particle, and this is an inherent dynamical feature of  nonextensivity \cite{Biro-1,Biro-2,Biro-3}. To some extent the nonextensive statistics replaces the nonextensive dynamics (and vice versa), therefore in many cases we can talk interchangeably about {\it dynamical nonextensivity} or {\it nonextensive dynamics}. In this paper we shall take a closer look at this feature. To facilitate this task we limit ourselves to the simplest possible implementation of the dynamics in the form of some specific quasi-particle model of interactions proposed in \cite{QPM,QPMu} in which interactions are modelled by only one (albeit temperature dependent) parameter for each type of particle (here quarks, $q$, and gluons, $g$), and the effective fugacities,$z_{i=q,g}$, cf. \cite{Symmetry} for more references). In such an approach a single nonextensivity parameter, $q$, can be confronted with a set of individual dynamical parameters, $z_{i=q,g}$.

\section{The quasi-particle model in extensive and nonextensive environments}
\label{sec:form-qz}

For completeness of presentation we start with a short reminder of the essentials of the quasi-particle model in extensive, $z$-QPM, and non-extensive, $qz$-QPM, environments (cf. \cite{Symmetry} for details and further references).

In the case of an extensive environment the $z$-QPM \cite{QPM} is based on the effective equilibrium distribution function for quasi-partons ($i=q,~s,~g$ for, respectively, $u$ and $d$ massless quarks, strange quarks of mass $m$ and gluons; $e(x) = \exp(x)$, $\xi = +1$ for bosons and $-1$ for fermions and $\beta = 1/T$):
\begin{eqnarray}
n\left[ x^{(i)}\right] &=& \frac{1}{\frac{1}{z^{(i)}}e\left[ x^{(i)} \right] - \xi} = \frac{1}{e\left[ \tilde{x}^{(i)}\right] - \xi}, \label{fex-d}\\
x^{(i)} &=& \left\{
 \begin{array}{lll}
  \beta \left[ E_i - \mu^{(i)}\right]~~ &~~ {\rm if}&~~ i=q, s,\\
  \beta E_i ~~&~~ {\rm if} &~~ i= g.
  \end{array}\right.\qquad {\rm and}\qquad \tilde{x}^{(i)} = x^{(i)} - \ln z^{(i)}(\tau),
  \label{x(i)}
\end{eqnarray}
$E_{i=q,g}=p$ and $E_s =\sqrt{m^2 + p^2}$. Note that $e(x)\cdot e(-x) =1$, a consequence of which is that $n(x) + n(-x) = \xi$. The $z^{(i)} \leq 1$ denote the effective fugacities which describe the interactions and, by assumption, depend only on the scaled temperature, $\tau = T/T_{c}$ ($T_{c}$ is the temperature of transition to the deconfined phase of QCD). They are obtained by fits to the lattice QCD data, for  $z^{(i)} = 1$ one has free particles. In its original version $z$-QPM is formulated for zero chemical potential $\mu^{(i)}$ \cite{QPM} (reflecting difficulties with accounting for it in lattice QCD calculations; nonzero $\mu$, not connected with $z$, was introduced recently in \cite{QPMu} anticipating expected future developments in lattice QCD calculations).

In the case of $qz$-QPM one takes some {\it extensive system of quasi-particles}, the interactions of which are described by fugacities $z^{(i)}$ (and, possibly, also by some chemical potentials $\mu^{(i)}$, not connected with these fugacities), and immerses it in a nonextensive environment characterised by a nonextensivity parameter $q\neq 1$. Note that, as mentioned before, this means that we are dealing now with particles which are not really free, even when the dynamics is switched off. As a result one gets a {\it nonextensive system of interacting quasi-particles} and, assuming that the external dynamical information encoded in the results of the lattice QCD simulations remains intact, one has to find a new set of fugacities, $z^{(i)}_q$, which, together with modifications in the distribution caused by $q\neq 1$, will reproduce this information. Technically speaking, one simply replaces in Eq. (\ref{fex-d}) the exponential function $e(x) = \exp(x)$ by its nonextensive equivalent, the nonextensive exponent $e_q(x) = \exp_q(x)$, and its dual, $e_{2-q}(-x) = \exp_{2-q}(-x)$, defined as:
\begin{equation}
e_q(x) = [ 1 + (q-1)x ]^{\frac{1}{q-1}},\qquad e_{2-q}(-x) = [1 + (1-q)(-x)]^{\frac{1}{1-q}}, \label{eqe}
\end{equation}
$e_q(x) \to e(x)$ and $e_{2-q}(-x) \to e(-x)$ for $q\to 1$, the correspondig $q$ and $(2-q)$-logarithm functions are:
\begin{equation}
\ln_q X = \frac{X^{q-1} - 1}{q - 1} \stackrel{q \rightarrow 1}{\Longrightarrow} \ln X\quad {\rm and}\quad \ln_{2-q}X = \frac{X^{1-q}-1}{1-q} \stackrel{q \rightarrow 1}{\Longrightarrow} \ln X.  \label{lnq}
\end{equation}
As a result one obtains the following nonextensive particle occupation numbers:
\begin{equation}
n_q\left[ \tilde{x}^{(i)} \right] = \frac{1}{ e_q\left[ \tilde{x}^{(i)}\right] - \xi} = \frac{e_{2-q}\left[- \tilde{x}^{(i)}\right]}{1 - \xi e_{2-q}\left[- \tilde{x}^{(i)}\right]}, \label{enq}
\end{equation}
with $\tilde{x}^{(i)}$ given by Eq. (\ref{x(i)}) (with energy, $E_i$, and chemical potential, $\mu$, remaining unchanged). Note that now  $e_q(x)\cdot e_{2-q} =1$ and, respectively, $n_q(x) + n_{2-q}(-x) = \xi$. To preserve thermodynamic consistency in the nonextensive environment, one has to use effective occupation numbers in the form of $\left[ n_q(x) \right]^q$ \cite{Symmetry}. In all calculations one has always to ensure that the corresponding $q$-exponents are nonnegative and real valued, cf. \cite{Symmetry} for details.

\section{The interplay between the fugacity $z$ and the nonextensivity $q$ - dynamics vs nonextensivity}
\label{sec:Interplay}

In such a QPM the investigation of the interplay between the dynamics and nonextensivity comes to the study of the interplay between the fugacities $z_{i=q,g}$, representing the dynamics, and the nonextensivity $q$, representing the action of the environment (always keeping in mind that $q\neq 1$ combines the action of all factors causing nonextensivity, even in the absence of any interaction, i.e., when $z_{i=q,g} = 1$, and that the nonextensive particles cannot be considered as being completely free \cite{Biro-1,Biro-2}). Note first that both $z$ and $q$ deform the original Fermi-Dirac or Bose-Einstein distributions of noninteracting particles, but they do it in two different, incompatible ways. The $z$-deformation is {\it local}, it is supposed to depend only on the scaled temperature $\tau = T/T_{crit}$. As can be seen in Eq. (\ref{fex-d}) the action of $z = z\left( \tau=T/T_{crit}\right)$ is, technically, the same as action of some kind of "artificial chemical potential", $ \mu_z = T\ln z(\tau)$ (which adds to the usual chemical potential). This means that the form of the distribution remains intact, only its argument changes (therefore adding some true chemical potential $\mu$ to the original $z$-QPM can be combined with fugacity $z$ and results in some new "effective" $\tilde{z}$). In contrast to this the $q$-deformation is {\it global} and, in principle, the parameter $q$ is assumed to be independent of temperature\footnote{At least as long as there is no energy exchange between the heat bath and the environment, this would result in the replacement $T \to T=T_{effective} = T(q)$ (so far this has been shown only for the Boltzmann statistics \cite{WW_Tq,WW_Tq1}).}. This means that one cannot fully replace $z$ by $q$ (and vice versa). To see it better let us identify the effective extensive (i.e., calculated for $q=1$) occupation number considered as a function of fugacity $z$, with its nonextensive counterpart defined mainly by a nonextensivity parameter $q$ (with no other dynamical effects, i.e., with $z=1$, note that it is given not by $n_q$, but by $n^q_q$):
\begin{equation}
n(x;z,q=1) = \frac{1}{\frac{1}{z}e(x) - \xi} = \left[ n_q(x;z=1,q)\right]^q =\left[ \frac{1}{e_q(x) - \xi}\right]^q.   \label{equality}
\end{equation}
The immediate result is the relation between $z$ and $q$,
\begin{equation}
z = z(q) = \frac{ e(x) }{ \left[ e_q(x) - \xi\right]^q + \xi} \stackrel{q \to 1}{\Longrightarrow} 1.\label{z-funq}
\end{equation}
The reverse relation, $q=q(z)$, is a transcedental function, not accessible analytically. However, relation (\ref{equality}) would imply that the resulting fugacity must depend not only on the nonextensivity $q$, but also on the energy, chemical potential and temperature, i.e., on $x$, $z = z(q;E,\mu,T)$. This is unacceptable because fugacity was introduced to model the dynamics which cannot depend on $x$.

The same argument also precludes the seemingly more general formulation of the $qz$-QPM discussed in \cite{JRGW3}. One could, for example, start with a system of  particles the interactions of which are described by some chemical potentials $\mu^{(i)}$ (not connected with any fugacities, $z=1$) and which are immersed in some nonextensive environment characterized by a nonextensive parameter $q\neq 1$. Adding now to such a system some additional interaction by means of fugacities $z^{(i)}_q$, results in
\begin{eqnarray}
n_q\left[ x^{(i)} \right] &=& \frac{1}{ \frac{1}{z_q^{(i)}}e_q\left[ x^{(i)}\right] - \xi} = \frac{1}{e_q\left[ x_q^{(i)}\right] - \xi}, \label{zenq}\\
x_q^{(i)} &=& \beta \left[ E^{(i)}_q - \tilde{\xi}\mu_q^{(i)}\right]\! -\! \zeta_q^{(i)} =
x^{(i)} \left[ z_q^{(i)}\right]^{1-q}\! -\! \zeta_q^{(i)} = \ln_q\left[ \frac{e_q(x)}{z_q}\right], \label{xenq}
\end{eqnarray}
where $E_q^{(i)} = E_i\cdot\left[ z_q^{(i)}\right]^{1-q}$, $\mu_q^{(i)} = \mu^{(i)}\cdot\left[ z_q^{(i)}\right]^{1-q}$ and $\zeta_q^{(i)} = \ln_{2-q}\left[z_q^{(i)}\right]$, i.e., the energy and chemical potentials become $q$-dependent quantities. However, such induced $q$-dependence of the initial energy and chemical potential is rather unphysical and precludes further applications of this approach.

In \cite{Symmetry,JRGW3} we have discussed in detail how the already known dynamics, represented by effective fugacities $z_{i=q,g}$ introduced in \cite{QPM}, is modified by immersing the system under consideration in some nonextensive environment characterized by the parameter $q\neq 1$. Fugacities $z_{i=q,g}$ were obtained from the results of lattice QCD simulations (which were treated as input data) by comparing the pressures of the gluons and quarks (expressed as functions of the fugacities) with the corresponding pressures obtained from the lattice calculations \cite{QPM}. As a result the effective fugacities were received as functions of scaled temperature, $\tau = T/T_{cr}$, with $T_{cr}$ being the critical temperature. To describe the lattice QCD data over the whole range of $\tau$ considered (i.e., for $\tau < 4$), the $\tau$ range had to be divided in two sectors, each of which was parameterized by a different functional form with the cross-over point at $\tau_g = 1.68$ for gluons and $\tau_q=1.7$ for quarks:
\begin{equation}
z^{(q,g)}(\tau) = a_{(q,g)}\exp\left[ - \frac{b_{(q,g)}}{\tau^5}\right]\Theta\left( \tau_{(q,g)} - \tau\right) + a'_{(q,g)}\exp\left[ - \frac{b'_{(q,g)}}{\tau^2}\right]\Theta\left( \tau - \tau_{(q,g)}\right).\label{zzz}
\end{equation}
The fit parameters are: $a_{(q,g)} = (0.810,~0.803)$, $a'_{(q,g)} = (0.960,~0.978)$, $b_{(q,g)} = (1.721,~1.837)$, $b'_{(q,g)} = (0.846,~0.942)$. The resulting $z^{(q,g)}(\tau)$ are shown in Fig. \ref{Rzqg}) (black curve) where we have also shown the positions of the points obtained from lattice QCD (by arrows) \cite{QPM}. Note that $z^{(q,g)}(\tau \to \infty) = a'_{(q,g)} < 1$. This result, when treated seriously, indicates that with increasing temperature $\tau$ the system of quarks and gluons considered in the QCD lattice simulations never becomes a gas of free streaming non-interacting particles. From the nonextensive point of view, presented below, this would mean that lattice QCD keeps a memory of the interaction and describes a quark-gluon system which is intrinsically nonextensive, with $q<1$, as will be shown below.

The same pressures of the gluons and quarks obtained from the lattice QCD results must now be reproduced in the nonextensive circumstances defined by the nonextensivity parameter $q$, i.e., by using nonextensive particle occupation numbers $n_q\left[ \tilde{x}^{(i)}\right]$ given by Eq. (\ref{enq}) with nonextensive effective fugacities $z_{q\neq 1}$. As in \cite{Symmetry}, we used for this purpose the same parameterization of $z_q(\tau)$ as was used for $z(\tau)$ in Eq. (\ref{zzz}), but with the $q$-dependent values of parameters $(a,b)$ and $(a',b')$ obtained from the same lattice QCD results. The results turned out to be very sensitive to the amount of nonextensivity imposed, limiting our considerations to $|q-1| <<1$. Because the corresponding changes in fugacities, $\delta z_q = z_{q\neq 1} - z_{q=1}$, for these values of $q$ were also small, $\delta z_q <1 $, the exact formulas practically coincided with a linear in $(q-1)$ approximation \cite{Symmetry} (for clarity of presentation we suppress indices $i=q,g$):
\begin{eqnarray}
z_q &\simeq& z_{q=1} + \delta z_q = z_{q=1} \left[ 1 + (1-q) \cdot F\left(q=1,z_{q=1}\right)\right], \label{zq}\\
&&F = \frac{\int_0^{\infty} dp p^2 \left\{ \ln^2[ 1 - \xi e\left( -x;z_{q=1}\right)] + n\left(x;z_{q=1}\right) x^2\right\}}{2\int_0^{\infty} dp p^2 n\left(x;z_{q=1}\right)}. \label{F}
\end{eqnarray}
\begin{figure}[b]
\hspace{-3mm}
\resizebox{0.58\textwidth}{!}{%
  \includegraphics{zquarksgas.eps}
}
\hspace{-12mm}
\resizebox{0.58\textwidth}{!}{%
  \includegraphics{zgluonsgas.eps}
}
\caption{(Color online) Fugacities $z_{(q,g})$ (left and right panels, respectively) obtained in the $z$-QPM \cite{QPM} (with $q=1$, black curves; the lattice QCD results which were fitted are also shown by arrows \cite{QPM}) as functions of the relative temperature $\tau$. They are compared with the nonextensive fugacities $z_q$ (red and blue curves) obtained from Eqs. (\ref{zq}) and (\ref{F}) for, respectively, $q=0.985$ and $0.981$ for quarks and for $q=0.989$ and $0.985$ for gluons. They were chosen in such a way as to get $z_{(q,g)} = 1$ for, respectively, $\tau = 4$ and $\tau=5$, which are also the values of $\tau$ onwards from which we assume that $z(\tau)$ remains constant.}
\label{Rzqg}
\end{figure}
As $F > 0$, this means that we always have $z_{(q>1)} < z_{(q=1)} < z_{(q<1)}$.

We shall concentrate now on the $q<1$ case, the example of which is displayed in Fig. \ref{Rzqg}. As already mentioned before, it looks like that for $q=1$ (the extensive case) the fugacity never reaches the value $z=1$ corresponding to a noninteracting gas of quarks and gluons, at least not in the considered range of $\tau < 4$ (which is otherwise quite big), one always has an interacting (confined?) system of quasiparticles.  A closer look reveals that this is because for $z=1$ the pressure in such a gas would exceed the pressure obtained from the lattice calculations. From our experience with a nonextensive version of the Nambu-- Jona-Lasinio model of QCD matter \cite{JRGW3} we know that this pressure is reduced when the system becomes nonextensive with $q<1$. When applying it here, using Eqs. (\ref{zq}) and (\ref{F}), one indeed observes that the respective $z_{(q<1)}$ exceeds $z_{(q=1)}$ for the same values of $\tau$ in such a way as to keep up with the original lattice QCD pressure, and very soon they reach, for a given $\tau$, the value $z_{(q<1)}=1$, and it will exceed it for larger values of $\tau$. This means that our system becomes a system of non-interacting nonextensive quasiparticles (noninteracting in the sense that the only interaction is that provided not by the dynamics but by the nonextensivity $q$ \cite{Biro-1,Biro-2,Biro-3}). Decreasing $q$ further for this $\tau$ (or increasing $\tau$ while keeping the same $q<1$) would result in $z_{(q<1)}$ exceeding unity, which we consider unrealistic (because it would correspond to interactions not present in QCD). Therefore, once the $z_{(q<1)}(\tau)$ reaches unity we assume, in what follows, that it remains unity for larger values of $\tau$ (and for the same $q$, as in Fig. \ref{Rzqg})).

In what concerns the examples shown in Fig. (\ref{Rzqg}), note that the range of $\tau$ considered in the lattice QCD simulations was limited, with the last point located at $\tau= 3.7$ \cite{QPM}. This means that extrapolation of $z(\tau)$ obtained from fits to lattice QCD results, cf. Eq. ({\ref{zzz}), to $\tau >> 3.7$, is highly uncertain. Because of this, we assume that at some $\tau$ the fugacity stops increasing and remains constant thereafter. As an example, we choose for further consideration two choices: $\tau =4$ and $\tau =5$. In other words, we tacitly assume that, starting from these values of $\tau$, the QCD interactions remain essentially constant (or decrease very slowly) therefore the corresponding fugacities remain virtually the same: $z_q = z_q(\tau=4)$ or $z_q = z_q(\tau=5)$, respectively, cf. Fig. (\ref{Rzqg}) (as we will see in a moment, this will allow the determination of the non-extensiveness of our QCD system of quarks and gluons). For each of these two values of $\tau$ the respective nonextensivities  $q$ were found (separately for quarks and gluons) such that the values of the corresponding $z_q(\tau)$ reach unity (corresponding to "free, nonextensive gas"). As was the case of the nonextensive $z(\tau)$, we assume that for $\tau > 4$ (or $\tau >5$)  the corresponding values of $z_q(\tau)$ remain unity. As see in Fig. \ref{Rzqg}) the respective values of $q$ are, for quarks, $q=0.981$ for $\tau = 4$ and $q = 0.985$ for $\tau = 5$ (for gluons they are only slightly larger). We therefore argue that the nonextensivity of quarks and gluons is in the range $0.98 < q < 0.99$.

For the case of $q>1$, not shown in Fig. (\ref{Rzqg}), the corresponding $z_q(q,g)$ are always below the $q=1$ curve for the full range of $\tau$. This is because in this case the pressure is increased, therefore the strength of the dynamical interactions (given by $z_q$) has to be reduced. This is a reflection of the fact that in both cases considered the non-extensiveness comes from different sources which work towards confinement for $q<1$ but work against it for $q>1$ (for example, by introducing some extra intrinsic fluctuations). In other words: because in lattice QCD one always has $z<1$ (which reflects the fact that it describes, in fact, more or less confined systems), in a possible equivalent nonextensive approach this would correspond to a nonextensivity parameter smaller than unity, $q<1$ (which in the colloquial understanding of non-extensiveness corresponds to broadly understood correlations). The $q>1$ type of nonextensivity corresponding (again, broadly speaking) to some kind of effective, intrinsic, fluctuations manifesting themselves as fluctuations of the temperature $T$, would demand $z>1$, which is forbidden in lattice QCD. On the other hand, if taken seriously, one could probably argue that this comes not so much from the QCD dynamics, but rather from the environment (and would, figuratively speaking, come from the heterogeneity of the heat bath and from possible energy transfers to and from it \cite{WW_Tq,WW_Tq1}.

\section{Summary}
\label{sec:SC}

Generally speaking, systems for which the fugacities $z = 1$ in the limit of large $T$ are extensive systems where a single-particle description works without an additional field persisting for high temperatures (without long range interactions). If, for $T\to \infty$, the fugacity $z<1$ our system is nonextensive with reduced pressure, whereas for $z>1$ in this limit it is nonextensive with increased pressure (always in comparison with the situation for the $z=1$ case). This means that we can adjust the system pressure in models with fugacity to determine its model non-extensivity.

Because, as already mentioned, the action of nonextensivity $q$ has a {\it global} character whereas that of the fugacity $z$ is {\it local}, therefore, in principle they are complementary and are not substitutable for each other. To put it differently, usually fugacity $z$ models phenomenologically dynamics of the mean field theory type (MFT) in which there are no correlations and intrinsic fluctuations, whereas these two features are exactly those described by the nonextensivity $q$. Therefore, the results obtained here should be considered as some illustration of how to introduce effects of correlation or fluctuations into Finally,MFT models, or, vice versa, how to add to the system with correlations and/or fluctuations a certain smoothing factor for these effects.

Finally, let us mention a recent work \cite{HJ} in which the authors use numerical field theoretical simulations to calculate particle yields. It turned out that in the model of local particle creation one observes rather natural deviations from the pure exponential distributions towards Tsallis like $q$-exponents. Interestingly, for a quantum SU(3) Yang-Mills gauge theory applied to gluons (which could perhaps be a good replacement for our QCD) they obtained similar values of $q<1$ and similar dependence on the temperature as we observe in our case. The $q>1$ was obtained for a toy model of classical $\Phi^4$ theory. From our point of view  the natural explanation would be the presence of some kind of confinement in the first case (which would correspond to our lattice QCD situation) and its lack in the second case (in our case it would correspond to immersing our system in a heat bath which is so nonuniform that it overcomes the confinement forces).
\vspace{1cm}

The research of GW was supported in part by the National Science Center (NCN) under contract 2016/22/M/ST2/00176. We would like to thank warmly Dr. Nicholas Keeley for reading the manuscript.


\section{Introduction}
Since the discovery of the first  quasar at redshift $z\geq4$ \citep{1987Natur.325..131W}, the number of active galactic nuclei (AGN) known at the highest redshifts is continuously increasing. As a result of extensive observing campaigns and surveys \citep[e.g.,][]{2000AJ....120.1579Y,2010AJ....140.1868W,2011AJ....142...72E, 2016arXiv161205560C,2016MNRAS.460.1270D, 2017AJ....154...28B}, there are nearly $3000$ sources identified between $z=4$ and $z=7.54$ to date, yet only $\sim6.5$ per cent of these exhibit radio emission \citep{2017FrASS...4....9P} detected at 1.4~GHz by either the Very Large Array (VLA) Faint Images of the Radio Sky at Twenty-Centimeters\footnote{\url{http://sundog.stsci.edu/}} \citep[FIRST,][]{1995ApJ...450..559B,2015ApJ...801...26H} survey or the  National Radio Astronomy Observatory (NRAO) VLA Sky Survey\footnote{\url{https://www.cv.nrao.edu/nvss/}} \citep[NVSS,][]{1998AJ....115.1693C}. The remaining sources are either outside the footprint of both FIRST and NVSS, or are below the detection limit of the two surveys.

Several studies are available in the literature on the overall radio luminosity function (LF) for AGN  \citep[e.g.,][]{1991ASPC...18..113C,1998MNRAS.300..625W,1999ApJ...511..612G,2014MNRAS.445..955B,2016MNRAS.460....2P,2017ApJ...842...87M,2017A&A...602A...6S}. Many authors found evidence for bimodality in the distribution of the AGN population \citep[e.g.,][]{1980AnA....88L..12S,1989AJ.....98.1195K,1999ApJ...511..612G,2007ApJ...654...99W}, dividing the population to distinct samples of radio-loud and radio-quiet objects. The validity of this dichotomy was recently challenged and found to be the result of observational and mathematical bias effects \citep{2008MNRAS.387..856Z,2012ApJ...754...12M,2012ApJ...759...30B,2013ApJ...768...37C}.
The weak radio emitters with flux densities below the detection threshold ($\sim1$~mJy for FIRST and $\sim2.5$~mJy for NVSS) might also exhibit radio jets at parsec scales \citep[e.g.,][]{2017ApJ...835L..20W}, but generally remain silent and hidden in flux density limited surveys due to their low radio power and large distance. Faint objects with radio flux densities below the FIRST detection limit can either host radio-quiet AGN, or harbour radio-loud quasars that reside at higher redshifts.

The lurking low- and modest-power members of the high-redshift AGN (hAGN) population can be uncovered by going below the survey threshold and estimating their typical flux densities by means of radio image stacking \citep[e.g.,][]{2007ApJ...654...99W}. Numerous works with stacking analysis are available in the literature on various samples of radio-weak quasars. Analysis of stacked FIRST images at $\sim 8000$ radio-quiet quasar positions from the 2dF QSO redshift survey at medium redshifts ($z\lesssim2.3$) resulted in median flux density levels between $20$ and $40~\mu$Jy  \citep{2005MNRAS.360..453W}. \citet{2008AJ....136.1097H} also found 10s of $\mu$Jy flux densities for a mixed set of Sloan Digital Sky Survey (SDSS) and Luminous Red Galaxy samples. Stacking of extremely red quasars at $2\leq z\leq 4$ using VLA maps at $1.4$ and $6.2$~GHz yielded similar values \citep{2018MNRAS.477..830H}.

There is a debate on the nature of physical processes responsible for the radio emission at the faint end ($P_\mathrm{1.4GHz}\leq10^{22.5-23}$~W~Hz$^{-1}$) of the continuous radio LF, with possible contribution of AGN jets \citep{1997ARAnA..35..607Z,2018ApJ...869..117R}, AGN winds  \citep{2010ApJ...711..125J,2014MNRAS.442..784Z,2016MNRAS.455.4191Z,2018ApJ...869..117R}, and enhanced star formation or starburst events in the host galaxy \citep{2004MNRAS.352..399J,2013ApJ...768...37C,2013ApJ...778...94R,2018MNRAS.477..830H}.  The origin of radio emission in radio-quiet AGN was also hypothesized to be driven by magnetically heated accreation disk coronae \citep{2008MNRAS.390..847L,2019MNRAS.482.5513L}.
Upper limits on star formation rate (SFR) determined in stacking studies vary between a few to $10$~M$_\odot$~yr$^{-1}$ \citep[e.g.,][]{2007AJ....134..457D,2008AJ....136.1097H}. In other works, low-luminosity AGN jets and  circumnuclear star-forming regions in the host galaxies together are considered responsible for the radio emission \citep[e.g.,][]{2011ApJ...730...61K,2011ApJ...742...45P}.
It was also found that the radio LF usually peaks at higher redshifts for radio-loud sources than for radio-quiet AGN \citep{1999ApJ...511..612G,2017ApJ...842...87M}.
A multi-wavelength study of  a radio-selected Cosmic Evolution Survey (COSMOS) field sample ($S_\mathrm{1.4GHz}\geq37~\mu$Jy) at $z\leq6$ suggested star formation as the dominant process for the sub-mJy radio population, explaining the dichotomy with AGN activity modes and connection to AGN--host galaxy feedback \citep{2017AnA...602A...3D}. Another study using a sample with $S_\mathrm{1.4GHz}\geq11.5~\mu$Jy applied spectral energy distribution fitting and led to similar conclusions: the peak of the radio LF is at $z\sim2$ \citep{2018AnA...620A.192C}. The radio loudness dichotomy was addressed also at low radio frequencies \citep[$120-168$~MHz,][]{2019AnA...622A..11G}, concluding that low-power quasars are dominated by star formation.

In this paper, we utilised $2229$ empty-field FIRST radio maps centred at positions of hAGN (quasars) identified in the optical or near-infrared. We analysed them with mean and median stacking methods to uncover the underlying radio quasar population, and to examine the possibility of the radio emission originating from either AGN activity or star formation. The aspects of sample selection are detailed in Section~\ref{sec:sample}. We describe the stacking procedure in Section~\ref{sec:stacking}. Results and derived properties are discussed in Section~\ref{sec:discussion}. We summarise our findings and conclude the paper in Section~\ref{sec:summary}.

Throughout this work, we assumed a standard $\Lambda$CDM cosmological model with $\Omega_\mathrm{m}=0.3, \Omega_{\Lambda}=0.7$, and $H_0=70$~km~s$^{-1}$~Mpc$^{-1}$ for calculations.

\section{Sample selection and data}\label{sec:sample}
The sample for the stacking analysis was defined by using the high-redshift ($z \geq 4$) AGN catalogue by \citet{2017FrASS...4....9P}{\footnote{The regularly updated version of the catalogue can be accessed at \url{http://astro.elte.hu/~perger/catalog.html}}}. We selected optically identified AGN from the catalogue with positions falling into the FIRST survey footprint. The flux density of the objects selected is below the detection limit of the survey (typically $\sim1$~mJy). At the time of the analysis, the total number of such AGN was $2232$. Intensity  fluctuations (rms noise) at the optical AGN positions in the FIRST images vary in the range from tens of $\mu$Jy to $\lesssim 1$~mJy. To test an underlying high-redshift radio AGN population, we stacked  FIRST radio maps centred at these positions. The images in Flexible Image Transport System (FITS) format \citep{1981AnAS...44..363W} were obtained from the FIRST image cutout service\footnote{\url{https://third.ucllnl.org/cgi-bin/firstcutout}}. We downloaded images of $4\farcm5\times4\farcm5$ size for each individual AGN position. After excluding the largely incomplete or entirely ragged maps, we obtained $2229$ images for the stacking process. The redshift distribution of the stacked sample is shown in Fig.~\ref{fig:histogram}, and the first 5 entries of the list of $2229$ AGN are given in Table~\ref{tab:sample}. The full list is available in the supplementary material.

\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{histogram.eps}
    \caption{Redshift distribution of the $2229$ high-redshift active galactic nuclei used for the stacking analysis.}
    \label{fig:histogram}
\end{figure}

\begin{table*}
    \centering\caption{The first 5 entries of the AGN sample used for the stacking analysis.}
    \begin{tabular}{lllll}
         \hline\hline
        Name                        &   R.A.  (h~m~s)                 &   Dec  ($\degr~'~''$)                   &   Redshift & Reference\\
    \hline
SDSS~J000046.69+010951.2	&	$00~00~46.69$	&	$+01~09~51.24$	&	4.25	&	\citet{2012AnA...548A..66P}	\\
SDSS~J000124.23+111212.6	&	$00~01~24.23$	&	$+11~12~12.69$	&	4.30	&	\citet{2017AnA...597A..79P}	\\
SDSS~J000404.71+000039.0	&	$00~04~04.71$	&	$+00~00~39.08$	&	4.31	&	\citet{2012AnA...548A..66P}	\\
SDSS~J000457.11$-$000538.7	&	$00~04~57.11$	&	$-00~05~38.78$	&	4.05	&	\citet{2012AnA...548A..66P}	\\
SDSS~J000527.14+025813.2	&	$00~05~27.15$	&	$+02~58~13.29$	&	4.11	&	\citet{2012AnA...548A..66P}	\\
    \hline
    \end{tabular}\\
   {\it Notes.} Column 1 -- object name, Column 2 -- right ascension, Column 3 -- declination, Column 4 -- spectroscopic redshift, Column 5 -- literature reference of discovery. The full table of $2229$ AGN is available in the electronic version of the journal.
    \label{tab:sample}
\end{table*}

\section{Stacking analysis}\label{sec:stacking}

\subsection{Stacking of catalogue sources}

The sample of $2229$ objects were binned based on the redshift of each source into four subsamples with approximately equal number of images in each bin. There are $554$, $559$, $559$, and $557$ objects in the four bins with redshift boundaries of $4.0\leq z_1<4.1,4.1\leq z_2<4.3,4.3\leq z_3<4.7,$ and $4.7\leq z_4<7.6$, respectively.

We performed the stacking in both mean and median procedures. Both the mean  \citep[e.g.][]{2007ApJ...654...99W,2007AJ....134..457D,2008AJ....136.1097H} and median \citep[e.g.][]{2007ApJ...654...99W,2007AJ....134..457D,2008AJ....136.1097H,2018MNRAS.477..830H,	2018ApJ...869..117R} stacking methods are commonly used in the literature. Mean stacking determines the arithmetic average of intensities in the image pixels within the sample. However, this method is very sensitive to outlying values in the sample, as well as to the threshold set to avoid the contaminating point sources. In turn, median stacking deals with the median of intensity values corresponding to the image pixels in the sample.

To evaluate the applicability of maps resulted from both methods, we calculated root-mean-square (rms) noise, maximum of the intensity, and signal-to-noise ratio (SNR) for both the mean and median stacked images in each bin and for the full sample.
Typical rms noise values in the separate redshift bins span the range $10-20~\mu$Jy~beam$^{-1}$ for mean stacked images, while median stacking resulted in $\sim 7\mu$Jy~beam$^{-1}$ noise levels. The full-sample maps have $8~\mu$Jy~beam$^{-1}$ and $3~\mu$Jy~beam$^{-1}$ rms noise for the mean and median methods, respectively.

\begin{table}
\centering
\caption{Properties of median-stacked maps.}\label{tab:data}
\begin{tabular}{cccc}
\hline\hline
Bin & rms [$\mu$Jy~beam$^{-1}$] & $I_\mathrm{max} $ [$\mu$Jy~beam$^{-1}$] & SNR \\
\hline
1	&		7		&		27		& 	4		\\
2	&		7		&		52		& 	7		\\
3	&		7		&		38		& 	5		\\
4	&		7		&		30		& 	4		\\
All	&		3		&		35		& 	11		\\

\hline
\end{tabular}\\
{\it Notes.} Column 1 -- redshift bin, Column 2 -- image noise, Column 3 -- maximum intensity, Column 4 -- signal-to-noise ratio
\end{table}

Despite the high SNR values (up to $50$) found in each mean-stacked map, the method provides insignificant results in the search of a hidden central source. The maps are contaminated and dominated by peaks at random off-centre locations, caused by strong intensity peaks from the stacked individual maps. The presence of off-centre sources leads to an increase of the rms image noise in the field compared to median stacking. Therefore mean-stacked maps are not considered in the further study and calculations.

Median maps on the contrary are not sensitive to occasional bright off-centre sources. They show a nearly uniform noise level in each subsample ($7~\mu$Jy~beam$^{-1}$), revealing a protruding radio peak at the image centre for the entire sample with SNR exceeding 10. Intensity maxima for the median-stacked images are found to be $27, 52, 38$ and $30~\mu$Jy~beam$^{-1}$ for the redshift-binned data, and $35~\mu$Jy~beam$^{-1}$ for the full sample. Radio maps obtained by median stacking are shown in Fig.~\ref{fig:bins}, for which the calculated image properties are listed in Table~\ref{tab:data}, while the radial profile of the full-sample median-stacked image is illustrated in Fig.~\ref{fig:radial}.

The image rms noise reached by stacking is expected to decrease with respect to the single-image rms by $\sqrt N$, where $N$ is the number of stacked images. The original FIRST image noise levels are $\approx 0.15$~mJy. For the sample of $N=2229$ objects, this predicts the rms of $150~\mu$Jy~beam$^{-1}$ divided by $\sqrt{2229}$, i.e. $\approx 3~\mu$Jy~beam$^{-1}$ to be reached with stacking. Indeed, this is equal to the actual rms values obtained for the full-sample median-stacked images.

\begin{figure*}

\centering
    \subfloat[Median-stacked maps for the four redshift bins.]{\includegraphics[width=0.85\linewidth]{medians_w.eps}}\\
    \subfloat[Median-stacked map derived from all images in the sample. The inner $30\arcsec \times 30\arcsec$ area is shown.\label{fig:bin_all}]{\includegraphics[width=0.85\linewidth]{allmedian_w.eps}}
\caption{The central $2\arcmin \times 2\arcmin$ area of the median-stacked FIRST image cutouts of $z\geq 4$ optically identified radio-non-detected AGN. For the binned data, the inner $20''\times20''$ area is also shown in the insets. Maximum intensities in the central 5 arcsec radius circle areas are indicated. The first contours are drawn at $\pm3$ times the rms noise level. Additional contours are at $6$ and $9$ times the rms. }\label{fig:bins
\end{figure*}

\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{radial_all_in_unbiased.eps}
    \caption{Pixel-by-pixel radial plot of the full-sample median-stacked intensity values for the inner $15''$ radius area. 1 pixel equals $1.8''$. For the data points, the unbiased values are shown, after applying the bias correction factor $1.4$ \citep{2007ApJ...654...99W}.}
    \label{fig:radial}
\end{figure}

\subsection{Stacking of fake sources and the image noise}
To check the reliability of our results, we created four samples of the same size, with arbitrary source coordinates in the manner of adding or subtracting $1\degr$ to/from either or both the right ascensions and declinations of objects in the  the original sample. Binning was based on the redshifts of AGN in the original (true) positions.

We did not find any radio emission in either of the four samples of fake sources, neither with mean, nor with median stacking methods. The mean images show outlier sources outside the central pixels, similarly to the results obtained for the true sample. The rms noise levels of the binned subsample maps are $10-20~\mu$Jy~beam$^{-1}$ and $\sim7~\mu$Jy~beam$^{-1}$ for the mean and median stacking, respectively. Stacking of all images resulted in rms values comparable to that of the real sample and consistent with the expected values from the stacking procedure, $8~\mu$Jy~beam$^{-1}$ for mean and $3~\mu$Jy~beam$^{-1}$ for median maps. The SNR values for the mean images vary from $16$ to $50$ in the separate redshift bins, and are below $\approx 30$ when the full fake samples were included in mean stacking. The image peaks are caused by occasional bright off-centre sources appearing in some fields, similarly to the real sample. For median stacking, the typical values in the fabricated samples and in their subsets are in the range of $4\lesssim \mathrm{SNR}\lesssim5$, and indicate no significant detection at the central location, as expected. Due to the similar SNR values in the samples of fake sources, the binned data were not considered for further analysis, only the full sample.

\subsection{Model fit to the stacked data}

For further analysis, we fitted a circular Gaussian model component to the central region of the full-sample median-stacked image, to characterise the brightness distribution using the {\sc imfit} task in the U.S. National Radio Astronomy Observatory (NRAO) Astronomical Image Processing System\footnote{http://www.aips.nrao.edu/index.shtml} \citep[{\sc aips},][]{2003ASSL..285..109G}. This resulted in an unresolved point source. By multiplying the fitted model component flux density with the bias correction factor of $1.4$ suggested by \citet{2007ApJ...654...99W}, we obtained an unbiased value of $S_\mathrm{1.4GHz}=52\pm1\,\mu{\mathrm Jy}$.

\subsection{Characteristic 1.4-GHz radio power }

To derive a characteristic value of the monochromatic rest-frame 1.4-GHz radio power for the hAGN population, we co-added the individual FIRST image cutouts around the positions of the full sample of $2229$ objects that are individually not detected in the FIRST survey. Motivated by the model-fitting results of the median stacking, we assumed that the radio emission originates from point sources unresolved in FIRST. Thus the flux density was derived from the cumulative brightness value of the central pixel for the following analysis. We consider that the co-added flux density of the $2229$ hAGN is the sum of the flux densities of the individual sources, and that all hAGN included in the stacking process have the same radio power at $1.4$~GHz, $P_\mathrm{1.4GHz}$. This latter assumption is obviously not true, but the value derived this way is characteristic to the hAGN sample as a whole.

We used the following relationship between flux density and power (including $K$-correction):
\begin{equation}
    S_\mathrm{1.4GHz,sum}=\sum_{i=1}^{N}\frac{P_\mathrm{1.4GHz}}{4\pi D_{\mathrm{L}i}^2(1+z_i)^{-\alpha-1}},
\end{equation}
where $S_\mathrm{1.4GHz,sum}=77$~mJy is the cumulative flux density derived from co-adding the $N=2229$ FIRST map centres, $\alpha$ is the radio spectral index (using the convention $S_\nu\propto\nu^\alpha$, where $\nu$ is the frequency), and $D_{\mathrm{L}i}$ and $z_i$ are the luminosity distance and the redshift of the $i$-th individual source, respectively. Since the individual power-law spectral indices are unknown for the sources, we assumed a common value for the whole sample. Various radio spectral indices appropriate for galaxies are found in the literature \citep[e.g.][]{2003ApJ...599..971H,2007AnA...462..525H,2007ApJ...654...99W,2018MNRAS.477..830H}. Therefore we calculated the characteristic power considering several different values in the range $-0.5\leq\alpha\leq-1$, corresponding to steep radio spectra. The estimated values of the $1.4$-GHz power for each spectral index are listed in Col.~2 in Table~\ref{tab:powsfr}.

\subsection{Upper limit on star formation}

Radio emission from galaxies may originate from AGN activity in the nuclear region where synchrotron radiation is produced in powerful relativistic plasma jets driven by accretion onto the central supermassive black hole. On the other hand, synchrotron emission from relativistic electrons and free-free emission from ionized hydrogen regions may be widespread in the AGN host galaxy as a consequence of recent start formation \cite[e.g.][]{1992ARA&A..30..575C}.

We can assume that the radio emission in our stacked hAGN sample originates solely from star formation taking place around the central regions of the host galaxy, rather than from AGN activity. This way we can calculate an upper limit on the star formation rate (SFR) using the estimated characteristic $1.4$-GHz power value, applying various correlations between the radio power and SFR from the literature \citep{2003ApJ...586..794B,2006ApJ...643..173S,2017MNRAS.466.4917D,2019MNRAS.482..560M}.  The SFR values obtained are listed in Cols.~3--6 in Table~\ref{tab:powsfr}. The upper limits on the SFR span about an order of magnitude, depending on the radio spectral index assumed and the correlation used. But we can conclude that the typical SFR upper limits are in the order of $10^3$~M$_\odot$~yr$^{-1}$.

\begin{table}
    \centering
    \caption{Estimated characteristic 1.4-GHz radio powers and star formation rates for the radio sources associated with the central pixels of the co-added $2229$ FIRST maps.}
    \begin{tabular}{cccccc}
    \hline\hline
        \multirow{2}{*}{$\alpha$}  &\multirow{2}{*}{$P_\mathrm{1.4GHz}$ [$10^{24}$~W~Hz$^{-1}$]}	&\multicolumn{4}{c}{SFR [M$_\odot$~yr$^{-1}$]}	 \\
            	&	        	&	(1)	    &	(2)	    &	(3)	    &	(4) 	\\ \hline
        $-0.5$  &  $2.9$   &  $1600$   & $1800$   & $370$   &  $1700$  \\
        $-0.6$  &  $3.4$   &  $1900$   & $2100$   & $420$   &  $2000$  \\
        $-0.7$  &  $4.1$   &  $2300$   & $2500$   & $480$   &  $2400$  \\
        $-0.8$  &  $4.8$   &  $2700$   & $3000$   & $540$   &  $2800$  \\
        $-0.9$  &  $5.7$   &  $3200$   & $3500$   & $610$   &  $3400$  \\
        $-1 $   &  $6.8$   &  $3700$   & $4200$   & $700$   &  $4000$ \\

    \hline
    \end{tabular}\label{tab:powsfr}\\
    {\it Notes.} Column 1 -- assumed spectral index, Column 2 -- 1.4-GHz radio power, Column 3--6 -- star formation rates calculated using various relationships by (1)~\citet{2003ApJ...586..794B}, (2)~\citet{2006ApJ...643..173S}, (3)~\citet{2017MNRAS.466.4917D}, and (4)~\citet{2019MNRAS.482..560M}
\end{table}

\section{Discussion}\label{sec:discussion}

\subsection{Flux densities}

The flux density of $52~\mu$Jy found in our median stacking analysis is comparable but somewhat higher than the values obtained in previous stacking works \citep[e.g.][]{2005MNRAS.360..453W,2007ApJ...654...99W,2008AJ....136.1097H,2018MNRAS.477..830H}. Considering that our sources are located at higher redshifts, this suggests that the hAGN are generally intrinsically more powerful than the sources studied in those samples.

For comparison, the ultraluminous quasar SDSS~J010013.02+280225.8 at $z=6.3$ --  not included in our stacking analysis due to its position outside the FIRST survey footprint -- was detected with the VLA at a higher frequency, 3~GHz ($0.65''$ resolution and $\sim3~\mu$Jy~beam$^{-1}$ sensitivity), with a flux density of $\sim100~\mu$Jy \citep{2016ApJ...830...53W}. Assuming a spectral index of $-0.7$, this scales to $\sim170~\mu$Jy at $1.4$~GHz. Further high-resolution observations performed with the Very Long Baseline Array (VLBA) at $1.5$~GHz revealed a radio structure partially resolved at $\sim 10$ milli-arcsecond (mas) level, with $\sim90~\mu$Jy flux density \citep{2017ApJ...835L..20W}. The detection of a mas-scale radio structure in this individual source raises the possibility that the FIRST-undetected sources in our stacked sample could also be targeted with sensitive very long baseline interferometry (VLBI) observations in the future, in a hope of revealing weak compact AGN-related radio emission. It also suggests that the radio emission in at least some of the sources in our $2229$-element sample does have synchrotron AGN jet contribution and is not associated with star formation only.

\subsection{Radio power and AGN radio emission}

All estimated values for the $1.4$-GHz radio power are $P_\mathrm{1.4GHz}\sim10^{24}$~W~Hz$^{-1}$ (Table~\ref{tab:powsfr}) in this study. If the radio emission in these objects would originate solely from AGN activity, the FIRST-undetected hAGN would be among the radio AGN population ($P_\mathrm{1.4GHz}\geq2\times10^{22.5}$~W~Hz$^{-1}$), since powers exceed the threshold found in empirical studies between supernova-related and AGN activity  \citep[$P_\mathrm{1.4GHz}=2\times10^{21}$~W~Hz$^{-1}$,][]{2000ApJ...530..704K, 2011AnA...526A..74M}.
The estimated radio powers also surpass values measured for prominent starburst galaxies, e.g. Arp\,220, Arp\,229A, and Mrk\,273 \citep[$P_\mathrm{1.4GHz}=2-4.5\times10^{22}$~W~Hz$^{-1}$,][]{2012MNRAS.423.1325A}.

Monochromatic powers found in this work are 2 to 5 orders of magnitude higher than those found for low-power sources at low redshifts ($z<0.3$), e.g. Lyman-break analogues (LBAs) including an LBA with dominant star formation (J0150+1308), an AGN (J1029+4829) and an AGN--SFR composite source (J0921+4509) \citep{2012MNRAS.423.1325A}. Also, low-redshift AGN SDSS\,J1155+1507,  SDSS\,J2104$-$0009, SDSS\,J2304$-$0933 \citep{2016ApJ...826..106G}, and  NGC\,3147 \citep{2004ApJ...603...42A} exhibit radio powers in the range of $5\times10^{21}$~W~Hz$^{-1}$ to a few times $10^{22}$~W~Hz$^{-1}$, or even lower powers between $10^{19}$~W~Hz$^{-1}$ to a few times $10^{20}$~W~Hz$^{-1}$ in the AGN Henize\,2-10 \citep{2012ApJ...750L..24R}, NGC\,4203, NGC\,4535 \citep{2002ApJ...581..925U,2004ApJ...603...42A}, NGC\,864, and NGC\,4123 \citep{2002ApJ...581..925U}. It was discussed in e.g. \citet[]{2005MNRAS.362...25B,2008ASPC..399..413S,2013MNRAS.433.2647S} and \citet{2016MNRAS.455.2731R} that the most luminous AGN have radio powers in the range from $10^{24}$~W~Hz$^{-1}$ to a few times $10^{26}$~W~Hz$^{-1}$, which is comparable with our results. These studies of low-redshift AGN found that host galaxies with $1.4$-GHz powers above $10^{23}$~W~Hz$^{-1}$ usually host radio-loud AGN in their centres. This outlines that AGN activity represents a significant contribution to the 1.4~GHz radio flux densities in the stacked sample.

Deep radio surveys at high redshift, such as the Extended Chandra Deep Field-South \citep[E-CDFS,][]{2008ApJS..179..114M,2013ApJS..205...13M} survey and the VLA-COSMOS 3 GHz Large Project \citep{2017A&A...602A...1S} can provide additional information about the most distant galaxy centres.
Our estimated radio powers ($\sim 10^{24}$~WHz$^{-1}$) are supported by the results achieved by e.g. \citet{2017AnA...602A...3D}, \citet{2017A&A...602A...2S}, and \citet{2017A&A...602A...6S}. In a multiwavelength analysis of $z\leq6$ radio sources in the COSMOS field at 3~GHz, \citet{2017AnA...602A...3D} divided the sample into three populations: radio-quiet (RQ) and radio-loud (RL) AGN, and star-forming galaxies (SFGs). The majority of RQ AGN hosts show enhanced star formation. Similar results were obtained by \citet{2017A&A...602A...2S}, revealing that radio sources at $1.4$~GHz with flux densities above $\sim200~\mu$Jy are predominantly AGN, and that with decreasing flux densities the SFGs take over as the dominant population in up to $\sim60$~per cent of the sample.
However, the $1.4$~GHz luminosity functions determined for AGN with radio excess (with respect to the expected contribution of SF to the radio emission) using different procedures by \citet{2017A&A...602A...6S} showed that rest-frame radio powers span the range of $\sim 10^{24}-10^{26}~$WHz$^{-1}$ for AGN at the highest redshifts, when relying solely on observational data, and
$10^{22}-10^{27}~$WHz$^{-1}$ with the application of evolution models. The fractional distribution found in the COSMOS field studies would imply an RQ--SFG dominance in the sample. Our estimated radio powers suggest that the objects in the stacking analysis are mostly radio-loud AGN.
Another study of radio sources below $\sim100~\mu$Jy based on the E-CDFS survey concluded that besides the dominance of SFGs and the decrase of the RL AGN towards the lowest flux densities, the number of RQ AGN increases \citep{2013MNRAS.436.3759B}. They also determined radio power distributions for all three populations, and found median values of both RQ and RL AGN which coincide with our estimated radio powers, while the power at the peak of the distribution for SFGs is an order of magnitude lower.

\subsection{Star formation rate}

SFR upper limits derived in this paper (Table~\ref{tab:powsfr}) are two to three orders of magnitude higher than SFR values found for individual AGN host galaxies in the literature, e.g.  $0.5-2$~M$_\odot$~yr$^{-1}$ for SDSS\,J2104$-$0009 and  SDSS\,J2304$-$0933 \citep{2016ApJ...826..106G}, $5-8$~M$_\odot$~yr$^{-1}$ for J0150+1308, J0921+4509, and J1029+4829  \citep{2012MNRAS.423.1325A}, and $\sim1$~M$_\odot$~yr$^{-1}$ with H$_\alpha$ star formation rates predicted by \citet{2002ApJ...581..925U}, as well as SFR with upper limits found in previous works using image stacking \citep[$\leq10$~M$_\odot$~yr$^{-1}$, e.g.][]{2007AJ....134..457D,2008AJ....136.1097H}. Even ultraluminous infrared galaxies were found to show star formation with rates varying between a couple of tens to a few hundreds of M$_\odot$~yr$^{-1}$  \citep[e.g.][]{2010ApJ...715..572H}, suggesting that our upper limits appreciably overestimate the contribution of star formation to the 1.4-GHz radio power.

However, turning to high-redshift objects, CO and [C II] line detection and dust emission in individual $z\sim6$ quasars indicate significant star formation in the central few kpc region in some host galaxies, with SFR in the order of $\sim1000$~M$_\odot$~yr$^{-1}$ \citep[e.g.][]{2003AnA...409L..47B,2013ApJ...770...13W,2013ApJ...773...44W,2019ApJ...876...99S}. Note that those AGN with an estimated SFR of a few thousand M$_\odot$~yr$^{-1}$ have not been included in our stacked sample because they are outside the FIRST survey coverage. For some individual objects that are included in our FIRST image stacking, molecular and atomic line observations determined a wide range of SFR values.
These start from a few tens of M$_\odot$~yr$^{-1}$, e.g. $48$~M$_\odot$~yr$^{-1}$ for CFHQS\,J0210$-$0456, $<40$~M$_\odot$~yr$^{-1}$ for CFHQS\,J2329$-$0301 \citep{2013ApJ...770...13W}, and $\sim80$~M$_\odot$~yr$^{-1}$ for CFHQS\,J0055+0146 \citep{2015ApJ...801..123W}.  SFR values of a couple of hundreds M$_\odot$~yr$^{-1}$ were found for VIKING\,J2348$-$3054 and VIKING\,J0109$-$3047 \citep[$\sim700$ and $\sim900$~M$_\odot$~yr$^{-1}$, respectively,][]{2016ApJ...816...37V},  ULAS\,J1120+0641  \citep[$\sim200$~M$_\odot$~yr$^{-1}$,][]{2017ApJ...837..146V}, and SDSS\,J0100+2800 \citep[$\sim650$~M$_\odot$~yr$^{-1}$,][]{2016ApJ...830...53W}. There are also examples of sources with SFR up to thousands of M$_\odot$~yr$^{-1}$ \citep[][]{2004ApJ...615L..17W,2019ApJ...874L..30V}.

The difference between the SFR upper limits we estimated as characteristic for the stacked hAGN sample ($\lesssim4000$~M$_\odot$~yr$^{-1}$, Table~\ref{tab:powsfr}) and the values measured for individual sources using independent methods is naturally explained if we relax the assumption that the entire radio emission comes from star forming activity in the host galaxies, without AGN contribution. Moreover, we calculated with a single characteristic radio power for all sources. In reality, the physical conditions in the individual objects are clearly more complex than just assuming a constant power, and one or the other process solely responsible for the radio emission. In fact, radio emission originating from both AGN and star formation must be present in most objects, in various proportions.

In a hope to refine the results, one may take the AGN radio luminosity function (LF) into account, up to the highest redshifts \citep[e.g.][]{2017A&A...602A...6S, 2018AnA...620A.192C}. Improvements in determining the radio LF at the earliest cosmological epochs are expected from the currently ongoing $2-4$~GHz VLA Sky Survey \citep[VLASS,][]{2015fers.confE...6M}. It will have higher sensitivity ($\sim70~\mu$Jy for the combined 3-epoch observations) and angular resolution ($2\farcs5$) compared to the FIRST survey. Also, repeating a stacking analysis similar to the one reported here, the deeper VLASS radio maps would allow for determining redshift-dependent properties of the FIRST-undetected AGN. Our study with binned sub-samples did not provide conclusive results because of the insufficient SNR in the stacked images.

\subsection{Origin of the sub-mJy radio emission}
Ascertaining the role of SF and AGN contribution to the  sub-mJy radio emission could be aided by the application of methods independent from the FIRST observations.
Based on mid-infrared (MIR) polycyclic aromatic hydrocarbon (PAH) detections in low-redshift galaxies, PAH emission lines can be found very close to the galaxy core, at $\sim1-20$~kpc distances \citep[e.g.,][]{2007ApJ...669..841S,2019ApJ...871..190M}. Since star formation indicators are found this close to the core, and our stacking analysis showed that the fitted Gaussian model component is unresolved in FIRST (corresponding to $26-35$~kpc linear size depending on the redshift), the angular resolution provided by FIRST is insufficient to distinguish between star formation and AGN related emission.
Observations with the upcoming Square Kilometre Array (SKA), in cooperation with globally distributed VLBI arrays could provide an adequate mas or sub-mas resolution \citep{2012PASA...29...42G,2015aska.confE.143P} and the thermal sensitivity of a few $\mu$Jy \citep{2015aska.confE.143P} sufficient for direct detection of compact AGN-related radio emission in the most powerful members of the radio-hAGN population.

The importance and contribution of SF to the radio emission can also be tackled using MIR observations. Considering the correlation between the MIR and $1.4$-GHz radio powers, and applying the fit parameters in \citet{2005ApJ...632L..79W}, we calculated the $8~\mu$m and $24~\mu$m powers that correspond to the characteristic $P_\mathrm{1.4GHz}$ values found in our study. Consistency of the MIR--radio correlation is considered invariant for 5 orders of magnitude of flux densities up to $z\sim3.5$  \citep{2001ApJ...554..803Y,2008MNRAS.386..953I} and was found reliable even at $\mu$Jy levels \citep{2009MNRAS.394..105G}. Characteristic MIR flux densities were estimated for our hAGN sample, based on the characteristic power derived from the co-added 1.4-GHz FIRST radio maps. Values of $10-15$~mJy  and $30-50$~mJy were found for the $8~\mu$m and $24~\mu$m flux densities, with mean values $14~$mJy and $45~$mJy, respectively. To estimate the level of AGN contamination contributing to the calculated values of SFR, we determined $q_{24}$ values using FIRST upper limits and  $24~\mu$m flux densities derived using $22~\mu$m emission measured by the Wide-field Infrared Survey Explorer \citep[WISE,][]{2010AJ....140.1868W} from the AllWISE catalog\footnote{http://wise2.ipac.caltech.edu/docs/release/allwise/} \citep{2014yCat.2328....0C}. We followed the analysis described by e.g. \citet{2013MNRAS.436.3759B}. On the one hand, all $124$ sources detected by WISE are above the theoretical limit separating RL AGN from RQ AGN and SFGs, so no further constraints could be obtained. On the other hand, the $24~\mu$m flux densities are in the range of $\sim2-11~$mJy with a mean value of $4~$mJy, which is an order of magnitude lower than the values calculated from the MIR--radio correlation, assuming all radio emission is SF-related. This indicates radio excess for the stacked objects (at least those with WISE detection), implying that the radio emission is AGN-related.

\section{Summary}\label{sec:summary}

We applied median stacking on 1.4-GHz VLA FIRST survey image cutouts centred on $2229$ optically identified but individually radio-undetected quasar positions. These objects populate the redshift range $4\leq z < 7.6$. Stacking of the full sample resulted in an unresolved point source with a bias-corrected flux density $52$~$\mu$Jy. Co-adding the radio map central pixels revealed a moderately radio-loud AGN population, with a characteristic 1.4-GHz radio power $P_\mathrm{1.4GHz}\sim10^{24}$~W~Hz$^{-1}$. Under the simplifying assumption that the entire radio emission in the sample is produced by star forming activity in the quasar host galaxies, we obtained upper limits of the star formation rate in the order of a few $1000$~M$_\odot$~yr$^{-1}$. Based on literature data on individual AGN, we argue that the source of radio emission in the sample is rather a mixture of star formation and AGN-related activity. The spatial resolution of FIRST images is not sufficient to distinguish between the different mechanisms responsible for the radio emission. Future measurements with VLBI and SKA-VLBI could help determining the relative importance of the two emission types in individual objects. Stacking studies similar to the one presented here will benefit from the improved sensitivity and angular resolution of the ongoing VLASS whose radio images could provide data sufficient for determining redshift-dependent properties of high-redshift radio quasars.

\section*{Acknowledgements}
We are grateful for the constructive comments by the anonymous referee that led to an improved discussion of our results.
This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration.
K\'{E}G was supported by the J\'{a}nos Bolyai Research Scholarship of the Hungarian Academy of Sciences, and by the Ministry of Human Capacities within the framework of the \'UNKP (New National Excellence Program).