# On self-similar solutions of a multi-phase Stefan problem

Evgeny Yu. Panov
Yaroslav-the-Wise Novgorod State University,
Veliky Novgorod, Russian Federation

###### Abstract

We study self-similar solutions of a multi-phase Stefan problem, first in the
case of one space variable, and then in the radial multidimensional case. In
both these cases we prove that a nonlinear algebraic system for determination
of the free boundaries is gradient one and the corresponding potential is an
explicitly written coercive function. Therefore, there exists a minimum point
of the potential, coordinates of this point determine free boundaries and
provide the desired solution. Moreover, in one-dimensional case the potential
is proved to be strictly convex and this implies the uniqueness of the
solution. In contrary, in the multidimensional case the potential is not
convex but the uniqueness of our solution remains true, it follows from the
general theory. Bibliography$:$ $3$ titles.

##  1 Stefan problem. One-dimensional case

In a half-plane $\Pi=\\{\ (t,x)\ |\ t>0,x\in{\mathbb{R}}\ \\}$ we consider the
multi-phase Stefan problem for the heat equation

$u_{t}=a_{i}^{2}u_{xx},\quad u_{i}<u<u_{i+1},$ (1.1)

where $u_{-}=u_{0}<u_{1}<\cdots<u_{m}<u_{m+1}=u_{+}$, $u_{i}$, $i=1,\dots,m$
being the temperatures of phase transitions, $a_{i}>0$, $i=0,\dots,m$, are the
diffusivity constants. On the unknown lines $x=x_{i}(t)$ of phase transitions
where $u=u_{i}$ the following Stefan condition

$d_{i}x_{i}^{\prime}(t)+k_{i}u_{x}(t,x_{i}(t)+)-k_{i-1}u_{x}(t,x_{i}(t)-)=0$
(1.2)

is postulated, where $k_{i}>0$ is the thermal conductivity of the $i$-th
phase, while $d_{i}\geq 0$ is the Stefan number (the latent specific heat) for
the $i$-th phase transition. In (1.2) the unilateral limits
$u_{x}(t,x_{i}(t)+)$, $u_{x}(t,x_{i}(t)-)$ on the line $x=x_{i}(t)$ are taken
from the domain corresponding to the warmer/colder phase, respectively. By the
physical reasons, the Stefan numbers $d_{i}$ should be positive. We will study
even more general case $d_{i}\geq 0$. In this case the problem (1.1), (1.2) is
well-posed for $u_{-}<u<u_{+}$ and reduces to a degenerate nonlinear diffusion
equation, see [2, Chapter 5]. We will study the Cauchy problem with the
Riemann initial data

