# Broadband energy squeezing and tunneling based on unidirectional modes

Lujun Hong1, Yazhou Wang2, Yun Shen1,6, Xiaohua Deng 1, Kai Yuan3, Sanshui
Xiao2, and Jie Xu 4,5,7 1 Institute of Space Science and Technology, Nanchang
University, Nanchang 330031, China  2 DTU Fotonik, Department of Photonics
Engineering, Technical University of Denmark, DK-2800 Kgs. Lyngby, Denmark  3
Information Engineering School of Nanchang University, Nanchang University,
Nanchang 330031, China  4 School of Medical Information and Engineering,
Southwest Medical University, Luzhou 646000, China  5 Medicine & Engineering &
Information Fusion and Transformation Key Laboratory of Luzhou City, Luzhou
646000, China  6<EMAIL_ADDRESS>7<EMAIL_ADDRESS>
###### Abstract

Energy squeezing attracts many attentions for its potential applications in
electromagnetic (EM) energy harvesting and optical communication. However, due
to the Fabry-Perot resonance, only the EM waves with discrete frequencies can
be squeezed and, as far as we know, in the previous energy-squeezing devices,
stringent requirements of the materials or the geometrical shape are needed.
We note that the structures filled with epsilon-near-zero (ENZ) mediums as
reported in some works can squeeze and tunnel EM waves at frequencies (e.g.
plasma frequency). However, the group velocity is usually near zero which
means few EM information travel through the structures. In this paper, low-
loss energy squeezing and tunneling (EST) based on unidirectional modes were
demonstrated in YIG-based one-way waveguides at microwave frequencies.
According to our theoretical analysis and the simulations using finite element
method, broadband EST was achieved and the EM EST was observed even for
extremely bended structures. Besides, similar EM EST was achieved in a
realistic three-dimensional remanence-based one-way waveguide as well. The
unidirectional modes-based EST paving the way to ultra-subwavelength EM
focusing, enhanced nonlinear optics, and designing numerous functional devices
in integrated optical circuits such as phase modulator.

##  I Introduction

One-way Electromagnetic (EM) waves propagating in systems in only one
direction are similar to those one-way edge modes founded in the quantum Hall
effectPrange and Girvin (1987). The one-way EM modes was experimentally
observed in 2009 by using a microwave magneto-optical (MO) material, i.e.
yttriun-iron-garnet (YIG)Wang et al. (2009). This kind of one-way EM modes
sustained at the MO interface are named surface magnetoplasmon (SMP) and
utilizing an external magnetic field to break the time-reversal symmetry is a
fundamental way to induce such one-way SMPBrion et al. (1972); Wallis et al.
(1974). Recently, Tsakmakidis and his colleagues proposed an one-way waveguide
consisting of one layer of semiconductor (InSb) and in the waveguide, one
metal layer was set to stop the EM wave, and it was reported that the time-
bandwidth limit in such one-way waveguide can be brokenTsakmakidis et al.
(2017). More recently, we proposed several similar one-way waveguides and
interesting phenomenons such as slowing wave and truly rainbow trapping were
reported as wellXu et al. (2019, 2021).

Squeezing electromagnetic waves at subwavelength or even deep subwavelength
scale is important for near-field microscopyAtkin et al. (2012); Berweger et
al. (2015), semiconductor laserMachida et al. (1987), optical sensingFrascella
et al. (2021) and nonlinear opticsTucker and Millea (1978); Thapliyal et al.
(2014). Epsilon-near-zero (ENZ) materials refers to those materials with a
near-zero relative permittivity and owing to the special physical characters,
they are widely used in researches of metamaterial and metadevices. In the
past two decades, researchers reported several types of devices based on ENZ
materials, in which the EM energy squeezing and tunneling (EST) were
observedSilveirinha and Engheta (2006); Liu et al. (2008, 2009); Jingjing et
al. (2010); Li et al. (2017). Moreover, utilizing metamaterials and
metasurface are believed to be an efficient way to realize EM ESTXu et al.
(2018); Liang et al. (2015). Recently, it was reported that, by utilizing two-
dimensional materials (graphene and hexagonal boron nitride), an ultimate
confinement limit of the length scale of only one atom was achieved based on a
graphene-insulator-metal structureIranzo et al. (2018).

However, the EM EST always suffers from the Fabry-Perot resonance and
reflections, and the coupling efficiency of the guided modes is always low.
Therefore, to achieve EM EST, researchers are required to carefully optimize
the geometric shape and sift the materials. On the other hand, one-way
waveguide can support unidirectionally propagating EM waves without
backscattering and those one-way EM waves, as reported in many works, can
nearly perfectly bypass disorders or bends. In this paper, we theoretically
study the propagation characters in metal-dielectric-YIG-metal (MDYM)
structures and by utilizing COMSOL software, EM EST are achieved and verified
in the YIG-based one-way waveguides. Moreover, the EM EST in such waveguides
are proved to be low-loss. Furthermore, similar EM EST are achieved in a
remanence-based three-dimensional (3D) MDYM waveguide as well. This kind of
unidirectional modes-based EM EST is useful for the oncoming optical
communication system.

##  II Broadband Energy Squeezing based on unidirectional modes

As reported in many works, the effect of energy squeezing is usually achieved
in cavities surrounded by metallic materials. Here, we propose a novel
waveguide performed like a energy compressor based on unidirectional EM modes.
The inset of Fig. 1(a) shows the schematic of that waveguide, which is
composed of two layers of metal layers (treated as perfect electric
conductor(PEC) in microwave regime), one layer of dielectric (air) and one
layer of YIG under an external magnetic field $H_{0}$, i.e. the MDYM
structure. The asymmetric relative permeability tensor of YIG in this case
takes the following form

$\bar{\mu}_{\mathrm{h}}=\left[\begin{array}[]{ccc}\mu_{\mathrm{r}}&-i\mu_{\mathrm{k}}&0\\\
i\mu_{\mathrm{k}}&\mu_{\mathrm{r}}&0\\\ 0&0&1\end{array}\right]$ (1)

with
${\mu_{\mathrm{r}}}=1+\frac{{\omega_{\mathrm{m}}}{(\omega_{0}-i\alpha{\omega})}}{{{(\omega_{0}-i\alpha{\omega})}^{2}-{\omega}^{2}}}$
and
${\mu_{\mathrm{k}}}=\frac{{{\omega_{\mathrm{m}}}{\omega}}}{{{(\omega_{0}-i\alpha{\omega})}^{2}}-{\omega}^{2}}$.
${\omega}$, ${\alpha}$, ${\omega_{\mathrm{m}}}$,
${\omega_{0}}=2{\pi}{\gamma}{\mathrm{H}_{0}}$ (${\gamma}=2.8\times 10^{6}$
rad/(sG), gyromagnetic ratio) are the angular frequency, the damping
coefficient, the characteristic circular frequency and the precession angular
frequency, respectivelyDeng et al. (2015). In the MDYM model, only the
transverse electric (TE) modes ($k_{\mathrm{z}}=0$) are sustained on the air-
YIG interface, and deriving from the Maxwells’ equations and two boundary
conditions in y direction, the dispersion relation of the TE modes (SMPs) can
be written as below

$\displaystyle{\alpha_{\mathrm{d}}}{\mu_{v}}+[\frac{\alpha_{\mathrm{h}}}{\tanh({\alpha_{\mathrm{h}}}h)}+\frac{\mu_{\mathrm{k}}}{\mu_{\mathrm{r}}}k]\tanh({\alpha_{\mathrm{d}}}d)=0$
(2)

where
${\alpha_{\mathrm{d}}}=\sqrt{{k}^{2}-{\varepsilon_{\mathrm{d}}}{k_{0}}^{2}}$
(${k_{0}=\omega/c}$) and
${\alpha_{\mathrm{h}}}=\sqrt{{k}^{2}-{\varepsilon_{\mathrm{h}}}{\mu_{v}}{k_{0}}^{2}}$
(${\mu_{v}}={\mu_{\mathrm{r}}}-{\mu_{\mathrm{k}}}^{2}/{\mu_{\mathrm{r}}}$, the
Voigt permeability) demonstrate the attenuation coefficients of the SMPs in
the air and YIG layersShen et al. (2019). In Eq. (2), by letting $k\to+\infty$
and $k\to-\infty$, we calculated two asymptotic frequencies of the SMPs, i.e.
$\omega_{\mathrm{AF}}^{+}=\omega_{0}+\omega_{\mathrm{m}}$ and
$\omega_{\mathrm{AF}}^{-}=\omega_{0}+\omega_{\mathrm{m}}/2$. As an example,
the dispersion diagram of the MDYM is shown in Fig. 1(a) for
$d=0.1\lambda_{\mathrm{m}}$ ($\lambda_{\mathrm{m}}=2\pi
c/\omega_{\mathrm{m}}$) and $h=0.15\lambda_{\mathrm{m}}$. The red line
indicates the dispersion curves of SMPs and there is a clear one-way
propagation (OWP) band (shaded rectangle area) limited by two AFs. Moreover,
the cyan shaded zones and the dot-dashed line respectively represent the bulk
zones in the YIG layer and the light cone of air, respectively. Different with
the SMPs in infinite YIG based one-way waveguide, the SMPs in the MDYM
waveguide strongly coupled with bulk modes when the frequency is near
$\omega_{\mathrm{AF}}^{+}$, making the horizontal dispersion curve branch in
the $k>0$ region. To explain the theory of energy squeezing in MDYM
structures, in Fig. 1(b), the zoomed-in dispersion curves of SMPs are plotted
for three different values of $d$, i.e. $d=0.1\lambda_{\mathrm{m}}$ (red
line), $d=0.03\lambda_{\mathrm{m}}$ (black line) and
$d=0.01\lambda_{\mathrm{m}}$ (cyan line). As a result, when reducing the
thickness of the air layer, the dispersion curves rose up and the dispersion
branch of $\omega<\omega_{\mathrm{AF}}^{-}$ gradually disappeared, in another
word, the dispersion valley in thick air case ($d=0.1$) gradually rose up and
in the thin air case ($d=0.01\lambda_{\mathrm{m}}$), no SMPs are found in the
frequency band $\omega<\omega_{\mathrm{AF}}^{-}$. Most importantly, although
the dispersion rose up or dropped down when changing the values of $d$, the
whole OWP band remain the same, which implies that the EM modes with
frequencies falling within the OWP band can always propagate to only one
direction even in ultra-thin or ultra-thick air-based MDYM structure.
Therefore, as shown the inset of Fig. 1(b), we believe that the one-way
propagating EM waves in the left thicker part can always nealy perfectly
couple with the sustained EM modes in the right part of the waveguide for the
one-way propagation character in both parts.

