# Efficient Optimisation of Physical Reservoir Computers using only a Delayed
Input

Enrico Picco1  Lina Jaurigue2  Kathy Lüdge2 and Serge Massar1  1 Laboratoire
d’Information Quantique CP224, Université libre de Bruxelles (ULB), Av. F. D.
Roosevelt 50, 1050 Bruxelles, Belgium
2 Institute of Physics, Technische Universität Ilmenau, Weimarer Straße 25,
98693 Ilmenau, Germany<EMAIL_ADDRESS>
††journal: opticajournal††articletype: Research Article

We present an experimental validation of a recently proposed optimization
technique for reservoir computing, using an optoelectronic setup. Reservoir
computing is a robust framework for signal processing applications, and the
development of efficient optimization approaches remains a key challenge. The
technique we address leverages solely a delayed version of the input signal to
identify the optimal operational region of the reservoir, simplifying the
traditionally time-consuming task of hyperparameter tuning. We verify the
effectiveness of this approach on different benchmark tasks and reservoir
operating conditions.

##  1 Introduction

In recent years, the quest for efficient and powerful machine learning
paradigms has led researchers to explore unconventional computational
architectures inspired by the dynamics of physical systems. Reservoir
Computing (RC) [1] is a possible promising approach. RC diverges from
conventional neural networks by introducing a distinct architecture
characterized by a fixed, randomly initialized recurrent layer, called
"reservoir", coupled with a simple linear readout layer: its characteristic
structure reduces the training computational complexity, and consequently, the
energy consumption. RC possesses a remarkable capacity to process temporal
data, distinguishing itself across a large variety of tasks, such as
equalization of distorted nonlinear communication channels [1, 2], audio
processing [3, 4], weather prediction [5], image [6, 7] and video [8, 9]
classification. RC has also gathered attention for its feasibility to being
implemented on a wide range of physical substrates [10, 11]. In particular,
photonic RC stands out as one of the most promising hardware platforms, thanks
to its advantages in parallelization [12, 13], high speed [14] and minimal
hardware requirements [15].

Hyperparameter optimisation is a crucial aspect of tuning reservoir computers,
aiming to enhance their performance across various tasks. Various optimization
techniques are employed to search for the most effective hyperparameter
values. Grid search, random search, and more sophisticated methods like
Bayesian optimization [16, 17] and genetic algorithms [18, 19] have been
proposed. Nonetheless, the optimization of hyperparameters remains a critical
challenge to the widespread adoption of RC [20]. The challenge lies in
exploring the vast hyperparameter space to find configurations that yield
optimal results for a given task, and the sensitivity of the reservoir to
hyperparameters. The high-dimensional nature of the hyperparameter space can
make the search process computationally expensive and time-consuming, even
more so in the case of slow experimental reservoirs that require more
hyperparameters and for which one iteration of the grid search takes up to
several hours [8].

In this work we study a novel approach, first introduced in [21] in numerical
simulations, based on the use of a delayed input. In [21] and later in [22],
Jaurige at al. show that adding a time-delayed version of the input can
considerably improve the performance of a reservoir whose hyperparameters are
not adjusted. This powerful approach requires to tune only two parameters
instead of the many hyperparameters usually associated to a reservoir, thus
allowing the optimisation with little computational cost. More recently, the
use of time delayed inputs was investigated numerically in [23] in the context
of a RC based on silicon microrings.

We test this approach for the first time on an experimental RC. We use the
architecture based on time multiplexing first introduced in [24], but in an
optoelectronic version as introduced in [25] and since studied extensively,
see e.g. [26, 27, 28]. We assess the impact of this technique on different
tasks involving time-series forecasting and audio classification, to verify
its effectiveness using data of various nature with different intrinsic
timescales. We show the superiority, in terms of performance, of this
procedure versus the standard hyperparameters optimisation.

We hope with this work to introduce a new technique in RC community for the
optimisation of physical reservoir computers; as discussed in the conclusion,
this approach could be particularly relevant for hardware reservoir computing
where it may not be possible to optimise all hyperparameters, or the
hyperparameters may drift over time.

##  2 Time-delay Reservoir Computing

A reservoir computer consists of three main components: an input layer, a
reservoir and an output layer. The input layer maps the input signal into the
reservoir. The reservoir is a recurrent neural network with fixed
interconnections. The reservoir layer exploits the rich dynamics of its
internal state, acting as a complex computational substrate that temporally
transforms input signals into high-dimensional representations. During the
training phase, the reservoir’s states are used to evaluate the parameters of
the linear output layer.

In this work, we use the well-known delay-based reservoir computer based on a
single nonlinear node. The reader can refer to [29, 15] for reviews on time-
delay RC. In this implementation, the reservoir nodes are arranged in a
topology similar to a ring structure where each node is connected to its
neighboring nodes. The dynamics of a time-delay reservoir of size $N$ in
discrete time is described by:

$\displaystyle x_{0}(n+1)=f(\alpha x_{N-1}(n-1)+\beta M_{0}u(n))$ (1)
$\displaystyle x_{i}(n+1)=f(\alpha x_{i-1}(n)+\beta M_{i}u(n))$ $\displaystyle
i=1,..,N-1,$

and the reservoir’s output $y$ is evaluated through a simple linear operation
$y(n)=W_{out}x(n)$ and it is used to evaluate the reservoir computer’s
performance, selecting the task-specific figure or merit. $x(n)$ is the
$N$-size vector of the reservoir states at timestep $n$; $u(n)$ is the
temporal dependent input signal of size $K$; $f$ is a nonlinear activation
function, in our experiments is a sinusoidal function (cf. sec. 4); $\alpha$
is the feedback attenuation, which represents the strength of the
interconnections between reservoir’s nodes; the input mask $M$ maps the input
signal in the reservoir, and its coefficient are usually drawn at random from
a uniform distribution in the range $[-1,+1]$; $\beta$ is the scaling of the
input signal, usually called input strength. After running the reservoir, all
the state vectors $x(n)$ are collected in a $N$x$K$ state matrix $X$ which is
used for the training, i.e. the to obtain the $N$ weights $W_{out}$ of the
linear readout layer. In this work we use the regularized linear regression
[43]:

$W_{out}=(X^{T}X+\lambda I)^{-1}X^{T}\tilde{y},$ (2)

where $\tilde{y}$ is the target signal and $\lambda$ is the regularization
parameter. In the case of a classification task with $C$ output classes,
$W_{out}$ becomes a $C$x$N$ matrix.

##  3 Reservoir Computing with Delayed Input

The idea of this work is to used a delayed version of the input time-series
$u(n)$ as the input signal to the reservoir. We use two $N$x$K$ input masks
$M_{1}$ and $M_{2}$ with values sampled from uniform distributions in the
range $[0,1]$. The value of the new input $J$ at timestep $n$ is then defined
as:

$J(n)=\beta_{1}u(n)M_{1}(n)+\beta_{2}u(n-d)M_{2}(n)+J_{0},$ (3)

where $d$ is the delay, $\beta_{1}$ and $\beta_{2}$ are scaling parameters,
and $J_{0}$ is a constant bias. The dynamics of the reservoir of equation (1)
become then:

$\displaystyle x_{0}(n+1)=f(\alpha x_{N-1}(n-1)+J(n+1))$ (4) $\displaystyle
x_{i}(n+1)=f(\alpha x_{i-1}(n)+J(n+1))$ $\displaystyle i=1,..,N-1.$

The new input $J$ contains two different components: one related to the input
at the present timestep, and one related to the input value $d$ timesteps in
the past. The masking process is now made with the two different masks $M_{1}$
and $M_{2}$. The idea is to optimise the reservoir tuning only the values of
$\beta_{2}$ and $d$, while keeping fixed the parameters specific to the
reservoir (such as $\alpha$ and $\beta_{1}$). A detailed explanation of this
procedure can be found in section 6.

##  4 Experimental Setup

Our optoelectronic setup, shown in Figure 1, implements the time-delay RC
introduced in section 2 and it is similar to experimental systems introduced
in previous works [26, 27].

Figure 1: Optoelectronic reservoir computer. In orange the optic part and in
blue the electronic part. SLD: superluminescent diode. MZ: Mach-Zender
intensity modulator. OA: Optical Attenuator. PD: Photo-Detector. FPGA: Field
Programmable Gate Array. Amp: amplifier.

A superluminescent diode (Thorlabs SLD1550P-A40) generates broadband light at
1550 nm. A Field Programmable Gate Array (FPGA) board generates the
reservoir’s input signal. This electrical signal coming from the FPGA drives
an electro-optic Mach-Zender (MZ) intensity modulator (EOSPACE AX2X2-0MSS-12):
the input signal is thus represented in the optical domain. Moreover, the
sinusoidal transfer function of the MZ modulator works in intensity only, so
that the MZ acts as
$I_{0}sin^{2}(V+\frac{\pi}{4})=I_{0}(\frac{1}{2}+\frac{1}{2}sin(2V))$. Because
the photodetector and the amplifier have a low pass filter, the constant term
has no role in the dynamics, leading to the equations (1), with $f$ as a
sinusoidal nonlinearity. The MZ modulator takes roughly 5 ns to reach a steady
state response, while the mask step duration is in the order of hundreds of
ns: it is therefore possible to describe the system using time discrete
equations.

The light passes trough an Optical Attenuator (JDS HA9) which attenuates the
light intensity in the loop by a fixed configurable factor. Then, the light
travels in a fiber spool, whose length of 1.7 km corresponds to a delay of
7.94 $\mu$s. This is the delay used to implement the delay-based RC: at this
point, after the spool, the light represents the time-multiplexed values of
reservoir’s nodes. Part of this light is collected by a Photo-Detector (PD),
electrically summed to the input coming from the FPGA, and amplified (with a
Mini Circuits ZHL-32A+ coaxial amplifier) to drive the MZ modulator over its
full $V_{\pi}$ range . The other part of the light is collected by another PD,
stored by the FPGA and saved on a PC for offline computations of the output
weights.

##  5 Tasks

We test our optoelectronic reservoir computer on four different tasks:
NARMA10, Mackey-Glass system, Spoken Digit Recognition and Speaker
Recognition. These are all widely used benchmark tasks in RC community. The
first two are time-series predictions tasks, whose goal is to forecast future
values in a sequence based on historical information; the last two are
classification tasks, aiming to accurately identify and categorize audio-based
temporal signals.

###  5.1 NARMA10

The NARMA10 task is a widely used benchmark in RC, and involves predicting the
next value in a sequence generated by a NARMA (Nonlinear Auto-Regressive
Moving Average) process of order 10. This task evaluates the capability of a
RC to capture complex temporal dependencies and make accurate predictions in
nonlinear dynamic systems. The input signal $u(n)$ is randomly drawn from a
uniform distribution between $[0,0.5]$. The output $q(n)$ of the NARMA10
system is defined as:

$q(n+1)=0.3q(n)+0.05q(n)(\sum_{i=0}^{9}q(n-1))+1.5u(n-9)u(n)+0.1.$ (5)

The aim is, given $u(n)$, to predict $q(n)$. We evaluate the accuracy of the
system using the normalised mean square error (NMSE), defined as

$NMSE=\frac{\langle(y(n)-\tilde{y}(n))^{2}\rangle}{\langle\tilde{y}(n)-\langle\tilde{y}(n)\rangle)^{2}\rangle},$
(6)

where $y(n)$ is the reservoir’s output and $\tilde{y}(n)$ is the target signal
($\tilde{y}(n)=q(n)$ for NARMA10).

###  5.2 Mackey-Glass system

The Mackey-Glass system describes a time-delay differential equation that
exhibits chaotic behavior [30]:

$\frac{dx}{dt}=\beta\frac{x(t-\tau)}{1+x(t-\tau)^{n}}-\gamma x.$ (7)

With $\beta=0.2$, $\tau=17$, $n=10$, $\gamma=0.1$ and timestep $dt=1$. In this
work we try to forecast the values of the series ten time steps into the
future. Similarly to NARMA10, also for this task we use the NMSE (eq. 6) as a
figure of merit.

###  5.3 Spoken Digit Recognition

The Spoken Digit Recognition dataset consists of 500 total utterances of the
10 (from 0 to 9) digits. Each digit is repeated ten times by five different
individuals. We use a dataset with a 3 dB Signal-To-Noise (SNR) ratio babble
noise, to increase the difficulty of the task and to evaluate the system’s
performance in noisy environments. The dataset is split 450/50 for train/test.
The audio signals are pre-processed using the Lyon Passive Ear model [31].
This model emulates the auditory canal’s biological response and transforms
the original time-domain utterances in a 86-channels frequency representation,
which is used as input signal. The length of the utterances ranges from 22 to
95 timesteps, with an average length of approximately 68 timesteps each.

In order to address the small size of our dataset containing 500 digits, we
implement k-fold cross-validation with k=10. This technique involves dividing
the dataset into 10 equal parts of 50 digits each. The training process is
then repeated 10 times, where each time a different subset of 50 digits is
used for testing while the remaining 450 are used for training. Since this is
is a classification task with 10 output classes, we use 10 distinct linear
classifiers: each classifier is designed to output "+1" if the input matches
the corresponding digit and "-1" otherwise. The system selects one of the
output classes as the winning one for every timestep of an utterance, and the
most voted class during the duration of an utterance is selected as the
prediction for that particular utterance: this approach is usually referred to
as winner-takes-all. Since this is a classification task and not a time-series
prediction task, we don’t use use anymore the NMSE as a figure of merit but
the Error Rate. The Error Rate is simply defined as the ratio between the
utterances predicted correctly and the total amount of utterances.

###  5.4 Speaker Recognition

The Japanese Vowels dataset [32] consists of 640 utterances of the Japanese
vowel ‘ae’, pronounced by 9 individuals. The goal is to correctly identify the
speaker for each utterance. To keep consistency with previous works, the
database is split in 270 train sequences and 370 test sequences. The audio
samples are pre-processed using the Mel-frequency cepstral coefficients
(MFCCs) to obtain their frequency representation, used as input signal to the
reservoir. Each sample consists thus of 12 MFCC coefficients for every
timestep. The length of the utterances ranges from 7 to 29 timesteps, with an
average length of approximately 15 timesteps each. For statistical purpose,
the training procedure is repeated 10 times: each time we select at random 270
audio samples for the training and use the rest for testing. Also here, as
with the Spoken Digit Recognition task, the Error Rate and the winner-takes-
all approach are used to evaluate the accuracy of the system.

