\section{Introduction }
When discussing quantum field theory in the eternal Schwarzschild black hole background, one usually considers the Boulware, Unruh, or Hartle-Hawking states \cite{Boulware:1974dm}, \cite{Unruh:1976db}, \cite{Hartle:1976tp}, \cite{DeWitt:2003pm}, \cite{Candelas:1980zt}. While the Boulware state is defined by taking the exact modes to be positive--definite w.r.t. $\partial / \partial t$ Killing vector, which would give the usual Poincare--invariant ground state in flat--spacetime QFT, the Unruh and Hartle-Hawking states are defined by introducing Kruskal coordinates (canonical affine parameters on the past and future horizons) and taking modes to be positive--definite w.r.t. these Kruskal coordinates. The outcome of the two latter definitions of positive-definite modes is that either both or one of the out--going or in--going modes are in a thermal state with the temperature that is equal to the Hawking one. Which state is more suitable to a certain physical situation is analyzed by the calculation of the energy-momentum tensor on the horizon and at spatial infinity \cite{Candelas:1980zt}, \cite{Akhmedov:2015xwa}. However, there was no consideration of the situation in which both out--going and in--going modes are in thermal states with temperatures which generically do not coincide and which are not equal to the Hawking temperature. It would be natural to assume that this state is the best approximation for a physical situation when one considers a black hole in a box of a gas with a temperature which is not equal to the Hawking one. The same question can be posed analogously in the Rindler and de Sitter space--times with respective canonical temperatures of the horizons.

The de Sitter space, Rindler space, and two-dimensional black hole examples were analysed in \cite{Akhmedov:2020ryq}. It was shown that when one considers fields in a thermal state with temperature which is not equal to the canonical temperature of the horizon of the background geometry, then propagators, both points of which are located at the horizon, acquire an anomalous singularity. This singularity is anomalous in a sense that the standard light--like separation divergence comes from high frequencies and does not depend on the temperature of the state, while at the horizon the singularity comes from the infrared region in the integral over frequencies and there is an explicit dependence on the temperature of the state in which the propagator is evaluated.

We consider the following quantum field theory
\begin{align*}
    S=\frac{1}{2}\int d^4 x \sqrt{-g}\Big[ \partial_\mu \phi \partial^\mu \phi -\mu^2 \phi^2\Big]
\end{align*}
on the Schwarzschild and Reissner-Nordström black hole backgrounds.

In \cite{Akhmedov:2020ryq} only the two--dimensional analogue of the Schwarzschild black hole was considered. In this paper these results are extended to the four--dimensional black holes. We start in the section \ref{modes_wight_sch} by constructing modes and the Wightman function in the Schwarzschild black hole background. In the section \ref{canonical_sing} it is shown that the singularity of the Wightman function is canonical when its points are light-like separated outside the horizon, i.e. the divergent term in the Wightman function has the standard coefficient. Then, in the section \ref{anom_sing} we demonstrate that the Wightman function acquires anomalous singularity when both of its points are located on the horizon, i.e. the coefficient before the divergent term differs from the canonical one and explicitly depends on the temperature of the out-going modes. Finally, we conclude in section \ref{RNsection} by showing that the anomalous singularity of the Wightman function also occurs in the Reissner-Nordström black hole background.

\section{Schwarzschild black hole}

In this section we consider both massive and massless fields in Schwarzschild spacetime. First, we set the notations and define modes and thermal Wightman functions. For more detailed review the reader may refer to \cite{Boulware:1974dm}, \cite{Unruh:1976db}, \cite{Hartle:1976tp}, \cite{DeWitt:2003pm}, \cite{Candelas:1980zt}. Then, the Wightman function is computed when both of its points are located at the event horizon.

\subsection{Modes and Wightman function} \label{modes_wight_sch}
\begin{figure}[!h]
\centering
\includestandalone[width=0.7\textwidth]{picSCHW}
\caption{Penrose diagram of the Schwarzschild black hole. Here and further "PNI" and "FNI" stand for "past null infinity" and "future null infinity" correspondingly.} \label{picschw}
\end{figure}

The Schwarzschild black hole background is given by the metric:
\begin{align}
    &ds^2 = \bigg( 1-\frac{2M}{r} \bigg) dt^2 - \frac{dr^2}{1-\frac{2M}{r}}-r^2(d \theta^2 + \sin^2 \theta d \varphi^2).
\end{align}
Here Eddington-Finkelstein coordinates are defined as:
\begin{align}
    &u = t- r_*, \qquad v=t+r_*, \nonumber \\
    &r_* = r+2M \log \big( r/2M-1 \big).
\end{align}
In Sec. \ref{canonical_sing} it will be shown that the Wightman function has a standard divergence when its points are lightlike separated near the horizon. Then, in Sec. \ref{anom_sing} I will consider the situation in which both points of the Wightman function are located on the future horizon, which is denoted as $H_{future}$ in Fig. \ref{picschw}.

\subsubsection{Massless case}
Solutions of the corresponding d'Alembertian equation are as follows:
\begin{equation} \label{split}
    u_{\omega lm} (x) = (4 \pi |\omega|)^{-1/2} r^{-1} f_{\omega l} (r,t) Y_{lm} (\Omega),
\end{equation}
where $\Omega$ denotes the angular coordinates.

To define the set of Boulware modes, we impose the condition that the modes are positive-definite w.r.t. ${\partial}/{\partial t}$ Killing vector, so one can represent the function $f_{\omega l} (r,t)$ in the form
\begin{equation} \label{f_split}
    {f}_{\omega l} (t,r) = e^{-i|\omega| t} {R}_l (\omega | r) ,
\end{equation}
in which the radial function $R_l (\omega | r)$ solves the equation
\begin{align} \label{radial}
    &\frac{d^2R}{dr_*^2}+\big[ \omega^2 - V_l (r) \big] R = 0,
\end{align}
with the potential
\begin{equation} \label{potential}
    V_l (r) = \bigg( 1 - \frac{2M}{r} \bigg) \bigg(  \frac{l(l+1)}{r^2} + \frac{2M}{r^3} \bigg).
\end{equation}
The complete set of modes in the outer region of the Schwarzschild black hole is given by the out-going and in-going modes, denoted accordingly as:
\begin{eqnarray} \label{set1}
&\overrightarrow{u}_{\omega lm} (x) = (4 \pi |\omega|)^{-1/2} e^{-i |\omega| t} r^{-1} \overrightarrow{R}_l(\omega | r) Y_{lm} (\theta, \varphi), \nonumber \\
&\overleftarrow{u}_{\omega lm} (x) = (4 \pi |\omega|)^{-1/2} e^{-i |\omega| t} r^{-1} \overleftarrow{R}_l(\omega | r) Y_{lm} (\theta, \varphi),
\end{eqnarray}
where the radial functions satisfy the boundary conditions as follows
\begin{equation}
                \overrightarrow{R}_l(\omega| r) = \begin{cases}
                 e^{i \omega r_*} + \overrightarrow{A}_l(\omega) e^{-i \omega r_*},
                & r \rightarrow 2M \\
                {B}_l (\omega) e^{i \omega r_*}, & r \rightarrow \infty
                \end{cases} \nonumber
            \end{equation}
    and
    \begin{equation} \label{boundary_cond}
                \overleftarrow{R}_l(\omega |r) = \begin{cases}
                {B}_l ({\omega}) e^{ - i \omega r_*}, & r \rightarrow 2M \\
                e^{-i \omega r_*} + \overleftarrow{A}_l ({\omega}) e^{i \omega r_*}, & r \rightarrow \infty
                \end{cases}
            \end{equation}
As we have only two types of modes, the mode decomposition of the field operator is as follows
\begin{equation}
    \phi (x) = \sum_{lm} \int_0^{\infty} d \omega \bigg( a_{\omega lm} \overrightarrow{u}_{\omega lm} (x) + b_{\omega lm} \overleftarrow{u}_{\omega lm} (x) + h.c. \bigg).
\end{equation}
Then for the Wightman function in a generic state with zero anomalous quantum averages (i.e. with $\langle a_{p} a_{q}\rangle = \langle a^{\dagger}_{p} a^{\dagger}_{q}\rangle =\langle b_{p} b_{q}\rangle = \langle b^{\dagger}_{p} b^{\dagger}_{q}\rangle = 0 $), but which is not necessary Fock ground state, one has:

\begin{multline} \label{wightman_general}
    W (x,x') \equiv \langle \phi (x) \phi(x') \rangle
= \int^{+\infty}_{0} {d \omega} \int^{+\infty}_{0} {d \omega'} \bigg[ \langle a_{\omega} a_{\omega'}^{\dagger} \rangle \overrightarrow{u}_{\omega}(x) \overrightarrow{u}^*_{\omega'}(x')
 + \langle a_{\omega}^{\dagger} a_{\omega'} \rangle \overrightarrow{u}^*_{\omega}(x) \overrightarrow{u}_{\omega'}(x')
 + \\
 +\langle b_{\omega} b_{\omega'}^{\dagger} \rangle\overleftarrow{u}_{\omega}(x) \overleftarrow{u}^*_{\omega'}(x') + \langle b_{\omega}^{\dagger} b_{\omega'} \rangle
  \overleftarrow{u}^*_{\omega}(x) \overleftarrow{u}_{\omega'}(x') \bigg].
\end{multline}
Here we have omitted spherical harmonic indices to simplify the equation.
Allowing the temperatures of out--going and in--going modes to differ from each other, in this paper we will study states in which
\begin{equation}
     \langle a_{\omega}^{\dagger} a_{\omega'} \rangle = \frac{ \delta(\omega -\omega')}{e^{\beta_R |\omega|}-1}, \qquad \langle b_{\omega}^{\dagger} b_{\omega'} \rangle = \frac{\delta(\omega -\omega')}{e^{\beta_L |\omega|}-1},
