\section{Introduction}
In order to deal with the challenges arising from the progress in autonomous driving, classical approaches such as the Lighthill-Whitham-Richards (LWR) model \cite{lighthill1955kinematic, richards1956shockwaves} have been extended to include more information on the surrounding traffic, see for example~\cite{BlandinGoatin2016,friedrich2018godunov,KeimerPflug2017}.
These models are {nonlocal} traffic flow models.
Here, the flux function depends on an integral evaluation of the density or velocity.
In case of autonomous vehicles the integration area allows for an interpretation as a connection radius.

Nonlocal traffic flow models have been studied in various research directions over the years, e.g. well-posedness~\cite{chiarello2018global,friedrich2018godunov,KeimerPflug2017}, numerical schemes~\cite{ ChalonsGoatinVillada2018,friedrich2019maximum,friedrich2018godunov}, its singular limit behavior \cite{coclite2020general,colombo2019singular,keimer2019local} or modeling extensions~\cite{chiarello2019multiclass, keimer2019nonlocal, bayen2022multilane,friedrich2020nonlocal,friedrich2020onetoone}.
Nevertheless, there are only a few works concerning control problems \cite{bayen2021boundary, friedrich2022lyapunov, huang2020stability, karafyllis2022control}.
The works \cite{huang2020stability, karafyllis2022control} consider Lyapunov functions on a ring road.
In \cite{bayen2021boundary} the authors prove the exact controllability towards a target state together with explicit rates of convergence.
Instead of a ring road they consider a bounded domain and apply the control at the entrance and exit point of the road.
Recently, \cite{friedrich2022lyapunov} established Lyapunov stabilization of a second order traffic model on the micro- and macroscopic scale by controlling the leading vehicle.

Here, we consider as in \cite{bayen2021boundary,huang2020stability} the nonlocal LWR model of \cite{BlandinGoatin2016} for general velocity functions and want to obtain a Lyapunov function and an exponential decay rate.
In particular, in contrast to \cite{bayen2021boundary,huang2020stability}, the Lyapunov function considers the $L^2$ distance in the nonlocal velocity and we apply, as in \cite{friedrich2022lyapunov}, the control on the leading vehicle.

The upcoming work is structured as follows: Section \ref{sec:models} introduces the traffic flow model of \cite{BlandinGoatin2016}. Section \ref{sec:stable} contains our main result with the explicit rate on the decay rate. In the last section we present numerical examples which demonstrate the theoretical results.
Further, the numerical results probably hold under less restrictive conditions than we assume in the theory and can also be obtained for the microscopic scale.

\section{A nonlocal LWR traffic flow model}\label{sec:models}
We base our analysis of the Lyapunov stabilization on the following first order nonlocal traffic flow model:
\begin{align}
\label{eq:ND}
&\partial_t \rho(t,x)+\partial_x \left(\rho V(t,x)\right)=0,\qquad (t,x)\in \mathbb{R}_+\times \mathbb{R}\\
&V(t,x):=v\left(W_{\ndt} \ast \rho\right)=v\left(\int_x^{x+{\eta}} W_{\ndt}(y-x)\rho(t,y)dy\right).\label{eq:ND:conv}
\end{align}
Here, $\rho$ is the traffic density, $v$ a suitable velocity function, ${\eta}>0$ the nonlocal reach and $W_{\ndt}$ a suitable kernel function.

This model was introduced in \cite{BlandinGoatin2016} and the well-posedness results were further extended in \cite{chiarello2018global,KeimerPflug2017}.
In particular, \cite{KeimerPflug2017} shows that no entropy condition is necessary as weak solutions are already unique.

The model \eqref{eq:ND:conv} is accompanied by initial conditions $\rho_0\in \operatorname{BV}(\mathbb{R};[0,\rho_{\max}])$ and we impose the following assumptions on the kernel function and the velocity function:
\begin{assumption}\label{ass:vprime}
Let $v\in C^2([0,\rho_{\max}];\mathbb{R}_+)$ with $v'<0$ and in addition we denote the upper bound on the derivative of $v$ by $v'_{\max}$, i.e.
        \[ v'(\rho)\leq v'_{\max}<0\qquad \forall \ \rho\in[0,\rho_{\max}].\]
\end{assumption}
\begin{assumption}\label{ass:kernelcons}
We set the kernel function to
\begin{equation}
W_{\ndt}(x)=\frac{1}{{\eta}}.
\end{equation}
\end{assumption}
A constant convolution kernel can model, e.g., connected autonomous vehicles which have the same degree of accuracy on information about the downstream traffic, independent of the distance.
\begin{remark}
For the theoretical result we need to restrict ourselves to constant kernel functions.
In general, the kernel needs to fulfill
\begin{align} \label{eq:ass:kernel}
W_{\ndt} \in C^1([0,{\eta}];\mathbb{R}^+)\  & \text{with} \ W_{\ndt}'\leq 0\ \text{and}\ \int_0^{\eta} W_{\ndt}(x) \dx=1 \ \forall {\eta}>0
\end{align}
to guarantee existence and uniqueness, see \cite{BlandinGoatin2016, chiarello2018global}.
\end{remark}
In order to control the velocity to a stationary solution, i.e. $V(t,x)=\bar v$, for $t$ sufficiently large, we apply specific initial data:
\begin{assumption}\label{ass:initcond}
We consider strictly positive initial data $\rho_0\in \operatorname{BV}(\mathbb{R};(0,\rho_{\max}])$ with
\begin{align}\label{eq:initcondcontrol}
\rho_0(x)=\bar \rho \qquad\text{for }x \geq b,
\end{align}
 $b\in\mathbb{R}$ and  $\bar \rho=v^{-1}(\bar v)$ for a given equilibrium velocity $\bar v$.
\end{assumption}
Note that $\bar \rho$ needs to be strictly positive such that $\bar v\in [0,v(0))$ holds.

As discussed in \cite[Remark 4.5]{KeimerPflug2017} the well-posedness of \eqref{eq:ND} with initial data $L^\infty(\mathbb{R})\ni \rho_0 \not\in L^1(\mathbb{R})$ is guaranteed.
In particular, the solution satisfies
\begin{align}\label{eq:maxprinciple}
0<\rho_{\min}:=\inf_{x\in\mathbb{R}}\rho_0(x)\leq \rho(t,x)\leq \sup_{x\in\mathbb{R}}\rho_0(x)\leq \rho_{\max},
\end{align}
see \cite{BlandinGoatin2016,chiarello2018global}, and hence we can set $v_{\max}'=v'(\rho_{\min})$.

\section{Lyapunov stabilization}\label{sec:stable}
The aim of this work is to control the nonlocal velocity $V(t,x)$ on the road towards a steady state $\bar v$.

In the following, we view the problem \eqref{eq:ND} with initial conditions fulfilling assumption \ref{ass:initcond} as an initial boundary value problem on the semi-infinite interval $(-\infty,\beta(t)]$ with $\beta(t):=b+\bar v t$ and we set the density for $x\geq\beta(t)$ to $\bar \rho$.
This induces the same waves as the initial value problem fulfilling assumption \ref{ass:initcond}, since for $x\geq\beta(t)$ the solution stays constant and $\beta(t)$ is the characteristic curve originating in $(0,b)$ in the $x$-$t$-plane, see \cite{KeimerPflug2017} for further details on the characteristics.
This allows to view the left boundary as a leading vehicle which controls the traffic and moves with speed $\bar v$, see \cite{friedrich2022lyapunov} for further details.
This leading vehicle has a direct influence on the interval of length ${\eta}$ behind it, i.e. $[\beta(t)-{\eta},\beta(t)]$.
Hence, our main result concerns this interval:
\begin{proposition}\label{prop:macroLyapunovDensity}
Let $0\leq \bar v<v(0)$ and the assumptions \ref{ass:vprime}, \ref{ass:kernelcons} and \ref{ass:initcond} hold.
We define
\begin{equation}\label{eq:macroLyapunov}
    L(t):=\int_{\beta(t)-{\eta}}^{\beta(t)}(V(t,x)-\bar v)^2 dx.