$u(0,x)=\left\\{\begin{array}[]{lr}u_{-},&x<0,\\\
u_{+},&x>0.\end{array}\right.$ (1.3)

By the invariance of our problem under the transformation group
$(t,x)\to(\lambda^{2}t,\lambda x)$, $\lambda\in{\mathbb{R}}$, $\lambda\not=0$,
it is natural to seek a self-similar solution of problem (1.1), (1.2), (1.3),
which has the form $u(t,x)=v(\xi)$, $\xi=x/\sqrt{t}$. For the heat equation
$u_{t}=a^{2}u_{xx}$ a self-similar solution must satisfy the linear ODE
$a^{2}v^{\prime\prime}=-\xi v^{\prime}/2$, the general solution of which is

$v=C_{1}F(\xi/a)+C_{2},\ C_{1},C_{2}=\mathrm{const},\mbox{ where
}F(\xi)=\frac{1}{2\sqrt{\pi}}\int_{-\infty}^{\xi}e^{-s^{2}/4}ds.$

This allows to write our solution in the form

$\displaystyle
v(\xi)=u_{i}+\frac{u_{i+1}-u_{i}}{F(\xi_{i+1}/a_{i})-F(\xi_{i}/a_{i})}(F(\xi/a_{i})-F(\xi_{i}/a_{i})),$
(1.4) $\displaystyle\xi_{i}<\xi<\xi_{i+1},\ i=0,\ldots,m,$

where $-\infty=\xi_{0}<\xi_{1}<\cdots<\xi_{m}<\xi_{m+1}=+\infty$ and we agree
that $F(-\infty)=0$, $F(+\infty)=1$. The parabolas $\xi=\xi_{i}$,
$i=1,\ldots,m$, where $u=u_{i}$, are free boundaries. They must be determined
by conditions (1.2). In the variable $\xi$ these conditions have the form (cf.
[3, Chapter XI])

$d_{i}\xi_{i}/2+\frac{k_{i}(u_{i+1}-u_{i})F^{\prime}(\xi_{i}/a_{i})}{a_{i}(F(\xi_{i+1}/a_{i})-F(\xi_{i}/a_{i}))}-\frac{k_{i-1}(u_{i}-u_{i-1})F^{\prime}(\xi_{i}/a_{i-1})}{a_{i-1}(F(\xi_{i}/a_{i-1})-F(\xi_{i-1}/a_{i-1}))}=0,$
(1.5)

$i=1,\ldots,m$. To investigate this nonlinear system, we notice that it is a
gradient one and coincides with the equality $\nabla E(\bar{\xi})=0$, where
the function

$\displaystyle
E(\bar{\xi})=-\sum_{i=0}^{m}k_{i}(u_{i+1}-u_{i})\ln(F(\xi_{i+1}/a_{i})-F(\xi_{i}/a_{i}))+\sum_{i=1}^{m}d_{i}\xi_{i}^{2}/4,$
(1.6) $\displaystyle\bar{\xi}=(\xi_{1},\ldots,\xi_{m})\in\Omega,$

the open convex domain $\Omega\subset{\mathbb{R}}^{m}$ is given by the
inequalities $\xi_{1}<\cdots<\xi_{m}$. Observe that $E(\bar{\xi})\in
C^{\infty}(\Omega)$. Since the function $F(x)$ takes values in the interval
$(0,1)$ all the terms in expression (1.6) are nonnegative while some of them
are strictly positive. Therefore, $E(\bar{\xi})>0$.

###  1.1 Coercivity of $E$

Let us introduce the sub-level sets

$\Omega_{c}=\\{\ \bar{\xi}\in\Omega\ |\ E(\bar{\xi})\leq c\ \\},\quad c>0.$

######  Theorem 1.1 (coercivity).

The sets $\Omega_{c}$ are compact for each $c>0$. In particular, the function
$E(\bar{\xi})$ reaches its minimal value.

###### Proof.

If $\bar{\xi}=(\xi_{1},\ldots,\xi_{m})\in\Omega_{c}$ then

$-k_{i}(u_{i+1}-u_{i})\ln(F(\xi_{i+1}/a_{i})-F(\xi_{i}/a_{i}))\leq
E(\bar{\xi})\leq c,\quad i=0,\ldots,m.$ (1.7)

It follows from (1.7) with $i=0$ that $F(\xi_{1}/a_{0})\geq
e^{-c/(k_{0}(u_{1}-u_{0}))}$, which implies the low bound $\xi_{1}\geq
r_{1}=a_{0}F^{-1}(e^{-c/(k_{0}(u_{1}-u_{0}))})$. Similarly, we derive from
(1.7) with $i=m$ and from the identity $1-F(\xi_{m}/a_{m})=F(-\xi_{m}/a_{m})$
that $\xi_{m}\leq r_{2}=-a_{m}F^{-1}(e^{-c/(k_{m}(u_{m+1}-u_{m}))})$. Since
all other coordinates of $\bar{\xi}$ are situated between $\xi_{1}$ and
$\xi_{m}$ we conclude that

$|\bar{\xi}|_{\infty}\doteq\max_{i=1,\ldots,m}|\xi_{i}|\leq
r=\max(|r_{1}|,|r_{2}|)$

and the set $\Omega_{c}$ is bounded.

Further, it follows from (1.7) that for all $i=1,\ldots,m-1$

$F(\xi_{i+1}/a_{i})-F(\xi_{i}/a_{i})\geq\delta\doteq\exp(-c/\alpha)>0,$ (1.8)

where $\displaystyle\alpha=\min_{i=1,\ldots,m-1}k_{i}(u_{i+1}-u_{i})>0$. Since
$F^{\prime}(\xi)=\frac{1}{2\sqrt{\pi}}e^{-\xi^{2}/4}<1$, the function $F(\xi)$
is Lipschitz with constant $1$, and it follows from (1.8) that

$(\xi_{i+1}-\xi_{i})/a_{i}\geq
F(\xi_{i+1}/a_{i})-F(\xi_{i}/a_{i})\geq\delta,\quad i=1,\ldots,m-1,$

and we obtain the estimates $\xi_{i+1}-\xi_{i}\geq\delta_{1}=\delta\min
a_{i}$. Thus, the set $\Omega_{c}$ is contained in a compact

$K=\\{\ \bar{\xi}=(\xi_{1},\ldots,\xi_{m})\in{\mathbb{R}}^{m}\ |\
|\bar{\xi}|_{\infty}\leq r,\ \xi_{i+1}-\xi_{i}\geq\delta_{1}\ \forall
i=1,\ldots,m-1\ \\}.$

Since $E(\bar{\xi})$ is continuous on $K$, the set $\Omega_{c}$ is a closed
subset of $K$ and therefore is compact. For $c>N\doteq\inf E(\bar{\xi})$, this
set is not empty and the function $E(\bar{\xi})$ reaches on it a minimal
value, which is evidently equal $N$. ∎

We have established the existence of minimal value $E(\bar{\xi}_{0})=\min
E(\bar{\xi})$. At the point $\bar{\xi}_{0}$ the required condition $\nabla
E(\bar{\xi}_{0})=0$ is satisfied, and $\bar{\xi}_{0}$ is a solution of system
(1.5). The coordinates of $\bar{\xi}_{0}$ determine the solution (1.4) of our
Stefan-Riemann problem. Thus, we establish the following existence result.

######  Theorem 1.2.

There exists a self-similar solution (1.4) of problem (1.1), (1.2), (1.3).

###  1.2 Convexity of the Function $E$ and Uniqueness of the Solution

In this section we prove that the function $E(\bar{\xi})$ is strictly convex.
Since a strictly convex function can have at most one critical point (and it
is necessarily a global minimum) the system (1.5) has at most one solution,
that is, a self-similar solution (1.4) of problem (1.1), (1.2), (1.3) is
unique. We will need the following simple lemma.

######  Lemma 1.1.

The function $P(x,y)=-\ln(F(x)-F(y))$ is strictly convex in the half-plane
$x>y$.

###### Proof.

The function $P(x,y)$ is infinitely differentiable in the domain $x>y$. To
prove the lemma, we need to establish that the Hessian $D^{2}P$ is positive
definite at every point. By the direct computation we find

$\displaystyle\frac{\partial^{2}}{\partial
x^{2}}P(x,y)=\frac{(F^{\prime}(x))^{2}-F^{\prime\prime}(x)(F(x)-F(y))}{(F(x)-F(y))^{2}},$
$\displaystyle\frac{\partial^{2}}{\partial
y^{2}}P(x,y)=\frac{(F^{\prime}(y))^{2}-F^{\prime\prime}(y)(F(y)-F(x))}{(F(x)-F(y))^{2}},\
\frac{\partial^{2}}{\partial x\partial
y}P(x,y)=-\frac{F^{\prime}(x)F^{\prime}(y)}{(F(x)-F(y))^{2}}.$

We have to prove positive definiteness of the matrix
$Q=(F(x)-F(y))^{2}D^{2}P(x,y)$ with the components

$\displaystyle Q_{11}=(F^{\prime}(x))^{2}-F^{\prime\prime}(x)(F(x)-F(y)),$
$\displaystyle Q_{22}=(F^{\prime}(y))^{2}-F^{\prime\prime}(y)(F(y)-F(x)),\
Q_{12}=Q_{21}=-F^{\prime}(x)F^{\prime}(y).$

Since $F^{\prime}(x)=e^{-x^{2}/4}$, then
$F^{\prime\prime}(x)=-\frac{x}{2}F^{\prime}(x)$ and the diagonal elements of
this matrix can be written in the form

$\displaystyle Q_{11}=F^{\prime}(x)(\frac{x}{2}(F(x)-F(y))+F^{\prime}(x))=$
$\displaystyle
F^{\prime}(x)(\frac{x}{2}(F(x)-F(y))+(F^{\prime}(x)-F^{\prime}(y)))+F^{\prime}(x)F^{\prime}(y),$
$\displaystyle
Q_{22}=F^{\prime}(y)(\frac{y}{2}(F(y)-F(x))+(F^{\prime}(y)-F^{\prime}(x)))+F^{\prime}(x)F^{\prime}(y).$

By Cauchy mean value theorem there exists such a value $z\in(y,x)$ that

$\frac{F^{\prime}(x)-F^{\prime}(y)}{F(x)-F(y)}=\frac{F^{\prime\prime}(z)}{F^{\prime}(z)}=-z/2.$

Therefore,

$\displaystyle
Q_{11}=F^{\prime}(x)(F(x)-F(y))(x-z)/2+F^{\prime}(x)F^{\prime}(y),$
$\displaystyle
Q_{22}=F^{\prime}(y)(F(x)-F(y))(z-y)/2+F^{\prime}(x)F^{\prime}(y),$

and it follows that $Q=R_{1}+F^{\prime}(x)F^{\prime}(y)R_{2}$, where $R_{1}$
is a diagonal matrix with the positive diagonal elements
$F^{\prime}(x)(F(x)-F(y))(x-z)/2$, $F^{\prime}(y)(F(x)-F(y))(z-y)/2$ while
$R_{2}=\left(\begin{smallmatrix}1&-1\\\ -1&1\end{smallmatrix}\right)$. Since
$R_{1}>0$, $R_{2}\geq 0$, then the matrix $Q>0$, as was to be proved. ∎

######  Remark 1.1.

In addition to Lemma 1.1 we observe that the functions $P(x,-\infty)$,
$P(+\infty,x)=P(-x,-\infty)$ of single variable $x$ are strictly convex. In
fact, it is sufficient to prove the strict convexity of the function
$P(x,-\infty)=-\ln F(x)$. By Lemma 1.1 in the limit as $y\to-\infty$ we obtain
that this function is convex, moreover,

$(F(x))^{2}\frac{d^{2}}{dx^{2}}P(x,-\infty)=F^{\prime}(x)(\frac{x}{2}F(x)+F^{\prime}(x))=\lim_{y\to-\infty}Q_{11}\geq
0.$

Since $F^{\prime}(x)>0$, we find, in particular, that
$\frac{x}{2}F(x)+F^{\prime}(x)\geq 0$. If $\frac{d^{2}}{dx^{2}}P(x,-\infty)=0$
at some point $x=x_{0}$ then $0=\frac{x_{0}}{2}F(x_{0})+F^{\prime}(x_{0})$ is
the minimum of the nonnegative function $\frac{x}{2}F(x)+F^{\prime}(x)$.
Therefore, its derivative $(\frac{x}{2}F+F^{\prime})^{\prime}(x_{0})=0$. Since
$F^{\prime\prime}(x)=-\frac{x}{2}F^{\prime}(x)$, this derivative

$(\frac{x}{2}F+F^{\prime})^{\prime}(x_{0})=F(x_{0})/2+\frac{x_{0}}{2}F^{\prime}(x_{0})+F^{\prime\prime}(x_{0})=F(x_{0})/2>0.$

But this contradicts our assumption. We conclude that
$\frac{d^{2}}{dx^{2}}P(x,-\infty)>0$ and the function $P(x,-\infty)$ is
strictly convex.

Now we are ready to prove the expected convexity of $E(\bar{\xi})$.

######  Theorem 1.3.

The function $E(\bar{\xi})$ is strictly convex.

###### Proof.

We introduce the functions

$E_{i}(\bar{\xi})=-k_{i}(u_{i+1}-u_{i})\ln(F(\xi_{i+1}/a_{i})-F(\xi_{i}/a_{i})),\quad
i=0,\ldots,m.$

By Lemma 1.1 and Remark 1.1 all these functions are convex. Since

$E(\bar{\xi})=\sum_{i=0}^{m-1}E_{i}(\bar{\xi})+E_{m}(\bar{\xi})+\sum_{i=1}^{m}d_{i}\xi_{i}^{2}/4$

and all the functions in this sum are convex, it is sufficient to prove the
strong convexity of the sum

$\tilde{E}(\bar{\xi})=\sum_{i=0}^{m-1}E_{i}(\bar{\xi}).$

By Lemma 1.1 and Remark 1.1 all the terms in this sum are convex function.
Therefore, the function $\tilde{E}$ is convex as well. To prove the strict
convexity, we assume that for some vector
$\zeta=(\zeta_{1},\ldots,\zeta_{m})\in{\mathbb{R}}^{m}$.

$D^{2}\tilde{E}(\bar{\xi})\zeta\cdot\zeta=\sum_{i,j=1}^{m}\frac{\partial^{2}\tilde{E}(\bar{\xi})}{\partial\xi_{i}\partial\xi_{j}}\zeta_{i}\zeta_{j}=0$
(1.9)

Since

$0=D^{2}\tilde{E}(\bar{\xi})\zeta\cdot\zeta=\sum_{i=0}^{m-1}D^{2}E_{i}(\bar{\xi})\zeta\cdot\zeta$

while all the terms are nonnegative, we conclude that

$D^{2}E_{i}(\bar{\xi})\zeta\cdot\zeta=0,\quad i=0,\ldots,m-1.$ (1.10)

By Lemma 1.1 for $i=1,\ldots,m-1$ the function $E_{i}(\bar{\xi})$ is strictly
convex as a function of two variables $\xi_{i},\xi_{i+1}$ and it follows from
(1.10) that $\zeta_{i}=\zeta_{i+1}=0$, $i=1,\ldots,m-1$. Observe that in the
case $m=1$ there are no such $i$. In this case we apply (1.10) for $i=0$.
Taking into account Remark 1.1, we find that $E_{0}(\bar{\xi})$ is a strictly
convex function of the single variable $\xi_{1}$, and it follows from (1.10)
that $\zeta_{1}=0$. In any case we obtain that the vector $\zeta=0$. Thus,
relation (1.9) can hold only for zero $\zeta$, that is, the matrix
$D^{2}\tilde{E}(\bar{\xi})$ is (strictly) positive definite, and the function
$\tilde{E}(\bar{\xi})$ is strictly convex. This completes the proof. ∎

Theorems 1.1,1.2 imply our main result.

######  Theorem 1.4.

There exists a unique self-similar solution (1.4) of problem (1.1), (1.2),
(1.3), and it correspond to the minimum of strictly convex and coercive
function (1.6).

######  Remark 1.2.

In recent paper [4] the problem (1.1), (1.2), (1.3) was studied in the case of
arbitrary (possibly negative) latent specific heats $d_{i}$. It was found a
necessary and sufficient condition for the coercivity of $E(\bar{\xi})$, as
well as a stronger sufficient condition of its strict convexity.

##  2 Multidimensional Stefan problem

Now we consider the multidimensional case, $x\in{\mathbb{R}}^{n}$, $n>1$. In
the half-space $\Pi=\\{\ (t,x)\ |\ t>0,x\in{\mathbb{R}}^{n}\ \\}$ we study the
Stefan problem

$u_{t}=a_{i}^{2}\Delta u,\quad u_{i}<u<u_{i+1},\ i=0,\ldots,m,$ (2.1)

where $u_{0}<u_{1}<\cdots<u_{m}<u_{m+1}=+\infty$. On a phase transition hyper-
surface $S$, where $u=u_{i}$, $i=1,\ldots,m$, the Stefan condition reads

$-d_{i}\nu_{t}+(k_{i}\nabla_{x}u_{+}-k_{i-1}\nabla_{x}u_{-})\cdot\nu_{x}=0,$
(2.2)

where $\nu=(\nu_{t},\nu_{x})\in{\mathbb{R}}\times{\mathbb{R}}^{n}$,
$\nu\not=0$ is a normal vector to $S$ directed from the $(i-1)$-th phase to
the $i$-th phase (so that the temperature increases in the direction of
$\nu$), $\nabla_{x}u_{-}$, $\nabla_{x}u_{+}$ are unilateral limits of the
gradient $\nabla_{x}u$ from the colder (respectively, the warmer) side. As in
section 1, the constants $k_{i}>0$, $i=0,\ldots,m$, are the thermal
conductivities, $d_{i}\geq 0$, $i=1,\ldots,m$, are the latent specific heats.

We are interesting in self-similar radial solutions

$u=v(r/\sqrt{t}),\ r=|x|=(x_{1}^{2}+\cdots+x_{n}^{2})^{1/2},$

defined for $r>0$. Putting the function $v(r/\sqrt{t})$ into the heat equation
$u_{t}=a^{2}\Delta u$, we find that this equation reduces to the linear ODE

$a^{2}v^{\prime\prime}+\left(a^{2}\frac{n-1}{\xi}+\frac{\xi}{2}\right)v^{\prime}=0.\quad\xi=r/\sqrt{t}.$

The general solution of this equation is the following:

$v=C_{1}G(\xi/a)+C_{2},\quad C_{1},C_{2}=\mathrm{const},\quad
G(y)=\int_{y}^{+\infty}s^{1-n}e^{-s^{2}/4}ds,\ y>0.$ (2.3)

Remark that

$G(+\infty)\doteq\lim_{y\to+\infty}G(y)=0,\ \lim_{y\to 0+}G(y)=+\infty.$

Moreover, as $y\to 0+$

$G(y)\sim{\mathcal{E}}_{n}(y)=\left\\{\begin{array}[]{lcr}\frac{1}{(n-2)y^{n-2}}&,&n\geq
3,\\\ -\ln y&,&n=2,\end{array}\right.$ (2.4)

so that $-\frac{1}{\omega_{n}}{\mathcal{E}}_{n}(|x|)$ is the fundamental
solution of the Laplace operator in ${\mathbb{R}}^{n}$, $\omega_{n}$ being the
surface area of a unit sphere in ${\mathbb{R}}^{n}$. We are going to construct
a solution to problem (2.1), (2.2) with the constant initial data

$u(0,x)=u_{0},$ (2.5)

and with the following growth condition at $x=0$:

$u(t,x)\sim A{\mathcal{E}}_{n}(|x|/(a_{m}\sqrt{t}))\ \mbox{ as }x\to 0,$ (2.6)

where $A$ is a positive constant (we will demonstrate below that this
condition corresponds to a positive heat source at the point $x=0$). It
follows from (2.3), (2.4), (2.6) that in a vicinity of zero our solution has
the form

$v(\xi)=u_{m}+A(G(\xi/a_{m})-G(\xi_{m}/a_{m})),\quad 0<\xi<\xi_{m},$ (2.7)

where $\xi_{m}$ corresponds to the unknown phase transition surface
$|x|=\xi_{m}\sqrt{t}$ with temperature $u_{m}$. Notice that the function
$u=v(|x|/\sqrt{t})$ satisfies the inhomogeneous heat equation

$u_{t}-a_{m}^{2}\Delta u=A\omega_{n}a_{m}^{n}t^{\frac{n}{2}-1}\delta(x)$ (2.8)

in the sense of distribution on $\Pi$ (in ${\mathcal{D}}^{\prime}(\Pi)$), with
the positive point heat source
$A\omega_{n}a_{m}^{n}t^{\frac{n}{2}-1}\delta(x)$ (in the case $n=2$ it does
not depend on time). Here, as usual, $\delta(x)$ denotes the Dirac
$\delta$-function. In fact, applying the distribution $u_{t}-a_{m}^{2}\Delta
u$ to a test function $f=f(t,x)\in C_{0}^{\infty}(\Pi)$ and integrating by
parts in the domain $|x|>\delta>0$, we obtain that

$\displaystyle-\int_{|x|>\delta}u[f_{t}+a_{m}^{2}\Delta
f]dtdx=\int_{|x|>\delta}[u_{t}-a_{m}^{2}\Delta u]fdtdx+$ $\displaystyle
a_{m}^{2}\int_{|x|=\delta}(u\nabla f-f\nabla u)\cdot\nu
dtd\sigma(x)=a_{m}^{2}\int_{|x|=\delta}(u\nabla f-f\nabla u)\cdot\nu
dtd\sigma(x),$ (2.9)

where $\nu=x/|x|$ is a unit normal vector on the sphere $|x|=\delta$, and
$d\sigma(x)$ is the surface Lebesgue measure on this sphere. Since for
$|x|=\delta$

$u=AG(\delta/(a_{m}\sqrt{t}))+\mathrm{const},\quad\nabla
u\cdot\nu=-A\delta^{1-n}(a_{m}\sqrt{t})^{n-2}e^{-\frac{\delta^{2}}{4a_{m}^{2}t}},$

then, taking into account (2.4), we find that in the limit as $\delta\to 0$

$\int_{|x|=\delta}(u\nabla f-f\nabla u)\cdot\nu d\sigma(x)\to
A\omega_{n}(a_{m}\sqrt{t})^{n-2}f(t,0).$

It now follows from (2) in the limit as $\delta\to 0$ that

$-\int_{\Pi}u[f_{t}+a_{m}^{2}\Delta
f]dtdx=A\omega_{n}a_{m}^{n}\int_{0}^{+\infty}t^{\frac{n}{2}-1}f(t,0)dt=\langle
A\omega_{n}a_{m}^{n}t^{\frac{n}{2}-1}\delta(x),f\rangle$

for all test functions $f\in C_{0}^{\infty}(\Pi)$. Hence, relation (2.8) holds
in ${\mathcal{D}}^{\prime}(\Pi)$.

Let $\xi_{i}$ be parameters of remaining phase transitions
$|x|=\xi_{i}\sqrt{t}$ corresponding to the temperatures $u_{i}$,
$i=1,\ldots,m-1$. Then $\xi_{m}<\xi_{m-1}<\cdots<\xi_{1}<\xi_{0}\doteq+\infty$
and in accordance with (2.3) for $i=0,\ldots,m-1$

$v(\xi)=u_{i}+\frac{(u_{i+1}-u_{i})(G(\xi/a_{i})-G(\xi_{i}/a_{i}))}{G(\xi_{i+1}/a_{i})-G(\xi_{i}/a_{i})},\quad\xi_{i+1}<\xi<\xi_{i},$
(2.10)

where the value $G(\xi_{0})=G(+\infty)=0$. The function $v(\xi)$, defined by
relations (2.7), (2.10), is a solution of our Stefan problem. Notice that this
function decreases from $+\infty$ at $\xi=0$ to $u_{0}$ at $\xi=+\infty$ and,
in particular, the initial condition (2.5) is satisfied. The unknown
parameters $\xi_{i}$ are determined by Stefan condition (2.2) on the surfaces
$|x|=\xi_{i}\sqrt{t}$, which can be written as

$d_{i}\xi_{i}/2+k_{i}v^{\prime}(\xi_{i}-0)-k_{i-1}v^{\prime}(\xi_{i}+0)=0,\
i=1,\ldots,m.$

Computing the derivatives $v^{\prime}(\xi_{i}\pm 0)$ from expressions (2.10),
(2.7), we arrive at the system

$\displaystyle
d_{i}\xi_{i}/2+k_{i}\frac{(u_{i+1}-u_{i})G^{\prime}(\xi_{i}/a_{i})}{a_{i}(G(\xi_{i+1}/a_{i})-G(\xi_{i}/a_{i}))}-$
$\displaystyle
k_{i-1}\frac{(u_{i}-u_{i-1})G^{\prime}(\xi_{i}/a_{i-1})}{a_{i-1}(G(\xi_{i}/a_{i-1})-G(\xi_{i-1}/a_{i-1}))}=0,\quad
i=1,\ldots,m-1,$ (2.11) $\displaystyle
d_{m}\xi_{m}/2+k_{m}\frac{A}{a_{m}}G^{\prime}(\xi_{m}/a_{m})-$ $\displaystyle
k_{m-1}\frac{(u_{m}-u_{m-1})G^{\prime}(\xi_{m}/a_{m-1})}{a_{m-1}(G(\xi_{m}/a_{m-1})-G(\xi_{m-1}/a_{m-1}))},\quad
i=m.$ (2.12)

Like in one-dimensional case this system turns out to be gradient one, it
coincides with the equality $\nabla E=0$, where the function

$\displaystyle
E(\bar{\xi})=-\sum_{i=0}^{m-1}k_{i}(u_{i+1}-u_{i})\ln(G(\xi_{i+1}/a_{i})-G(\xi_{i}/a_{i}))$
$\displaystyle+k_{m}AG(\xi_{m}/a_{m})+\frac{1}{4}\sum_{i=1}^{m}d_{i}\xi_{i}^{2},\quad\bar{\xi}=(\xi_{1},\ldots,\xi_{m})\in\Omega,$
(2.13)

where $\Omega$ is an open convex cone in ${\mathbb{R}}^{m}$ consisting of
vectors with strictly decreasing positive coordinates. We are going to
demonstrate that, like in the one-dimensional case, the function
$E(\bar{\xi})$ is coercive (but is not convex anymore). This implies the
existence of a point of its global minimum. Coordinates of this point
determine a solution (2.7), (2.10) of our problem.

We first prove that the function $E(\bar{\xi})$ is bounded from below.

######  Lemma 2.1.

There exists a constant $E_{0}$ such that $E(\bar{\xi})\geq E_{0}$ for all
$\bar{\xi}\in\Omega$.

###### Proof.

Since the function $G(y)$ decreases,

$G(\xi_{i+1}/a_{i})-G(\xi_{i}/a_{i})\leq G(\xi_{i+1}/a_{i})\leq
G(\xi_{m}/a_{i}),\quad i=0,\ldots,m-1.$

Taking also into account that $k_{i}(u_{i+1}-u_{i})>0$, $d_{i}\geq 0$, we
obtain that $E(\bar{\xi})\geq f(\xi_{m})$, where the function

$f(y)=-\sum_{i=0}^{m-1}k_{i}(u_{i+1}-u_{i})\ln G(y/a_{i})+k_{m}AG(y/a_{m})$
(2.14)

is continuous on $(0,+\infty)$. As follows from (2.4), the term
$k_{m}AG(y/a_{m})\to+\infty$ as $y\to 0+$ faster than all the terms $\ln
G(y/a_{i})$. Therefore, $\displaystyle\lim_{y\to 0+}f(y)=+\infty$. It is also
clear that $\displaystyle\lim_{y\to+\infty}f(y)=+\infty$. By these limit
relations we see that $f(y)$ has a global minimum. Taking
$E_{0}=\min\limits_{y>0}f(y)$, we complete the proof. ∎

Now we prove the coercivity of the function $E(\bar{\xi})$.

######  Theorem 2.1.

For each $c\in{\mathbb{R}}$ the set $\Omega_{c}=\\{\ \bar{\xi}\in\Omega\ |\
E(\bar{\xi})\leq c\ \\}$ is compact.

###### Proof.

If $\bar{\xi}=(\xi_{1},\ldots,\xi_{m})\in\Omega_{c}$ then $f(\xi_{m})\leq
E(\bar{\xi})\leq c$, where $f(y)$ is the function (2.14). Since
$f(y)\to+\infty$ as $y\to 0+$ there exists such $r_{1}>0$ that $f(y)>c$ for
all $y\in(0,r_{1})$. Then $\xi_{m}$ cannot be less than $r_{1}$:

$\xi_{m}\geq r_{1}.$ (2.15)

In particular, this implies that

$k_{i}(u_{i+1}-u_{i})\ln(G(\xi_{i+1}/a_{i})-G(\xi_{i}/a_{i}))\leq
M=\max_{i=0,\ldots,m-1}k_{i}(u_{i+1}-u_{i})\ln G(r_{1}/a_{i}).$

It follows from this bound that for each $i=0,\ldots,m-1$

$\displaystyle-
k_{i}(u_{i+1}-u_{i})\ln(G(\xi_{i+1}/a_{i})-G(\xi_{i}/a_{i}))\leq$
$\displaystyle E(\bar{\xi})+(m-1)M\leq c_{1}\doteq c+(m-1)M.$ (2.16)

Taking $i=0$, we find $-k_{0}(u_{1}-u_{0})\ln G(\xi_{1}/a_{0})\leq c_{1}.$
This implies that $G(\xi_{1}/a_{0})\geq e^{-c_{1}/(k_{0}(u_{1}-u_{0}))}>0$ and

$\xi_{1}\leq r_{2}\doteq a_{0}G^{-1}(e^{-c_{1}/(k_{0}(u_{1}-u_{0}))}).$ (2.17)

From relations (2) with $i=1,\ldots,m-1$ it follows that

$G(\xi_{i+1}/a_{i})-G(\xi_{i}/a_{i})\geq
q\doteq\min_{i=1,\ldots,m-1}e^{-c_{1}/(k_{i}(u_{i+1}-u_{i}))}>0$ (2.18)

Since $|G^{\prime}(y)|\leq y^{1-n}$, the function $G(y)$ is Lipschitz with
constant $L=\alpha^{1-n}$ on the interval $[\alpha,+\infty)$. Choosing such
$\alpha>0$ that $\alpha\max\limits_{i=1,\ldots,m-1}a_{i}\leq r_{1}$, so that
$\xi_{i}/a_{i}>\xi_{i+1}/a_{i}\geq r_{1}/a_{i}\geq\alpha$, we derive from
(2.18) the inequalities

$L(\xi_{i}-\xi_{i+1})/a_{i}\geq G(\xi_{i+1}/a_{i})-G(\xi_{i}/a_{i})\geq q,$

which implies the bounds

$\xi_{i}-\xi_{i+1}\geq\delta=\frac{q}{L}\min_{i=1,\ldots,m-1}a_{i}>0.$ (2.19)

In view of (2.15), (2.17), (2.19) we find that $\bar{\xi}\in K$, where

$K=\\{\ \bar{\xi}=(\xi_{1},\ldots,\xi_{m})\in{\mathbb{R}}^{m}\ |\
r_{2}\geq\xi_{1}\geq\cdots\geq\xi_{m}\geq r_{1}>0,\
\xi_{i}-\xi_{i+1}\geq\delta,\ i=1,\ldots,m-1\ \\}$

is a compact subset of $\Omega$. By the continuity of $E(\bar{\xi})$, we
conclude that the set $\Omega_{c}$ is a closed subset of the compact $K$ and,
therefore, is compact. ∎

As we have already demonstrated in the one-dimensional case, the coercivity of
$E(\bar{\xi})$ implies that there exists a global minimum
$E(\bar{\xi}_{0})=\min_{\Omega}E(\bar{\xi})$. Coordinates of the point
$\bar{\xi}_{0}$ provide the unknown parameters $\xi_{i}$, $i=1,\ldots,m$ of
the free boundaries in solution (2.7), (2.10). We establish the existence of
solution.

######  Theorem 2.2.

There exists a self-similar solution (2.7), (2.10) of the Stefan problem
(2.1), (2.2), (2.5), (2.6).

Observe that in the general case the function $E(\bar{\xi})$ is not convex. In
fact, assume that $m=1$, so that

$E(\bar{\xi})=E(\xi_{1})=-k_{0}(u_{1}-u_{0})\ln
G(\xi_{1}/a_{0})+Ak_{1}G(\xi_{1}/a_{1})+d_{1}\xi_{1}^{2}/4.$

Since the function

$f(\xi_{1})=-k_{0}(u_{1}-u_{0})\ln
G(\xi_{1}/a_{0})+d_{1}\xi_{1}^{2}/4\to-\infty\ \mbox{ as }\xi_{1}\to 0+,$

it cannot be convex on $(0,+\infty)$. This implies that
$E(\xi_{1})=f(\xi_{1})+Ak_{1}G(\xi_{1}/a_{1})$ is not convex either for all
sufficiently small $Ak_{1}>0$.

Nevertheless, a solution (2.7), (2.10) of Stefan problem (2.1), (2.2), (2.5),
(2.6) is unique. This can be proved like in [2, Chapter V, § 9]. For the sake
of completeness, we provide the details. First we notice that
$u=v(|x|/\sqrt{t})$ is a solution of (2.1), (2.2), (2.5), (2.6) if and only if
it is a weak solution to the problem

$\beta(u)_{t}-\Delta_{x}\alpha(u)=\omega_{n}k_{m}a_{m}^{n-2}t^{\frac{n}{2}-1}\delta(x),\quad
u(0,x)\equiv u_{0},$ (2.20)

where $\alpha(u)$, $\beta(u)$ are strictly increasing functions on
$(u_{0},+\infty)$ linear on each interval $(u_{i},u_{i+1})$, $i=0,\ldots,m$,
with slopes $\alpha^{\prime}(u)=k_{i}$, $\beta^{\prime}(u)=k_{i}/a_{i}^{2}$,
and such that

$\alpha(u_{i}+)-\alpha(u_{i}-)=0,\quad\beta(u_{i}+)-\beta(u_{i}-)=d_{i},\
i=1,\ldots,m.$

Remark that for all $u,v>u_{0}$, $u\not=v$

$0<\frac{\alpha(u)-\alpha(v)}{\beta(u)-\beta(v)}\leq\max a_{i}^{2}.$ (2.21)

Assume that $u_{i}=v_{i}(|x|/\sqrt{t})$, $i=1,2$, are two solutions. Then it
follows from (2.20) that

$(\beta(u_{1})-\beta(u_{2}))_{t}-\Delta_{x}(\alpha(u_{1})-\alpha(u_{2}))=0\
\mbox{ in }{\mathcal{D}}^{\prime}(\Pi)$ (2.22)

As follows from expressions (2.7), (2.10),
$P(t,x)\doteq\beta(v_{1}(\xi))-\beta(v_{2}(\xi))=\mathrm{const}$ for small
$\xi=|x|/\sqrt{t}>0$ and $|P(t,x)|\leq\mathrm{const}\cdot G(\xi/a_{0})$ for
large $\xi$. By the L’Hôpital’s rule

$\displaystyle\lim_{y\to+\infty}-\frac{2y^{-1}G^{\prime}(y)}{G(y)}=2\lim_{y\to+\infty}\frac{y^{-2}G^{\prime}(y)-y^{-1}G^{\prime\prime}(y)}{G^{\prime}(y)}=$
$\displaystyle 2\lim_{y\to+\infty}[y^{-2}+y^{-1}((n-1)y^{-1}+y/2)]=1,$

and, therefore, $G(y)\sim-2y^{-1}G^{\prime}(y)=2y^{-n}e^{-y^{2}/4}$ as
$y\to+\infty$. This implies that for large $\xi=|x|/\sqrt{t}$

$|P(t,x)|\leq\mathrm{const}\cdot\xi^{-n}e^{-\xi^{2}/(4a_{0}^{2})}.$

We conclude that the function $P$ is bounded and $P(t,\cdot)\in
L^{2}({\mathbb{R}}^{n})$, $\|P(t,\cdot)\|_{2}\leq\mathrm{const}\cdot t^{n/4}$.
The similar statements hold for the function $Q=\alpha(u_{1})-\alpha(u_{2})$
because $Q=CP$, where
$C=C(t,x)=(\alpha(u_{1})-\alpha(u_{2}))/(\beta(u_{1})-\beta(u_{2}))$ (if
$u_{1}=u_{2}$ we set $C=0$) is a nonnegative bounded function, in view of
(2.21).

Hence, we can apply (2.22) to a test function $f=f(t,x)$ from the Sobolev
space $W_{2}^{1,2}(\Pi_{T})$, $\Pi_{T}=(0,T)\times{\mathbb{R}}^{n}$ (so that
$f,f_{t},\nabla_{x}f,D_{x}^{2}f\in L^{2}(\Pi_{T})$) such that $f(T,x)=0$. As a
result, we obtain the relation

$\int_{\Pi_{T}}P(t,x)[f_{t}+C\Delta_{x}f]dtdx=0.$ (2.23)

Let $F(t,x)\in C_{0}^{1}(\Pi_{T})$, $\varepsilon>0$, and
$f^{\varepsilon}=f^{\varepsilon}(t,x)\in W_{2}^{1,2}(\Pi_{T})$ be a solution
of the backward Cauchy problem

$f_{t}+(C+\varepsilon)\Delta_{x}f=F,\quad f(T,x)=0.$ (2.24)

As is demonstrated in [2], such a solution exists and satisfies the estimate

$|\Delta_{x}f^{\varepsilon}|_{2}\leq C_{0}/\sqrt{\varepsilon},$ (2.25)

where $C_{0}$ is a constant independent of $\varepsilon$. It follows from
(2.23) with $f=f^{\varepsilon}$ that

$\int_{\Pi_{T}}P(t,x)F(t,x)dtdx=\varepsilon\int_{\Pi_{T}}P(t,x)\Delta_{x}f^{\varepsilon}(t,x)dtdx.$
(2.26)

By the Cauchy–Bunyakovsky inequality and (2.25)

$\left|\int_{\Pi_{T}}P(t,x)\Delta_{x}f^{\varepsilon}(t,x)dtdx\right|\leq\|P\|_{2}\|\Delta_{x}f^{\varepsilon}\|_{2}\leq
C_{0}\|P\|_{2}/\sqrt{\varepsilon},$

and the right-hand side of (2.26) vanishes as $\varepsilon\to 0$. We conclude
that

$\int_{\Pi_{T}}P(t,x)F(t,x)dtdx=0$

for all $F(t,x)\in C_{0}^{1}(\Pi_{T})$ and all $T>0$. This means that $P=0$
a.e. on $\Pi$, that is, $u_{1}=u_{2}$.

Thus, we have established the uniqueness.

######  Theorem 2.3.

A self-similar solution (2.7), (2.10) of the Stefan problem (2.1), (2.2),
(2.5), (2.6) is unique. In particular, the function $E(\bar{\xi})$ has only
one critical point, the point of its global minimum.

######  Remark 2.1.

In the case when $u_{0}$ coincides with temperature of phase transition and
the corresponding latent specific heat $d_{0}>0$ there exists another self-
similar solution $u_{1}=v_{1}(\xi)$ of problem (2.1), (2.2), (2.5), (2.6),
which equals the constant $u_{0}$ in the domain $\xi_{0}<\xi<+\infty$ and has
the form (2.7), (2.10) for $\xi<\xi_{0}$ (where $\xi_{0}$ is now finite). This
solution contains the extra phase transition surface $\xi=\xi_{0}$. It can be
easily verified that this new solution corresponds to the point of global
minimum of the function

$\displaystyle
E_{1}(\bar{\xi})=-\sum_{i=0}^{m-1}k_{i}(u_{i+1}-u_{i})\ln(G(\xi_{i+1}/a_{i})-G(\xi_{i}/a_{i}))$
$\displaystyle+k_{m}AG(\xi_{m}/a_{m})+\frac{1}{4}\sum_{i=0}^{m}d_{i}\xi_{i}^{2},\quad\bar{\xi}=(\xi_{0},\ldots,\xi_{m})\in\Omega,$
(2.27)

where the cone $\Omega\subset{\mathbb{R}}^{m+1}$ consists of points with
strictly decreasing positive coordinates. We notice that the only difference
between formulas (2.1) and (2) is the extra term $d_{0}\xi_{0}^{2}/4$. This
term guarantees coercivity of $E_{1}$ and the existence of solution $u_{1}$.
Moreover, the uniqueness of this solution (of the prescribed above form) can
be established as in the proof of Theorem 2.3.

Observe that both the functions $u,u_{1}$ are solutions of (2.20), where the
increasing functions $\alpha(u),\beta(u)$ are now defined up to $u_{0}$ and
satisfy the conditions $\alpha(u_{0})=\alpha(u_{0}+)$,
$\beta(u_{0})=\beta(u_{0}+)-d_{0}$. This seems surprising because of the
uniqueness statement of Theorem 2.3. But there is no contradiction here. In
fact, as is easy to verify,

$P=\beta(v(\xi))-\beta(v_{1}(\xi))\mathop{\to}_{\xi\to+\infty}d_{0}>0,$

therefore, $P(t,\cdot)\notin L^{2}({\mathbb{R}}^{n})$, and the reasoning used
in the proof of Theorem 2.3 is not applicable.

## References

  * [1]
  * [2] O. A. Ladyzhenskaya, V. A. Solonnikov and N. N. Ural’tseva, _Linear and Quasi-Linear Equations of Parabolic Type_ , AMS, Providence, 1968.
  * [3] H. S. Carslaw, J. C. Jaeger, _Conduction of heat in solids, Second edition_ , Oxford University Press, 1959.
  * [4] E. Yu. Panov, Solutions of an Ill-Posed Stefan Problem, _J. Math. Sci._ 274 (2023), no. 4, 534–543.


Naseer Ahmad Asif, PhD
Existence and Multiplicity Results for Systems of Singular Boundary Value Problems
– Monograph –


# Establishing limits for the monopole production in pp accelerators

B. M. Carlos<EMAIL_ADDRESS>M. B. Gay Ducati<EMAIL_ADDRESS>High Energy Physics Phenomenology Group, Physics Institute, Federal University
of Rio Grande do Sul,
Caixa Postal 15051, CEP 91501-970, Porto Alegre, RS, Brazil

(May 2021)

###### Abstract

The spin $1/2$ magnetic monopole pair production and the spin $0$ monopolium
production are studied in proton-proton collisions. The velocity-dependent
coupling model and an additional effective coupling with a magnetic moment
parameter are used to calculate the photon fusion and Drell Yan pair
production cross sections. The monopole mass is unknown, hence the mass range
employed here is based on the last results of ATLAS and MoEDAL experiments,
which set a minimal value around 2 TeV. The cross sections are calculated for
the LHC center of mass energies, its successor, the HE-LHC, and the future
collider FCC. The estimated mass limits for observation in the LHC and the
other accelerators are obtained and the advantages in using both effective
couplings and the monopolium as an indirect measure are discussed.

monopoles, pp collisions, effective coupling

##  I Introduction

One of the first formulations for a magnetic monopole was proposed by Dirac in
the years 1931[1] and 1948[2], and consisted of a point-like particle with
non-specified mass and spin, carrying a magnetic charge $g$ and a non-physical
string that carries the magnetic flux of the particle. This particle would be
responsible for bringing symmetry to Maxwell’s equations and explaining the
charge quantization through the Dirac quantization condition (in SI natural
units)

$eg=4\pi\frac{n}{2},$ (1)

where $e$ is the electron charge and $n$ a positive integer. New models for
this unknown particle were proposed later on, like the dyon, a particle with
both electric and magnetic charges theorized by Schwinger [3], and those in
electroweak [4] and unification theories [5][6][7][8][9], but even with the
support of several experiments[10][11] none of them were discovered yet. A
possible explanation for this lack of experimental evidence was first given in
[12], where it was assumed that due to the strong magnetic coupling the
magnetic poles would always be in a bound state, called monopolium. These
particles could have a much lower mass than the monopole alone and be produced
with higher rates in accelerators.

The large value of the magnetic charge $g$ has lead to the development of
effective models that alleviate the perturbative limits in order to make
approximated predictions. The simplest one assumes a velocity dependent
coupling, where the moving monopole is treated as an electric charge and
couples to the photon just like the electron. This idea was used to determine
the first limits on Dirac’s monopole mass[13][10], and has been followed since
then in some theoretical works [14],[15] and in the current experimental
search for monopoles in pp colliders [16][17][18]. A more recent work [19]
proposed the addition of a magnetic moment term to the usual velocity-
dependent coupling, and this new parameter increases the limits where
perturbative methods can be used.

Recent experiments dedicated to the monopole search [18],[16] have given lower
bounds on the monopole mass for different spin and charges, considering the
monopole production by photon fusion and Drell Yan in pp collisions. In this
work, the monopole pair production is studied considering the two processes
analyzed in the current experiments, assuming a Dirac’s monopole with spin
$1/2$ and employing the two effective coupling models cited above. The
monopolium production is investigated for the photon fusion process. The mass
bounds from recent results are considered for the estimates, with the
objective to establish upper bounds for production in the LHC and the future
colliders HE-LHC[20] and FCC[21].

In the next section the effective couplings are presented, followed by the
formalism used in the processes of monopole and monopolium production, and the
one for total cross sections in proton-proton collisions.

##  II The Monopole Production

###  II.1 Pair production

In Dirac’s formalism the coupling force between a monopole and an antimonopole
is about $10^{3}$ times stronger than the one between an electron and a
positron, which turns the usual perturbative methods of QED not directly
applicable to monopoles. Facing the scenario where no formal theory to treat
magnetic monopole interactions is fully developed, some effective couplings
have been proposed, like the minimalistic model where the coupling is given by

$\alpha_{m}=\frac{\beta^{2}g^{2}}{4\pi},$

with $\beta$ the monopole velocity. In this formulation the moving magnetic
monopole is treated as an electric charge in analogy with the behaviour of a
moving electron producing a magnetic field. With this coupling the monopole
can be considered in the place of an electron, or other charged lepton, in
many processes by the replacement $e\to g\beta$. The velocity dependence
decreases the production rates, but now the coupling $\alpha_{m}$ can be
perturbatively expanded in the limit $\beta\ll 1$, and more predictions can be
obtained.

An alternative to the velocity-dependent coupling is to add a consistent
magnetic moment term dependence[19]. As the monopole itself generates a
magnetic field, it will not be necessary an electromagnetic scattering to
generate a magnetic moment, like in the electron case. Also, due to the large
value of $g$, the monopole can be expected to have a great magnetic moment
that would be relevant already in a tree-level diagram. A parameter $\kappa$
can be defined in terms of the magnetic moment of the monopole

$\mathbf{\mu_{m}}=\frac{g\beta}{2m}2(1+2\hat{\kappa})\mathbf{S},$ (2)

where $\hat{\kappa}=\kappa m$ and $\mathbf{S}$ is the particle spin. Now, the
photon-monopole coupling will be proportional to $g\beta$ plus a term with
$\kappa$ dependence, and the perturbative expansions are valid in the limit
$\kappa\gg 1$, $\hat{\kappa}\gg 1$ and $\beta\ll 1$[19].

The results obtained with these two couplings are only indicative, for they
are based mostly in effective calculations. However, they can be used to make
estimations about whether the monopoles are detectable in the energies of
current accelerators or not.

##  III The Monopole Production

###  III.1 Pair production

Although in Dirac’s theory the monopole does not have a defined spin, here the
case $\hat{S}=1/2$ is going to be considered in order to build a symmetry with
the electron (for spin $0$ and $1$ monopoles, we refer to [22] and [19]). The
monopole-antimonopole pair production can be studied in two processes: photon
fusion (u and t channels) and Drell Yan, both in leading order to avoid the
issues that come up with the expansion of the magnetic coupling.

Due to the ultrarelativistic energies reached at LHC, the flux of photons
emitted by protons plays an important role in the production of many
particles, such as leptons or charged scalar particles that have been widely
studied [23] and used to investigate the Higgs boson production[24][25] and
particles beyond the Standard Model[26]. In collisions involving protons, the
Drell Yan mechanism[27] can also be used to investigate such particles, as it
has been indicated[28] that for lepton and Higgs production those cross
sections are $\sim 10^{2}$ times larger than those for photon fusion. However,
when considering the magnetic coupling, the photon fusion cross section has an
enhancement of $\sim g^{2}\beta^{2}$ as compared to Drell Yan, which makes
this process a more relevant candidate to consider for the monopole production
in virtue of the quantization condition.

The production cross sections can be obtained by replacing the electron charge
and the electromagnetic coupling by $g\beta$ and $\alpha_{m}$, respectively,
in the expressions for electron-positron production for photon fusion[28]

$\hat{\sigma}_{\gamma\gamma\to
m\bar{m}}(\hat{s})=\frac{\pi\beta^{5}}{4\alpha^{2}\hat{s}}\left[\frac{3-\beta^{4}}{2\beta}\ln\frac{1+\beta}{1-\beta}-(2-\beta^{2})\right]$
(3)

and for Drell Yan[14]

$\hat{\sigma}_{q\bar{q}\to
m\bar{m}}(\hat{s})=\frac{\pi\eta_{q}^{2}\beta^{3}}{12\hat{s}}\left(2-\frac{2}{3}\beta^{2}\right).$
(4)

Once including the magnetic term, the cross sections can be recalculated[19]
for photon fusion

$\begin{split}\hat{\sigma}^{\kappa}_{\gamma\gamma\rightarrow
m\bar{m}}(\hat{s},\kappa)&=\frac{\pi\alpha^{2}_{m}(\beta)}{3\hat{s}}\left\\{\ln\left(\frac{1-\beta}{1+\beta}\right)\left[\beta^{2}\kappa^{2}\hat{s}(3\beta^{2}\kappa^{2}\hat{s}-6\kappa^{2}\hat{s}+6)+6\beta^{4}\right.\right.\\\
&\left.-(36\beta^{2}-72\beta)\kappa\sqrt{(1-\beta^{2})\hat{s}}-9\kappa^{4}\hat{s}^{2}-60\kappa^{2}\hat{s}-18\right]\\\
&\left.-\beta\kappa^{2}\hat{s}(7\beta^{2}\kappa^{2}\hat{s}^{2}+15\kappa^{2}\hat{s}+132)+12\beta^{3}-24\beta-36\kappa\sqrt{(1-\beta^{2})\hat{s}}\right\\}\end{split}$
(5)

and Drell Yan

$\hat{\sigma}^{\kappa}_{q\bar{q}\to
m\bar{m}}(\hat{s},\kappa)=\frac{\pi\eta^{2}\beta^{3}}{18\hat{s}}\left[3-\beta^{2}-(2\beta^{2}-3)\kappa^{2}\hat{s}+6\kappa\sqrt{\hat{s}-\beta^{2}\hat{s}}\right].$
(6)

In the next section the formalism for the monopole-antimonopole bound state
will be presented.

###  III.2 Monopolium production

The monopolium is a theorized bound state between a monopole and an
antimonopole, first proposed in [12] as a possible relic of magnetic monopoles
produced in the early universe. Due to the large magnetic coupling, a
monopole-antimonopole pair would probably annihilate into a pair of
photons[15][29] or form a monopolium, which can have a lower mass and more
stability. The monopolium can assume two spin values, $0$ and $1$, but in this
work the spin $0$ monopolium is choosen to be studied in order to consider the
lowest energy case.

The binding energy is written in terms of the monopole and monopolium masses
$m$ and $M$, respectively, $M=2m+E_{binding}$, and to write the cross section
for the photon fusion production the knowledge of the interaction potential in
the pair is required. Considering the large coupling between the monopole and
its antiparticle, it is possible to argue that both have some spatial
extension[31], meaning that the interaction is non-singular when their
separation ($r$) goes to zero. This can be described with the potential[32]

$V(r)=-g^{2}\left(\frac{1-e^{-\mu r}}{r}\right),$ (7)

where $\mu=2m/g^{2}$ is a cut-off parameter that describes the interaction
when $r\rightarrow 0$. The wave function can then be obtained by solving the
Schrödinger equation with the interaction potential, and the cross section
takes the form as [33]

$\hat{\sigma}_{\gamma\gamma\to
M}(\hat{s})=\frac{2\sqrt{2}[R(R-1)]^{3/2}}{\alpha^{2}\epsilon^{6}M^{2}}\frac{\bar{\Gamma}_{M}(\epsilon^{2}-1)^{2}}{(\epsilon^{2}-1)^{2}+\bar{\Gamma}_{M}^{2}},$
(8)

where $R=2m/M$, $\epsilon=\sqrt{\hat{s}}/M$ and
$\bar{\Gamma}_{M}=\Gamma_{M}/M$ is the normalized decay rate of the
monopolium. This cross section again can be used for processes that involve
photon fluxes, like proton and nucleus collisions. However, it is not yet
understood how the monopolium can be directly detected, and some models for
the decay into two photons have been studied[15][29] and give similar cross
sections.

The collision of two ultrarelativistic protons can generate several processes
involving photons, quarks and gluons and has a cleaner sign than a collision
between a proton and a nucleus or two nucleus. Proton beams at LHC are also
more energetic than the nucleus ones and can produce heavy particles up to a
few TeV, which is the expected mass range of the monopole. A study of monopole
production by photon fusion in nucleus collisions can be seen in [34], and as
noted there, despite the enhancement of $Z^{2}$ due to the nucleus, the cross
sections are very small ($\sim 10^{-10}$ fb) for monopole masses in the range
$400$ to $1000$ GeV.

The usual formalism for proton-proton collisions can be done as presented in
[28], where the cross sections are written in a factorized form containing the
cross section of the subprocess and the necessary distribution functions and
photon fluxes. The total cross section for a photon fusion process is written
as the sum of the elastic, semielastic and inelastic contribution, while in
Drell Yan the cross section is given by only one term, summed over all quark
flavors. For the elastic photon flux of the proton it was used the expression
given in [35], and the equivalent photon spectrum of quarks for the non-
elastic terms and Drell Yan is given in [28]. For the proton structure
function, it was used the Cteq6-1L parametrization [36] with scale
$Q^{2}=\hat{s}/4$, to compare with the results in previous works[14][15][33].

###  III.3 Experimental search for monopoles

The current LHC experiments dedicated to the search of magnetic monopoles are
the MoEDAL[37] experiment and the search for HIPs (highly ionizing particles)
in ATLAS[18]. The last results from ATLAS[18] give a lower limit of $m\leq
2370$ GeV for a spin $1/2$ monopole with one or two units of magnetic charge
($n=1$ or $2$ in the quantization condition) produced via Drell-Yan,
considering that the monopoles couple to photon with $\alpha_{m}\sim g^{2}$.
The MoEDAL experiment evaluated the limits considering both Drell Yan and
photon fusion production, and obtained a lower mass limit of $m\leq 2420$
GeV[17] for the same type of coupling, and a limit of $m\leq 1760$ TeV using
the velocity dependent coupling. So far there is no prediction regarding the
magnetic moment dependent coupling. Taking these high limits into account, it
is very probable that the search for magnetic monopoles will continue in
future accelerators with higher collision energies and luminosities. For this
reason the simulations are extended to the HE-LHC[20] and FCC[21] colliders,
expected to begin their operations after 2035, and to the high luminosity LHC
(HL-LHC)[38], which will begin its operations in 2027 with the same energies
of the LHC and a higher luminosity.

In Table 1 the main parameters of the pp colliders considered in the
calculations are presented. The luminosity per year refers to the total
luminosity of allocated physics time in a year, and corresponds to $160$ days
for all accelerators.

Table 1: Main parameters[20][39] of the LHC, HL-LHC, HE-LHC and FCC colliders. The beam energy is given in TeV, the peak luminosity in $10^{-5}$ fb${}^{-1}/$s, and the luminosity per year in $fb^{-1}$. Parameters | LHC | HL-LHC | HE-LHC | FCC
---|---|---|---|---
Beam Energy | 14 | 14 | 27 | 100
Peak Luminosity | 1 | 5 | 16 | 5-30
Luminosity/year | 55 | 350 | 500 | 250-1000

##  IV Results and Comments

###  IV.1 Monopole Production

The cross sections for the monopole-antimonopole pair production by photon
fusion and Drell Yan ($\hat{\kappa}=0$), at center of mass energy
$\sqrt{s}=14$ TeV are presented in Fig. 1. The results corroborate with those
in [14], [15] and [19], meaning that the two-photon process has higher cross
sections and it is a good candidate for simulations in LHC and new
accelerators. For $m>5500$ GeV the Drell Yan cross section overcomes the
photon fusion, and this phenomenon also occurs in the HE-LHC and FCC
calculations for $m>10$ TeV and $m>40$ TeV, respectively, as shown in Fig. 2.
However, considering the luminosities in Table 1, the cross sections for these
masses are not relevant in any of the accelerators.

Figure 1: Monopole pair production via photon fusion and Drell Yan in pp
collisions with $\hat{\kappa}=0$ and $\sqrt{s}=14$ TeV.

Figure 2: Monopole pair production via photon fusion and Drell Yan in pp
collisions with $\hat{\kappa}=0$ at future accelerators.

In Table 2, the expected number of events per year for monopole production
considering both Drell Yan and photon fusion are displayed. It can be seen
that even for the LHC successor, the HL-LHC, monopoles with $m>3$ Tev would
probably not produce enough data to be confirmed. With the HE-LHC and FCC,
monopoles with masses $m<5$ TeV and $m<18$ TeV, respectively, would have a
higher probability to be detected.

Table 2: Number of events of monopole production (Drell Yan $+$ photon fusion) per year, for different monopole masses (in TeV). Mass | LHC | HL-LHC | HE-LHC | FCC
---|---|---|---|---
3 | $<10$ | $<40$ | $<3\cdot 10^{4}$ | $<2\cdot 10^{7}$
5 | $<2\cdot 10^{-5}$ | $<8\cdot 10^{-5}$ | $<150$ | $<2\cdot 10^{6}$
9 | 0 | 0 | $<4\cdot 10^{-4}$ | $<3\cdot 10^{4}$
20 | 0 | 0 | 0 | $<10$
30 | 0 | 0 | 0 | $<2\cdot 10^{-3}$

The total cross sections for pair production considering both photon fusion
and Drell Yan are compared in Fig. 3 for different values of the magnetic
moment parameter $\hat{\kappa}$, again for $\sqrt{s}=14$ TeV. The results in
[19] point out that the cross sections for the subprocesses and the total
cross sections increase with the parameter $\kappa$, and the same behavior is
achieved for the entire mass range. The cross section for $\hat{\kappa}=3$ is
up to $10^{2}$ times higher than the one for $\hat{\kappa}=0$ only in photon
fusion and around $10$ times higher in Drell Yan. It can then be concluded
that the magnetic moment dependence, besides providing more applicability to
the perturbation methods, can also increase the monopole detection chances.

Figure 3: Monopole pair production considering both photon fusion and Drell
Yan production in pp collisions for different values of $\hat{\kappa}$ and
$\sqrt{s}=14$ TeV.

###  IV.2 Monopolium production

The results for Monopolium production are in Fig. 4 and Fig. 5, and it can be
seen that the cross section decreases with a lower rate when the monopolium
mass is raised, compared to the monopole pair production. The production is
also increased for higher values of monopole mass, supporting the results in
[15], [33] and [34].

Figure 4: Monopolium production for different values of $M$ and $R=2m/M$ in pp
collisions with $\sqrt{s}=14$ TeV.

Figure 5: Monopolium production for fixed monopole masses in pp collisions at
different accelerators (FCC: straight line, HE-LHC: dotted line, HL-LHC:
dashed line).

Considering a minimum of 1 event per year, the limits of detection in the LHC
for monopolium production are $M<5$ TeV, for a fixed monopole mass of $m=3$
TeV. For the HE-LHC and FCC energies and luminosities this limit is close to a
maximum possible mass, $M=6$ TeV. An analysis of the possible decay channels
of the monopolium would improve the estimation of production and detection in
a next step.

##  V Conclusions

Based on these results, the production of magnetic monopoles under the
processes discussed here, makes it reasonable to expect that their
experimental evidence can take a while. The lack of a well built perturbative
theory to study the magnetic monopole in QED is still one of the greatest
difficulties in obtaining predictions for their interactions. Our results
confirm that the production by photon fusion is more relevant in the cross
section range that allows detection in the LHC and future accelerators.
Although also preliminary, the study of the magnetic moment term could lead to
more applicability and new results to increase the chances of detection of
such particles. The predictions for the mass range here obtained indicate that
the production of monopoles in this model are at least one order of magnitude
higher in comparison to the velocity dependent model.

The focus of this work is to restrict the limits of monopole production in
view of the accessible energies of the current and next generation
accelerators. Hence, it was opted to let GUT and other monopoles beyond the
Standard model, such as those in [44], aside in the presented predictions,
although they may be taken into consideration for future works. If the lower
bounds in future results given experiments in accelerators continue to grow,
this may indicate that one has to look for other possible monopole sources.

## Acknowledgments

The authors would like to thank the Brazilian funding agency Conselho Nacional
de Desenvolvimento Científico e Tecnológico (CNPq) that partially financed
this work.

## References

  * [1] P. A. M. Dirac, Proc. Roy. Soc. Lon. A 133, 60 (1931).
  * [2] P. A. M. Dirac, Phys. Rev. 74, 817 (1948).
  * [3] J. S. Schwinger, Science 165, 757 (1969).
  * [4] Y. Nambu, Nucl. Phys. B 130, 505 (1977).
  * [5] G. ’t Hooft, Nucl. Phys. B 79, 276 (1974).
  * [6] A. M. Polyakov, JETP Lett. 20, 194 (1974); Pisma Zh. Eksp. Teor. Fiz. 20, 430 (1974).
  * [7] Y. Nambu, Phys. Rev. D, 10, 4262 (1974).
  * [8] B. Julia and A. Zee, Phys. Rev. D 11, 2227 (1975).
  * [9] T. T. Wu and C. N. Yang, Phys. Rev. D 12, 3843 (1975).
  * [10] M. J. Mulhearn, Ph.D. thesis, Fermilab, 2004.
  * [11] A. Abulencia et al. (CDF Collaboration), Phys. Rev. Lett 96, 201801 (2006).
  * [12] C. T. Hill, Nucl. Phys. B 224, 469 (1983).
  * [13] G. R. Kalbfleisch, K. A. Milton, M. G. Strauss, L. Gamberg, E. H. Smith, and W. Luo, Phys. Rev. Lett 85, 5292 (2000).
  * [14] T. Dougall and S. D. Wick, Eur. Phys. J. A 39, 213 (2009).
  * [15] L. N. Epele, H. Fanchiotti, C. A. G. Canal, V. A. Mitsou and V. Vento, Eur. Phys. J. Plus 127, 60 (2012).
  * [16] B. Acharya et al. (MoEDAL Collaboration), Int. J. Mod. Phys. A 29, 1430050 (2014).
  * [17] B. Acharya et al. (MoEDAL Collaboration), Phys. Rev. Lett. 123, 021802 (2019).
  * [18] G. Aad et al. (ATLAS Collaboration), Phys. Rev. Lett. 124, 031802 (2020).
  * [19] S. Baines, N. E. Mavromatos, V. A. Mitsou, J. L. Pinfold and A. Santra, Eur. Phys. J. C 78, 966 (2018); 79, 166(E) (2019).
  * [20] A. Abada et al. (HE-LHC collaboration), Collider. Eur. Phys. J. Spec. Top. 228, 1109 (2019).
  * [21] A. Abada et al. (FCC collaboration), Eur. Phys. J. C 79, 474 (2019).
  * [22] Y. Kurochkin, I. Satsunkevich, D. Shoukavy, N. Rusakovich and Y. Kulchitsky, Mod. Phys. Lett. A 21, 2873 (2006).
  * [23] V. Budnev, I. Ginzburg, G. Meledin and V. Serbo, Phys. Rept. 15, 181 (1975).
  * [24] J. F. Gunion, H. E. Haber and G. L. Kane and S. Dawson, Front. Phys. 80, 1 (2000).
  * [25] M. B. Gay Ducati and G. G. Silveira, Phys. Rev. D 82, 073004 (2010).
  * [26] T. D. Lee, Phys. Rev. D 8, 1226 (1973); H. Georgi, Hadronic J. 1, 155 (1978).
  * [27] S.D. Drell and T.M. Yan, Phys. Rev. Lett. 25, 316 (1970).
  * [28] M. Drees, R. M. Godbole, M. Nowakowski and S. D. Rindani, Phys. Rev. D 50, 2335 (1994).
  * [29] H. Fanchiotti, C. A. García Canal and V. Vento, Int. J. Mod. Phys. A 32, 1750202 (2017).
  * [30] M. E. Peskin and D. V. Schroeder, An introduction to Quantum Field Theory (Westview, Boulder, 1995).
  * [31] L. I. Schiff, Phys. Rev. 160, 1257 (1967).
  * [32] L. N. Epele, H. Fanchiotti, C. A. García Canal and V. Vento, Eur. Phys. J. C 56, 87 (2008).
  * [33] L. N. Epele, H. Fanchiotti, C. A. García Canal and V. Vento, Eur. Phys. J. C 62, 587 (2009).
  * [34] J. T. Reis and W. K. Sauter, Phys. Rev. D 96, 075031 (2017).
  * [35] M. Drees and D. Zeppenfeld, Phys. Rev. D 39, 2536 (1989).
  * [36] J. Pumplin, D. R. Stump, J. Huston, H. L. Lai, P. M. Nadolsky and W. K. Tung, JHEP 0207, 012 (2002).
  * [37] B. Acharya et al. (MoEDAL Collaboration), Int. J. Mod. Phys. A 29, 1430050 (2014).
  * [38] High-Luminosity Large Hadron Collider (HL-LHC), Preliminary Design Report, edited by G. Apollinari, I. Béjar Alonso, O. Brüning and M. Lamont, L. Rossi, CERN-2015-005 (CERN, Geneva, 2015).
  * [39] M. Benedikt, D. Schulte and F. Zimmermann, Phys. Rev. ST Accel. Beams 18, 101002 (2015).
  * [40] I. K. Affleck and N. S. Manton, Nucl. Phys. B 194, 38 (1982).
  * [41] O. Gould and A. Rajantie, Phys. Rev. Lett. 119, 241601 (2017).
  * [42] O. Gould and A. Rajantie, Phys. Rev. D 96, 076002 (2017).
  * [43] A. Rajantie, Phil. Trans. Roy. Soc. Lond. A 377, 20190333 (2019).
  * [44] M. Eto, Y. Hamada and M. Nitta, Phys. Rev. D 102, no.10, 105018 (2020)


# Uniform quasi-convex optimisation via Extremum Seeking

Nicola Mimmo, Lorenzo Marconi and Giuseppe Notarstefano This research was
partially supported by the European Project ”AerIal RoBotic technologies for
professiOnal seaRch aNd rescuE” (AirBorne), Call: H2020, ICT-25-2016/17, Grant
Agreement no: 780960.The authors are with the Department of Electrical and
Information Engineering, ”Guglielmo Marconi”, University of Bologna, 40126
Bologna, Italy, e-mail: {nicola.mimmo2, lorenzo.marconi,
<EMAIL_ADDRESS>
###### Abstract

The paper deals with a well-known extremum seeking scheme by proving
uniformity properties with respect to the amplitudes of the dither signal and
of the cost function. Those properties are then used to show that the scheme
guarantees the global minimiser to be semi-global practically stable despite
the presence of local saddle points. To achieve these results, we analyse the
average system associated with the extremum seeking scheme via arguments based
on the Fourier series.

##  1 Introduction

The early researches on the Extremum Seeking (ES) date bake to 1920s [1] and
since then this strategy has been extensively exploited to solve several
optimisation problems in electronics [2], mechatronics [3], mechanics [4],
aerodynamics [5], thermohydraulics [6], and thermoacoustic [7]. Some of the
most popular ES schemes are those proposed in [8, 9], which represent the
subject of the proposed analysis, although a remarkable variety of schemes
were proposed, such as the adoption of an integral action in [10], the use of
a cost function’s parameter estimator [11, 12], the introduction of an
observer [13], the extension to fractional derivatives in [3], the use of a
predictor to compensate output delays in [14], the implementation of a Newton-
based algorithm avoiding the Hessian matrix inversion [15, 16], and the
concurrent use of a simplex-method to find the global minimiser [17].

All the methods, in a way or another, share the common philosophy of
perturbing the system subject to optimisation by means of the so-called dither
signal, which is periodic in most of the proposed solutions, to unveil in
which direction the associated cost function decreases (in the case of a
minimisation problem). The analysis tool that is typically adopted to prove
practical convergence to the minimiser exploits the averaging theory [18]. It
is shown that the average system asymptotically converges to the minimiser and
that the trajectories of the original system remain closed to the average ones
if a design parameter, namely $\gamma$, is kept sufficiently small [19, 20, 9,
21, 22].

For sake of completeness, it is worth mentioning the results in [23, 24] in
which the value of the cost function is directly used to perturb the phase of
the dither signal rather than to estimate the local gradient. Lie-derivative
arguments, instead of averaging techniques, are then adopted in the analysis.

Within the previous research context, this paper proposes two contributions.

First, We propose to study the average systems presented in [8] via Fourier
series arguments. This alternative approach allows to handle a class of non-
convex cost functions with a unique global minimiser. We show that the average
system trajectories converge to a neighbourhood of the minimiser. The size of
this neighbourhood is proportional to $\delta$, which is not required to be
small. As second contribution, it is shown that the addition of an high-pass
filter makes $\gamma$ independent (uniform) on the value of the cost function.
In fact, we prove that semiglobal and practical convergence to the minimiser
is achieved with the parameter $\gamma$ only depending on the Lipschitz
constant of the cost function in the domain of interest and not on its value.
This allows for a tuning of the algorithm that preserves good convergence
speed even in presence of large domain of attractions.

The rest of this paper is organised as follows. Section 2 provides the
formulation of the problem and reviews the basic ES scheme proposed by [8]. In
this section, we show that the averaging based on the Taylor expansion is not
a genuine representative of the ES scheme whose application can be extended to
non strictly convex cost functions via a different averaging analysis. Section
3 describes the ES scheme improved with the high-pass filter and states the
result about the uniformity of $\gamma$ with respect to the cost amplitude.
The performance of the investigated ES schemes are tested through the
simulations detailed in Section 4. Finally, Section 5 closes this paper with
some concluding comments. For lack of space, all the proofs of Lemmas,
Propositions, and Theorems claimed in this paper are reported in [25].

##  2 Problem Formulation

The ES problem consists in the optimisation of an unknown cost function
$h\,:\,\mathbb{R}\to\mathbb{R}$ satisfying the following two assumptions.

######  Assumption 1.

The function $h$ is smooth and there exists a $x^{\star}\in\mathbb{R}$ such
that

$h(x)-h(x^{\star})>0\quad\forall\,x\in\mathbb{R}\,:\,x\neq x^{\star}.$

######  Assumption 2.

There exist a locally Lipschitz and strictly quasi-convex function
$m\,:\,\mathbb{R}\to\mathbb{R}$, a class-${\cal K}_{\infty}$ function
$\alpha(\cdot)$, and a $A\geq 0$ such that for all
$x_{1},x_{2}\in\mathbb{R}\,:\,(x_{1}-x^{\star})(x_{2}-x^{\star})\geq 0$

$|m(x_{2})-m(x_{1})|\geq\alpha(|x_{2}-x_{1}|).$

######  Remark 1.

As for Assumption 2, it asks a (not necessarily strict) monotone behaviour of
$h$, with the latter that could have isolated saddle points. In other words,
$h(\cdot)$ belongs to the class of the so called strictly quasi-convex
functions [26]. Clearly, Assumption 2 is weaker than the common assumption
$(\partial h(x)/\partial x)x>0$ for any $x\neq x^{\star}$ typically present in
literature, (see, among the others, [[8], Assumptions 3 and 4]) ruling out the
existence of local saddle points.

The problem of semiglobal extremum seeking can be formulated in the following
way. For any $\epsilon>0$ and $r_{0}>0$, design a system of the form

$\dot{x}=\varphi_{\epsilon,r_{0}}(x,h(x),t)\quad x(0)=x_{0},$

so that for all $x_{0}$ satisfying $|x_{0}-x^{\star}|\leq r_{0}$ the resulting
trajectories $x(t)$ are bounded and satisfy
$\lim_{t\to\infty}\sup|x(t)-x^{\star}|\leq\epsilon$.

Among the different ES schemes proposed in literature to solve the previous
problem, a common one is given by (see [8])

$\dot{x}=-\gamma\,y_{\delta}(x,t)\,u(t)\qquad x(0)=x_{0}$ (1)

in which $y_{\delta}(x,t):=h(x+\delta u(t))$, $u(t):=\sin(2\pi t)$ is the
dither signal and $\gamma,\delta>0$ are tunable parameters111 The general case
of a dither takes the form $\tau\mapsto\sin(\omega\tau)$ with $\omega>0$,
$\tau\in\mathbb{R}$ (as considered in [8]) can be always obtained by rescaling
the time as $t=\tau\,2\pi/\omega$.. The block diagram of this algorithm is
represented in Figure 2. In the next part we briefly comment the main
properties, as available in the literature, of this algorithm and strengthen
them.

Since the right hand side of (1) is 1-periodic, the average system linked to
(1) is given by ([27], §10.4)

$\dot{{x}}_{a}=-\gamma\int_{0}^{1}y_{\delta}({x}_{a},t)u(t)dt\,.$ (2)

In the remaining part of the section we present a different route for the
analysis of (2) based on a Fourier serier expansion rather than on the Taylor
one used in the available literature. This analysis allows one to claim
stronger results on (1).

$x$$h(x)$$h(s)$$t$$y_{\delta}(s,t)$$s+\delta u(t)$$t$$s$

Figure 1: Graphical representation of the periodic behaviour of
$y_{\delta}(s,t)$. Conceiving $s$ as a parameter, we see that the oscillation
induced by $u(t)$ is elaborated by the non linear map $h(\cdot)$ and results
in a periodic function of time (top-right). This latter can be described
through its Fourier coefficients $a_{k}(s)$ and $b_{k}(s)$.

The smoothness of $h(\cdot)$ in Assumption 1 is essentially asked to guarantee
the existence of the Fourier series of the function $y_{\delta}$222Milder
regularity properties guaranteeing the existence of the series could be
assumed.. Then, since $y_{\delta}(x_{a},t)$ and its time derivatives are
continuous and periodic, $y_{\delta}(x_{a},t)$ can be expressed in terms of
its Fourier series as

$\displaystyle y_{\delta}(x_{a},t)=\dfrac{a_{0,\delta}(x_{a})}{2}+$ (3)
$\displaystyle\sum_{k=1}^{\infty}a_{k,\delta}(x_{a})\cos\left(k2\pi
t\right)+b_{k,\delta}(x_{a})\sin\left(k2\pi t\right)$

where

$\displaystyle a_{k,\delta}(x_{a})$
$\displaystyle:=2\int_{0}^{1}y_{\delta}(x_{a},t)\cos(k2\pi t)\,dt$ (4)
$\displaystyle b_{k,\delta}(x_{a})$
$\displaystyle:=2\int_{0}^{1}y_{\delta}(x_{a},t)\sin(k2\pi t)\,dt\,.$

Embedding (3) in (2) it is immediately seen that the average system linked to
(1) reads as

$\dot{x}_{a}=-{\gamma\over 2}b_{1,\delta}(x_{a})\,.$ (5)

For this system the following result holds. In the result we refer to the
class-${\cal K}_{\infty}$ function $\underline{\delta}^{\star}(\cdot)$ defined
as

$\underline{\delta}^{\star}(s):=2\int_{0}^{1/2}\alpha(s\sin(2\pi t))\,dt$ (6)

######  Lemma 1.

Let $h(\cdot)$ be such that Assumptions 1-2 are satisfied. Then:

  * a)

for all positive $\delta$ there exists a compact set ${{\cal
A}_{\delta}}\subseteq[x^{\star}-\delta,x^{\star}+\delta]$ that is globally
asymptotically and locally exponentially stable for (5).

  * b)

There exists a $\bar{\delta}^{\star}>0$ such that, for all positive $\delta$
such that $\delta\leq\bar{\delta}^{\star}$, there exists an equilibrium point
$x_{a\delta}^{\star}\in\mbox{\em int }{{\cal A}_{\delta}}$ that is locally
exponentially stable for system (5). If, in addition, the function
$\bar{h}(s):=h(x^{\star}+s)-h(x^{\star})$ is even, then
$x_{a\delta}^{\star}=x^{\star}$.

Item a) of the previous lemma states that the trajectories of the average
system reach a compact set ${\cal A}_{\delta}$ that is contained in a $\delta$
neighbourhood of $x^{\star}$ for all possible $\delta$. This, in particular,
implies that there exists a class $\cal KL$ function
$\beta(\cdot,\cdot):\mathbb{R}\times\mathbb{R}_{+}\to\mathbb{R}_{+}$ such that