##  6 Results

In this section we show the experimental results obtained with our
optoelectronic reservoir on the tasks listed in section 5. We use $N=50$
reservoir nodes for NARMA10 and Mackey-Glass, and $N=100$ nodes for the audio
tasks. For the NARMA10 and Mackey-Glass tasks, we use a time-series length of
$K=21000$. We divide these data points in four sets: the first 500 points are
removed for washout, the next 10000 points are used for training, then 500
points are again removed for washout, and the last 10000 points are used for
testing. For washout we mean the removing of initial or transient states in a
system, to ensure that the model focuses on the system’s steady-state behavior
rather than transient effects that might not represent the system’s long-term
dynamics. For the Digits and the Speakers tasks, the input length is given by
the dataset itself (cf. section 5). The frequency-encoded audio samples are
sent to the reservoir sequentially: the state of the reservoir are not reset
between every digit, meaning that the delayed input can combine data points
coming from the present and the previous digit.

To test the effectiveness of the optimisation using the delayed input, we use
the same approach for all the tasks. The input scaling $\beta_{1}$ and the
bias $J_{0}$ (cf. equation (3)) are selected for each task so that the range
of the input $J(n)$ falls in the range [0.2, 1.4] when $\beta_{2}$ is set to
zero. These values, though not optimized, are considered reasonable inputs for
the reservoir. The feedback attenuation is set with the OA of Figure 1 to 2 dB
(which roughly corresponds to a $\alpha=0.15$ in equation (1)), a reasonable
operating value for our optoelectronic reservoir, and the ridge regression
parameter is set to $\lambda=10^{-5}$ for all the tasks. We emphasize that we
do not optimise the value of $\alpha$, $\beta_{1}$ and $\lambda$ as
hyperparameters for every task, as it is usually done in RC. Another important
parameter that is not optimised is the feedback delay time. In fact, in this
work we want to highlight the reliability of using the delayed input for
optimisation without having to rely on the traditional tuning of
hyperparameters. Moreover, we test the NARMA10 and Mackey-Glass tasks also
with a value of feedback parameter of 15 dB (which roughly corresponds to
$\alpha=10^{-4}$), very far from reasonable values, to assess the
effectiveness of this approach when the hyperparameters are off from
acceptable values. Table 1 contains in a more compact form all the values
above mentioned, for all the tasks.

| $\beta_{1}$ | $J_{0}$ | Opt. att. | $\lambda$ | $N$
---|---|---|---|---|---
NARMA10 | 1.8 | 0.4 | 2 dB & 15 dB | $10^{-5}$ | 50
Mackey-Glass | 1 | 0 | 2 dB & 15 dB | $10^{-5}$ | 50
Spoken Digit | 10 | 0 | 2 dB | $10^{-5}$ | 100
Speaker | 0.3 | 0.8 | 2 dB | $10^{-5}$ | 100
Table 1: Parameters of the input and the reservoir, for all the considered
tasks. The optical attenuation (opt. att.) is set by the OA in Figure 1, and
is related to the hyperparameter $\alpha$ of equation (1).

The only two parameters that we need to tune for each task are then the
scaling of the delayed part of the input $\beta_{2}$ and the delay $d$ (cf.
equation (3)). We include in the search also the case where $\beta_{2}$ is
scanned while $d=0$, which gives a rough idea of how the reservoir behaves
when scanning the input strength with no delayed input (cf. section 6.2 for a
more detailed comparison of the optimisation with and without delayed input).

###  6.1 Effectiveness of optimisation with delayed input

Figure 2 show results on the NARMA10 task. We scan the parameters $\beta_{2}$
and $d$ looking for the best configuration. Our best results is a NMSE of 0.4:
it corresponds to a delay of 9 timesteps, most likely due to the intrinsic
timescales of the NARMA10 task (cf. equation (5)). As a comparison, in [27]
authors reach a best NMSE of 0.32 with a similar optoelectronic system, in all
likelihood because of the higher number of virtual nodes used ($N=80$ while we
use $N=50$). Moreover, the results in [27] are obtained through a more
advanced training method based on physical backpropagation, while we use a
simple scanning of two parameters. In [26] authors use an FPGA-based
optoelectronic reservoir similar to ours and reach a NMSE of 0.2, but using a
more complex training scheme based on an analog output layer and an online
learning algorithm. In [21] authors reach a NMSE of 0.3 using the delayed
input approach in numerical simulations.

Additionally, we wanted to test how the optimisation with delayed input works
in the case where the reservoir hyperparameters are rather poor. Figure 2(a)
and Figure 2(b) show results corresponding to reservoirs where a different
feedback attenuation $\alpha$ is applied: the attenuation is 2 dB in Figure
2(a) and 15 dB in Figure 2(b). As explained previously, with 2 dB our setup
can work reasonably on various tasks; instead, 15 dB is a high attenuation
value that brings the system in a operating region of poor performance.
Nonetheless, it is clear from Figure 2(b) that using a delayed input, for this
task, allows to reach the best performance even while operating in this sub-
optimal hyperparameter region. We also find quite poor performance without
delay ($\beta_{2}=0$) in both cases, most likely because the hyperparameters
of the reservoir (such as $\alpha$, $\beta_{1}$, $\lambda$) are not optimised.

Figure 2: Experimental results on the NARMA10 task as a function of
$\beta_{2}$ and the delay $d$ of equation (3), using different optical
feedback attenuation: 2 dB (a) and 15 dB (b).

Similar results, shown in Figure 3(a), are obtained on the prediction of the
Mackey-Glass system. Without delay ($\beta_{2}=0)$ the performance are already
acceptable (NMSE = 0.20), but they get further improved with the delay,
reaching an NMSE = 0.063 with $d=11$. Similar results were obtained on the
same task in numerical simulations in [21], using the delayed input approach.
Our results are also comparable to previous works [33, 34] which use similar
time-delay reservoirs with thousands of nodes. In the case of an attenuation
of 15 dB (Figure 3(b)), without delay the system performs very poorly, with
the NMSE exceeding 2; instead, by tuning the delayed input we reach an NMSE =
0.141. This proves that, even when very inadequate hyperparameters are
selected, the use of a delayed input can still bring the reservoir in a
valuable operating region. Note that the optimal delay $d$ is different from
the delay $\tau$ in equation (7) (while one would naively suppose they should
be similar).

Figure 3: Experimental results on the prediction of Mackey-Glass system 10
steps in the future, as a function of $\beta_{2}$ and the delay $d$ of
equation (3), using different optical feedback attenuation: 2 dB (a) and 15 dB
(b).

The results on "noisy" Spoken Digit Recognition are reported in Figure 4. For
this task the error rate is considerably improved using the delayed input,
decreasing from 0.250 with no delay ($\beta_{2}=0$) to 0.156 with a delay
$d=23$. In [35], on the same task and with the same optoelectronic setup, we
reported similar performance but using a deep architecture with 2
interconnected reservoir of size $N=100$ each, or using a single reservoir of
size $N=600$, while in this work we use only one reservoir of size $N=100$.

Figure 4: Experimental results on Spoken Digit Recognition with noise, as a
function of $\beta_{2}$ and the delay $d$ of equation (3).

Figure 5 shows results on the Speaker Recognition. Also in this case the use
of a delayed input increase the performance. The error rate decreases from
0.0384 with no delay ($\beta_{2}=0$) to 0.0170 error with a delay $d=5$. As a
comparison, in [36] authors reach a best error of roughly 0.020 with an
optical reservoir of $N=200$ nodes, whereas in [37] authors report an error of
roughly 0.035 with a reservoir size comparable to ours.

For these two latter tasks involving audio processing, the explanation for the
optimum values of the delay is not straightforward. For both tasks, the
optimum delay value is shorter than the average duration of an utterance: the
Digits task has an average sample length of 68 timesteps and an optimum delay
$d=23$, while the Speakers task has an average sample length of 15 timesteps
and an optimum delay $d=5$. A delay shorter than the utterance duration means
that most of the times the delayed input term and the non-delayed input term
refer to the same utterance. Thus, these results suggest that the performance
of the reservoir can be increased when different temporal components of the
same utterance are combined. While $d$ and the sample duration seem somewhat
correlated, their relationship does not look immediately quantifiable and it
is not possible so far to estimate the optimum delay a priori. This may be due
to the different audio preprocessing algorithms used for the two datasets.

Figure 5: Experimental results on Speaker Recognition, as a function of
$\beta_{2}$ and the delay $d$ of equation (3).

###  6.2 Optimisation with delayed input vs standard hyperparameters
optimisation

The previous section reports how, for all the tasks, the use of a delayed
input outperforms the case where no delay is used, i.e. $\beta_{2}=0$.
However, all the results presented so far refer to a reservoir whose
hyperparameters are not optimised, but rather set to a reasonable value. What
happens instead if the hyperparameters are optimized for each task ? Would the
superiority of the optimisation with delayed input still hold true? We try to
answer this question with Figure 6. The purple curves represent the "standard"
approach in RC with no delay, i.e. tuning the hyperparameters (in this case
$\alpha$, $\beta_{1}$ and $\lambda$) for every task. The orange curves
represent the optimisation using the delayed input, by setting the parameters
to the reasonable values of Table 1 and tuning only the delay parameters
$\beta_{2}$ and $d$. In both cases, we explore the behaviour of the system
when sweeping the optical attenuation (expressed in dB), which is related to
the feedback attenuation $\alpha$.

For all the tasks, it is clear that the optimisation with a delayed input
prevails: its accuracy is higher for every value of attenuation considered.
Especially for an optical attenuation larger than 10 dB, where performance
without delay are very bad, the use of a delay ensures accuracy remains high.

Figure 6: Experimental results on all the tasks, to compare the standard
approach (without delay) and the approach proposed in this work (with delay).
The purple curves are obtained using no delay, optimising the hyperparameters
$\beta_{1}$ and $\lambda$ and then sweeping the optical feedback attenuation.
The orange curves are obtained using the delayed input, optimising the
parameters $\beta_{2}$ and $d$ of equation (3), and then sweeping the optical
feedback attenuation.

##  7 Conclusions

In this work we investigate a novel approach for the optimisation of physical
reservoir computers using a delayed input, introduced numerically in [21]. We
verify its effectiveness using an optoelectronic setup on different benchmark
tasks, involving time-series prediction and audio signal classification.

We emphasize the advantages of this approach. It is a very simple technique
that does not require complex search algorithms. In fact, the reservoir’s
performance can be optimised only by tuning two parameters: the strength of
the delayed component ($\beta_{2}$) of the input, and the delay itself ($d$).
Therefore, it is not necessary to carefully select the reservoir’s
hyperparameters for every task: it is sufficient to set them to some
reasonable value and fine tune the delayed input. In addition, the performance
is not very sensitive on the value of $\beta_{2}$: values in the range 1-3 are
good, and need a bit of refinement to get the best values. Since one only
needs to act on the input signal, this approach can be useful for physical
reservoirs whose hyperparameters cannot be tuned due to real-world
constraints. Furthermore, even in the case where the hyperparameters are
considerably off from acceptable values, the use of a delayed input can still
bring the reservoir to satisfactory performance. Moreover, this optimisation
procedure outperforms the traditional hyperparameters tuning in terms of
performance.

Some open questions which we leave open for future research concern the
relationship between the timescale of a temporal task and the optimum amount
of delay to employ; the use of multiple delays (in which case improved
optimisation schemes, such as Bayesian optimisation, rather than the simple
grid scan used here, would be useful); and on the experimental side the
implementation of the delayed input with a physical delay line.

Funding E.P and S.M. acknowledge financial support from the H2020 Marie
Skłodowska-Curie Actions (Project POSTDIGITAL Grant number 860360), and from
the Fonds de la Recherche Scientifique - FNRS, Belgium under funding EOS n°
40007536 and CDR n° 40021243. L.J. and K.L. acknowledge funding from the
Deutsche 369 Forschungsgemeinschaft (DFG), Grant No. LU 1729/3-1 370
(Projektnummer 445183921).