Figure 1: (Color online) (a) Dispersion diagram of the MDYM configuration as
$h=0.15\lambda_{\mathrm{m}}$ and $d=0.1\lambda_{\mathrm{m}}$. (b) Dispersion
curves of SMPs for three different $d$, i.e. $d=0.1\lambda_{\mathrm{m}}$ (red
line), $d=0.03\lambda_{\mathrm{m}}$ (black line) and
$d=0.01\lambda_{\mathrm{m}}$ (cyan line), as $h=0.15\lambda_{\mathrm{m}}$. The
inset of (a) shows a straight MDYM waveguide while the inset of (b)
demonstrates a EM compressor based on MDYM configuration. The other parameters
are $\omega_{0}=\omega_{\mathrm{m}}$, $\lambda_{\mathrm{m}}\approx 60$ mm,
$\varepsilon_{\mathrm{d}}=1$, $\varepsilon_{\mathrm{h}}=15$, $\alpha=0$ and
$\mathrm{H}_{0}\approx 1785$ G. Figure 2: (Color online) Simulated (a)
electric field and (b) magnetic field distributions as
$d_{1}=0.1\lambda_{\mathrm{m}}$, $d_{2}=0.1\lambda_{\mathrm{m}}$ and
$h=0.15\lambda_{\mathrm{m}}$. (c) The amplitudes of the electric field and
magnetic field along the air-YIG interface. $L_{1}=L_{2}=60$ mm and the source
was located at $x=L_{1}/3$. Figure 3: (Color online) Simulations for multi-
frequencies one-way propagating EM waves for (a)
$d_{2}=0.01\lambda_{\mathrm{m}}$ and (b) $d_{2}=0.03\lambda_{\mathrm{m}}$.
Four working frequencies from top to bottom are respectively
$1.6f_{\mathrm{m}}$, $1.7f_{\mathrm{m}}$, $1.8f_{\mathrm{m}}$ and
$1.9f_{\mathrm{m}}$. The other parameters are the same as in Fig. 2. Figure 4:
(Color online) Analysis of the transmission efficiencies in simulations shown
in (a) Fig. 3(a) and (b) Fig. 3(b).

Note that the loss of MO materials, e.g. YIG, has no significant impact on the
one-way propagation character of the SMPs in one-way waveguidesTsakmakidis et
al. (2017); Xu et al. (2019). In this paper, the loss effects are considered
by setting $\alpha=0.001$ in the simulations using finite element method
(FEM). In Fig. 2, we performed a full wave simulation to reveal the high
performance energy squeezing based on our proposed MDYM configuration. Figs.
2(a) and 2(b) are the simulated electric field and magnetic field
distributions, respectively. One can clearly see that the excited EM wave
unidirectionally propagate to the end of part 1 of the waveguide and further
propagate to the narrow part (part 2) until it reaches to right most end. In
addition, no interference was observed in both Figs. 2(a) and 2(b), which is
in agreement with above analysis. More clear results of such supercoupling
between the EM modes in part 1 and part 2 are shown in Fig. 2(c),
demonstrating the electric field (blue line) and the magnetic field (red line)
distributions along the air-YIG interface. Different with the decreased
amplitude of the electric field in part 2, the magnetic field is significantly
enhanced in the narrow region and the magnetic field reaches the maximum value
on the right terminal end which is set to be a PEC wall. Due to the ability of
immune to back-scattering, this kind of energy squeezing should be high
performance and here we emphasis that, compared to the SPP- or ENZ-based
structures, our proposed thick-thin MDYM waveguides can squeeze SMPs with
frequencies falling within a continuous broadband frequency band, i.e. the OWP
band.

In Fig. 3, several FEM simulations were performed to show the broadband
squeezing property of our proposed thick-thin MDYM waveguide and in this case
we focused on the distributions of the electric field component
($E_{\mathrm{z}}$). Four working frequencies, i.e. $1.6f_{\mathrm{m}}$,
$1.7f_{\mathrm{m}}$, $1.8f_{\mathrm{m}}$ and $1.9f_{\mathrm{m}}$ were chose in
the simulations. As a result, complete one-way propagating EM waves and clear
energy squeezing were observed in the eight simulations in Fig. 3. Moreover,
all of the excited EM waves propagates to $+x$ direction and, in part 1,
according to the dispersion relation shown in Fig. 1(b) (cyan line), the
higher the frequency, the larger the values of $|k|$ (the wavenumber), which
makes the phase of $E_{\mathrm{z}}$ reaches $2\pi$ easier for the case of
$f=1.9f_{\mathrm{m}}$ than the case of $f=1.6f_{\mathrm{m}}$. In part 2, the
phase changing become complicate since the values of $k$ can be negative or
positive. In a word, Fig. 3(a) indicates that our designed MDYM structure has
the ability of squeezing EM waves with frequencies in a broadband OWP region
(in this case the bandwidth of the OWP band equals $0.5\omega_{\mathrm{m}}$).
On the other hand, as shown in Fig. 3(b), the thickness of air in part 2 of
the MDYM waveguide is three times of the one in simulations shown in Fig.
3(a), and the EM waves were squeezed in both cases, which proves that changing
the vertical thickness of the air cannot break the process of energy squeezing
in the MDYM structure. The transmission efficiencies in +x direction of the
FEM simulation in Fig. 3 are illustrated in Fig. 4 for $L_{1}\leq x\leq
L_{1}+L_{2}$. It is obvious that, for both cases, i.e.
$d_{2}=0.01\lambda_{\mathrm{m}}$ case and $d_{2}=0.03\lambda_{\mathrm{m}}$
case, the higher the frequencies of the propagating EM waves, the lesser the
propagation loss. In addition, the loss in the case of thicker air
($d_{2}=0.03$) is much lesser than the one in the case of $d_{2}=0.01$ since
more energy decayed in the YIG layer. It is worth to note that even the losses
seem like huge and only near $20\%$ energy remains for the case of
$f=1.6f_{\mathrm{m}}$ as shown in Fig. 4(a), one should notice that if we set
$L_{2}=10$ mm, more than $80\%$ energy will remain and the corresponding
transmission efficiency can be more than $93\%$ when considering
$f=1.9f_{\mathrm{m}}$. Therefore, the thick-thin MDYM structure can perfectly
squeeze unidirectionally propagating EM waves with high efficiency regardless
of the thickness of the air at subwavelength scale.

##  III Robust energy squeezing and tunneling

Figure 5: (Color online) Two kinds of U-type MDYM channel. The schematic and
simulated $E_{\mathrm{z}}$ for (a) straight channel and (b) bended channel
with $L_{2}=30$ mm and $r=7$ mm. For both (a) and (b), the working frequency
was set to be $f=1.75f_{\mathrm{m}}$. The other parameters are the same as in
Fig. 3(a).

One of the most important advantage of one-way waveguide is the robustness of
the unidirectional propagation character of SMPs, which in our opinion make it
possible to squeeze and recover the waves in a U type MDYM channel. In Fig. 5,
we designed two kinds of U-type channel and the first channel is straight with
the parameters are the same as the one in Fig. 3(a) except for the values of
$L_{2}$, while the second channel is bended with the radius of the corner
$r=7$ mm. For the first channel, we investigate the energy squeezing and
tunneling (EST) in cases of different $L_{2}$. According to the simulation
results illustrated in Fig. 5(a), in both simulations the excited waves
unidirectionally propagated to the right most end of the waveguide and were
completely squeezed in the squeezing part. More importantly, as we expect, the
EM waves recovered in the recovering part with no obvious reflection. Another
striking result shown in Fig. 5(a) is that the length of the squeezing part
only influence the phase of the EM waves at the end surface of the squeezing
part and it has no obvious impact on the supercoupling between the SMPs in the
thick and the thin MDYM waveguide. Since the supercoupling effects remain the
same regardless the values of the thickness of the air and the length of the
squeezing part, our proposed unidirectional SMPs-based EST are untouched by
the Fabry-Perot resonance. In Fig. 5(b), we designed a extraordinarily bended
U-type MDYM channel and the simulated electric field component
$E_{\mathrm{z}}$ is shown in the second diagram. According to the simulation
results, the excited unidirectional EM wave was squeezed in the squeezing part
and channelled through the half-circle shape bend, and recovered in the
recovering part of the MDYM structure. Again the loss effects can be greatly
eliminated by carefully designing the waveguide parameters $L_{1}$ and $r$.

Besides, the dispersion characters of the SMPs in the MDYM structures can be
engineered by controlling the external magnetic field, which will further
correction the OWP band. Therefore, as shown in Fig. 6, two AFs (red and black
lines) are plotted as functions of $\omega_{0}$ ($\omega_{0}=2\pi\gamma
H_{0}$) and the colored zone indicated the OWP regions for
$0<\omega_{0}\leq\omega_{\mathrm{m}}$. As a result, according to our
calculation, the robust EST can be achieved in an appreciable band
($0.5\omega_{\mathrm{m}}<\omega<2\omega_{\mathrm{m}}$) by optimizing the
external magnetic field. In order to verify the exactness of the conclusion,
we performed FEM simulations in the straight U-type MDYM structure shown in
Fig. 5(a) with $\omega_{0}=0.5\omega_{\mathrm{m}}$. Two working frequencies
are chosen, i.e. $f=1.2f_{\mathrm{m}}$ (marked by the lower green star) and
$f=1.4f_{\mathrm{m}}$ (marked by the upper green star) and for both
frequencies, the EM waves were squeezed and recorvered just like above
simulations shown in Fig. 5. Consequently, owing to the capability of
unidirectionally squeezing and recovering EM energy, the MDYM configuration
(or one-way waveguide) is proved to be a powerful and high performance
candidate in optical communication and optical integrated circuit.

Figure 6: (Color online) The AFs (or the OWP band) as a function of
$\omega_{0}$. The insets indicate the simulated electric field distribution
for $f=1.4f_{\mathrm{m}}$ (upper inset) and $f=1.2f_{\mathrm{m}}$ (lower
inset) when $\omega_{0}=0.5\omega_{\mathrm{m}}$ (green dashed line) which
corresponds to $H_{0}\approx 892.5$ G.

##  IV A 3D electromagnetic compressor based on remanence

Figure 7: (Color online) Schematic of a 3D MDYM channel based on remanence.
Figure 8: (Color online) (a) Dispersion diagram of 2D (solid lines) and 3D
(marked by stars) MDYM structures for $d_{1}=0.1\lambda_{\mathrm{r}}$ (red
lines) and for $d_{1}=0.01\lambda_{\mathrm{r}}$ (yellow line). (b) Simulated
electric field distribution in the 3D simulation as $f=0.75f_{\mathrm{r}}$
along the intersection line of a XZ plane ($y=0.005\lambda_{\mathrm{r}}$) and
a XY plane ($z=-0.5h_{0}$). The other parameters are
$\omega_{\mathrm{r}}=2\pi\times 3.587\times 10^{9}$ rad/s,
$d_{1}=0.1\lambda_{\mathrm{r}}$, $d_{2}=0.01\lambda_{\mathrm{r}}$, $h_{0}=4$
mm and $\alpha=0.001$.

YIG-based one-way waveguides, as reported in many works, sustain the
transverse electric (TE) modes rather transverse magnetic (TM) modes, making
it possible to cover the structure in the direction of $H_{0}$ with no
influence to the one-way propagation character of the SMPsXu et al. (2020).
Considering the realistic 3-dimensional (3D) condition, in this subsection, we
design the 3D MDYM channel based on remanence. Fig. 7 demonstrates the
schematic of our designed 3D channel, in which the air and YIG are surrounded
by four layers of metal and the channel is divided into three parts, i.e. the
incident part, the squeezing part and the recovering part. The lengths of
three parts are respectively $L_{1}$, $L_{2}$ and $L_{3}$ while the height of
the channel ($h_{0}$) was set to be 4 mm in this paper.

The relative permeability of YIG in this subsection takes the formPozar (2011)

$\bar{\mu}=\left[\begin{array}[]{ccc}\mu_{1}&-i\mu_{2}&0\\\
i\mu_{2}&\mu_{1}&0\\\ 0&0&1\end{array}\right]$ (3)