$|x_{a}(t)|_{{\cal A}_{\delta}}\leq\beta(|x_{a}(0)|_{{\cal A}_{\delta}},t)\,$

where $|\cdot|_{\mathcal{A}_{\delta}}$ denote the distance to the set
$\mathcal{A}_{\delta}$. Item $b)$ claims that the set ${\cal A}_{\delta}$
collapses to an equilibrium point if $\delta$ is also taken sufficiently
small. Moreover, the last point of item b) shows that $x^{\star}$ represents
the equilibrium point only for cost functions that are locally symmetric
around the optimum.

Standard averaging results can be then used to show that the same property is
preserved also for the trajectories of the original system (1) for
sufficiently small $\gamma$ but in a semiglobal and practical way. This is
detailed in the next Proposition 1 where we refer to the class ${\cal
K}_{\infty}$ function $\chi(s)$ defined as

$\chi(s):=\beta^{-1}(s,0).$

In the following analysis we denote by $L_{r}>0$ and $M_{r}>0$, respectively
the local Lipschitz constant and the upper bound of the function $h(\cdot)$ on
a closed interval of length $r$. In particular, regularity of $h$ implies that

  * •

for all $r>0$ there exists $L_{r}>0$ such that for all
$x_{1},x_{2}\in[x^{\star}-r,\,x^{\star}+r]$