## References

  * [1] H. Jaeger and H. Haas, “Harnessing nonlinearity: Predicting chaotic systems and saving energy in wireless communication,” Science 304, 78–80 (2004).
  * [2] I. Estébanez, S. Li, J. Schwind, _et al._ , “56 gbaud pam-4 100 km transmission system with photonic processing schemes,” Journal of Lightwave Technology 40, 55–62 (2022).
  * [3] D. Verstraeten, B. Schrauwen, D. Stroobandt, and J. Van Campenhout, “Isolated word recognition with the liquid state machine: a case study,” Information Processing Letters 95, 521–528 (2005).
  * [4] F. Triefenbach, A. Jalalvand, B. Schrauwen, and J.-P. Martens, “Phoneme recognition with large hierarchical reservoirs,” Advances in Neural Information Processing Systems 23 (2010).
  * [5] Y. D. Mammedov, E. U. Olugu, and G. A. Farah, “Weather forecasting based on data-driven and physics-informed reservoir computing models,” Environmental Science and Pollution Research pp. 1–14 (2022).
  * [6] N. Schaetti, M. Salomon, and R. Couturier, “Echo state networks-based reservoir computing for mnist handwritten digits recognition,” in _2016 IEEE Intl Conference on Computational Science and Engineering (CSE) and IEEE Intl Conference on Embedded and Ubiquitous Computing (EUC) and 15th Intl Symposium on Distributed Computing and Applications for Business Engineering (DCABES),_ (2016), pp. 484–491.
  * [7] P. Antonik, N. Marsal, and D. Rontani, “Large-scale spatiotemporal photonic reservoir computer for image classification,” IEEE Journal of Selected Topics in Quantum Electronics 26, 1–12 (2019).
  * [8] P. Antonik, N. Marsal, D. Brunner, and D. Rontani, “Human action recognition with a large-scale brain-inspired photonic computer,” Nature Machine Intelligence 1, 530–537 (2019).
  * [9] E. Picco, P. Antonik, and S. Massar, “High speed human action recognition using a photonic reservoir computer,” Neural Networks 165, 662–675 (2023).
  * [10] G. Tanaka, T. Yamane, J. B. Héroux, _et al._ , “Recent advances in physical reservoir computing: A review,” Neural Networks 115, 100–123 (2019).
  * [11] K. Nakajima and I. Fischer, _Reservoir Computing_ (Springer, 2021).
  * [12] A. Liutkus, D. Martina, S. Popoff, _et al._ , “Imaging with nature: Compressive imaging using a multiply scattering medium,” Scientific Reports 4, 1–7 (2014).
  * [13] M. Rafayelyan, J. Dong, Y. Tan, _et al._ , “Large-scale optical reservoir computing for spatiotemporal chaotic systems prediction,” Physical Review X 10, 041037 (2020).
  * [14] L. Larger, A. Baylón-Fuentes, R. Martinenghi, _et al._ , “High-speed photonic reservoir computing using a time-delay-based architecture: Million words per second classification,” Physical Review X 7, 011015 (2017).
  * [15] G. Van der Sande, D. Brunner, and M. C. Soriano, “Advances in photonic reservoir computing,” Nanophotonics 6, 561–576 (2017).
  * [16] J. Yperman and T. Becker, “Bayesian optimization of hyper-parameters in reservoir computing,” arXiv preprint arXiv:1611.05193 (2016).
  * [17] P. Antonik, N. Marsal, D. Brunner, and D. Rontani, “Bayesian optimisation of large-scale photonic reservoir computers,” Cognitive Computation pp. 1–9 (2021).
  * [18] A. A. Ferreira, T. B. Ludermir, and R. R. de Aquino, “An approach to reservoir computing design and training,” Expert Systems with Applications 40, 4172–4182 (2013).
  * [19] S. Basterrech and T. K. Sharma, “Re-visiting reservoir computing architectures optimized by evolutionary algorithms,” in _International Conference on Soft Computing and Pattern Recognition,_ (Springer, 2022), pp. 832–843.
  * [20] H. Joy, M. Mattheakis, and P. Protopapas, “Rctorch: a pytorch reservoir computing package with automated hyper-parameter optimization,” arXiv preprint arXiv:2207.05870 (2022).
  * [21] L. Jaurigue, E. Robertson, J. Wolters, and K. Lüdge, “Reservoir computing with delayed input for fast and easy optimisation,” Entropy 23, 1560 (2021).
  * [22] L. Jaurigue and K. Lüdge, “Reducing reservoir computer hyperparameter dependence by external timescale tailoring,” Neuromorphic Computing and Engineering (2024).
  * [23] B. J. G. Castro, C. Peucheret, and F. Da Ros, “Wavelength-multiplexed delayed inputs for memory enhancement of microring-based reservoir computing,” arXiv preprint arXiv:2312.04204 (2023).
  * [24] L. Appeltant, M. C. Soriano, G. Van der Sande, _et al._ , “Information processing using a single dynamical node as complex system,” Nature Communications 2, 468 (2011).
  * [25] Y. Paquot, F. Duport, A. Smerieri, _et al._ , “Optoelectronic reservoir computing,” Scientific Reports 2, 287 (2012).
  * [26] P. Antonik, F. Duport, M. Hermans, _et al._ , “Online training of an opto-electronic reservoir computer applied to real-time channel equalization,” IEEE Transactions on Neural Networks and Learning Systems 28, 2686–2698 (2016).
  * [27] M. Hermans, P. Antonik, M. Haelterman, and S. Massar, “Embodiment of learning in electro-optical signal processors,” Physical Review Letters 117, 128301 (2016).
  * [28] M. Nakajima, K. Inoue, K. Tanaka, _et al._ , “Physical deep learning with biologically inspired training method: gradient-free approach for physical hardware,” Nature Communications 13, 7847 (2022).
  * [29] T. Hülser, F. Köster, L. Jaurigue, and K. Lüdge, “Role of delay-times in delay-based photonic reservoir computing,” Optical Materials Express 12, 1214–1231 (2022).
  * [30] M. C. Mackey and L. Glass, “Oscillation and chaos in physiological control systems,” Science 197, 287–289 (1977).
  * [31] R. Lyon, “A computational model of filtering, detection, and compression in the cochlea,” in _ICASSP’82. IEEE International Conference on Acoustics, Speech, and Signal Processing,_ vol. 7 (IEEE, 1982), pp. 1282–1285.
  * [32] M. Kudo, J. Toyama, and M. Shimbo, “Multidimensional curve classification using passing-through regions,” Pattern Recognition Letters 20, 1103–1111 (1999).
  * [33] S. Ortín, M. C. Soriano, L. Pesquera, _et al._ , “A unified framework for reservoir computing and extreme learning machines based on a single time-delayed neuron,” Scientific reports 5, 14945 (2015).
  * [34] P. Antonik, M. Gulina, J. Pauwels, and S. Massar, “Using a reservoir computer to learn chaotic attractors, with applications to chaos synchronization and cryptography,” Physical Review E 98, 012215 (2018).
  * [35] E. Picco, A. Lupo, and S. Massar, “Deep photonic reservoir computer for speech recognition,” (2023).
  * [36] U. Paudel, M. Luengo-Kovac, J. Pilawa, _et al._ , “Classification of time-domain waveforms using a speckle-based optical reservoir computer,” Optics Express 28, 1225–1237 (2020).
  * [37] M. Dale, R. F. Evans, A. Sebald, and S. Stepney, “Reservoir computing with nanowire exchange-coupled spin torque oscillator arrays,” in _International Conference on Unconventional Computation and Natural Computation,_ (Springer, 2023), pp. 64–78.


# GeV muon beams with picometer-class emittance from electron-photon
collisions

C. Curatolo  Università degli Studi di Padova, via Marzolo 8, 35131 Padova,
Italy  INFN-Padova, via Marzolo 8, 35131 Padova, Italy  L. Serafini  INFN-
Milan, via Celoria 16, 20133 Milano, Italy

(6 June 2021)

###### Abstract

One of the challenge of future muon colliders is the production of muon beams
carrying high phase space densities. In particular the muon beam normalised
transverse emittance is a relevant figure of merit to meet luminosity
requests. A typical issue impacting the achieved transverse emittance in muon
collider schemes so far considered is the phase space dilution caused by
coulomb interaction of primary particles propagating into the target where
muons are generated. In this study we present a new scheme for muon beam
generation occurring in vacuum by interactions of electron and photon beams.
Setting the center of mass energy at about twice the threshold (i.e. around
$350$ MeV) the normalised emittance of the muon beam generated via muon pair
production reaction ($e^{-}+\gamma\rightarrow e^{-}+\mu^{+}/\mu^{-}$) is
largely independent on the emittance of the colliding electron beam and is set
basically by the excess of transverse momentum in the muon pair creation. In
absence of any other mechanism for emittance dilution, the resulting muon
beam, with energy in the range of few tens of GeV, is characterised by an
ultra-low normalised transverse rms emittance of a few nm rad, corresponding
to a geometrical emittance below $10$ pm rad. This opens the way to a new muon
collider paradigm based on muon sources conceived with primary colliding beams
delivered by $100$ GeV-class energy recovery linacs interacting with hard-X
ray free electron lasers. The challenge is to achieve the requested luminosity
of the muon collider adopting a strategy of low muon fluxes/currents combined
to ultra-low emittances, so to largely reduce also the levels of muon beam-
induced background.

##  I Introduction

Muon beam generation for muon colliders has been traditionally conceived using
hadronic interactions that go through the pion production channel with
consequent decay into muons: the total cross section for these reactions is
quite large, assuring a wealth of muon population when intense proton beams
are impinging on targets made of proper material MAP . Unfortunately the
coulomb interaction of primary and secondary beams into target largely
disrupts the phase space of the generated muon beam, diluting its emittance
well above the upper limits set by collider luminosity requests, so that a
challenging muon cooling process must be implemented via ionisation with the
aim to restore the emittance levels down to values compatible for collider
luminosity ref:mice . Alternatively, positron-electron annihilation into muon
pairs has been proposed as a promising scheme to achieve directly at the muon
source low emittance values LEMMA . Investigations on this scheme are still in
progress to assess its full potentiality: nevertheless the primary beam of
$45$ GeV positrons has to interact with a solid target in order to achieve
sufficient muon pair production rate, which in turns poses issues of target
handling. An exotic scheme of Hadron Photon Collider (HPC) was proposed to
generate muon pairs through the reaction $p^{-}+\gamma\rightarrow
p^{-}+\mu^{+}/\mu^{-}$: in this case the emittance of the muon beam was quite
good, basically conserving the emittance of colliding proton and photon beams
HPC ; HPC2 ; HPC3 . The very small total cross section typical of muon photo-
production reactions sets a severe limit to the muon flux achievable. In this
paper we revisit the basic concepts of that study, implementing the idea of
muon photo-production in a different environment - that of electron-photon
collisions. The reason why lays not really in a larger cross section of the
$e,\,\gamma$ reaction with respect to $p,\,\gamma$ (the increase is quite
modest), but in the possibility to run an electron-photon collider at much
larger luminosity than a proton-photon collider as in the HPC scheme. As a
matter of fact the focal spot size of a $100$ GeV-class electron beam
colliding in-vacuum with a counterpropagating photon beam can be much smaller
than that of a TeV proton beam stored in a machine like LHC, as well as that
of an electron-positron collider final focus. The former is actually limited
by the minimum beta achievable in a proton storage ring due to cm-long proton
bunches, the latter is actually limited by the beam-beam interaction mutuated
by the presence of colliding charged particles focused into a tight focal
region (the so called beamstrahlung effect). In case of electron-photon
collisions none of these effect are actually present - electrons and photons
interact just through 2-particle QED reactions without any significant
collective e.m. field effect impacting the beam propagation through the
interaction region. So the maximum focusing achievable in the collision region
is just dictated by final focus beam optics, electron and photon bunch lengths
and transverse emittances of the two colliding beams. In order to reach the
Center of Mass (CM) available reaction energies in excess of $300$ MeV (where
the total cross section exceeds few hundreds of nano barns) we are forced to
consider electrons with energies higher than $100$ GeV and photons with energy
in the $100-300$ keV range. Following the analysis pursued in HPC studies, in
order to maximise the luminosity of an electron-photon collider in this range
of particle-photon energies we are naturally addressed to consider a hard
X-ray Free Electron Laser (FEL) as the brightest radiation source and an
Energy Recovery Linac (ERL) to obtain the maximum electron beam current
(A-class) combined with high brightness - namely short electron bunches (fs-
class) with very small transverse normalised emittances (below $1$ mm mrad).
Since no concern is raised by beam-beam interaction effects, maximum
luminosity is achieved running round beams at the collision aiming at minimum
beta values at the final focus. We will show in the following discussion that
close to state of the art performances both of FELs and ERLs are indicating
the possibility to achieve an outstanding luminosity in the $e^{-}-\gamma$
collider, up to $10^{41}$ cm-2s-1, which allows to generate up to $10^{11}$
muon pairs per second with geometrical emittances down to 5 pm rad. The paper
is organised as follows: in Section II we summarize the theoretical aspects
regarding the electron-photon beams collision and the main formulas we will
use to analyse the various cases reported in Section III. In the latter we
discuss the performances requests to primary beams in different scenarios of
$e^{-}-\gamma$ colliders, either using parasitically or semi-parasitically the
layout of other projects under proposal in the High Energy Physics (HEP)
future scenarios (CLIC CLIC , FCC FCC ), or conceiving a dedicated scenario of
colliding ERLs (Section III.1). In Section III.2 we discuss how to achieve the
requested FEL performances, running dedicated ERLs to drive the hard X-ray
FEL, while in Section III.3 the power consumption demands are summarized. In
Section IV we describe the results of Monte Carlo dedicated simulations to
assess the phase space distributions and quality of the generated muon beams,
its flexibility and dynamic range. Conclusions roughly depict the way towards
a TeV muon collider based on $e^{-}-\gamma$ collider muon source.

##  II Theory

We analyse here the most important aspects concerning the electron-photon
beams collision at center of mass energies around 250-450 MeV.
The center of mass energy is given by

$E_{CM}=\sqrt{2E_{e}h\nu-2(\underline{p}_{e}\cdot\underline{k})+M_{e}^{2}}$
(1)

where $E_{e},\,\underline{p}_{e}$ and $h\nu,\,\underline{k}$ are electron and
photon energies and momenta respectively and $M_{e}=0.511$ MeV/c2 is the
electron mass (natural units are used i.e. c=1). Assuming a head-on collision
within the electron and the incident photon with energy $h\nu<<E_{e}$ and
supposing $M_{e}<<E_{e}h\nu$, the Lorentz factor of the center of mass is

$\gamma_{CM}=\frac{E_{tot}^{LAB}}{E_{CM}}\simeq\frac{E_{e}+h\nu}{\sqrt{4E_{e}h\nu+M_{e}^{2}}}\simeq\frac{E_{e}}{2\sqrt{E_{e}h\nu}}=\frac{\sqrt{E_{e}}}{2\sqrt{h\nu}}$
(2)

The CM threshold energy for the Muon Pair Production (MPP) is given by

$E_{CM}^{th}=M_{e}+2M_{\mu}$ (3)

If for example $E_{e}=200$ GeV, the equality
$E_{CM}^{th}=M_{e}+2M_{\mu}=\sqrt{4E_{e}h\nu+M_{e}^{2}}$ implies that
$h\nu^{th}=56.6$ keV. The photon energy in the electron rest frame is
$h\nu^{th^{\prime}}=2\gamma_{e}h\nu^{th}=44.3$ GeV while the photon energy in
the CM frame is $h\nu^{th*}=2\gamma_{CM}h\nu^{th}=106.5$ MeV.
Beside MPP, the other predominant reactions in the mentioned CM energy range
are Triplet Pair Production (TPP - $e^{-}+\gamma\rightarrow
e^{-}+e^{+}/e^{-}$) and Inverse Compton Scattering (ICS -
$e^{-}+\gamma\rightarrow e^{-}+\gamma$).

