\section{Introduction} \label{section-introduction}
\input{section/introduction}

\section{Gradient Descent} \label{section-gradient-descent}
\input{section/gradient-descent}

\section{Global Convergence Analysis without Summability} \label{section-no-summable}
\input{section/global-convergence-no-sum}

\section{Global Convergence Analysis with Summability} \label{section-summable}
\input{section/global-convergence-sum}

\section{An Example for Divergence} \label{section-divergence-example}
\input{section/divergence-example}

\section{Spuriousness of the Continuous Approximation} \label{section-continuous-approximation}
\input{section/continuous-approximation}

\section{Conclusion}
\label{section-conclusion}
\input{section/conclusion}

\bibliographystyle{siamplain}

\section{A detailed example}

Here we include some equations and theorem-like environments to show
how these are labeled in a supplement and can be referenced from the
main text.
Consider the following equation:
\begin{equation}
  \label{eq:suppa}
  a^2 + b^2 = c^2.
\end{equation}
You can also reference equations such as \cref{eq:matrices,eq:bb}
from the main article in this supplement.

\lipsum[100-101]

\begin{theorem}
An example theorem.
\end{theorem}

\lipsum[102]

\begin{lemma}
An example lemma.
\end{lemma}

\lipsum[103-105]

Here is an example citation: \cite{KoMa14}.

\section[Proof of Thm]{Proof of \cref{thm:bigthm}}
\label{sec:proof}

\lipsum[106-112]

\section{Additional experimental results}
\Cref{tab:foo} shows additional
supporting evidence.

\begin{table}[htbp]
\footnotesize
  \caption{Example table.}  \label{tab:smfoo}
\begin{center}
  \begin{tabular}{|c|c|c|} \hline
   Species & \bf Mean & \bf Std.~Dev. \\ \hline
    1 & 3.4 & 1.2 \\
    2 & 5.4 & 0.6 \\ \hline
  \end{tabular}
\end{center}
\end{table}

\bibliographystyle{siamplain}

\subsection{A Recursive Relationship for the Optimality Gap} \label{subsection-recursion-gap}

In order to establish a recursive relationship for the optimality gap between two sequential iterates, we study the behavior of the iterates prior to a stopping time. To define this stopping time, first let
\begin{equation}
\chi_{k}^z (R) = \begin{cases}
1 & \inlinenorm{ x_j - z }_2 \leq R,~ \forall j \leq k \\
0 & \text{otherwise},
\end{cases}
\end{equation}
where $k \in \mathbb{N}$; $R > 0$; and $z \in \mathbb{R}^p$. That is, $\chi_k^z(R)$ is one if all iterates up to $k$ remain with a radius $R$ of the point $z \in \mathbb{R}^p$, and is zero thereafter. Thus, $\lbrace \chi_k^z(R) \rbrace$ implicitly define a stopping time: the iterate at which $\chi_k^z(R)$ stops being $1$ (note, this iterate could be $k=0$).

To understand how $\chi_k^{z}(R)$ will be used, consider what happens if $\chi_k^z(R) = 1$ for all $k +1 \in \mathbb{N}$, for some $z\in\mathbb{R}^p$, and for some $R > 0$. Then, \cref{assumption-gradient} would imply that there is a universal Lipschitz constant for the gradient and we can proceed with the standard analysis to show that the gradient function goes to zero~\cite{zoutendijk1960}. However, owing to the loss of monotonic decay of the objective function and boundedness of level sets, we cannot guarantee that $\chi_k^z(R) = 1$ for all $ k + 1 \in \mathbb{N}$ and any $R > 0$ or any $z \in \mathbb{R}^p$; thus, we have no simple way of controlling the Lipschitz constant. To account for this, rather than establishing a relationship between $F(x_{k+1}) - F_{l.b.}$ and $F(x_k) - F_{l.b.}$, we establish a relationship between $[F(x_{k+1}) - F_{l.b.}]\chi_{k+1}^0(R)$ and $[F(x_k) - F_{l.b.}]\chi_{k}^0(R)$, which requires a careful proof.

\begin{lemma} \label{result-recursion-gap}
Suppose $F:\mathbb{R}^p \to \mathbb{R}$ satisfies \cref{assumption-lower-bound,assumption-gradient}. Let $x_0 \in \mathbb{R}^p$ and let $\lbrace x_k : k \in \mathbb{N} \rbrace$ be generated by \cref{eqn-recursion-iterate} satisfying \cref{property-spd}. Then, for all $ R \geq 0$ there exists a constant $C_R > 0$ such that
\begin{equation*}
[F(x_{k+1}) - F_{l.b.}]\chi_{k+1}^0(R) \leq \left[F(x_k) - F_{l.b.} - \dot F(x_k)^\intercal M_k \dot F(x_k) + C_R \norm{ M_k \dot F(x_k)}_2^2\right]\chi_{k}^0(R).
\end{equation*}
\end{lemma}
\begin{proof}
First, for any $R \geq 0$ define $L_R$ to be the Lipschitz constant for the gradient in the closed ball of radius $R$ around the point $0 \in \mathbb{R}^p$, which is well defined by \cref{assumption-gradient}. It readily follows that $L_R \leq L_{R+1}$. Second, let $L(x,y)$ be the Lipschitz constant for the gradient in a closed ball of radius $\inlinenorm{ x - y}_2$ around the point $y$. Finally, for any $R \geq 0$, define $G_R$ to be the maximum  $\inlinenorm{ \dot F(x)}_2$ for all $x$ in a closed ball of radius $R$ around the point $0 \in \mathbb{R}^p$.

By Taylor's remainder theorem,
\begin{equation*}
\begin{aligned}
F(x_{k+1}) - F_{l.b.} &= F(x_k) - F_{l.b.} + \dot F(x_k)^\intercal (x_{k+1} - x_k) \\
&\quad + \int_0^1 \left[ \dot F(x_k + t(x_{k+1} - x_k)) - \dot F(x_k) \right]^\intercal (x_{k+1} - x_k) dt.
\end{aligned}
\end{equation*}
By applying \cref{assumption-gradient} to the last term and substituting $x_{k+1} - x_k$ with \cref{eqn-recursion-iterate},
\begin{equation}\label{eqn-no-lipschitz-control}
\begin{aligned}
    F(x_{k+1}) - F_{l.b.} &\leq F(x_k) - F_{l.b.} + \dot F(x_k)^\intercal (x_{k+1} - x_k) + \frac{L(x_{k+1},x_k)}{2} \norm{ x_{k+1} - x_k}_2^2\\
    &= F(x_k) - F_{l.b.} - \dot F(x_k)^\intercal M_k \dot F(x_k) + \frac{L(x_{k+1},x_k)}{2} \norm{M_k \dot F(x_k)}_2^2.
\end{aligned}
\end{equation}
Note, to understand why we must keep going at this point in the proof, see \cref{discussion-no-lipschitz-control}. We now introduce $\chi_{k+1}^0(R)$ into \cref{eqn-no-lipschitz-control}. That is,
\begin{equation*}
\begin{aligned}
&[F(x_{k+1}) - F_{l.b.}] \chi_{k+1}^0(R) \\
&\quad \leq \left[F(x_k) - F_{l.b.} - \dot F(x_k)^\intercal M_k \dot F(x_k) + \frac{L(x_{k+1},x_k)}{2} \norm{M_k \dot F(x_k)}_2^2\right] \chi_{k+1}^0(R).
\end{aligned}
\end{equation*}
If $\chi_{k+1}^0(R) = 1$, then $\inlinenorm{x_{k+1}}_2 \leq R$ and $\inlinenorm{x_k}_2 \leq R$. Thus, $L(x_{k+1},x_k) \leq L_R \leq L_{R+1}$. When $\chi_{k+1}^0(R) = 0$, then both sides are trivially zero. Therefore,
\begin{equation*}
\begin{aligned}
&[F(x_{k+1}) - F_{l.b.}] \chi_{k+1}^0(R) \\
& \quad \leq \left[F(x_k) - F_{l.b.} - \dot F(x_k)^\intercal M_k \dot F(x_k) + \frac{L_{R+1}}{2} \norm{M_k \dot F(x_k)}_2^2\right] \chi_{k+1}^0(R).
\end{aligned}
\end{equation*}
However, we want $\chi_{k}^0(R)$ on the right hand side. So, we simply add and subtract $\chi_{k}^0(R)$ and then study the difference term. That is,
\begin{equation} \label{eqn-add-subtract-event}
\begin{aligned}
&[F(x_{k+1}) - F_{l.b.}] \chi_{k+1}^0(R) \\
&\quad \leq \left[F(x_k) - F_{l.b.} - \dot F(x_k)^\intercal M_k \dot F(x_k) + \frac{L_{R+1}}{2} \norm{M_k \dot F(x_k)}_2^2\right] \chi_{k}^0(R) \\
&\qquad  + \left[F(x_k) - F_{l.b.} - \dot F(x_k)^\intercal M_k \dot F(x_k) + \frac{L_{R+1}}{2} \norm{M_k \dot F(x_k)}_2^2\right] [\chi_{k+1}^0(R) - \chi_{k}^0(R)].
\end{aligned}
\end{equation}
We now have two cases to upper bound the last term of \cref{eqn-add-subtract-event}.
\paragraph{Case 1} Since $\chi_{k+1}^0(R) \leq \chi_{k}^0(R)$, if
\begin{equation*}
F(x_k) - F_{l.b.} - \dot F(x_k)^\intercal M_k \dot F(x_k) + \frac{L_{R+1}}{2} \norm{M_k \dot F(x_k)}_2^2 > 0,
\end{equation*}
then
\begin{equation*}
\left[F(x_k) - F_{l.b.} - \dot F(x_k)^\intercal M_k \dot F(x_k) + \frac{L_{R+1}}{2} \norm{M_k \dot F(x_k)}_2^2\right] [\chi_{k+1}^0(R) - \chi_{k}^0(R)] \leq 0.
\end{equation*}
Hence, in this case, we can upper bound the last term in \cref{eqn-add-subtract-event} by any non-negative term.
\paragraph{Case 2} If
\begin{equation} \label{eqn-add-subtrct-event-case-2}
F(x_k) - F_{l.b.} - \dot F(x_k)^\intercal M_k \dot F(x_k) + \frac{L_{R+1}}{2} \norm{M_k \dot F(x_k)}_2^2 \leq 0,
\end{equation}
then, since $\chi_{k+1}^0(R) \leq \chi_{k}^0(R)$,
\begin{equation*}
\left[F(x_k) - F_{l.b.} - \dot F(x_k)^\intercal M_k \dot F(x_k) + \frac{L_{R+1}}{2} \norm{M_k \dot F(x_k)}_2^2\right] [\chi_{k+1}^0(R) - \chi_{k}^0(R)] \geq 0.
\end{equation*}
Thus, we need only find a lower bound for the first term in the product to upper bound the entire term. Specifically,
\begin{equation*}
- \norm{ \dot F(x_k)}_2 \norm{ M_k \dot F(x_k)}_2 \leq F(x_k) - F_{l.b.} - \dot F(x_k)^\intercal M_k \dot F(x_k) + \frac{L_{R+1}}{2} \norm{M_k \dot F(x_k)}_2^2.
\end{equation*}
Now, when $\chi_{k+1}^0(R) < \chi_k^0(R)$, $\inlinenorm{x_k}_2 \leq R$. Moreover, \cref{eqn-add-subtrct-event-case-2} is only true if $R+1 \leq \inlinenorm{ x_{k+1}}_2$. Using \cref{eqn-recursion-iterate} and the triangle inequality, $R+1 \leq \inlinenorm{x_k}_2 + \inlinenorm{M_k \dot F(x_k)}_2 \leq R + \inlinenorm{M_k \dot F(x_k)}_2$. That is, $1 \leq \inlinenorm{ M_k \dot F(x_k)}_2 \leq \inlinenorm{ M_k \dot F(x_k)}_2^2$.

Hence,
\begin{align*}
&\left[F(x_k) - F_{l.b.} - \dot F(x_k)^\intercal M_k \dot F(x_k) + \frac{L_{R+1}}{2} \norm{M_k \dot F(x_k)}_2^2\right] [\chi_{k+1}^0(R) - \chi_{k}^0(R)] \nonumber \\
 & \quad \leq \norm{\dot F(x_k)}_2 \norm{ M_k \dot F(x_k)}_2^2 \left[ \chi_{k}^0(R) - \chi_{k+1}^0(R)\right] \\
 & \quad \leq G_R  \norm{ M_k \dot F(x_k)}_2^2 \chi_{k}^0(R).
\end{align*}

Putting these two cases together in \cref{eqn-add-subtract-event}, we conclude,
\begin{equation*}
\begin{aligned}
&[F(x_{k+1}) - F_{l.b.}] \chi_{k+1}^0(R) \\
 &\quad \leq \left[F(x_k) - F_{l.b.} - \dot F(x_k)^\intercal M_k \dot F(x_k) + \left(\frac{L_{R+1}}{2} + G_R \right) \norm{M_k \dot F(x_k)}_2^2\right] \chi_{k}^0(R).
\end{aligned}
\end{equation*}
Letting $C_R = L_{R+1}/2 + G_R$, the conclusion follows.
\end{proof}

If \cref{property-diminishing} holds, then we can simplify \cref{result-recursion-gap} as follows.