$|h(x_{1})-h(x_{2})|\leq L_{r}|x_{1}-x_{2}|\,.$ (7)

  * •

for all $r>0$, there exists $M_{r}>0$ such that for all
$x\in[x^{\star}-r,\,x^{\star}+r]$

$|h(x)|\leq M_{r}\,.$ (8)

######  Proposition 1.

Let $h(\cdot)$ be such that Assumptions 1-2 hold and let $r,\delta,d$ be
arbitrary positive numbers such that $r-d-2\delta>0$. Let
$r_{0}:=\chi(r-d-2\delta)$. There exist a $\bar{t}(r_{0},d)$ and a
$\gamma^{\star}(M_{r},L_{r},\delta,d)>0$ such that for any
$\gamma\in(0,\,\gamma^{\star})$, any
$x_{0}\in\mathbb{R}\,:\,|x_{0}-x^{\star}|\leq r_{0}$, the trajectories of (1)
are bounded and

$|x(t)|_{{\cal A}_{\delta}}\leq d\qquad\forall\,t\geq{\bar{t}\over\gamma}\,.$

An immediate consequence of the previous result is the next corollary showing
that under Assumptions 1-2 system (1) solves the problem of semiglobal
extremum seeking formulated before is fulfilled with a sufficiently small $A$.