Figure 1: Total cross sections of possible reactions in electron-photon
collision as function of $E_{cm}$ (MeV).

Their total cross sections $\sigma_{tot}$ as a function of $E_{CM}$ are
depicted in Fig.1. The TPP cross sections calculated as in Ref. motz  is

$\sigma^{TPP}_{tot}=\sigma^{TPP}_{Bethe-Heitler}-\sigma^{TPP}_{Borsellino}$
$\begin{split}\begin{cases}\sigma^{TPP}_{B-H}=\alpha
r_{e}^{2}&(3.11\,\ln(2k_{e})-8.07)\\\ \sigma^{TPP}_{Bors}=\frac{\alpha
r_{e}^{2}}{k_{e}}&\left(\frac{4}{3}\,\ln(2k_{e})^{3}-3\,\ln(2k_{e})^{2}+\right.\\\
&\left.\mbox{ }+6.84\,\ln(2k_{e})-21.51\right)\end{cases}\end{split}$ (4)

where $\alpha=1/137$ is the fine-structure constant, $r_{e}^{2}=0.079$ barn is
the classical electron radius squared and $k_{e}=2\gamma_{e}h\nu/M_{e}$.
In a similar way, the MPP cross section can be obtained from the one above by
substituting $r_{e}$ with $r_{\mu}=r_{e}(M_{e}/M_{\mu})$ and $k_{e}$ with
$k_{\mu}=2\gamma_{e}h\nu/M_{\mu}$. Its value is well-defined if
$E_{CM}^{2}/M_{\mu}^{2}>9$ and we also find that
$E_{CM}=\sqrt{2k_{\mu}M_{\mu}M_{e}}$.
The ICS cross section, as function of the electron recoil
$X=2h\nu^{\prime}/M_{e}=4\gamma_{e}h\nu/M_{e}$, is given by ics

$\begin{split}\sigma&{}^{ICS}_{tot}=2\pi r_{e}^{2}\frac{1}{X}\cdot\\\
&\cdot\left[\left(1-\frac{4}{X}-\frac{8}{X^{2}}\right)\ln(1+X)+\frac{1}{2}+\frac{8}{X}-\frac{1}{2(1+X)^{2}}\right].\end{split}$
(5)

The luminosity of the $e^{-}-\gamma$ collider is defined as

$\mathcal{L}=\frac{N_{e}\,N_{ph}\,r}{2\,\pi\,\sqrt{\sigma_{x_{e}}^{2}+\sigma_{y_{e}}^{2}}\sqrt{\sigma_{x_{ph}^{2}}+\sigma_{y_{ph}}^{2}}}\,[\mbox{cm}^{-2}\mbox{
s}^{-1}]$ (6)

being $N_{e},N_{ph}$ the number of electron and photons per bunch, $r$ the
repetition rate of the collisions and
$\sigma_{x_{e}},\sigma_{y_{e}},\sigma_{x_{ph}},\sigma_{y_{ph}}$ the transverse
dimensions of the electron and the photon beams respectively. Therefore the
number of muon and electron pairs and ICS events per second are

$\begin{split}\mathcal{N}_{\mu^{\pm}}=\mathcal{L}\cdot\sigma^{MPP}_{tot}(E_{CM})\,[\mbox{s}^{-1}]\\\
\mathcal{N}_{e^{\pm}}=\mathcal{L}\cdot\sigma^{TPP}_{tot}(E_{CM})\,[\mbox{s}^{-1}]\\\
\mathcal{N}_{ICS}=\mathcal{L}\cdot\sigma^{ICS}_{tot}(E_{CM})\,[\mbox{s}^{-1}]\end{split}$
(7)

The transverse normalized emittance of the produced muons is determined by the
intrinsic thermal contribution of the reaction and by the features of the
incoming electron beam. It can be described by

$\epsilon_{\mu}^{n}\simeq\frac{2}{3}\sigma_{0}\frac{\sqrt{E_{e}\,h\nu}-M_{\mu}}{M_{\mu}}+\frac{\langle\gamma_{\mu}\rangle\epsilon_{e}^{n}}{\gamma_{e}}$
(8)

where $\epsilon_{e}^{n}$ is the incoming electron beam transverse normalized
emittance, $\langle\gamma_{\mu}\rangle$ the mean energy of the muon beam,
$\sigma_{0}=\sqrt{\left(\sqrt{\sigma_{x_{e}}^{2}+\sigma_{y_{e}}^{2}}\sqrt{\sigma_{x_{ph}}^{2}+\sigma_{y_{ph}}^{2}}\right)/2}$
and the first addendum represents the normalized thermal emittance of the muon
beam. In the cases we will consider, the emittance of the incoming electron
beam has basically no impact on the muon beam one: the second part of Eq. 8 is
negligible because $\langle\gamma_{\mu}\rangle\ll\gamma_{e}$.

##  III $e^{-}-\gamma$ collider scenarios

We have conceived a few scenarios for $e^{-}-\gamma$ colliders achieving ultra
high luminosities, considering various technologies for the electron
accelerators involved in the generation of the primary electron beam and in
the generation of the FEL hard X-ray photon beam. An overarching goal is the
quest of very high sustainability: such levels of luminosities typically
require beam power levels that are hardly compatible with reasonable impacts
on society and governmental approvals when power consumption of the whole
scenario start to exceed the so called GW scale, that should be considered as
an unsurpassable borderline. The resulting criterion for the choice of the
electron accelerators is the following: since the generation of muon pairs in
the collision between the primary electrons and the FEL photons is
perturbative on the primary beams, with negligible losses of beam particles
and very small energy loss in the beams, a strategy for beam power
recirculation or recovery is mandatory. The main source of beam power loss in
these scenarios is set by the amount of synchrotron radiation power generated
at any electron beam deflection or recirculation. Four scenarios have been
analysed and described in detail in Table 1: two of them are based on
parasitic use of machines that are already in a proposal stage in the HEP
future strategy, i.e. CLIC and FCC-ee. Two further scenarios are based on
PERLE-like ERLs perle1 ; perle2 ): one, named FCC-twin Linac, is integrated in
FCC tunnels system, the other, named Twin Linac, is based on an ILC schematics
ILC ; shil . Only these scenarios satisfy the luminosity requirements but we
report for sake of comparison and point of reference also the performances of
parasitic operations of CLIC and FCC-ee.

Figure 2: FCC-twin Linac scheme: FCC ring with ERL insertions. Primary $e^{-}$ accelerated up to $100$ GeV collide with the counterpropagating FEL producing $\mu^{\pm}$. Both $e^{-}$ beams (primary and FEL) energies are recovered. $\mu^{\pm}$ accumulated and accelerated in a separate ring. Figure 3: Twin Linac scheme. Primary $e^{-}$ accelerated up to $200$ GeV collide with the counterpropagating FEL producing $\mu^{\pm}$. Both $e^{-}$ (primary and FEL) are decellerated in the opposite linacs and the energy recovered. A selected fraction of $\mu^{\pm}$ are injected in the opposite linac and accelerated before stogare and collision in a ring. Table 1: Parameters table of the four schemes: FCC-twin Linac, CLIC parasitic, FCC-ee parasitic, Twin Linac. Primary beams features and peak luminosities are reported. | FCC-twin Linac | CLIC parasitic | FCC-ee parasitic | Twin Linac
---|---|---|---|---
Energy $e^{-}$ beam (GeV) | $100$ | $1500$ | $200$ | 200
Bunch charge (pC) | 250 | 580 | 43200 | 250
Electrons per bunch | 1.6 $10^{9}$ | 3.7 $10^{9}$ | 2.7 $10^{11}$ | 1.6 $10^{9}$
Repetition rate (MHz) | 800 (CW) | 0.005 | 0.1 | 800 (CW)
Average Current (mA) | 200 | 0.003 | 5.4 | 2$\times$200
Nominal beam power (GW) | 20 | 0.0045 | n.a. | 2$\times$40
Energy loss 1 turn @ 15 km r (MeV) | 593 | n.a. | 9500 | n.a.
Synch. rad. 1 turn @ 15 km r (MW) | 119 | n.a. | 51 | n.a.
Beam power recovery fraction | 99.9% | n.a | n.a. | 99.9%
Beam power loss (MW) | 20 | 4.5 | n.a. | 2$\times$40
Bunch length (psec) | 0.3 | 0.15 | 12 | 0.3
$\epsilon^{n}_{x}$ (m rad) | 4 $10^{-7}$ | 6.6 $10^{-7}$ | 520 $10^{-6}$ | 4 $10^{-7}$
$\epsilon^{n}_{y}$(m rad) | 4 $10^{-7}$ | 3 $10^{-8}$ | $10^{-6}$ | 4 $10^{-7}$
$\beta_{x}$ (mm) | 0.2 | 0.6 | 1.6 | 0.2
$\beta_{y}$ (mm) | 0.2 | 12 | 1.6 | 0.2
$\sigma_{x}$ (nm) | 20 | 10 | 1442 | 14
$\sigma_{y}$ (nm) | 20 | 10 | 63 | 14
FEL photon energy (keV) | 300 | 30 | 75 | 150
Photons per pulse | $10^{12}$ | 1.6 $10^{14}$ | 2 $10^{13}$ | 5 $10^{12}$
$\epsilon_{x,y}$ (m rad) | 3 $10^{-13}$ | 3 $10^{-12}$ | 1.2 $10^{-12}$ | 6 $10^{-13}$
Focal spot size (nm) | 20 | 10 | (1400 $\times$ 60) | 14
Repetition rate (MHz) | 800 | 0.005 | 0.1 | 800
FEL beam power (MW) | 40 | 0.00375 | 0.025 | 2$\times$100
FEL $\rho$ | 3 $10^{-4}$ | $10^{-3}$ | 7 $10^{-4}$ | 5 $10^{-4}$
FEL efficiency (tapering) | 0.5% | 5% | 2% | 1%
FEL $e^{-}$ beam av. curr. (mA) | 200 | 0.0025 | 0.025 | 200
FEL $e^{-}$ beam bunch ch. (pC) | 250 | 500 | 250 | 250
FEL $e^{-}$ beam energy (GeV) | 70 | 30 | 50 | 50
FEL $e^{-}$ beam power (MW) | 14000 | 0.075 | 1.25 | 2$\times$10000
Beam power recovery fraction | 99.9% | n.a. | n.a. | 99.9%
Beam power loss (MW) | 14 | 0.075 | 1.25 | 2$\times$10
Total beam power loss (MW) | 202 | n.a. | n.a. | 200
Peak Luminosity $e^{-}-\gamma$ (cm-2s-1) | 2.5 $10^{40}$ | 2 $10^{38}$ | 3 $10^{35}$ | 2$\times$2.5 $10^{41}$

The final goal of generating a suitable beam of muons to be used in a TeV-
scale muon collider is accomplished by the Twin Linac scenario based on a twin
$200$ GeV ERL system coupled to a twin $50$ GeV FEL ERL system, with residual
beam power loss (after recovery, and taking into account all losses of
electron beam and FEL photon beam) of about $200$ MW, together with an
outstanding $4.5$ nm rad normalised emittance of the muon beam at $10^{11}$
muon pairs per second.

###  III.1 The primary electron beam

We illustrate here the rationale of the first and fourth column listed in
Table 1 for the primary electron beam parameter case. As a matter of fact the
intermediate cases, based on CLIC and FCC-ee, are just taken from electron
beam parameter sets reported in the literature and they are listed just to
show how far they are from the target value for luminosity. Instead, the first
case, based on a ERL integrated in the FCC tunnel (with proper tunnel
reshaping to achieve a race-track layout), and the fourth case, based on a
dedicated twin Linac layout, are both capable to meet the tight demands on
luminosity.

An important consideration concerns the peak properties fo the primary
electron beam for both these FCC-tL and TL cases: the parameters listed in the
table are quite consistent with present state of the art for electron beams,
with a rms normalised transverse emittance ($0.4$ mm mrad, round beam) quite
consistent with the accelerated bunch charge ($250$ pC). Also the value chosen
for the beta function at the collision point ($0.2$ mm) is very close to state
of the art performances. That allows to match the spot size of the FEL photon
beam down to $14$ nanometers.

The power budget for the primary electron beams is as follows: in case FCC-tL
the power loss is dominated by the synchrotron emission along the two
recirculating arcs (basically one turn of the FCC ring with about $15$ km
curvature radius), summing up to $119$ MW. In addition, considering a recovery
efficiency of $99.9\%$ in the ERL, further $20$ MW of electron beam power is
lost, bringing to a total beam power loss of about $139$ MW. In the Twin Linac
scenario no significant synchrotron radiation is expected, therefore the
primary beam power loss is set by the ERL efficiency, bringing up to $2\times
40=80$ MW.

###  III.2 High efficiency hard X-ray FELs

The photon beam needed to achieve ultra high luminosity in the $e^{-}-\gamma$
collider is unique: it must carry an outstanding number of photons per pulse
at the same repetition rate of the primary electron beam. It must also match
the ultra tight focus spot size at collision of the primary electron beam, set
by its very short beta function value at the focus (hundreds of microns), in
the range of a few tens of nanometers. The only radiation source capable to
meet these demanding requirements is a FEL driven by a dedicated ERL, and
operated in SASE mode with tapering, as illustrated in Ref. emma . Efficiency
in the range of a few percent is achievable, bringing to a number of photons
per pulse as listed in Table 1, according to photon energy. The partial
coherence of the amplified FEL radiation makes also possible to focus its
photon beam down to nanometer spot sizes, as discussed in Ref. matsu . This is
the second crucial property of FELs that makes possible to foresee a
luminosity for the $e^{-}-\gamma$ collider best scenario (based on a system of
two twin ERLs at $200$ GeV electron energy, combined to two twin $50$ GeV ERLs
driving the FELs) of up to $2\times 2.5\,10^{41}$ cm-2s-1.