\end{equation}
so the Wightman function can be written as
\begin{multline}
    W(x,x') = \sum_{lm} \int^{+\infty}_{0} {d \omega} \bigg[  \frac{\overrightarrow{u}_{\omega lm}(x) \overrightarrow{u}^*_{\omega lm}(x')}{1-e^{-\beta_R |\omega|}}
 + \frac{\overrightarrow{u}^*_{\omega lm}(x) \overrightarrow{u}_{\omega lm}(x')}{e^{\beta_R |\omega|}-1} +
  \frac{\overleftarrow{u}_{\omega lm}(x) \overleftarrow{u}^*_{\omega lm}(x')}{1-e^{-\beta_L |\omega|}}+ \\ + \frac{\overleftarrow{u}^*_{\omega lm}(x) \overleftarrow{u}_{\omega lm}(x')}{e^{\beta_L |\omega|}-1} \bigg] = \\
  =\sum_{lm} \int^{+\infty}_{0} \frac{d \omega}{4\pi |\omega|} \frac{1}{rr'}\bigg[ \frac{e^{i |\omega| (t-t')} \overrightarrow{R}^*_l (\omega|r) \overrightarrow{R}_l (\omega|r') }{e^{\beta_R | \omega|}-1}
  +\frac{e^{i |\omega| (t-t')} \overleftarrow{R}^*_l (\omega|r) \overleftarrow{R}_l (\omega|r') }{e^{\beta_L |\omega|}-1} \bigg] Y_{lm} (\Omega) Y^*_{lm} (\Omega') + \\+
  \sum_{lm} \int^{0}_{-\infty} \frac{d \omega}{4\pi |\omega|} \frac{1}{rr'}\bigg[ \frac{e^{-i |\omega| (t-t')} \overrightarrow{R}^*_l (\omega|r) \overrightarrow{R}_l (\omega|r') }{1-e^{-\beta_R | \omega|}}
  +\frac{e^{-i |\omega| (t-t')} \overleftarrow{R}^*_l (\omega|r) \overleftarrow{R}_l (\omega|r') }{1-e^{-\beta_L |\omega|}} \bigg] Y_{lm} (\Omega) Y^*_{lm} (\Omega'),
\end{multline}
where in the last relation we have made the substitution $\omega \to -\omega$ for the first and the third terms in the first line. Also we have used the relation $\overrightarrow{R}_l(-\omega|r) = \overrightarrow{R}^*_l(\omega|r)$, which can be deduced from \eqref{near_horizon_mode} or \eqref{pt_reflection} with \eqref{boundary_cond}. Then we put $|\omega| = +\omega$ in the first integral of the last equality and $|\omega| = -\omega$ in the second integral to obtain the following form of the Wightman function:
\begin{multline} \label{wightman}
    W(x,x') = \sum_{lm} \int^{+\infty}_{-\infty} \frac{d \omega}{4\pi \omega} \frac{1}{rr'}\bigg[ \frac{e^{i \omega (t-t')} \overrightarrow{R}^*_l (\omega|r) \overrightarrow{R}_l (\omega|r') }{e^{\beta_R  \omega}-1}
  +\\+\frac{e^{i \omega (t-t')} \overleftarrow{R}^*_l (\omega|r) \overleftarrow{R}_l (\omega|r') }{e^{\beta_L \omega}-1} \bigg] Y_{lm} (\Omega) Y^*_{lm} (\Omega').
\end{multline}
For the case $\beta_R = \beta_L = \infty$ we obtain the Boulware state, while $\beta_R=\frac{2\pi}{\kappa}, \beta_L = \infty$ corresponds to the Unruh state, and $\beta_R=\beta_L=\frac{2\pi}{\kappa}$ results in the Hartle-Hawking state. Note that in the limit $\beta=\infty$ Bose-Einstein distribution reduces to Heaviside theta-function $\theta(-\omega)$ and the Wightman function reduces to what one would have obtained from \eqref{wightman_general} by taking the state to be the Fock vacuum. Here $\kappa = (4M)^{-1}$ is the surface gravity of the black hole.
\subsubsection{Massive case}

In the case of massive field ($\mu \neq 0$), the equation for the radial function is
\begin{align} \label{radial}
    &\frac{d^2R}{dr_*^2}+\big[ \omega^2 - V_l (r) \big] R = 0,
\end{align}
with the potential
\begin{equation} \label{potential}
    V_l (r) = \bigg( 1 - \frac{2M}{r} \bigg) \bigg(  \frac{l(l+1)}{r^2} + \frac{2M}{r^3} + \mu^2 \bigg),
\end{equation}
which tends to $\mu^2$ in the $r^* \to \infty$ limit. So as in \cite{Akhmedov:2020ryq}, \cite{Akhmedov:2015xwa}, \cite{Akhmedov:2016uha} there are modes with $\omega^2 < \mu^2$ that are localised near the horizon and which are exponentially decreasing in the classically inaccessible domain.
Rewriting the equation \eqref{radial} in terms of coordinate $r$, taking the near-horizon limit and introducing the new coordinate $\xi^2 = \frac{r}{2M}-1$ leads to the equation
\begin{equation}
    \frac{d^2 R}{d\xi^2} + \frac{1}{\xi} \frac{dR}{d\xi} + \bigg[ \frac{(4M\omega)^2}{\xi^2} -4l(l+1)-(4M\mu)^2 \bigg] R = 0,
\end{equation}
solutions of which are given by the modified Bessel functions of the second kind. Keeping in mind the condition of exponential decrease in the classically inaccessible domain, one obtains
\begin{equation} \label{radial_massive_sch}
    R \approx  \sqrt{\frac{\omega \sinh({4\pi M \omega})}{\pi M}} K_{4iM\omega} \big( 2 \sqrt{(2M\mu)^2+l^2+l} \, \xi \big).
\end{equation}
Then, the modes are denoted as
\begin{align}
    &\varphi_{\omega lm} (x) = \frac{1}{\sqrt{ \pi |\omega|}} e^{-i |\omega| t} R_{ l} (\omega | r) Y_{lm} (\Omega),
\end{align}
where the radial part is given by \eqref{radial_massive_sch}, so that near the horizon one has
\begin{gather} \label{delta_sch}
    R_{l} (\omega | r) \approx \frac{1}{r} \cos (\omega r_* + \delta_{\omega l}), \\
    \delta_{\omega l} \approx \frac{\pi}{2} +2M\omega \log (4M^2 \mu^2 +l^2+l) - 2M\omega - \text{arg} \, \Gamma(1+4iM\omega).
\end{gather}
The modes with $\omega^2>\mu^2$ are similar to the massless case:
\begin{align} \label{f_modes}
    &\overrightarrow{F}_{\omega lm} (x) = (4 \pi |\omega|)^{-1/2} e^{-i |\omega| t} r^{-1} \overrightarrow{F}_l(\omega | r) Y_{lm} (\Omega), \nonumber \\
&\overleftarrow{F}_{\omega lm} (x) = (4 \pi |\omega|)^{-1/2} e^{-i |\omega| t} r^{-1} \overleftarrow{F}_l(\omega | r) Y_{lm} (\Omega),
\end{align}
where their radial parts obey the following boundary conditions:
\begin{align}
                &\overrightarrow{F}_l ({\omega} | r) = \begin{cases}
                e^{i \omega r_*} + \overrightarrow{C}_{\omega} e^{-i \omega r_*},
                & \qquad r \rightarrow 2M \\
                \sqrt{\frac{\omega}{p}} {D}_{\omega} e^{i p r_*}, & \qquad r \rightarrow \infty
                \end{cases} \\
                &\overleftarrow{F}_l ({\omega} |r) = \begin{cases}
                {D}_{\omega} e^{ - i \omega r_*}, & r \rightarrow 2M \\
                \sqrt{\frac{\omega}{p}} [ e^{-ipr_*} + \overleftarrow{C}_{\omega} e^{ipr_*}], & r \rightarrow \infty
                \end{cases}
            \end{align}
            with $p=\text{sgn} (\omega) \cdot \sqrt{\omega^2 - \mu^2}$, \, $\omega^2 > \mu^2$.

The mode decomposition of the field operator for the massive case is as follows
\begin{multline} \label{mode_massive_sch}
    \phi(x) = \sum_{lm} \int_0^{\mu} d\omega \big( \alpha_{\omega lm} \varphi_{\omega lm} + \alpha_{\omega lm}^{\dagger} \varphi_{\omega lm}^* \big) + \\
    +\sum_{lm} \int_{\mu}^{\infty} d \omega \big( \beta_{\omega lm} \overrightarrow{F}_{\omega lm} + \beta_{\omega lm}^{\dagger} \overrightarrow{F}^*_{\omega lm} + \gamma_{\omega lm} \overleftarrow{F}_{\omega lm} + \gamma_{\omega lm}^{\dagger} \overleftarrow{F}_{\omega lm}^* \big),
\end{multline}
and the Wightman function has the following form:
\begin{multline}
    W(x,x') = \sum_{lm} \int^{+\mu}_{-\mu} \frac{d \omega}{\pi \omega}  \frac{e^{i \omega (t-t')} R^*_l(\omega|r) R_l (\omega|r') }{e^{\beta_0 \omega}-1} Y_{lm}(\Omega) Y^*_{lm}(\Omega') + \\
    +\sum_{lm} \int_{|\omega|>\mu} \frac{d \omega}{4\pi \omega} \frac{1}{rr'} \bigg[ \frac{ e^{i \omega (t-t')} \overrightarrow{F}^*_{ l}(\omega | r) \overrightarrow{F}_{l}(\omega|r')}{e^{\beta_R \omega}-1}
 +
  \frac{ e^{i \omega (t-t')} \overleftarrow{F}^*_{l}(\omega|r) \overleftarrow{F}_{l}(\omega|r')}{e^{\beta_L \omega}-1} \bigg] Y_{lm}(\Omega) Y^*_{lm}(\Omega') .
\end{multline}
Note that we have added here the inverse temperature $\beta_0$ for the modes with $\omega^2 \leq \mu^2$, which does not have to be equal neither to $\beta_R$ nor to $\beta_L$.

\subsection{Canonical singularity} \label{canonical_sing}
First we want to explicitly show in what sense the divergences derived further in this section are anomalous. To do so, we consider the Wightman function \eqref{wightman} when both of it's points are located in the vicinity of the horizon, but not exactly on it. We will take light-like separation of these two points, keeping their angular coordinates equal. We restrict ourselves to the massless fields in this subsection.

The metric near the horizon can be written as
\begin{equation}
    ds^2 \approx \rho^2 d \eta^2 - d \rho^2 -r^2(\rho) \, d\Omega^2,
\end{equation}
where
\begin{align}
    \rho = 4M \sqrt{\frac{r}{2M}-1} ,\qquad \eta = \frac{t}{4M}.
\end{align}
Then, the square of the geodesic distance between the two points under consideration is written as
\begin{equation} \label{geodesic}
    L^2 \approx 2 \rho_1 \rho_2 \cosh (\eta_1-\eta_2) - \rho_1^2-\rho_2^2.
\end{equation}
 One can expand the potential near the horizon and find a solution (in terms of modified Bessel function of the second kind) that satisfies boundary condition \eqref{boundary_cond}:

\begin{equation} \label{near_horizon_mode}
    \overrightarrow{R}_l (\omega|r) \approx \frac{2}{\Gamma (-i 4 M \omega)} e^{-2iM \omega \log (l^2+l+1)+2iM\omega} K_{4iM\omega} (2 \sqrt{l^2+l+1} e^{\frac{r_*-2M}{4M}}),
\end{equation}
from which one can also find that
\begin{equation}
    \overrightarrow{A}_l (\omega) \approx - e^{4iM\omega - 4iM\omega \log(l^2+l+1) + 2i \text{arg} \Gamma(1+4iM \omega) }.
\end{equation}

Now, plugging \eqref{near_horizon_mode} into \eqref{wightman}, one obtains that
\begin{multline} \label{prop_non_horizon}
    W(x,x') \approx \frac{1}{\pi^2} \frac{1}{(2M)^2} \sum_{lm} Y_{lm} (\Omega) Y_{lm}^* (\Omega) \int^{\infty}_{-\infty} d\omega \frac{e^{\frac{i \omega (t-t')}{4M}}}{e^{\frac{\beta_R \omega}{4M}}-1} \sinh( \pi \omega) \cdot \\ \cdot
    K_{i\omega} (2 \sqrt{l^2+l+1} e^{\frac{r_*-2M}{4M}}) K_{i\omega} (2 \sqrt{l^2+l+1} e^{\frac{r_*'-2M}{4M}}).
\end{multline}
The simple formula is available if we take $\beta_R = \frac{8\pi M}{n}$ with $n$ being integer \cite{Akhmedov:2019esv}, \cite{Akhmedov:2020qxd}. Namely, we employ the fact that in this case
\begin{equation}
    \frac{\sinh(\pi \omega)}{e^{\frac{\beta_R \omega}{4M}}-1} = \frac{e^{-\pi \omega}}{2} \sum_{k=1}^{n} e^{\frac{2\pi \omega (k-1)}{n}}.
\end{equation}
Plugging this expressions into \eqref{prop_non_horizon} and evaluating the integral over $\omega$, we obtain that
\begin{multline} \label{wightman_discrete_sum}
    W(x,x') \approx \frac{1}{2 \pi} \frac{1}{(2M)^2} \sum_l \frac{(2l+1)}{4\pi} \times \\ \times
    \sum_{k=1}^{n} K_0 \bigg( 2 \Tilde{l} \sqrt{e^{\frac{r_*-2M}{2M}} + e^{\frac{r_*'-2M}{2M}} +2 e^{\frac{r_*-2M}{4M}} e^{\frac{r_*'-2M}{4M}} \cosh \bigg[ -\frac{t-t'}{4M} +i \pi - i \frac{2\pi (k-1)}{n}  \bigg] } \bigg),
\end{multline}
where $\Tilde{l} = \sqrt{l^2+l+1}$.
Here only the $k=1$ term depends on the geodesic distance between the two points \eqref{geodesic}, so detaching this terms from the sum and writing the Wightman function as:
\begin{equation} \label{wightman_k}
    W(x,x') = W_{k = 1} (x,x') + W_{k \neq 1} (x,x'),
\end{equation}
we obtain that for a light-like separated $x$ and $x'$, i.e. with $L=0$ (but not yet sitting at the horizon):
\begin{multline}
    W_{k=1}(x,x') \approx \frac{1}{2\pi} \frac{1}{(2M)^2} \sum_l \frac{(2l+1)}{4\pi}
    K_0 \bigg( \frac{\sqrt{l^2+l+1}}{2M} \sqrt{-{L^2}} \bigg) \approx \\
    \approx \frac{1}{4\pi^2} \frac{1}{(2M)^2} \int_0^{\infty} dl \, l K_0 \bigg( \frac{l}{2M} \sqrt{-L^2} \bigg) = -\frac{1}{4 \pi^2 L^2},
\end{multline}
and as the terms in $W_{k \neq 1} (x,x')$ do not depend on the geodesic distance between $x$ and $x'$, they are negligible when we are considering light-like separation of the points out of the horizon. Hence, we get that
\begin{equation}
    W(x,x') \approx  -\frac{1}{4 \pi^2 L^2},
\end{equation}
i.e. we have the standard coefficient of $(4\pi^2)^{-1}$ in front of the divergence.

However, when the points are located on the horizon, in \eqref{wightman_k} $k \neq 1$ terms also have divergent contributions. We can place the two points on the future horizon by taking the following limit:
\begin{align} \label{limit}
    &t=-\lambda, \qquad r_*=\lambda, \nonumber \\
    &t'=-\lambda, \qquad r_*'=\lambda+c, \qquad \lambda \to -\infty.
\end{align}
Now we demonstrate that in this limit the Wightman function acquires additional divergent terms. In the case when $n=2$, the Wightman function has the following form:
\begin{multline}
    W(x,x') \approx \frac{1}{2\pi} \frac{1}{(2M)^2} \sum_l \frac{(2l+1)}{4\pi} \bigg\{
    K_0 \bigg( \frac{\sqrt{l^2+l+1}}{2M} \sqrt{-{L^2}} \bigg) + \\ + K_0 \bigg( \frac{\sqrt{l^2+l+1}}{2M} \sqrt{-{L^2 \Lambda}} \bigg) \bigg\} \approx -\frac{1}{4 \pi^2 L^2} \frac{1+\Lambda}{\Lambda},
\end{multline}
where we have expressed the argument of the second term in the curly bracket through the geodesic distance in the limit \eqref{limit} via the coefficient $\Lambda$, which has the form
\begin{equation}
    \Lambda = \frac{1+e^{c/2M}+2 e^{c/4M}}{1+e^{c/2M}-2e^{c/4M}}.
\end{equation}
However, in this approach it is not clear how to obtain the answer for the Wightman function even in the case of the discrete inverse temperature $\beta_R = \frac{8\pi M}{n}$ with an arbitrary $n$, let alone the case of the arbitrary $\beta_R$. This difficulty comes from the fact that as opposed to the static patch of de Sitter space, in this case one cannot express the Wightman function with the non-canonical temperature as the sum of the canonical temperature Wightman functions. In the following section we calculate the Wightman function using the asymptotic form of the modes near the horizon, and the result is that the coefficient of the singularity of the Wightman function depends explicitly on the temperature (as it is shown in  \eqref{sing_schw_massive} and \eqref{schw_res} below). Furthermore, by taking high-frequency limit in \eqref{prop_non_horizon}, it can be shown that the singularity outside the horizon is an UV effect, while at the horizon the singularity comes from the IR region of the frequencies. This justifies the fact that everywhere outside the horizon the divergence does not depend on the temperature, while at the horzion there is an explicit dependence on the temperature (i.e. on the low laying state).

\subsection{Anomalous singularity} \label{anom_sing}

In this section the anomalous singularity in the Wightman function on the horizon for generic thermal state is derived.

\subsubsection{Massless case}

Using the asymptotic behavior of the modes \eqref{boundary_cond} in the limit \eqref{limit},
the Wightman function \eqref{wightman} can be written as:
\begin{multline} \label{prop}
    W(x,x') \approx \frac{1}{(2M)^2} \sum_{lm} Y^*_{lm} (\Omega) Y_{lm} (\Omega')  \int_{-\infty}^{\infty} \frac{d \omega }{4 \pi \omega} \bigg[ \frac{1}{e^{\beta_R \omega }-1} \bigg( e^{-i\omega c} + \overrightarrow{A}_l (\omega) e^{-2 i \lambda \omega } + \\ + \overrightarrow{A}_l^* (\omega) e^{2 i \lambda \omega } + | \overrightarrow{A}_l (\omega) |^2 e^{i\omega c} \bigg) + \frac{1}{e^{\beta_L \omega}-1} | B_l (\omega) |^2 e^{i \omega c} \bigg].
\end{multline}
To integrate over $\omega$, analytic expression for reflection coefficient $\overrightarrow{A}_l (\omega)$ is needed. To compute it, one needs to know transition coefficients between solutions of confluent Heun equation near the horizon and infinity. However, there is no knowledge about these coefficients \cite{Fiziev:2005ki}, but we can approximate potential for radial function by P{\"o}schl-Teller potential \cite{Chirenti:2017mwe},
\cite{Berti:2009kk}, which has the form
\begin{equation}
    \Tilde{V}_l(r) = \frac{V_0}{\cosh^2 \alpha (r_* - \Tilde{r}_*) },
\end{equation}
where $\Tilde{r}_*$ is the position of the peak and
\begin{align}
    V_0 = V_l (\Tilde{r}_*), \qquad \alpha^2 = - \frac{1}{2V_0} \frac{d^2 V_l}{dr_*^2}_{|r_*=\Tilde{r}_*}.
\end{align}
In this potential we can easily compute the reflection coefficient $\overrightarrow{A}_l (\omega)$, which is found to be
\begin{equation} \label{pt_reflection}
    \overrightarrow{A}_l (\omega) = \frac{\Gamma({i \omega}/{\alpha}) \Gamma(-i {\omega}/{\alpha}-s) \Gamma(1+s-i {\omega}/{\alpha})}{\Gamma(-i {\omega}/{\alpha}) \Gamma(1+s) \Gamma(-s) }, \qquad s = -\frac{1}{2} + \frac{i}{2} \sqrt{\frac{4V_0}{\alpha^2}-1}.
\end{equation}
So we have to evaluate the integral
\begin{multline}
    \int_{-\infty}^{\infty} \frac{d \omega }{4 \pi \omega} \frac{ \overrightarrow{A}_l (\omega) e^{-2 i \lambda \omega }}{e^{\beta_R \omega}-1} = \int_{-\infty}^{\infty} \frac{dz}{4 \pi z} \frac{1}{e^z-1}  \frac{\Gamma(iz/\alpha \beta_R) \Gamma(-iz/ \alpha \beta_R-s) \Gamma(1+s-iz/ \alpha \beta_R)}{\Gamma(1+s) \Gamma(-s) \Gamma(-iz/ \alpha \beta_R)} e^{2 i \frac{| \lambda |}{\beta_R} z} = \\
    = \pi i \text{Res}(f(z),0) + 2\pi i \text{Res}(f(z), \overline{0}),
\end{multline}
where the second term $\text{Res}(f(z), \overline{0})$ denotes the sum over all residues on the upper half-plane poles.

The residue at $z=0$ in the limit $\lambda \to -\infty$ is computed to be
\begin{equation}
    \pi i \text{Res}(f(z),0) =  \frac{|\lambda|}{2 \beta_R}.
\end{equation}
Here an additional term of the form $\Gamma'(0)/\Gamma(0)$ could arise from derivatives of gamma functions, but it is discarded in the limiting procedure as it is of the form
\begin{equation} \label{gammaterms}
    \frac{\Gamma(-iz / \alpha \beta_R) \Gamma'(iz/ \alpha \beta_R)}{\Gamma^2 (iz/ \alpha \beta_R)} + \frac{\Gamma(iz / \alpha \beta_R) \Gamma'(-iz/ \alpha \beta_R)}{\Gamma^2 (iz/ \alpha \beta_R)}.
\end{equation}
In fact, using that
\begin{equation}
    \Gamma(\epsilon) = - \Gamma(-\epsilon) \frac{\Gamma(1+\epsilon)}{\Gamma(1-\epsilon)}
\end{equation}
one concludes that \eqref{gammaterms} is zero.

Furthermore
\begin{equation}
    2\pi i \text{Res}(f(z), \overline{0})= -\frac{i}{2} \sum_{n>0} \frac{(-1)^n}{(n!)^2} \frac{1}{e^{i \alpha \beta_R n}-1} \frac{\Gamma(n-s) \Gamma(1+s+n)}{\Gamma(1+s) \Gamma(-s)} e^{-2 |\lambda| \alpha n},
\end{equation}
which tends to zero as $\lambda \to -\infty$.

So in the limit when both points of the propagator are located on the horizon, we obtain
\begin{equation}
    \int_{-\infty}^{\infty} \frac{d \omega }{4 \pi \omega} \frac{\overrightarrow{A}_l (\omega) e^{-2 i \lambda \omega }}{e^{\beta_R \omega}-1} \approx  \frac{|\lambda|}{2 \beta_R}, \qquad \lambda \to -\infty.
\end{equation}
The other complex conjugated term in \eqref{prop} gives the same contribution. Then, using that
\begin{align}
    P_l (\vec{x} \cdot \vec{y} ) = \frac{4 \pi}{2l+1} \sum_{m=-l}^l Y_{lm} (\vec{x}) Y_{lm}^*(\vec{y}), \qquad \sum_{l=0}^{\infty} \frac{2l+1}{2} P_l (x) P_l(y) = \delta(x-y),
\end{align}
the Wightman function can be rewritten as
\begin{equation} \label{schw_res}
    W(x,x') \approx  \frac{|\lambda|}{ \beta_R} \frac{\delta(\Omega,\Omega')}{(2M)^2} , \qquad \lambda \to - \infty.
\end{equation}
Actually, one can obtain the same answer without approximating the potential and employing exact solution in fitted potential. Namely, one can use the approximate near-horizon solution \eqref{near_horizon_mode},
but ultimately the result \eqref{schw_res} remains the same. However, the near-horizon approximation suffers from the fact that the potential for the radial function in this case is infinitely growing, so that one can only introduce the out-going modes that are fully reflected back from the potential, i.e. in this situation one has $|\overrightarrow{A}_l (\omega)|=1$. The approximation with the P{\"o}schl-Teller potential shows that deformation of the potential from the infinitely growing with no in-going modes to the bump-like with the presence of the transmitted in-going modes does not change the answer for the anomalous singularity.
This shows that this anomalous singularity is not a peculiarity for a given solution of Einstein's equations, but actually its presence does not depend on the form of the potential.

Also, one can see that for the discrete temperature $\beta = \frac{8 \pi M}{n}$, the result \eqref{schw_res} can be obtained from \eqref{wightman_discrete_sum} by taking the asymptotic form of the $K_0(x)$ as $x \to 0$. Dependence on $l$ splits off from the rest of the argument of this function because of the logarithm function in this asymptotic, so the leading terms in the horizon limit are given by the terms that depend on the geodesic distance, and the sum over $l$ gives the angular delta-function in \eqref{schw_res}.

\subsubsection{Massive case}

As the anomalous singularity on the horizon comes from the low-frequency modes (see previous subsection or \cite{Akhmedov:2020ryq}), we neglect the contribution from high-frequency modes in \eqref{mode_massive_sch}. Then,
the Wightman function in the horizon limit is as follows:
\begin{equation}
    W(x,x') \approx \sum_{lm} \int^{+\mu}_{-\mu} \frac{d \omega}{\pi \omega}  \frac{e^{i \omega (t-t')} R^*_l(\omega|r) R_l (\omega|r') }{e^{\beta_0 \omega}-1} Y_{lm}(\Omega) Y^*_{lm}(\Omega')
\end{equation}
Assuming that $\mu \beta_0$ is sufficiently large for the integrand to be zero at the endpoints of the integration domain, we extend the domain of integration to the whole real axis, and the integral evaluated in complete analogy with the two-dimensional case \cite{Akhmedov:2020ryq}. Namely, in the limit \eqref{limit}, the Wightman function is written as
\begin{multline}
    W(x,x') \approx \sum_{lm} \frac{Y_{lm}(\Omega) Y_{lm}^*(\Omega')}{(2M)^2} \int^{\infty}_{-\infty} \frac{d\omega}{4\pi \omega} \frac{1}{e^{\beta_0 \omega}-1} \bigg(e^{-i\omega c} +  e^{2 i \omega \lambda+2i \delta_{\omega l}}
    + e^{i\omega c} + e^{-2i \omega \lambda-2i \delta_{\omega l}} \bigg),
\end{multline}
where $\delta_{\omega l}$ is given by \eqref{delta_sch}.
Evaluating the contribution of the double pole at the $\omega=0$ as in the massless case,
we obtain that
\begin{equation} \label{sing_schw_massive}
    W(x,x') \approx \frac{| \lambda |}{ \beta_0} \frac{\delta(\Omega,\Omega')}{(2M)^2}, \qquad \lambda \to -\infty,
\end{equation}
where $\delta(\Omega,\Omega')$ is the delta function on the two-sphere.

\section{Reissner-Nordström black hole} \label{RNsection}
In this section first we construct modes and Wightman functions outside of the outer horizon. We investigate behavior of the Wightman function on the outer horizon and evaluate the anomalous singularity when both points are located there. The resulting singularity is analogous to the one in the Schwarzschild black hole.

\subsection{Modes and Wightman function}
\begin{figure}[!h]
\centering
\includegraphics[width=0.6\textwidth]{picRN}
\caption{Penrose diagram of the Reissner-Nordström black hole.} \label{picRN}
\end{figure}
Geometry of the Reissner-Nordström black hole is given by the metric
\begin{equation}
    ds^2 = \bigg(1-\frac{2M}{r}+\frac{Q^2}{r^2}\bigg) dt^2 - \frac{dr^2}{\big(1-\frac{2M}{r}+\frac{Q^2}{r^2}\big)}-r^2(d \theta^2 + \sin^2 \theta d \varphi^2).
\end{equation}
This geometry has two horizons at
\begin{equation}
    r_{\pm} = M \pm \sqrt{M^2-Q^2},
\end{equation}
where the surface $r=r_+$ is the event horizon and the surface $r=r_-$ is a Cauchy horizon (we consider non-extremal case $Q<M$). Surface gravity corresponding to these horizons is equal to:
\begin{align}
    \kappa_{\pm} = \frac{r_+-r_-}{2r^2_{\pm}},
\end{align}
and the tortoise coordinate may be defined as
\begin{align}
    r_* = r+\frac{1}{2\kappa_+} \log \bigg( \frac{|r-r_+|}{r_+-r_-} \bigg) -\frac{1}{2\kappa_-} \log \bigg( \frac{|r-r_-|}{r_+-r_-} \bigg).
\end{align}
Also, it will be useful to define the Eddington-Finkelstein coordinates, which are written in the exterior region as
\begin{align}
    u = t-r_*, \qquad v = t+r_*,
\end{align}
In the following I consider the situation in which both points of the Wightman function are located on the outer horizon $r=r_+$, which is denoted as $H_R$ in Fig. \ref{picRN}.

\subsubsection{Massless case}

In the exterior region we can proceed analogously as in the Schwarzschild spacetime. Namely, we split variables as in \eqref{split}

and represent these functions in the form \eqref{f_split}, where the radial functions $R(r)$ satisfy the equation as follows
\begin{align} \label{exterior_radial}
    &\qquad \qquad \qquad  \frac{d^2R}{dr_*^2}+\big[ \omega^2 - V_l (r) \big] R = 0, \nonumber \\
    &V_l(r) = \bigg(1-\frac{2M}{r}+\frac{Q^2}{r^2}\bigg) \bigg( \frac{l(l+1)}{r^2} + \frac{2M}{r^3} - \frac{2Q^2}{r^4} \bigg).
\end{align}
Again we distinguish two types of solutions \eqref{set1} of this equation with boundary conditions \eqref{boundary_cond}, corresponding to out-going and in-going modes. Then, as in the Schwarzschild spacetime, the Wightman function in the generic thermal state is written as
\begin{multline} \label{outer_hh}
    W(x,x') = \sum_{lm} \int^{+\infty}_{-\infty} \frac{d \omega}{4\pi \omega} \frac{1}{rr'}\bigg[ \frac{e^{i \omega (t-t')} \overrightarrow{R}^*_l (\omega|r) \overrightarrow{R}_l (\omega|r') }{e^{\beta_R  \omega}-1}
  +\\+\frac{e^{i \omega (t-t')} \overleftarrow{R}^*_l (\omega|r) \overleftarrow{R}_l (\omega|r') }{e^{\beta_L \omega}-1} \bigg] Y_{lm} (\Omega) Y^*_{lm} (\Omega').
\end{multline}

\subsubsection{Massive case}
In the massive case the equation for the radial function is
\begin{align} \label{radial_rn}
    &\frac{d^2R}{dr_*^2}+\big[ \omega^2 - V_l (r) \big] R = 0,
\end{align}
with the potential
\begin{equation} \label{potential}
    V_l (r) = \bigg( 1 - \frac{2M}{r} + \frac{Q^2}{r^2} \bigg) \bigg(  \frac{l(l+1)}{r^2} + \frac{2M}{r^3} - \frac{2Q^2}{r^4} + \mu^2 \bigg).
\end{equation}
Again, the modes with $\omega^2 \leq \mu^2$ are
\begin{align}
    &\varphi_{\omega lm} (x) = \frac{1}{\sqrt{ \pi |\omega|}} e^{-i |\omega| t} R_{ l} (\omega | r) Y_{lm} (\Omega), \nonumber \\
    &R_{ l} (\omega | r) \approx \frac{1}{r} \cos (\omega r_* + \delta_{\omega l}),
\end{align}
where the expression for the phase acquires the form
\begin{align}
    &\delta_{\omega l} \approx \frac{\pi}{2} +\frac{\omega}{2 \kappa_+} \log \big( \Tilde{l}^2+\mu^2 \big) -\omega r_+ - \text{arg} \, \Gamma(1+\frac{i\omega}{\kappa_+}), \nonumber \\
    &\qquad \qquad \qquad \qquad  \Tilde{l}^2 = l^2+l+1-\frac{r_-}{r_+}.
\end{align}
The modes with $\omega^2 > \mu^2$ are defined as in \eqref{f_modes}. Then, the Wightman function is as follows:
\begin{multline}
    W(x,x') = \sum_{lm} \int^{+\mu}_{-\mu} \frac{d \omega}{\pi \omega}  \frac{e^{i \omega (t-t')} R^*_l(\omega|r) R_l (\omega|r') }{e^{\beta_0 \omega}-1} Y_{lm}(\Omega) Y^*_{lm}(\Omega') + \\
    +\sum_{lm} \int_{|\omega|>\mu} \frac{d \omega}{4\pi \omega} \frac{1}{rr'} \bigg[ \frac{ e^{i \omega (t-t')} \overrightarrow{F}^*_{ l}(\omega | r) \overrightarrow{F}_{l}(\omega|r')}{e^{\beta_R \omega}-1}
 +
  \frac{ e^{i \omega (t-t')} \overleftarrow{F}^*_{l}(\omega|r) \overleftarrow{F}_{l}(\omega|r')}{e^{\beta_L \omega}-1} \bigg] Y_{lm}(\Omega) Y^*_{lm}(\Omega') .
\end{multline}

\subsection{Anomalous singularity}

\subsubsection{Massless case}
As in the previous section, we take the near outer horizon limit $r \to r_+$ \eqref{limit} of the expression \eqref{outer_hh}, evaluate the reflection coefficients $\overrightarrow{A}_l (\omega)$ by solving the radial equation near the outer horizon to obtain
\begin{align}
    &\overrightarrow{R}_l (\omega |r) \approx 2 \Tilde{l}^{-\frac{i\omega}{\kappa_+}}\frac{e^{i \omega r_+}}{\Gamma(-i \omega / \kappa_+)} K_{\frac{i\omega}{\kappa_+}} \big( 2 \Tilde{l} e^{\kappa_+ (r_*-r_+)} \big), \nonumber \\
    &\overrightarrow{A}_l (\omega) \approx - \frac{\kappa_+}{\pi \omega} e^{- \frac{2i\omega}{\kappa_+} \log \Tilde{l} +2i \omega r_+} \sinh \bigg( \frac{\pi \omega}{\kappa_+} \bigg) \Gamma \bigg( 1+\frac{i\omega}{\kappa_+} \bigg).
\end{align}
Substituting this reflection coefficient into the expression for the Wightman function in the outer horizon limit
\begin{multline}
    W(x,x') \approx \frac{1}{(2M)^2} \sum_{lm} Y^*_{lm} (\theta, \varphi) Y_{lm} (\theta',\varphi')  \int_{-\infty}^{\infty} \frac{d \omega }{4 \pi \omega} \bigg[ \frac{1}{e^{\beta_R \omega }-1} \big( e^{-i\omega c} + \overrightarrow{A}_l (\omega) e^{-2 i \lambda \omega } + \\ + \overrightarrow{A}_l^* (\omega) e^{2 i \lambda \omega } + | \overrightarrow{A}_l (\omega) |^2 e^{i\omega c} \big) + \frac{1}{e^{\beta_L \omega}-1} | B_l (\omega) |^2 e^{i \omega c} \bigg],
\end{multline}
and evaluating the integrals as in the Schwarzschild spacetime, one obtains
\begin{equation}
    W(x,x') \approx \frac{|\lambda| }{\beta_R} \, \frac{ \delta(\Omega,\Omega')}{r_+^2}   \, , \qquad \lambda \to -\infty.
\end{equation}

\subsubsection{Massive case}
As is the Schwarzschild spacetime, the leading contribution in the horizon limit \eqref{limit} comes from the low frequencies. Then, for the Wightman function in the horizon limit one has:
\begin{multline}
    W(x,x') \approx \sum_{lm} \int_{-\mu}^{\mu} d \omega \, \frac{\varphi_{\omega lm}^* (x) \, \varphi_{\omega lm} (x')}{e^{\beta_0 \omega}-1} \approx \\
    \approx \sum_{lm} \frac{Y_{lm}(\Omega) Y_{lm}^*(\Omega')}{(2M)^2} \int^{\infty}_{-\infty} \frac{d\omega}{4\pi \omega} \frac{e^{i \omega (t-t')}}{e^{\beta_0 \omega}-1} \bigg( e^{i \omega (r_*+r_*')+2i \delta_{\omega l}} + e^{-i \omega (r_*+r_*')-2i \delta_{\omega l}} \bigg).
\end{multline}
Evaluating the integrals as in the Schwarzschild spacetime, one obtains the result
\begin{equation}
    W(x,x') \approx \frac{| \lambda |}{ \beta_0} \frac{\delta(\Omega,\Omega')}{r_+^2}, \qquad \lambda \to -\infty.
\end{equation}
Thus, again we encounter the anomalous singularity at the outer horizon.

\section{Conclusions and discussions}

It is shown that the anomalous singularity also occurs in the four--dimensional Schwarzschild and Reissner-Nordström black hole backgrounds.
Furthermore, as in Rindler and in de Sitter $d$--dimensional examples \cite{Akhmedov:2020ryq}, it is shown that anomalous singularity does not occur when the coordinates additional to $(t,r)$ do not coincide, and the divergence is amplified if they do. As we can represent delta function of zero argument as a linear divergence, we obtain that the power of the divergence is quadratic, as is expected from the analyisis in the Rindler and de Sitter spacetimes where in $d$ dimensions the divergence was of power $(d-2)$. These anomalous singularities lead to the explosive behavior of the expectation value of the stress-energy tensor, which means that with non-canonical temperature the backreaction of the QFT on the background geometry is strong. Does that mean that one cannot place a black hole in a bath with non-canonical temperature? This is the question for the further study.

One more possible outcome of this anomalous singularity is that it can affect the beta function as the coefficient before the light-like divergence is different than the canonical ones.
It is interesting to see whether this anomalous singularity affects loop contributions to the two-point function in an interacting theory. Also, it is tempting to see how the results of this paper affect computation of the entropy of the fields outside the black hole. However, these problems are out of the scope of this paper and will be considered elsewhere.

\section{Acknowledgements}

I would like to thank E.T. Akhmedov, K.V. Bazarov, and D.V. Diakonov for valuable discussions.

This work was supported by the grant from the Foundation for the Advancement of Theoretical Physics and Mathematics “BASIS” and by Russian Ministry of education and science.


\section{Introduction}

When it comes to MPC strategy for continuous-time systems, we periodically compute an input signal from solving an optimization problem and apply this signal to the system over a bounded horizon, see e.g. \cite{Rawlings2009model,xi2019predictive}.
With the improving of control theory and technology, the new problems which involve large transients often require Nonlinear MPC (NMPC), see e.g. \cite{allgower2012nonlinear,grune2017nonlinear}.
NMPC for continuous-time systems with input constraints is proposed, see e.g. \cite{chen1998quasi}, which considers a nonlinear system given by
\begin{align}\label{eq:system:chen1997}
\dot{x}(t) & = f(x(t), u(t)), \quad x(0) = \bar{x}, \quad t \in \mathbb{R}_{\ge 0}, \\
& \quad  x(t) \in \mathbb{R}^n, \quad u (t) \in \mathbb{U} \subset \mathbb{R}^m,
\end{align}
where ${x}$ and $u$ are the state vector and the control input, respectively; $\mathbb{U}$ is the input constraint set and the mapping $f(\cdot, \cdot )$ satisfy $f(0, 0) = 0$. Besides, $f: \mathbb{R}^n \times \mathbb{U} \to \mathbb{R}^n$ is assumed to be twice continuously differentiable, which leads to unique solutions under the piece-wise right-continuous control inputs, and $\mathbb{U}$ is assumed to be compact, convex, and contains $0$ in its interior.

The quasi-infinite horizon NMPC developed in \cite{chen1998quasi} consists of designing a control invariant region $\Omega$, which is called as \emph{terminal region}, and a control Lyapunov function.
The set $\Omega$ must be chosen such that it is invariant for \eqref{eq:system:chen1997} controlled by a local auxiliary controller. It can be determined off-line such that invariance property of $\Omega$ holds and all the constraints are satisfied. The actual implementation of the NMPC developed in \cite{chen1998quasi} would not use the local auxiliary controller, but rather guarantee recursive feasibility of the closed-loop system. As a contrast, the dual-mode NMPC developed in \cite{Michalska:1993} is actually a receding horizon implementation of a switching control mechanism; that is, it applies the optimal control input from solving a optimization problem over a finite horizon so as to drive the state to the terminal region $\Omega$; once the state enters $\Omega$, it switches the control input to the local feedback law so as to steer the state to the origin.
We highlight that the terminal region $\Omega$ is a positively invariant set for \eqref{eq:system:chen1997} (see e.g.~\cite{bertsekas71}) and the analytic properties of positively invariant sets have been widely studied, see e.g.~\cite{blanchini:2007}. The design method of the set $\Omega$ developed in \cite{chen1998quasi} is based on the results presented in \cite{Michalska:1993}, which is calculated off-line by linearizing the system \eqref{eq:system:chen1997}. This method is quite conservative, especially for systems with high nonlinearities. The terminal region $\Omega$ is very important for stability analysis of the closed-loop NMPC and affects the feasible region of NMPC.
Therefore, this paper considers the data-driven approximations for feasible regions and the terminal regions of the nonlinear systems.

In this paper, we develop a data-driven learning framework for approximating feasible region and invariant set of a perturbed nonlinear systems using feasibility information of data samples.
This work can be seen as an extension of the results developed in \cite{ong2006enlarging} and \cite{chakrabarty2016support,chakrabarty2020active}.
Here, we consider the continuous-time nonlinear system with bounded disturbance and polyhedral state and input constraints.
The main contributions of this paper are as follows: (i) we develop a computation method to obtain the invariant set for the nonlinear system given by \eqref{eq:system:chen1997}; (ii) we develop a data-driven learning method to obtain an inner and outer approximation of the feasible region of the NMPC using support vector machine; (iii) we extend the above two results to perturbed nonlinear system, where disturbance is bounded and the state and input space are polyhedral.

\emph{Notation}:
We denote $\mathbb{R}^n$ for the $n$-dimensional Euclidean space and $\mathcal{C}^1(\mathbb{R}^n)$ for the space of continuously differentiable functions on $\mathbb{R}^n$. Besides, we let $\mathbb{ R}$ and $\mathbb{ R}_{\ge 0}$ denote the real and the non-negative real numbers. For a set $\mathbb{S}$, we let $\mathrm{Card}(\mathbb{S})$ denote the cardinality of $\mathbb{S}$, $\mathrm{Vol}(\mathbb{S})$ denote the Lebesgue measure of $\mathbb{S}$, and $\mathrm{Int}(\mathbb{S})$ denote the set of interior points of $\mathbb{S}$. Furthermore, a continuous function $\alpha: \mathbb{R}_+ \rightarrow \mathbb{R}_+$ belongs to class $\mathcal{K}$ if it is continuous, $\alpha(0) = 0$, and it is strictly increasing. It belongs to class $\mathcal{K}_{\infty}$ if it belongs to class $\mathcal{K}$ and it is unbounded. Finally, given a point $x \in\mathbb{R}^{n}$ and a bounded and non-empty set $\mathcal{A}\subseteq\mathbb{R}^{n}$, the distance of  $x$ to the set $\mathcal{A}$ is defined as $|x|_{\mathcal{A}} = \mathrm{dist}(x, \mathcal{A})   \triangleq    \mathrm{inf}  \{ \|x-\xi\|_{2}: \xi \in\mathcal{A}\}$. A closed ball centralized at a given point $\xi_{0}\in\mathbb{R}^{n}$ with a radius of $\rho \ge 0$ is denoted as $ \mathbb{B}^n(\xi_{0}, \rho) \triangleq \{ \xi \in \mathbb{R}^n : \| \xi - \xi_0\| \le \rho \}$. If $\xi_0 = 0$, we let $\mathbb{B}^n(\rho)$ denote $ \mathbb{B}^n(0, \rho)$ as a shorthand.

\section{Problem formulation}
We extend \eqref{eq:system:chen1997} to a nonlinear perturbed system described as
\begin{align}
\dot{x} {(t)} & = {F} \left(x (t), u (t), w (t) \right), \quad x(0) = \bar{x},  \quad t \in \mathbb{R}_{\ge 0}, \label{eq:system:description} \\
& \quad  x(t) \in \mathbb{X} \subset \mathbb{R}^n, \quad u (t) \in \mathbb{U} \subset \mathbb{R}^m, \label{eq:system:description:2}
\end{align}
where ${x} \in \mathbb{R}^{n}$, $u\in \mathbb{R}^{m}$, and $w \in \mathbb{R}^{n}$ are the state vector, the control input, and the additive disturbance. Besides, we confine the state space $\mathbb{X}$ and the input space $\mathbb{U}$ to be polytopes. We give the following assumptions for \eqref{eq:system:description}:

\begin{assum} \label{eq:J:l}
\begin{description}
	\item[A1] The set $\mathbb{U} \in \mathbb{R}^{m}$ is compact, convex and contains $0 \in \mathbb{R}^{m} $ in its interior. The set $\mathbb{X}$ is a $n$-dimensional hypercube and contains $0 \in \mathbb{R}^{n}$ in its interior.
	\item[A2] The disturbance $w(\cdot) \in \mathbb{W}$, where $\mathbb{W} = \left\{w: |w| \le \bar w \right\}$.
	\item[A3] $F: \mathbb{X} \times \mathbb{U} \times \mathbb{W} \to \mathbb{X}$ is locally Lipschitz in $x$ such that $\|F(x_1, u, w)-F(x_2, u, 0) \|  \le L_F  \|x_1 - x_2 \| + \beta(|w|_{\infty})$	for all ${x_1, x_2} \in \mathbb{X}$, ${u} \in \mathbb{U}$, and $w \in \mathbb{W}$, where $L_F$ is local Lipschitz constant and $\beta(\cdot)$ is a $\mathcal{K}$ function.
\end{description}
\end{assum}

The nominal model of the system \eqref{eq:system:description} is given by
\begin{align}\label{eq:system:nominal}
\dot{\hat{x}} (t) & = {F} \left(\hat x (t), u(t), 0 \right) \triangleq f(\hat x (t), u (t)), \notag \\
& \quad \qquad \qquad \hat x (0) = \bar{x},  \quad t \in \mathbb{R}_{\ge 0},
\end{align}
and the additive disturbance $w(t), t\ge 0$, in \eqref{eq:system:description} satisfies
\begin{align}\label{eq:5}
w(t) & = \dot{x} (t) - f(x (t), u (t)).
\end{align}
Assumption A2 indicates that $w(t)$ is bounded on $\mathbb{W}$. This kind of model uncertainty has been studied in previous papers about robustness in MPC as in \cite{zou2019event} and \cite{lu2019stochastic}.

Given a state measurement $x_k$, a prediction horizon $T$, a state cost $\Psi(\cdot, \cdot)$, a terminal penalty $\Phi(\cdot)$, and a terminal set $\Omega$, an optimization problem of NMPC for \eqref{eq:system:description} is defined as follows:
	\begin{subequations}\label{eq:7}
	\begin{align}\label{FHOCP}
	\hat{u}^\star = \underset{ u } { \arg\min }  ~ \mathcal{J}(x_k,  \hat{u}(\cdot)),
	\end{align}
	with
	\begin{align}\label{Perf}
	\mathcal{J}(x_k,  \hat{u}(\cdot))  \buildrel \vartriangle \over  = \int_{0}^{T} \Psi(\hat x(\tau), \hat u(\tau) ) \mathrm{d} \tau + \Phi(\hat x(T)),
	\end{align}
	subject to
		\begin{align}
		&  \dot{\hat x} (\tau) =   f(\hat{x} (\tau), \hat{u} (\tau)),  \quad 0 \le \tau \le T,   \label{eq:7A}  \\
		&   \hat{x} (\tau) \in \mathbb{X}, \quad \hat{u} (\tau) \in \mathbb{U}, \quad 0 \le \tau \le T,   \label{eq:7B} \\
		&  \hat{x} (0) = x_k,  \quad  \hat x (T) \in \Omega. \label{eq:7C}
		\end{align}
	\end{subequations}
For \eqref{eq:7}, we give the following two assumptions.
\begin{assum} \label{eq:J:II}
	The stage cost $\Psi(x, u)$ in \eqref{Perf} satisfies $\Psi(0, 0) = 0$ and $\Psi(x, u) \ge \alpha_\Psi (\|x\|) + \rho_\Psi$, where $\rho_\Psi\in \mathbb{R}$ and $\alpha_\Psi(\cdot)$ is a $\mathcal{K}_{\infty}$ function. Moreover, $\Psi(x, u)$ is locally Lipschitz in $x$ in $\mathbb{X} \times \mathbb{U}$ with constant $L_\Psi$; that is, $\| \Psi(x_1, u) - \Psi(x_2, u) \| \le L_\Psi \|x_1 - x_2\|$ for all $x_1, x_2 \in \mathbb{X}$ and all $u \in \mathbb{U}$.
\end{assum}
\begin{assum} \label{eq:J:F}
	The terminal cost $\Phi(x)$ in \eqref{Perf} satisfies
	\begin{enumerate}
		\item $\underline{\alpha}_{\Phi} (\|x\|) \le \Phi(x) \le \overline{\alpha}_{\Phi} (\|x\|)$, where $\underline{\alpha}_{\Phi} (\cdot)$ and $\overline{\alpha}_{\Phi} (\cdot)$ are $\mathcal{K}_{\infty}$ functions;
		\item $\Phi(f(x,\kappa(x))) - \Phi(x) \le - \Psi(x, \kappa(x))$ for all $x \in \Omega$, where $\kappa(\cdot): \mathbb{ R}^n \to \mathbb{ R}^m$ is continuous and $\kappa(0) = 0$;
		\item $\Phi(\cdot)$ is locally Lipschitz in $\Omega$ with a Lipschitz constant $L_\Phi$.
	\end{enumerate}
\end{assum}
Assumptions~\ref{eq:J:II} and \ref{eq:J:F} are fairly standard assumptions to guarantee feasibility and stability for nonlinear systems under MPC framework, see e.g. \cite{allgower2012nonlinear,grune2017nonlinear}. Specifically, several numerical methods to obtain $\Omega$ and the local auxiliary controller $\kappa(x)$ satisfying Assumption \ref{eq:J:F} were proposed in \cite{Michalska:1993,chen1998quasi}. The feature of $\kappa(x)$ makes it desirable to construct the feasible control input and further to prove the closed-loop stability. Furthermore, several ways to design $\Psi(\cdot, \cdot)$ and $\Phi(\cdot)$ satisfying Assumptions~\ref{eq:J:l} and \ref{eq:J:II} were proposed in \cite{Hashimoto2016self}.  In practice, $\Psi(x, u)$ and $\Phi(x)$ can be taken as quadratic functions.

Now, we assume that for \eqref{eq:7}, Assumptions~\ref{eq:J:l} and \ref{eq:J:II} hold. Then we introduce the following definitions.
\begin{defn}[Feasible State]  \label{Feasible:state}
	A sampling state $x_k \in \mathbb{X}$ is feasible if there exists a admissible control function $\kappa(\cdot) \in \mathbb{U}$, $t \ge 0$, satisfying all the conditions  \eqref{eq:7A}--\eqref{eq:7C} except the minimization of \eqref{Perf}.
\end{defn}
Using Definition~\eqref{Feasible:state}, we define the feasible region $\mathbb{F}$:
\begin{equation} \label{definition:F}
\mathbb{F} \buildrel \vartriangle \over  = \big\{ x_k \in \mathbb{X} :  x_k ~\mathrm{is}~\mathrm{feasible}  ~\mathrm{for}~ (7) \big\}.
\end{equation}
Moreover, we give another two definitions for the $\mathbb{F}$ in \eqref{definition:F}.
\begin{defn}[Full Region]  \label{Full:region}
	The feasible region $\mathbb{F}$ is full if $\mathbb{F} = \mathbb{X} \subset \mathbb{R}^n$.
\end{defn}
\begin{defn}[Globally Full Region]  \label{Globally: Full:region}
	The feasible region $\mathbb{F}$ is globally full if $\mathbb{F} = \mathbb{X} = \mathbb{R}^n$.
\end{defn}
Ideally we would want that the feasible region $\mathbb{F}$ is a globally full region. But it is impossible for both theoretical and application aspects. The best case is that the feasible region $\mathbb{F}$ is a full region such that, we can execute the receding horizon NMPC all over the state space $\mathbb{X}$ with guarantee of the closed-loop stability. Therefore, the main objective of this paper is to learn the feasible region of a nonlinear system under the MPC framework, as well as the nonlinear perturbed case.
Our method provides freedom in the design of the feasible region needed for closed-loop stability and does not need accurate information of the additively perturbed nonlinear system.
Furthermore, we give some considerations for enlarging the feasible region so as to make it to be a full region.

\section{Basic Results}
In this section, we use the support vector machine learning method to estimate the feasible region for the nonlinear system given by \eqref{eq:system:nominal} based on the feasibility information of a low-discrepancy sequence. Specifically, using deterministic sampling data, the SVM learning method is employed to learn the boundary function of the feasible region.

To that end, we let $\partial \mathbb{F}$ denote the boundary of the feasible region $\mathbb{F}$ given by \eqref{definition:F}. Then we make the following assumption.

\begin{assum}\label{zero:superlevel:set}
	There exists a continuous function $\phi(\cdot) \in \mathcal{C}^1 (\mathbb{X})$ such that, the feasible region $\mathbb{F}$ can be represented as the zero superlevel-set of $\phi(\cdot)$; that is
	\begin{align} \label{zero:superlevel:set1}
		\mathrm{Int}(\mathbb{F}) = \bigl\{ x \in \mathbb{X}: \phi(x) >0 \bigr\}.
	\end{align}
\end{assum}
Now using Assumption~\ref{zero:superlevel:set}, to characterize the feasible region $\mathbb{F}$, we just need to construct the function $\phi(\cdot)$.

\subsection{Point-wise sampling data and SVM learning}
To construct the function $\phi(\cdot)$, we need to sample the state space $\mathbb{X}$ and obtain the point-wise sampling data. To that end, we borrow the standard notion of low-discrepancy sequences from \cite{Niederreiter1988low} and construct a data set on a multilevel sparse-grid.
\begin{defn} [\cite{chakrabarty2016support}]
	The discrepancy of a sequence $\{x_{k}\}_{k=1}^{N} \subset \mathbb{X}$ is defined as
	\begin{align*}
	\mathscr{D}_{N}\left(  \{x_{k}\}_{k=1}^{N} \right)  \triangleq   \underset{X\in \mathscr{T}} {\mathrm{sup}} \bigg|\frac{ \mathrm{Card}_{\mathbb{X}} (\mathbb{S})}{N} - \frac{\mathrm{Vol}(\mathbb{S})}{\mathrm{Vol}(\mathbb{X})} \bigg|,
	\end{align*}
	where $\mathscr{T}$ is a set of state dimensional intervals and defined as $\mathscr{T}  \buildrel \vartriangle \over  = \prod_{k=1}^{M} \bigl[a_k,b_k\bigr)=\bigl\{ x_k \in \mathbb{X}:a_k\le x_k<b_k \bigr\} \subset\mathbb{X}$, and $M = \mathrm{Card}_{\mathbb{X}} (\mathbb{S}) \triangleq \mathrm{Card} \bigl\{ k\in\left\{ 1,\ldots,N\right\} : x_k \in \mathbb{S} \bigr\}$. For a low-discrepancy sequence $\{x_{k}\}_{k=1}^{N}$, $\underset{N\to\infty}{\mathrm{lim}} \mathscr{D}_{N} \bigl( \left\{ x_k\right\} _{k=1}^{N} \bigr)  =  0$.
\end{defn}

Over the state space $\mathbb{X}$, we generate  a point-wise data set
\begin{align}\label{data:set}
\varPi   \buildrel \vartriangle \over  = \left\{ x_k\right\} _{k=1}^{N}  = \bigl\{ x_1, x_2, \ldots, x_N \bigr\},
\end{align}
and assume that it satisfies the following assumption.

\begin{assum} \label{low:discrepancy:sequence}
The sequence $\left\{ x_k\right\} _{k=1}^{N}$ is a low-discrepancy sequence on $\mathbb{X}$.
\end{assum}

Now, using Assumption~\ref{low:discrepancy:sequence}, we describe the training data in $\varPi$ following the results presented in \cite{chakrabarty2016support}.
For each point $x_k$, $k=1, \ldots, N$ in $\varPi$, we solve the optimization problem \eqref{eq:7} over a finite-time horizon $T$.
If there is a feasible solution for $x_k$, then we label it as $\mathcal{O}(x_k) = +1$; otherwise, we label it as $\mathcal{O}(x_k) = -1$. Then we have the following two cases:
\begin{align}\label{Bi:classifier}
\mathcal{O}(x_k)  = \begin{cases}
+1, \quad x_k \in \mathbb{F}; \\
-1, \quad x_k \notin \mathbb{F}, \\
\end{cases}
\end{align}
for each point $x_k$ in $\varPi$.  Note that this is a two-class pattern recognition problem, which divides the data set $\varPi$ into two classes, that is,
\begin{align*}
\varPi^{+} & =\big\{ x_{k} \in \varPi: ~ \mathcal{O}(x_k)  = +1 \big\};  \\
\varPi^{-} & =\big\{ x_{k} \in \varPi: ~ \mathcal{O}(x_k)  = -1 \big\}.
\end{align*}
This problem can be solved using the SVM bi-classifier. Here, we use the kernel-based SVM classification method, see e.g. \cite{vapnik2013nature,burges1998tutorial}. Now, we can construct a decision function $\phi_N(\cdot) \in \mathcal{C}^1 (\mathbb{X})$,
\begin{align}\label{decison:function}
\phi_N(x): \varPi \to y_k \in \{+1, -1\},
\end{align}
 that accurately classifies an arbitrary state in $\varPi$ as feasible or infeasible.

Next, to facilitate the classification, we introduce a map: $\Gamma: \mathbb{X} \to \mathbb{H}$, which maps the data in $\varPi$ to a higher-dimensional Hilbert space $\mathbb{H}$ where the data is linearly separable, see e.g. \cite{burges1998tutorial}. In this case, we can let $\phi_N(x) = \rho \Gamma(x) + v$, where $\rho$ and $v$ are vector parameters that determine the orientation of the separating hyperplane. Referring to \cite{ong2006enlarging}, our problem can be formulated as the constrained optimization problem
\begin{subequations}\label{non:constrained:OP}
	\begin{align}\label{non:constrained:OP:target}
	\rho^\star = \underset{ \rho } { \arg\min }  ~  \frac{1}{2} \rho^\mathrm{T} \rho + L \sum_{k =1}^{N} s_k,
	\end{align}
	subject to
	\begin{align}
	&  \rho  \Gamma(x_k) + v \ge 1 - s_k, \quad  s_k \ge 0, \quad \forall k \in \varPi^{+},   \label{non:constrained:OP:1}  \\
	&  \rho  \Gamma(x_k) + v \le -1 + s_k, \quad \forall k \in \varPi^{-},   \label{non:constrained:OP:2}
	\end{align}
\end{subequations}
where $L$ is a regularization parameter and $s_k$ is a slack variable to relax separability constraints.

The optimization problem \eqref{non:constrained:OP} is convex and can be solved via its dual formulation. Letting $\alpha_i$ and $\alpha_j$ denote the Lagrange multiplier vector for \eqref{non:constrained:OP:1} and \eqref{non:constrained:OP:2}, \eqref{non:constrained:OP} can reformulated as the dual problem
\begin{subequations}\label{non:constrained:OP2}
	\begin{align}\label{non:constrained:OP2:target}
\alpha^{\star} = \underset{\alpha}{\mathrm{arg~min}}  &  ~\frac{1}{2}\sum_{i=1}^{N}\sum_{k=1}^{N}\alpha_i \alpha_{k} y_{i}y_{k} \Gamma(x_i) \Gamma(x_k) - \sum_{k=1}^{N}\alpha_k, \\
\mathrm{subject\;to}~~
	&  \sum_{k=1}^{N}\alpha_k y_{k}=0, \label{non:constrained:OP2:1}  \\
	&   0 \le \alpha_{k} \le L,  \quad  k = 1, \ldots, N. \label{non:constrained:OP2:2}
	\end{align}
\end{subequations}
After solving \eqref{non:constrained:OP2}, we obtain the SVM decision function as
\begin{equation}\label{SVM:function}
\phi_N(x)  = \sum_{k=1}^{N}\alpha_{k}^{\star} y_{k} \Gamma(x_k) \Gamma(x) + v,
\end{equation}
and the estimated feasibility region boundary is given by
\begin{equation} \label{SVM:functionregion}
\mathcal{S}_N  \triangleq   \left\{ x\in\mathbb{X}:  \phi_N(x) = 0 \right\}.
\end{equation}

To improve accuracy of the SVM learning when effective data in $\varPi^{+}$ are not enough, we iteratively select new points satisfying $\mathcal{O}(x_k)  = +1$ and expand into the data set $\varPi^{+}$ to ensure richness of data.
This ensures that a non-trivial decision function \eqref{SVM:function} exists.
Furthermore, for the inner product $\Gamma(x_k) \Gamma(x)$ in \eqref{SVM:function}, it can be replaced by using the kernel function. A common choice of kernel function is the Gaussian kernel:
	\[
\Gamma(x_k) \Gamma(x) = \frac{{\|x_{k}-x\|}^2}{2\sigma^{2}},
\]
where $\sigma^2$ is the kernel variance, see e.g. \cite{steinwart2001influence}.

\subsection{Inner and outer approximations of $\partial\mathbb{F}$}
The outer and inner approximations of a bounded region were introduced in \cite{darup2012low,deffuant2007approximating}. Using \eqref{zero:superlevel:set1}, we define the
strict outer ($\mathbb{F}^{-}$) and inner ($\mathbb{F}^{+}$) region of the feasible region $\mathbb{F}$ as
\begin{align}
\mathbb{F}^{+} & \triangleq \left\{ x\in\mathbb{X}:\phi(x)>\varepsilon\right\};  \label{F+}\\
\mathbb{F}^{-}  & \triangleq  \left\{ x\in\mathbb{X}:\phi(x)< - \varepsilon\right\}.  \label{F-}
\end{align}
for some $\varepsilon >0$.
In this subsection, we construct strict outer and inner approximations of the feasible region $\mathbb{F}$, which ensures that they contain no feasible (infeasible) samples.

Now, based on the super-level sets of $\phi_N(x)$, we can find a positive constant $\varepsilon^{+}>0$ which defines the inner approximation $\mathbb{F}_N^{+}$ by
\begin{align}
\mathbb{F}_N^{+} & \triangleq \bigl\{ x\in \mathbb{X}: \phi_N(x) > \varepsilon^{+}\bigr\}. \label{eq:outer}
\end{align}
Using the data in $\varPi^+$, we determine the parameter $\varepsilon^{+}$ by solving the optimization problem
\begin{subequations}\label{eq:outer:OP}
	\begin{align}
	\varepsilon^{+}  & =\underset{\varepsilon>0}   {\mathrm{arg~min}} ~ \varepsilon,  \label{eq:outer:OP:1}  \\
	\mathrm{subject\;to}~ & \phi_N (x_{k} ) > \varepsilon,  \quad\forall x_{k}\in \varPi^+. \label{eq:outer:OP:2}
	\end{align}
\end{subequations}

Analogously, there exists a negative constant $\varepsilon^{-} <0$ which defines the outer approximation $\mathbb{F}_N^{-}$ by
\begin{align}
\mathbb{F}_N^{-} & \triangleq \bigl\{ x\in  \mathbb{X}: \phi_N(x) < \varepsilon^{-}\bigr\}.  \label{eq:inner}
\end{align}
This constant $\varepsilon^{-}$ can be determined by using the data in $\varPi^-$ and solving the optimization problem
\begin{subequations}\label{eq:inner:OP}
	\begin{align}
	\varepsilon^{-}  & =\underset{\varepsilon<0}   {\mathrm{arg~max}} ~ \varepsilon,  \label{eq:inner:OP:1}  \\
	\mathrm{subject\;to}~ & \phi_N (x_{k} ) < \varepsilon,   \quad\forall x_{k}  \in \varPi^-.  \label{eq:inner:OP:2}
	\end{align}
\end{subequations}

We highlight that novel algorithms for constructing strict inner and outer approximations of the feasible region $\mathbb{F}$ are developed by solving \eqref{eq:outer:OP} and \eqref{eq:inner:OP}. Besides, the inner approximation $\mathbb{F}_N^{+}$ is based only on the feasible samples in $\varPi^+$ and the outer approximation $\mathbb{F}_N^{-}$ is based only on the infeasible samples in $\varPi^-$.

\subsection{Convergence analysis}
In this subsection, we provide convergence analysis for the approximation error.
We show that as the number of samples $N$ from the low-discrepancy sequence increases, the inner and outer approximations with universal kernels will converge to a strict approximation of the actual feasible region $\mathbb{F}$.

First, we introduce definition of universal kernels adopted from \cite{steinwart2001influence} and give a proposition.
\begin{defn} \label{universal:kernel}
	A continuous kernel $\mathcal{K}(\cdot, \cdot)$ is universal
	if the space of all functions induced by $\mathcal{K}(\cdot, \cdot)$ is dense in $\mathcal{C}(\mathbb{X})$.
\end{defn}
\begin{prop} \label{universal:kernel:property}
	Every universal kernel separates all compact subsets of $\mathbb{X}$.
\end{prop}
Definition~\ref{universal:kernel} indicates that for each continuous function $\phi$ and $\varepsilon>0$,
there exists a function $\phi_{\mathcal{K}}$ induced by a universal kernel $\mathcal{K}$ such that $|\phi-\phi_{\mathcal{K}}|<\varepsilon$.
Proposition~\ref{universal:kernel:property} indicates that the kernel-based SVM classification method can be adopted to solve the feasibility boundary learning problem.
In this case, if we design
\begin{align} \label{kernel:function}
	\mathcal{K}(x_k, x) = \Gamma(x_k) \Gamma(x),
\end{align}
then we can accurately separate the pairwise disjoint compact sets $\varPi^+$ and $\varPi^-$ in $\mathbb{X}$.

We highlight that there is not computational method for determining such a separating function in \cite{steinwart2001influence}.
By modifying the arguments presented in \cite{steinwart2001influence}, we collect the training data from a low-discrepancy sequence.
It guarantees that increasing the number of samples $N$ for
the SVM with a universal kernel results in an increasing accuracy of the feasible region learning, and reduces the conservativeness of the inner and outer approximations of $\mathbb{F}$.
We summarize our discussions as the following theorem, which is slightly similar to the results given in \cite{chakrabarty2016support}.

\begin{thm}\label{theprem:main}
	Consider the system \eqref{eq:system:nominal}, suppose Assumptions~\ref{eq:J:l}, \ref{eq:J:II}, \ref{eq:J:F}, \ref{zero:superlevel:set}, and \ref{low:discrepancy:sequence} hold, and let $\mathcal{K}(\cdot, \cdot)$ be a universal kernel on $\mathbb{X}$ in \eqref{kernel:function} and $\phi_{N, \mathcal{K}} (x)$ be the kernel-based SVM decision function from \eqref{SVM:function}. Then for any $L>0$ in the optimization problem \eqref{non:constrained:OP2} and for two compact subsets $\mathbb{F}^{+} \subset \mathbb{X}$ and $\mathbb{F}^{-} \subset \mathbb{X}$ in \eqref{F+} and \eqref{F-},
	\begin{align}
	& \underset{N \to \infty}{ \mathrm{lim}} \phi_N(x) = \phi(x), \\
	&	\underset{N \to \infty}{ \mathrm{lim}}  \mathbb{F}_N^+ = \mathbb{F}^+, \quad \underset{N \to \infty}{ \mathrm{lim}}  \mathbb{F}_N^- = 	\mathbb{F}^-, \\
	&	\underset{N \to \infty}{ \mathrm{lim}} \mathcal{S}_N  = \partial \mathbb{F}.
	\end{align}
\end{thm}

Theorem~\ref{theprem:main} indicates that there exists a sufficiently large number of samples $N_{0}$ such that if $N\ge N_{0}$,
\[
\begin{cases}
\phi(x) > 0,  &  \quad  \mathrm{for} \; x\in\mathbb{F}^{+};  \\
\phi(x) < 0,  &  \quad  \mathrm{for} \; x\in\mathbb{F}^{-},
\end{cases}
\]
and the strict inner and outer approximations of the feasible region $\mathbb{F}$ can be obtained by sampling $\mathbb{X}$ using low-discrepancy sequences.
Therefore, the SVM classifier using universal kernels is able to approximate our feasible region boundaries with sufficiently high accuracy.

\section{Main Results}
In Section~3, we present a method to estimate the feasible region where the states of the nominal system \eqref{eq:system:chen1997} can be driven to the terminal set $\Omega$ while satisfying the constraints \eqref{eq:system:description:2}. In this section, we extend this result to the perturbed nonlinear system given by \eqref{eq:system:description}. However, the region $\mathbb{F}$ obtaining from Section 3 may not be feasible for the perturbed system when there exists unknown disturbance.

Using Assumption A2, \eqref{eq:system:description}, and \eqref{eq:system:nominal}, it follows that
\begin{align}
| F(x(t), u(t), w(t)) - f(x(t), u(t)) | \le \bar{w}.
\end{align}
So for the perturbed system \eqref{eq:system:description}, the set $\Omega$ may not be positively invariant with the local feedback law $u(t)$. In this case, we introduce the robust control invariant (RCI) set for \eqref{eq:system:description} following the results presented in \cite{sun2018robust}.
\begin{defn}[RCI set]  \label{RCIset}
	Consider the perturbed system \eqref{eq:system:description} with constraints $x \in \mathbb{X}$, $u\in \mathbb{U}$ and $w \in \mathbb{W}$. The set $\bar{\Omega} \subset \mathbb{X} \subseteq \mathbb{R}^n$  is an RCI set, if there exists an ancillary controller $\kappa(x)$ with $\bar{u} = u + \kappa(x) \in \mathbb{U}$ such that $\forall \bar{x} \in \bar{\Omega}$ and $\forall w \in \mathbb{W}$, $x(t) \in \bar{\Omega}$.
\end{defn}

For the design of the RCI set $\bar{\Omega}$ (terminal region) for \eqref{eq:system:description}, the set-theoretic methods  and tube-based MPC scheme are very useful, see e.g. \cite{blanchini:2007} and \cite{Rakovic:2012}. However, the method developed in \cite{Rakovic:2012} can not be applied to our problem, since the perturbed system given by \eqref{eq:system:description} is continuous-time and nonlinear case; the authors in \cite{ong2006enlarging} and  \cite{chakrabarty2016support}  only consider the nonlinear systems, not for the perturbed system.

\subsection{RCI set design based on set-theoretic method}

Recalling that the disturbance in \eqref{eq:5} is actually an additive disturbance and satisfies
\begin{align} \label{eq:w}
w (t) = F(x(t), u(t), w(t)) - f(x(t), u(t)), ~ t \in \mathbb{ R}_{\ge 0}.
\end{align}
Besides, we note that the state and control do not depend on the disturbance $w$.

Now, using Assumption~A2, we introduce the following set
\begin{align} \label{eq:Set}
	\Omega_0 \triangleq \bigl\{ x \in \Omega: |x|_{\partial  \Omega  } \ge \bar{w} \bigr\},
\end{align}
where $\Omega$ is a terminal region for \eqref{eq:system:nominal} and $\partial  \Omega $ is the boundary of the set $\Omega$.
Then for \eqref{eq:Set}, we have the following lemma.

\begin{lem}  \label{Lem:Omega}
	Consider the set \eqref{eq:Set} and assume that $\Omega \supset \mathbb{R}^n (\bar{w})$ holds. Then, $\Omega_0$ is a RCI set for the perturbed system \eqref{eq:system:description}.
\end{lem}
\begin{pf}
	Note that $\Omega$ is a control invariant set for the nonlinear system \eqref{eq:system:nominal} and $\Omega_0 \subset \Omega$. Since $\Omega \supset \mathbb{R}^n (\bar{w})$, $\Omega \setminus \Omega_0 \neq \varnothing$. Then using \eqref{eq:Set}, we have that for any $x \in \Omega_0$, $ |x + w |_{ \partial  \Omega_0 } \ge |x|_{ \partial  \Omega } - \bar{w} \ge 0 $ for all $ w \in \mathbb{W}$. That is, for any $x \in \Omega_0 $, $x + w \in \Omega$ for all $ w \in \mathbb{W}$. According to Definition~\ref{RCIset}, $\Omega_0$ is a RCI set for \eqref{eq:system:description}.
\end{pf}

Lemma~\ref{Lem:Omega} indicates that the RCI set $\bar{\Omega}$ for \eqref{eq:system:description} can be obtained using the terminal set $\Omega$ for \eqref{eq:system:nominal}, while the set $\Omega$ can be computed using the results developed in \cite{chen1998quasi} for perturbed nonlinear  systems.

\subsection{Feasible region estimation for robust NMPC}
In this subsection, we develop a feasible region estimation method for \eqref{eq:system:description} under the NMPC scheme using the feasibility information of the low-discrepancy data samples.

First, using the RCI set $\Omega_0$ in \eqref{eq:Set} for the perturbed system given by \eqref{eq:system:description}, we modify the problem \eqref{eq:7} to a new optimization problem of robust NMPC which is given by
	\begin{subequations}\label{eq:robust:NMPC}
	\begin{align}\label{robust:NMPC}
	\hat{u}^\star = \underset{ u } { \arg\min }  ~   \mathcal{J} (x_k,  \hat{u}(\cdot)),
	\end{align}
	subject to
	\begin{align}
	&  \dot{\hat x} (\tau) =  f(\hat{x} (\tau), \hat{u} (\tau)),  \quad 0 \le \tau \le T,   \label{eq:robust:NMPC:A}  \\
	&   \hat{x} (\tau) \in \mathbb{X}, \quad \hat{u} (\tau) \in \mathbb{U}, \quad 0 \le \tau \le T,   \label{eq:robust:NMPC:B} \\
	&  \hat{x} (0) = x_k,   \quad \hat x (T) \in \Omega_0,    \label{eq:robust:NMPC:C}
	\end{align}
\end{subequations}
where $\mathcal{J}(x_k,  \hat{u}(\cdot))$ is given by \eqref{Perf} and $\Psi(\cdot, \cdot)$ and $\Phi(\cdot)$ satisfy Assumptions~\ref{eq:J:l} and \ref{eq:J:F}.

Next, using the same procedure developed in Section 3, we generate the data set and construct an approximate function of the boundary for the feasible region. Then similar to Theorem~\ref{theprem:main}, we have the following theorem which shows the approximation of the feasible region $\bar{\mathbb{F}}$ for the perturbed systems in \eqref{eq:system:description}.

\begin{thm}\label{theprem:main:two}
		Consider the system \eqref{eq:system:description}, suppose Assumptions~\ref{eq:J:l}, \ref{eq:J:F}, \ref{zero:superlevel:set}, and \ref{low:discrepancy:sequence} hold, and let $\bar{\Omega}$ be the maximal RCI set in $\mathbb{X}$
	and $\bar{\mathbb{F}}$ be the maximal feasible region for \eqref{eq:system:description}. Then,
	$\Omega_0 \subset \bar{\Omega}$ and
	\begin{align} \label{eq:Region}
    \bar{\mathbb{F}}_0 \triangleq \bigl\{ x \in {\mathbb{F}}: |x|_{\partial  {\mathbb{F}}  } \ge \bar{w} \bigr\} \subset \bar{\mathbb{F}}.
	\end{align}
\end{thm}
Theorem~\ref{theprem:main:two} indicates that based on the feasible region $\mathbb{F}$ for \eqref{eq:system:nominal}, we can obtain the conservative terminal region $\Omega_0$ and the conservative feasible region $\bar{\mathbb{F}}_0$ for the perturbed systems. Since the disturbance $w$ in \eqref{eq:w} is unknown, it is very difficult to compute the strict feasible region for the perturbed system \eqref{eq:system:nominal}. Therefore,  we have adopted the low-discrepancy sequences, optimal kernel-based SVM classification, and  set-theoretic methods to increase the accuracy of the data-driven feasible region estimation. Although the feasible region $\bar{\mathbb{F}}$ can be enlarged with a longer prediction horizon, it needs more computation resources for solving the optimization problems \eqref{eq:robust:NMPC}.

\section{illustrative numerical example}\label{sec:examples}
To illustrate the key ideas presented in this paper, we consider a perturbed nonlinear system representing the cart-and-spring system adopted from \cite{magni2003robust,li2014distributed} given by
\begin{equation*}\begin{cases}
\dot{x}_{1} (t) =  x_{2}(t), \\
\dot{x}_{2} (t) =  - \frac{k_0}{M} e^{- x_1(t)} x_1(t)  - \frac{h_d}{M} x_2(t) + \frac{u(t)}{M} + w(t),
\end{cases}
\end{equation*}
where $x = \left[x_1; x_2 \right] \in \mathbb{R}^2$, $u \in \mathbb{R}$,  and $w \in \mathbb{R}$. The parameters are given as: $M = 1.8~\mathrm{kg}$, $k_0 = 1.2~\mathrm{N/m}$, $h_d =0.25~\mathrm{N\cdot s/m}$,
The constrained state and control space associated with the system is given by $\mathbb{X} = \{ x \in \mathbb{R}^2: |x| \le 2 \}$ and $\mathbb{U} = \{ u \in \mathbb{R}: |u| \le 3\}$ and let $w(t) \in \mathbb{B}^1(0.01)$.
	%
	To satisfy Assumptions~\ref{eq:J:l} and \ref{eq:J:F}, we set the stage cost as $\Psi(x, u) = \|x\|_{Q}^2 + \|u\|_R^2$ with $Q =I_2$ and $R = 0.1$ and the terminal cost as $\Phi(x) = \|x\|_{P}^2$ with $P = \left[5.0511, -2.2731; -2.2731, 2.4586\right]$. The local controller is designed as $\kappa(x) = -Kx$ with $K = \left[4.2291~4.8551\right]$ using the procedure presented in \cite{chen1998quasi}.

	Next, the terminal region $\Omega$ is parametrized by $\Omega(\mu) = \bigl\{ x \in \mathbb{R}^2: x^\mathrm{T} P x \le \mu \bigr\}$, where $\mu \in (0, 1)$.
	%
	Then, we use a kernel-based SVM with $\sigma = 0.8$ for the feasible region approximation. Using the feasibility information of data samples, approximations of the feasible  region boundaries with different prediction horizons $T$ and different $\mu$  are illustrated in Figures.~\ref{figure:feasible:region:SVM:T10s}--\ref{figure:feasible:region:SVM:T20s}.
	Furthermore, it can be seen from Figures~\ref{figure:feasible:region:SVM:T10s} and \ref{figure:feasible:region:SVM:T15s} that larger terminal region results in larger feasible region. From Figures~\ref{figure:feasible:region:SVM:T10s} and \ref{figure:feasible:region:SVM:T20s}, we can see that longer prediction horizon results in larger feasible region.

\begin{figure}[ht]
	\includegraphics[width=.99\linewidth]{feasible_region_SVM_T1005.eps}
	\caption{The boundary of the feasible region~$\bar{\mathbb{F}}_0$ resulting from the optimization  \mbox{(\ref{eq:robust:NMPC})} with the terminal region $\Omega_0(0.5)$ and the prediction horizon $T = 1.0~s$ in \textcolor{red}{red}.
				%
		The infeasible samples in $\varPi^{-}$ are shown in \textcolor{black}{black};
		the feasible samples in $\varPi^{+}$ are shown in \textcolor{blue}{blue};
		the feasible samples selected using \eqref{SVM:functionregion} are shown in \textcolor{green}{green};
	   the boundary of the terminal region $\Omega_0(0.5)$ is shown in \textcolor{magenta}{magenta}.
   }
	\label{figure:feasible:region:SVM:T10s}
	\vspace{4mm}\end{figure}
\begin{figure}[ht]
	\includegraphics[width=.99\linewidth]{feasible_region_SVM_T1009.eps}
	\caption{The boundary of the feasible region~$\bar{\mathbb{F}}_0$ resulting from the optimization  \mbox{(\ref{eq:robust:NMPC})} with the terminal region $\Omega_0(0.9)$ and the prediction horizon $T = 1.0~s$ in \textcolor{red}{red}.
	The other notations are the same as Fig.~\ref{figure:feasible:region:SVM:T10s}.
}
	\label{figure:feasible:region:SVM:T15s}
\end{figure}
\begin{figure}
\includegraphics[width=0.99\linewidth]{feasible_region_SVM_T2005.eps}
\caption{The boundary of the feasible region~$\bar{\mathbb{F}}_0$ resulting from the optimization  \mbox{(\ref{eq:robust:NMPC})} with the terminal region $\Omega_0(0.5)$ and the prediction horizon $T = 2.0~s$ in \textcolor{red}{red}.
The other notations are the same as Fig.~\ref{figure:feasible:region:SVM:T10s}.
}
\label{figure:feasible:region:SVM:T20s}
\end{figure}

\section{Conclusions and future work}\label{sec:conclusion}
This paper develops a data-driven learning framework for designing the terminal region and approximating the
feasible region offline using the feasibility information of low-discrepancy data samples and support vector machine
learning. Our approach provides the freedom in the design of the feasible region needed for closed-loop stability, when
disturbance is subject to deterministic bounds. Future work will combine the offline feasible region learning with
online NMPC design.


\section{Conclusion}\label{sec:conclusion}

In this paper we have demonstrated an asynchronous, self-timed
inference datapath design with area and power of equal orders of magnitude to the
synchronous equivalent. Early propagation enables $10\times$ lower inference latency
than the equivalent synchronous circuit on average. The savings
are enabled by a reduced \ac{CD} scheme which can be applied to any dual-rail{}
asynchronous circuit. The new scheme introduces a timing assumption which can be
incorporated into the \ac{CD} circuit, so that the circuit environment does not need to be adapted.
This type of low-latency circuit can have applications in speech
recognition for wearables and other low-power applications where inference
latency is of particular importance.

In future work we will apply asynchronous design styles to the training
datapath of the \ac{TM} algorithm in order to enable a fully-asynchronous
\ac{ML} hardware capable of on-chip learning.

\section{Inference Datapath}\label{sec:infer}

The inference datapath of the \ac{TM} is derived
from the full \ac{TM} diagram (\cref{fig:tsetlin}). The \acp{TA} and their feedback
are not required for inference.
Only the \emph{exclude} action output is required from the \acp{TA} teams. In the
diagram this is abstracted to the \ac{PI} \signal{e}.
We split the majority voting of the \ac{TM} into two sections. Firstly we
distinctly count all positive votes and negative votes by means of population count{}s.
Secondly the two counts are compared using a magnitude comparator to determine
the winner. The result of the comparison is taken as the classifier outcome.

There are several ways to construct the circuit architecture. This architecture
has been chosen due to the simplicity and efficiency of the asynchronous
magnitude comparator as will become clear in \cref{sec:comparator}.

All \acp{PI} and \acp{PO} of the circuit are dual-rail{} encoded. These can interface
natively with other dual-rail{} signals, or with synchronous circuits using
converters~\cite{Sokolov2006}.

\subsection{Clause Calculation}

The \signal{e} input to the inference datapath controls whether the
corresponding feature input{} (\signal{f}) will be excluded from a clause computation.
We use \gates{or} to form a mask of each feature input{} in each clause.
The partial clause
values, \signal{pc}, must be aggregated using an \gate*{and} tree in order
to evaluate the entire clause comprising input from all \signal{f} and their
associated automaton actions. The \emph{exclude} signals (\signal{e}) from the
\acp{TA} mask \signal[m]{f} and \nsignal[m]{f} \finput{}s{} causing logic-1 at the
\gate{and} inputs. If \signal[2m]{e} (\resp{} \signal[2m+1]{e}) is logic-0
(\ie{} the feature input{} should be \emph{included} in the clause calculation), the
value of \signal[m]{f} (\resp{} \nsignal[m]{f}) is passed through to the
\gate{and} to be evaluated.
The partial clause evaluation circuit is replicated as many times as there are
\finput{}s{} to the \ac{TM}.

Since \signal[m]{f} will be dual-rail{} encoded in our system, we do not need to
generate \nsignal[m]{f} internally. By performing direct mapping of a single-rail{}
circuit, and along with negative gate optimization~\cite{Sokolov2006}, we arrive
at the optimized dual-rail{} circuit in \cref{fig:clause-dr-opt2}. All signal paths in
this circuit have a single inversion---satisfying spacer requirements and giving
the block an \emph{inverting spacer} overall.

\subsection{Population Count{}}

We base our population count{} circuit on the optimized design of
Dalalah~\cite{Dalalah2006}. The eight-input design comprises nine half-adders,
two full-adders and two \gates{or} and is illustrated in \cref{fig:popcount-sr}.
Each wire in the diagram represents two signals which form the dual-rail{} encoding.
The dual-rail{} \gate{or} is internally constructed from one \gate{or} and one
\gate{and}. The dual-rail{} half-adders are constructed using two complex gates and
two simple gates each. There is no spacer inversion within the half-adders as
all signal paths have an even number of inversions. The dual-rail{} full-adder is
constructed from six complex gates, two simple gates and four
inverters~\cite{Wheeldon2019}. It has inverted spacers on carry-in and carry-out
with respect to the other inputs and outputs, therefore we must accommodate for
these in the population count{} design by adding spacer inverters:
\begin{inparaenum}
	\item between \bname[8]{ha} and \bname[0]{fa};
	\item between \bname[1]{fa} and the \signal[3]{y} output.
\end{inparaenum}
The resulting dual-rail{} population count{} circuit has no spacer inversion overall,
therefore the output spacer will have the same polarity as the input spacer.

\begin{figure}[h]
	\tikzinput{infer-diagram}
	\caption{%
		Block diagram of the \acf{TM} inference path. Including dual-rail{}
		circuits for partial clause evaluation; and population count{}
		architecture, where each wire represents two signals with dual-rail{}
		encoding.
		spinv: spacer inverter.
	}\label{fig:dr-circuits}\label{fig:clause-dr-opt2}\label{fig:popcount-sr}
\end{figure}

\subsection{Magnitude Comparator}\label{sec:comparator}

The magnitude comparator compares the number of votes from the positive and
negative \ac{TA} teams. A larger number of positive votes indicates that the
input pattern belongs to the class in question, and conversely, a larger number
of negative votes indicates that the input pattern \emph{does not} belong to the
class in question.

The magnitude comparator is based on a request architecture~\cite{Wheeldon2019}
and compares the operands in bit-pairs, starting from the most significant bit.
Once a difference is found, the answer is known, and the remainder of the bits
need not be compared. This architecture enables huge average-case latency
improvement over a synchronous counterpart. Energy savings are also made by
due to saved switching power on the lower bits when the operands differ by a large
magnitude.

Since the comparator's outputs (less, equal, and greater) are mutually exclusive,
we take advantage of this in the asynchronous design. We use a
1-of-3 encoding on the output instead of the usual dual-rail{}---1-of-n encoding being
a superset of dual-rail{}. Provided a spacer{} seperates the valid{}s, the
switching of 1-of-n codes is monotonic~\cite{Bainbridge2003},
therefore satisfying \cref{ass:monotonic-circuit}.
Without this trick, three sets of dual-rail{} signals would be required at the comparator outputa at the expense of more logic to drive these signals.
The inputs to the comparator are dual-rail{} encoded.

\subsection{Inference Datapath Results}

\begin{table*}[!htb]
	\centering
	\caption{%
		Comparison of single-rail{} and dual-rail{} circuits after synthesis.
	}\label{tab:results}
	\input{results-table}
\end{table*}

The inference datapath was synthesized using \pnoun{Synopsys Design Compiler}
for two different \SI{65}{\nano\metre} silicon libraries. \processlibname{umc65ll} is a
commercially available, low-leakage library which we use with nominal
\SI{1.2}{\volt} supply and TT corner. \processlibname{pipistrelle4} is a custom
library aimed at high performance subthreshold
operation~\cite{Morris2017}. It uses a full diffusion sizing strategy with
non-minimum-length transistors in order to mitigate subthreshold effects. For
this silicon library the circuit is first synthesized at TT corner for nominal
\SI{1.2}{\volt} supply and results are shown with supply voltage in the range
\SIrange{0.25}{1.2}{\volt}.

Results in \cref{tab:results} show similar cell areas for both single-rail{} and dual-rail{}
designs for each silicon library. This is possible due to the careful choice
of dual-rail{} circuit architecture and the reduced completion detection scheme.
The dual-rail{} clause computation and magnitude comparator are more area efficient
than their single-rail{} conterparts due to exploitation of dual-rail{} encoding and clever
use of 1-of-3 encoding respectively.

For
the area of the sequential cells we count flip-flop area for the single-rail{} designs
and C-element{} area for dual-rail{} designs.
The sequential area is similar between designs, despite the dual-rail{} design
having twice as many sequential cells due to the doubled input rails. The
dual-rail{} circuit uses C-elements{} as latches. These comprise four simple gates in the
\processlibname{pipistrelle4} library (due to lack of \gate*{AOI32} cells) and a single complex gate
in the \processlibname{umc65ll} library. Note that the cell area
varies dramatically between the libraries due to transistor sizing%
---\processlibname{umc65ll} being minimally-sized for superthreshold
and \processlibname{pipistrelle4} larger for subthreshold operation. The
\emph{number} of cells does not vary significantly.

Latency is measured from \textsc{s\trans{}v}{} in the dual-rail{} designs, and the clock period
defines the latency for single-rail{} designs. The dual-rail{} circuit enables $10\times$ reduction in
average latency thanks to early propagation. Average throughput is worsened
however, due to the lengthened logic path and the need for the additional
\textsc{v\trans{}s}{} transition. Although the dual-rail{} switching power is greater due to
higher inherent activity factor, the computation energy is reduced due to
increased throughput.

Throughput period is defined by the single-rail{} circuit's clock period. For the dual-rail{}
design, throughput period is determined by $\tspcw + \tspcw$ so that
the \acp{PI} are ready for the next operand. \tcwsp{} has the same magnitude as
$\max (\tspcw)$.

\Cref{fig:infer-avg-pipi4} shows the effects of supply voltage on datapath
latency. The latency increases exponentially as the supply voltage is reduced
from \SIrange{0.6}{0.25}{\volt}. The key point is that the circuit functionality
is guaranteed across the whole supply voltage range thanks to the requirements
in \cref{sec:method} and without any alteration to the hardware.

\begin{figure}
	\centering
	\tikzinput{infer-avg-pipi4}
	\caption{%
		Scaling of dual-rail{} datapath latency with supply voltage for the
		\processlibname{pipistrelle4} library.
	}\label{fig:infer-avg-pipi4}
\end{figure}

\section{Introduction}

There is an accelerating demand for connected devices in the
\ac{IoT}~\cite{Capra2019}. Such devices often comprise a sensing aspect,
collecting environmental or personal data, for providing useful monitoring and decisions for transforming our everyday life. The
sensors collect vast amounts of data which must be processed into a usable or
more manageable form. Traditionally this was done by offloading the data into
cloud compute servers, usually over a wireless medium. However this paradigm is
quickly becoming unmaintainable as \ac{IoT} devices expand well into the
billions~\cite{shafik2018real}. Generated data sizes become overwhelming, wireless data transmissions
violate power budgets, and we see a shift towards data processing at the
edge~\cite{arm-trillium}.
Designers of \ac{IoT} products are turning to \ac{ML} in order to extract
meaningful features from the sensed data. Such products are often powered by
batteries or energy harvesters which demand low power and energy efficiency, as
well as robustness to supply variations~\cite{shafik2018real}.

There are several \ac{ML} algorithms which may be suited to such applications,
with \acp{NN} in widespread usage thanks to their often state-of-the-art
accuracy and powerful hardware\slash software ecosystem. \Ac{HDC} has also emerged in recent years with applications in low
power systems~\cite{Burrello2020}. Recently the \ac{TM} algorithm has been proposed as a
promising \ac{ML} algorithm based on \acp{TA}---specialized \acp{LA}. The
\acp{TA} use reinforcement learning locally, together creating an ensemble
learning effect on the global scale which is used to compose logic clauses.
Existing hardware based on \ac{TM} offers a new direction for \ac{ML} whose
inference engine is based on logic with little arithmetic~\cite{Wheeldon2020b}.
The logic-based underpinnings of the \ac{TM} algorithm provide opportunities
for low power and energy efficient \ac{ML} hardware design in the \ac{IoT}.

In this work we apply an asynchronous circuit design
methodology~\cite{Wheeldon2019} to the \ac{TM} algorithm. By removing the pairing between clock and supply voltage as in the synchronous digital designs, it enables an aggressive
voltage scaling~\cite{Diamant2015} for reduced energy per inference and also
adds robustness to environmental variations. Although we use \ac{ML} as the key
application driver, it is possible these techniques can also be applied in other
application areas. Our method is built on \emph{dual-rail{}} circuits with early
propagation~\cite{Brej2006}. Dual-rail{} is an asynchronous circuit design style in the family of
\ac{QDI} circuits. It is inherently robust to circuit delay variations which
means it can operate across a wide range of supply voltages and temperatures.
This usually comes at the cost of duplicated logic and \ac{CD}
overhead~\cite{Yakovlev2013}. In our design we carefully select circuit topology
to minimize such duplications. Additionally we use timing optimizations to reduce
overhead from \ac{CD}.

\subsubsection*{Nomenclature}
Positive- and negative-rail signals are denoted \prail{x} and \nrail{x}
respectively.
\signal[m]{x} denotes the \nth{m} signal in the bit vector \signal{x}.
\textsc{V\trans{}s}{} denotes a transition on a dual-rail{} signal from a valid
codeword to a spacer. Vice versa for \textsc{s\trans{}v}{}.

\subsubsection*{Major Contributions of this paper}
\begin{enumerate}
        \item application of early-propagative, reduced-overhead self-timed
                dual-rail{} circuits to \ac{ML} inference; and
        \item analysis of operand and delay probability distributions in the
                \ac{ML} inference circuit.
\end{enumerate}

\subsubsection*{Paper Organization}
\Cref{sec:tm} introduces the concepts of the \ac{TM} algorithm.
\Cref{sec:method} first briefly introduces the principles of dual-rail{} circuits
before describing our reduced \ac{CD} scheme. \Cref{sec:infer} presents our
dual-rail{} inference datapath design with in-depth analyses. We finally conclude our
findings in \cref{sec:conclusion}.

\section{Self-timing Methodology}\label{sec:method}

In dual-rail{} logic two wires are used to encode a codeword. For a single bit
$\signal{x}$, the dual-rail{} encoding consists of the positive and negative rails
$\{\prail{x},\nrail{x}\}$. $\signal{x}=0$ is encoded as $\{0,1\}$, and
$\signal{x}=1$ is encoded as $\{1,0\}$. One of the remaining states---$\{0,0\}$
or $\{1,1\}$---is chosen to represent the empty state, referred to as a
\emph{spacer}, which separates valid codewords temporally so they can be
distinguished from each other. Care must be taken to correctly handle spacer{}
in the design, otherwise data hazards could occur where one valid{}
overtakes another~\cite{Sokolov2006}. The remaining state is forbidden and must
be avoided by design.

In our design we abide by the following requirements to ensure correct circuit
operation:
\begin{enumerate}
	\item Monotonic switching at the \acsp{PI}.\label[assumption]{ass:monotonic-pi}
	\item Monotonic switching within the circuit.\label[assumption]{ass:monotonic-circuit}
	\item Acknowledgment of \textsc{s\trans{}v}{} on \acsp{PO}.\label[assumption]{ass:cd}
	\item \textsc{V\trans{}s}{} on \acsp{PO} and internal signals before new
		\acsp{PI} applied.\label[assumption]{ass:timing}
	\item \acsp{PI} must transition \textsc{s\trans{}v}{} and \textsc{v\trans{}s}{} for each operand.\label[assumption]{ass:determinism}
	\item \acsp{PI} transition \textsc{v\trans{}s}{} only after \textsc{s\trans{}v}{} on \acsp{PO}.\label[assumption]{ass:pi-cwtosp}
\end{enumerate}

\Cref{ass:monotonic-pi,ass:pi-cwtosp,ass:determinism} are assumed as part of the
circuit's environment. To ensure \cref{ass:monotonic-circuit}, the circuit must
be constructed solely from unate logic gates. To maintain monotonicity we must
exclude non-unate logic gates (\eg{} \gate*{xor} and \gate*{xnor}) from our
library when generating the dual-rail{} netlist.
\Cref{ass:cd} is taken care of by \ac{CD} insertion.
\Cref{ass:timing} can either be assumed as part of the environment or a delay
can be added to the falling edge of \ac{CD} assertion. The latter will be
discussed in \cref{sec:timing-assumption}.

\subsection{Reduced Completion Detection Scheme}\label{sec:timing-assumption}

\ac{CD} which acknowledges both \textsc{s\trans{}v}{} and \textsc{v\trans{}s}{} at the \acp{PO} is expensive to implement due to the
vast amount of complex C-elements{} required~\cite{Sparso2001}. By indicating only \textsc{s\trans{}v}{}
transitions we can significantly reduce the overhead of \ac{CD} by using a small
number of simple gates.

Full \ac{CD} on \emph{internal signals} is even more costly and
removes the possiblity of early propagation. Its job is to ensure \textsc{s\trans{}v}{} and \textsc{v\trans{}s}{} occurs on internal nets for each operand. Internal \ac{CD} can be safely
omitted by giving a grace period for the internal signal to reset to spacer
before applying new \acp{PI}. Codeword validity and correct operation can still
be guaranteed as long as \cref{ass:pi-cwtosp,,ass:determinism,ass:timing} are
met.

In order to meet \cref{ass:timing} there must be a sufficient grace period from
application of spacer{} at the \acp{PI} until application of the next
valid{} at the \acp{PI}. The grace period can be determined by using \ac{STA}
to find the maximum possible \textsc{s\trans{}v}{} time on all nodes of the circuit.
Consequently the grace period can be guaranteed by either
\begin{inparaenum}
	\item the circuit environment waiting for the required grace period; or
	\item an appropriate delay built into the done signal of the
		\ac{CD}.\label[method]{meth:delay}
\end{inparaenum}
The required delay can be calculated as
\newcommand*\td{\ensuremath{t_\mathrm{d}}}%
\newcommand*\tnet{\ensuremath{t_\mathrm{int}}}%
\newcommand*\tio{\ensuremath{t_\mathrm{io}}}%
\newcommand*\tdone{\ensuremath{t_\mathrm{done}}}%
$ \td = \tnet - \tio $,
where \tnet{} is the maximum internal net \textsc{v\trans{}s}{} time, and \tio{}
is the maximum \textsc{v\trans{}s}{} time from the \acp{PI} to \acp{PO}.
\tnet{} must include false paths. It is these false path which lead to the
distinction between \tnet{} and \tio{}. Since there may be some margin added to
\td{}, or due to implementation of the delay \td{} may be greater than the
requirement, the actual timing of the 1\texorpdfstring{\ignorespaces\,$\to$\,\ignorespaces}{ to }{}0 transition of \signal{done} can be calculated by
$ \tdone{}_{1\rightarrow{}0} = \tio + \td $.

\section{Tsetlin Machine Overview}\label{sec:tm}

The main inference component of the \ac{TM} is the conjunctive clause which uses
propositional logic expressions to produce a vote. The composition of each
clause (determined by inclusion of literals) is controlled by the action outputs of a team of
\aclp{TA}.
For inference, the \acp{TA} themselves are not required.
Following a number of reinforcement steps, the automata decide whether their associated literal should be excluded from (action 1) or included in (action 2) the clause.
\Cref{fig:tsetlin} illustrates a \ac{TM} classifier with automaton teams and
conjunctive clauses as one block for brevity.

Each clause
can produce a vote for its class. Half of the clauses can vote
positively, while the other half of the clauses can vote negatively. The inclusion of
inhibition in the voting system enables non-linearity in the inference
process. The votes are summed in a majority vote to produce a collective result
which gives an indication of confidence. This confidence is used to influence future decisions of the automata~\cite{Granmo2018}.

\begin{figure}
    \centering
    \tikzinput{tm-overview}
    \caption{%
        Simplified overview of a single \acf{TM} classifier.
    }\label{fig:tsetlin}
\end{figure}

A simple thresholding function can be used
to generate the final classification output. If the votes are positive (or
zero), the input data is determined to belong to the class. For a negative sum
the input data is determined to be \emph{not} in the class.

For purposes of studying inference, we abstract the \ac{TA} action outputs
to the circuit's environment and concentrate on only the clauses
calculation and majority voting.