######  Corollary 1.

Let $h(\cdot)$ be such that Assumptions 1-2 are fulfilled and let $r_{0}$ and
$\epsilon$ be arbitrary positive numbers. Then, there exist
$\bar{t}(r_{0},\epsilon)>0$, $\bar{\delta}^{\star}(\epsilon)>0$ and
$\gamma^{\star}(r_{0},\bar{\delta}^{\star},\epsilon)>0$ such that for any
$\delta\in(0,\bar{\delta}^{\star})$, any $\gamma\in(0,\,\gamma^{\star})$ and
any $x_{0}\in\mathbb{R}\,:\,|x_{0}-x^{\star}|\leq r_{0}$, the trajectories of
(1) are bounded and

$|x(t)-x^{\star}|\leq\epsilon\qquad\forall\,t\geq{\bar{t}\over\gamma}\,.$

By going through the proof of Proposition 1, it is immediately seen that
$\gamma^{\star}$ is inversely proportional to $M_{r}$. As a consequence, the
higher the cost function is within the set where $x(t)$ ranges, the lower the
value of $\gamma$ and, in turn, the slower the convergence rate of $x$ to the
neighbourhood of the optimum. Section 3 presents an improvement of (1)
overtaking this limitation.

###  2.1 Comments on Taylor expansion-based averaging analyses