The FEL beams considered in this study carry an impressive amount of photon
beam power: running at $800$ MHz in CW mode, the photon number per second
exceeds $10^{21}$. With photon energies in the range of tens to hundreds keV
that means up to $100$ MW radiation beam power. Since the FEL efficiency,
considering the special mode of FEL operation as illustrated in Ref. emma , is
about $1\%$, the power carried by the electron beam driving the FEL must be of
the order of $100$ MW/$0.01=10$ GW, as listed in Table 1. Therefore, the power
budget for the Linac driving the FEL is as follow: in case FCC-tL we have a
synchrotron radiation accounting for $29$ MW, a recovery beam efficiency for
$14$ MW and a FEL beam power loss of $40$ MW. Assuming that $50\%$ of the
radiation power may be recovered back as AC power, the total power budget is
$29+14+20=63$ MW of beam power loss. In case of Twin Linac we have twice the
recovery beam power loss of $10$ MW, and twice $50\%$ of $100$ MW, summing up
to $120$ MW.

###  III.3 Power budget

Summing up the power budget of primary electron beam to that of FEL driving
electron beam we get a total beam power loss of $202$ MW for the FCC-tL case,
and $200$ MW for the Twin Linac case. An expected efficiency beam-to-plug not
smaller than $20\%$, actually in the range $20-40\%$, would set the AC power
bill in the $250$ MW - $1$ GW range.

One should note that the power transferred from the primary colliding
electron/photon beams into the secondary beams of muon pairs, photons, and
electron-positron pairs is quite negligible compared to the stored power into
the colliding beams at the collision point. As a matter of fact, in case of
Twin Linac scenario, the muon pair beams ($10^{11}$ per second at an average
energy of 50 GeV) are taking out only $1$ kW of beam power. The power taken by
the back-scattered Compton gammas ($8\,10^{12}$ per second at an average
energy of $200$ GeV) is $250$ kW, and the power taken by the electron-positron
pair produced ($10^{16}$ per second at an average energy of $1$ GeV) is about
$1.6$ MW. All these are quite negligible with respect to the total power loss
quoted in Table 1 ($200$ MW), that is dominated by ERL efficiency and FEL
photons (and synchrotron emission in case of FCC-tL). Furthermore, in case of
Twin Linac scenario, the real estate on top of the underground whole Linac
bunker, that extends over about $6$ km2 (nearly $30$ km length times $200$ m
wide), if covered by solar panels would make the whole system almost self-
sustained by about $900$ MW AC power generated by sun light.

##  IV Simulation results for muon beams

The significant values concerning the emitted particles for the four analysed
schemes are reported in Table 2. The best option is represented by the Twin
Linac scheme that produces $2\times 5.4\,10^{10}$ muon pairs per second.

Table 2: Emitted particles characteristics. ∗ TPP and MPP cross section values for the FCC-twin Linac and Twin Linac schemes are calculated as in Section II while Whizard returns values lower by a factor $\sim 4$. ∗∗ The MPP total cross section in the FCC-ee parasitic scheme is not well-defined since $E_{CM}^{2}/M_{\mu}^{2}=5.4$. | FCC-tL | CLIC par | FCC-ee par | TL
---|---|---|---|---
$E_{CM}$ (MeV) | $346$ | $424$ | $244$ | $346$
$\sigma^{TPP}_{tot}$ (mb) | $19^{*}$ | $19.7$ | $17.7$ | $19^{*}$
$\sigma^{ICS}_{tot}$ ($\mu$b) | $14.7$ | $10$ | $27.9$ | $14.7$
$\sigma^{MPP}_{tot}$ (nb) | $216^{*}$ | $235$ | $183^{**}$ | $216^{*}$
$\mathcal{N}_{e^{\pm}}$ (s-1) | $4.8\,10^{14}$ | $2.9\,10^{11}$ | $4.2\,10^{10}$ | $2\times 4.8\,10^{15}$
$\mathcal{N}_{ICS}$ (s-1) | $3.7\,10^{11}$ | $1.5\,10^{8}$ | $6.6\,10^{7}$ | $2\times 3.7\,10^{12}$
$\mathcal{N}_{\mu^{\pm}}$ (s-1) | $5.4\,10^{9}$ | $3.4\,10^{6}$ | $4.3\,10^{5}$ | $2\times 5.4\,10^{10}$
$\epsilon^{n}_{x_{\mu}}$ (nm rad) | $6.7$ | $4.1$ | $873.9$ | $4.6$
$\mathcal{N}_{\mu^{\pm}}/\epsilon^{n}_{x_{\mu}}$ | $806$ | $0.8$ | $0.0005$ | $2\times 11739$
($10^{15}$ m-1s-1)
Figure 4: MPP features. Top line: FCC-Twin Linac scheme; middle line: CLIC
parasitic scheme; bottom line: FCC-ee parasitic scheme.

The MPP has been simulated by means of the Whizard event generator whizard ;
whizard1 , run in such a way to take into account the incoming beams features.
The muon beam transverse normalized emittance for all the schemes is also
reported in Table 2 and the number of muons per second divided by the
emittance, crucial figure of merit as suggested in frank , has been
calculated. The best calculated value for muon beam normalized transverse
emittance is $4.6$ nm rad: this value compares to the analytical prediction of
Eq. 8 giving $5.9$ nm rad. The emitted muon beams features for the FCC-twin
Linac, the CLIC parasitic and the FCC-ee parasitic schemes are reported in
Fig. 4. The energy, angle and transverse momentum distrubutions are shown. The
Twin Linac option has been analysed in detail. The incoming electron beam, the
FEL photon pulse and the MPP features are displyed in Fig. 5.

Figure 5: Twin Linac scheme: incident $e^{-}$ beam, FEL photon pulse and MPP
features.

The outstanding value of the transverse normalized emittance combined with the
number of muon pair per second returns the value
$\mathcal{N}_{\mu^{\pm}}/\epsilon^{n}_{x_{\mu}}=2\times 1.17\,10^{19}$ m-1s-1.
If we consider the muons around the energy distribution peak of $50$ GeV
corresponding to a $10\%$ rms relative energy spread, the $20\%$ of the
produced muons are selected (with a longituadinal emittance value of $\sim
4.5$ mm). The above mentioned coefficient corresponding to this selection is
$\mathcal{N}_{\mu^{\pm}}/\epsilon^{n}_{x_{\mu}}=2\times 2.34\,10^{18}$ m-1s-1,
comparable with the best option of the Gamma Factory GF1 ; GF2  combined with
LEMMA LEMMA  analysed in Ref. frank . Top line of Fig. 6 shows the energy
spectrum of the ICS photons produced in the collision: as expected in a high
electron recoil regime, the energy is peaked around $200$ GeV. The TPP, which
is the most probable collateral reaction, would in this case involve $2\times
4.8\,10^{15}$ primary electrons, still a small fraction of the total ($2\times
10^{18}$ per second). Moreover, the bottom line of Fig. 6 shows that most of
the primary electrons involved in TTP would face a negligible energy loss
since the $e^{+}/e^{-}$ pairs are mainly generated at very low energy (see
middle line).

Figure 6: Twin Linac. Emitted ICS photon, TPP positron and initial electron
beam after TPP reaction beam features.

A useful characteristic of the Twin Linac layout shown in Fig. 3 is the
possibility to accelerate muons in the same Linacs used to accelerate the two
twin primary electron beams. This is possible thanks to the very small
emittance of muons and the non-interacting muon-electron beams. Beam optics
could rely on RF focusing that would be effective both on muons and on
electrons, as further analysed in a future work. Linac acceleration of muons
would allow to bring them up rapidly to the TeV kinetic energy range requested
by muon collider physics, just in a few passes (each Twin Linac pass is $400$
GeV energy gain) through the Linac twin system, using proper muon
recirculation arcs.

##  V Conclusions

We described a muon source based on electron-photon collisions at ultra-high
luminosity, capable to reach muon fluxes up to a few $10^{11}$ muon pairs per
second at an outstanding normalized transverse emittance of a few nm rad, with
muon beam energies peaked at $50$ GeV. The electron-photon collider is based
on a primary electron Linac with energy in the $100-200$ GeV range and a FEL
Linac driver with FEL photon energy in the $100-300$ keV range. Extremely
large beam power (both electron and photon beams) are requested to achieve the
ultra-high luminosity needed, in excess of $10^{40}$ cm-2s-1. Since the
electron-photon collisions transfer only a very small amount of power from the
primary beams into the secondary beams, an efficient energy recovery must be
implemented in the scenario, so to reduce the amount of beam power loss down
to hundreds MW level, from $100$ GW beam power stored in primary beams at
collision. This is the main challenge of such a muon source, together with
challenging beam collision spot size, in the few ten nanometers range, and
handling of an extremely large FEL photon beam power. Two scenarios have been
analysed in details: one adopting the future FCC tunnel, with foreseen $15$ km
radius of curvature, so to host in a slightly modified race-track geometry the
two Linacs generating the primary electron beam and the FEL driving electron
beam. This scenario is limited in performances by the huge amount of power
lost due to synchrotron radiation emission in the arcs. The second scenario is
based on a twin array of Linacs arranged face-to-face, providing both the
primary electron beam and the FEL driving beam: the performances of this
scenario are well much higher in terms of muon flux, for the same beam power
loss. Nevertheless, in this case it is of paramount importance to find a
strategy of beam energy recovery with counterpropagating electron beams in the
same Linac. This is the original ERL layout conceived by M. Tigner long time
ago (see Ref. tigner ), subject to R&D on twin super conducting cavities
accomodating counterpropagating beams, but not yet demonstrated marix ; uh .
Further studies on the feasibility of these scenarios are necessary to assess
the achievable luminosity of a muon collider based on this muon source,
depending on the kind of accumulator scheme to combine with such muon source
scenario. The promise is to achieve the requested luminosity using much lower
muon beam currents with respect to other schemes, thanks to the very low
emittance, which would alleviate significantly the issue of muon beam-induced
background. Additional potentialities are also to be studied and carefully
analysed, e.g. the simultaneous acceleration of muon beams and primary
electron beams in the main Linacs, possible polarization of muon beams if
polarized primary electron beams are used (FEL photons are naturaly
polarized), parasitic use of intense GeV-class positron beams (up to $10^{16}$
pairs per second) generated in the electron-photon collisions, and of $100$
GeV-class intense (in excess of $10^{12}$ per second) monochromatic photons.

## References

  * (1) Muon Accelerator Program, https://map.fnal.gov
  * (2) MICE collaboration, _Demonstration of cooling by the Muon Ionization Cooling Experiment_ , _Nature_ 578 (2020) pg. 53-59
  * (3) D. Alesini _et al._ , _Positron driven muon source for a muon collider_ , arXiv:1905.05747v2 (2019)
  * (4) L. Serafini, C. Curatolo and V. Petrillo, _Low emittance pion beams generation from bright photons and relativistic protons_ , arXiv:1507.06626v2 (2015)
  * (5) L. Serafini, F. Broggi, and C. Curatolo, _Study of Hadron-Photon Colliders for Secondary Beam Generation_ in _Proc. 8th Int. Particle Accelerator Conf. (IPAC’17)_ , Copenaghen, Danimarca (2017) pg. 2865-2867
  * (6) C. Curatolo, F. Broggi and L. Serafini, _Phase space analysis of secondary beams generated in hadron-photon collisions_ , _Nucl. Instrum. Methods Phys. Res., Sect. A_ 865 (2017) pg. 128-132
  * (7) Compact Linear Collider, https://clic.cern/
  * (8) Future Circular Collider, https://fcc-cdr.web.cern.ch/
  * (9) J. W. Motz, H. A. Olsen and H. W. Koch, _Pair Production by Photons_ , _Rev. Mod. Phys._ 41 (1969) pg. 632
  * (10) C. Curatolo _et al._ , _Analytical description of photon beam phase spaces in inverse Compton scattering sources_ , _Phys. Rev. Accel. Beams_ 20 (2017) 080701
  * (11) D. Angal-Kalinin _et al._ , _PERLE. Powerful energy recovery linac for experiments. Conceptual design report_ , _Journal of Physics G_ 45 (2018) 065003
  * (12) M. Klein and A. Stocchi, _PERLE: A High Power Energy Recovery Facility for Europe A contribution to the Update of the European Strategy on Particle Physics_ , (2018)
  * (13) International Linear Collider, https://linearcollider.org/
  * (14) V. Shiltsev and F. Zimmermann, _Modern and future colliders_ , _Reviews of Modern Physics_ 93 (2021) pg. 015006
  * (15) C. Emma, K. Fang, J. Wu and C. Pellegrini, _High efficiency, multiterawatt x-ray free electron lasers_ , _Phys. Rev. Accel. Beams_ 19 (2016) 020705
  * (16) S. Matsuyama _et al._ , _Nanofocusing of X-ray free-electron laser using wavefront-corrected multilayer focusing mirrors_ , _Sci Rep_ 8 (2018) 17440
  * (17) Whizard, https://whizard.hepforge.org/
  * (18) W. Kilian, T. Ohl, J. Reuter, _WHIZARD: Simulating Multi-Particle Processes at LHC and ILC_ , _Eur.Phys.J._ C71 (2011) pg. 1742
  * (19) F. Zimmermann, _LHC/FCC-based muon colliders_ , _J. Phys.: Conf. Ser._ 1067 (2018) 022017
  * (20) M. W. Krasny _et al._ , _The CERN Gamma Factory initiative: an ultra-high intensity gamma source_ in _Proc. 9th Int. Particle Accelerator Conf. (IPAC’18)_ , Vancouver, BC, Canada (2018)
  * (21) W. Placzek _et al._ , _The Gamma Factory Project at CERN: a New Generation of Research Tools Made of Light_ , _Acta Phys. Polon. Supp. 13_ (2020) pg. 645-652
  * (22) M. Tigner, _A possible apparatus for electron clashing-beam experiments_ , _Il Nuovo Cimento (1955-1965)_ 37 (1965) pg. 1228-1231
  * (23) MariX, http://www.marix.eu/
  * (24) I. Konoplev _et al._ , _UH-FLUX: Compact, Energy Efficient Superconducting Asymmetric Energy Recovery LINAC for Ultra-high Fluxes of X-ray and THz Radiation_ in _Proc. of International Particle Accelerator Conference (IPAC’16)_ , Busan, Korea (2016)