\begin{corollary} \label{result-recursion-gap-simplified}
Suppose the conditions of \cref{result-recursion-gap} hold. If \cref{property-diminishing} holds, then there exists $ K \in \mathbb{N}$ such that for all $k \geq K$,
\begin{equation*}
[F(x_{k+1}) - F_{l.b.}]\chi_{k+1}^0(R) \leq \left[F(x_k) - F_{l.b.} - \frac{1}{2}\lambda_{\min}(M_k) \norm{\dot F(x_k)}_2^2 \right]\chi_{k}^0(R).
\end{equation*}
\end{corollary}
\begin{proof}
From \cref{result-recursion-gap}, we need to show that $\exists K \in \mathbb{N}$ such that $\forall k \geq K$,
\begin{equation} \label{eqn-recursion-gap-simplified-sufficient}
- \dot F(x_k)^\intercal M_k \dot F(x_k) + C_R \norm{ M_k \dot F(x_k)}_2^2 \leq -\frac{1}{2} \lambda_{\min}(M_k) \norm{ \dot F(x_k)}_2^2.
\end{equation}
To this end, we prove, if $M$ is symmetric positive definite with $\lambda_{\max}(M) < 2/C_R$ then
\begin{equation*}
\max_{ v : \inlinenorm{v}_2 = 1} - v^\intercal M v + C_R v^\intercal MM v \leq - \frac{1}{2} \lambda_{\min}(M).
\end{equation*}
Let $0 < \lambda_{\min}(M) = \lambda_p \leq \lambda_{p-1} \leq \cdots \leq \lambda_2 \leq \lambda_1 = \lambda_{\max}(M) < 2/C_R$. Using the Schur Decomposition, there exists an orthogonal matrix $Q$ such that
\begin{equation*}
- v^\intercal M v + C_R v^\intercal MM v = \sum_{\ell=1}^p ( - \lambda_\ell + C_R \lambda_{\ell}^2 ) w_\ell^2,
\end{equation*}
where $w_\ell$ is the $\ell^\mathrm{th}$ component of $Qv$ (note, $\inlinenorm{w}_2 = \inlinenorm{Qv}_2 = \inlinenorm{v}_2 = 1$). Since $\lambda_{\ell} < 2/C_R$, it follows that $C_R \lambda_{\ell}^2 < \lambda_{\ell}/2$. Subtracting $\lambda_{\ell}$ from both sides, $-\lambda_{\ell} + C_R \lambda_{\ell}^2 < - \lambda_{\ell}/2 \leq -\lambda_{\min}(M)/2$. Thus,
\begin{equation*}
- v^\intercal M v + C_R v^\intercal MM v \leq - \lambda_{\min}(M)/2 \sum_{\ell=1}^p w_\ell^2 = - \lambda_{\min}(M)/2.
\end{equation*}
Since $\lambda_{\max}(M_k) \to 0$, there exists a $K \in \mathbb{N}$ such that $\forall k \geq K$, \cref{eqn-recursion-gap-simplified-sufficient} holds.
\end{proof}