\end{equation}
Then, we obtain the following bound
\begin{align*}
    L(t)\leq L(0) \exp\left(\frac{2}{{\eta}} v'_{\max} \rho_{\min}t\right).
\end{align*}
\end{proposition}
\begin{proof}
First note that the boundary condition for $\beta(t)$ simplifies the nonlocal term in $\eqref{eq:ND:conv}$ and its partial derivatives for $x\in [\beta(t)-{\eta},\beta(t)]$ and a constant kernel as in the assumption \ref{ass:kernelcons}.
For simplicity, we define:
\begin{align*}
R(t,x):=\frac{1}{{\eta}}\left( \int_x^{\beta(t)} \rho(t,y)\dy+(x+{\eta}-\beta(t))\bar \rho\right).
\end{align*}
Next, we have
\begin{align}
    V(t,x)&=v\left(R(t,x)\right)\notag\\
    \partial_x V(t,x)&= \frac{1}{{\eta}}\left(\bar \rho- \rho(t,x) \right)v'\left(R(t,x)\right)\label{eq:partialxV}
\end{align}
and hence, a direct relationship between the local density and the equilibrium density is obtained.
This derivative will play a key role for proving our stabilization result.
Next, we exploit the temporal derivative of the nonlocal term:
\begin{align}\notag
 \partial_t V(t,x)&=\frac1{\eta}\left(\bar v \bar \rho+\int_x^{\beta(t)}\partial_t \rho(t,y)\dy-\bar v \bar \rho\right)v'\left(R(t,x)\right)\\
\notag &=-\frac{1}{{\eta}}\int_x^{\beta(t)}\partial_x \left(\rho(t,y)V(t,y)\right)\dy v'\left(R(t,x)\right)\\
\label{eq:partialtV} &=\frac{1}{{\eta}}\left(\rho(t,x)V(t,x)-\bar v\bar\rho \right) v'\left(R(t,x)\right).
\end{align}
Here, we used $\partial_t \rho(t,x)=-\partial_x (\rho(t,x)V(t,x))$.
Since $\rho$ is bounded and the integral evaluations are on a finite interval, we obtain that $R,\ V,\ \partial_x V,\ \partial_t V$ are all bounded.
In particular, this ensures that the integrals we consider in the following and the Lyapunov function $L(t)$ are bounded, too.

Turning now to the Lyapunov function \eqref{eq:macroLyapunov}, we directly obtain
\begin{align*}
    \frac{d}{dt}L(t)=&\bar v (V(t,\beta(t))-\bar v)^2 - \bar v(V(t,\beta(t)-{\eta})-\bar v)^2\\
    &+2\int_{\beta(t)-{\eta}}^{\beta(t)}(V(t,x)-\bar v)\partial_t V(t,x) dx.
\end{align*}
Due to the boundary conditions, we have $V(t,\beta(t))=\bar v$ and the first term equals zero.
Next, by using \eqref{eq:partialtV} and adding a zero we obtain
\begin{align*}
    \frac{d}{dt}L(t)=& - \bar v(V(t,\beta(t)-{\eta})-\bar v)^2+2\int_{\beta(t)-{\eta}}^{\beta(t)}\frac{ v'\left(R(t,x)\right)}{{\eta}}\rho(t,x)(V(t,x)-\bar v)^2 dx\\
    &+2\int_{\beta(t)-{\eta}}^{\beta(t)}\frac{ v'\left(R(t,x)\right)}{{\eta}}\bar v \left(\rho(t,x)-\bar\rho \right)(V(t,x)-\bar v)  dx.
\end{align*}
Using \eqref{eq:partialxV} we can rewrite the last term as
\begin{align*}
 &2\int_{\beta(t)-{\eta}}^{\beta(t)}\frac{ v'\left(R(t,x)\right)}{{\eta}}\bar v \left(\rho(t,x)-\bar\rho \right)(V(t,x)-\bar v)  dx \\
& =-\bar v \int_{\beta(t)-{\eta}}^{\beta(t)}\partial_x (V(t,x)-\bar v)^2  dx   =\bar v\left(V(t,\beta(t)-{\eta})-\bar v\right)^2.
\end{align*}
Hence, we obtain
\begin{align*}
        \frac{d}{dt}L(t)=2\int_{\beta(t)-{\eta}}^{\beta(t)}\frac{ v'\left(R(t,x)\right)}{{\eta}}\rho(t,x)(V(t,x)-\bar v)^2 dx
        \leq \frac{2}{{\eta}}\rho_{\min}v_{\max}'L(t).
\end{align*}
The latter inequality holds due to assumption \ref{ass:vprime} and the maximum principle for $\rho$ given by \eqref{eq:maxprinciple}.

Applying Gr\"onwall's lemma yields the result.
\end{proof}

\begin{remark}
In the case of a linear velocity function a stabilization result concerning the $L^2$ distance of the density and the stationary density is given by \cite[Proposition 3.8]{friedrich2022lyapunov}.
In this case, the Lyapunov function decreases only on a specific subinterval of $[\beta(t)-{\eta},\beta(t)]$.
\end{remark}
\begin{remark}
In order to estimate $\rho(t,x)$ from below and obtain a decreasing Lyapunov function, we need to ensure $\rho_{\min}>0$.
Hence, we cannot prove the exponential stability for $\bar v=v(0)$.
\end{remark}
As already mentioned in the introduction, \cite{huang2020stability} considers the stability of \eqref{eq:ND:conv} on a ring road.
In addition, they restrict themselves to a linear velocity function.
One of the main results in \cite{huang2020stability} is that for a constant kernel as in assumption \ref{ass:kernelcons}, specific initial data and a specific nonlocal reach ${\eta}$ traveling waves can be created, such that no control is possible.
Nevertheless, this is no contradiction to our results as we do not consider a ring road and such phenomena do not occur.

\section{Numerical simulations}
\subsection{Numerical scheme}
In order to illustrate the theoretical results we will use the numerical scheme presented in \cite{chiarello2019multiclass,friedrich2018godunov}.
For the spatial step size $\Delta x>0$ the numerical flux is given by
\[F_{j+\frac12}^n=V_j^n \rho_j^n,\]
where
\[V_j^n = v\left(\sum_{k=0}^{\lfloor{\eta}/\Delta x\rfloor-1}\gamma_k\rho^n_{j+k+1}\right)\quad\text{ and }\quad\gamma_k=\int_{k\Delta x}^{(k+1)\Delta x} W_{\ndt}(y-x)dy.\]
The complete scheme is then given by
\begin{equation}\label{eq:macroscheme}
     \rho_j^{n+1}= \rho_j^{n}-\frac{\Delta t^n}{\Delta x}\left(F^n_{j+\frac12}-F^n_{j-\frac12}\right).
\end{equation}
At the right boundary we apply ghost cells equal to $\bar \rho$.
As CFL condition we choose an adaptive step size control determined by the maximum nonlocal velocity, i.e. $\Delta t^n\leq \Delta x/\max_j V_j^n $.
\subsection{Numerical example}
As a numerical example we consider the situation of a bulk of cars standing bumper-to-bumper and starting to move.
For example such situations occur when a red light turns green or a completely stopped traffic jam begins to resolve.
This leads to the initial condition $\rho_0(x)=1$ for $x\leq 0=b$.
The speed law is given by $v(\rho)=1-\rho$.
Without a control the leading vehicle would move at maximum speed $v(0)=1$.
Instead we want to stabilize the traffic towards the velocity $\bar v=0.5$ (and $\bar \rho=0.5$).
In the numerical simulations the space step size is given by $\Delta x=5\cdot 10^{-3}$.
The constant ${\eta}$ is ${\eta}=1$.
We compare the numerical stabilization for the following kernels
\begin{align*}
\wt^{\text{const.}}(x)=\frac{1}{{\eta}}, \qquad \wt^{\text{lin.}}(x)=\frac{2({\eta}-x)}{{\eta}^2}, \qquad \wt^{\text{conc.}}(x)=\frac{3({\eta}^2-x^2)}{2{\eta}^3}.
\end{align*}
The theoretical results only ensure the convergence of the Lyapunov function \eqref{eq:macroLyapunov} for the kernel $\wt^{\text{const.}}$.
Nevertheless, the kernels $\wt^{\text{lin.}}$ and $\wt^{\text{conc.}}$ fulfill \eqref{eq:ass:kernel} such that the well-posedness of the model \eqref{eq:ND} is satisfied.
Figure \ref{fig:lyapunov} displays the logarithm of the Lyapunov function \eqref{eq:macroLyapunov} for all kernels together with the theoretical upper bound of Proposition \ref{prop:macroLyapunovDensity}.
As can be seen, even though not covered by our theoretical results, the bound is valid for all considered kernels.

\begin{figure}
\setlength{\fwidth}{0.5\textwidth}
\centering
\input{Lyapunov_diffKernels}
\captionsetup{margin=4cc}
\caption{Logarithm of the Lyapunov function \eqref{eq:macroLyapunov} over time together with the exponential upper bound for different kernel functions.}
\label{fig:lyapunov}
\end{figure}
\begin{remark}
In order to prove the Lyapunov stabilization on the complete interval $[\beta(t)-{\eta},\beta(t)]$ we cannot consider the density, i.e.
\begin{equation}
\label{eq:lyapunovdensity}
\tilde{L}(t)=\int_{\beta(t)-{\eta}}^{\beta(t)}(\rho(t,x)-\bar\rho)^2 dx.
\end{equation}
A numerical example with initial data
\[\rho_0(x)=\begin{cases}
0.01,\quad&\text{for }x\leq-0.5,\\
0.35,\quad&\text{for }x\in(-0.5,0],\\
0.5,\quad&\text{for }x>0,
\end{cases}\]
and the rest of the parameters as before can be seen in Figure \ref{fig:lyapunovdensity}.
In particular, the numerical Lyapunov function \eqref{eq:lyapunovdensity} has an increasing part which demonstrates that a stabilization result for \eqref{eq:lyapunovdensity} cannot be expected.
In order to get a decreasing Lyapunov function one needs to consider a specific subinterval of $[\beta(t)-{\eta},\beta(t)]$ as in \cite{friedrich2022lyapunov}.
\begin{figure}
\setlength{\fwidth}{0.5\textwidth}
\centering
\input{Lyapunov_density}
\captionsetup{margin=4cc}
\caption{Logarithm of the Lyapunov function \eqref{eq:lyapunovdensity} over time with the kernel $\wt^{\text{const.}}$.}
\label{fig:lyapunovdensity}
\end{figure}
\end{remark}
\begin{remark}
In \cite{friedrich2022lyapunov} the Lyapunov stabilization of a corresponding microscopic model is derived with similar results as on the macroscopic scale.
For the model \eqref{eq:ND:conv} a microscopic formulation is considered in, e.g. \cite{ridder2019traveling}.
A numerical example of the situation above with $\rho_0(x)=1$ for $x\leq 0$ which considers the corresponding formulation of the Lyapunov function \eqref{eq:macroLyapunov} on the microscopic scale with $\wt^{\text{const.}}$ can be seen in Figure \ref{fig:lyapunovmicro}.
The latter figure demonstrates that a stabilization to the equilibrium velocity can be obtained.
Nevertheless, for an analytic treatment one has to take care of discontinuities which can be seen in Figure \ref{fig:lyapunovmicro}.
These arise due to cars entering or leaving the domain $[\beta(t),\beta(t)-{\eta}]$.
\begin{figure}
\setlength{\fwidth}{0.5\textwidth}
\centering
\input{Lyapunov_discmicro}
\captionsetup{margin=4cc}
\caption{Logarithm of a microscopic equivalent to the Lyapunov function \eqref{eq:macroLyapunov} over time together with the corresponding exponential upper bound.}
\label{fig:lyapunovmicro}
\end{figure}
\end{remark}

\section{Conclusion}
In this work we have presented a Lyapunov function and the explicit rate such that the nonlocal velocity of a LWR-type traffic model on a single road tends to its equilibrium velocity.
For the theoretical analysis, we had to restrict ourselves to a constant kernel function.
Nevertheless, numerical examples suggest that the asymptotic stabilization effect can be obtained for general kernels and the corresponding microscopic model.
Future work may include extending the results to those cases.

Furthermore, by using the stabilization in the nonlocal velocity we are able to obtain our result on an interval of length ${\eta}$ behind the leading vehicle.
This is more intuitive than the stabilization result of \cite{friedrich2022lyapunov} which needs a specific subinterval.
Hence, future work may additionally include studying the Lyapunov stabilization of the nonlocal velocity for the second order traffic flow model considered in \cite{friedrich2020micromacro, friedrich2022lyapunov}.

\section*{Acknowledgement}
J.F. was supported by the German Research Foundation (DFG) under grant HE 5386/18-1, 19-2, 22-1, 23-1.

\vspace{\baselineskip}
\bibliographystyle{plain}
\providecommand{\othercit}{}
\providecommand{\jr}[1]{#1}
\providecommand{\etal}{~et~al.}


\section{Introduction and basic notation} \label{secc01}
Robust optimization arises as a deterministic approach when addressing an optimization problem under uncertainty data. This paper revisites
the following robust optimization problem:
\begin{equation}\label{prob:00}
\min\Big\{\dfrac{1}{2}x^\top Ax+a^\top x :~\alpha\leq
\dfrac{1}{2}x^\top Bx+b^\top x\leq\beta,~\forall~(B,b)\in{\mathcal B}_b\Big\},
\end{equation}
where $\mathcal{B}_b\doteq\{B_1+\mu B_2:\mu\in[\mu_1,\mu_2]\}\times\{b_1+\delta b_2:\delta\in[\delta_1,\delta_2]\}$, with all the matrices being real symmetric, $a,b\in\mathbb{R}^n$ and
$\alpha,\beta,\delta_1,\delta_2,\mu_1,\mu_2$ are given real numbers.
Optimization problems that can be  modeled by quadratic functions appear, for instance, in \cite{Beck,mat,slemma,Polyak,sid}.

Problem \eqref{prob:00} includes that examined in \cite{jeya2013}:
\begin{equation}\label{prob:01}
\min\Big\{\dfrac{1}{2}x^\top Ax+a^\top x :~
\dfrac{1}{2}x^\top Bx+b^\top x\leq\beta,~\forall~(B,b)\in{\mathcal B}_b\Big\}.
\end{equation}
Theorem 5.1 in \cite{jeya2013} provides a characterization of robust optimality for
problem \eqref{prob:01} under the convexity of the set
 $$
\Big\{(x^\top H_0x,x^\top H_{1}x,x^\top H_{2}x):x\in\mathbb{R}^{n+1}\Big\},
$$
where
$$
H_0=\begin{pmatrix}A & a\\
a^\top  &  2\gamma \end{pmatrix},
~
H_{1}=\begin{pmatrix}B_1+\mu_1 B_2 & b_1+\delta_1 b_2\\
(b_1+\delta_1b_2)^\top  &  -2\beta \end{pmatrix},~
H_{2}=\begin{pmatrix}B_1+\mu_2 B_2 & b_1+\delta_2 b_2\\
(b_1+\delta_2b_2)^\top  &  -2\beta \end{pmatrix}
$$
with $\gamma=-f(\overline{x})$. Here, $\overline{x}$ is a feasible point of the robust optimization problem \eqref{prob:01}, which is either to be supposed becoming an optimal solution, or to be optimal for deriving optimality conditions.

Certainly the presence of the matrices $H_0, H_1, H_2$ is because the authors in \cite{jeya2013}  homogenize problem \eqref{prob:01} in order to apply  the Dines convexity theorem (\cite{Dines})  valid for quadratic forms. Finally, we realize there is a gap in the proof of Theorem 5.1 in \cite{jeya2013}, but we were unable to find a counterexample to such a result under their assumptions. This is discussed in detail after Theorem \ref{minimocorregido} in Section \ref{secc:uni}. Notice that the approach employed in  \cite{jeya2013} was also applied in \cite{nuevo}.

We have to point out that problem \eqref{prob:00} (and so problem \eqref{prob:01}) was studied without homogenizing the problem thanks to the convexity result established in \cite[Theorem 4.19]{Flores} valid for inhomogeneous quadratic functions. This  allows us to impose the convexity of
a set being the image of $\mathbb{R}^n$ via the  inhomogeneous quadratic functions.

Associated to problem \eqref{prob:00}, purposes of the present paper are to establish an alternative robust result (Theorem \ref{teoralt}), a  robust S-lemma (Theorem \ref{slema}) and a characterization of robust optimality (Theorem \ref{minimonohomoalt}), for problem \eqref{prob:00}. Finally, we provide a counterexample (Example \ref{ex100}) to the argument employed in the proof of Theorem 5.1 in \cite{jeya2013} related to problem \eqref{prob:01}.

Thus, the structure of the present paper is as follows. Section \ref{sect02} establishes the convexity of images for quadratic mappings by applying the Ramana-Goldman criterion \cite{rg1995} (see also \cite[Theorem 2.1 ]{Flores2}). The main results are presented in  Section \ref{sect003}, and Section \ref{secc:uni} revisites problem
\eqref{prob:01} discussed in \cite{jeya2013}.

\section{A preliminary result: convexity of images}\label{sect02}

By ${\mathcal S}^n$ we denote the set of symmetric matrices of order $n\in\mathbb{N}$ with real entries; ${\mathcal S}_+^n$ denotes the subset of
${\mathcal S}^n$ whose elements are positive semidefinite matrices, and we write
$A\succeq 0$ if $A\in {\mathcal S}_+^n$; and ${\mathcal S}_{++}^n$ stands for the matrices in ${\mathcal S}^n$ that are positive definite, and in this case we write
$A\succ 0$ if $A\in {\mathcal S}_{++}^n$.

It is our purpose to prove the convexity of images for quadratic mappings under the Ramana-Goldman criterion \cite{rg1995} (see also \cite[Theorem 2.1 ]{Flores2}). To that end, we are
given $A_i\in{\mathcal S}^n$, $b_i\in\mathbb{R}^n$, $c_i\in\mathbb{R}$ for $i=0, 1,\ldots,m$, we set
$$
M_i=\begin{pmatrix}A_i & b_i\\
b_i^\top  &  2c_i\end{pmatrix}, f_i(x)= x^\top A_i x + 2b_i^\top x,\; \overline{f}_i(x)= x^\top A_i x,~x\in\mathbb{R}^n.
$$
Furthermore, let us consider the function
$$
G(x,t)= (g_0(x,t),  g_1(x,t),  \ldots,  g_m(x,t)), ~~(x,t)\in\mathbb{R}^n\times\mathbb{R},
$$
where $g_i$ is defined by
$g_i(x,t)= \begin{pmatrix} x\\t \end{pmatrix}^\top M_i\begin{pmatrix} x\\t \end{pmatrix} $.

\begin{lemma} \label{convex_R_5} Let $A_i\in {\mathcal S}^n$,  $c_i \in \mathbb{R}$,
$b_i\in\mathbb{R}^n$,  $i= 0, 1, \dots, m$ be as above. Set
$$
F(x)= (f_0(x),  f_1(x),  \ldots,  f_m(x)),~~\overline{F}(x)= (\overline{f}_0(x),  \overline{f}_{1}(x), \ldots,  \overline{f}_{m}(x)).
$$
If  $F(\mathbb{R}^n)$ and  $\overline{F}(\mathbb{R}^n)$ are convex then  $G(\mathbb{R}^{n+1})$ is convex.
\end{lemma}
\begin{proof} Set
	 $\Lambda:= F(\mathbb{R}^n)$ and  $\overline{\Lambda}:= \overline{F}(\mathbb{R}^n)$ and $\Omega:= G(\mathbb{R}^{n+1})$.

	As $\Lambda$ is convex, by the convexity criterion due to Ramana-Goldman  (see also \cite[Theorem 2.1 ]{Flores2}),
	$\Lambda+\overline{\Lambda}= \Lambda$.  We easily get that for $(x,t) \in \mathbb{R}^{n+1}$,
	\[
	\begin{pmatrix} x\\t \end{pmatrix}^\top M_i\begin{pmatrix} x\\t \end{pmatrix} = x^\top A_ix+2tb_i^{\top} x+2c_i t^2,~~i= 0,1, \dots, m.
	\]

	By setting $\overline{\gamma} = 2(c_0, c_1, \dots,  c_m)$, we obtain
	\begin{equation}\label{contain}
	\Lambda +\overline{\gamma} \subseteq \Omega.
	\end{equation}
	Let ${z}_1 = G(x_1, t_1)$,  ${z}_2 = G(x_2, t_2)$ be any elements in $\Omega$ and let $\lambda \in ~]0,1[$. We distinguish three cases.  \\
	$(i)$: $t_1 \neq 0$ and $t_2 \neq 0$. Then,
	$
	{z}_1= t_1^2(F(x_1/t_1)+ \overline{\gamma})
	$
	and
	$
	{z}_2= t_2^2(F(x_2/t_2)+ \overline{\gamma}).
	$\\
	The convexity of  $\Lambda$ implies that

	\begin{align*}
	\frac{\lambda}{\lambda t_1^2 +(1-\lambda)t_2^2} &{z}_1+ \frac{1-\lambda}{\lambda t_1^2 +(1-\lambda) t_2^2} {z}_2 \\
	& =
	\frac{\lambda t_1^2}{\lambda t_1^2 +(1-\lambda)t_2^2}F(x_1/t_1)+\frac{(1-\lambda)t_2^2}{\lambda t_1^2 +(1-\lambda)t_2^2}F(x_2/t_2) + \overline{\gamma} \in  \Lambda + \overline{\gamma}.
	\end{align*}
	Taking into account \eqref{contain} and the fact the $\Omega$ is a cone,  we obtain
	\[
	\lambda {z}_1+ (1-\lambda){z}_2 \in \mathbb{R}_{++}( \Lambda +\overline{\gamma})  \subseteq \Omega.
	\]
	$(ii)$: $t_1=t_2= 0$.  Then,  ${z}_1, {z}_2 \in \overline{\Lambda}$, and because of the convexity of $\overline{\Lambda}$,  we get
	\[
	\lambda \left({z}_1+(1-\lambda){z}_2\right)  \in \overline{\Lambda} \subseteq \Omega.
	\]
	$(iii)$: $t_1\neq 0$ and $t_2=0$. Then,  since $\overline{\Lambda}$ is a cone,
	\[
	\lambda {z}_1+(1-\lambda) {z}_2 \in   \lambda t_1^2( \Lambda +\overline{\gamma}) +  \lambda t_1^2 \overline{ \Lambda}  \subseteq \lambda t_1^2( \Lambda +\overline{\gamma}) \subseteq \Omega.
	\]
	This completes the proof  that $\Omega$ is convex.
	\qed
\end{proof}

Part $(a)$ of the following result is exactly Theorem 2.3 $(i)$ in \cite{Beck}, and $(b)$ is a consequence of the previous lemma.

\begin{corollary} \label{cor_R_5}
	Let the same  hypotheses  of Lemma \ref{convex_R_5} be satisfied.  Let  $\rho_i\in \mathbb{R}$,  for $i=1, \dots, m$.
	If $n\geq m+1$,  $A_0\in {\mathcal S}_{++}^n$,  $A_i= \rho_i A_0$ for $i=1, \dots, m$,  then
	\begin{itemize}
		\item[$(a)$] $F(\mathbb{R}^n)$ and  $\overline{F}(\mathbb{R}^n)$ are convex.
		\item[$(b)$] $G(\mathbb{R}^{n+1})$ is convex.
	\end{itemize}
\end{corollary}
\begin{proof} By assumption on $A_i$, we can apply \cite[Theorem 2.3 (i)]{Beck} to obtain the convexity of $F(\mathbb{R}^n)$ and  $\overline{F}(\mathbb{R}^n)$. Then, $(b)$ follows from Lemma \ref{convex_R_5}.
\qed
\end{proof}

\section{The main results}\label{sect003}

Denote the function:
$$
f(x)\doteq\dfrac{1}{2}x^\top Ax+a^\top x
$$
and let us define the following  matrices in ${\mathcal S}^{n+1}$:

$$
H_0\doteq\begin{pmatrix}A & a\\
a^\top  &  2\gamma \end{pmatrix},~
W(\delta,\lambda)\doteq\begin{pmatrix}B_1 & b_1+\delta b_2\\
(b_1+\delta b_2)^\top  &  -2\lambda \end{pmatrix},~
W_2\doteq\begin{pmatrix}B_2 & 0\\
0^\top  &  0 \end{pmatrix},
$$
and set
$$
W_{1_{\beta}}=W(\delta_1,\beta);~
W_{2_{\beta}}=W(\delta_2,\beta);~
W_{1_{\alpha}}=W(\delta_1,\alpha);~
W_{2_{\alpha}}=W(\delta_2,\alpha).
$$

The following set will play an important role in the following.
$$
\Omega_W\doteq\Big\{\Big(\dfrac{1}{2}y^\top H_0y
,\max_{\mu\in[\mu_1,\mu_2]}\dfrac{1}{2}y^\top (W_{1_{\beta}}+\mu W_2)y),\max_{\mu\in[\mu_1,\mu_2]}\dfrac{1}{2}y^\top (W_{2_{\beta}}+\mu W_2)y,
$$
\begin{equation
-\min_{\mu\in[\mu_1,\mu_2]}\dfrac{1}{2}y^\top (W_{1_{\alpha}}+\mu W_2)y,-\min_{\mu\in[\mu_1,\mu_2]}\dfrac{1}{2}y^\top (W_{2_{\alpha}}+\mu W_2)y\Big):
y\in\mathbb{R}^{n+1}\Big\}+\rm{int} \ \mathbb{R}^5_+.
\end{equation}
By Corollary 1 in \cite{fbp2022}, $\Omega_W$ is convex if the set
$$\Omega_\mu\doteq
\Big\{\Big(y^\top H_0y,y^\top(W_{1_{\beta}}+\mu_1 W_2)y,
y^\top(W_{2_{\beta}}+\mu_1 W_2)y, y^\top (W_{1_{\beta}}+\mu_2 W_2)y,
$$
$$
y^\top (W_{2_{\beta}}+\mu_2 W_2)y,
-y^\top(W_{1_{\alpha}}+\mu_1 W_2)y,
-y^\top(W_{2_{\alpha}}+\mu_1 W_2)y, -y^\top (W_{1_{\alpha}}+\mu_2 W_2)y,
$$
$$
-y^\top (W_{2_{\alpha}}+\mu_2 W_2)y\Big):y\in\mathbb{R}^{n+1} \Big\}+\rm{int} \ \mathbb{R}^9_+
$$
is so.

\begin{theorem} $($A robust alternative result$)$\label{teoralt}
 Let $A, B_1,B_2\in {\mathcal S}^n$, $a, b_1, b_2\in\mathbb{R}^n$ and $\gamma, \alpha, \beta, \mu_1, \mu_2,\delta_1,\delta_2\in\mathbb{R}$, with $\mu_1\leq\mu_2$, $\delta_1\leq\delta_2$ and $\alpha<\beta$. Assume that $\Omega_W$ is convex. Then, exactly one of the two following assertions hold:
	\begin{itemize}
		\item[$(a)$] $\exists~ x\in\mathbb{R}^n: \dfrac{1}{2}x^\top Ax+a^\top x + \gamma<0$, $\alpha<\dfrac{1}{2}x^\top (B_1+\mu B_2)x+(b_1+\delta b_2)^\top x<\beta$, $\forall~\mu\in [\mu_1,\mu_2]$, $\forall~\delta\in [\delta_1,\delta_2]$.
		\item[$(b)$] $\exists~(\lambda_0,\lambda_1,\lambda_2)\in\mathbb{R}^3_+\backslash\{0\}, \exists~\mu_{\alpha},\mu_{\beta}\in [\mu_1,\mu_2 ], \exists~\delta_{\alpha},\delta_{\beta}\in [\delta_1,\delta_2 ]:~ \forall~ x\in \mathbb{R}^n$
		$$
		\lambda_0\Big(\dfrac{1}{2}x^\top Ax+a^\top x+\gamma\Big) +\lambda_1\Big(	\dfrac{1}{2}x^\top (B_1+\mu_{\beta} B_2)x+(b_1+\delta_{\beta} b_2)^\top x-\beta\Big)+
		$$
		$$ \lambda_2\Big(\alpha-\Big(\dfrac{1}{2}x^\top (B_1+\mu_{\alpha} B_2)x+(b_1+\delta_{\alpha} b_2)^\top x \Big)\Big)\geq 0,
		$$
	\end{itemize}
	where $\mu_{\alpha}+\mu_{\beta}=\mu_1+\mu_2$.\\
	\\
	In addition, we observe that $(b)$ may be written equivalently
	as
	$$
	\lambda_0 A+\lambda_1(B_1+\mu_{\beta}B_2)-\lambda_2(B_1+\mu_{\alpha}B_2)\succeq 0 \ \rm{and}$$
	$$ \exists~ \overline{x}\in\mathbb{R}^n: \Big( \lambda_0A+\lambda_1(B_1+\mu_{\beta}B_2)-\lambda_2(B_1+\mu_{\alpha}B_2)\Big)\overline{x}+\lambda_0 a + \lambda_1(b_1+\delta_{\beta}b_2)+\lambda_2(b_1+\delta_{\alpha} b_2)=0.$$
\end{theorem}
\begin{proof} It is obvious that both statements $(a)$ and
	$(b)$ cannot be fulfilled simultaneously. Thus, we must
	check that if $(a)$ does not hold, $(b)$ does.\\
	\\
	\textbf{Step 1: The homogenization system.} If $(a)$ does not hold, then there exists no  $x\in\mathbb{R}^n$ such that for all $\mu\in[\mu_1,\mu_2]$ and all $\delta\in[\delta_1,\delta_2]$
	$$
	\dfrac{1}{2}x^\top Ax+a^\top x + \gamma<0, \quad \dfrac{1}{2}x^\top (B_1+\mu B_2)x+(b_1+\delta b_2)^\top x<\beta,
	$$
	$$
	-\left(\dfrac{1}{2}x^\top (B_1+\mu B_2)x+(b_1+\delta b_2)^\top x\right)<-\alpha.
	$$
	By setting $\mathcal{B}_b\doteq\{B_1+\mu B_2:\mu\in[\mu_1,\mu_2]\}\times\{b_1+\delta b_2:\delta\in[\delta_1,\delta_2]\}$, the previous is equivalent to the nonexistence of $x\in\mathbb{R}^n$ such that
	$$
	\dfrac{1}{2}x^\top Ax+a^\top x + \gamma<0,~~	\max\left\{\dfrac{1}{2}x^\top Bx+b^\top x-\beta: (B,b)\in\mathcal{B}_b\right\}<0,
	$$
	$$
	-\min\left\{\dfrac{1}{2}x^\top Bx+b^\top x-\alpha: (B,b)\in\mathcal{B}_b\right\}<0.
	$$
	We claim that the following homogeneous system in $\mathbb{R}^{n+1}$:

	\begin{equation}\label{sis:hom1}
	\dfrac{1}{2}x^\top Ax+ta^\top x + t^2\gamma<0,~ \max\left\{\dfrac{1}{2}x^\top Bx+tb^\top x-t^2\beta:(B,b)\in\mathcal{B}_b\right\}<0,
	\end{equation}
	\begin{equation}\label{sis:hom2}
	-\min\left\{\dfrac{1}{2}x^\top Bx+tb^\top x-t^2\alpha: (B,b)\in\mathcal{B}_b\right\}<0
	\end{equation}
	has no solution. If, on the contrary, there was  a solution $(\overline{x},\overline{t})\in\mathbb{R}^{n+1}$ such that
	$$
	\dfrac{1}{2}\overline{x}^\top A\overline{x}+\overline{t}a^\top \overline{x} + \overline{t}^2\gamma<0, \quad \max\left\{\dfrac{1}{2}\overline{x}^\top B\overline{x}+\overline{t}b^\top \overline{x}-\overline{t}^2\beta:(B,b)\in\mathcal{B}_b\right\}<0,
	$$
	$$
-\min\left\{\dfrac{1}{2}\overline{x}^\top B\overline{x}+\overline{t}b^\top \overline{x}-\overline{t}^2\alpha: (B,b)\in\mathcal{B}_b\right\}<0,
	$$
	we immediately reach a contradiction in case $\overline{t}\neq 0$.

	So, suppose that $\overline{t}=0$. Then, the system \eqref{sis:hom1}-\eqref{sis:hom2} reduces to
	$$
	\dfrac{1}{2}\overline{x}^\top A\overline{x}<0, \quad	\max\left\{\dfrac{1}{2}\overline{x}^\top (B_1+\mu B_2)\overline{x}:\mu\in [\mu_1,\mu_2]\right\}<0,
	$$
	$$
	-\min\left\{\dfrac{1}{2}\overline{x}^\top (B_1+\mu B_2)\overline{x}: \mu\in [\mu_1,\mu_2]\right\}<0,
	$$
	which is impossible to hold. Thus, the claim is proved. \\
	On the other hand, observe that for every $(x,t)\in\mathbb{R}^n\times\mathbb{R}$, the minimum and maximum values in \eqref{sis:hom1}-\eqref{sis:hom2} are achieved in, at least, one of the extreme points of the rectangle $[\mu_1,\mu_2]\times[\delta_1,\delta_2]$, that is, in one of the elements  $(B_1+\mu_1B_2,b_1+\delta_1b_2)$, $(B_1+\mu_1B_2,b_1+\delta_2b_2)$, $(B_1+\mu_2B_2,b_1+\delta_1b_2)$ or  $(B_1+\mu_2B_2,b_1+\delta_2b_2)$. Then, the nonexistence of solution to the system \eqref{sis:hom1}-\eqref{sis:hom2} is equivalent to the nonexistence of $y\in\mathbb{R}^{n+1}$ solution to the system
	$$
	\dfrac{1}{2}y^\top H_0y<0
	$$
	$$
	\max\left\{\dfrac{1}{2} y^\top(W_{1_{\beta}}+\mu W_2)y
	:~\mu\in[\mu_1,\mu_2]\right\}<0,
	$$
	$$
	\max\left\{\dfrac{1}{2} y^\top(W_{2_{\beta}}+\mu W_2)y
	:~\mu\in[\mu_1,\mu_2]\right\}<0.
	$$
	$$
	-\min\left\{\dfrac{1}{2} y^\top(W_{1_{\alpha}}+\mu W_2)y:~\mu\in[\mu_1,\mu_2]\right\}
	<0,
	$$
	$$
	-\min\left\{\dfrac{1}{2} y^\top(W_{2_{\alpha}}+\mu W_2)y:~\mu\in[\mu_1,\mu_2]\right\}
	<0.
	$$

This means that	 $(0,0,0,0,0)\notin \Omega_W$.\\
\textbf{Step 2: A first use of a separation result.}
Since $\Omega_W$ is convex,
there exist  $(\lambda_0,\lambda_1,\lambda_2,\lambda_3,\lambda_4)\in\mathbb{R}^5_+\backslash\{0\}$ such that for all $y\in\mathbb{R}^{n+1}$
$$
\lambda_0 \left(\dfrac{1}{2}y^\top H_0y\right)+\lambda_1\left(\max_{\mu\in[\mu_1,\mu_2]}\dfrac{1}{2}y^\top (W_{1_{\beta}}+\mu W_2)y\right)+\lambda_2\left(\max_{\mu\in[\mu_1,\mu_2]}\dfrac{1}{2}y^\top (W_{2_{\beta}}+\mu W_2)y\right)+
$$
$$
\lambda_3\left(-\min_{\mu\in[\mu_1,\mu_2]}\dfrac{1}{2}y^\top (W_{1_{\alpha}}+\mu W_2)y\right)+\lambda_4\left(-\min_{\mu\in[\mu_1,\mu_2]}\dfrac{1}{2}y^\top (W_{2_{\alpha}}+\mu W_2)y\right)\geq 0.
$$
Since each of the minimum or maximum values are achieved in either $\mu_1$ or $\mu_2$, we obtain for all $y\in\mathbb{R}^{n+1}$,
$$
\lambda_0 \left(\dfrac{1}{2}y^\top H_0y\right)+\lambda_1\left(\max\left\{\dfrac{1}{2}y^\top (W_{1_{\beta}}+\mu_1 W_2)y, \dfrac{1}{2}y^\top (W_{1_{\beta}}+\mu_2 W_2)y\right\}\right)+
$$
$$
\lambda_2\left(\max\left\{\dfrac{1}{2}y^\top (W_{2_{\beta}}+\mu_1 W_2)y, \dfrac{1}{2}y^\top (W_{2_{\beta}}+\mu_2 W_2)y\right\}\right)+
$$
$$
\lambda_3\left(-\min\left\{\dfrac{1}{2}y^\top (W_{1_{\alpha}}+\mu_1 W_2)y, \dfrac{1}{2}y^\top (W_{1_{\alpha}}+\mu_2 W_2)y\right\}\right)+
$$
$$
\lambda_4\left(-\min\left\{\dfrac{1}{2}y^\top (W_{2_{\alpha}}+\mu_1 W_2)y, \dfrac{1}{2}y^\top (W_{2_{\alpha}}+\mu_2 W_2)y\right\}\right)\geq 0.
$$
Thus, there is no $y\in\mathbb{R}^{n+1}$ solution to the system:
$$
\dfrac{1}{2}y^\top \Big(\lambda_0 H_0+\lambda_1 (W_{1_{\beta}}+\mu_1 W_2)+\lambda_2 (W_{2_{\beta}}+\mu_1 W_2)-\lambda_3 (W_{1_{\beta}}+\mu_2 W_2)-\lambda_4 (W_{2_{\beta}}+\mu_2 W_2)\Big)y<0;
$$
$$
\dfrac{1}{2}y^\top \Big(\lambda_0 H_0+\lambda_1 (W_{1_{\beta}}+\mu_2 W_2)+\lambda_2 (W_{2_{\beta}}+\mu_2 W_2)-\lambda_3 (W_{1_{\beta}}+\mu_1 W_2)-\lambda_4 (W_{2_{\beta}}+\mu_1 W_2)\Big)y<0.
$$
This means that $(0,0)\not\in \Omega_2$, where
$$
\Omega_2\doteq
\Big\{\Big(y^\top (\lambda_0 H_0+\lambda_1
H_{1,1_{\beta}}+\lambda_2H_{1,2_{\beta}}
-\lambda_3 H_{2,1_{\beta}}-\lambda_4 H_{2,2_{\beta}})y,
$$
$$
y^\top(\lambda_0 H_0+\lambda_1 H_{2,1_{\beta}}
+\lambda_2 H_{2,2_{\beta}}
-\lambda_3 H_{1,1_{\beta}}
-\lambda_4 H_{1,2_{\beta}})
y\Big):~y\in\mathbb{R}^{n+1}\Big\},
$$
with
	$$
W_{1_{\beta}}+\mu_1W_2=H_{1,1_{\beta}},\quad W_{1_{\beta}}+\mu_2W_2=H_{2,1_{\beta}},
$$
$$
W_{2_{\beta}}+\mu_1W_2=H_{1,2_{\beta}},\quad W_{2_{\beta}}+\mu_2W_2=H_{2,2_{\beta}},
$$
$$
W_{1_{\alpha}}+\mu_1W_2=H_{1,1_{\alpha}},\quad W_{1_{\alpha}}+\mu_2W_2=H_{2,1_{\alpha}},
$$
$$
W_{2_{\alpha}}+\mu_1W_2=H_{1,2_{\alpha}}\ \quad W_{2_{\alpha}}+\mu_2W_2=H_{2,2_{\alpha}}.
$$
\\
\textbf{Step 3: A second use of a separation result and conclusion.}
By the Dines theorem, $\Omega_2$ is convex. Thus, there exists $(\xi_1,\xi_2)\in\mathbb{R}^2_+\backslash\{0\}$ such that for all $y\in\mathbb{R}^{n+1}$,
$$
\xi_1 y^\top \Big(\lambda_0 H_0+\lambda_1 H_{1,1_{\beta}}+\lambda_2H_{1,2_{\beta}}
-\lambda_3 H_{2,1_{\beta}}-\lambda_4 H_{2,2_{\beta}}\Big)y
$$
$$
+\xi_2y^\top\Big(\lambda_0 H_0+\lambda_1 H_{2,1_{\beta}}
+\lambda_2 H_{2,2_{\beta}}
-\lambda_3 H_{1,1_{\beta}}
-\lambda_4 H_{1,2_{\beta}}\Big)y\geq 0.
$$
In particular, for $y=(x,1)$ with $x\in\mathbb{R}^n$, and by setting
$$\overline{\lambda}_0=\lambda_0(\xi_1+\xi_2),~~\overline{\lambda}_1=(\lambda_1+\lambda_2)(\xi_1+\xi_2),~~\overline{\lambda}_2=(\lambda_3+\lambda_4)(\xi_1+\xi_2),
$$
$$
\mu_{\beta}=\dfrac{\xi_1\mu_1+\xi_2\mu_2}{\xi_1+\xi_2},~
\mu_{\alpha}=\dfrac{\xi_1\mu_2+\xi_2\mu_1}{\xi_1+\xi_2},~\delta_\beta=\dfrac{\lambda_1\delta_1+\lambda_2\delta_2}{\lambda_1+\lambda_2},~{\rm and}~\delta_\alpha=\dfrac{\lambda_3\delta_2+\lambda_4\delta_1}{\lambda_3+\lambda_4},
$$
one gets  $(\overline{\lambda}_0,\overline{\lambda}_1,\overline{\lambda}_2)\in\mathbb{R}^3_+\backslash\{0\}$, $\mu_{\beta},\mu_{\alpha}\in[\mu_1,\mu_2]$, $\delta_{\beta},\delta_{\alpha}\in[\delta_1,\delta_2]$ and
$$
\overline{\lambda}_0\Big(\dfrac{1}{2}x^\top Ax+a^\top x+\gamma\Big) +\overline{\lambda}_1\Big(\dfrac{1}{2}x^\top (B_1+\mu_{\beta} B_2)x+(b_1+\delta_{\beta} b_2)^\top x-\beta\Big)+
$$
$$
\overline{\lambda}_2\left(\alpha-\Big(\dfrac{1}{2}x^\top (B_1+\mu_{\alpha} B_2)x+(b_1+\delta_{\alpha} b_2)^\top x \Big)\right)\geq 0 \quad \forall~ x\in\mathbb{R}^n,$$
where $\mu_{\alpha}+\mu_{\beta}=\mu_1+\mu_2$. This proves that $(b)$ holds.
\qed
\end{proof}

\begin{theorem} $($A robust S-lemma$)$\label{slema}
	Let $A, B_1,B_2\in {\mathcal S}^n$, $a, b_1, b_2\in\mathbb{R}^n$ and $\gamma, \alpha, \beta$, $\mu_1, \mu_2,\delta_1,\delta_2\in\mathbb{R}$, with $\mu_1\leq\mu_2$, $\delta_1\leq\delta_2$ and $\alpha<\beta$. Assume that $\Omega_W$ is convex and that there exists
	$x_0\in\mathbb{R}^n$ satisfying
 \begin{equation}\label{eq3}  \alpha<\dfrac{1}{2}x_0^\top (B_1+\mu B_2)x_0+(b_1+\delta b_2)^\top x_0 <\beta,~\forall~\mu\in [\mu_1,\mu_2], \forall~ \delta\in [\delta_1,\delta_2].
 \end{equation}
Then, the following two assertions are equivalent:
	\begin{itemize}
		\item[$(a)$] $\alpha\leq \dfrac{1}{2}x^\top (B_1+\mu B_2)x+(b_1+\delta b_2)^\top x \leq\beta$, $\forall~ \mu\in [\mu_1,\mu_2]$, $\forall ~\delta\in[\delta_1,\delta_2]$, $\Rightarrow \dfrac{1}{2}x^\top Ax+a^\top x+\gamma\geq 0$.
		\item[$(b)$] $\exists~(\lambda_1,\lambda_2)\in\mathbb{R}^2_+, \exists~\mu_{\alpha},\mu_{\beta}\in [\mu_1,\mu_2 ], \exists~\delta_{\alpha},\delta_{\beta}\in [\delta_1,\delta_2 ]:~ \forall~ x\in \mathbb{R}^n$
		$$
		\dfrac{1}{2}x^\top Ax+a^\top x+\gamma +\lambda_1\Big(	\dfrac{1}{2}x^\top (B_1+\mu_{\beta} B_2)x+(b_1+\delta_{\beta} b_2)^\top x-\beta\Big)+
		$$
		$$ \lambda_2\Big(\alpha-\Big(\dfrac{1}{2}x^\top (B_1+\mu_{\alpha} B_2)x+(b_1+\delta_{\alpha} b_2)^\top x \Big)\Big)\geq 0,
		$$
	\end{itemize}
	where $\mu_{\alpha}+\mu_{\beta}=\mu_1+\mu_2$.
\end{theorem}
\begin{proof} Clearly $(b)\Rightarrow(a)$.\\
Assume now that $(a)$ is satisfied. Then $(a)$ in Theorem \ref{teoralt} does not hold. Thus $(b)$ of the same theorem fulfills, but then $\lambda_0$ is strictly positive because of \eqref{eq3}, which implies the desired result.
\qed
\end{proof}

We are now ready to establish a characterization of optimality
for the problem \eqref{prob:00}.

\begin{theorem}\label{minimonohomoalt}
	$($Characterizing robust optimality$)$   Let $A, B_1,B_2\in {\mathcal S}^n$, $a, b_1, b_2\in\mathbb{R}^n$
	and $\alpha, \beta, \mu_1, \mu_2,$ $\delta_1,\delta_2\in\mathbb{R}$, with $\mu_1\leq\mu_2$, $\delta_1\leq\delta_2$ and $\alpha<\beta$.

	 Let
$\overline{x}$ be feasible for problem \eqref{prob:00} and put $\gamma=-f(\overline{x})$. Assume that $\Omega_W$ is convex, and the
Slater-type condition \eqref{eq3} holds. Then, $\overline{x}$ is optimal if, and only if there exist $(\lambda_1,\lambda_2)\in\mathbb{R}^2_+$, $\mu_{\alpha},\mu_{\beta}\in[\mu_1,\mu_2]$, $\delta_{\alpha}$, $\delta_{\beta}\in[\delta_1,\delta_2]$
	such that the following statements are satisfied:
	\begin{itemize}
		\item[$(a)$] $\Big(A+\lambda_1(B_1+\mu_{\beta}B_2)-\lambda_2(B_1+\mu_{\alpha}B_2)\Big)\overline{x}=-\Big(a+\lambda_1(b_1+\delta_\beta b_2)-\lambda_2(b_1+\delta_{\alpha}b_2)\Big);$
		\item[$(b)$] $\lambda_1\left(\dfrac{1}{2}\overline{x}^\top (B_1+\mu_{\beta}B_2)\overline{x}+(b_1+\delta_{\beta}b_2)^\top \overline{x}-\beta\right)=0;$\\
		$\lambda_2\left(\alpha-\Big(\dfrac{1}{2}\overline{x}^\top (B_1+\mu_{\alpha}B_2)\overline{x}+(b_1+\delta_{\alpha}b_2)^\top \overline{x}\Big)\right)=0;$
		\item[$(c)$] $A+\lambda_1(B_1+\mu_{\beta}B_2)-\lambda_2(B_1+\mu_{\alpha}B_2)\succeq 0$.
	\end{itemize}
\end{theorem}

\begin{proof} The necessary condition follows from Theorem \ref{slema} where $\gamma$ is substituted by $-f(\overline{x})$. Indeed, if $\overline{x}$ is optimal then $(a)$ in Theorem \ref{slema} holds, which means that $(b)$ of the same theorem is also satisfied. This finally implies the desired statements.\\
The sufficiency part is already standard since $(c)$ implies the convexity of the function
	$$
	h(x)\doteq\dfrac{1}{2}x^\top Ax+a^\top x-f(\overline{x}) +{\lambda}_1\Big(\dfrac{1}{2}x^\top (B_1+\mu_{\beta} B_2)x+(b_1+\delta_{\beta} b_2)^\top x-\beta\Big)+
	$$
	$$
	{\lambda}_2\left(\alpha-\Big(\dfrac{1}{2}x^\top (B_1+\mu_{\alpha} B_2)x+(b_1+\delta_{\alpha} b_2)^\top x \Big)\right),$$
	and $(a)$ and $(b)$ allow us to prove that $\overline{x}$ is in fact a solution to problem \eqref{prob:00}.
	\qed
\end{proof}

\section{Revisiting the case $\alpha=-\infty$}
\label{secc:uni}

We consider the problem:

\begin{equation}\label{uni0}
\begin{split}
\min &~ \dfrac{1}{2} x^\top Ax+a^\top x  \\
\text{s.t.}~ & \dfrac{1}{2}x^\top (B_1+\mu B_2)x+(b_1+\delta b_2)^\top\leq\beta,~ \forall~ \mu\in [\mu_1,\mu_2],\forall~\delta\in[\delta_1,\delta_2],
\end{split}
\end{equation}
where $A$, $B_1$, $B_2$ are real symmetric matrices, $\mu_1,\mu_2,\beta\in\mathbb{R}$ satisfying
$\mu_1< \mu_2$, $\delta_1<\delta_2$.
This problem was also discussed in \cite{jeya2013}, and actually this paper motivated our study.

By looking at carefully the proof of Theorem \ref{teoralt}, we immediately realize that in case there is no lower bound in the inequality constraint, all the terms where $\alpha$ appears, actually dissapear: they are superfluous. Hence, the set
$\Omega_W$ reduces to
$$
\Omega^\beta_W\doteq\Big\{\Big(\dfrac{1}{2}y^\top H_0y
,\max_{\mu\in[\mu_1,\mu_2]}\dfrac{1}{2}y^\top (W_{1_{\beta}}+\mu W_2)y),\max_{\mu\in[\mu_1,\mu_2]}\dfrac{1}{2}y^\top (W_{2_{\beta}}+\mu W_2)y\Big):
$$
\begin{equation
y\in\mathbb{R}^{n+1}\Big\}+\rm{int} \ \mathbb{R}^3_+.
\end{equation}
Thus, by Corollary 1 in \cite{fbp2022}, $\Omega_W^\beta$ is convex if the set
$$\Omega^\beta_\mu\doteq
\Big\{\Big(y^\top H_0y,y^\top(W_{1_{\beta}}+\mu_1 W_2)y,
y^\top(W_{2_{\beta}}+\mu_1 W_2)y, y^\top (W_{1_{\beta}}+\mu_2 W_2)y,
$$
$$
y^\top (W_{2_{\beta}}+\mu_2 W_2)y\Big):y\in\mathbb{R}^{n+1} \Big\}+\rm{int} \ \mathbb{R}^5_+.
$$

\begin{theorem}\label{minimocorregido} Let the data be as described above. Let $\overline{x}$ be feasible for  problem \eqref{uni0} and put $\gamma=-f(\overline{x})$. Assume that $\Omega_W^\beta$ is convex, and the
Slater-type condition:  there exists $x_0\in\mathbb{R}^n$ such that
\begin{equation
\dfrac{1}{2}x_0^\top (B_1+\mu B_2)x_0+(b_1+\delta b_2)^\top x_0<\beta, \ \forall~ \mu\in[\mu_1,\mu_2], \ \forall ~\delta\in[\delta_1,\delta_2]
\end{equation}
is satisfied. Then, $\overline{x}$ is optimal if, and only if there exist $\lambda\geq 0$, $\mu\in[\mu_1,\mu_2]$, $\delta\in[\delta_1,\delta_2]$
such that the following statements are satisfied:
\begin{itemize}
	\item[$(a)$] $\Big(A+\lambda(B_1+\mu B_2)\Big)\overline{x}=-\Big(a+\lambda(b_1+\delta b_2)\Big);$
	\item[$(b)$] $\lambda\left(\dfrac{1}{2}\overline{x}^\top (B_1+\mu B_2)\overline{x}+(b_1+\delta b_2)^\top \overline{x}-\beta\right)=0;$
	\item[$(c)$] $A+\lambda (B_1+\mu B_2)\succeq 0$.
\end{itemize}
\end{theorem}

We recall that
$$
H_0=\begin{pmatrix}A & a\\
a^\top  &  2\gamma \end{pmatrix},
$$
with $\gamma=-f(\overline{x})$ as in the previous theorem.
Set
$$
H_1\doteq H_{1,1_{\beta}},~H_2\doteq H_{2,2_{\beta}}, H_3\doteq H_{1,2_{\beta}},~H_4\doteq H_{2,1_{\beta}}.
$$
The authors in \cite[Theorem 5.1]{jeya2013} proved the same result expressed in
Theorem \ref{minimocorregido}  under the convexity of the set
\begin{equation}\label{setw00}
\Big\{\Big(y^\top H_0y,y^\top H_1y, y^\top H_2y\Big):~y\in\mathbb{R}^{n+1}\Big\}+
{\rm int}~\mathbb{R}_+^3;
\end{equation}
whereas ours requires the convexity of
\begin{equation*}
\Omega_\mu^\beta=
\Big\{\Big(y^\top H_0y,y^\top H_1y, y^\top H_2y, y^\top H_3y, y^\top H_4y
\Big):~y\in\mathbb{R}^{n+1}\Big\}+{\rm int}~\mathbb{R}_+^5.
\end{equation*}
We believe that there is a gap in the proof of Theorem 5.1 in \cite{jeya2013}. More precisely, the authors assert in page 221 of the same paper that  the nonexistence of solution to the system (notice that our $\beta$ is $-\beta$ in \cite{jeya2013})
\begin{equation}\label{sis:hom10}
\dfrac{1}{2}x^\top Ax+ta^\top x + t^2\gamma<0,~ \max\left\{\dfrac{1}{2}x^\top Bx+tb^\top x-t^2\beta:(B,b)\in\mathcal{B}_b\right\}<0,
\end{equation}
implies that
\begin{equation}\label{ecua00}
\dfrac{1}{2} y^\top  H_0 y<0\quad \text{and}\quad \forall~ \mu\in[\mu_1,\mu_2], \quad \dfrac{1}{2}y^\top (W_1+\mu W_2)y<0,
\end{equation}
has no solution. This is not necessarily true as Example \ref{ex100} below shows. We recall that
$$
W_1=\begin{pmatrix}B_1 & b_1+\dfrac{\delta_1\mu_2-\delta_2\mu_1}{\mu_2-\mu_1}b_2\\
(b_1+\dfrac{\delta_1\mu_2-\delta_2\mu_1}{\mu_2-\mu_1}b_2)^\top  &  -2\beta
\end{pmatrix}
\quad \text{and}\quad
$$
$$
W_{2}=\begin{pmatrix}B_2 & \dfrac{\delta_2-\delta_1}{\mu_2-\mu_1}b_2\\
\dfrac{\delta_2-\delta_1}{\mu_2-\mu_1}b_2^\top  &  0 \end{pmatrix}.
$$
\vskip0.3truecm
\begin{example} \label{ex100} Let  $n \geq 5$.  Taking $A=I_n$, $B_1=B_2=2I_n$,  $a=b_1=b_2 =s:=(1,1,\ldots,1) \in \mathbb{R}^n$,
$\mu_1=-1$,  $\mu_2=1$,  $\delta_1=-1, $  $\delta_2=1$ and $\beta =1$,  Problem \eqref{uni0} takes the form
\begin{equation}\label{example_j}
	\begin{split}
	\min ~~ & \dfrac{1}{2} \|x\|^2 + \sum_{i=1}^n x_i = \frac{1}{2} \|x+s\|^2 -\frac{n}{2}  \\
	\text{s.t.}~~ & g(x, \mu, \delta):=(1+\mu)\|x\|^2 + (1+\delta) \sum_{i=1}^n x_i -1 \leq 0~~\forall ~\mu,\delta\in[-1,1].
	\end{split}
	\end{equation}
	Let us define the functions
$$
f(x):=\dfrac{1}{2} \|x\|^2 + \sum_{i=1}^n x_i = \frac{1}{2} \|x+s\|^2 -\frac{n}{2};
$$
$$
g(x, \mu, \delta):=(1+\mu)\|x\|^2 + (1+\delta) \sum_{i=1}^n x_i -1.
$$
Then, problem \eqref{example_j} is equivalent to
	\begin{align}
	\min ~~& f(x)\notag\\
	\hbox{ s.t.}~~ &\|x\|^2 -\frac{1}{2} \leq 0; \label{eq1_example_j}\\
	& \sum_{i=1}^n x_i -\frac{1}{2} \leq 0;  \label{eq2_example_j}\\
	&\|x\|^2 +  \sum_{i=1}^n x_i -\frac{1}{2} \leq 0. \label{eq3_example_j}
	\end{align}
    Let $C$ be the set of constraints, that is, those $x$ satisfying
    \eqref{eq1_example_j}-\eqref{eq3_example_j}. It is clear that $C$ is convex and compact. Thus, the unique solution to problem \eqref{example_j} is the projection
    of $-s$ on $C$, which is  $\overline{x} = -\dfrac{1}{\sqrt{2n}}s$. Hence,
   $\overline{x}$ is a robust solution to problem \eqref{example_j}.  In
	 this case, $\gamma = -f(\overline{x}) = \sqrt{\dfrac{n}{2}} -\dfrac{1}{4}$. By identifying the matrices involved in Theorem \ref{minimocorregido}, we get
	\[
	H_0=\begin{pmatrix}
	I_n & s \\
	s^\top & 2\gamma
	\end{pmatrix};
	H_1= \begin{pmatrix}
	\overline{0} & 0 \\
	0^{\top} & -2
	\end{pmatrix};
	H_2= \begin{pmatrix}
	4I_n & 2s \\
	2s^{\top} & -2
	\end{pmatrix};
	\]
	\[
	H_3=\begin{pmatrix} \overline{0} &2s\\
	2s^\top  &  2\beta \end{pmatrix};\quad H_4=\begin{pmatrix}4I_n  & 0\\
	0^\top  &  2\beta \end{pmatrix}.
	\]
	Since
	 \begin{equation}
	-2H_0+ (-2-2\gamma) H_1+ H_2= \begin{pmatrix}
	2I_n & 0 \\
	0^{\top} & 2
	\end{pmatrix} \succ 0,
	\end{equation}
  applying \cite[Theorem 2.1 ]{Polyak} we have that
\begin{equation}\label{setw:0}
\Big\{\Big(y^\top H_0y,y^\top H_1y, y^\top H_2y\Big):~y\in\mathbb{R}^{n+1}\Big\}
\end{equation}
is convex.  This means, according to \cite {jeya2013}, that  \eqref{example_j} is regular with respect to $\overline{x}$.
Also,  by taking $x_0=0$,   we have  $g(x_0, \mu, \delta)  <0$, for all  $\mu, \delta \in [-1, 1]$.  So,  we have all the conditions of   \cite[Theorem 5.1]{jeya2013} are satisfied. From the first part of the proof of \cite[Theorem 5.1]{jeya2013} we have that the following homogeneous system in $\mathbb{R}^{n+1}$, \eqref{sis:hom10} has no solution.\\
Coming back to our example, we obtain
\[W_1=\begin{pmatrix}
	2I_n & s \\
	s^{\top} & -2
	\end{pmatrix} \hbox{ and }
	W_2=\begin{pmatrix}
	2I_n & s \\
	s^{\top} & 0
	\end{pmatrix}.
	\]
Then, the following homogeneous system in $\mathbb{R}^{n+1}$ (see \eqref{ecua00})
\begin{equation*}
\frac{1}{2}\begin{pmatrix} x\\t\end{pmatrix}^{\top} H_0\begin{pmatrix} x\\t\end{pmatrix}  <0 \hbox{ and }
\forall ~\mu \in [-1, 1],~\frac{1}{2}\begin{pmatrix} x\\t\end{pmatrix}^{\top}(W_1+\mu W_2)\begin{pmatrix} x\\t\end{pmatrix}<0,
\end{equation*}
becomes
\begin{equation}\label{ecua01}
\frac{1}{2}\begin{pmatrix} x\\t\end{pmatrix}^{\top}\begin{pmatrix}
	I_n & s \\
	s^{\top} & 2\gamma
	\end{pmatrix}\begin{pmatrix} x\\t\end{pmatrix}  <0 \hbox{ and }
\end{equation}
\begin{equation}\label{ecua02}
\forall~ \mu \in [-1, 1],~\frac{1}{2}\begin{pmatrix} x\\t\end{pmatrix}^{\top}\begin{pmatrix}
	2(1+\mu)I_n & \; & (1+\mu)s \\
	(1+\mu)s^{\top} & \;& -2
	\end{pmatrix}\begin{pmatrix} x\\t\end{pmatrix}<0.
\end{equation}
We will see that such a system admits a solution, contradicting the assertion made in page 221 of \cite{jeya2013} about the nonexistence of solution to the same system. Indeed, for $(-s,1) \in \mathbb{R}^{n+1}$, \eqref{ecua01} reduces to
\[
\frac{1}{2}\begin{pmatrix} -s\\1\end{pmatrix}^{\top}\begin{pmatrix}
	I_n & s \\
	s^{\top} & 2\gamma
	\end{pmatrix} \begin{pmatrix} -s\\1\end{pmatrix} = \frac{1}{2}(-n+2\gamma )= \frac{1}{2}\big(-n+  \sqrt{\frac{n}{2}} -\frac{1}{4}\big) < 0;
\]
whereas \eqref{ecua02} becomes: for all $\mu \in [-1, 1]$,
	 \[
	\frac{1}{2}\begin{pmatrix} -s\\1\end{pmatrix}^{\top}\begin{pmatrix}
	2(1+\mu)I_n & \;& (1+\mu)s \\\
	(1+\mu)s^{\top} & \; & -2
	\end{pmatrix}\begin{pmatrix} -s\\1\end{pmatrix}= 	\frac{1}{2}(n(1+\mu)-(1+\mu)n -2) =-1 <0.
	\]
This proves our claim.

Observe that taking  in Corollary \ref{cor_R_5}: $m=4$, $A_0=I_n$, $\rho_1=\rho_3=4$,  $\rho_2=\rho_4=0$,
$a_0= s, a_1=a_3=0$,  $a_2=a_4=2s$,  $c_0= \gamma$, $c_1=c_2= -1$, $c_3=c_4= \beta$ we obtain that
$$
\Big\{\Big(y^\top H_0y,y^\top H_1y, y^\top H_2y, y^\top H_3y, y^\top H_4y\Big):~y\in\mathbb{R}^{n+1}\Big\} \subset \mathbb{R}^5
$$
is a convex set, which is required in our Theorem \ref{minimocorregido}, providing the characterization of robust optimality for our example.
\qed\end{example}

Clearly, the previous example only shows there is a gap in the proof of Theorem 5.1 in \cite{jeya2013}. We were unable to construct a real counterexample to that result under the convexity either of the set given in \eqref{setw00} or in \eqref{setw:0}.


\section{Introduction}

Deep Learning (DL) is a subfield of Machine Learning (ML) that focuses on developing Deep Neural Network (DNN) models \cite{shahroudnejad2021survey}. DNNs are complex models that are capable of achieving high performance on a variety of tasks. Many deep neural network models are uninterpretable black-boxes, which often leads to users having less trust in them \cite{miller2019but}. This lack of interpretability and trust can have negative consequences – people might use an AI that makes mistakes, or not use an AI that could improve the chances of a desired result. To improve the interpretability and trustworthiness of black-box models, Explainable Artificial Intelligence (XAI) research focuses on developing methods to explain the behaviour of these models in terms that are comprehensible to humans \cite{molnar2020interpretable}.

Explanations provided by XAI can impact human behaviour and cognition \cite{donadello2020explaining, dragoni2020explainable}. This paper aims to apply frameworks from research on behaviour change within the broader field of behavioural science \cite{ruggeri_behavioral_2018}. Two notable frameworks for thinking about and categorising behaviour change techniques are ‘nudges’ and ‘boosts’ \cite{grune2016nudge}. A nudge is any aspect of choice architecture – the context in which people make decisions - aimed at influencing the behaviours of individuals without limiting or forcing options \cite{sunstein2014nudging}. Boosts are interventions that foster people’s competencies through changes in knowledge and skills, so that they can make their own choices more effectively \cite{hertwig2017consider}. By doing so this paper seeks to put forward a framework and a paradigm for evaluating whether within the context of human-AI interaction an XAI is nudging performance or boosting capability. It will discuss the implications of each in relation to human-machine collaboration as well as the ethics of influence.

\section{Different explanations, different outcomes}

A great deal of XAI research has gone into developing methods to improve the explainability of DNN models \cite{molnar2020interpretable}. As it stands it is not fully clear how different methods may influence human behaviour.

Feature importance methods (such as saliency methods) generate scores that reveal how important a feature (like a word vector or pixel) is to the AI's decision-making process \cite{bhatt2020explainable}. Explanations generated by these methods can be either global or local in scope \cite{lundberg2020local}. Global explanations, like those provided by SHAP (SHapley Additive exPlanations) models, give a quantitative indication of the importance of each input variable on the model's output \cite{lubo2020machine}. Local explanation methods, such as LIME (Local Interpretable Model Agnostic Explanations), generate a numeric score showing the importance of an input variable in relation to the outcome variable \cite{lee2019developing}.

Counterfactual explanations show what the model's output would have been if one or more of the inputs had been different \cite{verma2020counterfactual}. This can be helpful in understanding why the model arrived at a particular output. Finally, concept-based explanations attempt to explain a model's output by referencing pre-defined or automatically generated sets of concepts that are comprehensible to humans \cite{kazhdan2020now}.

Even though there has been progress in the development of XAI models, it is still not completely clear which model should be used for human-machine collaboration, and for what purpose. There are still many open question. For example, do local explanations nudge performance in the short term but not necessarily provide enough information to educate a user and boost their capability in the long term? Do global explanations provide enough information to educate? Can counterfactual explanations teach people how their AI works or allow them to identify errors?

A systematic literature review of 241 papers looked at how the validity and usefulness of explanations have been evaluated by the authors of XAI methods \cite{anjomshoae2019explainable}. Most studies only conducted user studies in simple scenarios or completely lacked evaluations. The results show that 32\% of the research papers did not include any type of evaluation. Furthermore, 59\% of the research papers conducted a user study to evaluate the usefulness of the explanation (with a small minority also evaluating user trust towards the AI system). Finally, 9\% of the papers used an algorithmic evaluation, which did not involve any empirical user research. These initial findings suggest that different explanations will lead to variations in performance on a task \cite{lage2019evaluation, narayanan2018humans}, and will not always necessarily improve performance and understanding \cite{kindermans2019reliability}. If viewed as a behaviour change intervention, when does an XAI serve as a nudge, changing behaviour, and when does it serve as a boost, improving capability?

\section{Local and concept-based explanations as nudges}

 The common aim of nudges is to predictively change targeted behaviours. A nudge is any aspect of choice architecture aimed at influencing people’s behaviour, without limiting or forcing options, or significantly changing their economic incentives \cite{thaler2008nudge}. All environments influence behaviour to some extent, even when people are not aware of it. Intentionally changing choice architecture is nudging. Nudges take many shapes \cite{sunstein2014nudging}. Default rules such as automatic enrollment in programs, automate decision-making for people by setting a default. Simplification nudges reduce the amount of information presented to people to avoid information overload. The use of descriptive social norms - telling people what most other people are doing - influences behaviour. As a policy tool, nudging has been used in over 80 countries worldwide, and by major supranational institutions such as the World Bank and UN \cite{oecd2017behavioural}.

Nudges have been heavily influenced by Daniel Kahneman's dual-process account of reasoning \cite{kahneman2003perspective}. He proposed that people have "two systems" in their mind - System 1 and System 2. System 1 thinking is heuristic. It reacts intuitively and effortlessly, without analysing all available information. System 2 is an analytical and effortful, rationalising process. System 1 thinking is fast, and thus accounts for most behaviour. System 2 can re-evaluate System 1 thinking, thus using System 2 thinking leads to fewer erroneous decision. However, this is difficult, as it requires more cognitive effort. Importantly, some factors and contexts are more likely to trigger System 1 or System 2 thinking than others. Per \citet{sunstein2016ethics}, nudges work by targeting either System 1 thinking, thus influencing behaviour without the awareness of the decision maker, or System 2 thinking, thus promoting deliberative thinking.

A famous example of nudging is using disclosure nudges. Disclosure nudges disclose decision-relevant information \cite{sunstein2014nudging}. They are educative, because they provide a learning experience, and target System 2, by promoting deliberative thinking. Disclosure nudges are rooted in three insights. First, as uncertainty promotes erroneous decision making \cite{kochenderfer2015decision}, disclosure nudges seek to reduce uncertainty with decision-relevant information. Second, when too-much decision-irrelevant information is present, people find decision-making more challenging \cite{rogers2013infobesity}; disclosing only decision-relevant information can, therefore, reduce error. Finally, disclosure nudges create an emphasis frame, making relevant information more salient \cite{chong2007framing}.

Viewed through this framework, local feature importance explanations and concept-based explanations can be viewed as a type of disclosure nudge. They will provide a small amount of decision-relevant information at the time a person needs to decide whether or not to trust an AIs advice or prediction. It thus could be the case that findings from the extensive body of research on disclosure nudges may generalise to local and concept-based explanations. Such comparisons can guide future practice and research directions.

\section{Global and counterfactual explanations as boosts}

Boosts are interventions that aim to promote people’s competencies, so they can make better decisions \cite{grune2016nudge}. Proponents of boosts aim to foster skills and decision heuristics that can persist over time, throughout different decision contexts \cite{hertwig2017nudging}. An example of a boost is teaching people better decision-making skills with the use of decision trees \cite{hertwig2009fast}, or by teaching them how to calculate the expected value of a prospect \cite{franklin2019optimising}. Unlike nudges, a boosts effectiveness requires people’s awareness of the boost and motivation to improve their competencies. Furthermore, boosts differ from educational nudges, in that they do not guide people towards a decision, but rather they assume that better-skilled people will make advantageous decisions. In the context of a decision, boosts provide people with the right capability for the task at hand.

Boosts can either promote abilities that are specific to a single domain (e.g., decision-making under risk) or generalise across many domains (e.g., statistical literacy; \citet{hertwig2017nudging}). Furthermore, boosts can change the choice architecture to foster competencies, directly teach people skills, or do both. Boosts have been used to improve people’s decision-making capabilities by promoting people’s understanding of statistical information. This competence has been previously achieved through: using graphical representations of statistical information \cite{lusardi2017visual}, training math skills \cite{berkowitz2015math}, representing information that avoids framing effects \cite{gigerenzer2007helping}, and giving people simulation based representations of statistical information \cite{hogarth2015providing}.

Through the lens of behaviour change theories, global feature importance explanations and counterfactual explanations can be viewed as a type of boost. Counterfactual explanations can be used to train people to understand when their AI tool is more likely to err. This may generalise as a more broader knowledge about how AI tools "work" more broadly. Global explanations may increase people's knowledge about what factors are relevant for a given context. This can improve people's own predictions and competencies within a certain domain. As was the case with nudges, findings from the research literature on boosts may generalise to global and counterfactual explanations. These comparisons could guide both practice and future research on improving human-AI collaboration.

\section{Evaluating XAI Influence}

In light of the previous discussion on certain XAIs as performance boosting and others as capability boosting, this section will propose  methods for testing whether a explanation is nudging performance or boosting capability. It will outline human-AI interaction paradigms for evaluating an XAI's influence.

Measures of performance will be highly context-specific \cite{hoffman2018metrics}. It is also possible to identify whether explanations improve performance by helping people do better, or by reducing the number of mistakes people make \cite{franklin2022human}. Ideally, performance on the same task without an AI will be compared to performance on a task with an AI, and with an XAI. Multiple different XAIs can be used in this comparison. It could either between-subject - comparing the performance of people across different groups - or within-subject - tracking the performance of an individual across time.

To give an example, one could compare the effects of local and global explanations for improving a doctor's diagnostic ability. Specifically, the study would evaluate the benefits of an AI diagnostic tool on performance, and see how this improves when people receive explanations for how that diagnostic tool makes decisions. In this example, doctors will be tested on their ability to detect pneumonia from x-ray images. In a between-subject paradigm, doctors could be randomly placed in four separate groups: 1) doctors with no AI tools, 2) doctors with AI tools, 3) doctors with AI tools and a local explanation, 4) doctors with AI tools and a global explanation. The difference between the first and second group would be a pure measure of how the AI tool improves pneumonia detection. The difference between the second and third, or the second and fourth group would allow one to look at the benefits of the explanation on performance. Finally, the difference between the third and fourth group would serve as a measure of what XAI is more effective in improving performance.

A similar paradigm can be used to measure capability. An AI without an explanation can allow people to understand the AI's decision-making process by detecting patterns \cite{peltola2019machine}. XAI can provide a more direct learning experience for users. This learning could be \textit{procedural}, for example a change in a certain ability, or \textit{semantic}, for example an increase in factual knowledge \cite{berliner2013handbook}. Procedural knowledge can be measured as people's performance in the absence of an AI tool after using an AI tool. The pause between initially using the AI tool and then measuring performance without it can be more or less longitudinal. Longer time frames would make for a more valid measure as they would reduce the confounding impact memory may have on performance. Back to the previous example, each doctor can be brought back for testing 2 weeks after the initial study. In this case they would be performing the same activity - pneumonia detection from x-ray images - but this time with no AI tools or explanations. A lack of differences between their previous and current performance would suggest that their interaction with the AI or XAI models led to an increase in capability. Measuring semantic knowledge with a similar paradigm  would involve direct, context-specific questions to test a doctor's knowledge.

Finally, measuring changes in people's mental models can provide useful information for whether an XAI has led to new knowledge and capability. Mental models are representations of a person's understanding of some system or object \cite{lagnado2021explaining}. An XAI could change the person's mental models about the task, domain or the AI tool. Given that people are able to infer causal structures from explanations \cite{kirfel2021inference}, explanations can both establish the presence and change the direction of perceived cause and effect relationships. Mental models can be measured using a nearest neighbor task, where participants select the explanation or diagram that best fits their beliefs, or with concept mapping, where users create a diagram which outlines their knowledge \cite{hoffman2018metrics}.

\section{Issues with XAI influence}

Not understanding how an XAI's explanation influences behaviour and cognition raises three broader issues. First, there may be individual differences in the way explanations influence specific users on certain tasks \cite{tomsett2018interpretable}. A review of 137 articles on XAI applied to different domains shows that although there is preliminary evidence for visual explanations being more acceptable to everyday users, most studies have been directed at expert users \cite{islam2022systematic}. More research evaluating everyday users' reactions to XAI is needed \cite{franklin2022human}. The person receiving the information matters.

Second, it is not evident what amount of and type of information a user should receive. Detailed explanations might conceal more information than they reveal. This phenomenon is known as \textit{information overload} in psychology – when a person receives too much information, or more specifically when the amount of input to a system exceeds its processing capacity \cite{sutcliffe2009information}. This is also true even when all of the information is directly relevant to the task at hand. Recent research has found that information overload effects can occur in response to explanations provided by XAI methods \cite{ferguson2022explanations}. In general, it was found that the more detailed the explanation, the less useful and trustworthy it was considered to be. Specifically, non-numerical explanations in plain English that listed the variables involved in a model were rated as more understandable and trust-worthy than SHAP models.

Other research has found the opposite to be true for highly technical individuals. A 2021 study found that people with backgrounds in AI have different requirements for explanations \cite{ehsan2021explainable}. Specifically, individuals with a background in AI get more use out of numerical explanations compared to less technical individuals. This means that it is not possible to make generalizations, such as that more information will always lead to better decisions.

Finally, persuasive explanations will influence behaviour and preferences in ways that are currently difficult to predict. Current research suggests that an XAI tool, as opposed to a regular AI tool, is more likely to produce a change in behaviour \cite{donadello2020explaining, dragoni2020explainable}. Researchers have also developed XAI methods which can generate misleading explanations that can both increase trust in the AI and mislead domain experts \cite{lakkaraju2020fool}. It is therefore essential to understand the direction of the behaviour change, whether it is desirable, and how this will differ from one XAI method to another. This is especially relevant given the fact that preference and behaviour have a bidirectional causal relationship \cite{ariely2008actions}, which means that AI systems can influence preferences by changing behaviours \cite{ashton2022problem, franklin2022recognising}. By pushing behaviour in a certain direction, over time, people's preferences for what to do and how to do it can change.

\section{Conclusions}

This paper proposed a framework for understanding XAI influence, and a method for measuring it. It proposed that local feature importance explanations and concept-based explanations are adjacent to disclosure nudges, and global feature importance explanations and counterfactual explanations are adjacent to boosts. It is thus possible that the extensive research conducted on nudges and boosts can provide insights for how XAI will influence behaviour. Translational research using existing insights from AI, Behavioural Science and Human-Computer Interaction  would be beneficial. Such research could expand our understanding of influence dynamics in human-AI interaction.

There are at least three reasons for why understanding XAI influence is important for better human-AI collaboration and teaming. First, the adequate use of explanations allows people to identify AI errors; namely due to spurious correlations \cite{sagawa2020investigation} or runaway feedback loops \cite{ensign2018runaway}. This is especially relevant in light of recent work on making models robust to spurious correlations by leveraging humans’ common sense knowledge of causality \cite{srivastava2020robustness}. Identifying XAIs for promoting this error-detecting ability would be beneficial.

Second, it is important to understand which XAIs result in unwanted or negative (i.e., capability and performance decreasing) influence or \textit{sludge} \cite{thaler2018nudge}. Sludge can encourage self-defeating behaviours or discourage a person's best interest. A common form of sludge is friction by making something slower or unnecessarily complicated \cite{sunstein2018sludge}. Research showing that certain forms of XAI result in information overload is an example of this \cite{ferguson2022explanations}. \citet{sunstein2020sludge} argues for sludge audits whereby one identifies the sources of sludge and eliminates them. Similar practices should be used in XAI.

Finally, understanding the influence of XAI can allow for a more optimal distribution and selection of XAI methods. Some people performing certain tasks will want performance boosts. Others will want to learn. Understanding XAI influence would allow us to match these preferences.


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

    $W^{\pm}$ vector bosons are heavy particles with masses of $ m_{W^\pm}
\approx 80$~GeV/$c^{2}$~\cite{ParticleDataGroup:2020ssz,CDF:2022hxs}. They are
mainly produced in the hard partonic scattering processes with large momentum
transfer at the early stage of the (ultra-)relativistic nuclear
collisions. In comparison with evolution time of the heavy-ion collision
system, 10 to 100~fm/$c$ for instance, the $W^{\pm}$ decay time
($\sim 0.0922~{\rm fm/}c$, estimated with full decay width, $t=h/\Gamma$
\cite{ParticleDataGroup:2020ssz}) is very short. The $W^{\pm}$ leptonic decays
\begin{equation*}
    W^{\pm}\rightarrow l^{\pm}\nu_{l},~(l{\rm :}~e,\mu,\tau)
\end{equation*}
are nearly instantaneous. As the produced leptons weakly interact with the
partonic and hadronic matter created in nuclear collisions, the $W^{\pm}$ is a
powerful probe for investigating the properties of the initial state in the
nuclear collision system. The main production processes of $W^\pm$ are
\begin{equation*}
    u\bar{d} \rightarrow W^{+}, \hspace{0.4cm} d\bar{u} \rightarrow W^{-}
\end{equation*}
in the leading-order (LO) approximation~\cite{ParticleDataGroup:2020ssz,Martin:1999ww}.
With the involvement of the valence $u$- and $d$-quarks inside the nuclear
system, the isospin difference between proton and neutron is important for the
production of $W^{+}$ and $W^{-}$. It is therefore expected that the $W^\pm$
boson charge asymmetry and the corresponding decayed lepton charge asymmetry
are related to the relative abundance of protons and neutrons (to the relative
abundance of valence $u$- and $d$-quarks) of the collision system. The lepton
charge asymmetry $A_{l}$ could be defined as the difference between the
multiplicity of positive and negative leptons divided by their sum:
\begin{equation}
    A_{l} = \frac{N_{l^{+} \leftarrow W^{+}} - N_{l^{-} \leftarrow W^{-}}}
	    {N_{l^{+} \leftarrow W^{+}} + N_{l^{-} \leftarrow W^{-}}}.
\label{eq:A_l}
\end{equation}
Similarly, the collision system valence $u$ ($d$) quark number asymmetry
$A_v$ is defined as
\begin{equation}
A_v = \frac{N_u - N_d}{N_u + N_d},
\end{equation}
where $N_u$ ($N_d$) refers to the number of valence quark $u$ ($d$) in
the nuclear collision system.

    The lepton charge asymmetry have been measured in $p\bar p$ collisions at
$\sqrt{s}=1.8$ and $1.96$~TeV by CDF~\cite{CDF:1991igh,CDF:1994cas,CDF:1998uzn,CDF:2005cgc,CDF:2009cjw}
and D0~\cite{D0:2007pcy,D0:2008cgv,D0:2013xqc} collaborations on Fermilab
Tevatron. They give early complementary constrains on Parton Distribution
Function (PDF) relative to the deep inelastic scattering experiments.
More precise measurements in $pp$ collisions at $\sqrt{s}$=7 and 8~TeV are
reported by ATLAS~\cite{ATLAS:2011pph,ATLAS:2011qdp,ATLAS:2016nqi,ATLAS:2019fgb},
CMS~\cite{CMS:2011aa,CMS:2012ivw,CMS:2012fgk,CMS:2013pzl,CMS:2016qqr} and
LHCb~\cite{LHCb:2012lki,LHCb:2014liz,LHCb:2015mad,LHCb:2016zpq} collaborations
on the Large Hadron Collider (LHC) at CERN. ATLAS also give the measurements
in $pp$ collisions at $\sqrt{s}=2.76$~\cite{ATLAS:2019fyu} and $5.02$~TeV
~\cite{ATLAS:2018pyl}. Beyond the elementary $pp$ collisions, the lepton charge
asymmetry has been measured in the p+Pb collisions at
$\sqrt{s_{\rm NN}}=5.02$~TeV by ALICE~\cite{ALICE:2016rzo} and at
$\sqrt{s_{\rm NN}}=8.16$~TeV by ALICE~\cite{ALICE:2022cxs} and CMS~\cite{CMS:2019leu}.
The results of lepton charge asymmetry in Pb+Pb collisions at
$\sqrt{s_{\rm NN}}=2.76$~TeV are first given by ATLAS~\cite{ATLAS:2014sic} and
CMS~\cite{CMS:2012fgk}. Recently, ATLAS and ALICE published their results
of $W^\pm$ production cross-section and lepton charge asymmetry in Pb+Pb
collisions at $\sqrt{s_{\rm NN}}=5.02$~TeV~\cite{ATLAS:2019ibd,ALICE:2022cxs}.
All those measurements are declared to be well reproduced by the LO and/or
next-to-leading-order (NLO) perturbative Quantum ChromoDynamics (pQCD)
calculations~\cite{Dulat:2015mca,Eskola:2009uj,Kusina:2016fxy,Eskola:2016oht}
using the CT14 PDF set~\cite{Dulat:2015mca} with or without the nuclear
modified PDF (nPDF).

    In this paper, the parton and hadron cascade model PACIAE \cite{Sa:2011ye} is
employed to reproduce the ATLAS and ALICE lepton charge asymmetry data
in Pb+Pb at $\sqrt{s_{\rm NN}}=5.02$~TeV~\cite{ATLAS:2019ibd,ALICE:2022cxs}. We
then calculate the lepton charge asymmetry as a function of the collision
system valence quark number asymmetry in both the elementary nucleon-nucleon
(NN) and the nucleus-nucleus (AA) collisions at $\sqrt{s_{\rm NN}}=5.02$~TeV. An
interesting linear scaling behavior in the lepton charge asymmetry as the
function of collision system valence quark number asymmetry is observed
among the different NN and AA collision systems.

    The paper is organized as follows. In Sec.~\ref{sec:mod}, we describe the
PACIAE model and the setups for $W^\pm$ production in \textsc{Pythia} 6.4,
with a precision equivalent to the NLO pQCD. The simulated results and
discussions are presented in Sec.~\ref{sec:res}. At last, a summary is
given in Sec.~\ref{sec:summ}.

\section{\label{sec:mod} Model}

    A parton and hadron cascade model PACIAE~\cite{Sa:2011ye} is employed to
simulate the elementary nucleon-nucleon (NN: $pp$, $pn$, $np$ and $nn$)
and the nucleus-nucleus (AA: Cu+Cu, Xe+Xe, Au+Au, Pb+Pb and U+U)
collisions at $\sqrt{s_{\rm NN}}=5.02$~TeV.

    The PACIAE model is based on \textsc{Pythia} event generator (version $6.4.28$)
\cite{Sjostrand:2006za}. For NN collisions, with respect to \textsc{Pythia},
the partonic and hadronic rescatterings are introduced before and after
the hadronization, respectively. The final hadronic state is developed
from the initial partonic hard scattering and parton showers, followed by parton
rescattering, string fragmentation, and hadron rescattering stages. Thus,
the PACIAE model provides a multi-stage transport description on the evolution
of the NN collision system.

    For AA collisions, the initial positions of nucleons in the colliding nucleus
are sampled according to the Woods-Saxon distribution. Together with the
initial momentum setup of $p_{x} = p_{y} = 0$ and $p_{z} =p_{\rm beam}$ for
each nucleon, a list containing the initial state of all nucleons in a given
AA collision is constructed. A collision happened between two nucleons
from different nuclei if their relative transverse distance is less than or
equal to the minimum approaching distance:
$D\leq\sqrt{\sigma_{\rm NN}^{\rm tot}/\pi}$. The collision time is calculated
with the assumption of straight-line trajectories. All such nucleon pairs
compose an NN collision time list. An NN collision with least collision time is
selected from the list and executed by \textsc{Pythia} (PYEVNW subroutine) with
the hadronization temporarily turned-off, as well as the strings and diquarks
broken-up. The nucleon list and NN collision time list are then updated. A new
NN collision with least collision time is selected from the updated NN collision
time list and executed by \textsc{Pythia}. With repeating the aforementioned
steps till the NN collision list empty, the initial partonic state is
constructed for a AA collision.

    Then, the partonic rescatterings are performed, where the LO-pQCD parton-parton
cross section~\cite{Combridge:1977dm,field} is employed. After partonic
rescattering, the string is recovered and then hadronized with the Lund string
fragmentation scheme resulting in an intermediate hadronic state. Finally, the
system proceeds into the hadronic rescattering stage and produces the final
hadronic state observed in the experiments.

    The PACIAE Monte-Carlo simulation provides a complete description of the NN
and AA collisions, which includes the partonic initialization, partonic
rescattering, hadronization, and the hadronic rescattering stages. Meanwhile,
a PACIAE model simulation could be selected to stop at any stage desired
conveniently. In this work, the simulations are stopped at the partonic
initialization stage, after the initial hard scattering and parton shower
processes happened.

    In this study, the $W^\pm$ relevant production channels are activated in a
user-controlled approach by setting ${\rm MSEL}=0$ in $\textsc{Pythia}$,
together with the following subprocesses switched-on:
\begin{align*}
    f_{i}\overline{f}_{j} & \rightarrow W^{+}/W^{-},   \\
    f_{i}\overline{f}_{j} & \rightarrow g W^{+}/W^{-}, \\
    f_{i}\overline{f}_{j} & \rightarrow \gamma W^{+}/W^{-}.
\end{align*}
In the above equations, $f$ refers to fermions (quarks) and its subscript
stands for the flavor code.

    In addition, the parton shower is switched off by setting MSTP(61)=0 and
MSTP(71)=0, to prevent the double-counting.
The PDF and other parameters are kept the same as the default \textsc{Pythia}.
Thus the PACIAE (\textsc{Pythia}) model simulation is corresponding to the NLO
pQCD in precision.

\section{\label{sec:res} Results and Discussions}

    In the left panel of Fig.~\ref{fig:comp}, we compare the PACIAE model
calculated electron and muon combined lepton charge asymmetry (black solid
squares) with ATLAS data~\cite{ATLAS:2019ibd} (red solid circles) for the
0-80\% centrality class in Pb+Pb collisions at $\sqrt{s_{\rm NN}}=5.02$~TeV.
The bottom sub-panel shows the differences between PACIAE calculation and ATLAS
data, where the shaded boxes show the total uncertainties in experimental data.
The kinematic cuts in PACIAE simulation is the same as the ones in ATLAS
experiments~\cite{ATLAS:2019ibd}, i.e.
\begin{align*}
    p^{l}_T > 25 {\rm~GeV}/c, & \hspace{0.4cm} p^{\nu}_T > 25 {\rm~GeV}/c, \\
    |\eta_l| < 2.5, & \hspace{0.4cm} m_T > 40 {\rm~GeV}/c^2,
\end{align*}
where $p^{l}_T$ and $p^{\nu}_T$ stand for the transverse momentum of the
lepton and the corresponding neutrino, respectively.

\begin{figure*}[htbp]
    \includegraphics[width=0.8\textwidth]{PACIAE_ATLAS_ALICE.eps}
    \caption{\label{fig:comp} Left panel: the combined electron and muon
charge asymmetry as a function of absolute pseudorapidity for the 0-80\%
centrality class in Pb+Pb collisions at $\sqrt{s_{\rm NN}}=5.02$~TeV. The black
solid squares with vertical error bar (hard to see) show the results of PACIAE
simulations. The red solid circles with vertical error bars and shaded error
boxes represent the ATLAS data ~\cite{ATLAS:2019ibd}. The lower sub-panel shows
the differences between the PACIAE calculations and the ATLAS data, with the
shaded boxes representing the total experimental uncertainties. The
PACIAE result points are shifted horizontally for better visibility.
Right panel: the muon charge asymmetry for the 0-90\% centrality class in
Pb+Pb collisions at $\sqrt{s_{\rm NN}}=5.02$~TeV, with kinematic cuts of
$p^{\mu}_T > 10$~GeV/$c$ and $2.5 < y_{\mu} < 4.0$. Here the open square with
vertical error bar shows the result of PACIAE simulation, and the blue solid
circle with vertical error bar and shaded error box represents ALICE
darum~\cite{ALICE:2022cxs}.}
\end{figure*}

\noindent $|\eta_l|$ refers to the
absolute pseudorapidity of the lepton, and $m_T$ is the transverse mass of the
lepton and neutrino system.
The pseudorapidity bins are divided as follows:
\begin{align*}
    |\eta_l| = &~(0,0.21),~~~~~(0.21,0.42), ~(0.42,0.63), ~(0.63,0.84), \\
               &~(0.84,1.05), ~(1.05,1.37), ~(1.37,1.52), ~(1.52,1.74), \\
               &~(1.74,1.95), ~(1.95,2.18), ~(2.18,2.5).
\end{align*}

    We see in this panel,
the PACIAE calculated lepton charge asymmetry is well consistent with ATLAS
data~\cite{ATLAS:2019ibd}, except the most forward $|\eta_l|$ bin. This
exception is also existed in the pQCD predictions calculated with the CT14 NLO
PDF set, even in the EPPS16 and nCTEQ15 nPDF sets~\cite{ATLAS:2019ibd}.

    The PACIAE calculation of the muon charge asymmetry (open square) for the
0-90\%  centrality class in Pb+Pb collisions at $\sqrt{s_{\rm NN}}=5.02$~TeV is
compared with the ALICE datum ~\cite{ALICE:2022cxs} (blue solid circle) in the
right panel of Fig.~\ref{fig:comp}. The same kinematic cuts of
$p^{\mu}_T > 10$~GeV/$c$ and $2.5 < y_{\mu} < 4$~\footnote[1]
{The ALICE muon spectrometer covers negative pseudorapidity. However, in
symmetric Pb+Pb collisions, positive values of rapidity are conventionally
used for the muon coverage.}
are set in both the experiment and theory. One can see in this panel that, the
calculated result agrees well with the ALICE datum.

    In Fig.~\ref{fig:Amu_pp}, we give the muon charge asymmetry as a function of
the valence quark number asymmetry in the elementary NN collisions of $pp$
(open triangle-ups), $pn$ (triangle-lefts), $np$ (triangle-rights), and $nn$
(triangle-downs)

\begin{figure}[htbp]
    \includegraphics[width=0.45\textwidth]{Amu_pp_vs_Av.eps}
    \caption{\label{fig:Amu_pp} The muon charge asymmetry as a function of the
valence quark asymmetry in the elementary NN collisions ($pp$, open
triangle-ups; $pn$, open triangle-lefts; $np$, open triangle-rights;
and $nn$, open triangle-downs) at $\sqrt{s}=5.02$~TeV for three rapidity
intervals of $|y|<1$, $2.5<|y|<4.0$, and full rapidity phase-space. The dashed,
dotted, and dash-dotted lines are the corresponding linear fitting to the $A_\mu$
vs. $A_v$ functions, respectively. The kinematic cut is $p^\mu_T > 10$~GeV/$c$.
The red short-dashed zero-baseline and the shaded area are added for better
distinguishment between the positive and negative values of $A_\mu$.}
\end{figure}

\noindent at $\sqrt{s}=5.02$~TeV. The kinematic cut of $p^\mu_T >
10$~GeV/$c$ is implemented in the simulations for three $y_\mu$ intervals:
$|y_\mu|<1$, $2.5<|y_\mu|<4.0$, and full $y_\mu$ phase-space.
The dashed, dotted and dash-dotted lines represent the linear fitting to the
function of $A_\mu$ vs. $A_v$ in three $y_\mu$ intervals above, respectively.
One can see all of $A_\mu$ in $pp$, $pn$, $np$ and $nn$ collisions are positive
in mid-rapidity ($|y_\mu| < 1$), while negative in the forward rapidity
($2.5 < |y_\mu| <4$). This sign-changing phenomenon has been also observed in
LHCb measurements~\cite{LHCb:2012lki,LHCb:2014liz,LHCb:2015mad,LHCb:2016zpq}.
It comes from the different helicity dependence of the lepton couplings to the
boson~\cite{LHCb:2012lki}. For full rapidity phase-space, $A_\mu$ in $pn$ and
$np$ collisions is consistent with zero owing to the vanishing relative numbers
between $u$- and $d$-valence quarks, while the $A_\mu$ in $pp$ and $nn$ are
almost equal in magnitude but opposite in sign. On the other hand, we can see
the slope fitted in mid-rapidity is gentler than the one in the forward
rapidity, due to the isospin effect in $W^{\pm}$ production is more pronounced
at forward rapidity~\cite{ConesaDelValle:2007dza}.

    The correlation between $A_\mu$ and $A_v$ in symmetrical AA collisions of Cu+Cu
(solid circles), Xe+Xe (solid squares), Au+Au (solid stars), Pb+Pb (solid
pentagons) and U+U (solid diamonds) at $\sqrt{s_{\rm NN}}=5.02$~TeV for three
rapidity intervals is shown in Fig.~\ref{fig:Amu_AA}, with the kinematic cut of
$p^\mu_T > 10$~GeV/$c$. The corresponding linear fitting is also given.
In this figure, the valence quark number asymmetry ($A_v$) are all negative,
since the number of neutrons (the number of $d$ quarks) is always greater than
that of protons ($u$ quarks) in heavy nuclei. The $A_\mu$ is positive
in the mid-rapidity bin for all five AA collision systems but negative
in the forward and full rapidity regions. It can be attributed to the $V-A$
decays of $W$ boson. This structure also affects the slope of the linear
fitting. With the ratio of neutron to proton ($A_v$ magnitude) increasing, the
fraction of $W^- \rightarrow l^- \nu_l$ becomes larger in forward rapidity
resulting a steeper slope than the one in mid-rapidity~\cite{LHCb:2012lki,ConesaDelValle:2007dza}.
It is more clear in Fig.~\ref{fig:Amu_AApp}, where we draw the $A_\mu$ vs. $A_v$
function in all the NN and AA collisions together. Here we see the fitted slope
in the full rapidity is larger than the one in mid-rapidity, but is smaller than
that in forward rapidity. One more interesting feature is that, the linear
scaling behavior of $A_\mu$ vs. $A_v$ shown in Fig.~\ref{fig:Amu_pp} is also
observed in Fig.~\ref{fig:Amu_AA} and~\ref{fig:Amu_AApp}.

\begin{figure}[htbp]
    \includegraphics[width=0.45\textwidth]{Amu_AA_vs_Av.eps}
    \caption{\label{fig:Amu_AA} The muon charge asymmetry as a function of
the valence quark number asymmetry in the symmetric AA collisions (Cu+Cu, solid
circles; Xe+Xe, solid squares; Au+Au, solid stars; Pb+Pb, solid pentagons; and
U+U, solid diamonds) at $\sqrt{s_{\rm NN}}=5.02$~TeV for three rapidity
intervals of $|y_\mu|<1$, $2.5<|y_\mu|<4.0$, and full rapidity phase-space. The
dashed, dotted and dash-dotted lines are the corresponding linear fitting to the
$A_\mu$ vs. $A_v$ functions, respectively. The kinematic cut is
$p^\mu_T > 10$~GeV/$c$. The red short-dashed zero-baseline and the shaded area
are added for better distinguishment between the positive and negative values
of $A_\mu$.}

\end{figure}

\begin{figure}[htbp]
    \includegraphics[width=0.45\textwidth]{Amu_AApp_vs_Av.eps}
    \caption{\label{fig:Amu_AApp} The muon charge asymmetry as a function of
the valence quark number asymmetry in the elementary NN collisions and the symmetric
AA collisions at $\sqrt{s_{\rm NN}}=5.02$~TeV, obtained by putting
Fig.~\ref{fig:Amu_pp} and \ref{fig:Amu_AA} together.}
\end{figure}

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

    In summary, we studied the lepton charge asymmetry from the $W^\pm$ decays in
both the NN and AA nuclear collisions at $\sqrt{s_{\rm NN}}=5.02$~TeV in frame
of the Monte-Carlo simulation model PACIAE. The simulated results well
reproduce the corresponding ATLAS and ALICE data. Furthermore, the correlation
between the lepton charge asymmetry $A_\mu$ and the valence quark number
asymmetry $A_v$ is studied in the above nuclear collisions for three rapidity
intervals: $|y_\mu|<1$, $2.5 < |y_\mu| < 4$, and the full rapidity phase space.
In each rapidity interval, a linear scaling behavior of $A_\mu$ vs. $A_v$
function is observed in the elementary NN collision systems
($pp$, $pn$, $np$, and $nn$), the AA collision systems
(Cu+Cu, Xe+Xe, Au+Au, Pb+Pb, and U+U), and even in their combination.

    The universal linear scaling behavior between $A_\mu$ and $A_v$ explored in
both nucleon-nucleon and nucleus-nucleus collision systems may provide an
additional constraint in the extraction of PDF (nPDF). On the contrary, we
also note that the nuclear modifications on the vacuum nucleon PDF change the
in-going parton kinematics in the hard scattering processes and consequently
impact the $p_T$- and $y$-distribution of the out-going particles. This effect
should result in a breaking of the $A_\mu$ vs. $A_v$ scaling established if
looking into a given lepton $p_T$--$y$ phase space. Therefore, by measuring
$p_T$ and $y$ differentially, the deviation of the $A_\mu$ vs. $A_v$
correlations between NN and AA collision systems should provide a new tool to
constrain the nPDF in a wide Bjorken-$x$ region. This has to be studied
further.

\begin{acknowledgments}
    The authors thank helpful discussion with Ming-Rui Zhao. This work was
supported by the National Natural Science Foundation of China (11775094,
11805079, 11905188, 11775313, 11905163), the Continuous Basic Scientific
Research Project (No.WDJC-2019-16) in CIAE, National Key Research and
Development Project (2018YFE0104800) and by the 111 project of the foreign
expert bureau of China.
\end{acknowledgments}