where $\mu_{1}=1+\frac{i\alpha\omega_{\mathrm{r}}}{\omega(1+\alpha^{2})}$ and
$\mu_{2}=-\frac{\omega_{\mathrm{r}}}{\omega(1+\alpha^{2})}$.
$\omega_{\mathrm{r}}=2\pi\times 3.587\times 10^{9}$ rad/s represents the
characteristic circular frequency. The dispersion relation in the remanence-
based 2D MDYM structure can directly derived from Eq. (2) by replacing the
attenuation efficient $\alpha_{\mathrm{h}}$ by
$\alpha_{\mathrm{h1}}=\sqrt{k^{2}-\varepsilon_{\mathrm{h}}\mu_{\mathrm{v1}}k_{0}^{2}}$
($\mu_{\mathrm{v1}}=\mu_{1}-\mu_{2}^{2}/\mu_{1}$). As shown in Fig. 8(a), the
solid lines indicate the dispersion curves of SMPs for
$d=0.1\lambda_{\mathrm{r}}$ (red line) and $d=0.01\lambda_{\mathrm{r}}$
(yellow line) and
$\lambda_{\mathrm{r}}=2\pi\mathrm{c}/\omega_{\mathrm{r}}\approx 83.6$ mm. The
stars in Fig. 8(a) represent the corresponding solutions of the 3D dispersion
equation by using the mode analysis module of COMSOL and we found that all the
stars is in the solid lines. In the 3D simulation, we set $\alpha=0.001$,
$L_{1}=L_{3}=10$ mm, $L_{2}=30$ mm and $h=0.15\lambda_{\mathrm{r}}$. The
thicknesses of the air of the incident, squeezing and recovering part are
respectively $d_{1}=0.1\lambda_{\mathrm{r}}$, $d_{2}=0.01\lambda_{\mathrm{r}}$
and $d_{3}=0.1\lambda_{\mathrm{r}}$. The inset of Fig. 8(b) shows the
distribution of the simulated electric field and we found that, similar with
above 2D simulations, the EM wave excited by the source at $x=L_{1}/2$ and it
channel through the squeezing part and transmit into the recovering part and
the energy accumulating at the end surface. The purple line indicates the
amplitudes of the electric field along the intersection line of $y=d_{2}/2$
and $z=h_{0}/2$. Moreover, the smooth line also demonstrates that the EM wave
travel through the waveguide without reflection or backscattering and the EM
energy is truly squeezed and recovered in such a remanence-based MDYM channel.

##  V Conclusion

In conclusion, we have proposed a novel metal-dielectric-YIG-metal (MDYM)
configuration and theoretically investigated the capability of such MDYM
structure in energy squeezing and tunneling (EST). Based on the theoretical
analysis and the simulations using COMSOL software, the MDYM configuration,
owing to the one-way propagation character of the SMPs, has been proved to be
high-performance in squeezing the EM energy and the corresponding transmission
efficiencies can be more than $93\%$ when the parameters of the waveguides are
carefully engineered. Besides, continuous and broadband EST were achieved
based on the one-way SMPs in MDYM structures, which, to our knowledge, has
never been reported before. The vertical and lateral lengths of the squeezing
part has also be discussed and, interestingly, different with the energy
squeezing based on ENZ materials, Fabry-Perot resonance has no obvious impact
on the capability of energy squeezing in such MDYM waveguides. Moreover, in
two kinds of our designed U-type MDYM channel, due to the supercoupling
between the one-way EM modes in incident, squeezing amd recovering parts, the
EM waves traveled unidirectionally and channelled through the U-type MDYM
channel even when the squeezing part is greatly bended. Finally, we designed a
realistic 3D MDYM channel for EST based on remanence and similar EST were
observed in such 3D structure. Our findings pave a way to squeeze energy at
ultra-subwavelength scale and the theory proposed in this paper is highly
possible to be applied in terahertz or even visible regime.

## Fundings

National Natural Science Foundation of China (NSFC) (61927813, 61865009,
11904152); Start-up funding of Southwest Medical University (20/00040186);
Department of Science and Technology of Sichuan Province (14JC0153); Science
and Technology Strategic Cooperation Programs of Luzhou Municipal People’s
Government and Southwest Medical University (2019LZXNYDJ18).

## References

  * Prange and Girvin (1987) R. E. Prange and S. M. Girvin, _The Quantum Hall Effect_ (Springer, 1987).
  * Wang et al. (2009) Z. Wang, Y. Chong, J. D. Joannopoulos, and M. Soljačić, Nature 461, 772 (2009).
  * Brion et al. (1972) J. J. Brion, R. F. Wallis, A. Hartstein, and E. Burstein, Physical Review Letters 28, 1455 (1972).
  * Wallis et al. (1974) R. Wallis, J. Brion, E. Burstein, and A. Hartstein, 9, 3424 (1974).
  * Tsakmakidis et al. (2017) K. L. Tsakmakidis, L. Shen, S. A. Schulz, X. Zheng, J. Upham, X. Deng, H. Altug, A. F. Vakakis, and R. W. Boyd, Science 356, 1260 (2017).
  * Xu et al. (2019) J. Xu, S. Xiao, C. Wu, H. Zhang, X. Deng, and L. Shen, 27, 10659 (2019), ISSN 1094-4087.
  * Xu et al. (2021) J. Xu, P. He, D. Feng, K. Yong, L. Hong, Y. Shen, and Y. Zhou, 29, 11328 (2021), ISSN 1094-4087.
  * Atkin et al. (2012) J. M. Atkin, S. Berweger, A. C. Jones, and M. B. Raschke, 61, 745 (2012), ISSN 0001-8732, 1460-6976.
  * Berweger et al. (2015) S. Berweger, J. C. Weber, J. John, J. M. Velazquez, A. Pieterick, N. A. Sanford, A. V. Davydov, B. Brunschwig, N. S. Lewis, T. M. Wallis, et al., 15, 1122 (2015), ISSN 1530-6984, 1530-6992.
  * Machida et al. (1987) S. Machida, Y. Yamamoto, and Y. Itaya, Physical Review Letters 58, 1000 (1987).
  * Frascella et al. (2021) G. Frascella, S. Agne, F. Y. Khalili, and M. V. Chekhova, 7, 72 (2021), ISSN 2056-6387.
  * Tucker and Millea (1978) J. R. Tucker and M. F. Millea, 33, 611 (1978), ISSN 0003-6951, 1077-3118.
  * Thapliyal et al. (2014) K. Thapliyal, A. Pathak, B. Sen, and J. Peřina, Physical Review A 90, 013808 (2014).
  * Silveirinha and Engheta (2006) M. Silveirinha and N. Engheta, Physical Review Letters 97, 157403 (2006).
  * Liu et al. (2008) R. Liu, Q. Cheng, T. Hand, J. J. Mock, T. J. Cui, S. A. Cummer, and D. R. Smith, Physical Review Letters 100, 023903 (2008).
  * Liu et al. (2009) L. Liu, C. Hu, Z. Zhao, and X. Luo, Optics Express 17, 12183 (2009).
  * Jingjing et al. (2010) Z. Jingjing, L. Yu, and M. Niels Asger, Optics Express 18, 3864 (2010).
  * Li et al. (2017) Z. Li, Y. Sun, H. Sun, K. Wang, J. Song, L. Liu, X. Chen, and C. Gu, Journal of Physics D: Applied Physics 50, 375105 (2017).
  * Xu et al. (2018) Z. Xu, S. Liu, S. Li, H. Zhao, L. Liu, and X. Yin, Applied Physics Express 11, 042002 (2018).
  * Liang et al. (2015) G. Liang, C. Wang, Z. Zhao, Y. Wang, N. Yao, P. Gao, Y. Luo, G. Gao, Q. Zhao, and X. Luo, Advanced Optical Materials 3, 1248 (2015).
  * Iranzo et al. (2018) D. A. Iranzo, S. Nanot, E. J. Dias, I. Epstein, C. Peng, D. K. Efetov, M. B. Lundeberg, R. Parret, J. Osmond, J.-Y. Hong, et al., Science 360, 291 (2018).
  * Deng et al. (2015) X. Deng, L. Hong, X. Zheng, and L. Shen, Applied Optics 54, 4608 (2015).
  * Shen et al. (2019) Q. Shen, L. Shen, W. Min, J. Xu, C. Wu, X. Deng, and S. Xiao, Optical Materials Express 9, 4399 (2019).
  * Xu et al. (2020) J. Xu, X. Deng, H. Zhang, C. Wu, M. Wubs, S. Xiao, and L. Shen, 22, 025003 (2020), ISSN 2040-8978.
  * Pozar (2011) D. M. Pozar, _Microwave Engineering_ (John wiley & sons, 2011).


#  A Semi-Decentralized Tikhonov-based Algorithm for Optimal Generalized Nash
Equilibrium Selection