The analysis of (2) is typically approached, see [8], by using a Taylor
expansion of $h(\cdot)$ to obtain a system of the form

$\dot{{x}}_{a}=\,-\gamma c_{1}\delta\left.\dfrac{\partial h}{\partial
x}\right|_{{x}_{a}}-\gamma\sum_{k=2}^{\infty}c_{k}\,{\delta^{2k-1}}\,\left.\dfrac{\partial^{2k-1}h}{\partial
x^{2k-1}}\right|_{{x}_{a}}$ (9)

where $c_{k}$ are suitably defined positive coefficients. A key role in the
study of this system is played by the first order approximation

$\dot{{x}}_{a}=\,-\gamma c_{1}\delta\left.\dfrac{\partial h}{\partial
x}\right|_{{x}_{a}}\,.$ (10)

In fact, if our Assumption 2 is strengthen by asking that, for any $x\neq
x^{\star}$, it holds $(\partial h(x)/\partial x)(x-x^{\star})>0$ (respectively
$>\alpha(|x-x^{\star}|)$ with $\alpha(\cdot)$ a class-${\cal K}$ function, see
[8], Assumptions 3 and 4), then the Lyapunov arguments of [[8], eq. (45)]
demonstrate $x^{\star}$ to be a stable (respectively globally asymptotically
stable) equilibrium point for (10). These stability properties are transferred
to (2) for sufficiently small $\delta$. Then, averaging techniques [27] can be
used to prove that, for sufficiently small $\gamma$, the trajectories of (1)
and (2) remain arbitrarily close. Namely, the semiglobal extremum seeking
problem is solved.

