text:
\begin{document}

{\title{ Detecting nonlinear corrosion \\ by electrostatic
measurements\thanks{Work supported in part by MIUR, grant n.
2002013279.} }}
\author{G. Alessandrini\thanks{Dipartimento di  Matematica e Informatica, Universita' degli Studi di Trieste, Italy,
\tt{alessang@univ.trieste.it}} \ and \ E.
Sincich\thanks{S.I.S.S.A., Trieste, Italy, \tt { sincich@sissa.it}
}}

\date{}

\maketitle

\begin{abstract}We deal with an inverse problem arising in corrosion detection. We
prove a stability estimate for a nonlinear term on the
inaccessible portion of the boundary by electrostatic boundary
measurements on the accessible one.
\end{abstract}

{\small{\bf Keywords }:  \ inverse boundary problems, corrosion,
stability.}

{\small{\bf 2000 Mathematics Subject Classification }: 35R30, 35R25, 31B20 .}

\section{Introduction}
In this paper we deal with the inverse problem of determining the
nonlinear term $f(u)$ in the boundary value problem
\begin{equation}\label{P}
\left\{
\begin{array}
{lcl}
\Delta u=0\ ,& \mbox{in $\Omega$ ,}
\\
\dfrac{\partial u}{\partial \nu}=g\ ,   & \mbox{on $\G_2$ ,}
\\
\dfrac{\partial u}{\partial \nu}=f(u)\ , & \mbox{on $\G_1$ ,}
\\
u=0\ , & \mbox{on $\G_{D}$ ,}
\end{array}
\right.
\end{equation}
where $\G_1$ and $\G_2$
are two open, disjoint portions of  $\partial\O$
and $\G_{D}=\partial \O \setminus (\G_1\cup\G_2)$ and $f$
is a Lipschitz
 function such that $f(0)=0$.

\noindent Such a problem with a specific choice of the profile of $f$, has been introduced and discussed in recent years by M.Vogelius and others for the modeling of the electrochemical phenomenon of surface corrosion in metals \cite{Vog1},\cite{Vog2},\cite{Vog3}.
 In the above boundary value problem, $\O$ represents the metal specimen, $\G_1$ represents the corroded part of the boundary, which is not accessible to direct inspection, $\G_2$ is the portion of the boundary which is accessible to measurements and $\G_D$ is the remaining part of $\partial \O$ which is assumed to be grounded.

\noindent The inverse problem thus consists in the determination of $f$ when one pair of Cauchy data $u|_{\G_2},{\frac{\partial u}{\partial \nu}}|_{\G_2}$ is known for one non-trivial solution $u$ to \eqref{P}.

\noindent It must be noted that the models of nonlinearity which have been discussed in \cite{Vog2},\cite{Vog3} are of the form
$$f(u)=\lambda(\exp(\a u)-\exp(-(1-\a)u))$$
and it turns out that are such that in the direct problem the
existence and the uniqueness of the solution are not granted. For
this reason, we found it necessary to require some additional a
priori assumptions. First an energy bound of the measured
electrostatic potential $u$
\begin{eqnarray}\label{E}
\int_{\O}|\nabla u(x)|^2\le E^2\ .
\end{eqnarray}
Next, an a priori bound of the Lipschitz continuity of $f$, namely
\begin{eqnarray}\label{Lip}
|f(u)-f(v)|\le L|u-v|\ , \ \mbox{for every}  \ u, v \ \in
\mathbb{R}\ .
\end{eqnarray}
Moreover, we found it necessary to assume the knowledge of some
additional information on the measured current density $g$ on the
accessible part of the boundary $\G_2$. More precisely, we assume
a bound on the H\"older continuity of $g$, with
\begin{eqnarray}\label{G}
 \|g\|_{C^{0,\a}(\G_2)}\le G\ .
\end{eqnarray}
Also, we shall require a lower bound on the same current density
$g$. Namely, we shall prescribe that, for a given inner portion
$\G_{2,2r_0}$ of $\G_2$ (see \eqref{gi} below, for a precise definition)
and a given number $m > 0$, we have
\begin{equation}\label{lower b}
\|g\|_{L^{\infty}(\G_{2,2r_0})}\ge m>0 \ .
\end{equation}

We note that, indeed, a lower bound of this sort appears to be
necessary. In fact, knowing the Cauchy data on $\G_2$ for a single
potential $u$, one can expect to identify $f$ only on the range of
values taken by $u$ on the inaccessible boundary $\G_1$. Thus, as
a preliminary step of the treatment of this inverse problem, it is
necessary to evaluate the amplitude of the range of $u$ on $\G_1$,
see Theorem \ref{teo1}, in which we prove that the oscillation of
$u$ on $\G_1$ is bounded from below by
$\exp(-\big(\frac{m}{c}\big)^{-\gamma})$, where $c>0,\gamma>1$ are
constants depending on the a priori data only. Next, as the main
result of this paper, we prove a stability estimate,  Theorem
\ref{teo2}, that is we show that if $u_1, u_2$ are two potentials
corresponding to nonlinearities $f_1,\ f_2,$ whose Cauchy data are
close
\begin{eqnarray*}
&&\|{\psi}_1-{\psi}_2\|_{L^2(\G_2)}\le \varepsilon\ ,\\
&&\|g_1-g_2\|_{L^2(\G_2)}\le \varepsilon\ ,
\end{eqnarray*}
where $\psi_i=u_1\big|_{\G_2}$ and $\displaystyle{\frac {\partial
u_i}{\partial \nu}}\bigg|_{\G_2}=g_i$ with $i=1,2$ and if $g=g_1$
satisfies the a priori bounds (\ref{G}), (\ref{lower b}), then the
ranges of $u_1, u_2$ on $\G_1$ agree on a sufficiently large
interval $V$ such that
$$\mbox{length of}\  V\  \sim \ \ \exp\bigg[- \bigg(\frac{m}{c}\bigg)^{-\gamma}\bigg] \ .$$
On such an interval the nonlinearities $f_1, f_2$ agree up to an
error $\omega(\varepsilon)$ such that
$$\omega(\varepsilon) \sim \left|\log\bigg(\frac{1}{\varepsilon}\bigg)\right|^{-\theta}$$
where $0<\theta<1$.

The study of determination of nonlinear terms in elliptic inverse
boundary value problems has produced various different results.
Just to mention a few, Cannon \cite{Ca}, Beretta and Vogelius
\cite{BV}, Isakov and Sylvester \cite{IS}, Sun \cite{Sun}. Among
the specific peculiarities of the present study, in comparison
with the mentioned previous results, we emphasize that the
nonlinearity is part of a boundary condition and that existence
and uniqueness for the corresponding direct problem are not
assured, see, in this respect, Kavian and Vogelius \cite{Vog2} and
Medville and Vogelius \cite{MV}. For the determination of a
coefficient in a \emph{linear} boundary condition, we refer for
instance to Fasino and Inglese \cite{FI}, Alessandrini, Del Piero
and Rondi \cite{adr} and Chaabane, Fellah, Jaoua and Leblond
\cite{CJ}.

\subsection{Plan of the paper}
In Section 2 we formulate our main hypotheses and state the main
results, Theorem \ref{teo1} and Theorem \ref{teo2}, whose proofs
are deferred to Section 3 and 4 respectively. In Section 3 we
prove some preliminary lemmas, in which we observe that, taking
into account the a priori assumptions, the solution is H\"older
continuous with its first order derivatives in a suitable
neighborhood of the inaccessible part $\G_1$, (see Lemma
\ref{localbound}, Theorem \ref{teoreg}). In order to prove Theorem
\ref{teo1} we need also a stability estimate near the boundary for
a Cauchy problem (see Proposition \ref{tritten}). The proof of
this Proposition is based on a method developed by Payne and
Trytten \cite{payne}, \cite{tr}. We conclude the Section with the
proof of Theorem \ref{teo1}. In Section 4 we adapt the above
mentioned stability estimate for the Cauchy problem (see Theorem
\ref{stabilita'}), when both the Dirichlet and Neumann data are
affected by errors in the $L^2$ norm, and at the same time, we
obtain a stability result up to the boundary. In Proposition
\ref{inv} we show that the solution is locally invertible along a
suitable curve on $\G_1$ and we evaluate its length. Finally we
prove Theorem \ref{teo2}.

\section{Main assumptions and results}
\subsection{Preliminaries}
{\it A priori information on the domain}
\\
\\

We shall assume throughout that $\O$ is a bounded, simply
connected domain in $\Rn$, $n\ge2$ such that $\mbox{diam} (\O)\le
D$ and with Lipschitz boundary $\partial \O$ with constants $r_0,
M$. More precisely, for every $x_0 \in \partial \O$, exists a
rigid transformation of coordinates under which,
\begin{eqnarray}\label{gamma}
\O \cap B_{r_0}(x_0)=\{(x',x_n): x_n>\gamma(x')\}
\end{eqnarray}
where $x \in \Rn,\ x=(x',x_n)$, with $x'  \in \Rnu , \ x_n \in \mathbb{R}$ and
$$\gamma :B_{r_0}(x_0)\subset \Rnu \rightarrow \mathbb{R}\ ,$$
satisfying $\gamma(0)=0$ and $$\|\gamma\|_{C^{0,1}( B_{r_0}(x_0))}\le Mr_0\ ,$$ where we denote by
$$\|\gamma\|_{C^{0,1}( B_{r_0}(x_0))}=\|\gamma\|_{L^{\infty}( B_{r_0}(x_0))} +\ r_0\!\!\!\!\!\!\!\!\!\sup_{\substack {x,y  \in B_{r_0}(z_0)\\x\ne y }}
\frac{|\gamma (x)-\gamma (y)|}{|x-y|}\ . $$
Moreover, we assume that the portions of the boundary $\Gamma_i$ are contained respectively into surfaces $S_i$, $i=1,2$ which are $C^{1,\a}$ smooth with constants $r_0, M$.

More precisely,
 for a given $\a,\
 0<\a<1$ and for any $z_0 \in S_i,\  i=1,2 $, we have that up to a rigid change of coordinates,
\begin{eqnarray}\label{si}
 S_i \cap B_{r_0}(z_0)=\{(x',x_n): x_n=\varphi_i(x')\}\ ,
\end{eqnarray}
where
$$\varphi_i:B_{r_0}(z_0)\subset \Rnu \rightarrow \mathbb{R} $$
 are  $C^{1,\a}$ functions satisfying $\varphi_i(0)=|\nabla \varphi_i(0)|=0$ and $$\|\varphi_i \|_{C^{1,\a}(B_{r_0}(z_0))}\le Mr_0\ ,$$
where we denote
\begin{eqnarray}\nonumber
\|\varphi\|_{C^{1, \a}( B_{r_0}(z_0))}&=&\|\varphi\|_{L^{\infty}( B_{r_0}(z_0))}+r_0\|\nabla \varphi\|_{L^{\infty}( B_{r_0}(z_0))}+\\
& & +\ {r_0}^{1+\a}\!\!\!\!\!\!\!\!\!\sup_{\substack {x,y  \in B_{r_0}(z_0)\\x\ne y }}\frac{|\nabla \varphi (x)-\nabla \varphi (y)|}{|x-y|^{\a}}\ .
\end{eqnarray}

In particular it follows that if
$$x_0 \in \Gamma_i\ \ \mbox{and} \ \ \mbox{dist}(x_0,\G_D)>r_0\ ,$$
then
$$\O \cap B_{r_0}(z_0)=\{(x',x_n) \in  B_{r_0}(z_0): x_n>\varphi_i(x')\}\ ,$$
where $\varphi_i$ is the Lipschitz function whose graph locally represents $\partial \O$. Moreover, since $\O \cap B_{r_0}(z_0) \cap \G_D =\emptyset$, $\varphi_i$ must also be the $C^{1,\a}$ function whose graph locally represents $S_i$.
We also suppose that the boundary of $\G_i$, within $S_i$, is of $C^{1,\a}$ class with constants $r_0, M$, namely, for any $z_0\in \partial \G_i$, there exists a rigid transformation of coordinates under which
$$\partial \G_i \cap B_{r_0}(z_0)=\{(x',x_n)\in B_{r_0}(z_0):\ \ x_n=\varphi_i(x'),\  x_{n-1}=\psi_i(x'')\}$$
where $x'=(x'',x_{n-1})$, with $x'' \in \mathbb{R}^{n-2}$, $x_{n-1} \in \mathbb{R}$ and
$$\psi_i : B_{r_0}(z_0)\subset  \mathbb{R}^{n-2} \longrightarrow \mathbb{R} $$
satisfying $\psi_i(0)=|\nabla \psi_i(0)|=0$ and $\|\psi_i\|_{C^{1,\a} (B_{r_0}(z_0))}\le M$.

We introduce some notation that we shall use in the sequel
\begin{eqnarray}\label{ui}
U^i_{\rho}=\{x \in \bar{\O}:\mbox{dist}(x,\partial\O\setminus \G_i)>\rho \}\ \ ,
\end{eqnarray}
\begin{eqnarray}\label{gi}
\Gamma_{i,\rho}=U^1_{\rho}\cap \G_i\ ,
\end{eqnarray}
with i=1,2  and
\begin{eqnarray}\label{or}
\O_{\rho}=\{x\in \O: \mbox{dist}(x, \partial \O)>\rho \}\ .
\end{eqnarray}

 {\it A priori information on the boundary data}
\\
\\
The current flux $g$
is a prescribed function such that
\begin{eqnarray*}
 \|g\|_{C^{0,\a}(\G_2)}\le G\ .
\end{eqnarray*}

 {\it A priori information on the nonlinear term}
\\
\\
We assume that the function $f$ belongs to $C^{0,1}(\mathbb{R},\mathbb{R})$ and, in particular,
\begin{eqnarray}\label{f}
  f(0)=0\ \  \mbox{and}\ \  | f(u)-f(v)|\le L|u-v|\ \ \
 \mbox{for every } u, v\in \mathbb{R}\ .
\end{eqnarray}

We recall that a weak solution of problem \eqref{P} is a  function $ u \in \i$, such that $u|_{\G_D}=0$ and which satisfies
\begin{eqnarray}\label{soluzione}
\io \n u\cdot \n \rho  =\id g  \rho + \iu f(u) \rho \quad \mbox{for all } \rho \in \i, \, \,\rho _{|_{_{_{\G_D}}}}=0
\end{eqnarray}

where  $u_{|_{_{_{\G_D}}}}  $ denote the trace on $\G_D$ .
\\
From now on we shall refer to the \emph{a priori data} as to the set of quantities
$r_0, M, \a, L, G, E, D, m$.\\
In the sequel we shall denote with $\eta(t)$ and $\omega(t)$, two positive increasing functions defined on $(0, +\infty)$, that satisfy
\begin{eqnarray}\label{eta}
&&\eta(t)\ge\exp{\bigg[-\bigg(\frac{t}{c}\bigg)^{-\gamma}\bigg]}, \ \ \ \mbox{for every}\ \ 0<t\le G\ \ ,
\end{eqnarray}
\begin{eqnarray}\label{omega}
&&\omega(t)\le C \left|\log (t) \right|^{-\vartheta},\ \ \ \mbox{for every}\ \ 0<t<1\ \ ,
\end{eqnarray}
where $c>0$, $C>0$, $\gamma>1$, $0<\theta<1$ are constants depending on the \emph{a priori data} only.
\subsection{Main theorems}
The statements of the main results are the following.
\begin{theorem}[Lower bound for the oscillation]\label{teo1}
Let $\O, g $ satisfying the a priori assumptions. Let $u$ be a
weak solution of \eqref{P} satisfying the a priori bound \eqref{E}
then $$\mathop{\rm osc}\limits_{\G_1} u\ge
\eta(\|g\|_{L^{\infty}(\G_{2,2r_0})})$$ where $\eta$ satisfies
\eqref{eta}.
\end{theorem}
\begin{theorem}[Stability for the nonlinear term $f$]\label{teo2}
Let $u_i\ \in \i$, $i=1,2$ be weak solutions of the problem $\eqref{P}$, with $f=f_i$ and $g=g_i$
respectively and such that \eqref{E} holds for each $u_i$.
Let us also assume that, for some positive number $m$, the following holds
\begin{eqnarray}\label{m}
\|g_1\|_{L^{\infty}(\G_{2,2r_0})}\ge m>0\ .
\end{eqnarray}
Moreover, let $\psi_i=u_i\big|_{\G_2}, \ i=1,2$.
There exists $\eps_{0}>0$ only depending on the \emph{a priori data} and on $m$ such that, if, for some $\eps,\ 0<\eps<\eps_0$, we have
\begin{eqnarray*}
&& \|\psi_1 -\psi_2  \|_{L^{2}(\G_2)}\le \eps\ , \\
&& \|g_1 -g_2  \|_{L^{2}(\G_2)}\le \eps\ ,
\end{eqnarray*}
 then
$$\|f_1(u)-f_2(u) \|_{L^{\infty}(V)}\le \omega(\eps)\ , $$
where
 $$V=(\a,\beta)\subseteq [-CE,CE]\ ,$$
is  such that
 $$\beta -\a >\frac{\eta(m)}{2}$$
and $\eta, \omega$ satisfy \eqref{eta}, \eqref{omega} respectively.
\end{theorem}

\section{The lower bound for the oscillation}

\begin{proposition}[Stability for a Cauchy problem]\label{tritten}
Let $\O$ satisfies the a priori assumptions and let $v \in
H^1(\O)$ be a solution of the following Cauchy problem
\begin{equation}
\left\{
\begin{array}{lcl}
\Delta v=0\ ,& \mbox{in $\Omega$},
\\
v=\varphi\ ,  & \mbox{on $\G_1$},
\\
\displaystyle\frac{\partial v}{\partial \nu}=h, & \mbox{on $\G_1$},
\end{array}
\right.
\end{equation}
where $\varphi, h\in L^2(\G_1)$ and the boundary conditions are considered in the weak sense.\\
Then, for every $P_1 \in \G_{1,2r_0}$, $v$ satisfies the following estimate
\begin{eqnarray*}
\|v\|_{L^2(B_{\rho}(P_0)\cap U^1_{2r_0})}&\leq&\ C \left(\|\varphi\|_{L^2(\G_{1,\rho})} +\
\|h\|_{L^2(\G_{1,\rho})}\ +\|v\|_{\i}
\right)^{1-\delta}\cdot\\
&&\ \ \ \ \ \ \ \ \ \cdot\left(\|\varphi\|_{L^2(\G_{1,\rho})} +\
\|h\|_{L^2(\G_{1,\rho})}\right)^{\delta}
\end{eqnarray*}
where $ \rho \in \left(\frac{M}{4\sqrt{1+M^2}}r_0,\frac{3M}{4\sqrt{1+M^2}}r_0 \right)$, $P_0=P_1+\frac{M}{4\sqrt{1+M^2}}r_0\cdot \nu$, $\nu$ is the outer unit normal to $\O$ at $P_1$ and
 $C>0$, $0<\delta<1$ are constants depending on
 $ \rho,r_0, n, M $ only.
\end{proposition}

\begin{proof}By a well-known estimate of stability for a Cauchy problem, we have that for any $P_1 \in \G_{1,r_0}$
\begin{eqnarray}\label{grad}
\|v\|_{L^2(B_{\rho}(P_0)\cap U^1_{r_0})}&\leq&C_1 \left(\|\varphi\|_{L^2(\G_{1,\rho})} +\
\|\nabla v\|_{L^2(\G_{1,\rho})}\ +\|v\|_{\i}
\right)^{1-\delta}\cdot \nonumber\\
&&\ \ \ \ \ \ \ \ \ \cdot\left(\|\varphi\|_{L^2(\G_{1,\rho})} +\
\|\nabla v\|_{L^2(\G_{1,\rho})}\right)^{\delta}
\end{eqnarray}
where $\rho \in \left(\frac{M}{4\sqrt{1+M^2}}r_0,\frac{3M}{4\sqrt{1+M^2}}r_0 \right)$ and $P_0=P_1+\frac{M}{4\sqrt{1+M^2}}r_0\cdot \nu$, where $\nu$ is the outer unit normal to $\O$ at $P_1$ and $C_1>0$, $0<\eta<1$ are constants depending on
 $ \rho,r_0, n, M $ only.
For a  proof we refer to \cite{tr}, see also \cite{payne}.
Let us define the function $\tilde{h} \in L^2(\partial \O)$ as follows
$$
\tilde{h}(x)=
\left\{
\begin{array}{ll}
h(x), &\quad\mbox{for a.e. $x \in \G_{1,r_0}$}\ , \\
\displaystyle-\frac{1}{|\G_2|}\int_{\G_{1,r_0}}h\ , & \quad\mbox{for a.e. $x \in \G_2$}\ ,\\0\ , & \quad\mbox{otherwise .}
\end{array}
\right.
$$
Let us consider the following Neumann problem
\begin{equation}\label{z}
 \left\{
\begin{array}
{lcl}
\Delta z=0\ ,& \mbox{in $\Omega$}\ ,
\\
\dfrac{\partial z}{\partial \nu}=\tilde{h}\ , & \mbox{on $\partial \O$ .}
\end{array}
\right.
\end{equation}
Note that $\int_{\partial \O}\tilde{h}=0$, hence a weak solution $z\in H^1(\O)$ exists and it is unique up to an additive constant.
We select the solution $z$ of \eqref{z} with zero average, it is well-known that for such a $z$ the following holds
$$\|z\|_{\i}\le C_2\|\tilde{h}\|_{L^2(\partial \O)}\le C_3\|h\|_{L^2(\G_{1,r_0})}$$
where $C_2$ and $C_3$ are positive constants depending on the \emph{a priori data} only.
Let us set $w=v-z$, thus $w$ solves the following Cauchy problem
\begin{equation}\label{w}
 \left\{
\begin{array}
{lcl}
\Delta w=0\ ,& \mbox{in $\Omega$},
\\
w=\varphi-z\ , & \mbox{on $\G_{1,r_0}$},
\\
\dfrac{\partial w}{\partial \nu}=0\ , & \mbox{on $\G_{1,r_0}$.}
\end{array}
\right.
\end{equation}
By a standard boundary regularity estimate (see for instance
\cite[p.667] {adn}), we have that $w \in
C^{1,\beta}(U^1_{\frac{3}{2}r_0})$ and the following holds
\begin{eqnarray}\label{beta}
\|w\|_{C^{1,\beta}(U^1_{\frac{3}{2}r_0})}\le C_4\|w\|_{\i}
\end{eqnarray}
where  $0<\beta<1$ and $C_4>0$ are constants depending on
$r_0,M,\a$ only. By an interpolation inequality, (see for instance
\cite[p.777]{ABRV}) we have that
\begin{eqnarray}\label{interpol}
\|\nabla w\|_{{L^2}(\G_{1,2r_0})}\le C_5\|w\|^{1-\gamma}_{C^{1,\beta}(U^1_{\frac{3}{2}r_0})}\|w\|^{\gamma}_{L^2(\G_{1,2r_0})}
\end{eqnarray}

where $C_5>0$ and $0<\gamma<1$ are constants depending on $M,\a,r_0$ only.

Moreover,
\begin{eqnarray}\label{distriang}
\|w\|_{\i}\le \|v\|_{\i}+\|z\|_{\i}\le C_6\left(\|v\|_{\i}+\|h\|_{L^2(\G_{1,r_0})}\right)
\end{eqnarray}
where $C_6=\max\{1,C_3 \}$.
From \eqref{beta},\eqref{interpol} and \eqref{distriang} it follows that
\begin{eqnarray}\label{nonso}
\|\nabla w\|_{L^{2}(\G_{1,2r_0})}\le C_7\left(\|v\|_{\i}+\|h\|_{L^2(\G_{1,r_0)})}\right)^{1-\gamma}\cdot\\\left(\|\varphi\|_{L^2(\G_{1,2r_0})}+\|z\|_{L^2(\G_{1,2r_0})} \right)^{\gamma}\ .
\end{eqnarray}

Applying \eqref{grad} to $w$ and using \eqref{nonso} we obtain
\begin{eqnarray*}
\|v\|_{L^2(B_{\rho}(P_0)\cap U^1_{2r_0})}&\leq&C \left(\|\varphi\|_{L^2(\G_{1,r_0})} +\
\|h\|_{L^2(\G_{1,r_0})}\ +\|v\|_{\i}
\right)^{1-\gamma\eta}\cdot\\
&&\cdot\left(\|\varphi\|_{L^2(\G_{1,r_0})} +\
\|h\|_{L^2(\G_{1,r_0})}\right)^{\gamma\eta}\ .
\end{eqnarray*}
And the thesis follows with $\delta=\gamma\eta$.
\end{proof}

\begin{lemma}[H\"older regularity at the boundary]\label{localbound}
Let $u$ be a solution of \eqref{P}, satisfying the  a priori bound \eqref{E}
then there exists a constant $C>0$, depending on the \emph{a priori data} only, such that
\begin{eqnarray}\label{x}
 \|u\|_{L^{\infty}(B_{\frac{r_0}{4}}(z_0)\cap \O)}\le CE ,\  \  \mbox{for every }  \ z_0\in \G_1
\end{eqnarray}
and
\begin{eqnarray}\label{y}
\|u\|_{C^{0,\a}(\G_1)}\le C E
\end{eqnarray}
where $0<\a<1$ is a constant depending on $r_0,M,n$ only.
\end{lemma}

\begin{proof}
From the weak formulation of the problem \eqref{soluzione} we have
\begin{eqnarray*}
\int_{B_{\frac{r_0}{2}}}\nabla u\cdot \nabla \varphi=\int_{\G_{\frac{r_0}{2}}}f(u)\varphi
\end{eqnarray*}
where $B_{\frac{r_0}{2}}=B_{\frac{r_0}{2}}(z_0) \cap \O$ and $\G_{\frac{r_0}{2}}=\partial
B_{\frac{r_0}{2}} \cap \G_1$ and $\varphi$ is any test function in $\i$ such that $\mbox{supp}{\varphi}\subset B_{\frac{r_0}{2}} \cup \G_{\frac{r_0}{2}}$.\\
By \eqref{f} we have that
\begin{eqnarray*}
\bigg|\int_{B_{\frac{r_0}{2}}}\nabla u\cdot \nabla \varphi\bigg|\le L\int_{\G_{\frac{r_0}{2}}}|u\varphi|
\end{eqnarray*}
and by a trace inequality (see \cite[p.114]{Adams}) it follows that
\begin{eqnarray}\label{traccia}
\bigg|\int_{B_{\frac{r_0}{2}}}\nabla u\cdot \nabla \varphi\bigg|\le CL \int_{B_{\frac{r_0}{2}}}|\nabla (u\varphi)|.
\end{eqnarray}
From now on one can proceed by the standard iteration technique due to Moser (see for instance \cite {gt}) and \eqref{x} follows. By the local bound \eqref{x} and by applying again the method by Moser leading to the Harnack inequality (see \cite{gt}) we obtain \eqref{y}.
 \end{proof}

\begin{theorem}[$C^{1,\a}$ regularity at the boundary]\label{teoreg}
Let $u$ be a solution of \eqref{P}, satisfying the  a priori bound \eqref{E}, then $u \in C^{1,\alpha}(G)$ and there exists a constant $C_{\rho}>0$, depending on the \emph{a priori data} and on $\rho$ only, such that the following estimate holds
\begin{eqnarray}\label{stima2}
\|u \|_{C^{1,\alpha}(U^1_{\rho})}\le C_{\rho}E
\end{eqnarray}
where  $\rho \in (0,r_0)$ .
\end{theorem}
\begin{proof}
 Since, by Lemma \ref{localbound}, we know that $u \in  C^{0,\alpha}(\G_1)$, by the Lipschitz
 regularity of $f$ we have that
$$\frac{\partial u}{\partial \nu}(x)=f(u(x))\ \ \in  C^{0,\alpha}(\G_1)\ . $$
By well-known regularity bounds for the Neumann problem (see for instance \cite[p.667]{adn}) it follows that $u \in C^{1,\alpha}(U^1_{\rho})$ and the following estimate holds
\begin{eqnarray}\label{stimaadn}
\|u\|_{C^{1,\alpha}(U^1_{\rho})}&\le& C\left(\|u\|_{C^{0,\alpha}(\G_{1,\frac{\rho}{2}})}+\left\|\frac{\partial u}{\partial \nu}\right\|_{C^{0,\alpha}(\G_{1,\frac{\rho}{2}})}+\|\nabla u\|_{L^2(\O)}\right)\le\nonumber\\
&\le&C\left(\left\|\frac{\partial u}{\partial \nu}\right\|_{C^{0,\alpha}(\G_{1,\frac{\rho}{2}})}+E \right)
\end{eqnarray}
where $C>0$ depends on the \emph{a priori data} and on $\rho$ only.
Moreover, we can estimate the $C^{0,\alpha}$ norm of $\displaystyle\frac{\partial u}{\partial \nu}$ in terms of $E$, in fact

\begin{eqnarray}
\left\|\frac{\partial u}{\partial \nu}\right\|_{C^{0,\alpha}(\G_{1,\frac{\rho}{2}})}=
\sup_{x \in \G_{1,\frac{\rho}{2}}}\left|\frac{\partial u(x)}{\partial \nu}\right|+\left({\frac{\rho}{2}}\right)^{\a}\sup_{x,y \in \G_{1,\frac{\rho}{2}}}\frac{\left|\frac{\partial u(x)}{\partial \nu}-\frac{\partial u(y)}{\partial \nu}\right|}{|x-y|^{\alpha}}=\nonumber \\
=\sup_{x \in \G_{1,\frac{\rho}{2}}}\left|f(u(x))\right|+\left({\frac{\rho}{2}}\right)^{\a}\sup_{x,y \in \G_{1,\frac{\rho}{2}}}\frac{\left|f(u(x))-f(u(y))\right|}{|x-y|^{\alpha}}\ .\nonumber
\end{eqnarray}
 By the Lipschitz bound \eqref{Lip} on $f$ and by Lemma \ref{localbound} we obtain
\begin{eqnarray}
\left\|\frac{\partial u}{\partial \nu}\right\|_{C^{0,\alpha}(\G_{1,\frac{\rho}{2}})}&\le& L\sup_{x \in \G_{1,\frac{\rho}{2}}}|u(x)|+L\left({\frac{\rho}{2}}\right)^{\a}\sup_{x,y \in \G_{1,\frac{\rho}{2}}}\frac{|u(x)-u(y)|}{{|x-y|}^{\alpha}}\le\nonumber\\
&\le& L\|u \|_{C^{0,\alpha}(\G_1)}\le CE\ .
\end{eqnarray}
So inserting this estimate in \eqref{stimaadn} we have the thesis.
\end{proof}
\begin{corollary}\label{corollario}
Let $u$ be as above, then, for every $\rho>0$, the function $\frac{\partial u}{\partial \nu}$ belongs to $C^{0,1}(\G_{1,\rho})$, with Lipschitz constant $\tilde{L}$ depending on the a \emph{a priori data} and on $\rho$ only.
\end{corollary}
\begin{proof}
Let $x$ and $y$ be two points in $\G_{1,\rho}$ then, by the assumption \eqref{Lip} and by Theorem \ref{teoreg}, it follows that
\begin{eqnarray*}
\bigg|\frac{\partial u(x)}{\partial \nu}-\frac{\partial u(y)}{\partial \nu}\bigg|&=& |f(u(x))-f(u(y))|\le L |u(x)-u(y)|\le\\
&\le& LC_{\rho}E|x-y|\ .
\end{eqnarray*}
The thesis follows with $\tilde{L}= LC_{\rho}E$.
 \end{proof}

\begin{proofb}
Let  $\eps=\mathop{\rm osc\ u}\limits_{\G_1}>0 $, since $u=0$ on $\G_D$, we have that
\begin{eqnarray}\label{bu}
 \|u\|_{L^2({\G_{1,r_0})}}\le C_1\eps
\end{eqnarray}
where $C_1$ is a positive constant depending on the \emph{a priori
data} only. By the a priori assumption \eqref{f} on $f$, we have
that  $|f(u)|\le L\eps $, moreover, since $$\left|\frac{\partial
u(x)}{\partial \nu}\right|=|f(u(x))|\ \  \mbox{on   $\G_1$},$$
then
\begin{eqnarray}\label{bu2}
\left\|\frac{\partial u}{\partial \nu}\right \|_{L^2(\G_{1,r_0})}\le |\G_{1,r_0}|^{\frac{1}{2}}L\eps\ .
\end{eqnarray}
By Proposition \ref{tritten}, it follows
\begin{equation}\label {tri}
 \|u\|_{L^2(B_{\rho}(P_0)\cap U^1_{2r_0})}\le C (\eps+E)^{1-\d}\cdot \eps ^{\d}
\end{equation}
where $C$ is a constant depending on the \emph{a priori data} only.
Since the boundary of $\O$ is of Lipschitz class, then it satisfies the cone property.
More precisely, if $Q$ is a point of $\partial \O$, then there exists a rigid transformation of coordinates under which we have $Q=0$. Moreover, considering the finite cone
$$\mathcal{C}=\bigg\{x :|x|<r_0 ,\,\,\frac{x\cdot \xi}{|x|}>\cos\theta \bigg \} $$
with axis in the direction $\xi$ and width $2\theta$, where $\theta=\arctan \frac {1}{M}$, we have that $\mathcal{C}\subset \O$. Let us consider now a point $Q \in \G_{2,r_0}$  and let $Q_0$ be a point lying on the axis $\xi$ of the cone with vertex in $Q=0$ such that $d_0=\mbox{dist}(Q_0,0)<\frac{r_0}{2}$. Following Lieberman \cite{Li}, we introduce a regularized distance $\tilde{d}$ from the boundary of $\O$. We have that there exists $\tilde{d}$ such that $\tilde{d}\in C^2(\O)\cap C^{0,1}({\bar\O})$, satisfying the following properties
\begin{itemize}
\item $\gamma_0\le \displaystyle \frac{\mbox{dist}(x,\partial \O)}{\tilde{d}(x)}\le \gamma_1$,
\item $|\nabla {\tilde{d}}(x)|\ge c_1$,\ \ \ for every $x$ such that ${\mbox {dist}}(x,\partial \O)\le br_0$,
\item $\|\tilde{d}\|_{C^{1,\a}}\le c_2r_0 $,
\end{itemize}
where $\gamma_0, \gamma_1, c_1, c_2, b$ are positive constants depending on $M, \a$ only, (see also \cite[Lemma 5.2]{ABRV}).

Let us define for every $\rho>0$
$$\tilde\O_{\rho}=\{x\ \in \O\ \ : \ \ {\tilde{d}}(x)>\rho \}\ . $$
It follows that, there exists $a$, $0<a\le1$, only depending on $M, \a$ such that for every $\rho$, $0<\rho\le ar_0$, $\tilde\O_{\rho}$ is connected with boundary of class $C^1$
and
\begin{eqnarray}\label{dc}
{\tilde c_1}\rho\le \mbox{dist}(x,\partial\O)\le \tilde c_2\rho \ \ \ \ \mbox{for every}\ \ x \ \in \ \partial\tilde\O_{\rho} \cap \O
\end{eqnarray}
where $\tilde c_1, \tilde c_2$ are positive constants depending on $M,\a$ only.
By \eqref{dc} it follows that
$$\O_{\tilde c_2\rho}\subset \tilde \O_{\rho}\subset \O_{\tilde c_1\rho}\ . $$
Using the notation introduced in the Proposition \ref{tritten}, we define the point $P=P_0-\frac{1}{4\sqrt{1+M^2}}r_0\cdot\nu$ and  ${\rho}_0=\min\{\frac{1}{32M\sqrt{1+M^2}}r_0,\frac{r_0}{4}\sin \theta\}$.
Moreover, let  $\gamma$ be a path in $\tilde\O_{\frac{\rho_0}{\tilde c_1}} $ joining  $P$ to  $Q_0$ and let us define $\{y_i\}$, $ i=0,\ldots,s$
as follows $y_0=Q_0$, $y_{i+1}=\gamma(t_{i})$, where $t_i=\max\{\mbox{$t$ s.t. }|\gamma(t)-y_i|=2\rho_0\}$
if $|P-y_i|>2\rho_0$ otherwise let $i=s$
and stop the process.

Now, we will use the three spheres inequality for harmonic functions (see for instance \cite{KM} or \cite[Appendix E]{ab}) that is
$$ \int_{B_{3\rho_0}(y_0)} u^2 \le \left(\int_{B_{\rho_0}(y_0)} u^2\right)^{\tau}\cdot \left(\int_{B_{4\rho_0}(y_0)} u^2\right)^{1-\tau}$$

where $0<\tau<1$ is an absolute constant. Now since $B_{\rho_0}(y_0)\subset B_{3\rho_0}(y_1)$ and since, by hypothesis  $\|u\|_{\i} \le E$, then we have

$$ \int_{B_{\rho_0}(y_0)} u^2 \le \left(\int_{B_{3\rho_0}(y_1)} u^2\right)^{\tau}\cdot E^{1-\tau}\ .$$
An iterated application of the three spheres inequality leads to
$$  \int_{B_{\rho_0}(y_0)} u^2 \le   \left(\int_{B_{\rho_0}(y_s)} u^2\right)^{{\tau}^{s}}\cdot E^{1-{\tau}^s}\ .$$
Finally, since we have $B_{\rho_0}(y_s)\subset B_{\frac{3M}{4\sqrt{1+M^2}}r_0}(P_0)\cap U^1_{2r_0}$, then by the Proposition \ref{tritten} it follows
$$ \int_{B_{\rho_0}(y_0)} u^2 \le C\big\{ (\eps+E)^{1-\d}\cdot(\eps )^{\d}\big\}^{{\tau}^s}\ .$$
We shall construct a chain of balls $B_{\rho_k}(Q_k)$ centered on the axis of the cone, pairwise tangent to each other and all contained in the cone
$${\mathcal{C}}^{\prime}=\bigg\{x :|x|<r_0 ,\,\,\frac{x\cdot \xi}{|x|}>\cos{\theta}^{\prime} \bigg \} $$
where $\theta^{\prime}=\arcsin\big(\frac{\rho_0}{d_0}\big) .$
Let $B_{\rho_0}(Q_0)$ be the first of them, the following are defined by induction in such a way
\begin{equation*}
\begin{array}{l}
Q_{k+1}=Q_{k}-(1+\mu)\rho_k\xi\ ,
\\
\rho_{k+1}=\mu\rho_{k}\ ,
\\
d_{k+1}=\mu d_{k}\ ,
\end{array}
\end{equation*}
with
\begin{equation*}
\begin{array}{l}
 \mu=\dfrac{1-\sin\theta^{\prime}}{1+\sin\theta^{\prime}}\ .
\end{array}
\end{equation*}
Hence, with this choice, we have $\rho_{k}=\mu^{k}\rho_0$ and  $B_{\rho_{k+1}}(Q_{k+1})\subset B_{3\rho_k}(Q_k)$.\\
Let us now consider the following estimate obtained by a repeated application of the three spheres inequality
\begin{eqnarray}\label{s}
\| u \|_{L^2(B_{\rho_k}(Q_k))}&\le& \| u \|_{L^2(B_{3\rho_{k-1}}(Q_{k-1}))
}\le \nonumber\\
&\le&  \| u \|^{\tau}_{L^2(B_{\rho_{k-1}}(Q_{k-1}))}  \| u \|^{1-\tau}_{L^2(B_{4\rho_{l-1}}(Q_{k-1}))}\nonumber\\
&\le& C \| u  \|^{{\tau}^k}_{L^2(B_{{\rho}_0}(Q_0))}\le\nonumber\\
&\le& C\Big\{ \big[ (\eps+E)^{1-\d}\cdot(\eps )^{\d}\big]^{{\tau}^s}\Big\}^{{\tau}^k}\ .
\end{eqnarray}
For every $r$, $0<r<d_0$, let $k(r)$ be the smallest positive integer such that  $d_k\le r$, then since $d_k={\mu}^k d_0$, it follows
\begin{eqnarray}\label{r}
\dfrac{|\log(\frac{r}{d_0})|}{\log{\mu}}\le k(r)\le \dfrac{|\log(\frac{r}{d_0})|}{\log{\mu}} +1
\end{eqnarray}
and by $(\ref{s})$, we have
\begin{eqnarray}\label{k}
   \| u \|_{L^2(B_{\rho_k(r)}(Q_k(r)))}\le
  C  \Big\{ \big[ (\eps+E)^{1-\d}\cdot(\eps)^{\d}\big]^{{\tau}^s}\Big\}^
{{\tau}^{k(r)}}\ .
\end{eqnarray}
Since, by hypothesis, $\G_2$ is contained in a $C^{1,\a}$ surface and by the regularity assumption \eqref{G} on $g$, it follows, by the same argument used in Theorem \ref{teoreg}, that  $u \in  C^{1,\a}(U^2_{2r_0}) $.

Let $\bar{x} \in \G_{2,2r_0}$, $x \in B_{\frac{\rho_{k(r)-1}}{2}}(Q_{k(r)-1}) $,  since $u \in
 C^{1,\a}(U^2_{2r_0})$ we have
$$\left|\frac{\partial u(\bar{x})}{\partial \nu}\right| \le
 \left|\dfrac{\partial u({x})}{\partial \nu}\right| +C|x-\bar{x}|^{\a} \le
\left|\dfrac{\partial u({x})}{\partial \nu}\right| +C\bigg( {\frac{2}{\mu} r}\bigg)^{\a}\ .$$
Integrating over $ B_{\frac{\rho_{k(r)-1}}{2}}(Q_{k(r)-1}) $, we deduce that
\begin{eqnarray}
\left|\frac{\partial u(\bar{x})}{\partial \nu}\right|^2&\le&
\frac{2}{{\omega_n}{(\frac{\rho_{k-1}}{2})}^n}\int_{B_{\frac{\rho_{k(r)-1}}{2}}\big(Q_{k(r)-1}\big)}
\bigg| \frac{\partial u(x)}{\partial \nu} \bigg|^2 \mbox{d}x + 2 C^2{\bigg(\frac{4 r^2}{{\mu}^2}\bigg)}^{\a}\le \nonumber\\
&\le&
\frac{2}{{\omega_n}{(\frac{\rho_{k-1}}{2})}^n}\int_{B_{\frac{\rho_{k(r)-1}}{2}}\big(Q_{k(r)-1}\big)}
|\nabla u(x)|^2 \mbox{d}x + 2 C^2{\bigg(\frac{4 r^2}{{\mu}^2}\bigg)}^{\a}\ .\nonumber
\end{eqnarray}
Applying the Caccioppoli inequality, we have
$$\left|\frac{\partial u(\bar{x})}{\partial \nu}\right|^2\le \frac{C}{{\big(\rho_{k-1}\big)}
^{n+2}}\int_{ B_{\rho_{k(r)-1}}(Q_{k(r)-1})} u(x)^2 \mbox{d}x +C r^{2\a} $$
and since $k$ is the smallest integer such that $d_k\le r$, then $d_{k-1}>r$, it follows

$$\left|\frac{\partial u(\bar{x})}{\partial \nu}\right|^2\le\frac{C}{\big({r \sin {\theta}^{\prime}\big)}
^{n+2}}  \int_{ B_{\rho_{k(r)-1}}(Q_{k(r)-1})} u(x)^2  \mbox{d}x + C r^{2\a}\ .$$ From $\eqref {k}$, we deduce
$$\left|\frac{\partial u(\bar{x})}{\partial \nu}\right|^2\le
\frac{C}{{r}^{n+2}} \Big\{ \big[ (\eps+E)^{1-\d}\cdot(\eps)^{\d}\big]^{{\tau}^s}
\Big\}^
{{\tau}^{k(r)-1}} +Cr^{2\a}\ .
$$
Let us define
$$ \sigma (\eps)= \big[(\eps+E)^{1-\d}\cdot(\eps)^{\d}\big]^{{\tau}^s},$$
thus the previous inequality becomes
\begin{eqnarray*}
\left|\frac{\partial u(\bar{x})}{\partial \nu}\right|^2\le
\frac{C}{{r}^{n+2}} \big\{\sigma (\eps)
\big\}^
{{\tau}^{k(r)-1}} +C r^{2\a}\ .
\end{eqnarray*}

Now, using $\eqref{r}$, we have
$$\tau^{k(r)-1}\ge \bigg(\frac{r}{d_0}\bigg)^{\nu} $$
where $\nu= -\log\big(\frac{1}{\mu}\big)\log\tau$.
We have
$$ \left|\frac{\partial u(\bar{x})}{\partial \nu}\right|\le C \bigg\{r^{-{\frac{n+2}{2}}}\Big[\sigma (\eps)\Big]^{\frac{r^{\nu}}{2}}+r^{\a} \bigg\}\ .$$
Now minimizing the function on the right hand side, with respect to $r$, with $r \in (0,\frac{r_0}{4})$, we deduce
$$\left|\frac{\partial u(\bar{x})}{\partial \nu}\right|\le C{\bigg(\log \frac{1}{\sigma(\eps)}\bigg)}
^{-\frac{2\a}{\nu+2}}\ . $$
Since this estimate holds for every $\bar{x} \in \G_{2,2r_0}$, we infer
$$ \Big\|\frac{\partial u}{\partial \nu}\Big\|_{{L^\infty}(\G_{2,2r_0})}\le C{\Big(\log{
\frac{1}{\sigma(\eps)}}\Big)}^{-\frac{2\a}{\nu+2}}$$
where $C$ is a constant depending on the \emph{a priori data} only.
Hence, solving for $\eps$, we can compute
$$\eps\ge C\exp\Big\{ {-{\Big\|\frac{\partial u}{\partial \nu}\Big\|_{L^\infty
(\G_{2,2r_0})}}}^{\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!-\frac{\nu+2}{2\a}}\Big\}\ .
$$ Note that, recalling the a priori bound \eqref{G}, and choosing
$c=2(1-\log C G^{\gamma})$ and $\gamma=\frac{\nu+2}{2\a}$ one
trivially obtains
$$\eps\ge \exp\bigg[{-\bigg(\frac{t}{c}\bigg)^{-\gamma}}\bigg],\ \ \ \mbox{for every}\ \ t\in (0,G]\ . $$
\end{proofb}

\section{A stability result}

\begin{theorem}[Stability for a Cauchy problem]\label{stabilita'}
Let $\O$, $f_i\  i=1,2$ and $g_i$ satisfy the a priori assumptions
described above. Let $u_i \in \i$, $i=1,2$ be weak solutions of
the problem $\eqref{P}$, with $f=f_i$ and $g=g_i$
respectively and such that \eqref{E} holds for each $u_i$.\\
Moreover, let $\psi_i=u_i\big|_{\G_2},\ \ i=1,2$.
Suppose that
\begin{eqnarray*}
&& \|\psi_1 -\psi_2  \|_{L^{2}(\G_2)}\le \eps\ , \\
&& \|g_1 -g_2  \|_{L^{2}(\G_2)}\le \eps\ ,
\end{eqnarray*}
 then, for every $\rho \in (0,r_0)$
\begin{eqnarray}\label{stab}
  \| u_1 -u_2 \|_{C^1(\G_{1,\rho)}}\le \omega(\eps)
\end{eqnarray}
where $\omega$ is given by \eqref{omega} with a constant $C>0$ which depends on the \emph{a priori data} and  on $\rho$ only.
\end{theorem}

\begin{proof}
Arguing as in Theorem \ref{teo1}, we find the following estimate
\begin{eqnarray*}
&&\|u_1 -u_2 \|_{L^{\infty}(\G_{1,\frac{\rho}{2}})}\le
\tilde{\omega}(\eps)\ ,\\
&&\left\|\frac{\partial{u_1}}{\partial\nu} -\frac{\partial{u_2}}{\partial{\nu}}\right \|
_{L^{\infty}(\G_{1,\frac{\rho}{2}})}\le
\tilde{\omega}(\eps)
\end{eqnarray*}
where $\tilde{\omega}$ is a positive increasing function of the type \eqref{omega}, such that
$$\tilde
{\omega}(t)\le\tilde{C}|\log(t)|^{-\gamma_1}\ \ \mbox{for every}\ \ 0<t<1$$
where $\tilde{C}>0, 0<\gamma_1<1$ are constants depending on the \emph{a priori data} and on $\rho$ only.
By an interpolation inequality we have
$$\|\nabla_t(u_1-u_2)\|_{ L^{\infty}(\G_{1,\rho})}\le C\|u_1-u_2 \|_{ L^{\infty}(\G_{1,\frac{\rho}{2}})} ^{\beta}{\|u_1-u_2 \|_{C^{1,\a}(\G_{1,\rho})}}^{1-\beta}$$
where $\beta=\frac{\a}{\a+1}$ and $C>0$ depends on the \emph{a priori data} and on $\rho$ only,
thus by Theorem \ref{teoreg} it follows that
$$\|\nabla_t(u_1-u_2)\|_{ L^{\infty}(\G_{1,\rho})}\le C\|u_1-u_2 \|_{ L^{\infty}(\G_{1,\frac{\rho}{2}})} ^{\beta}{E}^{1-\beta}$$
where $C>0$ only depends on the \emph{a priori data} and on $\rho$ only.

It follows that for every $\eps<\eps_0$, with $\eps_0$ depending only on the \emph{a priori data}
\begin{eqnarray}
&&\|\nabla(u_1-u_2)\|_{ L^{\infty}(\G_{1,\rho})}\le\nonumber
\\
&&\le\left\| \frac{\partial{u_1}}{\partial{\nu}} -\frac{\partial{u_2}}{\partial{\nu}}\right \|
_{L^{\infty}(\G_{1,\rho})} +
\|\nabla_t(u_1-u_2)\|_{ L^{\infty}( \G_{1,\rho})}\le\nonumber
\\
&&\leq\tilde{C}\tilde{\omega}(\eps)^{\beta}
\end{eqnarray}
where $\tilde{C}>0$ depends on the \emph{a priori data} only.
Hence,
\begin{eqnarray}
\|u_1-u_2 \|_{ L^{\infty}(\G_{1,\rho})}+
\|\nabla(u_1-u_2)\|_{ L^{\infty}(\G_{1,\rho})}\le \omega(\eps)
\end{eqnarray}
where
$$\omega(t)\le C {\tilde{\omega}(t)}^{\b}\ \ \ \mbox {for every}\ \ 0< t<1\ .$$

\end{proof}
\begin{proposition}[Local monotonicity]\label{inv}
Let $u$
be a solution of \eqref{P} satisfying \eqref{E}, then there exist a point $\bar{x}
\in \G_{1,\tau}$
and a direction  $\xi \in {\mathbb{R}}^{n-1} ,\ |\xi|=1$ such that, in the representation \eqref{si} of $\G_{1}$ near $\bar{x}$, the following holds
\begin{eqnarray*}
\left|\nabla _{x'}u(x',\varphi(x'))\cdot \xi \right|\ge {\eta\left(\|g\|_{L^{\infty}(\G_{2,2r_0})} \right)} ,\
 \ x'\in U_{\bar{x}'}=\{x'=t\cdot\xi+\bar{x}',\ |t|\le\tau\}
\end{eqnarray*}
with
\begin{eqnarray}\label{tau}
\tau=\min\left\{\frac{r_0}{4} ,\frac{a \tilde{c_1}r_0}{4}, {\eta(\|g\|_{L^{\infty}(\G_{2,2r_0})})}\right\}
\end{eqnarray}
where $0<a<1, \tilde{c_1}>0$ are constants depending on the \emph{a priori data} only and $\eta$ satisfies \eqref{eta}.
\end{proposition}
\begin{proof}
Arguing as in Theorem \ref{teo1}, we can introduced a regularized distance, in the sense of Lieberman, on $S_1$ from the boundary of $\G_1$ and consequently construct connected sets $\tilde\G_{1,\rho}$ for every $\rho$, $0<\rho\le a r_0$, which satisfy
\begin{eqnarray}\label{distanza}
\G_{1,\tilde c_2h}\subset \tilde \G_{1,h}\subset \G_{1,\tilde c_1h}
\end{eqnarray}
where $0<a<1,\tilde c_2>\tilde c_1>0$ are constants depending on $M,\a$ only.

Since, by Lemma \ref{localbound}, $u \in C^{0,\a}(\G_1)$, we have that  by \eqref{distanza} it follows
$$\mathop{\rm osc}\limits_{\tilde\G_{1,\frac{\rho}{\tilde{c_1}}}} u\ge\mathop{\rm osc}\limits_{\G_{1,\frac{\tilde {c_2}\rho}{\tilde{c_1}}}} u\ge  \mathop{\rm osc}\limits_{\G_{1}} u - 2{C}E\left({\frac{\rho}{\tilde{c_1}}}\right)^{\a}{\tilde{c_2}}^{\a}\ .$$
Moreover by Theorem \ref{teo1}, we infer that
$$\mathop{\rm osc}\limits_{\tilde\G_{1,\frac{\rho}{\tilde{c_1}}}} u \ge \eta(\|g\|_{L^{\infty}(\G_{2,2r_0})})-2{C}E\left({\frac{\rho}{\tilde{c_1}}}\right)^{\a}{\tilde{c_2}}^{\a}\ .$$
Possibly replacing $c$ by a larger constant in \eqref{eta} and  taking
$$r_1=\min \left\{\eta(\|g\|_{L^{\infty}(\G_{2,2r_0})}),a\tilde{c_1}r_0, r_0\right\}$$
we have that
\begin{eqnarray}\label{osc}
\mathop{\rm osc}\limits_{\tilde\G_{1,\frac{r_1}{\tilde{c_1}}}} u \ge\eta(\|g\|_{L^{\infty}(\G_{2,2r_0})})\ .
\end{eqnarray}
Let us set, for simplicity,
$\eta=\eta\left(\|g\|_{L^{\infty}(\G_{2,2r_0})} \right)$. Since in
the a priori assumptions we have assumed that the portion $\G_1$
of the boundary is of $C^{1,\a}$ class, then we can locally
represent the restriction of $u$ (the solution to \eqref{P}) to
$\G_1$, as a function of $n-1$ variables, more precisely, for
every $x_0 \in \G_1$, up to a rigid change of coordinates, we
denote
\begin{eqnarray}\label{rapr}
w(x')= u(x',\varphi_1(x'))
 \ \ \mbox{for all } \ x\ \in \G_1 \cap B_{r_0}(x_0).
\end{eqnarray}
 By \eqref{osc}, it follows that exist two points
 $x$
and $y$ in  $ \tilde\G_{1,\frac{r_1}{\tilde{c_1}}}$,
 such that
\begin{eqnarray}\label{K}
\eta\le u(x)-u(y)\ .
\end{eqnarray}
Let us consider a continuous path $\sigma \subset \tilde\G_{1,\frac{r_1}{\tilde{c_1}}}$
joining $x$
to $y$
and let us define a sequence $\{x_i\}_{i=0,\dots,l}$ as follows
$x_0=x,\ x_i=\sigma(s_i)$ where $s_i=\max\{\mbox{$t$ },|\sigma(s)-x_i|
=\frac{r_1}{4}\}$ if $|y-x_i|>\frac{r_1}{4}$ otherwise let $i=l$ and otherwise stop the process.

The number $l$ of balls is bounded from above by $CM {\left(\frac{D}{r_1}\right)}^{n-1}$, where $C>0$ is a constant depending on $n$ only.

Let us define
$$M_i=\!\!\!\!\!\!\!\!\max_{{\overline{B_{\frac{r_1}{4}}} (x_i)
}\cap \G_1 }\!\!|\nabla_t u(x)| $$ where $\nabla_t$ denotes the
tangential gradient on $\G_1$. Let $\bar{M},\ \bar{\imath},\
\bar{x}$ be such that $\bar{x} \in B_{\frac{r_1}{4}}
(x_{\bar{i}})\cap \G_1$ and
\begin{eqnarray}\label{L}
\bar{M}=\max_{i=1,\dots,l}\left\{M_i\right\}=|\nabla_t u(\bar{x})|.
\end{eqnarray}
 By \eqref{K} and the mean value Theorem, it follows that
\begin{eqnarray*}
&\eta&\le |u(x)-u(x_1)|+\dots+|u(x_l)-u(y)|\le
\\&&\le \sum_{i=1,\dots,l} M_i\frac{r_1}{4}
\le \bar{M} C_{1}
\end{eqnarray*}
where $C_{1}>0$
is a constant depending on the \emph{a priori data} only.
Thus we have
\begin{eqnarray}\label{M}
 \bar{M}\ge \frac{\eta}{C_1}>0\ .
\end{eqnarray}
Now we use the local representation of $u$ as a function of $n-1$ variables \eqref{rapr}, within $\G_1\cap B_{\frac{r_1}{4}} (x_{\bar{i}})$.
Let us define the direction $\xi=\frac{\nabla _{x'}w}{|\nabla _{x'}w|}(\bar{x}')$.
We shall further restrict the function $w$ to the segment $t\cdot\xi+\bar{x}'$, with
$$v(t)= w(t\cdot\xi+\bar{x}')\ .$$
Now, we look for a neighborhood $U_{0}$  of $t=0$ such that

\begin{eqnarray}\label{A}
 |v'(t) |\ge \frac{\eta}{2C_1}\ \ \mbox{for every }\ \ t\in  U_{0}\ .
\end{eqnarray}

It follows that for every $|t|<\frac{r_1}{4} $
$$|v'(0)-v'(t)|\le C_2|t|^{\a}$$
where $C_2>0$ is a constant depending on the \emph{a priori data} only.

Thus  we have
$$ \bar{M}=|v'(0)|\le
 |v'(t)|  + C_2|t|^{\a}\ . $$

Hence by \eqref{M},
$$\frac{\eta}{C_1} -  C_2|t|^{\a}\le|v'(t)|\ .$$
Let us choose $t$ in such a way
$$ C_2|t|^{\a}\le\frac{\eta}{2C_1}\ . $$
Hence \eqref{A} holds with  $U_0=[-\tau,\tau]$,
where $\tau=\min\left\{\frac{r_1}{4} ,\bigg(\frac{\eta}{2C_1C_2} \bigg)^{\frac{1}{{\a}}}\right\}$.
The thesis follows, observing that $v'(t)=\frac{\partial w(x')}{\partial \xi}=\nabla_{x'}u(x',\varphi(x'))\cdot \xi$ and, possibly, by a further adjustment of the constant $c$ in \eqref{eta}.
\end{proof}

\begin{proofa}
Let $\bar{x}\in\G_{1,\tau_1},\tau_1, \xi\in \mathbb{R}^{n-1}$ be the point, the length and the direction introduced in Proposition \ref{inv}, with $u$ replaced with $u_1$. Up to a change of coordinates, we assume $\xi=e_1$.
Let $$v_i(t)=u_i(t\cdot\xi+\bar{x}, \varphi_1(t\cdot\xi+\bar{x}))\ ,\ \ \ i=1,2\ ,$$where $x=(x',\varphi_1(x'))$ is the local representation of $\G_1$ near $\bar{x}$.\\
By Proposition \ref{inv} and assumption \eqref{m}, we have that
\begin{eqnarray}\label{1}
|v_1^{\prime}(t)|\ge \eta(m)\ ,\ \ \ \  \ \mbox{for every } \ t \in U_{0}=[-\tau_1,\tau_1]\ .
\end{eqnarray}
We shall denote by $\eta_1=\eta(\|g_1\|_{L^{\infty}(\G_{2,2r_0})})$.
By the stability estimate \eqref{stab} of Theorem \ref{stabilita'}, we have that $$v'_2(t)\ge \eta_1-\omega(\eps)\ ,\ \mbox{for every}\ t\in U_0\ . $$
Thus choosing $\eps_0$ such that
$$\omega(\eps_0)\le \frac{\eta_1}{2}$$
we have
\begin{eqnarray}\label{2}
 |v_2^{\prime}(t)|\ge \frac{\eta_1}{2}\ ,\  \ \mbox{for every }  t \in U_{0}\ .
\end{eqnarray}

Thus the functions $v_i$
are invertible on $U_0$, let us denote by $V_i$ their respective images and by
\begin{eqnarray}
s^i:V_i\rightarrow U_{0}\ ,\ \ \ i=1,2\ ,
\end{eqnarray}
their inverse functions.
Let us observe that the intervals $V_1$
 and $V_2$ overlap on a sufficiently large interval $V$.
In fact, by \eqref{1} and \eqref{2} it follows that $v_i$ are monotone.Without loss of generality, let us assume they are both increasing. We have that, taken
$$a=-\frac{\tau_1}{2}\ , \ \ b=\frac{\tau_1}{2}\ ,$$
the following hold
\begin{eqnarray*}
v_i(a)<v_i(t)<v_i(b)\ ,\  \ \ \mbox{for every } t\in (a,b)\ ,\ \ i=1,2\ .
\end{eqnarray*}
Moreover, since by the Theorem \ref{stabilita'}
we have
$$ \| u_1 -u_2 \|_{L^{\infty}(\G_{1,\frac{r_1}{2}})}\le \omega(\eps) $$
then, it follows that, for $\varepsilon< \varepsilon_0$, setting $V=(v_1(a)+2\omega(\eps),v_1(b)-2\omega(\eps))$, for every $u \in  V$, there  exists $t \in(a,b) \ \mbox{such that}\ \  v_2(t)=u .$

Let us estimate from below the length of the interval $V$.
By the mean value Theorem, \eqref{1} and \eqref{tau}, it follows that
\begin{eqnarray*}
&&|v_1(a)-v_1(b)|=|v'_1(\xi)||b-a|\ge \eta_1\tau_1\ .
\end{eqnarray*}
Thus the length $\mathcal{L}$
of $V$ is bounded from below by
\begin{eqnarray*}
\mathcal{L}\ge \tau_1\eta_1-\omega(\eps)\ .
\end{eqnarray*}
 Hence, possibly adjusting the constant $c$ in the definition \eqref{eta} of $\eta$, we have that
$$\mathcal{L}\ge \eta(m)-\omega(\eps_0)\ge \frac{1}{2}\eta(m)>0\ . $$

Let us consider any value $u \in V$, then using the inverse function $s^i$,
we have
$$u=v_1(s^1(u))=v_2(s^2(u))\ .$$
Let us estimate
\begin{eqnarray*}
&&\hspace{-1cm}\!|f_1(u)-f_2(u)|=\ \nonumber\\
&&\hspace{-1cm}\!\left|\frac{\partial u_1}{\partial \nu}
(s^1(u)e_1,\varphi_1(s^1(u)e_1))\!-\!\frac{\partial u_2}{\partial \nu}
(s^2(u)e_1,\varphi_1(s^2(u)e_1))\right|\le\nonumber\\
&&\hspace{-1cm}\!\left|\frac{\partial u_1}{\partial \nu}
(s^1(u)e_1,\varphi_1(s^1(u)e_1))-\frac{\partial u_2}{\partial \nu}
(s^2(u)e_1,\varphi_1(s^2(u)e_1))\right|+ \\
&&\hspace{-1cm}\left|\frac{\partial u_2}{\partial \nu}
(s^1(u)e_1,\varphi_1(s^1(u)e_1))-\frac{\partial u_2}{\partial \nu}
(s^2(u)e_1,\varphi_1(s^2(u)e_1))\right|\hspace{+1cm}\textrm{\ \ }
\end{eqnarray*}
where $e_1=(1,0,\cdots,0)\ \in \mathbb{R}^{n-1}$.
By  Theorem
\ref{stabilita'} it follows that, for all $u \in V$,
\begin{eqnarray}\label{stima1}
&& \left|\frac{\partial u_1}{\partial \nu}
(s^1(u)e_1,\varphi_1(s^1(u)e_1))-\frac{\partial u_2}{\partial \nu}
(s^1(u)e_1,\varphi_1(s^1(u)e_1))\right|\le\omega(\eps)\ .\ \ \ \ \ \ \ \ \ \
\end{eqnarray}

 By Corollary \ref{corollario}, we infer that
\begin{eqnarray*}
&&\left|\frac{\partial u_2}{\partial \nu}
(s^1(u)e_1,\varphi_1(s^1(u)e_1))-\frac{\partial u_2}{\partial \nu}
(s^2(u)e_1,\varphi_1(s^2(u)e_1))\right|\le\\
&&\tilde{L}\left(|s^1(u)-s^2(u)|+|\varphi_1(s^1(u)e_1)-\varphi_1(s^2(u)e_1)|   \right)\le\nonumber\\
&& \tilde {L}(1+M)|s^1(u)-s^2(u)|\ .
\end{eqnarray*}
By the mean value Theorem, we find
\begin{eqnarray*}
v_2(s^2(u))=v_2(s^1(u)) + v_2'(\bar{s})(s^2(u)-s^1(u))
\end{eqnarray*}
where $\bar{s}$
is a point between $s^2(u)$ and $s^1(u)$.
Since
$$ v_2(s^2(u))=v_1(s^1(u))\ , $$
by \eqref{2} and by Theorem \ref{stabilita'}, it follows that
\begin{eqnarray}
|s^1(u)-s^2(u)|&\le& \frac{2}{\eta_1}|v_2(s^1(u))-v_1(s^1(u))|\le\nonumber\\
&\le& \frac{2}{\eta_1}\omega(\eps)\ , \ \ \ \ \ \mbox{for every}\ \ u \in V\ .\nonumber
\end{eqnarray}
Finally, we infer that
\begin{eqnarray*}
|f_1(u)-f_2(u)|\le\omega(\eps)\ ,\ \ \ \ \mbox{for every} \ \ u \in\  V\ ,
\end{eqnarray*}
possibly by a further adjustment of the constant $C$ in \eqref{omega}.
\end{proofa}

config: arxiv


TITLE: How can I find the equation of the tangent line to the curve at the point(−3,f(−3))?
QUESTION [0 upvotes]: Let f(t)=2/(t+8)
(a) Simplify the following difference quotient as much as possible:
[f(t+h)−f(t)]/h
MY ANSWER (CORRECT): −2/[(t+8)(t+h+8)]
(b) Use your result from (a) and the limit definition of the derivative to calculate
f'(t)=lim_h→zero [f(t+h)−f(t)]/h=
MY ANSWER (CORRECT): −2/[(t+8)^2]
(c) Use your answer from part (b) to evaluate
f'(−3)
MY ANSWER (CORRECT): −2/25
(d) Find the equation of the tangent line to the curve at the point
(−3,f(−3))
y= ?
I have tried the point-slope formula, y−y_1=m(x−x_1), but my answer was incorrect.

REPLY [1 votes]: Hint: Your equation has the form
$$y=-\frac{2}{25}x+n$$ you you can plug $$x=-3,y=\frac{2}{5}$$ in the equation above and you will get $n$
I have got
$$y=-\frac{2}{25}x+\frac{4}{25}$$


text:
\begin{document}
\preprint{AIP/123-QED}
\title[]{Multi-branched resonances, chaos through quasiperiodicity,
 and asymmetric states in a superconducting dimer
}
\author{J. Shena}

\affiliation{National University of Science and Technology ``MISiS'',
             Leninsky Prospekt 4, Moscow 119049, Russia}
             \author{N. Lazarides}
\affiliation{Department of Physics, University of Crete, 71003 Herakleio, Greece}

\author{J. Hizanidis}
\email{hizanidis@physics.uoc.gr}
\affiliation{Department of Physics, University of Crete, 71003 Herakleio, Greece}

\date{\today}

\begin{abstract}
A system of two identical SQUIDs (superconducting quantum interference devices)
symmetrically coupled through their mutual inductance and driven by a sinusoidal
field is investigated numerically with respect to dynamical properties such as
its multibranched resonance curve, its bifurcation structure, as well as its
synchronization behavior. The {\em SQUID dimer} is found to exhibit
a hysteretic resonance curve with a bubble connected to it through Neimark-Sacker
(torus) bifurcations, along with coexisting chaotic branches in their vicinity.
Interestingly, the transition of the SQUID dimer to chaos occurs through a
period-doubling cascade of a two-dimensional torus (quasiperiodicity-to-chaos
transition). The chaotic states are identified through the calculated Lyapunov
spectrum, and their basins of attraction have been determined. Bifurcation
diagrams have been constructed on the parameter plane of the coupling strength
and the driving frequency of the applied field, and they are superposed to maps
of the maximum Lyapunov exponent on the same plane. In this way, a clear
connection between chaotic behavior and torus bifurcations is revealed. Moreover,
asymmetric states that resemble localized synchronization have been detected
using the correlation function between the fluxes threading the loop of the
SQUIDs. The effect of intermittent chaotic synchronization, which seems to be
present in the SQUID dimer, is only slightly touched.
\end{abstract}
\maketitle

{\bf Networks of coupled limit cycle oscillators represent a class of systems
with special interest in physics, chemistry, and biological sciences. Arrays of
coupled nonlinear oscillators, in particular, which often exhibit complex
dynamical behavior, have attracted large amounts of theoretical, numerical, and
experimental efforts. However, the prominent features of a finite large number
of such oscillators can be sometimes understood by analysing just two coupled
oscillators. A unique in many aspects oscillator is the superconducting quantum
interference device (SQUID) that has been investigated extensively for many
years. The SQUID is a low-loss, highly nonlinear resonant element that responds
strongly to applied magnetic field(s) and exhibits rich dynamical behavior.
Moreover, it has recently been the elementary unit for the construction of
metamaterials in one and two dimensions. SQUIDs and SQUID arrays are
technologically important devices, and they also serve as a testbed for
exploring complex dynamics. Two SQUIDs in close proximity are coupled together
through magnetic dipole-dipole forces, and their dynamical complexity significantly
increases. Specifically, more complex bifurcation scenarios appear, along with
transitions to chaos through quasiperiodicity, the emergence of localized
synchronization, and intermittent chaotic synchronization. These dynamical
properties are explored numerically with a well-established model whose
parameters are acquired from recent relevant experiments.
}

\section{\label{sec:level1} Introduction}
There has been great interest on the behavior of systems with interacting and
forced nonlinear oscillations with applications in physics, biology, and
chemistry, which in addition are in abundance in the natural world
\cite{Awrejcewicz1991}. The interaction introduces qualitatively new behavior
to the network of oscillators, as compared to that of a single oscillator. For
example, synchronization of the oscillating units of the system
\cite{Rosenblum2003}, extreme multistability \cite{Wiesenfeld1989,Nishio1992},
multiple resonance and anti-resonance effects
\cite{Jothimurugan2016,Kominis2019,Sarkar2019},
complex bifurcation structure \cite{Yin1998,Kenfack2003}, localization
\cite{Lazarides2010}, emergence of exceptional points \cite{Kominis2018}, the
intriguing effect of amplitude death \cite{Herrero2000}, chaos synchronization
\cite{Tafo2009},
localized synchronization \cite{Hohl1997}, chaos to hyperchaos transitions
\cite{Kapitaniak1991,Kapitaniak1995,Kapitaniak2000,Grygiel2000}, and
quasi-periodicity with subsequent transition to chaos
\cite{Rand1982,Bishop1986,Kozlowski1995,Dixon1996,Elhadj2008,Borkowski2015},
are some of the experimentally or numerically observed behaviors.

Here we consider a pair of coupled  nonlinear oscillators (i.e., a nonlinear
dimer) that belong to the large class of externally driven and dissipative
systems. More specifically, we consider a pair of identical mutually coupled
SQUID oscillators \cite{Kleiner2004,Fagaly2006}, where the acronym stands for
``superconducting quantum interference device''. SQUIDs are mesoscopic
superconducting devices which are modeled efficiently by equivalent electrical
circuits, while their dynamics is governed by a second order nonlinear ordinary
differential equation. The simplest variant of a SQUID consists of a
superconducting loop interrupted by a Josephson junction (JJ), called rf SQUID,
which was found to exhibit very rich dynamic behavior, including complex
bifurcation structure and chaos \cite{Hizanidis2018}. Moreover, rf SQUIDs are
employed as elementary units for the fabrication of superconducting metamaterials,
whose investigation has revealed several extraordinary properties
\cite{Trepanier2013,Zhang2015}.
Recent theoretical works on SQUID metamaterials have also reported the existence of chimera states and pattern formation~\cite{LAZ15,HIZ16a,HIZ16b,HIZ20}.

\begin{figure}[!h]
   \includegraphics[scale=0.2]{Fig1.jpg}
   \caption{Schematic diagram of a SQUID dimer in a magnetic field where ``Mf''
            is the Magnetic field, ``Sr'' is the Superconducting ring, ``JJ'' is
            the Josephson Junction, and $I_1$, $I_2$ are the induced currents.
            Panels (a) and (b) correspond to negative and positive magnetic
            coupling strength, respectively, appropriate for the planar and the
            axial geometry.}
\label{fig1}
\end{figure}

Two SQUIDs in relatively close proximity are coupled together by magnetic
dipole-dipole forces due to their mutual inductance $\cal M$, while the sign and
the magnitude of their coupling strength coefficient depends strongly on their
relative positions. Specifically, the two SQUIDs may be arranged either in the
planar geometry (Fig. \ref{fig1}(a)), where they lie in the same plane with
their axes in parallel, or in the axial geometry, they lie the one on top of the
other with their axes lying on the same line (Fig. \ref{fig1}(b)). The SQUIDs are
driven by externally applied periodic and constant (bias) magnetic fields of
appropriate polarization, to which they respond resonantly. Although its
potential for technological applications, the dynamics of a SQUID dimer are not
properly investigated. Recently however, the existence of homoclinic chaos in a
pair of parametrically driven SQUIDs has been exposed
\cite{Agaoglou2015,Agaoglou2017}.

In the present work, particular dynamical properties of a SQUID dimer are
addressed, such as its highly unusual resonance curve that consists of a
bistable part with a bubble connected to it through Neimark-Sacker (NS) and
pitchfork bifurcations (PBs). These bubbles are of very different character compared to those reported
long ago by Bier and Bountis~\cite{BIE84} and they deserve to
be investigated on their own. Moreover, in addition to those solutions, there
are at least two more coexisting solutions in the frequency interval spanned by
the resonance curve, which end up becoming chaotic through quasiperiodicity,
i.e., through period doubling bifurcations of two-dimensional (2D) tori. The
basins of attraction leading to chaos are calculated for two parameter sets of
particular importance. Furthermore, asymmetric solutions either synchronized
(leading to localized synchronization) or anti-synchronized (i.e., having
opposite phases) are investigated by calculating the maximum Lyapunov exponent
and the correlation function over the entire control parameter space.

In the next section (II), the electrical circuit equivalent model for the SQUID
dimer is described in detail, and the dynamic equations are obtained. In section
III, the bifurcation structure for the SQUID dimer is discussed along with the
emergence of quasiperiodicity and chaos. Section IV is devoted to the transition
to chaos through quasiperiodicity, through a Neimark-Sacker (NS) bifurcation. The
existence of highly asymmetric solutions related to localized synchronization is
revealed using an appropriate correlation function in section V.
Conclusions are given in section VI.

\section{Modeling Two Coupled SQUIDs}
\label{sec:model}
Consider two identical rf SQUIDs in close proximity so that they are coupled by
magnetic dipole-dipole forces through their mutual inductance $\cal M$ (Fig.
\ref{fig1}). Each of the SQUIDs is modeled by an equivalent electrical circuit
that features a self-inductance $L$ due to the superconducting ring, which is
connected in series with a ``real'' Josephson junction \cite{Josephson1962}
characterized by a critical current $I_c$, capacitance $C$, and Ohmic resistance
$R$ \cite{Hizanidis2018}. The SQUIDs are subject to an externally applied
spatially constant and time-periodic (ac) magnetic field and a constant in time
and space (dc) magnetic field. Those fields add an electromotive force in series
to the SQUID's equivalent circuits. The magnetic flux threading the loops of the
SQUIDs includes supercurrents as well as normal (i.e., quasiparticle) currents
around the SQUID's rings through Faraday's law. In turn, the induced currents
produce their own magnetics field which counter-acts the applied ones. Then, the
flux $\Phi_1$ and $\Phi_2$ through the loop of the SQUID number $1$ and $2$,
respectively, is given by the following {\em flux-balance relations}
\begin{eqnarray}
\label{2.01}
   \Phi_1 =\Phi_{ext} +L\, I_1 +{\cal M} I_2,
\\
\label{2.02}
   \Phi_2 =\Phi_{ext} +L\, I_2 +{\cal M} I_1,
\end{eqnarray}
where $I_1$ and $I_2$ are the currents induced by the external magnetic fields,
and the external flux
\begin{equation}
\label{2.03}
   \Phi_{ext} =\Phi_{dc} +\Phi_{ac}\, \cos(\omega t),
\end{equation}
with $\Phi_{dc}$ the dc flux bias, $\Phi_{ac}$ the amplitude of the ac flux,
$\omega$ its frequency, and $t$ the temporal variable. Eqs. (\ref{2.01}) and
(\ref{2.02}) can be solved for the currents and written in matrix form as
\begin{eqnarray}
\label{2.04}
  \begin{bmatrix}
    I_1 \\
    I_2
  \end{bmatrix}
  =\frac{1}{L}
\begin{bmatrix}
    1       & \lambda \\
    \lambda & 1
  \end{bmatrix}^{-1}
  \begin{bmatrix}
    \Phi_1 -\Phi_{ext} \\
    \Phi_2 -\Phi_{ext}
  \end{bmatrix},
\end{eqnarray}
where $\lambda=\frac{\cal M}{L}$ is the magnetic coupling strength between the
SQUIDs. Note that the sign of $\lambda$ depends on the mutual position of the two
SQUIDs. In the planar geometry (Fig. \ref{fig1}(a)), that sign is negative,
while in the axial geometry (Fig. \ref{fig1}(b)), it is positive. The currents
flowing in the SQUIDs are given within the framework of the resistively and
capacitively shunted junction (RCSJ) model \cite{Likharev1986}, as
\begin{equation}
\label{2.05}
    I_n =-C \frac{d^2 \Phi_n}{dt^2} -\frac{1}{R} \frac{d \Phi_n}{dt}
     -I_c\, \sin\left( 2 \pi \frac{\Phi_n}{\Phi_0} \right),
\end{equation}
where $\Phi_0$ is the flux quantum, and $n=1,2$. Combining Eqs. (\ref{2.04}) and
(\ref{2.05}) we get the coupled equations for the fluxes through the loops of
the SQUIDs in natural units, as
\begin{eqnarray}
\label{2.06}
   L\, C \frac{d^2 \Phi_1}{dt^2} +\frac{L}{R} \frac{d \Phi_1}{dt}
   +L\, I_c\, \sin\left( 2 \pi \frac{\Phi_1}{\Phi_0} \right)
   +\frac{1}{1-\lambda^2} \Phi_1
\nonumber \\
   -\frac{\lambda}{1-\lambda^2} \Phi_2 =\frac{1}{1+\lambda} \Phi_{ext}, \\
\label{2.07}
   L\, C \frac{d^2 \Phi_2}{dt^2} +\frac{L}{R} \frac{d \Phi_2}{dt}
   +L\, I_c\, \sin\left( 2 \pi \frac{\Phi_2}{\Phi_0} \right)
   +\frac{1}{1-\lambda^2} \Phi_2
\nonumber \\
   -\frac{\lambda}{1-\lambda^2} \Phi_2 =\frac{1}{1+\lambda} \Phi_{ext},
\end{eqnarray}
In normalized form, Eqs. (\ref{2.06}) and (\ref{2.07}) can be written as
\begin{eqnarray}
\label{2.08}
   \ddot{\phi_1} +\gamma \dot{\phi_1}
   +\beta\, \sin\left( 2 \pi \phi_1 \right)
   +\frac{1}{1-\lambda^2} \phi_1
   -\frac{\lambda}{1-\lambda^2} \phi_2
\nonumber \\
   =\frac{1}{1+\lambda} \phi_{ext},
\\
\label{2.09}
   \ddot{\phi_2} +\gamma \dot{\phi_2}
   +\beta\, \sin\left( 2 \pi \phi_2 \right)
   +\frac{1}{1-\lambda^2} \phi_2
   -\frac{\lambda}{1-\lambda^2} \phi_1
\nonumber \\
   =\frac{1}{1+\lambda} \phi_{ext},
\end{eqnarray}
with
\begin{equation}
\label{2.10}
   \phi_{ext} =\phi_{dc} +\phi_{ac}\, \cos\left( \Omega \tau \right),
\end{equation}
where all the fluxes $\Phi_1, \Phi_2, \Phi_{ac}, \Phi_{dc}$ are in units of
the flux quantum $\Phi_0$ and the new temporal variable $\tau$ is defined as
$\tau =t\, \omega_{LC}$ with $\omega_{LC} =1/\sqrt{L\, C}$ being the
inductive-capacitive ($L\, C$) SQUID frequency. Note that the normalized driving
frequency $\Omega$ is in units of $\omega_{LC}$, i.e.,
$\Omega =\omega/\omega_{LC}$. The overdots in Eqs. (\ref{2.08}) and (\ref{2.09})
denote differentiation with respect to the normalized temporal variable $\tau$.
The rescaled SQUID parameter and the loss coefficient are given respectively by
\begin{equation}
\label{2.11}
   \beta =\frac{L\, I_c}{\Phi_0}, \qquad \gamma =\frac{1}{R} \sqrt{ \frac{L}{C} }.
\end{equation}

\begin{figure}[!t]
   \includegraphics[scale=0.22]{Fig2.jpg}
   \caption{(a) The potential $U_{SQ}(\phi_1,\phi_2)$ of the SQUID dimer for
            $\tau =T =2\pi/\Omega$ where $\Omega =0.891$.
            (b) The corresponding one-dimensional projection of
            $U_{SQ}(\phi_1,\phi_2)$ on the $\phi_2 =0$ plane.
            Other parameters: $\lambda=0.31$, $\phi_{ac}=0.02$, $\gamma=0.024$,
            and $\beta=0.1369$.}
\label{fig2}
\end{figure}

In the lossless case, i.e., when $\gamma=0$ ($R \rightarrow \infty$), the
(Newtonian) dynamics of the SQUID dimer is governed by the equations
\begin{equation}
\label{2.12}
   \ddot{\phi}_1 =-\frac{\partial U_{SQ} (\phi_1, \phi_2)}{\partial \phi_1},
\qquad
   \ddot{\phi}_2 =-\frac{\partial U_{SQ} (\phi_1, \phi_2)}{\partial \phi_2}
\end{equation}
in terms of the SQUID dimer potential
\begin{eqnarray}
\label{2.13}
   U_{SQ} (\phi_1,\phi_2) &=&-\frac{1}{1 +\lambda}
                         \left( \phi_1 +\phi_2 \right) \phi_{ext}(\tau)
\nonumber \\
                        & &-\frac{\beta}{2\pi}
        \left[ \cos( 2 \pi \phi_1 ) +\cos( 2 \pi \phi_2 ) \right]
\nonumber \\
                        & &+\frac{1}{2 ( 1 - \lambda^2 )}
        \left( \phi_1^2 -2\lambda \phi_1 \phi_2 +\phi_2^2 \right).
\end{eqnarray}

\begin{figure*}[!t]
   \includegraphics[scale=0.38]{Fig3.jpg}
   \caption{(a) The resonance curve of the SQUID dimer. Black lines correspond
            to stable branches and red and blue lines correspond to unstable
            ones. The latter are born through saddle-node bifurcations of limit
            cycles (SN), Neimark-Sacker bifurcations (NS), or pitchfork
            bifurcations of limit cycles (PB).
            (b) Maximum values of the magnetic flux $\phi_1$ (top) and all five
            Lyapunov exponents (bottom), as a function of the driving frequency
            in the interval corresponding to the left blue branches of (a).
            For each value of $\Omega$, $n_i =50$ different initial conditions
            were used.
            (c) The same as in (b) for the right blue branches of (a).
            Other parameters: $\lambda=0.31$, $\phi_{ac}=0.02$, $\gamma=0.024$,
            and $\beta=0.1369$.}
\label{fig3}
\end{figure*}

Since the SQUID dimer considered here consists of two identical non-hysteretic
SQUIDs ($2\pi\beta<1$), the potential $U_{SQ} (\phi_{1},\phi_{2})$ is a
two-dimensional corrugated parabola with a single minimum that depends
periodically on time $\tau$. A snapshot of the potential at $\tau =T$, with
$T =2\pi/\Omega$ being the period of the (driving) external periodic field, is
shown in Fig.~\ref{fig2}(a) while its one-dimensional projection on the
$\phi_2 =0$ plane is shown in Fig.~\ref{fig2}(b).

In what follows, the external dc flux is set to zero, i.e., $\phi_{dc} =0$.
For obtaining the values of the SQUID-dependent normalized parameters $\beta$
and $\gamma$ that go into the model equations (\ref{2.08}) and (\ref{2.09}), we
choose $L=120 ~pH$, $C=1.1 ~pF$, $R=500 ~\Omega$, and $I_c =2.35 ~\mu A$ for
their equivalent lumped-circuit elements. These values are typical for
non-hysteretic SQUIDs investigated recently in the context of SQUID
metamaterials \cite{Trepanier2013,Zhang2015}. By substituting them into  Eq.
(\ref{2.11}), we get $\beta_L =0.86$ ($\beta \simeq 0.1369$) and $\gamma =0.024$,
while the $L C$ or {\em geometrical} frequency is
$f_{LC} =\omega_{LC} /(2\pi) \simeq 13.9 ~GHz$. As mentioned earlier, the value
of the coupling strength $\lambda$  depends significantly on the relative
positions of the two SQUIDs, and may assume its value from the relatively broad
range of $[-0.06, +0.50]$ that corresponds to technologically feasible SQUID
dimer designs \cite{Anlage0000}. Besides, there are also the parameters
$\phi_{ac}$ and $\Omega$ ($\phi_{dc} =0$), which can be tuned externally. Here,
however, we prefer to fix the former to $\phi_{ac} =0.02$, which is well into
the range of the experimentally accessible values, and use $\lambda$ and $\Omega$
as control parameters.

\section{Bifurcations, Quasi-periodicity, and Chaos}
It is known that for relatively high $\phi_{ac}$, the single SQUID exhibits a
multistable response which is reflected in its corresponding ``snake-like''
resonance curve~\cite{Hizanidis2018}, i.e., the curve formed by the maximum flux
within one driving period $T=2\pi/\Omega$, $\phi_{max}$, as a function of the
driving frequency $\Omega$. For the value of $\phi_{ac}=0.02$ considered here,
the single SQUID is bistable. When two SQUIDs are coupled positively (i.e., when
$\lambda > 0$) the bistability is maintained but the dynamics presents additional
complexity. The resonance curve for the SQUID dimer is plotted in
Fig.~\ref{fig3}(a) in terms of the magnetic flux $\phi_1$ threading the loop of
the SQUID number $1$. Specifically, the maximum flux threading the loop of the
SQUID number $1$ within one driving period $T$, $\phi_{1,max}$, is plotted
against the driving frequency $\Omega$.
We see the typical tilted curve associated with bistability and hysteresis,
where two stable solutions (black curves) coexist with an unstable one (red
curve) below the geometric resonance frequency $\Omega =1$. Starting at lower driving frequencies and following the periodic
solution in $\Omega$ there is a saddle-node bifurcation of limit cycles (SN)
occurring at $\Omega=1.2$. The solution then becomes unstable and when it
reaches the peak of the resonance curve it turns stable again in a second SN
bifurcation at $\Omega=0.8036$. The stable branch that is born then undergoes a
pitchfork bifurcation of limit cycles (PB) at $\Omega=0.8684$, where in addition
to the unstable branch two more stable branches of periodic solutions are born.
These branches consecutively become unstable in a Neimark-Sacker bifurcation (NS)
on either side at $\Omega=0.87$, shown in blue color. The blue branches become stable again in a
second NS bifurcation at $\Omega=0.8944$, again on either side. The same scenario takes place at
$\Omega=1.12$ (right PB) and this family of solutions emanating from the two
``opposite to each other'' PB bifurcations, undergoes NS bifurcations (at $\Omega=1.083$ and $\Omega=1.114$) on
either side forming, thus, a closed loop.
The bifurcation lines have been obtained
using a very powerful software tool that executes a root-finding algorithm for
continuation of periodic solutions~\cite{ENG02}.

We focus now on the two $\Omega$ intervals corresponding to the blue unstable
branches of Fig.~\ref{fig3}(a) which lie within the two PB bifurcations on both right and left side.
In the upper panels of Figs.~\ref{fig3}(b)
and (c) the maxima of the solution for $\phi_1 (\tau)$ are plotted over $\Omega$
for the left and right interval, respectively. Note that at each frequency
$\Omega$, a large number of different initial conditions $n_i$ were used (here
$n_i =50$) in order to obtain all the coexisting solutions within the $\Omega$
interval shown, and also to obtain a good representation of possible chaotic
behavior. Apparently, chaotic behavior appears in both the upper panels of
Figs.~\ref{fig3}(b) and (c). These solution branches coexist with the blue
unstable branches of Fig.~\ref{fig3}(a) mentioned above; the reason for which
they are not superposed to them is merely clarity.

In the corresponding lower
panels, the Lyapunov exponents of the obtained solutions that are plotted for
the same intervals of $\Omega$, reveal chaotic behavior in the regions where
the maximum Lyapunov exponent $L_1$ is positive (positive Lyapunov exponents are
shown as blue points). Careful inspection of the Lyapunov exponents indicate
that transitions to chaos through quasiperiodicity and reversely take place in
Figs.~\ref{fig3}(b) and (c) (see next section). The Lyapunov exponents were calculated employing the algorithm from \cite{Geist1990} using
the Julia$^R$ software library. Similar multibranched resonance curve are
sometimes encountered in dissipative-driven oscillators
\cite{Warminski2015,Marchionne2018,Zang2019} where in addition, isolated
branches (``isolas'') may exist. However, to the best of our knowledge,
multibranched resonance curves whose branches are connected through
Neimark-Sacker bifurcations have never before been observed.

Since the system dynamics is very sensitive to initial conditions, we have
identified the basins of attraction of chaotic and periodic behavior for two
$\Omega$ values within the intervals shown in Figs.~\ref{fig3}(b) and (c)
close to the NS bifurcations of limit cycles.
These are displayed in Figs.~\ref{fig4}(a) and (b), respectively, both in the
$(\phi_1,\phi_2)$ plane (left panels) and $(\dot{\phi}_1, \dot{\phi}_2)$ plane
(right panels). Blue (dark) regions denote the set of initial conditions leading
to chaotic motion, while light yellow (light) regions the ones leading to
periodic or quasiperiodic oscillations. The other state variables are initially
fixed as $\dot{\phi}_1 (\tau=0) =\dot{\phi}_1 (\tau=0) =0$ and
${\phi}_1 (\tau=0) ={\phi}_1 (\tau=0) =0$, respectively. In Fig.~\ref{fig4}(a),
an interesting intertwining between the basins of attraction of chaotic and
periodic or quasiperiodic states can be seen. This is apparent from the presence
of yellow spots located erratically within the blue regions in both subfigures,
that resemble riddled basins investigated in various low-dimensional systems
\cite{Lai1994,Sathiyadevi2019}. The inset of Fig.~\ref{fig4}(a) reveals this intricate structure of the basin of attraction on a smaller scale. Note that a classification of basins of
attraction in several low-dimensional systems has been perormed in Ref.
\cite{Sprott2015}. The knowledge of the basins is essential for determining the
usefulness of the systems in practical applications.

The other control parameter in our system is the coupling strength $\lambda$,
which (as described in Section~\ref{sec:model}) can obtain both negative and
positive values within a certain range provided by calculations on feasible
SQUID dimer designs \cite{Anlage0000}. In Fig.~\ref{fig6}(a) the co-dimension 2
bifurcation diagram is shown in the $(\Omega,\lambda)$ parameter plane. Black
lines mark the continuation of the saddle-node bifurcations and red lines that
of the Neimark-Sacker (torus) bifurcations. In Fig.~\ref{fig6}(b) the maximum Lyapunov
exponent is plotted in the $(\Omega,\lambda)$ parameter space. Blue and yellow
regions correspond to chaotic (positive $L_{\text{max}}$) and non-chaotic
(negative or zero $L_{\text{max}}$) dynamics, respectively. We can see that the
chaotic regimes lie within the boundaries of the NS bifurcation curves. In the
following section we explore in detail the transition from quasiperiodicity to
chaos in our system.

\begin{figure}
   \includegraphics[scale=0.4]{Fig4.jpg}
   \caption{Basins of attraction for the SQUID dimer in terms of the maximum
            Lyapunov exponent $L_{\text{max}}$. Blue and yellow color
            corresponds to chaotic and non-chaotic dynamics, respectively.
	    Left: $\phi_1$ versus $\phi_2$ for
            $(\dot{\phi}_1(0),\dot{\phi_2}(0)) = (0,0)$.
            Right: $\dot{\phi}_1$ versus $\dot{\phi}_2$ for
            $(\phi_1(0),\phi_2(0)) = (0,0)$.
            (a) $\Omega=0.891$, and (b) $\Omega=1.1$. The inset in the left panel of (a)
            shows a blow-up of the region inside the red box. Other parameters:
            $\lambda=0.31$, $\phi_{ac}=0.02$, $\gamma=0.024$, and $\beta=0.1369$.}
\label{fig4}
\end{figure}

\section{Torus-doubling transition to chaos}
From Figs.~\ref{fig3}(b)-(c) and ~\ref{fig5}(b) it is clear that the SQUID dimer
described by Eqs. (\ref{2.08}) and (\ref{2.09}) becomes chaotic through
Neimark-Sacker bifurcations. This quasiperiodic transition
to chaos is a well-known scenario whereby a torus in low-dimensional dynamical
systems loses its stability and develops into chaos. A path from a limit cycle,
then a transition into 2D torus following by a 3D torus and eventually direct
transition to chaos has been presented numerically and experimentally is an
analog electrical circuit representing the ring of unidirectionally coupled
single-well Duffing oscillators \cite{Borkowski2015}. Another possible scenario
to chaos is a cascade of period doublings of the torus which has been
investigated numerically long ago using low-dimensional maps
\cite{Arneodo1983,Kaneko1984}, and experimentally in electrochemical reactions
\cite{Bassett1989}. According to the second scenario, a finite cascade of period
doubling of such invariant tori leads eventually to chaos. This scenario has
been reported numerically for a quintic complex Ginzburg-Landau equation
\cite{Kim1997} and bimodal laser model \cite{Letellier2007} as well as
experimentally, in Rayleigh-Benard convection \cite{Flesselles1994}, in a simple
thermoacoustic system \cite{Mondal2017}, and near the ferroelectric phase
transition of $KH_{2}PO_{4}$ crystals \cite{Shin1999}.

\begin{figure}[!t]
   \includegraphics[scale=0.22]{Fig5.jpg}
   \caption{(a) Bifurcation diagram in the ($\Omega,\lambda$) parameter plane of
            the two coupled SQUIDs. Black and red lines correspond to saddle-node
            (SN) bifurcations of limit cycles and Neimark-Sacker (NS)
            bifurcations, respectively.
            (b) Map of the maximum Lyapunov exponent $L_{\text{max}}$ in the
            ($\Omega, \lambda$) parameter space, on which the NS bifurcation
            curve from (a) is superposed (red lines). Blue and yellow color
            corresponds to chaotic and non-chaotic dynamics, respectively.
            Other parameters:
            $\phi_{ac}=0.02$, $\gamma=0.024$, and $\beta=0.1369$.}
\label{fig5}
\end{figure}

\begin{figure}[!b]
   \includegraphics[scale=0.24]{Fig6.jpg}
   \caption{(a) The resonance curve of the SQUID dimer for $\lambda=0.15$. Black lines correspond
            to stable branches and red and blue lines correspond to unstable
            ones. The latter are born through saddle-node bifurcations of limit
            cycles (SN), Neimark-Sacker bifurcations (NS), or pitchfork
            bifurcations of limit cycles (PB).
            (b) Maximum values of the magnetic flux $\phi_1$ of stable solutions as a function of the driving frequency
            in the same $\Omega$ interval.
            Other parameters: $\phi_{ac}=0.02$, $\gamma=0.024$,
            and $\beta=0.1369$.
            }
\label{fig6}
\end{figure}

 \begin{figure}[!b]
   \includegraphics[scale=0.22]{Fig7.jpg}
   \caption{The complete set of Lyapunov exponents $L_i$, $i=1,...,5$ for the SQUID dimer (top) and the corresponding Lyapunov dimension (bottom), as a function of the driving frequency $\Omega$. The
            bifurcation analysis depicts a transition from a stable
            two-dimensional torus to chaos via a sequence of four consecutive
            period-doubling bifurcations marked by $P_i$ ($i=2, 4, 8, 16$), where
            the period of the solution doubles.
            Other parameters:
            $\lambda=0.15$, $\phi_{ac}=0.02$, $\gamma=0.024$ and $\beta=0.1369$.
            }
\label{fig7}
\end{figure}

We examine this route to chaos through period doubling of a 2D torus in our system
by fixing the magnetic coupling strength to a certain value $\lambda=0.15$
(marked by a black horizontal line in Fig.~\ref{fig5}(b)) and by varying the
driving frequency $\Omega$ in a narrow interval that encloses chaotic states.

For this choice of parameters, the resonance curve (Fig.~\ref{fig6}(a)) exhibits the same sequence of bifurcations (PB and consequently NS) as in Fig.~\ref{fig3}(a), but presents no second chaotic band and therefore no ``bubble'' connected to it. The transition into chaos is shown in Fig.~\ref{fig6}(b) where the maximum values of the magnetic flux $\phi_1$ are plotted for the same $\Omega$ interval.

For each value of $\Omega$ the whole spectrum of the Lyapunov exponents is
calculated and plotted in Fig.~\ref{fig7} (top). Starting from the right end of the
frequency interval, it is observed that the two largest Lyapunov exponents are
zero ($L_1=L_2=0$), while the other three are less than zero. This indicates
quasiperiodic dynamics on a stable 2D torus.
When the
frequency $\Omega$ decreases, one more Lyapunov exponent reaches zero at
$\Omega=1.27032$ (marked as P2 on the figure), indicating a period doubling
bifurcation of the torus. By further decreasing $\Omega$, we find three more
points in which the three largest Lyapunov exponents are zero at
$\Omega=1.269505$, $1.26933$, and $1.269295$ indicating three more period
doupling bifurcations of the 2D torus (marked by P4, P8, and P16, respectively,
on the figure). For frequencies even lower than that at P16 the 2D torus
break down and the system enters into a chaotic state, as it is indicated by the
positivity of the largest Lyapunov exponent ($L_1 =L_{max} >0$).
This transition to chaos is also reflected in the
Lyapunov dimension $d_L$, which according to Kaplan and York~\cite{KAP79} is given by:
\begin{equation}
\label{01}
  d_L = k +\frac{1}{|L_{k+1}|} \sum_{i=1}^k L_i,
\end{equation}
where $k$ is defined by the condition that
\begin{equation}
\label{02}
  \sum_{i=1}^k L_i \geq 0 ~~~{\rm and} ~~~ \sum_{i=1}^{k+1} L_i < 0.
\end{equation}
As the 2D torus transitions into chaos, the Lyapunov dimension changes from the value $d_L=2$ to $d_L>3$, as shown in Fig.~\ref{fig7} (bottom).

\begin{figure*}[t!]
   \includegraphics[scale=0.45]{Fig8.jpg}
   \caption{Poincar{\'e} cross sections on the $\Phi-\dot{\Phi}$ plane (upper
            panels) and the corresponding frequency spectra plotted in
            a logarithmic scale (bottom panels), for (a) $\Omega=1.2705$;
            (b) $\Omega=1.27$; (c) $\Omega=1.2694$;(d) $\Omega=1.26931$;
            (e) $\Omega=1.26923$. In the upper panels, blue and red lines
            correspond to the first ($\Phi =\phi_1$, $\dot{\Phi} =\dot{\phi}_1$)
            and second ($\Phi =\phi_2$, $\dot{\Phi} =\dot{\phi}_2$) SQUID,
            respectively.
            {\em Insets:} Enlargement of a small part of Poincar{\'e} sections
            of the trajectories of the first SQUID to make clear the period of
            the tori. Other parameters as in Fig. \ref{fig7}.
\label{fig8}}
\end{figure*}

Typical Poincar{\'e} sections of the period-2, -4, -8, and chaotic 2D tori
on the $\phi_{1,2} - \dot{\phi}_{1,2}$ phase plane are shown in the upper panels
of Figs.~\ref{fig8}(a) - (e). Blue and red lines correspond to flux threading
the loop of the first ($1$) and the second ($2$) SQUID, respectively. The insets
show enlargements of a small part of the Poincar{\'e} sections of the first
SQUID (SQUID number $1$), in order to see clearly the varying number of layers
of the tori as the frequency $\Omega$ decreases.
The Poincar{\'e} section of the period-16 torus is very similar to that of the
period-8 torus and thus it has been omitted. In the lower panels of
Fig.~\ref{fig8}, the corresponding Fourier power spectral densities in a
logarithmic scale are shown.
The 2D torus (upper panel of Fig.~\ref{fig8}(a)) has a dominant frequency which
in the present case coincides with the external (driving) frequency $\Omega$,
as well as four more frequencies denoted as $\Omega_{\pm 1}$, $\Omega_{\pm 2}$
(lower panel of Fig.~\ref{fig8}(a)). The exact values of these frequencies are
obtained from the locations of the sharp peaks. Then, as the period of the 2D
torus increases to $2$, $4$, and $8$ through period-doubling bifurcations, a
doubling of the sharp peaks occurs (as can be seen in Fig.~\ref{fig8}(b), (c),
and (d), respectively) that correspond to frequencies with significant spectral
content. In Fig.~\ref{fig8}(e), in which the SQUID dimer is in a chaotic state,
the power spectrum exhibits a substantial noisy background, along with several
sharp peaks that correspond to frequencies with significant spectral content.
Note that the driving frequency $\Omega$ is dominant in all five subfigures.

Thus, we have identified a novel route to chaos for the SQUID dimer, which
initiates from a single limit cycle. The latter then bifurcates into a 2D
torus, which follows a finite period doubling cascade which eventually leads to
a chaotic state. Note that in the case of two coupled Duffing oscillators with
dissipation and driving force, three period-doubling bifurcations of a
three-torus were identified for the transition to chaos through quasiperiodicity
\cite{Kozlowski1995}.

\section{Localized $\&$ Chaotic Synchronization}
The SQUID dimer considered here consists of two symmetrically coupled identical
nonlinear oscillators. Such systems were once believed to support only totally
synchronous or asynchronous states. However, as it is known today, there are more
possible forms of synchronization such as the so called localized
synchronization. This form of synchronization has been originally investigated
in a slightly dissimilar pair of solid state \cite{Kuske1997} and semiconductor
\cite{Hohl1997} lasers. For two coupled oscillators, localized synchronization
means that one of the oscillators exhibits strong oscillations while the other
one exhibits weak oscillations. In the other limit, that of asynchronous dynamics,
extreme case of having the one oscillator exhibiting periodic motion and the other
chaotic motion has been recently reported \cite{Awal2019}. We explore below the
effect of localized synchronization for the SQUID dimer using the correlation
function between the time-dependent fluxes through the loop of the first and the
second SQUID $\phi_{1}(t)$ and $\phi_{2}(t)$, respectively, i.e.,
\begin{equation}
\label{999}
   C_{12} =\frac{\langle (\phi_{1}(t)-\mu_{1})(\phi_{2}(t)-\mu_{2}) \rangle}
                {\sigma_{1} \sigma_{2}},
\end{equation}
where $\mu_{1,2}$ and $\sigma_{1,2}$ are the mean value and the standard
deviation of the corresponding time series of $\phi_1$ and $\phi_2$,
respectively. Using this measure we can quantify the asymmetry in the amplitudes
of the two coupled SQUID magnetic fluxes. If $\phi_{1}$ and $\phi_{2}$ are
increasing or decreasing simultaneously then the correlation function will be
positive. In contrast, if the value of $\phi_{1}$ is increasing while $\phi_{2}$
is decreasing and vice versa, then $C_{12}$ is negative. The maximum/minimum
value of the correlation function is $\pm 1$. For $C_{12}\simeq 1$ the time
series of $\phi_1$ and $\phi_2$ are almost identical thereby the two SQUIDs
exhibit in-phase synchronization, whereas for $C_{12}\simeq -1$ they are in
anti-phase synchronization.

In Fig.~\ref{fig9}(a), the asymptotic value of the correlation function $C_{12}$
is mapped onto the $\Omega, \lambda$ parameter space. As mentioned earlier, the
value of $C_{12}$ reflects the dynamical behavior of the SQUID dimer as long as
its synchronization properties are concerned. As it can be observed, the value
of $C_{12}$ is equal to $1$ in a large area (yellow) of the parameter space but
there also exist regions where $C_{12} \neq 1$. For example, a dark-blue region
(\Romannum{1}) can be observed in which $C_{12}$ approaches $-1$, indicating
anti-synchronized (i.e., phase difference $\pi$ between $\phi_1 (t)$ and
$\phi_2 (t)$) dynamics. Moreover, we observe a ``horn''-shaped region of
intermediate to high values of $C_{12}$ which coincides with the corresponding
``horn''-shaped chaotic region of Fig.~\ref{fig6}(b). This particular region,
which is bounded by the Neimark-Sacker (NS) bifurcation lines, intermittent
chaotic synchronization \cite{Baker1998,Zhao2005} is observed. In the
case of the SQUID dimer, the effect of intermittent chaotic synchronization is
manifested by the entrainment of the two SQUID time series for $\phi_1(t)$ and
$\phi_2(t)$ in random temporal intervals of finite duration~\cite{Zhao2005}.
This will be further investigated in future works involving
SQUID trimers where more interesting synchronization phenomena can be observed,
like small chimeras~\cite{Banerjee2018}.

Next we focus on the behavior of the correlation function $C_{12}$ in the
non-chaotic region of the $(\Omega,\lambda)$ parameter space and reveal another
interesting behavior. For a fixed $\Omega$ we study the dependence of $C_{12}$
as the coupling strength $\lambda$ varies from region \Romannum{1} to
\Romannum{2} of Fig.~\ref{fig9}(a). A cross section of $C_{12}$ for
$\Omega =0.9822$ is plotted in Fig.~\ref{fig9}(b) in red color, for one random
set of initial conditions (top) and $n_i =30$ sets of different, randomly chosen
initial conditions (bottom). The difference between the maximum values of
the two SQUID magnetic fluxes $\phi_{1\text{max}}-\phi_{2\text{max}}$ is also
plotted in black color. We make the following observations: first of all, the
fully phase-synchronized state ($C_{12}=1$) exists for all $\lambda$ in this
interval.
Secondly, the correlation function appears to have a (co)sinusoidal dependence
on the coupling strength $\lambda$ as expected, since the correlation function
of two out-of-phase cosine signals is a cosine function of their phase
difference \cite{Oppenheim1996}.
In the low-coupling non-chaotic regions \Romannum{1} and \Romannum{2}
we can assume that the two SQUID solutions are also (co)sine-like. For
$\phi_{1}(t)$ and $\phi_{2}(t+\tau)$ or vice versa for $\phi_{1}(t+\tau)$ and
$\phi_{2}(t)$, $\tau$ is the phase difference between the SQUID fluxes and is
dependent on the set of initial conditions. From Fig.~\ref{fig9}(b) (bottom) we
may claim that the $\tau=0$ case corresponds to $C_{12}=1$ and the
(co)sinusoidal-like form of $C_{12}$ corresponds to the case where $\tau$ is a
function of $\lambda)$. In the former case the two SQUIDs are in phase and
$\phi_{1\text{max}}-\phi_{2\text{max}}$ is zero, while in the latter case their
phase difference varies from values close to $-1$ (almost anti-phase) to values
close to $+1$ (almost in-phase). The difference
$\phi_{1\text{max}}-\phi_{2\text{max}}$ attains a nonzero positive/negative
value depending on whether $\phi_1$ is larger/smaller than $\phi_2$, for all
$\lambda$ and all initial conditions. Therefore this case corresponds to an
asymmetric solution as far as the amplitudes of the two coupled SQUIDs are
concerned.

This is illustrated in Fig.~\ref{fig9}(c) in the phase portraits (top) and time
series (bottom) of the magnetic flux for the first SQUID (green line) and the
second SQUID (blue line) for parameters tha belong to region \Romannum{1}
$(\lambda=-0.0398, C_{12}=-0.9103)$ and region \Romannum{2}
$(\lambda=0.0672, C_{12}=0.9248)$. These plots show the phase portraits of the
corresponding stable periodic solutions where large asymmetry in the amplitudes of
the two coupled SQUIDs occurs. Due to that difference between the amplitudes,
these two states are illustrative examples of localized synchronization.
Moreover, through the time series, it can be clearly seen that the fluxes
$(\phi_{1}$  and $\phi_{2})$ are almost out of phase (I) and almost in phase
(II) because the correlation function has values close to $-1$ and $1$,
respectively.

\begin{figure*}[!t]
   \includegraphics[scale=0.4]{Fig9.jpg}
   \caption{(a) Dynamical behavior in the ($\Omega, \lambda$) parametric space
            of the correlation function $C_{12}$. For a completely in-phase
            synchronized state $C_{12} =1$ while for a completely anti-phase
            synchronized state $C_{12} =−1$.
            (b) The correlation function (red line) and maxima differences
            between $\phi_{1}$ and $\phi_{2}$ (black line), versus the coupling
            strength ($\lambda$), for one random set of initial conditions (top)
            and for $n_i =30$ different initial conditions (bottom) where
            $\Omega=0.9822$.
            (c) Phase portraits (top) and time series (bottom) of the magnetic
            flux for the first SQUID, $\phi_1$ (green line), and the second
            SQUID, $\phi_2$ (blue line), for (I) $\lambda=-0.0398$,
            $C_{12}=-0.9103$; (II) $\lambda=0.0672$, $C_{12}=0.9248$.
            Other parameters
            $\phi_{ac} =0.02$, $\gamma =0.024$, and $\beta =0.1369$.}
\label{fig9}
\end{figure*}

\section{Conclusions}
In summary, a dimer comprising two magnetically coupled SQUIDs is investigated
numerically with respect to its bifurcation structure and its transitions from
quasiperiodicity to chaos through a finite period-doubling cascade of tori.
The SQUID dimer is a prototype, dissipative-driven dynamical system of
considerably complex dynamic behavior, on which numerous nonlinear dynamic
effects can be investigated. In the present work, the bifurcation structure of
the SQUID dimer was obtained using the driving frequency $\Omega$ and the
coupling strength $\lambda$ as external control parameters that form a
two-dimensional parameter space. Bifurcation (stroboscopic) diagrams as a
function of $\Omega$ reveal multibranched resonance curves resulting through
different kinds of bifurcations, transitions to chaos after four period-doubling
bifurcations of a 2D torus, and riddle-like basins. Bifurcation diagrams for
SN and NS bifurcations are shown on the $\Omega - \lambda$ plane. In combination
with maps of the maximum Lyapunov exponent $\Lambda_{max}$, it is observed that
the areas of the $\Omega - \lambda$ plane where chaotic behavior is expected is
bounded by the NS bifurcation curves. Eventually, localized synchronization is
observed and quantified in the SQUID dimer using the correlation function
$C_{12}$. The possibility of intermittent chaotic synchronization is briefly
mentioned.

\begin{acknowledgments}
This work was supported by the General Secretariat for Research and Technology
(GSRT) and the Hellenic Foundation for Research and Innovation (HFRI) (Code: 203).
JS acknowledges support by the Ministry of Education and Science of the Russian
Federation in the framework of the Increase Competitiveness Program of NUST
``MISiS'' (Grant number K4-2018-049). The authors would like to thank
S. M. Anlage for fruitful discussions.
\end{acknowledgments}


TITLE: Is there a real continuous function of a Cube?
QUESTION [3 upvotes]: Sorry if this seems like a basic question, but I'm having troubles finding an answer in my research. I'm looking for a function that if given a test point (x, y, z) will return a real value describing whether the given point is

On the surface of the shape
$$f(x) = 0$$
Inside the volume of the shape
$$f(x) > 0$$
Or outside the volume of the shape
$$f(x) < 0$$

I've found such equations for a sphere and other shapes but cubic, and in general faceted shapes have eluded me. Unfortunately I'm a novice in this space with no formal training so any source or in depth explanations would be greatly appreciated.
Edit: As added bonus, is there an equation that handles rectangular cuboids?

REPLY [0 votes]: whist not an exact cube
$$x^n+y^n+z^n=1$$
approximates a cube of side length 2, for large n.
An image for the shape with n=40

for a general rectangular cuboids, of length=$l$, width=$w$, and height=$h$,
the equation becomes
$$\left(\frac{2x}{l}\right)^n+\left(\frac{2y}{w}\right)^n+\left(\frac{2z}{h}\right)^n=1$$
hopefully that helps someone


TITLE: What is the basis of Huygens principle?
QUESTION [8 upvotes]: When we were studying mechanical waves like sound waves, and waves on strings in class, we never studied Huygens' principle with these - and nor did we really derive the laws of reflection or refraction with them. We always in any case took them to travel linearly, and did not talk of any wavefronts then. But in studying geometric optics recently, we've begun by proving the laws of reflection and refraction by using Huygens' principle. In doing so we looked at the behaviour of wavefronts - not linear waves - which we had not with sound waves previously. In real life, I understand there must exist mostly wavefronts - there can hardly be any isolated linear wave disturbances, so reflection and refraction occur the way they do.
But if we were to talk of wavefronts of sound, will they still follow Huygens' principle? Alternately, will linear light waves behave like sound waves - for example, will they set up standing waves as well?
And is Huygens' principle any more than a simple method of analysis of how wavefronts propagate? It must just a consequence of the wave behaviour - and must be obtainable from the wave theory. Surely it isn't an independent principle?

REPLY [13 votes]: The grounding of Huygens's principle is essentially the observation that the Green's Function for the Helmholtz Wave Equation is the spherical wave source
$$\psi_G(\vec{r})=\frac{e^{i\,k\,r}}{r}$$
Since approximately monotonal sound waves also fulfill the Helmholtz equation, the reasoning below, and thus Huygens's principle, applies exactly to them as well.
A simple calculation to show that one can sum up the effect of such sources on a wavefront and get approximately the right answer runs as follows. We look at a semi-infinite region $V$ with boundary $\partial V$. Only part of that boundary - an aperture $A$ - has a significantly nonzero disturbance $\psi(\vec{r})$. We wish to find $\psi$ at some position $\vec{r}_0$ away from the aperture. We consider two functions $\psi(\vec{r})$ and the Green's function $\psi_G(\vec{r}-\vec{r}_0)$ and form the vector field $\psi(\vec{r})\,\nabla \psi_G(\vec{r}-\vec{r}_0)-\psi_G(\vec{r}-\vec{r}_0)\,\nabla \psi(\vec{r})$ and then plug this little beast into Gauss's divergence theorem for a surface comprising (1) the boundary $\partial V$, which, by assumption, for this calculation is the same as simply the aperture $A$ (because the field is by assumption small elsewhere) and (2) a small sphere of radius $\epsilon$ that excises the singularity in $\psi_G(\vec{r}-\vec{r}_0)$ at $\vec{r}_0$. So we are applying the divergence theorem to the part of the volume within $\partial V$ that is outside the little "quarantining" sphere. The divergence theorem yields:
$$\int_{\partial V}  (\psi(\vec{r})\,\nabla \psi_G(\vec{r}-\vec{r}_0)-\psi_G(\vec{r}-\vec{r}_0)\,\nabla \psi(\vec{r}))\cdot \hat{n} \,\mathrm{d} S = -4\,\pi\,\psi(\vec{r}_0) + \int_V (\psi(\vec{r})\,\nabla^2 \psi_G(\vec{r}-\vec{r}_0)-\psi_G(\vec{r}-\vec{r}_0)\,\nabla^2 \psi(\vec{r}))\,\mathrm{d} V + O(\epsilon) = -4\,\pi\,\psi(\vec{r}_0) + O(\epsilon)$$
where the $-4\,\pi\,\psi(\vec{r}_0)$ term (recall that $\psi(\vec{r}_0)$ is what we want to find) comes from the surface integral of the little quarantining sphere and the volume integral vanishes because both $\psi$ and $\psi_G(\vec{r}-\vec{r}_0)$ fulfill the Helmholtz equation outside the quarantining sphere. So you can see that what we're left with is that (after we take the limit as $\epsilon\to 0$):
$$\psi(\vec{r}_0) \propto \int_A  (\psi(\vec{r})\,\nabla \psi_G(\vec{r}-\vec{r}_0)-\psi_G(\vec{r}-\vec{r}_0)\,\nabla \psi(\vec{r}))\cdot \hat{n} \,\mathrm{d} S \approx  \int_A  \frac{\exp(i\,k|\vec{r}-\vec{r}_0|)}{|\vec{r}-\vec{r}_0|}(i\,k\,\psi(\vec{r})\,\cos(\theta(\vec{r}))-\nabla \psi(\vec{r})\cdot \hat{n}) \,\mathrm{d} A$$
which, under various approximations, leads to Huygen's principle (witness the first term, which sums up the spherical wave $\exp(i\,k|\vec{r}-\vec{r}_0|)/|\vec{r}-\vec{r}_0|$ weighted by the value of the field $\psi(\vec{r})$ on the aperture. $\cos(\theta(\vec{r}))$ is responsible for the so-called obliquity factor.

Can your explanation please by simplified to high school level? I apologise, but I cannot follow the maths.

Apologies on my part then. It's a good idea to state one's level in the question: your own question is well and thoughtfully put, so you actually come across as someone more advanced than high school. Describe your level in terms of what you know rather than age as we have several teenage physicists on this site who approach graduate level.
Anyhow, you will need to be content with the explanation that Huygens's principle was postulated by Huygens as a "guess" to explain the nature of waves. He understood the linear superposition principle: that the disturbance caused by the sum of waves is the sum of the individual disturbances and thus understood, for example, that a line of point radiators could sum together to yield a plane wavefront. In the 19th century, mathematicians came up with a rigorous description of these thoughts - the Green's function method is essentially the building up of general solutions to linear equations by linear superposition of "fundamental solutions". Lo and behold, if you work out what the "fundamental solution" for the monochromatic wave equation (Helmholtz's equation) is, it turns out to be a spherical wave point emitter, exactly as Huygens had guessed.
So Huygens's principle works for any phenomenon that is described by the Helmholtz equation. This includes light, sound and water waves (in certain approximations). The Helmholtz equation is another form of the D'Alembert wave equation that is valid when the waves are roughly monochromatic or monotonal. Even in modern quantum field theory, where the relevant equations aren't wave equations, certain path integrals are still calculated by Huygens-like methods.


TITLE: Subject GRE Exam 0568 Q.51
QUESTION [1 upvotes]: The question and its answer is given in the following pictures:

The first line in the solution is not clear for me, my questions are:
1- where is the interval from 0 to 1 ?
2- why we changed the x inside the floor function to n and did not change the x in the power of e ?
Could anyone illustrate this for me please?

REPLY [3 votes]: Our task here is nothing but calculating Laplace Transform of $\lfloor t\rfloor$ at $s=1$

$$\begin{align}F(s)&=\displaystyle\int_0^\infty \lfloor t \rfloor e^{-st}dt\\&=\displaystyle\int_0^1 \lfloor t \rfloor e^{-st}dt +\displaystyle\int_1^2 \lfloor t \rfloor e^{-st}dt +\displaystyle\int_2^3 \lfloor t \rfloor e^{-st}dt + \cdots\\&=\displaystyle\sum_{n=0}^\infty n\displaystyle\int_{n}^{n+1}  e^{-st}dt\\&=\dfrac{1}{s}\displaystyle\sum_{n=0}^\infty n \left[ e^{-sn}-e^{-s(n+1)} \right]\\&=\dfrac{1}{s}\left(1-e^{-s}\right)\displaystyle\sum_{n=0}^\infty n  e^{-sn}\\&=\dfrac{1}{s}\left(1-e^{-s}\right)\dfrac{e^{-s}}{\left(1-e^{-s}\right)^{2}}\\&=\dfrac{1}{s\left(e^s-1\right)}\\&=\dfrac{\coth\left(\dfrac{s}{2}\right)-1}{2s}\end{align}\tag*{}$$

Now just put $s=1$
Note that:
$$\displaystyle\sum_{n=0}^{\infty} e^{-n x} = \dfrac{1}{1-e^{-x}}$$
$$\displaystyle\sum_{n=0}^{\infty}n e^{-n x}=- \dfrac{d}{dx}\left(\dfrac{1}{1-e^{-x}}\right)$$