Emilio Benenati, Wicak Ananduta, and Sergio Grammatico The authors are with
the Delft Center for Systems and Control (DCSC), TU Delft, the Netherlands.
E-mail addresses: {e.benenati, w.ananduta<EMAIL_ADDRESS>This work
was partially supported by the ERC under research project COSMOS (802348).

###### Abstract

To optimally select a generalized Nash equilibrium, in this paper, we propose
a semi-decentralized algorithm based on a double-layer Tikhonov regularization
method. Technically, we extend the Tikhonov method for equilibrium selection
in non-generalized games to the generalized case by coupling it with the
preconditioned forward-backward splitting, which guarantees linear convergence
to the solutions of the inner layer problem and allows for a semi-
decentralized implementation. We then establish a conceptual connection and
draw a comparison between the proposed algorithm and the hybrid steepest
descent method, the other known distributed framework for solving the
selection problem.

##  I Introduction

Several multi-agent decision processes can be modelled as a game, that is, a
set of inter-dependent optimization problems. In particular, if the agents are
coupled not only through their respective objective functions, but also
through a shared constraint set, then we label the setting as a _generalized_
game. Application examples for generalized games include traffic routing [1],
peer-to-peer energy markets [2, 3] and cognitive radio networks [4]. A typical
solution paradigm is the generalized Nash equilibrium (GNE), that is, an
optimal situation for each agent given the decisions of the remaining agents,
and especially the sub-class of variational GNEs (v-GNES), which has recently
received widespread attention due to its stability properties [5].

Plenty efficient v-GNE seeking algorithms, e.g. [6, 7, 8, 9, 10, 11], have
been developed for games that satisfy a monotonicity condition. Crucially,
monotone games admit in general an infinite number of v-GNEs (unless a much
more restrictive _strong_ monotonicity condition is imposed). An appealing
method to deal with the non-uniqueness of the solution is to select a GNE that
optimizes some desirable objective, as proposed in [12, 13, 14, 15]. The
selection algorithms in [12, 13] use Tikhonov’s regularization method and
builds on the literature of Variational Inequalities (VIs) to cast the
selection problem as a VI-constrained VI, which is solved by finding a
sequence of approximate solutions to regularized games. The algorithms in our
previous works [14, 15], instead, rely on fixed point selection theory and use
the hybrid steepest descent method (HSDM) [16], which pairs an appropriate
nonexpansive operator with a gradient descent. While the latter is recently
proposed as distributed algorithms for GNE selection, the former works for
non-generalized games only.

In this paper, we close the gap between the two solution frameworks by
devising a Tikhonov-based algorithm for equilibrium selection in generalized,
monotone games. Similarly to [13], we cast the GNE selection problem as a VI-
constrained VI. Then, we design a decentralized algorithm based on the
preconditioned Forward Backward (pFB) [9] to solve each regularized sub-
problem with linear convergence rate. The proposed algorithm requires a
central coordinator to compute the gradient of the selection function and to
check the termination condition for the pFB algorithm. Compared with the
Tikhonov-based GNE seeking algorithms in [6, 7] that compute the minimum-norm
v-GNE, our proposed algorithm works for general convex selection functions.
Interestingly, we also find a theoretical connection between the proposed
Tikhonov method and the HSDM. Although neither method generalizes the other,
the HSDM can be cast as a forward-backward step towards the solution of the
Tikhonov regularized inner layer problem with an appropriate choice of the
nonexpansive operator. In Section V, we present a numerical comparison between
the two methods.

#### Notation

The Euclidean inner product and norm are denoted respectively by $\langle
x,y\rangle$ and $\|\cdot\|$. For a symmetric matrix $\Psi\succ 0$, we denote
the $\Psi$-induced norm by $\|\cdot\|_{\Psi}$ and define $\langle
x,y\rangle_{\Psi}=\langle x,\Psi y\rangle$. Nonlinear set-valued operators are
denoted in calligraphic letters, e.g.
$\mathcal{T}:\mathbb{R}^{n}\rightrightarrows\mathbb{R}^{n}$. For a matrix $A$,
$\|A\|$ denotes its spectral norm. We denote the vector of all $1$ ($0$) with
dimension $n$ by $\mathbf{1}_{n}$ ($\bm{0}_{n}$). The column stack operation
is denoted by $\operatorname{col}(\cdot)$.

#### Operator theory

We denote by $\mathrm{Id}$ the identity operator. We denote the normal cone of
a closed convex set $C$ as
$\mathrm{N}_{C}(\cdot):\mathbb{R}^{n}\rightrightarrows\mathbb{R}^{n}$ [17,
Def. 6.38] and $\operatorname{proj}_{C}(x)=\operatorname{argmin}_{z\in
C}\|x-z\|$. For an operator
$\mathcal{T}\colon\mathbb{R}^{n}\rightrightarrows\mathbb{R}^{n}$, we denote
its zero set by
$\operatorname{zer}(\mathcal{T}):=\\{x\in\operatorname{dom}(\mathcal{T})\mid
0\in\mathcal{T}(x)\\}$, its fixed point set by
$\operatorname{fix}(\mathcal{T}):=\\{x\in\operatorname{dom}(\mathcal{T})\mid
x\in\mathcal{T}(x)\\}$ and its resolvant by
$\mathcal{J}_{\mathcal{T}}:=(\mathrm{Id}+\mathcal{T})^{-1}$ [17, Def. 23.1].
We refer to an operator $\mathcal{T}:C\rightrightarrows\mathbb{R}^{n}$ as:
_( $\mu$-strongly) monotone_: if, for any two pairs
$(x,y),(x^{\prime},y^{\prime})\in\operatorname{gph}(\mathcal{T}):=\\{(x,y)|x\in
C,y\in\mathcal{T}(x)\\},$ it holds that $\langle
y-y^{\prime},x-x^{\prime}\rangle\geq\mu\|x-x^{\prime}\|^{2}$, with $\mu\geq 0$
($>0$);
_Monotone-plus_ : if it is monotone and, for any two pairs
$(x,y),(x^{\prime},y^{\prime})\in\operatorname{gph}(\mathcal{T})$ such that
$\langle y-y^{\prime},x-x^{\prime}\rangle=0$, it holds
$(x,y^{\prime})\in\operatorname{gph}(\mathcal{T})$ and
$(x^{\prime},y)\in\operatorname{gph}(\mathcal{T})$;
_Lipschitz continuous_ : if there exists $L>0$, such that, for all
$x,x^{\prime}\in C$, $\|\mathcal{T}(x)-\mathcal{T}(x^{\prime})\|\leq
L\|x-x^{\prime}\|$;
_Nonexpansive_ : if Lipschitz continuous with $L=1$;
_Contractive_ : if Lipschitz continuous with $L<1$;
_$\alpha$ -averaged nonexpansive_: if for $\alpha\in(0,1)$ there exists
$\mathcal{R}\colon C\to\mathbb{R}^{n}$ nonexpansive such that
$\mathcal{T}=(1-\alpha)\mathrm{Id}+\alpha\mathcal{R}$.

#### Variational inequalities (VIs)

For a closed convex set $K$ and $\mathcal{T}:K\to K$, we denote by
$\mathrm{VI}(K,\mathcal{T})$ the problem

$\text{find}~{}\bm{x}^{*}\in K~{}\text{s.t.}~{}\inf_{\bm{x}\in
K}\langle\mathcal{T}(\bm{x}^{*}),\bm{x}-\bm{x}^{*}\rangle\geq 0,$

and its solution set by

$\mathrm{SOL}(K,\mathcal{T}):=\\{\bm{x}^{*}\in
K|\bm{x}^{*}~{}\text{solves}~{}\mathrm{VI}(K,\mathcal{T})\\}.$

A point $\hat{\bm{x}}$ is an $\varepsilon$-approximate solution of
$\mathrm{VI}(K,\mathcal{T})$ if

$\varepsilon\geq\textstyle\inf_{\bm{x}\in\mathrm{SOL}(K,\mathcal{T})}\|\bm{x}-\hat{\bm{x}}\|,$

where $\varepsilon\geq 0$ denotes the approximation error.

##  II Equilibrium selection as a variational inequality problem

We consider the multi-agent decision process in which each of $N$ agents aims
at solving an optimization problem over the decision variables
$x_{i}\in\mathcal{X}_{i}\subset\mathbb{R}^{n_{i}}$, where
$i\in\\{1,...,N\\}=:\mathcal{I}$. Let us denote
$n:=\sum_{i\in\mathcal{I}}n_{i}$ and
$\mathcal{I}_{-i}:=\mathcal{I}\setminus\\{i\\}$. Crucially, the decision
problem associated with agent $i$ is coupled to the decision variables of the
remaining agents $\bm{x}_{-i}:=(x_{j})_{j\in\mathcal{I}_{-i}}$ both through
the objective function $J_{i}:\mathbb{R}^{n}\to\mathbb{R}$, where
$n=\sum_{{i}\in\mathcal{I}}n_{i}$, and some constraints. We consider the case
in which the coupling constraints are defined via $m\in\mathbb{N}$
inequalities of the form

$\textstyle\sum_{{i}\in\mathcal{I}}A_{i}x_{i}\leq b,$ (1)

where $A_{i}\in\mathbb{R}^{m\times n_{i}}$, for each $i\in\mathcal{I}$, and
$b\in\mathbb{R}^{m}$. This problem is commonly referred to as a _generalized
game_ and we formalize it as follows:

$\displaystyle\underset{x_{i}\in\mathcal{X}_{i}}{\operatorname{min}}\quad$
$\displaystyle J_{i}(x_{i},\bm{x}_{-i})$ (2a)
$\displaystyle\operatorname{s.t.}\quad$ $\displaystyle A_{i}x_{i}\leq
b-\textstyle\sum_{{j}\in\mathcal{I}_{-i}}A_{j}x_{j}.$ (2b)

We note that we can reformulate the more general case with separable convex
coupling constraint functions into the form of (2) by augmenting the decision
variables [8, Rem. 2]. Let us now denote
$\bm{\mathcal{X}}:=\Pi_{i\in\mathcal{I}}\mathcal{X}_{i}$ and the collective
feasible set of the game in (2) as

$\Gamma:={\bm{\mathcal{X}}}\cap\\{\bm{x}~{}|~{}\eqref{eq:coupling_constr}~{}\text{holds
true}\\}.$

As a solution paradigm for the problem in (2), we consider the set of
generalized Nash equilibria, that is, a set of points from which no agent has
any incentive to unilaterally deviate:

######  Definition 1.

A set of decision variables $\bm{x}^{*}\in\bm{\mathcal{X}}$ is a GNE for the
game in (2) if, for each $i\in\mathcal{I}$,

$J_{i}(x_{i}^{*},\bm{x}^{*}_{-i})\leq J_{i}(x_{i},\bm{x}^{*}_{-i}),$

for any
$x_{i}\in\mathcal{X}_{i}\cap\\{{y\in\mathbb{R}^{n_{i}}}~{}|~{}A_{i}{y}\leq-\sum_{j\in{\mathcal{I}_{-i}}}A_{j}x_{j}^{{*}}\\}$.

The following assumptions, which are standard for GNE seeking problems [10,
11], guarantee the existence of a GNE [18, Prop. 12.11]:

######  Assumption 1.

For each $i$, $J_{i}(\cdot,\bm{x}_{-i})$ is convex and continuously
differentiable for any $\bm{x}_{-i}$.

######  Assumption 2.

For all $i\in\mathcal{I}$, $\mathcal{X}_{i}$ is compact and convex; $\Gamma$
is non-empty and it satisfies Slater’s constraint qualification.

######  Assumption 3.

The pseudogradient mapping of the game in (2), defined as

$F(\bm{x}):=\operatorname{col}((\nabla_{x_{i}}J_{i}(x_{i},\bm{x}_{-i}))_{i\in\mathcal{I}}),$

is monotone and $L_{F}$-Lipschitz continuous.

For solving the problem in (2) via a decentralized algorithm, we assume that
the agents can exchange information over an undirected, connected
communication network. We denote the set of neighbours of agent $i$ in this
network by $\mathcal{N}_{i}$. For simplicity, we consider the case where, for
all $i\in\mathcal{I}$, $J_{i}$ depends on $x_{i}$ and the decision variables
of (a subset of) $\mathcal{N}_{i}$, so that each agent is able to evaluate its
cost function by receiving the decision variables implemented by some
neighbouring agents. Furthermore, each agent maintains a local estimate of the
dual variable associated with the shared constraints in (1), denoted by
$\lambda_{i}$, as well as an auxiliary variable $\nu_{i}$, which is
communicated to the other agents in order to achieve consensus over the dual
variable estimates. We then define the following extended Karush-Kuhn-Tucker
(KKT) operator, first proposed in [9], which includes both the optimality
conditions for the coupled optimization problems in (2) and the consensus
condition:

$\displaystyle\begin{split}\mathcal{T}^{\text{KKT}}(\bm{\omega})=&\mathcal{A}(\bm{\omega})+\mathcal{B}(\bm{\omega})+\mathcal{C}(\bm{\omega}),\\\
\mathcal{A}(\bm{\omega}):=&{\mathrm{N}_{\bm{\mathcal{X}}}(\bm{x})\times\mathrm{N}_{\mathbb{R}^{|\mathcal{I}|m}_{\geq
0}}(\bm{\lambda})\times\\{\bm{0}_{|\mathcal{I}|m}\\}},\\\
\mathcal{B}:=&\begin{bmatrix}{F}(\bm{x})\\\ -\bar{\mathcal{L}}\bm{\lambda}\\\
0\end{bmatrix},\
\mathcal{C}(\bm{\omega}):=\begin{bmatrix}\bm{A}^{\top}\bm{\lambda}\\\
b-\bm{A}\bm{x}-\bar{\mathcal{L}}\bm{\nu}\\\
\bar{\mathcal{L}}\bm{\lambda}\end{bmatrix},\end{split}$ (3)

where $\bm{\omega}=\operatorname{col}(\bm{x},\bm{\lambda},\bm{\nu})$ and
$\bar{\mathcal{L}}=\mathcal{L}\otimes I_{m}$ and $\mathcal{L}$ is the
Laplacian matrix of the communication graph. It is well known that a subset of
the GNEs of (2) can be found as the zero set of $\mathcal{T}^{\text{KKT}}$ [9,
Thm. 2]. Such GNEs are known as the variational GNEs (v-GNEs) of the game in
(2) [5, Def. 3.10]. This class of solutions has been widely studied and
multiple efficient algorithms for their computation exist. However, since
Assumptions 1–3 are not enough to guarantee that
$\operatorname{zer}(\mathcal{T}^{\text{KKT}})$ is a singleton, most of the
algorithms in the literature resort to computing an unspecified v-GNE among
the possibly infinitely many.

Instead, we develop an algorithm to compute an optimally selected v-GNE
according to the selection function $\phi$, in the sense that we aim at
solving the selection problem

$\displaystyle\min_{\bm{\omega}}\phi(\bm{\omega})\quad\operatorname{s.t.}\
\bm{\omega}\in\operatorname{zer}(\mathcal{T}^{\text{KKT}}).$ (4)

Let us postulate the following assumptions on the selection function:

######  Assumption 4.

The selection function $\phi$ is convex, continuously differentiable, and
coercive. Furthermore, its gradient is $L_{\nabla\phi}$-Lipschitz continuous.

By [17, Prop. 23.39], $\operatorname{zer}(\mathcal{T}^{\text{KKT}})$ is a
closed convex set under Assumptions 1 and 3. Therefore, (4) is a convex
optimization problem. However, $\operatorname{zer}(\mathcal{T}^{\text{KKT}})$
can seldom be written in a closed form and, thus, (4) cannot be solved by
standard optimization algorithms. To derive an algorithmic solution, we
further massage its formulation by noting that
$\operatorname{zer}(\mathcal{T}^{\text{KKT}})$ is the solution set to a
Variational Inequality (VI) [19, Sec. 1.3.1]. Specifically, by defining
$\bm{\Omega}:=\bm{\mathcal{X}}\times\mathbb{R}_{\geq
0}^{Nm}\times\mathbb{R}^{Nm}$ and by [19, Eq. (1.1.3)],

$\bm{\omega}^{*}\in\operatorname{zer}(\mathcal{T}^{\text{KKT}})\Leftrightarrow\bm{\omega}^{*}\in\mathrm{SOL}(\bm{\Omega},\mathcal{B}+\mathcal{C}).$
(5)

Following the equivalence between convex optimization problems and VIs [19,
Sec. 1.3.1], we then recast (4) as the following VI-constrained VI:

${\mathrm{VI}(\mathrm{SOL}(\bm{\Omega},\mathcal{B}+\mathcal{C}),\nabla\phi)}.$
(6)

As discussed in [13], we can solve (6) by finding the
$\varepsilon_{k}$-approximate solutions to a sequence of regularized sub-
problems, indexed by $k\in\mathbb{N}$:

$\displaystyle\mathrm{VI}(\bm{\Omega},\mathcal{B}+\mathcal{C}+\gamma_{k}\nabla\phi+\alpha(\mathrm{Id}-\bm{\omega})),~{}~{}\bm{\omega}\in\bm{\Omega},$
(7)

where the sequence of regularization weights $\alpha$,
$(\gamma_{k})_{k\in\mathbb{N}}$ and approximation errors
$(\varepsilon_{k})_{k\in\mathbb{N}}$ are chosen according to the following
criteria:

######  Assumption 5.

The parameter $\alpha$ is positive; $(\gamma_{k})_{k\in\mathbb{N}}$ and
$(\varepsilon_{k})_{k\in\mathbb{N}}$ are positive and non-negative sequences
of real numbers, respectively, such that
$\textstyle\sum_{k\in\mathbb{N}}\gamma_{k}=\infty$ and eventually
$\varepsilon_{k}=0.$

By denoting ${\bar{\gamma}}:=\sup_{k\in\mathbb{N}}\gamma_{k}$, we observe the
following properties of the operators that define (7):

######  Lemma 1.

Let Assumptions 3–5 hold true. For any $k\in\mathbb{N}$, $\alpha>0$, and
$\bm{\omega}\in\bm{\Omega}$:

  1. 1.

$\mathcal{B}+\mathcal{C}+{\gamma_{k}}\nabla\phi+\alpha(\mathrm{Id}-\bm{\omega})$
is $\alpha$-strongly monotone,

  2. 2.

$\mathcal{B}+{\gamma_{k}}\nabla\phi+\alpha(\mathrm{Id}-\bm{\omega})$ is
$L_{G}$\- Lipschitz continuous, where
$L_{G}:=\max(L_{F},2|\mathcal{N}_{i}|_{i\in\mathcal{I}})+{{\bar{\gamma}}}L_{\nabla\phi}+\alpha$.

By applying [19, Cor. 2.2.5] and [19, Thm. 2.3.3a], we conclude that the
regularized problem in (7) admits one and only one solution. The following
result, which follows directly from [13, Thm. 2], formalizes a prototypical
algorithmic solution to the problem in (6) based on a sequential solution of
the problem in (7).

######  Proposition 1.

Let Assumptions 1–5 hold. Let
$(\bm{\omega}^{(k)})_{k\in\mathbb{N}}\in\bm{\Omega}$ and, for every $k$,
$\bm{\omega}^{(k+1)}$ be the $\varepsilon_{k}$-approximate solution of the VI
in (7) with $\bm{\omega}=\bm{\omega}^{(k)}$. Then, the sequence
$(\bm{\omega}^{(k)})_{k\in\mathbb{N}}$ is bounded and each of its limit points
is a solution of (6).

##  III Distributed equilibrium selection

In view of Proposition 1, next we derive a semi-decentralized algorithm for
recursively generating a sequence $(\bm{\omega}^{(k)})_{k\in\mathbb{N}}$ such
that, for all $k$, $\bm{\omega}^{(k+1)}$ is a $\varepsilon_{k}$-approximate
solution to (7) with $\bm{\omega}=\bm{\omega}^{(k)}$. We note that the exact
solution, denoted by $\bm{\omega}^{*}_{k}$, equivalently satisfies the
monotone inclusion

$\displaystyle\begin{split}&0\in(\mathcal{A}+\mathcal{B}+\mathcal{C}+\gamma_{k}\nabla\phi+\alpha(\mathrm{Id}-\bm{\omega}^{(k)}))(\bm{\omega}_{k}^{*}).\end{split}$
(8)

We then apply the preconditioned forward-backward (pFB) method on (8),
proposed in [9] in the context of decentralized GNE seeking, for computing
$\bm{\omega}^{*}_{k}$. To this end, let us define the following matrix, which
collects the step sizes of the proposed algorithm:

$\Psi=\operatorname{diag}(\bm{\rho}^{-1},\bm{\tau}^{-1},\bm{\sigma}^{-1}),$
(9)

where $\bm{\rho}:=\operatorname{diag}(\rho_{i}I_{n_{i}})_{i\in\mathcal{I}}$,
$\bm{\tau}:=\operatorname{diag}(\tau_{i}I_{n_{i}})_{i\in\mathcal{I}}$,
$\bm{\sigma}:=\operatorname{diag}(\sigma_{i}I_{n_{i}})_{i\in\mathcal{I}}$
collect the step sizes associated with the primal, dual and auxiliary
variables, respectively. Let us make the following design choice on the matrix
$\Psi$:

######  Assumption 6.

Let

$\displaystyle
r^{\mathrm{x}}_{i}:=\textstyle\max_{j=1,\dots,n_{i}}\sum_{k=1}^{m}|[A_{i}]_{jk}|,$
$\displaystyle
r^{\lambda}_{i}:=\textstyle\max_{j=1,\dots,n_{i}}\sum_{k=1}^{m}|[A_{i}]_{jk}|+2|\mathcal{N}_{i}^{\lambda}|,$
$\displaystyle r^{\nu}_{i}:=2|\mathcal{N}_{i}^{\lambda}|.$

Furthermore, let
$r=\max_{i\in\mathcal{I}}(r^{\mathrm{x}}_{i},r_{i}^{\lambda},r_{i}^{\nu})$ and
$\delta>\max(\frac{L_{G}^{2}}{\alpha},2r)$, with $L_{G}$ defined as in Lemma
1. For all $i\in\mathcal{I}$:

  1. 1.

$(2\delta-r^{\mathrm{x}}_{i})^{-1}\leq\rho_{i}\leq(\delta+r^{\mathrm{x}}_{i})^{-1}$;

  2. 2.

$(2\delta-r^{\lambda}_{i})^{-1}\leq\tau_{i}\leq(\delta+r^{\lambda}_{i})^{-1}$;

  3. 3.

$(2\delta-r^{\nu}_{i})^{-1}\leq\sigma_{i}\leq(\delta+r^{\nu}_{i})^{-1}$.

Let us also define the preconditioning matrix

$\Phi=\Psi+\begin{bmatrix}\bm{0}&-\bm{A}^{\top}&\bm{0}\\\
-\bm{A}&\bm{0}&-\bar{\mathcal{L}}\\\
\bm{0}&-\bar{\mathcal{L}}&\bm{0}\end{bmatrix}.$ (10)

######  Lemma 2.

Under Assumption 6, $\Phi\succcurlyeq\delta I$ and
$\frac{\delta}{\|\Phi\|}\geq\frac{1}{2}$.

The pFB operator for the inclusion in (8) reads, for all $k\in\mathbb{N}$, as

$\displaystyle\begin{split}\mathcal{T}^{\text{pFB}}_{{k}}=&(\mathrm{Id}+\Phi^{-1}(\mathcal{A}+\mathcal{C}))^{-1}\\\
&(\mathrm{Id}-\Phi^{-1}(\mathcal{B}+\gamma_{k}\nabla\phi+\alpha(\mathrm{Id}-\bm{\omega}^{(k)})))).\end{split}$
(11)