The fact that the asymptotic properties of the average system (2) are just
ensured by the first order term (10), and thus by the gradient of $h$, implies
that isolated local saddle points, of $h$, cannot be handled by that proof
technique. This justifies why [[8], Assumption 3], which is stronger than our
Assumption 2, is needed. Furthermore, we observe that the previous analysis
requires that the dither amplitude is kept sufficiently small for the higher
order terms of the average dynamics be negligible.

##  3 The High-Pass Filer (HPF)-ES Algorithm

$h(s)$$+$Eq. (11b)$-\bar{y}$$\dot{{x}}=-\gamma
v$$\times$$+$$v$${x}$$s$$u$$\delta$ Figure 2: The HPF-ES scheme consists of an
extension of the classic ES algorithm with the addition of the dashed block.

In [8] the authors proposed the next modification of (1)

$\displaystyle\dot{x}=$
$\displaystyle\,-\gamma\left(y_{\delta}({x},t)-\bar{y}\right)\,u(t)$
$\displaystyle{x}(0)={x}_{0}$ (11a) $\displaystyle\dot{\bar{y}}=$
$\displaystyle\,\gamma\left(y_{\delta}({x},t)-\bar{y}\right)$
$\displaystyle\bar{y}(0)=\bar{y}_{0}$ (11b)

where $u(t)$ is the dither signal defined before, and
$(x,\bar{y})\in\mathbb{R}\times\mathbb{R}$. A block representation of (11) is
depicted in Figure 2.

The intuition behind the previous scheme is to interpret
$y_{\delta}({x},t)-\bar{y}$ as the output of a high pass filter of
$y_{\delta}$. Moreover, the difference $y_{\delta}({x},t)-\bar{y}$ represents
a proxy of the local mean variation of $h(x)$, directly proportional to the
mean local Lipschitz constant. In the following we show how this feature
guarantees that the upper bound for the value of $\gamma$ is not dependent on
$M_{r}$ (eq. (8)) but rather only on $L_{r}$ ((7)). As for $\gamma$, similarly
to (1), it must be small to let $x$ and $\bar{y}$ be sufficiently slow to
preserve the correlation between the oscillations of
$y_{\delta}({x},t)-\bar{y}$ and those of $u(t)$. The average system of (11) is
defined as

$\displaystyle\dot{{x}}_{a}$
$\displaystyle=-\gamma\int_{0}^{1}\left(y_{\delta}({x_{a}},\tau)-\bar{y}_{a}\right)\,u(\tau)\,d\tau$
(12a) $\displaystyle\dot{\bar{y}}_{a}$
$\displaystyle=-\gamma\,\bar{y}_{a}+\gamma\int_{0}^{1}y_{\delta}({x}_{a},\tau)\,d\tau.$
(12b)

By expanding $y_{\delta}(x_{a},t)$ in terms of Fourier series as in the
previous section, it turns out that

$\displaystyle\displaystyle\int_{0}^{1}y_{\delta}(x_{a},\tau)\,d\tau$
$\displaystyle=\dfrac{a_{0,\delta}(x_{a})}{2}$ (13a)
$\displaystyle\displaystyle\int_{0}^{1}\left(y_{\delta}({x_{a}},\tau)-\bar{y}_{a}\right)u(\tau)\,d\tau$
$\displaystyle=\dfrac{b_{1,\delta}(x_{a})}{2}$ (13b) where in the latter we
exploited $\int_{0}^{1}\bar{y}_{a}u(\tau)d\tau=0$ and (4).

Hence, the average system reads as

$\displaystyle\dot{{x}}_{a}$
$\displaystyle=-\gamma\dfrac{b_{1,\delta}(x_{a})}{2}$ (14a)
$\displaystyle\dot{\bar{y}}_{a}$
$\displaystyle=-\gamma\,\bar{y}_{a}+\gamma\,\dfrac{a_{0,\delta}({x}_{a})}{2}$
(14b)

which is a cascade where the first subsystem coincides with (5) and the second
subsystem is linear and asymptotically stable. From this, the next result
follows from Lemma 1.

######  Lemma 2.

Let $h(\cdot)$ be such that Assumptions 1-2 are satisfied. Then:

  * a)

for any positive $\delta>0$ there exist a compact set ${\cal
A}_{\delta}\subseteq[x^{\star}-\delta,x^{\star}+\delta]$ and a continuous
function $\tau:\mathbb{R}\to\mathbb{R}$ such that the set

$\mbox{\rm graph}\left.\tau\right|_{{\cal
A}_{\delta}}=\\{(x_{a},y_{a})\in{\cal
A}_{\delta}\times\mathbb{R}\;:\;y_{a}=\tau(x_{a})\\}$

is globally asymptotically and locally exponentially stable for (14).

  * b)

There exists $\bar{\delta}^{\star}>0$ such that, for any
$\delta\in(0,\bar{\delta}^{\star})$, there exists an equilibrium point
$(x_{a\delta}^{\star},\bar{y}_{a\delta}^{\star})\in\mbox{\rm
graph}\left.\tau\right|_{{\cal A}_{\delta}}$ that is globally asymptotically
and locally exponentially stable for system (5). If, in addition, the function
$\bar{h}(s):=h(x^{\star}+s)-h(x^{\star})$ is even, then
$x_{a\delta}^{\star}=x^{\star}$.

From this, the following result mimics the one of Proposition 1 with the
remarkable difference that the upper bound $\gamma^{\star}$ of $\gamma$ is
uniform with respect to $M_{r}$.

######  Theorem 1.

Let $h(\cdot)$ be such that Assumptions 1-2 are fulfilled and let $r,\delta,d$
be arbitrary positive numbers such that $r-d-2\delta>0$. Let
$r_{0}:=\chi(r-d-2\delta)$. There exist a $\bar{t}(r_{0},d)$ and a
$\gamma^{\star}(L_{r},\delta,d)>0$ such that for any
$\gamma\in(0,\,\gamma^{\star})$ and any $(x_{0},\,\bar{y}_{0})$ fulfilling
$|x_{0}-x^{\star}|\leq r_{0}$ and
$|\bar{y}_{0}-a_{0,\delta}(x_{0})/2|\leq\gamma^{\star}$, then the trajectories
of (11) are bounded and