much as possible in a way that is $\max$–$\min$ fair and $T_{i}\leq
T_{j_{0}}$, $\forall i\in\mathcal{J}$. Hence, we have to find
$\\{k_{i}\\}_{i\in\mathcal{J}}$ such that $\\{w_{i}\\}_{i\in\mathcal{J}}$ are
increased by $k_{i}$ each. The optimal way of doing this is by first assigning
$k_{i}=\frac{T_{j_{0}}}{\log_{2}\left(1+\gamma_{g,i}\right)}-w_{i}$, $\forall
i\in\mathcal{J}$ (step 6) and checking that
$\sum\limits_{i\in\mathcal{J}}k_{i}\leq
W_{\mathcal{G}}-\sum\limits_{i=1}^{|\mathcal{U}_{g}|}w_{i}$. In case such an
inequality is not satisfied, then assign
$k_{i}=\frac{W_{\mathcal{G}}-\sum\limits_{i\notin\mathcal{J}}w_{i}}{\log_{2}\left(1+\gamma_{g,i}\right)\sum\limits_{i\in\mathcal{J}}\frac{1}{\log_{2}\left(1+\gamma_{g,i}\right)}}-w_{i}$,
$\forall i\in\mathcal{J}$ (step 8).

Once $k_{i}$ is derived, we assign $w_{i}\leftarrow w_{i}+k_{i}$, $\forall
i\in\mathcal{J}$ (step 10). Now, in case that we have assigned
$k_{i}=\frac{T_{j_{0}}}{\log_{2}\left(1+\gamma_{g,i}\right)}-w_{i}$, $\forall
i\in\mathcal{J}$, we reassign set $\mathcal{J}$ as
$\mathcal{J}\leftarrow\mathcal{J}\cup\\{j_{0}\\}$ (step 12), and start all
over. Otherwise, the algorithm stops and the optimal solution is found.

Algorithm 4 Optimal solution of CP (43) (unbounded $\tau_{g}$).

1: Initialize: $w_{i}\\!\leftarrow\\!W_{\mathcal{G}}^{\min}$,
$T_{i}\\!\leftarrow\\!w_{i}\log_{2}\\!\left(1\\!+\\!\gamma_{g,i}\\!\right)$,
$\forall 1\\!\leq\\!i\\!\leq\\!|\mathcal{U}_{g}|$.

2: Assume, WLOG, $T_{i-1}\leq T_{i}$, $\forall 2\leq i\leq|\mathcal{U}_{g}|$.

3:
$\mathcal{J}\leftarrow\left\\{i\in\\{1,\ldots,|\mathcal{U}_{g}|\\}\,\mid\,T_{i}=\min\limits_{1\leq
j\leq|\mathcal{U}_{g}|}T_{j}\right\\}$.

4: while $\sum\limits_{i=1}^{|\mathcal{U}_{g}|}w_{i}<W_{\mathcal{G}}$ and
$|\mathcal{J}|\neq|\mathcal{U}_{g}|$ do

5: $j_{0}\leftarrow\arg\min\\{T_{i}\,\mid\,i\notin\mathcal{J}\\}$ and
$K_{b}\leftarrow 1$.

6:
$k_{i}\leftarrow\frac{T_{j_{0}}}{\log_{2}\left(1+\gamma_{g,i}\right)}-w_{i}$,
$\forall i\in\mathcal{J}$.

7: if
$\sum\limits_{i\in\mathcal{J}}k_{i}>W_{\mathcal{G}}-\sum\limits_{i=1}^{|\mathcal{U}_{g}|}w_{i}$
then

8:
$\\!\\!k_{i}\\!\leftarrow\\!\frac{W_{\mathcal{G}}-\sum\limits_{i\notin\mathcal{J}}w_{i}}{\log_{2}\\!\left(1+\gamma_{g,i}\right)\\!\sum\limits_{i\\!\in\\!\mathcal{J}}\\!\frac{1}{\log_{2}\left(1+\gamma_{g,i}\right)}}\\!-\\!w_{i}$,
$\forall i\\!\in\\!\mathcal{J}$ and $K_{b}\\!\leftarrow\\!0$.

9: end if

10: $w_{i}\leftarrow w_{i}\\!+\\!k_{i}$, $\forall i\\!\in\\!\mathcal{J}$ and
$T_{i}\leftarrow w_{i}\log_{2}\\!\left(1\\!+\\!\gamma_{g,i}\right)$, $\forall
i\\!\in\\!\mathcal{J}$.

11: if $K_{b}=1$ then

12: $\\!\\!\mathcal{J}\leftarrow\mathcal{J}\cup\\{j_{0}\\}$.

13: end if

14: end while

15: $T\leftarrow\min\limits_{1\leq i\leq|\mathcal{U}_{g}|}T_{i}$.

In Note 1 we prove that the assigned $\\{k_{i}\\}_{i\in\mathcal{J}}$ at each
iteration of Algorithm 4 provides the optimal $\max$–$\min$ fair distribution
of resources.

######  Note 1.

Given a distribution of resources $\\{w_{i}\\}_{i=1}^{|\mathcal{U}_{g}|}$ and
users throughput rates $\\{T_{i}\\}_{i=1}^{|\mathcal{U}_{g}|}$ such that
$T_{i}=w_{i}\log_{2}\left(1+\gamma_{g,i}\right)$, $\forall 1\leq
i\leq|\mathcal{U}_{g}|$, we define the set $\mathcal{J}$ as in Eq. (56).

Hence, we have that
$w_{i}\log_{2}\\!\left(1\\!+\\!\gamma_{g,i}\right)=w_{k}\log_{2}\\!\left(1\\!+\\!\gamma_{g,k}\right)$,
$\forall i,k\in\mathcal{J}$.

Now, given $j_{0}=\arg\min\\{T_{i}\,\mid\,i\notin\mathcal{J}\\}$, we want to
increase $\\{w_{i}\\}_{i\in\mathcal{J}}$ as much as possible by $k_{i}$ each
in a $\max$–$\min$ fair way so that
$(w_{i}+k_{i})\log_{2}\left(1+\gamma_{g,i}\right)\leq T_{j_{0}}$, $\forall
i\in\mathcal{J}$. Hence, we must solve the following convex program:

$\begin{cases}\max\min\limits_{i\in\mathcal{J}}\left(w_{i}+k_{i}\right)\log_{2}\left(1+\gamma_{g,i}\right);&\\\
\text{s.t.:}\\\ (w_{i}+k_{i})\log_{2}\left(1+\gamma_{g,i}\right)\leq
T_{j_{0}},&\\!\\!\forall i\\!\in\\!\mathcal{J};\\\
\sum\limits_{i=1}^{|\mathcal{U}_{g}|}w_{i}+\sum\limits_{i\in\mathcal{J}}k_{i}\leq
W_{\mathcal{G}};&\end{cases}$ (57)

Since CP (57) has a $\max$–$\min$ utility function, we reformulate this CP
into an equivalent CP to which we can directly apply KKT conditions:

$\begin{cases}\max L;&\\\ \text{s.t.:}\\\
(w_{i}+k_{i})\log_{2}\left(1+\gamma_{g,i}\right)\geq L,&\forall
i\\!\in\\!\mathcal{J};\\\
k_{i}\leq\frac{T_{j_{0}}}{\log_{2}\left(1+\gamma_{g,i}\right)}-w_{i},&\forall
i\\!\in\\!\mathcal{J};\\\ \sum\limits_{i\in\mathcal{J}}k_{i}\leq
W_{\mathcal{G}}-\sum\limits_{i=1}^{|\mathcal{U}_{g}|}w_{i};&\end{cases}$ (58)

Please, note that the first constraint of CP (58) can be rearranged as
$k_{i}\geq\frac{L}{\log_{2}\left(1+\gamma_{g,i}\right)}-w_{i}$, $\forall
i\in\mathcal{J}$. Also, the second and third constraints of CP (58) are
equivalent to the constraints of the $\max$–$\min$ CP (57).

Now, we make use of KKT conditions to solve CP (58). The decision variables
are gathered in vector
$x=\left[\\{k_{i}\\}_{i\in\mathcal{J}},L\right]\in\mathds{R}^{|\mathcal{J}|+1}$.
We define the following KKT functions:

$\displaystyle f(x)=$ $\displaystyle-L;$ (59) $\displaystyle g_{i}(x)=$
$\displaystyle\frac{L}{\log_{2}\left(1+\gamma_{g,i}\right)}-k_{i}-w_{i},$
$\displaystyle\forall i\in\mathcal{J};$ (60) $\displaystyle
g_{|\mathcal{J}|+i}(x)=$ $\displaystyle
k_{i}-\frac{T_{j_{0}}}{\log_{2}\left(1+\gamma_{g,i}\right)}-w_{i},$
$\displaystyle\forall i\in\mathcal{J};$ (61) $\displaystyle
g_{2|\mathcal{J}|+1}(x)=$
$\displaystyle\sum\limits_{i\in\mathcal{J}}k_{i}-W_{\mathcal{G}}-\sum\limits_{i=1}^{|\mathcal{U}_{g}|}w_{i},$
$\displaystyle\forall i\in\mathcal{J}.$ (62)

Hence, the KKT gradients are the following:

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\nabla\\!f(x)=$
$\displaystyle[0,\ldots,0,1];$ (63)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\nabla\\!g_{i}(x)=$
$\displaystyle\left[\\!0,...,0,-1\\!|_{i},0,...,0,\frac{1}{\log_{2}\\!\left(1\\!+\\!\gamma_{g,i}\\!\right)}\\!\right]\\!\\!,$
$\displaystyle\\!\\!\forall i\\!\in\\!\mathcal{J}\\!;$ (64)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\nabla\\!g_{|\mathcal{J}|\\!+\\!i}(x)=$
$\displaystyle\left[0,\ldots,0,1\\!|_{i},0,\ldots,0\right],$
$\displaystyle\\!\\!\forall i\\!\in\\!\mathcal{J}\\!;$ (65)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\nabla\\!g_{2|\mathcal{J}|\\!+\\!1}(x)=$
$\displaystyle\left[1,\ldots,1,0\right].$ (66)

We derive the following KKT conditions:

$\displaystyle-\mu_{i}+\mu_{|\mathcal{J}|+i}+\mu_{2|\mathcal{J}|+1}=0,$
$\displaystyle\forall i\in\mathcal{J};$ (67)
$\displaystyle-1+\sum\limits_{i=1}^{|\mathcal{J}|}\frac{\mu_{i}}{\log_{2}\left(1+\gamma_{g,i}\right)}=0;$
(68)
$\displaystyle\mu_{i}\cdot\left(\frac{L}{\log_{2}\left(1+\gamma_{g,i}\right)}-k_{i}-w_{i}\right)=0,$
$\displaystyle\forall i\in\mathcal{J};$ (69)
$\displaystyle\mu_{|\mathcal{J}|+i}\cdot\left(k_{i}-\frac{T_{j_{0}}}{\log_{2}\left(1+\gamma_{g,i}\right)}+w_{i}\right)=0,$
$\displaystyle\forall i\in\mathcal{J};$ (70)
$\displaystyle\mu_{2|\mathcal{J}|+1}\cdot\left(\sum\limits_{i\in\mathcal{J}}k_{i}-W_{\mathcal{G}}-\sum\limits_{i=1}^{|\mathcal{U}_{g}|}w_{i}\right)=0.$
(71)

First, we note that if assigning
$k_{i}=\frac{T_{j_{0}}}{\log_{2}\left(1+\gamma_{g,i}\right)}-w_{i}$, $\forall
i\in\mathcal{J}$ accomplishes that $\sum\limits_{i\in\mathcal{J}}k_{i}\leq
W_{\mathcal{G}}-\sum\limits_{i=1}^{|\mathcal{U}_{g}|}w_{i}$, we have the
optimal solution, as each $k_{i}$ receives the maximum possible value and
constraints are satisfied.

Hence, we can assume that $\exists
i_{0}\in\mathcal{J}\,\mid\,k_{i_{0}}<\frac{T_{j_{0}}}{\log_{2}\left(1+\gamma_{g,i_{0}}\right)}-w_{i_{0}}$
and hence, according to Eq. (70), $\mu_{|\mathcal{J}|+i_{0}}=0$.

The optimal $\\{k_{i}\\}_{i\in\mathcal{J}}$ and $L$ that solve the KKT
conditions are:

$\displaystyle k_{i}=$
$\displaystyle\frac{W_{\mathcal{G}}-\sum\limits_{i\notin\mathcal{J}}w_{i}}{\log_{2}\left(1+\gamma_{g,i}\right)\sum\limits_{i\in\mathcal{J}}\frac{1}{\log_{2}\left(1+\gamma_{g,i}\right)}}-w_{i},$
$\displaystyle\forall i\in\mathcal{J};$ (72) $\displaystyle L=$
$\displaystyle\frac{W_{\mathcal{G}}-\sum\limits_{i\notin\mathcal{J}}w_{i}}{\sum\limits_{i\in\mathcal{J}}\frac{1}{\log_{2}\left(1+\gamma_{g,i}\right)}}.$
(73)

It is simple to check that
$\sum\limits_{i\in\mathcal{J}}k_{i}=W_{\mathcal{G}}-\sum\limits_{i=1}^{|\mathcal{U}_{g}|}w_{i}$,
so that Eq. (71) is satisfied.