The following result formalizes the convergence of the fixed-point iteration
generated by $\mathcal{T}^{\text{pFB}}_{k}$ to the solution of (7).

######  Lemma 3.

Let Assumptions 1–6 hold true. Let $\bm{y}^{0}\in\bm{\Omega}$. Then, for all
$k\in\mathbb{N}$, the sequence $(\bm{y}^{(t)})_{t\in\mathbb{N}}$ generated by
the fixed-point iteration

$\bm{y}^{(t+1)}=\mathcal{T}^{\mathrm{pFB}}_{k}(\bm{y}^{(t)}),$ (12)

where $\mathcal{T}^{\text{pFB}}_{{k}}$ is defined in (11), converges linearly
in the $\Phi$-induced norm to $\bm{\omega}^{*}_{k}$ in (8) and

$\|\bm{y}^{(t)}-\bm{\omega}^{*}_{k}\|_{\Phi}\leq\|\bm{y}^{(t+1)}-\bm{y}^{(t)}\|_{\Phi}/(1-\beta),$
(13)

with $\beta=(1+\tfrac{L_{G}^{2}}{\delta^{2}}-\tfrac{2\alpha}{\|\Phi\|})$.

######  Remark 1.

Differently from the pFB operator in [9], that in (11) has additional
regularization terms, involving the gradient of the selection function, and
thus achieves the desired linear rate.

Algorithm 1 Tikhonov-regularized-pFB for optimal GNE selection

Initialization. Let $\alpha$, $(\varepsilon_{k})_{k\in\mathbb{N}}$ and
$(\gamma_{k})_{k\in\mathbb{N}}$ satisfy Assumption 5. Let
$\bm{\omega}^{(0)}\in\Pi_{i\in\mathcal{I}}({\mathcal{X}_{i}})\times{\mathbb{R}_{\geq
0}}^{|\mathcal{I}|m}\times\mathbb{R}^{|\mathcal{I}|m}$. Let
$\bm{\rho},\bm{\sigma},\bm{\tau}$ satisfy Assumption 6.
Outer iteration: for $k\in\mathbb{N}_{0}$

  1. 1.

Each agent $i\in\mathcal{I}$ sets:

$({x}^{(k,0)}_{i},{\lambda}^{(k,0)}_{i},{\nu}^{(k,0)}_{i})\leftarrow
y_{i}^{(k,0)}\leftarrow\omega_{i}^{(k)}.$ (14)

  2. 2.