$\|(x(t),\bar{y}(t))\|_{\mbox{\rm graph}\left.\tau\right|_{{\cal
A}_{\delta}}}\leq d\qquad\forall\,t\geq{\bar{t}\over\gamma}\,.$

##  4 Numerical Results

This section presents numerical results obtained adopting

the following cost function $h(\cdot)\,:\,\mathbb{R}\to\mathbb{R}$

$h(x)=h_{0}+\left\\{\begin{array}[]{cc}(x-\pi)^{2}-1&x<\pi\\\
\cos(x-\pi)-2&x\in[\pi,\,2\pi)\\\ (x-2\pi)^{2}-3&x\geq 2\pi\end{array}\right.$
(15)

where $h_{0}\in\mathbb{R}$. This function, which verifies Assumptions 1 and 2,
is depicted in Figure 3, for $h_{0}=10$, on $[-\pi,\,3\pi]$.

Figure 3: Cost function (15) for $h_{0}=0$. It presents multiple points with
null first derivative and it is also non-symmetric around the minimiser. The
saddle point is located at $x=\pi$ and the minimiser is at $x^{\star}=2\pi$.
Figure 4: The presence of saddle points stacks the Taylor-based averaging
(yellow) of the classic ES (blue). Vice versa, the Fourier-based averaging
(red) better represents the actual behaviour of the classic ES (blue). These
results are obtained for $h_{0}=0$ and $\gamma=0.1$. Figure 5: The average
based on the Taylor expansion (yellow) of the classic ES wrongly assess
$x^{\star}$ as an equilibrium point. Indeed, due to the asymmetry of
$h(\cdot)$ around $x^{\star}$, the classic ES (blue) converges to an
equilibrium point different from $x^{\star}$. Moreover, the average based on
the Fourier series (red) tracks the classic ES more accurately than the
average based on the Taylor expansion. These results are obtained for
$h_{0}=0$ and $\gamma=0.1$. Figure 6: The average of the classic ES, based on
the Fourier series (red), converges to an equilibrium point within the set
$[x^{\star}-\delta,\,x^{\star}+\delta]$ accordingly to what foreseen in Lemma
1. Thus, the smaller is $\delta$ the closer the equilibrium point is to
$x^{\star}$. These results are obtained for $h_{0}=0$ and $\gamma=0.1$.

Figures 4 and 5 show that the asymmetry of $h(\cdot)$ in the neighbourhood of
$x^{\star}$ implies that the equilibrium point of (2) does not correspond to
$x^{\star}$ as wrongly assessed by the averaging through Taylor expansion, see
yellow lines in Figures 5 and 6. Vice versa, Figures 5 and 6 show that the
trajectory of (5) (red) tracks, better than those of (10) (yellow), the
trajectory of (1) (blue) and, as claimed in Lemma 1, converges into a
$2\delta$-wide set centred at $x^{\star}$.

The second test is performed to show that the classic ES suffers of large
values of $|h(\cdot)|$. Indeed, as depicted in Figure 7 (blue lines in
subplots from (a) to (d)), while keeping $\gamma$ fixed, larger $M_{r}$ lead
to more oscillatory behaviours. To mathematically support this result we
observe that

$h(x+\delta u)=h(x)+R(x,\delta u)$ (16)

where $R(\cdot,\cdot)$ represents the remainder of the Taylor expansion around
$x$ of $h(\cdot)$. Exploit the definition of the Lipschitz constant of
$h(\cdot)$ and $|u(t)|_{\infty}\leq 1$ to bound the remainder from above as

$|h(x+\delta u)-h(x)|=|R(x,\delta u)|\leq L_{r}\delta.$ (17)

Substitute (16) into (1)

$\dot{x}=-\gamma h(x+\delta u)u=-\gamma(h(x)+R(x,\delta u))u$ (18)

and investigate the following support system

$\dot{x}_{1}=-\gamma h(x_{1})u\qquad x(0)=x_{10}$ (19)

conceivable as approximation of (18) for $|h(x)|\gg L_{r}\delta$. Let
$H(x):=\int h(x)^{-1}dx$ and solve (19) by parts as

$x_{1}(t)=H^{-1}(H(x_{10})-\gamma(\cos(t)-1))$ (20)

which is a pure oscillation whose amplitude is proportional to $|H(x_{10})|$,
evident in subplot (d) of Figure 7.

(a)

(b)

(c)

(d)

Figure 7: The behaviour of the classic ES changes for increasing $M_{r}$. In
these simulations the value of $h_{0}$ is increased from $-10$ to $10^{3}$
(subplots from (a) to (d)) demonstrating that, keeping fixed $\gamma$, the
classic ES (blue) becomes oscillatory for large values of the cost function.
These simulations are performed with $\gamma=\delta=0.1$.

The third group of simulations, reported in Figure 8, shows the performance of
the HPF-ES scheme whose main feature consists of better convergence
performance in case of large cost functions.

Figure 8: The behaviour of the HPF-ES does not change for increasing $M_{r}$.
In these simulations the value of $h_{0}$ is increased from $-10$ to $10^{3}$
demonstrating that, keeping fixed $\gamma$, the HPF-ES keeps the same speed of
convergence. These simulations are performed with $\gamma=\delta=0.1$. (left)
Behaviour of $x$ and $x_{a}$. (right) Behaviour of $\bar{y}$ and
$\bar{y}_{a}$.

Through the proposed simulations we confirm three results: a) the classic ES
is able to deal with cost functions with local saddle points where, at the
opposite, the average of the classic ES obtained through the Taylor expansion
gets stacked (Figure 4); b) the classic ES and its average via Taylor
expansion do not converge to the same equilibrium point (Figure 5); c) in case
of large cost functions the classic ES behaves as an oscillator whose
amplitude is proportional to the value of the cost function (Figure 7(d)).
Vice versa, the adoption of the high-pass filter makes the basic ES
convergence rate uniform with respect to the amplitude of the cost function.

##  5 Conclusions

This paper deals with two well-known extremum seeking schemes to show that
they work under less restrictive assumptions. Relying on averaging and
Fourier-series arguments, it is demonstrated that these schemes are able to
deal with strictly quasi-convex cost functions making the global minimiser
(assumed to be unique) semi-global practically stable. Moreover, it is shown
that the presence of a high-pass filter elaborating the cost function makes
the tuning of the parameters independent on the cost magnitude.

## References

  * [1] D. Dochain, M. Perrier, and M. Guay, “Extremum seeking control and its application to process and reaction systems: A survey,” Mathematics and Computers in Simulation, vol. 82, no. 3, pp. 369 – 380, 2011.  6th Vienna International Conference on Mathematical Modelling.
  * [2] S. F. Toloue and M. Moallem, “Multivariable sliding-mode extremum seeking control with application to mppt of an alternator-based energy conversion system,” IEEE Transactions on Industrial Electronics, vol. 64, no. 8, pp. 6383–6391, 2017.
  * [3] H. Malek and Y. Chen, “Fractional order extremum seeking control: Performance and stability analysis,” IEEE/ASME Transactions on Mechatronics, vol. 21, no. 3, pp. 1620–1628, 2016.
  * [4] C. Zhang and R. Ordonez, “Numerical optimization-based extremum seeking control with application to abs design,” IEEE Transactions on Automatic Control, vol. 52, no. 3, pp. 454–467, 2007.
  * [5] Hsin-Hsiung Wang, S. Yeung, and M. Krstic, “Experimental application of extremum seeking on an axial-flow compressor,” IEEE Transactions on Control Systems Technology, vol. 8, no. 2, pp. 300–309, 2000.
  * [6] D. J. Burns, C. R. Laughman, and M. Guay, “Proportional–integral extremum seeking for vapor compression systems,” IEEE Transactions on Control Systems Technology, vol. 28, no. 2, pp. 403–412, 2020.
  * [7] W. H. Moase, C. Manzie, and M. J. Brear, “Newton-like extremum-seeking for the control of thermoacoustic instability,” IEEE Transactions on Automatic Control, vol. 55, no. 9, pp. 2094–2105, 2010.
  * [8] Y. Tan, D. Nešić, and I. Mareels, “On non-local stability properties of extremum seeking control,” Automatica, vol. 42, no. 6, pp. 889–903, 2006.
  * [9] M. Krstić and H.-H. Wang, “Stability of extremum seeking feedback for general nonlinear dynamic systems,” Automatica, vol. 36, no. 4, pp. 595 – 601, 2000.
  * [10] M. Guay, “A perturbation-based proportional integral extremum-seeking control approach,” IEEE Transactions on Automatic Control, vol. 61, no. 11, pp. 3370–3381, 2016.
  * [11] M. Guay and D. Dochain, “A minmax extremum-seeking controller design technique,” IEEE Transactions on Automatic Control, vol. 59, no. 7, pp. 1874–1886, 2014.
  * [12] D. Nesic, A. Mohammadi, and C. Manzie, “A framework for extremum seeking control of systems with parameter uncertainties,” IEEE Transactions on Automatic Control, vol. 58, no. 2, pp. 435–448, 2013.
  * [13] L. Hazeleger, M. Haring, and N. van de Wouw, “Extremum-seeking control for optimization of time-varying steady-state responses of nonlinear systems,” Automatica, vol. 119, p. 109068, 2020.
  * [14] T. R. Oliveira, M. Krstić, and D. Tsubakino, “Extremum seeking for static maps with delays,” IEEE Transactions on Automatic Control, vol. 62, no. 4, pp. 1911–1926, 2017.
  * [15] C. Labar, E. Garone, M. Kinnaert, and C. Ebenbauer, “Newton-based extremum seeking: A second-order lie bracket approximation approach,” Automatica, vol. 105, pp. 356–367, 2019.
  * [16] T. R. Oliveira, J. Feiling, S. Koga, and M. Krstić, “Multivariable extremum seeking for pde dynamic systems,” IEEE Transactions on Automatic Control, vol. 65, no. 11, pp. 4949–4956, 2020.
  * [17] Y. Zhang, M. Rotea, and N. Gans, “Simplex guided extremum seeking control with convergence detection to improve global performance,” IEEE Transactions on Control Systems Technology, vol. 24, no. 4, pp. 1266–1278, 2016.
  * [18] F. V. J.A. Sanders, Averaging Methods in Nonlinear Dynamical Systems.  Applied Mathematical Sciences, Springer-Verlag, 1985.
  * [19] A. Teel, L. Moreau, and D. Nesic, “A unified framework for input-to-state stability in systems with two time scales,” IEEE Transactions on Automatic Control, vol. 48, no. 9, pp. 1526–1544, 2003.
  * [20] Y. Tan, D. Nešić, and I. Mareels, “On non-local stability properties of extremum seeking control,” IFAC Proceedings Volumes, vol. 38, no. 1, pp. 550 – 555, 2005.  16th IFAC World Congress.
  * [21] A. Teel and D. Nešić, “Averaging with disturbances and closeness of solutions,” Systems & Control Letters, vol. 40, no. 5, pp. 317 – 323, 2000.
  * [22] A. R. Teel, J. Peuteman, and D. Aeyels, “Semi-global practical asymptotic stability and averaging,” Systems & Control Letters, vol. 37, no. 5, pp. 329 – 334, 1999.
  * [23] A. Scheinker and M. Krstić, “Extremum seeking with bounded update rates,” Systems & Control Letters, vol. 63, pp. 25–31, 2014.
  * [24] Y. Zhu and E. Fridman, “Extremum seeking via a time-delay approach to averaging,” Automatica, vol. 135, p. 109965, 2022.
  * [25] N. Mimmo, L. Marconi, and G. Notarstefano, “Supplement to the paper: Uniform quasi-convex optimisation via exrtremum seeking,” arXiv, arXiv:submit/4239409, 2022.
  * [26] S. Karamardian, “Strictly quasi-convex (concave) functions and duality in mathematical programming,” Journal of Mathematical Analysis and Applications, vol. 20, no. 2, pp. 344–358, 1967.
  * [27] H. K. Khalil and J. W. Grizzle, Nonlinear systems, vol. 3.  Prentice hall Upper Saddle River, NJ, 2002.
  * [28] S. Z. Khong, C. Manzie, D. Nešić, and Y. Tan, “Multidimensional global extremum seeking via the direct method,” in 2012 2nd Australian Control Conference, pp. 277–282, 2012.
  * [29] J. K. Hale, L. T. Magalhães, and W. Oliva, Dynamics in infinite dimensions, vol. 47.  Springer Science & Business Media, 2006.
  * [30] L. Marconi, L. Praly, and A. Isidori, “Robust asymptotic stabilization of nonlinear systems with non-hyperbolic zero dynamics,” IEEE Transactions on Automatic Control, vol. 55, no. 4, pp. 907–921, 2010.
  * [31] L. Marconi and L. Praly, “Essential and redundant internal models in nonlinear output regulation,” in Analysis and Design of Nonlinear Control Systems, pp. 263–283, Springer, 2008.
  * [32] L. Marconi, L. Praly, and A. Isidori, “Output stabilization via nonlinear luenberger observers,” SIAM Journal on Control and Optimization, vol. 45, no. 6, pp. 2277–2298, 2007.