Now, we assign
$\mu_{i}=\frac{w_{i}}{\sum\limits_{i\in\mathcal{J}}w_{i}}\log_{2}\left(1+\gamma_{g,i}\right)$,
$\forall i\in\mathcal{J}$, and hence, since $\mu_{|\mathcal{J}|+i_{0}}=0$, and
according to Eq. (67), we have that
$\mu_{2|\mathcal{J}|+1}=\mu_{i_{0}}=\frac{w_{i_{0}}}{\sum\limits_{i\in\mathcal{J}}w_{i}}\log_{2}\left(1+\gamma_{g,i_{0}}\right)$.
Hence, according also to Eq. (67),
$\mu_{|\mathcal{J}|+i}=\mu_{i}-\mu_{2|\mathcal{J}|+1}=\frac{w_{i}}{\sum\limits_{i\in\mathcal{J}}w_{i}}\log_{2}\left(1+\gamma_{g,i}\right)-\frac{w_{i_{0}}}{\sum\limits_{i\in\mathcal{J}}w_{i}}\log_{2}\left(1+\gamma_{g,i_{0}}\right)$,
$\forall i\in\mathcal{J}$. Since $i,i_{0}\in\mathcal{J}$, we have that, by
definition on $\mathcal{J}$, $\mu_{|\mathcal{J}|+i}=0$, $\forall
i\in\mathcal{J}$.

As a result, since with such a solution all KKT conditions, are satisfied, we
have the optimal solution to CP (58).

With limited backbone capacity $\boldsymbol{\tau_{g}}$. Now, we assume that
$\tau_{g}$ is bounded. In order to find the optimal solution under this
assumption, we first solve the problem assuming unbounded $\tau_{g}$, as
detailed above. If such an output provides a feasible solution, i.e.,
$\sum\limits_{u\in\mathcal{U}_{g}}T_{u}\leq\tau_{g}$, we have the optimal
solution. Hence, we assume that such an output provides an unfeasible
solution, i.e., $\sum\limits_{u\in\mathcal{U}_{g}}T_{u}>\tau_{g}$.

Let $\\{w_{u}\\}_{u\in\mathcal{U}_{g}}$, $\\{T_{u}\\}_{u\in\mathcal{U}_{g}}$
be the unfeasible solution provided by the $\max$–$\min$ optimization with
unbounded $\tau_{g}$, and let $T=\min\limits_{u\in\mathcal{U}_{g}}\\{T_{u}\\}$
be the minimum achieved throughput by users so far (note that $T$ is the value
of the utility). Due to the $\max$–$\min$ fairness nature, every user $u$ such
that $T_{u}>T$ disposes of the minimum amount of resources,
$W_{\mathcal{G}}^{\min}$ and
$T_{u}=W_{\mathcal{G}}^{\min}\log_{2}\left(1+\gamma_{g,u}\right)$, $\forall
u\in\mathcal{U}_{g}\,\mid\,T_{u}>T$ (otherwise, if such users disposed of more
than $W_{\mathcal{G}}^{\min}$ resources, such exceeded resources could be
reallocated to those users with minimum throughput to increase the utility
function, which is not possible from the output $\max$–$\min$ fairness
optimization with unbounded $\tau_{g}$).

Hence, as no resources can be removed from any user $u$ such that $T_{u}>T$,
and $\sum\limits_{u\in\mathcal{U}_{g}}T_{u}>\tau_{g}$, we can reduce the rate
of these users to, for instance, $T_{u}=T$ (hence, now $T_{u}=T$, $\forall
u\in\mathcal{U}_{g}$). If now
$\sum\limits_{u\in\mathcal{U}_{g}}T_{u}\leq\tau_{g}$, we have found the
optimal solution. Otherwise, we need to reduce more individual throughput
rates. In order to be $\max$–$\min$ fair, we have to reduce every individual
rate the same amount until the $\tau_{g}$–constraint is satisfied. Hence, we
define $R=\sum\limits_{u\in\mathcal{U}_{g}}T_{u}-\tau_{g}>0$ and decrease
every individual rate by $\frac{R}{|\mathcal{U}_{g}|}$, i.e.,
$T_{u}=T-\frac{R}{|\mathcal{U}_{g}|}$, $\forall u\in\mathcal{U}_{g}$. Hence,
the $\max$–$\min$ fairness optimization of CP (43) is solved also under the
assumption of bounded $\tau_{g}$.

###  D.2 With One Drone

In this scenario, we consider that $g\\!B\\!S$ $g$ serves its users
$\mathcal{U}_{g}$ and also one $a\\!B\\!S$ $a$, i.e., $\exists
a\in\mathcal{A}\,\mid\,\mathcal{A}_{g}=\\{a\\}$. Hence, CP of Fig. 6
simplifies since backhaul resources do not need to be split over multiple
$a\\!B\\!S$s, but all of them are assigned to only one $a\\!B\\!S$, i.e.,
$w^{a}=W_{\mathcal{B}}$. Hence, we solve the following CP:

${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\begin{cases}\,\,\,\qquad\qquad\qquad\qquad\qquad\mathclap{\max\limits_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}w_{u},}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}T_{u}}}U_{\mathrm{cvx},g}^{\alpha}=\begin{cases}\sum\limits_{u\in\mathcal{U}_{g}\cup\,\mathcal{U}_{a}}\\!\left(T_{u}\right)^{1-\alpha}\cdot\frac{1}{1-\alpha},&\alpha\neq
1;\\\
\sum\limits_{u\in\mathcal{U}_{g}\cup\,\mathcal{U}_{a}}\\!\log\left(T_{u}\right),&\alpha=1;\end{cases}}\\\
\text{s.t.:}\\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}T^{a}\\!\leq\\!W_{\mathcal{B}}\log_{2}\\!\left(1\\!+\\!\gamma_{g,a}^{\mathcal{B}}\right)\\!;}&\\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}w_{u}\geq
W^{\min}_{\mathcal{G}},}&{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\forall
u\\!\in\\!\mathcal{U}_{g};}\\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sum\limits_{u\in\mathcal{U}_{g}}\\!w_{u}=W_{\mathcal{G}};}&\\!\\!\\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}w_{u}\geq
W^{\min}_{\mathcal{A}},}&{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\forall
u\\!\in\\!\mathcal{U}_{a};}\\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sum\limits_{u\in\mathcal{U}_{a}}\\!w_{u}=W_{\\!\mathcal{A}};}&\\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}T_{u}\\!\leq\\!w_{u}\log_{2}\\!\left(1\\!+\\!\gamma_{g,u}\right)\\!,}&{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\forall
u\in\mathcal{U}_{g};}\\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}T_{u}\\!\leq\\!w_{u}\log_{2}\\!\left(1\\!+\\!\gamma_{a,u}\right)\\!,}&{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\forall
u\in\mathcal{U}_{a};}\\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sum\limits_{u\in\mathcal{U}_{a}}\\!T_{u}\\!\leq
T^{a};}&\\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sum\limits_{u\in\mathcal{U}_{g}}\\!T_{u}\\!+\\!T^{a}\leq\tau_{g}.}&\\!\\!\end{cases}}$
(74)

In order to solve CP (74), we do not derive KKT conditions, as we will make
use of the analysis conducted in Subection D.1.

With unbounded backbone capacity $\boldsymbol{\tau_{g}}$. First, we assume
that $\tau_{g}$ is unbounded, and hence users throughput rates are only
limited by their Shannon capacity, according to the allocated bandwidth
resources. As $g\\!B\\!S$ $g$ disposes of two independent baskets of
resources, $W_{\mathcal{G}}$ for $g\\!B\\!S$-served users, and
$W_{\mathcal{B}}$ for the backhaul-served $a\\!B\\!S$, the distribution of
$g\\!B\\!S$-served users resources disregards from $a\\!B\\!S$-served users
resources (since $\tau_{g}$ is unbounded). Hence, the backhaul rate for
$a\\!B\\!S$ $a$, $T^{a}$, is limited only by the Shannon capacity, so that we
can assume that
$T^{a}=W_{\mathcal{B}}\log_{2}\left(1+\gamma_{g,a}^{\mathcal{B}}\right)$.

Hence, we solve optimal resource allocation for $g\\!B\\!S$-served users
assuming unbounded $\tau_{g}$, and solve also optimal resource allocation for
$a\\!B\\!S$-served users assuming that their aggregated throughput is limited
by the backhaul capacity $T^{a}$, i.e.,
$\sum\limits_{u\in\mathcal{U}_{a}}T_{u}\leq T^{a}$. Both cases can be solved
as detailed in Subsection D.1. As a result, CP (74) is solved under the
assumption of unbounded $\tau_{g}$, for any value of $\alpha$, including
$\alpha\rightarrow\infty$.

The solution provided is optimal although not necessarily unique, as the
backhaul capacity $T^{a}$ might be higher than the aggregated served
throughputs. Hence, to ease the understanding of upcoming sections, we assume
that the provided optimal solution satisfies that
$T^{a}=\sum\limits_{u\in\mathcal{U}_{a}}T_{u}$, which remains feasible and
optimal.

With limited backbone capacity $\boldsymbol{\tau_{g}}$. Now, we do not assume
that $\tau_{g}$ is unbounded. In order to solve this case, we first search the
optimal solution assuming that $\tau_{g}$ is indeed unbounded and checking if
$\sum\limits_{u\in\mathcal{U}_{g}}T_{u}+T^{a}=\sum\limits_{u\in\mathcal{U}_{g}}T_{u}+\sum\limits_{u\in\mathcal{U}_{a}}T_{u}\leq\tau_{g}$.
In case the inequality is satisfied, we have found the optimal solution.
Otherwise, users throughputs must be decreased in order to satisfy the
$\tau_{g}$–constraint.

Hence, according to what proven and described in Subsection D.1, we can
decrease the convenient users throughputs $T_{u}$ until we get to satisfy that
$\sum\limits_{u\in\mathcal{U}_{g}}T_{u}+T^{a}=\sum\limits_{u\in\mathcal{U}_{g}}T_{u}+\sum\limits_{u\in\mathcal{U}_{a}}T_{u}=\tau_{g}$
and hence provide the optimal solution.

###  D.3 Generic Case

In this scenario, we consider the generic case, formulated in Fig. 6. As in
the previous cases, we distinguish between unbounded and bounded $\tau_{g}$
assumptions. In Algorithm 5 we show how to find the optimal solution, as
detailed in what follows.

With unbounded backbone capacity $\boldsymbol{\tau_{g}}$. Here, we assume that
$\tau_{g}$ is unbounded, so that distribution of user resources at
$g\\!B\\!S$s is not affected by the presence of $a\\!B\\!S$s and their
associated users. However, now we cannot know a priori how many resources each
$a\\!B\\!S$ will get, as they share a common bandwidth of
$W_{\\!\mathcal{B}}$. Each $a\\!B\\!S$ must receive $w^{a}$ resources
according to the optimization output. We find the optimal solution as detailed
as follows.

First, the optimal resource allocation for $g\\!B\\!S$-served users is as
previously detailed in Section D.1. Then, for each $a\\!B\\!S$
$a\in\mathcal{A}_{g}$ we solve also optimal resource allocation assuming
unbounded backhaul capacity $T^{a}$, i.e., $a\\!B\\!S$-served users throughput
rates are limited only by their Shannon capacity, according to the bandwidth
allocated, using the scheme detailed in Section D.1. For convenience, we now
denote by $T_{a,u}$ the throughput of access link $(a,u)$, $\forall
a\in\mathcal{A}_{g}$ $\forall u\in\mathcal{U}_{a}$. For each
$a\in\mathcal{A}_{g}$, we need a backhaul capacity of
$T^{a}=\sum\limits_{u\in\mathcal{U}_{a}}T_{a,u}$. Hence, $a\\!B\\!S$ $a$ needs
$w^{a}=\max\left(W_{\\!\mathcal{B}}^{\min},\sum\limits_{u\in\mathcal{U}_{a}}\frac{T_{a,u}}{\log_{2}\left(1+\gamma_{g,a}^{\mathcal{B}}\right)}\right)$.

Now, in case $\sum\limits_{a\in\mathcal{A}_{g}}w^{a}\leq W_{\\!\mathcal{B}}$,
we choose $a_{0}\in\mathcal{A}_{g}$ arbitrarily and add resources to
$w^{a_{0}}$ so that we get to
$\sum\limits_{a\in\mathcal{A}_{g}}w^{a}=W_{\\!\mathcal{B}}$ (i.e.,
$w^{a_{0}}=w^{a_{0}}+W_{\\!\mathcal{B}}-\sum\limits_{a\in\mathcal{A}_{g}}w^{a}$).
Hence, in this case the optimization is solved.

Conversely, in case
$\sum\limits_{a\in\mathcal{A}_{g}}w^{a}>W_{\\!\mathcal{B}}$, we are assigning
to backhaul $a\\!B\\!S$s more resources than what available at the
$g\\!B\\!S$. Hence, some resources should be removed. In this case, we need to
build the optimal solution from scratch. First, we assign to each $a\\!B\\!S$
$a\in\mathcal{A}_{g}$ the minimum bandwidth:
$w^{a}=W_{\\!\mathcal{B}}^{\min}$, $\forall a\in\mathcal{A}_{g}$ (step 1 in
Algorithm 5). Second, we assign to each $a\\!B\\!S$ $a\in\mathcal{A}_{g}$ the
highest achievable throughput:
$T^{a}=w^{a}\log_{2}\left(1+\gamma_{g,a}^{\mathcal{B}}\right)$, $\forall
a\in\mathcal{A}_{g}$ (step 2). Third, we solve optimal resource allocation at
each $a\in\mathcal{A}_{g}$ and also at $g\\!B\\!S$ $g$. Hence, we dispose of
$\\{T_{a,u}\\}_{u\in\mathcal{U}_{a}}$, $\forall a\in\mathcal{A}_{g}$, and of
$\\{T_{u}\\}_{u\in\mathcal{U}_{g}}$ (step 3). Please note that $\forall
a\in\mathcal{A}_{g}$, $\forall u\in\mathcal{U}_{a}$ it might happen that
$T_{a,u}<w_{u}\log_{2}\left(1+\gamma_{a,u}\right)$ because of the backhaul
limitation $T^{a}$. Hence, if $T^{a}$ increases, such $T_{a,u}$’s can also
increase. The contribution of $\\{T_{a,u}\\}_{u\in\mathcal{U}_{a}}$, $\forall
a\in\mathcal{A}_{g}$ to utility $\mathcal{U}_{\mathrm{cvx},g}^{\alpha}$ is
either $\left(T_{a,u}\right)^{1-\alpha}$ or $\log\left(T_{a,u}\right)$. Hence,
the distribution of throughput contribution follows an increasing and concave
function. Since $\xi(x)=x^{1-\alpha}$ and $\xi(x)=\log\left(x\right)$ are both
increasing and concave functions, the way of increasing most utility
$\mathcal{U}_{\mathrm{cvx},g}^{\alpha}$ is by raising the lowest values of
$\\{T_{a,u}\\}_{u\in\mathcal{U}_{a}}$, $\forall a\in\mathcal{A}_{g}$ as much
as possible as long as constraints are not violated. Also in the $\max$–$\min$
case: the way of increasing as much as possible the utility is by equally
raising the lowest values of $\\{T_{a,u}\\}_{u\in\mathcal{U}_{a}}$, $\forall
a\in\mathcal{A}_{g}$, as long as constraints are not violated, due to the
well-known $\max$–$\min$ fairness nature. Let