Inner iteration: for $t\in\mathbb{N}_{0}$
-For each agent $i\in\mathcal{I}$:

    1. (a)

Receive $x_{j}^{(k,t)}$, $\lambda_{j}^{(k,t)},\nu_{j}^{(k,t)}$ from agent
$j\in\mathcal{N}_{i}$ and $\nabla\phi_{\omega_{i}}(\bm{y}^{(k,t)})$ from the
aggregator.

    2. (b)

Update:

$\displaystyle\begin{split}&x_{i}^{(k,t+1)}={\operatorname{proj}_{\mathcal{X}_{i}}}\left[{x}_{i}^{(k,t)}-\rho_{i}\big{(}\nabla_{x_{i}}{J}_{i}(\bm{x}^{(k,t)})+\right.\\\
&\left.~{}A_{i}^{\top}\lambda_{i}^{(k,t)}+\gamma_{k}\nabla_{x_{i}}\phi(\bm{y}^{(k,t)})+\alpha(x_{i}^{(k,t)}-{x}^{(k,0)}_{i})\right],\end{split}$
(15)
$\displaystyle\begin{split}\nu_{i}^{(k,t+1)}&=\nu_{i}^{(k,t)}-\sigma_{i}\left(\textstyle\sum_{j\in\mathcal{N}_{i}}(\lambda_{i}^{(k,t)}-\lambda_{j}^{(k,t)})\right.\\\
&\left.+\alpha(\nu_{i}^{(k,t)}-{\nu}^{(k,0)}_{i})+\gamma_{k}\nabla_{\nu_{i}}\phi(\bm{y}^{(k,t)})\right).\end{split}$
(16)

    3. (c)

Receive $\nu_{j}^{(k,t+1)}$ from agent $j\in\mathcal{N}_{i}$.

    4. (d)

Update:

$\displaystyle\begin{split}&\lambda_{i}^{(k,t+1)}=\operatorname{proj}_{\mathbb{R}^{m}_{\geq
0}}\Big{[}\Big{.}\lambda_{i}^{(k,t)}+\tau_{i}\Big{(}\Big{.}\alpha(\lambda_{i}^{(k,t)}-{\lambda}^{(k,0)}_{i})-b_{i}\\\
&+A_{i}(2x_{i}^{(k,t+1)}-x_{i}^{(k,t)})+\Big{.}\Big{.}\textstyle\sum_{j\in\mathcal{N}_{i}}\left(2\nu_{i}^{(k,t+1)}-2\nu_{j}^{(k,t+1)}\right.\\\
&\left.-\nu_{i}^{(k,t)}+\nu_{j}^{(k,t)}-\lambda_{i}^{(k,t)}+\lambda_{j}^{(k,t)}-\nabla_{\lambda_{i}}\phi(\bm{y}^{(k,t)})\right)\Big{)}\Big{]}.\end{split}$
(17)

-Coordinator:

    1. (a)

Set
$\bm{y}^{(k,t+1)}\leftarrow(\bm{x}^{(k,t+1)},\bm{\lambda}^{(k,t+1)},\bm{\nu}^{(k,t+1)})$.

    2. (b)

Communicate $\nabla_{\omega_{i}}\phi(\bm{y}^{(k,t+1)})$ to each agent
$i\in\mathcal{I}$.

    3. (c)

If the following is satisfied,

$\|\bm{y}^{(k,t+1)}-\bm{y}^{(k,t)}\|_{\Phi}\leq(1-\beta)\varepsilon_{k},$ (18)

terminate inner iteration. Each agent then sets

$\omega_{i}^{(k+1)}=y_{i}^{(k,t+1)}.$ (19)

The linear convergence of $\mathcal{T}^{\text{pFB}}_{{k}}$ guarantees that,
for $\varepsilon_{k}>0$, the iterate $\bm{\omega}^{(k+1)}$ in Proposition 1 is
found within a finite number of inner iterations, whose termination can be
based on a simple stopping criterion derived from (13). The proposed method,
which results from the expansion of the pFB operator, is illustrated in
Algorithm 1, where we use $t$ as the inner iteration index.

######  Proposition 2.

Let $(\bm{\omega}^{(k)})_{k\in\mathbb{N}}$ be generated by Algorithm 1. Under
Assumptions 1–6, for each $k$, $\bm{\omega}^{(k+1)}$ is an
$\varepsilon_{k}$-approximate solution of (7) with
$\bm{\omega}=\bm{\omega}^{(k)}$ and, if $\varepsilon_{k}>0$, the condition in
(18) is verified in a finite number of steps.

######  Remark 2.

By Proposition 1, Assumption 5 requires the algorithm to eventually seek an
exact solution of a VI. This is a stringent requirement, as the pFB algorithm
only achieves such a solution asymptotically. The authors of [13],
alternatively, consider the less restrictive condition
$\lim_{k\xrightarrow{}\infty}\tfrac{\varepsilon_{k}}{\gamma_{k}}=0$ in the
case where the definition set of the VI-constrained VI is compact. Although
$\bm{\Omega}$ does not satisfy this condition as the dual variables belong to
an unbounded set, in practice the dual variables are bounded in view of [20,
Prop. 3.3].

##  IV Theoretical connection and comparison with the Hybrid Steepest Descent
Method

In our previous works [14, 15], we take a different algorithm design path,
where we reformulate Problem (4) as a fixed point selection problem. This
reformulation allows one to use the hybrid steepest descent method (HSDM) [16]
to solve (4). Specifically, one has to find a quasi-shrinking mapping
$\mathcal{T}$ [15, Def. 1] such that
$\operatorname{fix}(\mathcal{T})=\operatorname{zer}(\mathcal{T}^{\text{KKT}})$.
Then, the limit point of the sequence $(\bm{z}^{(k)})_{k\in\mathbb{N}}$,
defined by

$\bm{z}^{(k+1)}=\mathcal{T}(\bm{z}^{(k)})-\gamma_{k}\nabla\phi(\mathcal{T}\bm{z}^{(k)}),$
(20)

converges to the solution of (4) if $(\gamma_{k})_{k\in\mathbb{N}}$ is square-
summable but non-summable [16, Thm. 5]. The vanishing weight on $\nabla\phi$
is reminiscent of the Tikhonov regularization introduced in Section II.
Indeed, the two methods are related, as shown next.

Let us consider the exact solution to the VI in (7) with
$\bm{\omega}=\bm{\omega}^{(k)}$. Then, from [19, Prop. 12.3.6],

$\bm{\omega}^{(k+1)}=\mathcal{J}_{\frac{1}{\alpha}(\mathcal{A}+\mathcal{B}+\mathcal{C}+\gamma_{k}\nabla\phi)}(\bm{\omega}^{(k)})=:\mathcal{T}^{\text{Tik}}_{k}(\bm{\omega}^{(k)}).$
(21)

The properties of the HSDM update in (20) depend on the choice of
$\mathcal{T}$. Let us consider the particular case
$\mathcal{T}=\mathcal{J}_{\mathcal{A}+\mathcal{B}+\mathcal{C}}$, which is
averaged and, thus, can be shown to be quasi-shrinking by [15, Lem. 1]. We
rewrite (20) with this particular choice of $\mathcal{T}$ as:

$\displaystyle\begin{split}&\bm{v}^{(k+1)}=\mathcal{J}_{\mathcal{A}+\mathcal{B}+\mathcal{C}}(\bm{z}^{(k)}),\\\
&\bm{z}^{(k+1)}=\bm{v}^{(k+1)}-\gamma_{k}\nabla\phi(\bm{v}^{(k+1)}),\end{split}$
(22)

where we introduced the auxiliary sequence $(\bm{v}^{(k)})_{k\in\mathbb{N}}$.
By rearranging (22), we note that this sequence evolves as

$\bm{v}^{(k+1)}=\mathcal{J}_{\mathcal{A}+\mathcal{B}+\mathcal{C}}\circ(\mathrm{Id}-\gamma_{k}\nabla\phi)(\bm{v}^{(k)})=:\mathcal{T}^{\text{HSDM}}_{k}(\bm{v}^{(k)}).$
(23)

With this particular choice of $\mathcal{T}$, the operator
$\mathcal{T}^{\text{HSDM}}_{k}$ in (23) corresponds to the forward-backward
(FB) splitting method applied to the inclusion

$0\in\mathcal{A}+\mathcal{B}+\mathcal{C}+\gamma_{k}\nabla\phi.$

From [17, Prop. 26.1iv] and [17, Prop. 23.38],

$\operatorname{fix}(\mathcal{J}_{\frac{1}{\alpha}(\mathcal{A}+\mathcal{B}+\mathcal{C}+\gamma_{k}\nabla\phi)})=\operatorname{fix}(\mathcal{J}_{\mathcal{A}+\mathcal{B}+\mathcal{C}}\circ(\mathrm{Id}-\gamma_{k}\nabla\phi)),$

which implies

$\operatorname{fix}(\mathcal{T}^{\text{Tik}}_{k})=\operatorname{fix}(\mathcal{T}^{\text{HSDM}}_{k}).$
(24)

Thus, we conclude that both the Tikhonov update in (21) and the HSDM step in
(20) (in terms of the auxiliary sequence $\bm{v}^{(k)}$) apply at each step
$k$ a single update of a fixed point iteration, and the two operators in (24)
have the same fixed point set. This analogy is only theoretical, as in
practice the operator $\mathcal{J}_{\mathcal{A}+\mathcal{B}+\mathcal{C}}$
cannot be implemented in a distributed fashion. Nevertheless, it outlines that
both the Tikhonov method and the HSDM function by the same underlying
principle of tracking the solutions to a sequence of regularized problems.
Moreover, we note that (24) does not hold for a generic choice of
$\mathcal{T}$, thus the HSDM includes algorithms that are not covered by the
Tikhonov method. On the other hand, clearly the operator
$\mathcal{T}^{\text{Tik}}_{k}$ cannot be rewritten in terms of an FB operator.
Therefore, we should conclude that neither method is a generalization of the
other.

The key differences between the two frameworks are summarized in Table I. In
addition, we remark that the Tikhonov framework can be paired with any
(splitting) methods for strongly monotone games to obtain a decentralized
algorithm. Meanwhile, the HSDM requires methods for monotone games that are
quasi-shrinking, such as the forward-backward-forward (FBF) splitting [17,
Sec. 26.6]. Thus, Tikhonov-based methods can benefit from a larger pool of
available algorithms for future development.

| Tikhonov | Tikhonov | HSDM
---|---|---|---
| (bounded set)
parameters |  $\sum\gamma_{k}=\infty$, |  $\sum\gamma_{k}=\infty$, |  $\sum\gamma_{k}=\infty$,
| $\varepsilon_{k}=0$ | $\lim_{k\to\infty}\frac{\varepsilon_{k}}{\gamma_{k}}=0$ | $\sum\gamma_{k}^{2}<\infty$
$t$ (#inner | $\infty$ | $t\to\infty$ | $1$
iterations) | eventually
$\phi$ coercive | yes | no | no
TABLE I: Theoretical property differences of the Tikhonov method and HSDM for
GNE selection.

##  V Numerical simulations

We test the proposed algorithm on $100$ game instances with $10$ agents, where
the pseudogradient is in the form

$F(\bm{x})=Q_{F}\bm{x}+\bm{c}_{F}.$

The parameters $Q_{F}$, which is singular and positive semi-definite, and
$\bm{c}_{F}$ are randomly generated. We define the selection function

$\phi(\bm{\omega})=\|\bm{x}\|^{2}_{Q_{\phi}}+\bm{c}_{\phi}^{\top}\bm{x}+\theta(\|\bm{\lambda}\|^{2}+\|\bm{\nu}\|^{2}),$

where $Q_{\phi}\succcurlyeq 0$ and $\bm{c}_{\phi}$ are as well randomly
generated and $\theta=10^{-3}$. For all $i$, we define the local constraint
set ${\mathcal{X}_{i}=}\\{\mathbb{R}^{5}:\|x_{i}\|_{\infty}\leq 1\\}.$
Furthermore, we let $A_{i}=I$, for all $i\in\mathcal{I}$, and
$b=2\cdot\mathbf{1}_{5}$.We then set $(\gamma_{k})_{k\in\mathbb{N}}$ and
$(\varepsilon_{k})_{k\in\mathbb{N}}$ in Algorithm 1 to

$\gamma_{k}=10^{-3}k^{-\xi};~{}~{}\varepsilon_{k}=\begin{cases}10^{-3}k^{-\xi\zeta}&\text{if}~{}10^{-3}k^{-\xi\zeta}\geq\underline{\varepsilon}\\\
0&\text{if}~{}10^{-3}k^{-\xi\zeta}<\underline{\varepsilon}\end{cases}$

where $\underline{\varepsilon}$ is set to the computer numerical precision.
The parameter $\xi$ controls the decay of the regularization weight
$\gamma_{k}$, while $\zeta$ controls the decay rate of $\varepsilon_{k}$ with
respect to $\gamma_{k}$. We evaluate Algorithm 1 in terms of the residual
computed for each outer iteration $k$ and inner iteration $t$ as

$\displaystyle\mathfrak{R}(t)=\|{\bm{y}^{(k,t)}}-\mathcal{(}\mathrm{Id}+\mathcal{A})^{-1}\circ(\mathrm{Id}-\mathcal{B}-\mathcal{C})({\bm{y}^{(k,t)}})\|,$

and in terms of the reduction of selection function $\phi$ with respect to the
value obtained by the standard FBF algorithm (without optimal selection) [10,
Alg. 2] in Figures 1–3. we observe that, for decreasing values of $\xi$, the
algorithm achieves a lower selection function value and a larger residual (cf.
Figure 1). This trade-off between convergence to a GNE (measured by the
residual) and convergence to a $\phi$-optimal point is expected, because a too
slow decay of the regularization weight $\gamma_{k}$ leads the algorithm to
disregard the GNE seeking in order to compute the unconstrained optimal value
of $\phi$. However, if $\xi$ is set too high, the algorithm achieves a value
of $\phi$ higher than the one computed without selection. Moreover, an
increasing value of $\alpha$ improves the algorithm performance (cf. Figure
2). Finally, we observe that for increasing values of $\zeta$, the algorithm
reaches a higher residual and a lower selection function value (cf. Figure 3).
In our experience, setting $\zeta$ too high might cause convergence failure as
$\varepsilon_{k}$ might become 0 before $\gamma_{k}$ reaches a negligible
value. In Figure 4 we compare Algorithm 1 with the HSDM paired with the FBF
algorithm [15, Alg. 1]. The parameters for Algorithm 1 are chosen among the
ones that performed reasonably well in both the performance metrics considered
in Figures 1–3. We find the two algorithms to have similar convergence speed.
One might find this surprising, as the Tikhonov method is double-layered;
thus, one could expect slower convergence compared to the single-layer HSDM.
This can be explained by noting that the slowdown caused by the double-layer
iterations is compensated by the linear convergence of the pFB. In contrast,
the HSDM uses the FBF, which only achieves sublinear convergence.
Nevertheless, as observed in the first set of simulations, the Tikhonov method
requires more careful parameter tuning than the HSDM.

Figure 1: Average performance of Algorithm 1 for $\zeta=2$ and $\alpha=1$.
Figure 2: Average performance of Algorithm 1 for $\xi=0.6$ and $\zeta=2$.
Figure 3: Average performance of Algorithm 1 for $\alpha=1$ and $\xi=0.6$.
Figure 4: Comparison between Algorithm 1 with $\xi=0.6$, $\zeta=2$, $\alpha=1$
and the HSDM-FBF method [15, Alg. 1]. The x-axis indicates the cumulative
number of inner iterations for Algorithm 1.

##  VI Conclusion

The generalized Nash equilibrium selection problem can be solved with a semi-
decentralized algorithm based on the Tikhonov regularization method combined
with the preconditioned forward backward algorithm, which achieves linear
convergence for the regularized sub-problems. Interestingly, both the Tikhonov
regularization method and a particular instance of the hybrid steepest descent
method seek at each iteration an approximate solution to the same regularized
problem, indicating a conceptual connection. Although theoretically less
practical (as shown in Table I), the Tikhonov algorithm demonstrates
comparable convergence performance to the state-of-the-art in our simulations.

###  -A Proof of Lemma 1

The result 1) is immediate from the definition of strong monotonicity by
applying the monotonicity of $\mathcal{B}+\mathcal{C}$ and $\nabla\phi$, which
follow from [9, Lem. 5] and Assumption 4, respectively. From Gerschgorin’s
theorem, $\|\bar{\mathcal{L}}\|\leq
2\max(|\mathcal{N}_{i}|_{i\in\mathcal{I}})$. Thus, $\mathcal{B}$ is
$\max(L_{F},2|\mathcal{N}_{i}|_{i\in\mathcal{I}})$-Lipschitz continuous from
Assumption 3. The result 2) then follows immediately. $\hfill\blacksquare$