\subsection{Applying Zoutendjik's Analysis Approach} \label{subsection-zoutendjik}

We now apply the recursive relationship established in \cref{result-recursion-gap-simplified} to study the gradient using Zoutendjik's analysis method \cite{zoutendijk1960}. As a minor bonus, we also add a statement about the behavior of the optimality gap.

\begin{theorem} \label{result-zoutendjik}
Suppose $F:\mathbb{R}^p \to \mathbb{R}$ satisfies \cref{assumption-lower-bound,assumption-gradient}. Let $x_0 \in \mathbb{R}^p$ and let $\lbrace x_k : k \in \mathbb{N} \rbrace$ be generated by \cref{eqn-recursion-iterate} satisfying \cref{property-spd,property-diverge,property-diminishing}. Then, for all $ R \geq 0$,
\begin{equation}
\lim_{k \to \infty} F(x_k)\chi_k^0(R) ~\text{exists and is finite, and}~
\liminf_{k \to \infty} \norm{ \dot F(x_k) }_2 \chi_k^0(R) = 0.
\end{equation}
If $\sup_{k} \inlinenorm{ x_k}_2 < \infty$ then $\lim_{k \to \infty} F(x_k)$ exists and is finite, and $\liminf_{k \to \infty} \inlinenorm{ \dot F(x_k) }_2 = 0$.
\end{theorem}
\begin{proof}
Let $R \geq 0$. The hypotheses of \cref{result-recursion-gap-simplified} are satisfied, and its conclusion is used freely herein. For the objective function, $\exists K \in \mathbb{N}$ such that for all $k \geq K$, $[ F(x_{k+1}) - F_{l.b.} ] \chi_{k+1}^0(R) \leq [F(x_k) - F_{l.b.}] \chi_k^0(R)$. As $\lbrace [ F(x_k) - F_{l.b.} ] \chi_k^0(R) : k \geq K \rbrace$ is a nonincreasing sequence bounded from below, it converges. Now, if we further assume that $\sup_k \inlinenorm{x_k}_2 < \infty$, then there exists an $R > 0$ such that $\chi_k^0(R) = 1$ $\forall k+1 \in \mathbb{N}$. Hence, $\lim_{k \to \infty} F(x_k) - F_{l.b.}$ exists and is finite.

For the gradient function, applying the conclusion of \cref{result-recursion-gap-simplified} and rearranging terms, for all $k \geq K$,
\begin{equation*}
\frac{1}{2}\lambda_{\min}(M_k) \norm{ \dot F(x_k) }_2^2 \chi_k^0(R) \leq [F(x_{k}) - F_{l.b.}] \chi_{k}^0(R) - [ F(x_{k+1}) - F_{l.b.} ] \chi_{k+1}^0(R).
\end{equation*}
Letting $j \geq K$ and using $F(x_{j+1}) - F_{l.b.} \geq 0$,
\begin{equation} \label{eqn-zoutendjik-telescoping-sum}
\sum_{k=K}^j \frac{1}{2}\lambda_{\min}(M_k) \norm{ \dot F(x_k) }_2^2 \chi_k^0(R) \leq [F(x_{K}) - F_{l.b.}] \chi_{K}^0(R).
\end{equation}
Now for a contradiction, suppose $\exists \; c > 0$ such that $\liminf_{k \to \infty} \inlinenorm{ \dot F(x_k)}_2^2 \chi_k^0(R) > c$. Then, there exists a $K' > K$ such that
\begin{equation*}
\frac{c}{2} \sum_{k=K'}^j \lambda_{\min}(M_k) \leq \sum_{k=K'}^j \frac{1}{2}\lambda_{\min}(M_k) \norm{ \dot F(x_k) }_2^2 \chi_k^0(R) \leq [F(x_{K}) - F_{l.b.}] \chi_{K}^0(R) < \infty.
\end{equation*}
By \cref{property-diverge}, we have a contradiction. This part of the result follows for any $R \geq 0$.

Now, if $\sup_{k} \inlinenorm{x_k}_2 < \infty$ then there exists an $R > 0$ such that $\sup_{k} \inlinenorm{ x_k}_2 < R$. Therefore, $\chi_k^0(R) = 1$ for all $k+1 \in \mathbb{N}$, and, thus, the final part of the result follows.
\end{proof}

\begin{remark} \label{discussion-no-lipschitz-control}
Suppose we directly attempt to use Zoutendjik's analysis approach in \cref{eqn-no-lipschitz-control}. We begin by rearranging \cref{eqn-no-lipschitz-control} and summing up to $j \in \mathbb{N}$ to conclude,
\begin{equation*}
\sum_{k=0}^j \dot F(x_k)^\intercal M_k \dot F(x_k) - \frac{L(x_{k+1},x_k)}{2} \norm{M_k \dot F(x_k)}_2^2 \leq  F(x_0) - F_{l.b.}.
\end{equation*}
Thus, we conclude
\begin{equation*}
\lim_{k \to \infty} \dot F(x_k)^\intercal M_k \dot F(x_k) - \frac{L(x_{k+1},x_k)}{2} \norm{M_k \dot F(x_k)}_2^2 = 0.
\end{equation*}
Unfortunately, this conclusion does not imply that $\dot F(x_k) \to 0$ as $k \to 0$. For instance, suppose as $k \to \infty$, $L(x_{k+1},x_k) \to \infty$. If $M_k = 2L(x_{k+1},x_k)^{-1} I$ for all $k$, then a straightforward substitution will show that the limit is satisfied yet $\dot F(x_k)$ does not have to be zero. Hence, using Zoutendjik's analysis method on this line of logic would not produce the desired conclusion. However, as shown in \cref{result-zoutendjik}, using Zoutendjik's analysis method on the conclusion of \cref{result-recursion-gap} is fruitful.
\end{remark}

\subsection{Convergence of the Gradient}

One of the limitations of \cref{result-zoutendjik} is that it only states that the limit infimum of the gradient function goes to zero. By using \cref{property-summable} in place of \cref{property-diminishing}, we can show that the limit supremum agrees with the limit infimum, which we establish in the next result.

\begin{theorem} \label{result-gradient-convergence}
Suppose $F:\mathbb{R}^p \to \mathbb{R}$ satisfies \cref{assumption-lower-bound,assumption-gradient}. Let $x_0 \in \mathbb{R}^p$ and let $\lbrace x_k : k \in \mathbb{N} \rbrace$ be generated by \cref{eqn-recursion-iterate} satisfying \cref{property-spd,property-diverge,property-summable}. Then, for all $ R \geq 0$,
\begin{equation}
\lim_{k \to \infty} \norm{ \dot F(x_k) }_2 \chi_k^0(R) = 0.
\end{equation}
Moreover, if $\sup_{k} \inlinenorm{ x_k}_2 < \infty$ then $\lim_{k \to \infty} \inlinenorm{ \dot F(x_k) }_2 = 0$.
\end{theorem}
\begin{proof}
By \cref{result-zoutendjik}, for every $R \geq 0$, $\liminf_{k} \inlinenorm{ \dot F(x_k) }_2 \chi_k^0(R) = 0$. For a contradiction, suppose $\exists \  \delta > 0$ such that $\limsup_{k} \inlinenorm{ \dot F(x_k) }_2 \chi_k^0(R) > \delta$. Then, for some $\gamma > 0$ sufficiently small, there exists a subsequence $\lbrace k_j : j \in \mathbb{N} \rbrace \subset \mathbb{N}$ such that $\inlinenorm{ \dot F(x_{k_j+1})}_2 \chi_{k_j+1}^0(R) \geq \delta + \gamma$ while $\inlinenorm{ \dot F(x_{k_j})}_2 \chi_{k_j}^0(R) \leq \delta$. Thus,
\begin{align*}
\delta + \gamma &\leq \norm{ \dot F(x_{k_j+1})}_2 \chi_{k_j+1}^0(R) \\
			   &\leq \norm{ \dot F(x_{k_j+1}) - \dot F(x_{k_j})}_2 \chi_{k_j+1}^0(R) + \norm{ \dot F(x_{k_j} }_2 \chi_{k_j+1}^0(R) \\
			   &\leq L(x_{k_j+1},x_{k_j}) \norm{ x_{k_j+1} - x_{k_j}}_2 \chi_{k_j + 1}^0(R) + \norm{ \dot F(x_{k_j})}_2 \chi_{k_j}^0(R) \\
			   &\leq L_R \norm{ M_{k_j} \dot F(x_{k_j})}_2 \chi_{k_j + 1}^0(R) + \delta \\
			   &\leq L_R \delta \lambda_{\max}(M_{k_j}) + \delta.
\end{align*}
Subtracting $\delta$ from both sides we have $\gamma < L_R \delta \lambda_{\max}(M_{k_j})$, from which it follows that $\gamma^q < L_R^q \delta^q \lambda_{\max}(M_{k_j})^q$, for $q\geq 0$. Summing over all $j$, $\sum_{j=1}^\infty \gamma^q < L_R^q \delta^q \sum_{j=1}^\infty \lambda_{\max}(M_{k_j})^q$. The sum on the left diverges while the sum on the right is finite by \cref{property-summable} for any appropriate choice of $q$, which leads to a contradiction. Hence, the limit supremum and the limit infimum must agree and the first conclusion of the statement follows.

Finally, when $\sup_{k} \inlinenorm{ x_k}_2 < \infty$, there exists an $R > 0$ such that $\chi_{k}^0(R) = 1$ for all $k+1 \in \mathbb{N}$, and thus the final conclusion of the statement follows.
\end{proof}

\subsection{Convergence of the Iterates}

Using a nearly identical argument as for the convergence of the gradient, we show that the iterates either converge to a point or they diverge. Again, \cref{property-summable} is essential to this conclusion.

\begin{theorem} \label{result-iterate-convergence}
Suppose $F:\mathbb{R}^p \to \mathbb{R}$ satisfies \cref{assumption-gradient}. Let $x_0 \in \mathbb{R}^p$ and let $\lbrace x_k : k \in \mathbb{N} \rbrace$ be generated by \cref{eqn-recursion-iterate} satisfying \cref{property-spd,property-summable}. Then, $\lbrace x_k : k+1 \in \mathbb{N} \rbrace$ converges to a fixed point or it diverges.
\end{theorem}
\begin{proof}
Let $z \in \mathbb{R}^p$. We study the sequence $\lbrace \inlinenorm{ x_k - z}_2 \rbrace$. Suppose for the purposes of a  contradiction that  $\exists \ R > 0$ such that $\liminf_k \inlinenorm{ x_k - z}_2 < R < \limsup_k \inlinenorm{ x_k - z}_2$. Then,  $\exists \ \gamma > 0$ and a subsequence $\lbrace k_j : j \in \mathbb{N} \rbrace$ such that $\inlinenorm{ x_{k_j} - z }_2 \leq R < R + \gamma \leq \inlinenorm{ x_{k_j + 1} - z }_2$. Notice,
\begin{align*}
R + \gamma
&\leq \inlinenorm{ x_{k_j + 1} - z }_2 \\
&\leq \inlinenorm{ x_{k_j + 1 } - x_{k_j} }_2 + \inlinenorm{ x_{k_j} - z}_2 \\
&\leq \inlinenorm{ M_{k_j} \dot F(x_{k_j})}_2 + R \\
&\leq \lambda_{\max}(M_{k_j}) G_{R,z} + R,
\end{align*}
where $G_{R,z} = \max_{ x : \inlinenorm{x - z}_2 \leq R} \inlinenorm{ \dot F(x)}_2$, which is finite by continuity and compactness. Hence, $\gamma \leq \lambda_{\max}(M_{k_j}) G_{R,z}$. In particular, $\sum_{j=1}^\infty \gamma^q \leq G_{R,z}^q \sum_{j=1}^\infty \lambda_{\max}(M_{k_j})^q$, for $q\geq0$. The sum on the left hand side is infinite, while the sum on the right hand side is finite by \cref{property-summable} for any appropriate choice of $q$, which is a contradiction. Therefore, the limit supremum and limit infimum must agree.

This leaves two possibilities: $\lim_{k \to \infty} \inlinenorm{ x_k - z}_2$ is either infinite or finite. In the case that the limit is infinite, there is nothing more that we need to do as this implies that $\lbrace x_k \rbrace$ must diverge. In the case that the limit is finite, we can replace $z$ by any $p+1$ distinct points, $z_1,\ldots,z_{p+1} \in \mathbb{R}^p$, whose affine hull is $\mathbb{R}^p$. Then $\lim_{k \to \infty} \inlinenorm{x_k - z_j }_2 < \infty$ for $j=1,\ldots,p+1$. By triangulation, we conclude that $\lbrace x_k \rbrace$ must converge to a fixed point.
\end{proof}

Combining the conclusions of  \cref{result-zoutendjik,result-gradient-convergence,result-iterate-convergence}, we have the following summary.
\begin{corollary} \label{result-zoutendjik-summable}
Suppose $F:\mathbb{R}^p \to \mathbb{R}$ satisfies \cref{assumption-lower-bound,assumption-gradient}. Let $x_0 \in \mathbb{R}^p$ and let $\lbrace x_k : k \in \mathbb{N} \rbrace$ be generated by \cref{eqn-recursion-iterate} satisfying \cref{property-spd,property-diverge,property-summable}. Then, either $\lbrace x_k : k+1 \in \mathbb{N} \rbrace$ converges to a stationary point or diverges.
\end{corollary}
\begin{proof}
By \cref{result-iterate-convergence}, $\lbrace x_k \rbrace$ either converges to a point, $x^*$, or diverges. Moreover, if $\lbrace x_k \rbrace$ converges to $x^*$, then the iterates remain bounded.  By \cref{result-gradient-convergence} and \cref{assumption-gradient}, $\dot F(x_k) \to 0$ and $\dot F(x_k) \to \dot F(x^*)$, which implies $\dot F(x^*) = 0$.
\end{proof}

\subsection{Optimization Problem}

Consider the optimization problem
\begin{equation} \label{eqn-problem}
\min_{x \in \mathbb{R}^p} F(x),
\end{equation}
under the following assumptions.
\begin{assumption} \label{assumption-lower-bound}
The objective function, $F:\mathbb{R}^p \to \mathbb{R}$, is bounded from below by a constant $F_{l.b.}$.
\end{assumption}
\begin{assumption} \label{assumption-gradient}
The gradient function, $\dot F(x) = \nabla F(z) \vert_{z = x}$, exists $\forall x \in \mathbb{R}^p$, and is locally Lipschitz continuous.
\end{assumption}

The optimization problem, \cref{eqn-problem}, that satisfies \cref{assumption-lower-bound,assumption-gradient} is relatively general, and applies to a variety of problems in data science \cite[see Appendix A]{patel2021global}. Importantly, any optimization problem whose objective function is twice continuously differentiable immediately satisfies \cref{assumption-gradient}, as is proven in the next well-known lemma.

\begin{lemma} \label{result-hessian-implies-lipschitz}
Suppose $F$ is twice continuously differentiable for all $x \in \mathbb{R}^p$. Then $\dot F(x)$ is locally Lipschitz continuous.
\end{lemma}
\begin{proof}
Let $\ddot F(x)$ denote the Hessian of $F$. Then, by assumption, $\inlinenorm{ \ddot F(x)}_2$ is a continuous function and it is bounded over any compact region. By Taylor's theorem, for any $x,y \in \mathbb{R}^p$,
$\dot F(x) - \dot F(y) = \int_{0}^1 \ddot{F}(y + t (x - y) ) (x -y) dt$.
Let $K \subset \mathbb{R}^p$ be compact. By continuity and compactness, there exists an $L$ for $K$ such that $\inlinenorm{ \ddot F(x)}_2 \leq L$ for all $x \in K$. Hence, by H\"{o}lder's inequality, for any $x, y \in K$,
$\inlinenorm{ \dot F(x) - \dot F(y) } \leq L\inlinenorm{x - y}_2$.
As $K$ is arbitrary, the result follows.
\end{proof}

\subsection{Gradient Descent with Diminishing Step Sizes}

Now, suppose we apply gradient descent with diminishing step-sizes to solve \cref{eqn-problem}. Specifically, given $x_0 \in \mathbb{R}^p$, we generate a sequence $\lbrace x_k : k \in \mathbb{N} \rbrace$ according to
\begin{equation} \label{eqn-recursion-iterate}
x_{k+1} = x_k - M_k \dot F(x_k),
\end{equation}
where $M_k$ satisfies some of the following properties:
\begin{property} \label{property-spd}
$\lbrace M_k : k + 1 \in \mathbb{N} \rbrace \subset \mathbb{R}^{p \times p}$ are symmetric positive definite matrices;
\end{property}
\begin{property} \label{property-diverge}
$\sum_{k=0}^\infty \lambda_{\min}(M_k)$ diverges, where $\lambda_{\min}(M_k)$ denotes the smallest eigenvalue of $M_k$; and \textbf{either}
\end{property}
\begin{property} \label{property-diminishing}
$\lim_{k \to \infty} \lambda_{\max}(M_k) = 0$ where $\lambda_{\max}(M_k)$ denotes the largest eigenvalue of $M_k$, \textbf{or}
\end{property}
\begin{property} \label{property-summable}
$\exists q > 1$ such that $\sum_{k=0}^\infty \lambda_{\max}(M_k)^q < \infty$.
\end{property}

\cref{property-spd,property-diverge,property-diminishing} are a natural generalization of classical diminishing step size requirements \cite[Proposition 1.2.4]{bertsekas2016}. Moreover, owing to its contrast with \cref{property-summable}, we refer to requiring properties \cref{property-spd,property-diverge,property-diminishing} as the \textit{without summability} case. \cref{property-spd,property-diverge,property-summable}  have their roots in the Robbins-Monro conditions \cite{robbins1951} and are closely related to the properties in \cite{patel2020,patel2021global,patel2021}. Despite its heritage, \cref{property-summable} is an important generalization to those of \cite{robbins1951,patel2020,patel2021global,patel2021} as we can choose any $q > 2$. Indeed, \cref{property-summable} is not enough for the stochastic gradient case, as summability for $q = 2$ is required.\footnote{This is the case when we assume the variances of the stochastic gradients exist. While we can permit even lower order moments to exist (i.e., less than the second moment of the stochastic gradient), we would then need more stringent requirements on the value of $q$ in \cref{property-summable}. See \cite{patel2021global,patel2021}.} Thus, \cref{property-summable} allows for a much larger selection of step sizes.

To be concrete, consider $M_k = (k+1)^{-1/4} I$, where $I \in \mathbb{R}^{p \times p}$ is the identity matrix. Then, $\lbrace M_k \rbrace$ satisfies \cref{property-diverge} and \cref{property-summable} with $q > 4$, but does not satisfy \cref{property-summable} for $q = 2$ as would be needed for stochastic gradient descent. As an informative counterexample, \cref{property-summable} does not apply to all decaying sequences of $\lbrace M_k \rbrace$, as the next proposition shows.
\begin{proposition}
For all $k + 1\in \mathbb{N}$, let $M_k = \log(k+1)^{-1} I$, where $I \in \mathbb{R}^{p \times p}$ is the identity matrix. Then, $\lbrace M_k : k+1 \in \mathbb{N} \rbrace$ satisfies \cref{property-diverge}, but fails to satisfy \cref{property-summable}.
\end{proposition}
\begin{proof}
It is enough to show that for all $q \geq 1$, $\sum_{k=1}^\infty \log(k)^{-q}$ diverges. For any $q \geq 1$, it is straightforward to verify that $\log(k)^q / k \to 0$ as $k \to \infty$. Therefore, there exists a $K \in \mathbb{N}$, such that $\log(k)^{-q} \geq k^{-1}$ for all $k \geq K$. Since $\sum_{k=K}^\infty k^{-1}$ diverges, the result follows.
\end{proof}

We now analyze the global convergence of gradient descent with diminishing step sizes on nonconvex functions specified by \cref{assumption-lower-bound,assumption-gradient}. We will split our analysis into two cases: the case without summability (\cref{property-spd,property-diverge,property-diminishing}; \cref{section-no-summable}) and the case with summability (\cref{property-spd,property-diverge,property-summable}; \cref{section-summable}).

\subsection{Construction of the Objective Function}
Let $\lbrace m_k : k+1 \in \mathbb{N}\rbrace$ be a sequence of scalars such that $m_k > 0$, $\sum_k m_k = \infty$, and $m_k \to 0$ as $k \to \infty$. Define $S_0 = 0$ and $S_{k+1} = \sum_{j=0}^k m_k$ for all integers $k \geq 0$.

For the objective function, define $F: \mathbb{R} \to \mathbb{R}$ by
\begin{equation} \label{eqn-divergence-objective}
F(x) = \begin{cases}
-x & x \leq 0 \\
f_j(x) & x \in (S_j, S_{j+1}], ~\forall j+1 \in \mathbb{N},
\end{cases}
\end{equation}
where $\lbrace f_j: (S_j, S_{j+1}] \to \mathbb{R} : j + 1\in\mathbb{N} \rbrace$ are defined iteratively as follows. Let
\begin{equation} \label{eqn-divergence-component-base}
f_0(x) = \begin{cases}
-x & x \in (0, \frac{m_0}{16}) \\
\frac{8}{m_0}(x- \frac{m_0}{8})^2 - \frac{3m_0}{32} & x \in [\frac{m_0}{16}, \frac{3m_0}{16}) \\
- \frac{5m_0}{16} \exp\left( \frac{5m_0/16}{x - m_0/2 } + 1 \right) + \frac{m_0}{4} & x \in [\frac{3m_0}{16}, \frac{ m_0}{2}) \\
\frac{m_0}{4} & x = \frac{m_0}{2} \\
\frac{5m_0}{16} \exp \left( -\frac{5m_0/16}{x - m_0/2} + 1 \right) + \frac{m_0}{4} & x \in (\frac{m_0}{2}, \frac{13m_0}{16}) \\
\frac{-8}{m_0}(x - \frac{7m_0}{8})^2 + \frac{19m_0}{32} & x \in [\frac{13 m_0}{16}, \frac{15m_0}{16}) \\
-x + \frac{3 m_0}{2} & x \in [\frac{15 m_0}{16}, m_0],
\end{cases}
\end{equation}
which is plotted for a particular choice of $m_0$ in \cref{plot-divergent-example-component}. Now, for $j \in \mathbb{N}$, let
\begin{equation} \label{eqn-divergence-component-general}
f_j(x) = \begin{cases}
-x + S_j + f_{j-1}(S_j) & x \in (S_j, S_j + \frac{m_j}{16}) \\
\frac{8}{m_j}(x - S_j - \frac{m_j}{8})^2 - \frac{3m_j}{32} + f_{j-1}(S_j) & x \in [S_j + \frac{m_j}{16}, S_j + \frac{3m_j}{16}) \\
- \frac{5m_j}{16} \exp\left( \frac{5m_j/16}{x - S_j - m_j/2} + 1\right) + \frac{m_j}{4} + f_{j-1}(S_j) & x \in [S_j + \frac{3m_j}{16}, S_j + \frac{m_j}{2}) \\
\frac{m_j}{4} + f_{j-1}(S_j) & x = S_j + \frac{m_j}{2} \\
\frac{5m_j}{16} \exp\left( \frac{-5m_j/16}{x - S_j - m_j/2} + 1\right) + \frac{m_j}{4} + f_{j-1}(S_j) & x \in (S_j + \frac{m_j}{2}, S_j + \frac{13m_j}{16}) \\
\frac{-8}{m_j}(x - S_j - \frac{7m_j}{8})^2 + \frac{19m_j}{32} + f_{j-1}(S_j) & x \in [S_j + \frac{13m_j}{16}, S_j + \frac{15m_j}{16}) \\
-x + S_j + \frac{3m_j}{2} + f_{j-1}(S_j) & x \in [S_j + \frac{15m_j}{16}, S_{j+1}].
\end{cases}
\end{equation}

\begin{figure}[h]
\centering
\input{figure/tikz-divergent-example-component}
\caption{Plot of $f_0(x)$ with $m_0 = 8.0$ with each component shown in a different color.}
\label{plot-divergent-example-component}
\end{figure}

\subsection{Properties of the Objective Function}

We show that $F: \mathbb{R} \to \mathbb{R}$ as defined in \cref{eqn-divergence-objective}--\cref{eqn-divergence-component-general} satisfies \cref{assumption-lower-bound,assumption-gradient}. We begin by proving that each component, $f_j: (S_j, S_{j+1}] \to \mathbb{R}$, satisfies \cref{assumption-lower-bound,assumption-gradient} on its domain.

\begin{remark}
Below, we define the continuous extension of $f_j$ on $[S_j, S_{j+1}]$ by the value of $f_j(x)$ on $(S_j, S_{j+1}]$ and by $\lim_{x \downarrow S_j} f_j(x)$ for the point $x= S_j$. Moreover, at the ends of the interval, we use differentiability and the corresponding notation, $\dot f_j$, to mean the one-sided derivatives.
\end{remark}
\begin{proposition} \label{result-divergence-component-properties}
The continuous extension $f_0: (S_0, S_1] \to \mathbb{R}$ $($as defined in \cref{eqn-divergence-component-base}$)$ to $[S_0,S_1]$ is continuous on its domain, bounded from below by $-3m_0/32$, differentiable on its domain with $\dot f_0(S_0) = \dot f_0(S_1) = -1$, and its derivative is locally Lipschitz continuous. Similarly, the continuous extension of $f_j: (S_j, S_{j+1}] \to \mathbb{R}$ $($as defined in \cref{eqn-divergence-component-general}$)$ to $[S_j,S_{j+1}]$ is continuous on its domain, bounded from below by $f_{j-1}(S_j) - 3m_j/32$, differentiable on its domain with $\dot f_j(S_j) = -1$ and $\dot f_j(S_{j+1}) = -1$, and its derivative is locally Lipschitz continuous.
\end{proposition}
\begin{proof}
We only look at an arbitrary $j \in \mathbb{N}$ as the proof is identical for $f_0$. Moreover, since $f_j$ is equal to its reflection across the vertical axis $x = S_j + m_j/2$ followed by a reflection over the horizontal axis $y = m_j/4 + f_{j-1}(S_j)$, it is enough to show continuity and differentiability on $[S_j, S_j + m_j/2]$.

To establish continuity, we need to show that the left sided limits of $f_j$ agree with the function value at $S_j + \delta m_j/16$ for $\delta = 1, 3, 8$. Starting with $\delta = 1$, $\lim_{x \uparrow S_j + m_j/16} -x + S_j + f_{j-1}(S_j) = -m_j/16 + f_{j-1}(S_j)$. By direct substitution,
\begin{equation*}
f_j(S_j + m_j/16) = \frac{8}{m_j}\left( \frac{m_j}{16} \right)^2 - \frac{3m_j}{32} + f_{j-1}(S_j) = \frac{m_j}{32} - \frac{3m_j}{32} + f_{j-1}(S_j).
\end{equation*}
Hence, the left limit agrees with the function value at $\delta = 1$. For $\delta = 3$, the symmetry and continuity of the quadratic function implies $\lim_{ x \uparrow S_j + 3m_j/16} f_j(x) = -m_j/16 + f_{j-1}(S_j)$. By direct substitution,
\begin{equation*}
f_j(S_j + 3m_j/16) = -\frac{5m_j}{16} \exp\left( \frac{5m_j/16}{-5m_j/16} + 1 \right) + \frac{4m_j}{16} + f_{j-1}(S_j) = -\frac{m_j}{16} + f_{j-1}(S_j).
\end{equation*}
Hence, the left limit agrees with the function value at $\delta = 3$. For $\delta = 8$,
\begin{equation*}
\lim_{ x \uparrow S_j + m_j/2} -\frac{5m_j}{16} \exp\left( \frac{5m_j/16}{x - S_j - m_j/2} + 1 \right) + \frac{m_j}{4} + f_{j-1}(S_j) = \frac{m_j}{4} + f_{j-1}(S_j),
\end{equation*}
which is just $f_j(S_j + m_j/2)$. Hence, the continuous extension of $f_j$ is continuous on $[S_j, S_{j+1}]$.

We now compute the derivatives of the components of $f_j$ on $[S_j, S_{j+1}]$ with the convention of assigning the one-sided derivative to the component function that includes its end point.
\begin{equation} \label{eqn-divergence-component-derivatives}
\dot f_j(x) =\begin{cases}
-1 & x \in [S_j, S_j + \frac{m_j}{16}) \\
\frac{16}{m_j}(x - S_j - \frac{m_j}{8}) & x \in [S_j + \frac{m_j}{16}, S_j + \frac{3m_j}{16}) \\
\left( \frac{5m_j/16}{x - S_j - m_j/2} \right)^2 \exp\left( \frac{5m_j/16}{x - S_j - m_j/2} + 1 \right) & x \in [S_j + \frac{3m_j}{16}, S_j + \frac{m_j}{2}) \\
\left( \frac{5m_j/16}{x - S_j - m_j/2} \right)^2 \exp\left( \frac{-5m_j/16}{x - S_j - m_j/2} + 1 \right) & x \in (S_j + \frac{m_j}{2}, S_j + \frac{13m_j}{16}) \\
-\frac{16}{m_j}(x - S_j - \frac{7m_j}{8}) & x \in [S_j + \frac{13m_j}{16}, S_j + \frac{15m_j}{16}) \\
-1 & x \in [S_j + \frac{15m_j}{16}, S_{j+1}].
\end{cases}
\end{equation}
In order to extend these component derivatives to the continuous extension of $f_j$, we need to verify that the left hand limits of the derivatives agree with the right side derivatives at $S_j + \delta m_j/16$ for $\delta=1,3$ of \cref{eqn-divergence-component-derivatives}, and we need to verify that the left hand limit of the derivative is $0$ at $\delta = 8$ of \cref{eqn-divergence-component-derivatives}. Starting with $\delta = 1$, the left hand limit is $-1$ and, by direct calculation,
\begin{equation*}
\frac{16}{m_j}\left(S_j + \frac{m_j}{16} - S_j - \frac{m_j}{8} \right) = -1.
\end{equation*}
For $\delta = 3$, the left hand limit is
\begin{equation*}
\lim_{ x \uparrow S_{j} + {3m_j/16} } \frac{16}{m_j}\left(x - S_j - \frac{m_j}{8} \right) = 1,
\end{equation*}
and a direct evaluation of the third component \cref{eqn-divergence-component-derivatives} is
\begin{equation*}
\left( \frac{5m_j/16}{S_{j} + {3m_j/16} - S_j - m_j/2} \right)^2 \exp\left( \frac{5m_j/16}{S_{j} + {3m_j/16} - S_j - m_j/2} + 1 \right) = 1.
\end{equation*}
For $\delta = 8$, we need to check that the left hand limit is zero, which can be confirmed by checking that the argument of the exponential term goes to $-\infty$ as $x \uparrow S_j + 8m_j/16$.

Overall, the derivative of the continuous extension of $f_j$ is well defined at every point on its interval, is continuous, and is given by
\begin{equation} \label{eqn-divergence-component-derivative-definition}
\dot f_j(x) = \begin{cases}
-1 & x \in [S_j, S_j + \frac{m_j}{16}) \\
\frac{16}{m_j}(x - S_j - \frac{m_j}{8}) & x \in [S_j + \frac{m_j}{16}, S_j + \frac{3m_j}{16}) \\
\left( \frac{5m_j/16}{x - S_j - m_j/2} \right)^2 \exp\left( \frac{5m_j/16}{x - S_j - m_j/2} + 1 \right) & x \in [S_j + \frac{3m_j}{16}, S_j + \frac{m_j}{2}) \\
0 & x = S_j + m_j/2 \\
\left( \frac{5m_j/16}{x - S_j - m_j/2} \right)^2 \exp\left( \frac{-5m_j/16}{x - S_j - m_j/2} + 1 \right) & x \in (S_j + \frac{m_j}{2}, S_j + \frac{13m_j}{16}) \\
-\frac{16}{m_j}(x - S_j - \frac{7m_j}{8}) & x \in [S_j + \frac{13m_j}{16}, S_j + \frac{15m_j}{16}) \\
-1 & x \in [S_j + \frac{15m_j}{16}, S_{j+1}].
\end{cases}
\end{equation}

Using this derivative, we can calculate the lower bound for the function. By the derivative of the extension of $f_j$, \cref{eqn-divergence-component-derivative-definition}, we see that the function is decreasing only on $[S_j, S_j + m_j/8]$ and $[S_j + 7m_j/8, S_{j+1}]$. Moreover, $f_j(S_j + m_j/8) = -3m_j/32 + f_{j-1}(S_j)$ and $f_j(S_{j+1}) = m_j/2 + f_{j-1}(S_j)$. Thus, the lower bound of the extension of $f_j$ is as stated.

Our last step is to verify the local Lipschitz continuity of the derivative. It is easy to verify that within its interval, the components are twice continuously differentiable. As a result, we can use \cref{result-hessian-implies-lipschitz}. Similarly, if we define the second derivative at $S_j + m_j/2$ to be $0$, we can verify that the objective is twice continuously differentiable at $S_j + m_j/2$. Then, we can use \cref{result-hessian-implies-lipschitz} again. To conclude, we need to examine what happens around the points $S_j + \delta m_j/16$ for $\delta = 1, 3$.

Starting with $\delta = 1$, consider the points $S_j + m_j/16 - \epsilon_1 m_j/16$ and $S_j + m_j/16 + \epsilon_2 m_j/16$ for $\epsilon_1, \epsilon_2 > 0$ sufficiently small. Then the difference in the derivatives at these points divided by the distance between the points is
\begin{align*}
\frac{\left\vert \frac{16}{m_j}\left(-\frac{m_j}{16} + \epsilon_2 \frac{m_j}{16}\right) + 1  \right\vert }{(\epsilon_2 + \epsilon_1) m_j/16} = \frac{\epsilon_2}{(\epsilon_2 + \epsilon_1) m_j/16} \leq \frac{16}{m_j}.
\end{align*}
Therefore, we conclude that the derivative is locally Lipschitz near $S_j + m_j/16$. For $\delta=3$, we compute the same ratio at the points $S_j + 3m_j/16 - \epsilon_1 m_j/16$ and $S_j + 3m_j/16 + 5 \epsilon_2 m_j/16$ for $\epsilon_1, \epsilon_2 \in [0,1/4]$ where at most either $\epsilon_1$ or $\epsilon_2$ is zero. The ratio of the difference in the derivatives and the points is
\begin{align*}
\frac{\left\vert \frac{1}{(1 - \epsilon_2)^2} \exp\left( \frac{\epsilon_2}{\epsilon_2 - 1}\right) -1 + \epsilon_1 \right\vert}{(5 \epsilon_2 + \epsilon_1) m_j/16} \leq \frac{\epsilon_2 + \epsilon_2^2/2 + \epsilon_1}{(5 \epsilon_2 + \epsilon_1) m_j/16} \leq \frac{16}{m_j}.
\end{align*}
Therefore, we conclude that the derivative is locally Lipschitz near $S_j + 3m_j/16$.
\end{proof}

With this calculation complete, we can now verify that $F$ satisfies \cref{assumption-lower-bound,assumption-gradient}.

\begin{proposition} \label{result-divergence-objective-properties}
The function $F : \mathbb{R} \to \mathbb{R}$ is continuous and differentiable on its domain, is lower bounded, and its derivative is locally Lipschitz continuous. Moreover, $F(S_j) = S_j/2$ and $\dot F(S_j) = -1$ for all $j + 1 \in \mathbb{N}$.
\end{proposition}
\begin{proof}
By \cref{result-divergence-component-properties}, we need only verify the continuity of $F$ and the differentiability of $F$ at $x = S_j$. Note, $F(S_j) = f_{j-1}(S_j)$. We need to verify that the right side limit converges to $f_{j-1}(S_j)$. That is,
\begin{equation*}
\lim_{ x \downarrow S_j } F(x) = \lim_{x \downarrow S_j} f_j(x) = \lim_{x \downarrow S_j} -x + S_j + f_{j-1}(S_j) = f_{j-1}(S_j).
\end{equation*}
Thus, $F$ is continuous at all $S_j$ for $j \in \mathbb{N}$. We check $x = S_0 = 0$ as well. $F(0) = 0$ by definition. Moreover,
\begin{equation*}
\lim_{ x\downarrow 0 } F(x) = \lim_{x \downarrow 0} f_0(x) = \lim_{ x \downarrow 0} -x = 0.
\end{equation*}
Hence, $F$ is continuous on its domain.

For differentiability, by  \cref{result-divergence-component-properties} it follows that the derivative at each $S_j$ is $-1$ from the left and the right for $j \in \mathbb{N}$. Hence, the derivative exists at each $S_j$ for $j \in \mathbb{N}$. Moreover, the left hand derivative at $S_0$ is $-1$, which agrees with the right hand derivative. Hence, $F$ is continuously differentiable on its domain. Since the derivative is constant in a small neighborhood of $S_j$, it is Lipschitz continuous in this region. Finally, $\dot F(S_j) = -1$ for all $j+1 \in \mathbb{N}$.

To show that $F$ is lower bounded, we will first calculate the values of $F(S_j)$. We proceed by induction. $F(S_0) = 0 = S_0/2$. Now suppose the statement holds up to $j$. Then $F(S_{j}) = S_{j}/2$. By construction, $f_{j-1}(S_j) = F(S_{j}) = S_j/2$. Now,
\begin{align*}
F(S_{j+1}) = f_j(S_{j+1}) &= -S_{j+1} + S_j + \frac{3m_j}{2} + f_{j-1}(S_j) \\
&= -S_j - m_j + S_j + \frac{3m_j}{2} + \frac{S_j}{2} = \frac{S_{j+1}}{2}.
\end{align*}
To show the lower bound property, recall that $f_j(x) \geq f_{j-1}(S_j) - 3m_j/32 = F(S_{j}) - 3m_j/32 = S_j/2 - 3m_j/32$. Since $S_j \to \infty$ and $m_j \to 0$ by construction, $F(x) \geq \inf_{j} S_j/2 - 3m_j/32 > -\infty$.
\end{proof}

\subsection{Properties of Gradient Descent on the Objective Function}

We are now ready to show that gradient descent with diminishing step sizes generates iterates that diverge, and whose objective function diverges and  gradient function remains bounded away from zero.

\begin{proposition} \label{result-divergence}
Let $\lbrace m_k : k+1 \in \mathbb{N} \rbrace$ be any positive sequence such that $\sum_{k} m_k$ diverges and $m_k \to 0$. Define $F : \mathbb{R} \to \mathbb{R}$ as in \cref{eqn-divergence-objective}. Suppose $x_0 = 0$ and let $\lbrace x_k : k \in \mathbb{N} \rbrace$ be generated according to \cref{eqn-recursion-iterate} with $M_k = m_k I$ for all $k+1 \in \mathbb{N}$. Then, $\lbrace M_k \rbrace$ satisfies \cref{property-spd,property-diverge,property-diminishing}. Moreover, (a) $\lim_k x_k = \infty$, (b) $\lim_k F(x_k) = \infty$, and (c) $\lim_k |\dot F(x_k)| = 1$.
\end{proposition}
\begin{proof}
To prove the result, we need only show that $x_k = S_k$ where we recall that $S_0 = 0$ and $S_k = \sum_{j=0}^{k-1} m_j$. For $k=0$, $x_0 = 0 = S_0$. Suppose this holds up to $k$. Then, by \cref{result-divergence-objective-properties},
\begin{equation*}
x_{k+1} = x_k - M_k \dot F(x_k) = S_k - m_k (-1) = S_{k+1}.
\end{equation*}
Now, since $S_k$ diverges, the iterates diverge (part $(a)$). Moreover, by \cref{result-divergence-objective-properties}, since $F(x_k) = F(S_k) = S_k/2$, the objective function also diverges (part $(b)$). Finally, by \cref{result-divergence-objective-properties}, $\dot F(x_k) = \dot F(S_k) = -1$ (part $(c)$).
\end{proof}

In summary, as the example from \cref{result-divergence} and the example of $F(x) = \exp(-x^2)$ show, under our assumptions about the objective function and properties of gradient descent, we cannot conclude anything additional about the objective behavior of the function or the  gradient in the regime where the iterates generated by gradient descent with diminishing step sizes diverge.

\subsection{Continuous Gradient Descent}
Let $x \in \mathbb{R}^p$ be arbitrary. Consider the initial value problem,
\begin{equation} \label{eqn-ivp}
\begin{cases}
\frac{dy}{dt} &= - \dot F(y(t)), ~ t \geq 0 \\
y(0) &= x.
\end{cases}
\end{equation}
To explicitly show the dependence on the initial value, $x$, we will write $\Phi_h(x)$ to be the solution (if it exists) to \cref{eqn-ivp} at time $h \geq 0$ that is initialized at $x \in \mathbb{R}^p$. If $F$ satisfies \cref{assumption-gradient}, we can apply the Picard-Lindel\"{o}f theorem to conclude that $\Phi_h(x)$ exists and is unique for any finite, non-negative choice of $h$. Therefore, we define continuous gradient descent initialized at $x \in \mathbb{R}^p$ to be $\lbrace \Phi_h(x) : h \geq 0 \rbrace$.

\subsection{Limit Cycles of Continuous Gradient Descent}
There are many possible choices for the asymptotics of continuous gradient descent; two particularly interesting cases are cycles and limit cycles. Recall, a cycle cannot be a solution for continuous gradient descent. To see why, suppose such a cycle exists and $x$ is a point on this cycle. Then there would exist a $t > 0$ such that $F(\Phi_t(x)) = F(x)$, but
\begin{equation*}
0 = F(\Phi_t(x)) - F(x) = \int_0^t \dot F(\Phi_s(x))^\intercal \left( \frac{d}{dr} \Phi_r(x) \vert_{r=s} \right) ds = \int_0^t - \norm{ \dot F(\Phi_s(x))}_2^2 ds < 0.
\end{equation*}
Hence, a cycle cannot be a solution for continuous gradient descent.

However, it is possible for continuous gradient descent to have a limit cycle. To see how, consider a small modification of the Palis-de Melo function \cite[\S 1, Example 3]{palis2012}, $F: \mathbb{R}^2 \to \mathbb{R}$, specified by
\begin{equation} \label{eqn-pdm}
F(x) = \begin{cases}
\exp\left(\frac{1}{\norm{x}_2^2 - 1}  \right) & \norm{x}_2 < 1 \\
0 & \norm{x}_2 = 1 \\
-\exp\left(\frac{-1}{\norm{x}_2^2 - 1} \right) \left[ \sin\left( \frac{1}{\norm{x}_2-1}  \right) \frac{\cmp{x}{1}}{\norm{x}_2} -  \cos\left( \frac{1}{\norm{x}_2-1}  \right) \frac{\cmp{x}{2}}{\norm{x}_2} \right] & \norm{x}_2 > 1,
\end{cases}
\end{equation}
where $\cmp{x}{j}$ is the $j^\mathrm{th}$ component of $x$. The zero set of the Palis-de Melo function is plotted in \cref{figure-palis-de-melo-zero}, and a profile of the function when $\cmp{x}{2} = 0$ is plotted in \cref{figure-palis-de-melo-profile}. Importantly, the Palis-de Melo function satisfies \cref{assumption-lower-bound,assumption-gradient}.

\begin{proposition} \label{result-pdm-assumptions}
The Palis-de Melo function, \cref{eqn-pdm}, is bounded from below, is differentiable at every point in $\mathbb{R}^2$, and its gradient is locally Lipschitz continuous at every point in $\mathbb{R}^2$.
\end{proposition}
\begin{proof}
To verify that the Palis-de Melo function is bounded from below, we look at the three components separately. When $\inlinenorm{x}_2 < 1$, $F(x) \geq 0$. When $\inlinenorm{x}_2 = 1$, $F(x) = 0$. When $\inlinenorm{x}_2 > 1$, define
\begin{equation}
\theta = \begin{cases}
\arcsin\left( \frac{\cmp{x}{2}}{\norm{x}_2} \right) & \cmp{x}{1} \geq 0 \\
\pi - \arcsin\left( \frac{\cmp{x}{2}}{\norm{x}_2} \right) & \cmp{x}{1} < 0
\end{cases}
\end{equation}
Then, when $\inlinenorm{x}_2 > 1$,
\begin{equation}
F(x) = - \exp\left(\frac{-1}{\norm{x}_2^2 - 1} \right) \sin\left( \frac{1}{\norm{x}_2 -1} - \theta \right) \geq - \exp\left(\frac{-1}{\norm{x}_2^2 - 1} \right) \geq -1.
\end{equation}
Thus, $F$ is bounded from below.

We now show that $F$ is continuously differentiable. When $\inlinenorm{x}_2 < 1$, we can differentiate to conclude
\begin{equation}
\dot F(x) = \exp\left( \frac{1}{\norm{x}_2^2 - 1} \right) \frac{-2}{\left( \norm{x}_2^2 - 1  \right)^2} x.
\end{equation}
Next, let $x \in \mathbb{R}^2$ such that $\inlinenorm{x}_2 = 1$. Then,
\begin{equation}
\dot F(x) = \lim_{x_n \to x} \frac{|F(x_n)|}{\norm{x_n - x}_2} \leq \lim_{x_n \to x } \frac{\min\left\lbrace \exp\left( \frac{1}{\norm{x_n}_2^2 -1 }  \right),\exp\left( \frac{-1}{\norm{x_n}_2^2 -1 }  \right) \right\rbrace}{\norm{x_n - x}_2} = 0.
\end{equation}
When $\inlinenorm{x}_2 > 1$, we can differentiate to conclude
\begin{equation*}
\begin{aligned}
\dot F(x) & =
\exp\left( \frac{-1}{\norm{x}_2^2 - 1} \right) \frac{2}{\left( \norm{x}_2^2 - 1  \right)^2}  \left[ \sin\left( \frac{1}{\norm{x}_2-1}  \right) \frac{\cmp{x}{1}}{\norm{x}_2} -  \cos\left( \frac{1}{\norm{x}_2-1}  \right) \frac{\cmp{x}{2}}{\norm{x}_2} \right]x \\
&\quad +
\exp\left( \frac{-1}{\norm{x}_2^2 - 1} \right)\left[
\frac{-\cmp{x}{1}}{(\norm{x}_2 - 1)^2 \norm{x}_2^2}  \cos\left( \frac{1}{\norm{x}_2-1}  \right) x
- \sin\left( \frac{1}{\norm{x}_2-1}  \right)\frac{\cmp{x}{1}}{\norm{x}_2^3}x \right. \\
&\quad + \sin\left( \frac{1}{\norm{x}_2-1}  \right) \frac{1}{\norm{x}_2}\begin{bmatrix} 1 \\ 0 \end{bmatrix}
-\frac{\cmp{x}{2}}{(\norm{x}_2 - 1)^2 \norm{x}_2^2}  \sin\left( \frac{1}{\norm{x}_2-1}  \right) x  \\
&\quad \left. + \cos\left( \frac{1}{\norm{x}_2-1}  \right)\frac{\cmp{x}{2}}{\norm{x}_2^3}x
-\cos\left( \frac{1}{\norm{x}_2-1}  \right) \frac{1}{\norm{x}_2}\begin{bmatrix} 0 \\ 1 \end{bmatrix}\right].
\end{aligned}
\end{equation*}
Finally, notice that when when $\norm{x}_2$ is near $1$, $\dot F(x)$ is dominated by the exponential term which goes to zero. Hence, $\dot F$ is continuously differentiable.

We can do the identical, yet more tedious calculation to conclude that the second derivative of $F(x)$ exists for all $x \in \mathbb{R}^2$ and is continuous. We can then apply \cref{result-hessian-implies-lipschitz} to conclude that $\dot F(x)$ is locally Lipschitz continuous.
\end{proof}

By \cref{result-pdm-assumptions}, continuous gradient descent is well-defined on the Palis-de Melo function. There is an appropriate initialization of continuous gradient descent that results in a path whose limit cycle is the circle of radius one, which we show next.

\begin{proposition} \label{result-pdm-limit-cycle}
Let $\mathcal{S} = \lbrace z \in \mathbb{R}^2 : \inlinenorm{z}_2 = 1 \rbrace$. Let $\lbrace \Phi_t(x): t \geq 0 \rbrace$ denote continuous gradient descent on \cref{eqn-pdm} when initialized at $x \in \mathbb{R}^2$. Let
\begin{equation}
\mathcal{U} = \left\lbrace \begin{bmatrix} r\cos(\theta) \\ r \sin(\theta) \end{bmatrix}: ~r\geq 0, ~\theta \geq 0, 1+ \frac{1}{2\pi+\theta} < r < 1 + \frac{1}{\pi + \theta}  \right\rbrace \subset \mathbb{R}^2.
\end{equation}
There exists an $x \in \mathbb{R}^2$ for which $\cmp{x}{1} \in (1+\frac{1}{2\pi}, 1 + \frac{1}{\pi})$ and $\cmp{x}{2} = 0$ such that $\Phi_h(x) \in \mathcal{U}$ for all $h \geq 0$. Moreover, for all $ z \in \mathcal{S}$,  $\exists \lbrace t_i : i \in \mathbb{N} \rbrace \subset \mathbb{R}_{\geq 0}$ such that $\lim_{i \to \infty} \Phi_{t_i}(x) = z$.
\end{proposition}
\begin{proof}
The proof is a small modification of that given in \cite[\S 1, Example 3]{palis2012}.
\end{proof}

The basic idea of the proof is best understood geometrically. Consider the initialization of continuous gradient descent between the purple points in \cref{figure-palis-de-melo-zero} (equivalently, \cref{figure-palis-de-melo-profile}). From \cref{figure-palis-de-melo-profile}, the initialization is in a region where the objective function is positive. Hence, as continuous gradient descent attempts to reduce the objective function value, continuous gradient descent initialized in this region should go towards the black spiral or orange spiral in \cref{figure-palis-de-melo-zero}, or towards another point between the purple points in \cref{figure-palis-de-melo-zero}. Hence, by the Picard-Lindel\"{o}f theorem (which can be applied by \cref{result-pdm-assumptions}), there is a point between the set of initializations going towards the black spiral and the other two alternatives. Thus, if continuous gradient descent is initialized at one of these points, it must remain between the black and orange spirals (i.e., stay within $\mathcal{U}$) and must have a limit cycle of $\mathcal{S}$.

\begin{figure}[hbt]
\centering
\input{figure/tikz-palis-de-melo}
\caption{A plot of the zero set of the Palis-de Melo function, \cref{eqn-pdm}. The lines indicate the sets where the function is zero. The green line corresponds to $\inlinenorm{x}_2 = 1$; the black line corresponds to one branch of the zero set of the sinusoidal term (when $\inlinenorm{x}_2 > 1$); and the orange line corresponds to the another branch of the zero set of the sinusoidal term with a phase shift of $\pi$. The pluses and minuses indicate where $F$ is positive and negative, respectively. The two purple points correspond to $(1+\frac{1}{2\pi},0)$ and $(1+\frac{1}{\pi},0)$.} \label{figure-palis-de-melo-zero}
\end{figure}

\begin{figure}[hbt]
\centering
\input{figure/tikz-palis-de-melo-profile}
\caption{A profile plot of the Palis-de Melo function, \cref{eqn-pdm}, when $\cmp{x}{2} = 0$.} \label{figure-palis-de-melo-profile}
\end{figure}

\subsection{Spurious Solution for the Palis-de Melo Function} \label{subsection-pdm-spurious}

Now, suppose we apply gradient descent to any initialization with scalar step sizes satisfying \cref{property-spd,property-diverge,property-summable}, and that the iterates generated by gradient descent remain in $\mathcal{U}$. Then, the iterates are finite and \cref{result-zoutendjik-summable} implies that they must converge to a stationary point. Thus, we have a contradiction. To conclude, continuous gradient descent has a spurious solution for gradient descent with diminishing, summable step sizes. In other words, the properties of continuous gradient descent, \textit{regardless of our initialization}, do not all apply to gradient descent with diminishing, summable step sizes.

While this proof is completely sufficient, we can provide some geometric insight into why we should expect this to hold. In particular, recall that we effectively need to spiral around the circle, $\mathcal{S}$, using a discrete method that only uses gradient information. As we get closer and closer to this circle, our gradient will be perpendicular to the circle and so, no matter our step size, we must deviate from the tracing of this circle.


\section{Introduction}

Recently, we have seen a surge of interest to bound states
in the continuum (BICs)~\cite{Hsu16, Koshelev19, sadreev2021interference} that have grown to an important
tool in nanophotonics paving a
way to optical devices with enhanced light-matter interaction.
BICs do not couple the incident light, however, if the symmetry of the system is broken the
BICs are observed as narrow Fano resonances in the scattering spectrum~\cite{Kim,Shipman,SBR,Blanchard16}.
In more detail, the BICs are spectrally surrounded by a leaky band of high{-}quality resonances (quasi-BICs) which can be excited
from the far-field~\cite{Yuan17}. The excitation of the strong resonances results in critical field enhancement~\cite{Yoon15, Mocella15a}
 with the near-field amplitude controlled by the frequency and the angle of incidence of the
incoming monochromatic wave.
The critical field enhancement triggers nonlinear optical effects even with a low intensity of the
incident light.
This resonant enhancement of nonlinear effects can lead to symmetry breaking~\cite{Bulgakov11},
channel dropping~\cite{Bulgakov13}, excitation of non-linear standing waves~\cite{Yuan16}
as well as self-adaptive robust~\cite{bulgakov2014robust} and tunable Fabry-Perot~\cite{Bulgakov2010} BICs.

In the field of nonlinear optics the BICs have been applied
for second harmonic  generation~\cite{ndangali2013resonant, koshelev2019meta, wang2018large}.
Quasi-BICs in subwavelength dielectric resonators~\cite{rybin2017high, bogdanov2019bound} are also shown
to demonstrate to enhance second harmonic generation~\cite{carletti2018giant, koshelev2020subwavelength}.
In~\cite{yuan2020excitation} it was found that, otherwise decoupled, BICs can be excited via second harmonic generation by illuminating
the structure from the far field.
At the same time it was shown theoretically that quasi-BICs allow for optical bistability~\cite{Yuan17,Krasikov18,Bulgakov19,maksimov2020optical}
due to the Kerr effect.

More recently, the research focus shifted towards BICs in lossy structures. In the presence of material absorption
the BIC are shown to acquire finite-life, albeit remain localized and decoupled from the outgoing channels~\cite{Hu20}.
Quasi-BIC in lossy periodic structures are found to be instrumental for enhancement of light absorption~\cite{Saadabad21, Tan22}
in the critical coupling regime even in low loss dielectrics. This opens novel opportunities for highly
efficient light absorbers~\cite{zhang2015ultrasensitive,wang2020controlling,sang2021highly,xiao2021engineering, cai2022enhancing}.

Lately, it has been suggested that thermo-optical effects can be the dominating nonlinear effects
in BIC supporting structures due to heating by absorbed radiation\cite{zograf2021all}.
In this work we investigate resonantly enhanced thermo-optical bistability~\cite{Sun10, Khandekar17, Gao17a, Pottier21, Ryabov22} in
a system supporting an optical BICs. Taking into account
radiative heat transfer as a major cooling mechanism we shall construct a non-linear model based on the temporal coupled mode theory (TCMT)~\cite{Fan03}
and theoretically demonstrate thermo-optic hysteresis.

\section{Scattering theory}

We consider an array of identical dielectric rods of radius  $R_0$ linearly arranged with
period $L$ in vacuum. The rods are made of amorphous silicon. The axes of the rods are collinear and aligned with the $z$-axis as shown in Fig.~\ref{fig1}~(a).
Such a system is known to support an abundance of BICs as demonstrated in~\cite{Bulgakov2015,Bulgakov17a}.
In this work we performed numerical simulations with application of FDTD Lumerical to examine the properties of the BIC induced optical response taking
into account both temperature and frequency dependence of the refractive index. In Fig.~\ref{fig1}~(b) we compare the eigenmode profile of an optical BIC with
vacuum wavelength $\lambda_{\scriptscriptstyle{\rm BIC}}=782~\rm{nm}$ against the scattering solution obtained under illumination by a TE plane wave with $\lambda=785~\rm{nm}$ at the incidence angle
$\theta=8.97~\mathrm{deg}$, and the incident wave vector on the $x0y$-plane as shown in Fig.~\ref{fig1}~(a). The results are obtained with the numerical values
of the refractive index from~\cite{Jellison94}.
In Fig.~\ref{fig1}~(b) one can see a striking similarity between the two field profiles.

By definition a BIC can not couple
to incident light. However, each BIC is a singular point on the dispersion sheet of a leaky band where the $Q$-factor diverges to infinity. The BIC is, thus,
spectrally surrounded by a family of high-$Q$ leaky modes (quasi-BICs) with the resonant $Q$-factor controlled by variation of the angle of incidence.
In our case the BIC occurs in the $\Gamma$-point being symmetrically mismatched from the single radiation channel of the zeroth diffraction order. Therefore,
the high-$Q$ resonant response is triggered by setting off the angle of incidence from zero. Thus, the similarity between
the BIC and the scattering solution is explained by both BIC and quasi-BIC sitting in the same dispersion band.

Our goal in this section is to set up the optimal regime for enhanced light absorption leading to the most significant change of the refractive index by heating.
The problem of enhanced absorption by quasi-BICs has been previously considered in the literature~\cite{Saadabad21}. The central result for upside-down symmetric
structures, such as the one in Fig.~\ref{fig1}~(a),  is that the maximal absorption of $50\%$ can be achieved in the critical coupling point
where the radiation and absorption loss rates are equal to one another.

\begin{figure}
\includegraphics[width=0.49\textwidth, height=0.25\textwidth,trim={0cm 0cm 0cm 0cm},clip]{fig1.pdf}
\caption{(a) Set-up of the array of dielectric rods. (b) Comparison between the eigenmode profile for a BIC with
$\lambda_{\scriptscriptstyle{\rm BIC}}=782~\rm{nm}$ and scattering solution under illumination by a plane wave with $\lambda=785~\rm{nm}$.
The solutions are visualized as the $z$-components of the electric field.
The geometric parameters: $R_0=128~{\rm{nm}}$, $L=428~\rm{nm}$.}
\label{fig1}
\end{figure}

The BIC induced optical response can be understood in the framework of temporal coupled mode theory (TCMT)~\cite{Fan03}. For the reader's convenience
we list the most important TCMT formulas below. Let us consider two-channel scattering for TE-polarized light.
The $S$-matrix is implicitly defined
through the following equation
\begin{equation}\label{Scat}
\left(
\begin{array}{c}
s_{1}^{\scs (\!-\!)} \\
s_{2}^{\scs (\!-\!)} \\
\end{array}
\right)
=\widehat{S}
\left(
\begin{array}{c}
s_{1}^{\scs (\!+\!)} \\
s_{2}^{\scs (\!+\!)} \\
\end{array}
\right),
\end{equation}
where $s_{m}^{\scs (\!\pm\!)}$ are the amplitudes of the plane waves
in the far-field with subscript $m=1,2$ corresponding
to the upper and lower half-spaces while superscripts $^{\scs (\!+\!)}, ^{\scs (\!-\!)}$ stand for incident and outgoing waves respectively. We assume that the system is illuminated
by a monochromatic wave of frequency $\omega$. In what follows we introduce the vectors of incident and outgoing amplitudes $|s^{\scs (\!\pm\!)}(t)\rangle$ which oscillate in time with the harmonic factor $e^{-i\omega t}$.
The TCMT equations~\cite{Fan03} take the following form
\begin{align}\label{CMT_linear}
& \frac{d a(t)}{d t}=-(i\omega_0+\gamma+\gamma_0)a(t)+\langle d^{*}|s^{\scs (\!+\!)}\rangle, \nonumber \\
&|s^{\scs (\!-\!)}\rangle=\widehat{C}|s^{\scs (\!+\!)}\rangle+
a(t)|d\rangle,
\end{align}
where $\widehat{C}$ is the matrix of direct (non-resonant) process, $\omega_0$ is the resonance center frequency,
$\gamma$ is the radiation decay rate, $\gamma_0$ is material loss decay rate, $a$ is the amplitude of the resonant eigenmode, and
$|d\rangle$ is the $2\times1$ vector of coupling constants satisfying
\begin{equation}\label{CMT1}
\langle d|d\rangle=2\gamma.
\end{equation}
The solution
for the $S$-matrix reads
\begin{equation}\label{S}
\widehat{S}=\widehat{C}+\frac{|d\rangle\langle d^*|}{i(\omega_0-\omega)+\gamma+\gamma_0}.
\end{equation}
In the case of the center-plane mirror symmetry
we have
\begin{equation}\label{direct}
\widehat{C}=
e^{i\phi}\left(
\begin{array}{cc}
\rho & i\tau\\
i\tau & \rho
\end{array}
\right),
\end{equation}
where $\phi, \rho, \tau$ are real valued parameters such as
\begin{equation}
\rho^2+\tau^2=1.
\end{equation}
The following equation also holds true~\cite{Fan03}
\begin{equation}\label{CMT2}
\widehat{C}|d^{*}\rangle=-|d\rangle.
\end{equation}
Using Eq.~(\ref{CMT1}) together with Eq.~(\ref{CMT2})  one finds
\begin{equation}
\label{d}
|d\rangle=
e^{i\frac{\phi}{2}}
\sqrt{\frac{\gamma }{2(1+\rho)}}
\left(
\begin{array}{l}
\tau-i(1+\rho) \\
\tau- i(1+\rho)
\end{array}
\right),
\end{equation}
for symmetric modes, and
\begin{equation}
\label{dm}
|d\rangle=
e^{i\frac{\phi}{2}}
\sqrt{\frac{\gamma }{2(1+\rho)}}
\left(
\begin{array}{l}
\tau+i(1+\rho) \\
-\tau- i(1+\rho)
\end{array}
\right),
\end{equation}
for anti-symmetric ones.
 The reflectance can be found from Eq.~(\ref{S}) as
\begin{equation}\label{R}
R=\frac{\rho^2(\omega-\omega_0)^2\mp 2\rho\tau\gamma(\omega-\omega_0)+\tau^2\gamma^2+\rho^2\gamma_0^2}{(\omega-\omega_0)^2+(\gamma+\gamma_0)^2},
\end{equation}
while the transmittance is given by
\begin{equation}\label{T}
T=\frac{\tau^2(\omega-\omega_0)^2\pm 2\rho\tau\gamma(\omega-\omega_0)+\rho^2\gamma^2+\tau^2\gamma_0^2}{(\omega-\omega_0)^2+(\gamma+\gamma_0)^2},
\end{equation}
where the top sign is used for symmetric modes.
Finally, the absorbance is given by
\begin{equation}\label{A}
A=\frac{2\gamma\gamma_0}{(\omega-\omega_0)^2+(\gamma+\gamma_0)^2}.
\end{equation}
The maximal absorbance is obtained at the critical coupling condition
$\gamma=\gamma_0$ and $\omega=\omega_0$. Substituting the above into Eq.~\eqref{A} one finds
\begin{equation}
A_{\mathrm{max}}=\frac{1}{2}.
\end{equation}

\begin{figure}
\includegraphics[width=0.49\textwidth, height=0.25\textwidth,trim={0cm 0cm 0cm 0cm},clip]{fig2.pdf}
\caption{Scattering spectra for the system shown in Fig.~\ref{fig1}~(a), $R_0=128~{\rm{nm}}$, $L=428~\rm{nm}$. (a) The reflectance, transmittance, and absorbance spectra found from
FDTD simulations at the critical coupling angle $\theta=8.97~\mathrm{deg}$. (b)
The TCMT approximation of the data in subplot (c); $L\omega/c=3.425$, $L\gamma/c=L\gamma_0/c=0.0056$, $\rho=0$.
}
\label{fig2}
\end{figure}
Since the system supports its antisymmetric BIC in the $\Gamma$-point, we expect that the resonant frequency $\omega_0$ and the radiation decay rate $\gamma$
are given by the following equations~\cite{Bulgakov17oe}
\begin{align}\label{Taylor1}
& \omega_0=\omega_{\scriptscriptstyle \rm {BIC}}+\kappa_{\omega}\theta^2+\mathcal{O}(\theta^4), \nonumber \\
& \gamma=\kappa_{\gamma}\theta^2+\mathcal{O}(\theta^4),
\end{align}
where $\kappa_{\omega,\gamma}$ are the coefficients of the Taylor expansion.
At the same time the absorption decay rate can be assessed in the perturbative manner~\cite{vainstein1988electromagnetic}  as
\begin{equation}\label{gamma0}
\gamma_0=\epsilon^{\prime\prime}\mkern-1.2mu \int\limits_{S_R} dS \frac{{\bf E}^{\dagger}{\bf E}+{\bf H}^{\dagger}{\bf H}}{4},
\end{equation}
where $\epsilon^{\prime\prime}\mkern-1.2mu$ and
 ${\bf E}$ and ${\bf H}$ are the electric and magnetic vectors of the BIC eigenfield
satisfying the normalization condition
\begin{equation}\label{normalization}
\int\limits_S dS \frac{\epsilon(x,y){\bf E}^{\dagger}{\bf E}+{\bf H}^{\dagger}{\bf H}}{4}=1,
\end{equation}
with $S$ as  the area of the elementary cell, and $S_R$ as the area of the rod cross section.
By a close examination of Eq.~(\ref{Taylor1}) and Eq.~(\ref{gamma0}) one concludes that the critical coupling
condition can always be achieved by increasing the angle of incidence. In our simulations we found that the critical coupling
condition is fulfilled at $\theta=8.97~\mathrm{deg}$ which is the case shown in Fig.~\ref{fig1}~(b). In Fig.~\ref{fig2}~(a) we show transmittance, reflectance,
and absorbance at the critical angle found via direct numerical simulations. In Fig.~\ref{fig2}~(b) we plotted the TCMT approximation with
Eqs.~(\ref{R},~\ref{T},~\ref{A}). The non-resonant reflection coefficient is taken $\rho=0$ according to~\cite{Bulgakov19}. The other approximation parameters
are listed in the caption to Fig.~\ref{fig2}. By comparing Fig.~\ref{fig2}~(a) against Fig.~\ref{fig2}~(b) one can see a certain deviation from the TCMT approach
which
is due to the perturbative nature of the solution given by Eqs.~(\ref{R},~\ref{T},~\ref{A}). In particular, the absorption in the direct process
is not accounted for by Eq.~\eqref{direct}. A better accuracy of the perturbative solution can be achieved in the near infrared where the imaginary part of
the refractive
index drops by two orders of magnitude~\cite{schinke2015uncertainty}. In Fig.~\ref{fig3} we compare the numerical data for dielectric array with
$L=580~\mathrm{nm}, R_0=
174~{\rm{nm}}$ supporting a BIC at wavelength $\lambda_{\scriptscriptstyle\rm{BIC}}=1027~\rm{nm}$. One can see an excellent coincidence between the TCMT and
numerical simulations. Also
notice that the position of the critical coupling spot in the absorption spectrum is shifted towards smaller angles of incidence due to a
smaller $\gamma_0$. Further
on in this work we stay in the visible relying on approximation shown in Fig.~\ref{fig2}~(b), and bearing in mind that the accuracy of the TCMT can be
improved by going
to the near infrared.
\begin{figure}[t]
\includegraphics[width=0.49\textwidth, height=0.63\textwidth,trim={0.0cm 0.0cm 0.0cm 0.0cm},clip]{fig3.pdf}
\caption{Absorbance, reflectance, and transmittance in the spectral vicinity of the BIC in the near infrared, $\lambda_{\scriptscriptstyle\rm{BIC}}=1027~\rm{nm}$, $L=580~\mathrm{nm}, R_0=
174~{\rm{nm}}$.
Column (a) shows the data obtained by direct FDTD simulations. Column (b) shows the TCMT approximation, $L\omega_{\scriptscriptstyle\mathrm{BIC}}/c=3.549$, $\rho=0$, $L\kappa_{\omega}/c=-0.00016$,
 $L\kappa_{\gamma}/c=0.000072$.  }
\label{fig3}
\end{figure}

\section{Heat transfer}

Since our system consists of dielectric rods in vacuum, the only mechanism of heat transfer is radiation from the rods.
The total energy lost per unit surface area across all wavelengths per unit time is given by the Stefan-Boltzmann law
\begin{equation}\label{SB}
    j=\alpha\sigma[(\Theta_0+\Delta \Theta)^4-\Theta_0^4],
\end{equation}
where $\alpha$ is the emissivity, $\sigma=
5.670\times 10^{-8}
~\mathrm{W}\cdot\mathrm{m}^{-2}\cdot\mathrm{K}^{-4}$ is the Stefan-Boltzmann constant, $\Theta_0$ is the temperature
of the surrounding medium and $\Delta \Theta$ is the temperature increase due to absorption of light. The emissivity, $\alpha$, is a quantity of key importance,
that describes the difference in the radiative properties between
the dielectric rod and the black body. In this work we assessed the emissivity by computing the
absorption coefficient of a single dielectric rod in the temperature range $25^{\circ}\mathrm{C}\leq\Theta\leq 300^{\circ}\mathrm{C}$.
The absorption coefficient was found as the average over all possible direction of the incident wave vector for both polarizations and, then,
averaged with the Planck's distribution
\begin{figure}
    \centering
    \includegraphics[width=0.30\textwidth, height=0.27\textwidth,trim={0cm 0cm 0.0cm 0cm},clip]{fig4.pdf}
    \caption{Emissivity of a single dielectric rod as a function of temperature; $L=428$~nm, $R_0=128$~nm. The average emissivity, $\tilde{\alpha}$, is shown
    by the horizontal line.}
    \label{fig4}
\end{figure}
\begin{equation}
    p(\omega)=\frac{\hbar \omega^3}{4 \pi^3 c^2}
    \frac{1}{e^{\hbar \omega/(k_\mathrm{B}\Theta )}-1}
\end{equation}
 peaked at about $10~\mu{\mathrm{m}}$. The numerical values of the refractive index were taken from~\cite{Chandler-Horowitz05}. The temperature dependance of
 the refractive index was neglected in the mid-infrared.
 The result is demonstrated in Fig.~\ref{fig4}. One can see that the emissivity
is dependant on temperature of the rods due to the shift of the peak of Planck's distribution. For further convenience we computed
the temperature average $\tilde{\alpha}=2.8\times10^{-5}$ which will be applied in our simulations. The average emissivity is shown in  Fig.~\ref{fig4} as a horizontal line.

The regime of heating can be determined by the value of the Biot number~\cite{bergman2011fundamentals}
\begin{equation}\label{Bi}
\mathrm{Bi}=\frac{h_S}{k}L_0,
\end{equation}
where $L_0$ is the characteristic length, $h_S$ is the surface heat transfer coefficient, and $k$ is the thermal conductivity of the body.
Values of the Biot number smaller than $0.1$ imply that the heat conduction inside the body is much faster than the heat lost
at the interface between the body and the surrounding medium. In this regime we can assume that the temperature is constant across the body.
By using Eq.~\eqref{SB} we find
\begin{equation}\label{transfer_coeff}
h_S=4\tilde{\alpha}\sigma\Theta_0^3\Delta\Theta.
\end{equation}
Taking $L_0$ as the perimeter of the rod, and $k=1.8~\mathrm{W}/(\mathrm{m}\cdot\mathrm{K})$ according to~\cite{Wada96} we find at the room temperature
$\mathrm{Bi}=6.5\cdot10^{-11}$. This allows us to assume that the temperature is uniform in the cross-section of the rod and avoid simulating heat transfer within the rods. Now our TCMT and heat transfer models can be linked by equating
the energy absorption rate in a unit length of the rod
\begin{equation}\label{energy1}
    \frac{dE}{dt}=\gamma_0(\Theta)|a|^2
\end{equation}
to the radiation rate from the surface per unit length of the same rod
\begin{equation}\label{energy2}
   \frac{dE}{dt}=2\pi R_0 \tilde{\alpha}\sigma[(\Theta_0+\Delta \Theta)^4-\Theta_0^4].
\end{equation}

\begin{figure}[t]
    \centering
    \includegraphics[width=0.37\textwidth, height=0.28\textwidth,trim={0cm 0cm 0.0cm 0cm},clip]{fig5.pdf}
    \caption{Real $n'$ and imaginary $n''$ parts of the refractive index of the amorphous silicon as a function of temperature at $\lambda=782~\mathrm{nm}$. Solid
    lines show linear approximation.}
    \label{fig5}
\end{figure}
\section{Non-linear TCMT equation and optical bistability}
\begin{figure}[ht]
\includegraphics[width=85mm]{fig6.pdf}
\caption{Thermo-optic hysteresis in the system shown in Fig.~\ref{fig1}~(a) with the same values of geometric parameters at incident angle $\theta=8.97^{\circ}$.
(a) Temperature increment against frequency of incident light, $I_0=0.54 \rm{W/m^2}$. (c) Scattering spectra, $I_0=0.54 \rm{W/m^2}$.
(b) Temperature increment against the intensity of incident light with vacuum wavelength $\lambda=795.5 \rm{nm}$. (d) Absorbance, reflectance, and transmittance
against the intensity, $\lambda=795.5 \rm{nm}$.}
\label{fig6}
\end{figure}
Frequency domain coupled mode equation reads
\begin{equation}\label{CMT_frequency}
[i\omega_0(\Theta)-i\omega+\gamma(\Theta)+\gamma_0(\Theta)]a=d_1(\Theta)\sqrt{LI_0},
\end{equation}
where $I_0$ is the intensity of the incident monochromatic light. Notice that according to Eq.~\eqref{normalization} the BIC eigenfield is
normalized to store a unit energy in the elementary cell per unit length of the rod. Thus, the intensity in Eq.~\eqref{CMT_frequency}
is multiplied by the period under the basic TCMT assumption
that the scattering channel function are normalized to supply (retract) a unit energy
in a unit of time~\cite{Fan03}.
All parameters including the resonant eigenfrequency, absorption and radiation decay rates are dependent on temperature. Hence we write
\begin{align}\label{Taylor2}
    & \omega_0(\Theta)=\omega_0^{\scriptscriptstyle{(0)}}
    +\omega_0^{\scriptscriptstyle{(1)}}\Delta\Theta, \nonumber \\
    & \gamma_0(\Theta)=\gamma_0^{\scriptscriptstyle{(0)}}
    +\gamma_0^{\scriptscriptstyle{(1)}}\Delta\Theta, \nonumber \\
    & \gamma(\Theta)_=\gamma^{\scriptscriptstyle{(0)}}
    +\gamma^{\scriptscriptstyle{(1)}}\Delta\Theta.
\end{align}
As it was mentioned the temperature dependence is due to the temperature dependance of the refractive index. By using
the experimental data from~\cite{Jellison94} we plotted the temperature dependance of both real $n'$ and imaginary $n''$ parts
of the refractive index of the amorphous silicon at $\lambda=782~\mathrm{nm}$ in Fig.~\ref{fig5}. One can see in Fig.~\ref{fig5} that
both quantities can be linearly approximated in the temperature range of interest.

By numerically computing the scattering spectra with
$\theta=8.97~\mathrm{deg}$ at elevated temperatures
and performing a TCMT fit we found the numerical values of the expansion coefficients in Eq.~\eqref{Taylor2} as follows
\begin{align}
& L\omega_0^{\scriptscriptstyle{(1)}}/c=-2.80\times 10^{-4},~\mathrm{K}^{-1}, \nonumber \\
& L\gamma_0^{\scriptscriptstyle{(1)}}/c=2.67\times10^{-5},~\mathrm{K}^{-1}, \nonumber \\
& L\gamma^{\scriptscriptstyle{(1)}}/c=6.67\times10^{-6},~\mathrm{K}^{-1}.
\end{align}
One can see from the above data that the dominating effect is a shift of the resonant frequency. The resonant frequency is at least one
order of magnitude more sensitive to heating than the other TCMT parameters. Also notice
that the absorption loss rate is much more sensitive to heating than the radiative loss rate. This can be understood from the explicit linear dependance
of $\gamma_0$ on the imaginary part of the refractive index given by Eq.~\eqref{gamma0}. In what follows we neglect
the temperature dependance of $\gamma$ and, consequently, of $|\bf{d}\rangle$.
From \eqref{CMT_frequency} and \eqref{CMT1} we find
\begin{equation}
    |a|^2=\frac{\gamma I_0L}{[\omega-\omega_0(\Theta)]^2
    +[\gamma+\gamma_0(\Theta)]^2},
\end{equation}
which in combination with Eq. \eqref{energy1} and Eq. \eqref{energy2} yields
\begin{align}\label{final}
   & \frac{ (\gamma_0^{\scriptscriptstyle{(0)}}
    +\gamma_0^{\scriptscriptstyle{(1)}}\Delta \Theta)\gamma I_0L}{(\omega-\omega_0^{\scriptscriptstyle{(0)}}
    -\omega_0^{\scriptscriptstyle{(1)}}\Delta \Theta)^2
    +(\gamma+\gamma_0^{\scriptscriptstyle{(0)}}
    +\gamma_0^{\scriptscriptstyle{(1)}}\Delta \Theta)^2}= \nonumber \\
    & 2\pi R_0 \tilde{\alpha}\sigma[(\Theta_0+\Delta \Theta)^4-\Theta_0^4].
\end{align}
The above formula constitutes a non-linear equation which has to be solved for $\Delta\Theta$.

By solving Eq.~\eqref{final} numerically in the frequency domain we find that a bistabily window of width $\Delta\lambda=1.6\rm{nm}$
exists at
\begin{equation}
\frac{2\pi R_0 \tilde{\alpha}\sigma}{I_0L}=6.15\cdot 10^{-12}
\end{equation}
at the incident wavelength $\lambda=795.9~\rm{nm}$. This in combination with our assessment for the emissivity $\tilde{\alpha}=2.80\cdot 10^{-5}$
yields the incident intensity $I_0=0.54~\rm{W/m^2}$. Our findings are illustrated in Fig.~\ref{fig6}. Figure~\ref{fig6}~(a), and Fig.~\ref{fig6}~(b)
demonstrate thermo-optical hysteresis in the frequency domain. One can see that at the peak of the resonance the system is heated up to $\Theta=187~^{\circ}\rm{C}$.
Figure~\ref{fig6}~(c), and Fig.~\ref{fig6}~(d)
demonstrate thermo-optical hysteresis in the intensity domain with the width of bistability window $\Delta I_0=1.6~\rm{W/m^2}$.
\section{Summary and Conclusions}

We considered thermo-optic hysteresis in a system supporting a symmetry protected bound state in the continuum. It is found that
bistability occurs at surprisingly low intensities of incident light, about $I_0=0.54~\rm{W/m^2}$. Two aspects are of key importance for
obtaining low bistability thresholds. (i) Critical coupling which is obtained by breaking the symmetry of the system by incident light.
With the increase of the angle of incidence the radiation decay rate grows quadratically until it becomes equal to the absorption decay rate,
$\gamma=\gamma_0$. In this situation, no matter how small the absorption decay rate is, $50\%$ of the incident light is absorbed and transformed to heat.
This approach paves a way for engineering effective absorbers with low loss dielectrics. (ii) Low emissivity which is due to a small optical path within the absorbing
dielectric at the peak wavelength of Planck's distribution at about $10~\rm{\mu m}$. The dielectric array is assembled of rods with a subwavelength cross section
in the visible. Thus, the optical path in the mid-infrared is at least an order of magnitude less than the wavelength of the thermal radiation. Given that the imaginary part
of the refractive index of silicon is vanishingly small in the mid-infrared~\cite{Chandler-Horowitz05}, $n''\approx 10^{-4}$, this results in a small value of
the emissivity $\tilde{\alpha}=
2.8\cdot 10^{-5}$.

The bistability occurs with intensities eleven orders of magnitude smaller than those needed for hysteresis due to Kerr nonlinearity in silicon
\cite{maksimov2020optical}. Thus, that thermo-optic effects should dominate in non-linear response of silicon nano-structures supporting bound states in
the continuum. In this work we have used a rough estimate for the emissivity. Yet, it is worth mentioning that even with $\tilde{\alpha}=1$, i.e. with
a mid-infrared
black body, the bistability thresholds would drop by less than five orders of magnitude and still be dominating in the non-linear response. Here we neglected
diffusive and convective heat transfer. In particular, the system under consideration has no substrate, which is necessary for realistic dielectric
metasurfaces. One would expect that the presence of substrate would further decrease the bistability thresholds as the heat would be leaking into the substrate.
This issue will be addressed in future studies.

\section*{Acknowledgments} This work received financial support through the grant of Russian Science Foundation and Krasnoyarsk
Regional Fund of Science No~22-22-20056, https://rscf.ru/project/22-22-20056/,
the authors acknowledge discussions with A.A. Bogdanov.


\section{Introduction}
Many systems have been proposed to improve the browsability and searchability of formalized mathematical libraries. The MathWiki Project\footnote{{https://www.nwo.nl/en/projects/612066825}} aims to improve the readability of formalized mathematical libraries and make them accessible to wider communities. A Wiki for Mizar~\cite{urban2010wiki}, Large Formal Wikis for Coq/CoRN~\cite{alama2011large}, and Agora System for Flyspeck Project~\cite{tankink2013formal} were proposed in the MathWiki project. However, since the MathWiki Project was terminated in 2014, its contents cannot follow the libraries' updates. ProofWiki~\cite{vyskocil2018disambiguating} and the Lean Mathematical Library~\cite{van2020maintaining} accumulate mathematical libraries and convert them into highly readable HTML documents. However, these two systems do not collaborate with advanced search engines or graphical tools to visualize library dependencies.

We developed the emwiki system, a Wiki service that hosts the Mizar Mathematics Library (MML)~\cite{bancerek2018role} while addressing issues with existing systems. The emwiki system is a Web platform based on the Django framework\footnote{https://www.djangoproject.com/}, featuring extensibility and interoperability. Vue.js\footnote{https://vuejs.org/} is used as the front-end framework, making it a partial single-page application. This eliminates extra rendering to reduce user stress. Vuetify\footnote{https://vuetifyjs.com/} is used as the user interface library, which allows for intuitive operation. Currently, this service is available at {https://em1.cs.shinshu-u.ac.jp/emwiki/release/}.

\section{Wiki Function}
The Wiki feature is implemented to embed additional comments to read and understand articles in the MML.
Although users cannot edit mathematical statements written in the formal language itself, they can add comments to theorems and definitions.
We reused the HTMLized MML~\cite{urban2013atp} to improve the convenience of the MML.
The HTMLized MML is a document in which the MML is converted to HTML format, and the reference relationships are expressed as hyperlinks.
We adopted the LaTeX syntax for mathematical expressions in the comments, which is familiar to mathematicians.
MathJax~\cite{cervone2012mathjax} is used as a mathematical expression rendering engine.
A real-time preview feature is also implemented to check rendering results while editing.

The version tracking feature manages the history of comments.
This function is achieved by using Git, a distributed version control system, as its backend.
Moreover, it is necessary to maintain the linkage between theorems/definitions and comments during library updates.
Our system exploits the Git merge function for this purpose.
As theorems/definitions do not have persistent identifiers in the Mizar language, the most effective way to identify theorems/definitions before and after a library update is to track text differences.
Our system embeds comments directly into the MML.
Since the comments are written just before the theorems/definitions, administrators can maintain the consistency of the linkage between the theorems/definitions and comments using the three-way merge function during library updates.

The user management function allows administrators to track and block malicious users.
The Wiki function stores a history of editors and their revisions, and when a comment is accidentally rewritten, it is possible to contact the user who wrote it and roll back it.
Also, if a malicious comment is found, it will be deleted, and the user who wrote it will be blocked.
The user management function is available in any component on our platform.

The screenshot of the Wiki feature is shown in Fig. \ref{fig:emwiki}.
\begin{figure}[hbt!]
        \includegraphics[width=\textwidth]{figures/emwiki.png}
    \caption{emwiki} \label{fig:emwiki}
\end{figure}

\section{Search Function}
The emwiki system has three types of search components.
The two incremental search components accept the name of articles and symbols as input.
As of 2018, the MML contains 1,290 articles and 8,852 symbols~\cite{bancerek2018role}, and being able to search for these articles and symbols efficiently is critical for users.
MML Reference~\cite{nakasho2015documentation} generates HTML documents from the MML to help users understand its symbols.
Each HTML page contains symbol definitions, referrers, and references.
It also has an incremental search function that enables users to search with a fast and intuitive operation.
We integrated MML Reference into emwiki and linked this search function with the Wiki function.

We also developed a flexible theorem search component.
For a long time, grep has been used for full-text search in the MML~\cite{matuszewski2005mizar}.
MML Query~\cite{bancerek2003information}, developed by G. Bancerek et al. in 2001, has dramatically improved the efficiency of the search of the MML.
However, MML Query is not easy for beginners, because it has its own query language.
MML Query is a pattern-matching-based search system.
Therefore, it is not good at searching modified (but logically equivalent) theorems.
The Alcor system~\cite{cairns2007integrating}, developed by P. Carins et al., provides an Latent Semantic Indexing (LSI)-based search function for the Mizar Mathematical Library.
We have implemented a similar LSI-based search component into the emwiki system.
The search component does not require its own query language, and the user can perform searches by inputting the desired theorems in the Mizar language.
We have also implemented a function to collect user evaluations by displaying a ``good'' button on the search results.
In the future, we plan to use the evaluations as training data for machine learning.

\section{Visualization of Library Dependencies}
The dependency graph component is designed to encourage library maintainers to refactor their libraries by visualizing the dependencies of files in the MML.
The MML has been maintained for the past 30 years by the University of Bialystok~\cite{grabowski2015four} and is continuously refactored using specialized tools to minimize the dependencies of the articles and eliminate the cyclical dependencies between groups~\cite{naumowicz2015tools}.
However, some of the Mizar language specifications hinder the refactoring of the MML~\cite{rudnicki2003integrity}.
The Mizar language does not have namespace or package features.
Also, the length of the Mizar file name is limited to 8 characters due to the legacy of MS-DOS constraints.
Therefore, all of the 1,290 articles with ambiguous file names are currently arranged in a flat structure in the MML.
Furthermore, the Mizar language has a constructor overload feature and a context-sensitive grammar in which the last imported constructor takes precedence.
These language specifications should be improved in the future, but we must continue to maintain our vast library of assets in the meantime.

J. Alama (2011)~\cite{alama2011mizar} created a database of reference relations among theorems, definitions, and notations in the MML and provided a way to access the dependency graph through a Web interface.
J. Heras et al. (2014)~\cite{heras2014hott} provides a tool to visualize the dependency graph of the HoTT library written in Coq.
J. Alama et al. (2012)~\cite{alama2012dependencies} extracted the dependencies in the Coq and Mizar libraries and performed a quantitative comparative analysis of the library features.
R. Marcus et al. (2020)~\cite{marcus2020tgview3d} proposed the TGView3D system, which renders the dependencies of formalized mathematical libraries as 3D graphs in a hybrid of force-directed and hierarchical layouts.

In this study, we adopted the hierarchical graph.
The hierarchical graph is intended to clarify the flow of dependencies and utilize refactoring, such as minimizing dependencies and eliminating circular references among groups.
In the Mizar language, reference relationships between articles are described in the environment part at the beginning of an article.
Since Mizar does not recursively load external files written in environment sections, a transitive reduction is performed to cut redundant edges before constructing the dependency graph.
The data structure is saved in dot and sfdp formats of Graphviz and is drawn in a Web browser using Cytoscape.js library.
The component also provides functions such as highlighting nearby nodes, searching nodes, moving nodes, hyperlinking to articles, and zooming in/out.
A detailed usage of the dependency graph component can be found in the Help menu in the upper right corner.

\section{Conclusions and Future Work}
In this paper, we have developed the emwiki system as a Web platform for hosting the MML.
This system is differentiated from existing services hosting mathematical libraries in consistency, extendability, and interoperability.
One of the most severe difficulties in the long-term use of a service that hosts formalized mathematical libraries is keeping up with library updates.
In this study, we attempted to solve this problem using Git and its merge function to maintain reference relationships during version updates.
Since this platform is built using Django, a new component development is stylized as adding Django applications.
Also, all components are allowed to access their version-controlled library, user management function, and database.
The components implemented on the platform share Web pages and data with each other.
For example, MML Reference, now integrated into the system, shares a Wiki function and Web pages.
Ease of operation is also improved through automated testing during development and deployment using containers.

For the dependency graph component, several issues remain.
While our hierarchical graph could properly place the dependencies of the articles, it did not classify the articles into meaningful groups.
Mizar articles have metadata such as the Mathematics Subject Classification (MSC)\footnote{{https://msc2020.org/}}.
We consider taking this information into account for graph layout and including it as additional information for article search.
We also think it will be effective to analyze the vocabulary contained in the articles and reflect the semantic classification in the layout algorithm.
Another beneficial enhancement is drawing theorem graphs.
However, handling as many nodes as the theorems in the MML requires a more powerful graph database and drawing library than the conservative library we currently use.

This project aims to improve the browsability, searchability, and comprehensive understandability of formalized mathematical libraries.
Although the system currently supports only the Mizar language, it could be extended to other languages.
It is also significant to incorporate other components, including a Web IDE.

It is also important to establish a cooperative framework for the long-term operation of the system. In particular, we need to seek the cooperation of the core team that maintains the MML and the Mizar system to create a flow to ensure consistency between the MML versions. Also, we would like to invite a wide range of developers from the Mizar community to work with us in defining requirements and developing new functions.

\section*{Acknowledgments}
This work was supported by JSPS KAKENHI Grant Number JP20K19863.

\bibliographystyle{splncs04}