$\displaystyle
T_{m}\\!=\\!\min\\!\left\\{T_{a,u}\\!\\!\,\mid\,\\!T_{a,u}\\!<\\!w_{u}\log_{2}\\!\left(1\\!+\\!\gamma_{a,u}\right),a\\!\in\\!\mathcal{A}_{g},u\\!\in\\!\mathcal{U}_{a}\\!\right\\}$
(75)

(step 6) be the minimum throughput rate that has not reached the corresponding
to the Shannon capacity (note that if such a $T_{m}$ does not exist, we are
done). Let

$\displaystyle\mathcal{L}_{m}=\left\\{(a,u)\in\mathcal{A}_{g}\times\bigcup\limits_{a\in\mathcal{A}_{g}}\mathcal{U}_{a}\,\mid\,T_{a,u}=T_{m}\right\\}$
(76)

be the set of those $a\\!B\\!S$–UE links such that the link rate is the same
as the minimum $T_{m}$ (step 7). Let

$\displaystyle
T_{M}=\min\left\\{T_{a,u}\,\mid\,T_{a,u}>T_{m},a\in\mathcal{A}_{g},u\in\mathcal{U}_{a}\right\\}$
(77)

be the minimum throughput rate among those rates that are not as the minimum
$T_{m}$ (step 8). Let

$\displaystyle
T_{M_{2}}=\min\left(T_{M},\min\limits_{(a,u)\in\mathcal{L}_{m}}w_{u}\log_{2}\left(1+\gamma_{a,u}\right)\right)$
(78)

(step 9). The goal now is to increase
$\\{T_{a,u}\\}_{(a,u)\in\mathcal{L}_{m}}$ as much as possible not exceeding
$T_{M_{2}}$, as long as those involved $a\in\mathcal{A}_{g}$ can request more
resources to increase $T^{a}$. Let $\beta\in[0,1]$ be an undetermined
parameter, $\\{T_{a,u}\\}_{(a,u)\in\mathcal{L}_{m}}$ will be increased by
$\beta\\!\cdot\\!(T_{M_{2}}-T_{m})$, i.e., at most, by $T_{M_{2}}-T_{m}$ (step
14). Parameter $\beta$ will be defined later. Let

$\displaystyle\overline{\mathcal{U}}_{a}=\\{u\in\mathcal{U}_{a}\,\mid\,(a,u)\in\mathcal{L}_{m}\\},\forall
a\in\mathcal{A}_{g}$ (79)

(step 10). Now, we set
$T_{a,u}^{{}^{\prime}}=T_{a,u}+\beta\\!\cdot\\!(T_{M_{2}}-T_{m})$,
$\forall(a,u)\in\mathcal{L}_{m}$ to increase the involved throughput rates.
Hence, we set $\forall a\in\mathcal{A}_{g}$

$\displaystyle T^{a}$ $\displaystyle\\!\\!\\!\\!=\\!\\!\\!\\!$
$\displaystyle\sum\limits_{u\notin\overline{\mathcal{U}}_{a}}T_{a,u}+\sum\limits_{u\in\overline{\mathcal{U}}_{a}}T_{a,u}^{{}^{\prime}}=$
(80) $\displaystyle\\!\\!\\!\\!=\\!\\!\\!\\!$
$\displaystyle\sum\limits_{u\notin\overline{\mathcal{U}}_{a}}T_{a,u}+\sum\limits_{u\in\overline{\mathcal{U}}_{a}}\left(T_{a,u}+\beta\\!\cdot\\!(T_{M_{2}}-T_{m})\right)=$
$\displaystyle\\!\\!\\!\\!=\\!\\!\\!\\!$
$\displaystyle\sum\limits_{u\notin\overline{\mathcal{U}}_{a}}T_{a,u}+\sum\limits_{u\in\overline{\mathcal{U}}_{a}}T_{a,u}+|\overline{\mathcal{U}}_{a}|\beta\\!\cdot\\!(T_{M_{2}}-T_{m}).$

Hence, we set $\forall a\in\mathcal{A}_{g}$

$\displaystyle w^{a}_{new}$ $\displaystyle\\!\\!\\!\\!=\\!\\!\\!\\!$
$\displaystyle\frac{T_{a}}{\log_{2}\\!\left(1\\!+\\!\gamma_{g,a}^{\mathcal{B}}\\!\right)}=\frac{\sum\limits_{u\in\mathcal{U}_{a}}\\!\\!T_{a,u}\\!+\\!|\overline{\mathcal{U}}_{a}|\beta\\!\cdot\\!(\\!T_{M_{2}}\\!\\!-\\!T_{m}\\!)}{\log_{2}\\!\left(1\\!+\\!\gamma_{g,a}^{\mathcal{B}}\\!\right)}=$
(81) $\displaystyle\\!\\!\\!\\!=\\!\\!\\!\\!$ $\displaystyle
w^{a}+\frac{|\overline{\mathcal{U}}_{a}|\beta\\!\cdot\\!(T_{M_{2}}\\!-\\!T_{m})}{\log_{2}\left(1+\gamma_{g,a}^{\mathcal{B}}\right)}$

(step 12). Now, the aggregation of the new backhaul bandwidth allocation has
to be lower than the total bandwidth, i.e.,

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\sum\limits_{a\in\mathcal{A}_{g}}\\!\\!w^{a}_{new}$
$\displaystyle\\!\\!\\!\\!=\\!\\!\\!\\!$
$\displaystyle\\!\sum\limits_{a\in\mathcal{A}_{g}}\\!\\!\left(w^{a}+\frac{|\overline{\mathcal{U}}_{a}|\beta\\!\cdot\\!(T_{M_{2}}-T_{m})}{\log_{2}\left(1+\gamma_{g,a}^{\mathcal{B}}\right)}\right)=$
(82) $\displaystyle\\!\\!\\!\\!=\\!\\!\\!\\!$
$\displaystyle\\!\sum\limits_{a\in\mathcal{A}_{g}}\\!\\!\\!w^{a}\\!+\\!\beta\\!\cdot\\!(T_{M_{2}}\\!-\\!T_{m})\\!\\!\sum\limits_{a\in\mathcal{A}_{g}}\\!\\!\\!\frac{|\overline{\mathcal{U}}_{a}|}{\log_{2}\\!\left(1\\!\\!+\\!\gamma_{g,a}^{\mathcal{B}}\\!\right)}$

has to be lower or equal than $W_{\\!\mathcal{B}}$. Hence, isolating $\beta$
we get that necessarily,

$\displaystyle\beta\leq\frac{W_{\\!\mathcal{B}}-\sum\limits_{a\in\mathcal{A}_{g}}w^{a}}{(T_{M_{2}}-T_{m})\sum\limits_{a\in\mathcal{A}_{g}}\frac{|\overline{\mathcal{U}}_{a}|}{\log_{2}\left(1+\gamma_{g,a}^{\mathcal{B}}\right)}}.$
(83)

Hence, we define $\beta$ as

$\displaystyle\beta=\min\left(1,\frac{W_{\\!\mathcal{B}}-\sum\limits_{a\in\mathcal{A}_{g}}w^{a}}{(T_{M_{2}}-T_{m})\sum\limits_{a\in\mathcal{A}_{g}}\frac{|\overline{\mathcal{U}}_{a}|}{\log_{2}\left(1+\gamma_{g,a}^{\mathcal{B}}\right)}}\right)$
(84)

(step 11). Once the parameter $\beta$ is derived, we assign
$w^{a}=w^{a}_{new}$ and
$T^{a}=w^{a}\log_{2}\left(1+\gamma_{g,a}^{\mathcal{B}}\right)$, $\forall
a\in\mathcal{A}_{g}$ (step 13). In case that $\beta=1$ (step 5), we repeat the
process defining $T_{m}$ again and increasing the corresponding throughput
rates. This yields the optimal solution.

With limited backbone capacity $\boldsymbol{\tau_{g}}$. Now, we do not assume
that $\tau_{g}$ is unbounded. Hence, we now find the optimal solution by
solving the problem as above, assuming that $\tau_{g}$ is indeed unbounded,
and hence progressively decreasing highest individual throughput rates
$T_{u}$, $\forall
u\in\mathcal{U}_{g}\cup\bigcup\limits_{a\in\mathcal{A}_{g}}\mathcal{U}_{a}$
until
$\sum\limits_{u\in\mathcal{U}_{g}}T_{u}+\sum\limits_{a\in\mathcal{A}_{g}}\sum\limits_{u\in\mathcal{U}_{a}}T_{u}=\tau_{g}$,
as done also in the analyzed cases in Subsections D.1 and D.2 (step 18 of
Algorithm 5).

Algorithm 5 Optimal solution to CP of Fig. 6. Generic case

1: $w^{a}\leftarrow W_{\mathcal{B}}^{\min}$, $\forall a\in\mathcal{A}_{g}$.

2: $T^{a}\leftarrow w^{a}\log_{2}\left(1+\gamma_{g,a}^{\mathcal{B}}\right)$,
$\forall a\in\mathcal{A}_{g}$.

3: Derive $\\{T_{u}\\}_{u\in\mathcal{U}_{g}}$ and
$\\{T_{a,u}\\}_{u\in\mathcal{U}_{a}}$, $\forall a\in\mathcal{A}_{g}$ by
solving optimal resource allocation for $g\\!B\\!S$ $g$ assuming unbounded
$\tau_{g}$ and also $\forall a\in\mathcal{A}_{g}$ assuming a backhaul
limitation of $T^{a}$, as detailed in Subsection D.1.

4: $\beta\leftarrow 1$.

5: while $\beta=1$ do

6:
$T_{m}\\!\leftarrow\\!\min\\!\left\\{T_{a,u}\\!\\!\,\mid\,\\!\\!T_{a,u}\\!\\!<\\!w_{u}\log_{2}\\!\left(\\!1\\!\\!+\\!\gamma_{g,a}\\!\right)\\!,a\\!\in\\!\mathcal{A}_{g},u\\!\in\\!\mathcal{U}_{a}\right\\}\\!$.

7:
$\mathcal{L}_{m}\\!\leftarrow\\!\left\\{(a,u)\in\mathcal{A}_{g}\times\bigcup\limits_{a\in\mathcal{A}_{g}}\mathcal{U}_{a}\,\mid\,T_{a,u}=T_{m}\right\\}$.

8:
$T_{M}\\!\leftarrow\\!\min\left\\{T_{a,u}\,\mid\,T_{a,u}>T_{m},a\in\mathcal{A}_{g},u\in\mathcal{U}_{a}\right\\}$.

9:
$T_{M_{2}}\\!\leftarrow\\!\min\left(T_{M},\min\limits_{(a,u)\in\mathcal{L}_{m}}w_{u}\log_{2}\left(1+\gamma_{a,u}\right)\right)$.

10:
$\overline{\mathcal{U}}_{a}\\!\leftarrow\\!\left\\{u\in\mathcal{U}_{a}\,\mid\,(a,u)\in\mathcal{L}_{m}\right\\}$,
$\forall a\in\mathcal{A}_{g}$.

11:
$\beta\\!\leftarrow\\!\min\left(1,\frac{W_{\mathcal{B}}-\sum\limits_{a\in\mathcal{A}_{g}}w^{a}}{\left(T_{M_{2}}-T_{m}\right)\cdot\sum\limits_{a\in\mathcal{A}_{g}}\frac{|\overline{\mathcal{U}}_{a}|}{\log_{2}\left(1+\gamma_{g,a}^{\mathcal{B}}\right)}}\right)$.

12:
$w^{a}\\!\leftarrow\\!w^{a}+\frac{|\overline{\mathcal{U}}_{a}|\cdot\beta\cdot\left(T_{M_{2}}-T_{m}\right)}{\log_{2}\left(1+\gamma_{g,a}^{\mathcal{B}}\right)}$,
$\forall a\in\mathcal{A}_{g}$.

13:
$T^{a}\\!\leftarrow\\!w^{a}\log_{2}\left(1+\gamma_{g,a}^{\mathcal{B}}\right)$,
$\forall a\in\mathcal{A}_{g}$.

14: $T_{a,u}\\!\leftarrow\\!T_{a,u}+\beta\cdot\left(T_{M_{2}}-T_{m}\right)$,
$\forall(a,u)\in\mathcal{L}_{m}$.

15: end while

16: $T^{a}\\!\leftarrow\\!\sum\limits_{u\in\mathcal{U}_{a}}T_{a,u}$, $\forall
a\in\mathcal{A}_{g}$.

17: if
$\sum\limits_{u\in\mathcal{U}_{g}}T_{u}+\sum\limits_{a\in\mathcal{A}_{g}}T^{a}>\tau_{g}$
then

18: Apply Algorithm 2 to
$n\\!=\\!|\mathcal{U}_{g}|+\\!\\!\sum\limits_{a\in\mathcal{A}_{g}}\\!\\!|\mathcal{U}_{a}|$,
get new line
$\\{x_{i}\\}_{i=1}^{n}\\!=\\!\\{T_{u}\\}_{u\in\mathcal{U}_{g}}\\!\cup\\!\bigcup\limits_{a\in\mathcal{A}_{g}}\\!\\!\\{T_{a,u}\\}_{u\in\mathcal{U}_{a}}$,
$X\\!=\\!\tau_{g}$.

19: end if