###  -B Proof of Proposition 1

Let us denote the solution set of (6) by $\mathcal{S}$. The proposition
follows from [13, Thm. 2] if $\bm{\Omega}$ is closed and convex;
$\mathcal{B}+\mathcal{C}$ is continuous and monotone; $\nabla\phi$ is
continuous and monotone plus; $\mathcal{S}$ is bounded and not empty; and the
set

$\mathfrak{L}_{1}:=\\{\bm{\omega}_{2}\in\bm{\Omega}|~{}\exists~{}\bm{\omega}_{1}\in\mathcal{S}~{}\text{s.
t.}~{}\nabla\phi(\bm{\omega}_{2})^{\top}(\bm{\omega}_{1}-\bm{\omega}_{2})>0\\}$

is bounded. The conditions on $\bm{\Omega}$ and the continuity of
$\mathcal{B},\mathcal{C},\nabla\phi$ follow directly from the assumptions. The
monotonicity of $\mathcal{B}+\mathcal{C}$ follows from that of skew-symmetric
linear operators [17, Ex. 20.35] as in [9, Lemma 5]. $\nabla\phi$ is monotone-
plus from the convexity of $\phi$ and [17, Ex. 22.4i]. From [9, Thm. 2],
$\operatorname{zer}(\mathcal{T}^{\text{KKT}})\neq\varnothing$. Let
$\bm{\omega}_{1}\in\operatorname{zer}(\mathcal{T}^{\text{KKT}})$ and consider
the set

$\mathfrak{L}_{2}:=\\{\bm{\omega}_{2}\in\bm{\Omega}|\nabla\phi(\bm{\omega}_{2})^{\top}(\bm{\omega}_{1}-\bm{\omega}_{2})\geq
0\\}.$

We apply the convexity inequality on $\phi$ to find

$\phi(\bm{\omega}_{1})-\phi(\bm{\omega}_{2})\geq\nabla\phi(\bm{\omega}_{2})^{\top}(\bm{\omega}_{1}-\bm{\omega}_{2}),\quad{\forall}\bm{\omega}_{2}\in\bm{\Omega}.$
(25)

From the coercivity of $\phi$, $\phi(\bm{\omega}_{1})-\phi(\bm{\omega}_{2})<0$
for sufficiently large $\|\bm{\omega}_{2}\|$. Thus, the set $\mathfrak{L}_{2}$
is bounded. Therefore, from [19, Prop. 2.2.3], $\mathcal{S}$ is non-empty and
compact. Consequently, as (25) holds for any $\bm{\omega}_{1}\in\mathcal{S}$,
we find for all $\bm{\omega}_{2}\in\bm{\Omega}$,

$\max_{\bm{\omega}_{1}\in\mathcal{S}}\nabla\phi(\bm{\omega}_{2})^{\top}(\bm{\omega}_{1}-\bm{\omega}_{2})\leq\phi^{*}-\phi(\bm{\omega}_{2}),$

where $\phi^{*}=\max_{\bm{\omega}_{1}\in\mathcal{S}}\phi(\bm{\omega}_{1})$. By
the coercivity of $\phi$, $\phi^{*}-\phi(\bm{\omega}_{2})<0$ for sufficiently
large $\|\bm{\omega}_{2}\|$. Therefore, $\mathfrak{L}_{1}$ is bounded.
$\hfill\blacksquare$

###  -C Proof of Lemma 2

$\Phi\succcurlyeq\delta I$ follows from [9, Lem. 6]. From Gerschgorin’s disc
theorem,

$\|\Phi\|\leq\max_{i\in\mathcal{I}}(\max(r^{\mathrm{x}}_{i}+\rho_{i}^{-1},r^{\lambda}_{i}+\tau_{i}^{-1},r^{\nu}_{i}+\nu_{i}^{-1}))\leq
2\delta.$

$\hfill\blacksquare$

###  -D Proof of Lemma 3

The operator $(\mathrm{Id}+\Phi^{-1}(\mathcal{A}+\mathcal{C}))$ is non-
expansive in the $\Phi$-induced norm, following [9, Lem. 7ii]. For
compactness, let us denote
$G=\mathcal{B}+\gamma_{k}\nabla\phi+\alpha(\mathrm{Id}-\bm{\omega}^{(k)})$. We
proceed by proving that $(\mathrm{Id}-\Phi^{-1}G)$ is contractive. First, we
observe

$\displaystyle\|\Phi\|I\succcurlyeq\Phi\succcurlyeq\delta
I\Rightarrow\|\Phi\|\|\bm{z}\|^{2}\geq\|\bm{z}\|^{2}_{\Phi}\geq\delta\|\bm{z}\|^{2},$
(26a) $\displaystyle\Phi\succcurlyeq\delta
I\Rightarrow\Phi^{-1}\preceq\delta^{-1}I\Rightarrow\|\bm{z}\|^{2}_{\Phi^{-1}}\leq\delta^{-1}\|\bm{z}\|^{2},$
(26b)

for any $\bm{z}\in\bm{\Omega}$. Furthermore, for any pair
$\bm{z},\bm{z}^{\prime}\in\bm{\Omega}$,

$\displaystyle\langle\bm{z}-\bm{z}^{\prime},\Phi^{-1}(G\bm{z}-G\bm{z}^{\prime})\rangle_{{\Phi}}$
$\displaystyle\geq\langle\bm{z}-\bm{z}^{\prime},\Phi^{-1}\gamma_{k}(\nabla\phi(\bm{z})-\nabla\phi(\bm{z}^{\prime}))+\alpha\Phi^{-1}(\bm{z}-\bm{z}^{\prime})\rangle_{\Phi}$
$\displaystyle=\langle\bm{z}-\bm{z}^{\prime},\gamma_{k}(\nabla\phi(\bm{z})-\nabla\phi(\bm{z}^{\prime}))+\alpha(\bm{z}-\bm{z}^{\prime})\rangle$
$\displaystyle\geq\alpha\|\bm{z}-\bm{z}^{\prime}\|^{2}\geq\tfrac{\alpha}{\|\Phi\|}\|\bm{z}-\bm{z}^{\prime}\|_{\Phi}^{2}.$
(27)

We use the fact that $\Phi^{-1}\mathcal{B}$ is cocoercive in the
$\Phi$-induced norm [9, Lem. 7i] (and, thus, monotone) in the first
inequality, the monotonicity of $\nabla\phi$ in Euclidean norm in the second
inequality, and (26a) in the third inequality. We then have that

$\displaystyle\begin{split}&\|\Phi^{-1}G(\bm{z})-\Phi^{-1}G(\bm{z}^{\prime})\|_{\Phi}^{2}=\|G(\bm{z})-G(\bm{z}^{\prime})\|_{\Phi^{-1}}^{2}\overset{\eqref{eq:norm_ineq_2}}{\leq}\\\
&\tfrac{1}{\delta}\|G(\bm{z})-G(\bm{z}^{\prime})\|_{2}^{2}{\leq}\tfrac{L_{G}^{2}}{\delta}\|\bm{z}-\bm{z}^{\prime}\|_{2}^{2}\overset{\eqref{eq:norm_ineq_1}}{\leq}\tfrac{L_{G}^{2}}{\delta^{2}}\|\bm{z}-\bm{z}^{\prime}\|_{\Phi}^{2},\end{split}$
(28)

where we use the Lipschitz continuity of $G$ (Lemma 1). By expanding the
square and from (27) and (28), we have that

$\displaystyle\|(\mathrm{Id}-\Phi^{-1}G)(\bm{z})-(\mathrm{Id}-\Phi^{-1}G)(\bm{z}^{\prime})\|^{2}_{\Phi}$
$\displaystyle=\|\bm{z}-\bm{z}^{\prime}\|_{\Phi}^{2}+\|\Phi^{-1}G(\bm{z})-\Phi^{-1}G(\bm{z}^{\prime})\|_{\Phi}^{2}$
$\displaystyle\quad\
{-2}\langle\bm{z}-\bm{z}^{\prime},\Phi^{-1}(G\bm{z}-G\bm{z}^{\prime})\rangle_{\Phi}$
$\displaystyle\leq(1+\tfrac{L_{G}^{2}}{\delta^{2}}-\tfrac{2\alpha}{\|\Phi\|})\|\bm{z}-\bm{z}^{\prime}\|_{\Phi}^{2}.$

From Assumption 6 and Lemma 2,

$\frac{L_{G}^{2}}{\delta^{2}}\leq\frac{L_{G}^{2}\alpha}{L_{G}^{2}\delta}\leq\frac{2\alpha}{\|\Phi\|},$

which implies that $\Phi^{-1}G$ is contractive in $\Phi$-norm. Finally, from
[17, Prop. 26.1iv] and [19, Eq. (1.1.3)],
$\operatorname{fix}(\mathcal{T}^{\mathrm{pFB}}_{k})=\operatorname{zer}(\mathcal{A}+\mathcal{C}+G)=\mathrm{SOL}(\bm{\Omega},\mathcal{C}+G)$
and the thesis follows immediately from [17, Thm. 1.50 (iii) and (v)].
$\hfill\blacksquare$

###  -E Proof of Proposition 2

For any $k\in\mathbb{N}$, denote
$G_{k}=\mathcal{B}+\gamma_{k}\nabla\phi+\alpha(\mathrm{Id}-\bm{\omega}^{(k)})$
and $\bm{y}^{(k,0)}=\bm{\omega}^{(k)}$. Let us consider the iteration

$\bm{y}^{(k,t+1)}=\mathcal{T}^{\mathrm{pFB}}_{k}(\bm{y}^{(k,t)}),$ (29)

and denote
$\bm{y}^{(k,t)}=(\bm{\lambda}^{(k,t)},\bm{y}^{(k,t)},\bm{\nu}^{(k,t)})$. From
(11),

$(\mathrm{Id}-\Phi^{-1}G_{k})(\bm{y}^{(k,t)})\in(\mathrm{Id}+\Phi^{-1}(\mathcal{A}+\mathcal{C}))(\bm{y}^{(k,t+1)}).$

By multiplying each side by $\Phi$ and by substituting the definition of
$\mathcal{A},\mathcal{C}$ and $\Phi$ in the latter, we obtain

$\displaystyle(\Phi-G_{k})(\bm{y}^{(k,t)})\in\Psi\bm{y}^{(k,t+1)}+$
$\displaystyle\begin{bmatrix}(\mathrm{N}_{\bm{\mathcal{X}}}(\bm{x}^{(k,t+1)}))\\\
\mathrm{N}_{\mathbb{R}^{|\mathcal{I}|m}_{\geq 0}}(\bm{\lambda}^{(k,t+1)})\\\
\bm{0}_{|\mathcal{I}|m}\end{bmatrix}+\begin{bmatrix}\mathbf{0}_{{n}}\\\
b-2\bm{A}\bm{x}^{(k,t+1)}-2\bar{\mathcal{L}}\bm{\nu}^{(k,t+1)}\\\
\mathbf{0}_{{|\mathcal{I}|m}}\end{bmatrix}.$

By multiplying both sides by $\Psi^{-1}$, rearranging the terms and by
applying $(\mathrm{Id}+\mathrm{N}_{C})^{-1}=\operatorname{proj}_{C}$ [17, Ex.
23.4] for any $C$ closed, convex, it is immediate to verify that the updates
in (LABEL:eq:primal_update)–(17) are equivalent to the iteration (29). From
(18), (19), and Lemma 3,

$\displaystyle\|\bm{\omega}^{(k+1)}-\bm{\omega}_{k}^{*}\|_{\Phi}$
$\displaystyle=\|\bm{y}^{(k,t+1)}-\bm{\omega}_{k}^{*}\|_{\Phi}$
$\displaystyle\leq\|\bm{y}^{(k,t+1)}-\bm{y}^{(k,t)}\|_{\Phi}/(1-\beta)\leq\varepsilon_{k}.$

Thus, $\bm{\omega}^{(k+1)}$ is an $\varepsilon_{k}$-approximate solution.
Next, by triangle inequality, we have that

$\displaystyle\|\bm{y}^{(k,t+1)}-\bm{y}^{(k,t)}\|_{\Phi}$
$\displaystyle\leq\|\bm{y}^{(k,t+1)}-\bm{\omega}_{k}^{*}\|_{\Phi}+\|\bm{\omega}_{k}^{*}-\bm{y}^{(k,t)}\|_{\Phi}$
$\displaystyle\overset{\\{1\\}}{\leq}(1+\beta)\|\bm{\omega}_{k}^{*}-\bm{y}^{(k,t)}\|_{\Phi}$
$\displaystyle\overset{\\{2\\}}{\leq}\beta^{t}(\beta+1)\|\bm{\omega}_{k}^{*}-\bm{y}^{(k,0)}\|_{\Phi},$

where $\\{1\\}$ follows from the contractivity of
$\mathcal{T}^{\mathrm{pFB}}_{k}$ (see proof of Lemma 3) and $\\{2\\}$ follows
from [17, Thm. 1.50iii]. Thus, if $\varepsilon_{k}>0$, (18) is satisfied for

$t\geq\log_{\beta}\left(\tfrac{\varepsilon_{k}}{(\beta+1)\|\bm{\omega}_{k}^{*}-\bm{y}^{(k,0)}\|_{\Phi}}\right).$

$\hfill\blacksquare$

## References

  * [1] B. G. Bakhshayesh and H. Kebriaei, “Decentralized equilibrium seeking of joint routing and destination planning of electric vehicles: A constrained aggregative game approach,” IEEE Transactions on Intelligent Transportation Systems, pp. 1–10, 2021.
  * [2] G. Belgioioso, W. Ananduta, S. Grammatico, and C. Ocampo-Martinez, “Operationally-Safe Peer-to-Peer Energy Trading in Distribution Grids: A Game-Theoretic Market-Clearing Mechanism,” IEEE Transactions on Smart Grid, 2022.
  * [3] W. Ananduta and S. Grammatico, “Equilibrium seeking and optimal selection algorithms in peer-to-peer energy markets,” Games, vol. 13, no. 5, 2022\.
  * [4] J. Wang, M. Peng, S. Jin, and C. Zhao, “A Generalized Nash Equilibrium Approach for Robust Cognitive Radio Networks via Generalized Variational Inequalities,” IEEE Transactions on Wireless Communications, vol. 13, no. 7, pp. 3701–3714, 2014.
  * [5] F. Facchinei and C. Kanzow, “Generalized Nash equilibrium problems,” Annals of Operations Research, vol. 175, pp. 177–211, Mar. 2010.
  * [6] H. Yin, U. V. Shanbhag, and P. G. Mehta, “Nash equilibrium problems with scaled congestion costs and shared constraints,” IEEE Transactions on Automatic Control, vol. 56, no. 7, pp. 1702–1708, 2011.
  * [7] A. Kannan and U. V. Shanbhag, “Distributed computation of equilibria in monotone nash games via iterative regularization techniques,” SIAM Journal on Optimization, vol. 22, no. 4, pp. 1177–1205, 2012.
  * [8] S. Grammatico, “Proximal dynamics in multiagent network games,” IEEE Transactions on Control of Network Systems, vol. 5, no. 4, pp. 1707–1716, 2018\.
  * [9] P. Yi and L. Pavel, “An operator splitting approach for distributed generalized Nash equilibria computation,” Automatica, vol. 102, pp. 111–121, 2019.
  * [10] B. Franci, M. Staudigl, and S. Grammatico, “Distributed forward-backward (half) forward algorithms for generalized Nash equilibrium seeking,” in Proceedings of the 2020 European Control Conference (ECC), pp. 1274–1279, IEEE, 2020.
  * [11] G. Belgioioso, A. Nedić, and S. Grammatico, “Distributed generalized Nash equilibrium seeking in aggregative games on time-varying networks,” IEEE Transactions on Automatic Control, vol. 66, no. 5, pp. 2061–2075, 2020\.
  * [12] G. Scutari, F. Facchinei, J.-S. Pang, and L. Lampariello, “Equilibrium selection in power control games on the interference channel,” in 2012 Proceedings IEEE INFOCOM, pp. 675–683, IEEE, Mar. 2012.
  * [13] F. Facchinei, J.-S. Pang, G. Scutari, and L. Lampariello, “VI-constrained hemivariational inequalities: distributed algorithms and power control in ad-hoc networks,” Mathematical Programming, vol. 145, pp. 59–96, June 2014\.
  * [14] E. Benenati, W. Ananduta, and S. Grammatico, “On the optimal selection of generalized Nash equilibria in linearly coupled aggregative games,” in 2022 IEEE 61st Conference on Decision and Control (CDC), pp. 6389–6394, IEEE, Dec. 2022.
  * [15] E. Benenati, W. Ananduta, and S. Grammatico, “Optimal selection and tracking of generalized Nash equilibria in monotone games,” Mar. 2022.  arXiv:2203.07765.
  * [16] I. Yamada and N. Ogura, “Hybrid steepest descent method for variational inequality problem over the fixed point set of certain quasi-nonexpansive mappings,” Numerical Functional Analysis and Optimization, vol. 25, no. 7-8, pp. 619–655, 2005.  Publisher: Taylor & Francis.
  * [17] H. H. Bauschke and P. L. Combettes, Convex Analysis and Monotone Operator Theory in Hilbert Spaces.  CMS Books in Mathematics, Springer international publishing, 2017.
  * [18] D. P. Palomar and Y. C. Eldar, Convex Optimization in Signal Processing and Communications.  Cambridge university press, 2009.
  * [19] F. Facchinei and J.-S. Pang, Finite-dimensional variational inequalities and complementarity problems.  Springer series in operations research, New York: Springer, 2003.
  * [20] A. Auslender and M. Teboulle, “Lagrangian duality and related multiplier methods for variational inequality problems,” SIAM Journal on Optimization, vol. 10, pp. 1097–1115, Jan. 2000.