Yan, Yang and Fang

*Huazhen Fang,

# High-Order Leader-Follower Tracking Control under Limited Information
Availability††thanks: This work was supported in part by the U.S. National
Science Foundation under Awards CMMI-1763093 and CMMI-1847651.

Chuan Yan  Tao Yang  Huazhen Fang  Department of Mechanical Engineering,
University of Kansas, Kansas, USA  State Key Laboratory of Synthetical
Automation for Process Industries, Northeastern University, Shenyang, China
<EMAIL_ADDRESS>
###### Abstract

[Abstract]Limited information availability represents a fundamental challenge
for control of multi-agent systems, since an agent often lacks sensing
capabilities to measure certain states of its own and can exchange data only
with its neighbors. The challenge becomes even greater when agents are
governed by high-order dynamics. The present work is motivated to conduct
control design for linear and nonlinear high-order leader-follower multi-agent
systems in a context where only the first state of an agent is measured. To
address this open challenge, we develop novel distributed observers to enable
followers to reconstruct unmeasured or unknown quantities about themselves and
the leader and on such a basis, build observer-based tracking control
approaches. We analyze the convergence properties of the proposed approaches
and validate their performance through simulation.

###### keywords:

Leader-follower tracking; multi-agent systems; distributed observer;
distributed control; high-order dynamics

††articletype: Article Type

##  1 Introduction

Cooperative autonomy based on multi-agent systems (MASs) is finding ever-
increasing application in different fields. This has driven a surge of
research on distributed cooperative control for different tasks, including
consensus, leader-follower tracking, synchronization, rendezvous, flocking,
optimization, learning, formation control, and coverage control 1, 2, 3, 4, 5,
6, 7, 8, 9, 10. Most of the current literature considers agents governed by
first- or second-order models. Although such low-order models are useful as
well as amenable to control design, they are often found inadequate to
characterize agents with more complex higher-order dynamics. It is also
neither trivial nor easy to extend low-order cooperative control designs to
high-order systems. Recent years have thus seen a growing amount of work on
high-order MAS control synthesis 11.

The study 12 takes a lead in investigating high-order MAS consensus,
presenting a distributed consensus control algorithm. This subject has since
attracted considerable research efforts, with many studies proposed to deal
with various challenges, e.g., directed communication topologies 13, 14,
switching topologies 13, 15, output feedback design 16, 17, 18, 19, 20, 21,
bipartite consensus 22, 23, external disturbances 24, 25, switching or
heterogeneous dynamics 26, quantization 27, and constrained energy budget 8,
to name a few. Besides consensus, high-order leader-follower tracking has
emerged as another problem of great interest. A basic form of the problem is
introduced in 12, which assumes that the leader agent continuously broadcasts
its state information to all the followers. A consensus-based control
algorithm is then developed therein to make the followers track the leader.
The study 28 considers a more general setting where only a subset of the
followers can receive information from the leader. It proposes a leader-
follower tracking control method and proves that followers with small degrees
must be informed by the leader to ensure tracking convergence. Further,
nonlinear dynamics constitutes a stronger challenge for high-order leader-
follower tracking. The investigations 29, 30 propose to adaptively estimate
the nonlinearity inherent in an agent’s dynamics using neural networks and
then offset it in the control run. A few well-known nonlinear control
techniques, including backstepping control and iterative nonlinear control,
have also found use in addressing tracking problems of nonlinear multi-agent
systems 31, 32, 33. In 34, a finite-time tracking control approach is
developed for a high-order nonlinear MAS with actuator saturation, and the
work in 35 studies the problem of finite-time higher-order tracking with
mismatched disturbances.

Despite the importance, the above studies generally assume that a follower can
obtain a large amount of information to make control decisions. For example,
it is required in 12, 28, 29, 34 that a follower must know all of its own
states, all of the states of its neighbors, and if connected with the leader,
all of the leader’s states. This requirement can be hardly met in a real world
where relevant sensors can be unavailable 28. One can also find similar
requirements in studies about high-order consensus control. This hence
motivates us to explore a more realistic setting—when only the first state of
every agent (leader and followers) is measured. It is unsurprising that the
low information availability will increase the difficulty for tracking
control. In addition, the present literature often requires that the leader’s
dynamics is input-free 29, 36, 30, 32, 31, 33, which comes as another
limitation in practice. To overcome the challenges, we propose to exploit the
notion of observer-based control design and present two major contributions.
First, we design an observer-based tracking control approach for linear high-
order MASs. We propose a set of distributed observers to compensate for the
limited information, allowing a follower to comprehensively estimate the
leader’s maneuver input and states along with its own states. These observers
are then combined with a nominal controller to form an observer-embedded
tracking controller. We further characterize the convergence of tracking when
the proposed controller is applied. As the second contribution, the study
extends to the more challenging case when the agents’ dynamics is not only
high-order but also nonlinear. Extrapolating the design for the linear case,
we develop an observer-based tracking control approach and analyze its
convergence. Compared to the literature, our work is distinct in two aspects.
First, the control design requires very limited information—only every agent’s
first state. Second, the leader’s dynamics is input-driven, and the input is
only known by the followers that communicate with the leader.

Our work is further related with two lines of research. 1) Leader-follower
tracking via output feedback, which generally considers a state-space model
and uses local state or parameter observers to estimate certain unknown
quantities 37, 38, 39, 40, 41, 42. Differing from them, our study designs
distributed observers to help the followers cooperate to infer the leader’s
input and states, removing the restriction that the leader’s maneuver input
must be either zero or known by all the followers. 2) Observer-based first-
and second-order tracking control. The literature includes various kinds of
observers designed to allow a follower to estimate its own velocity 43, 44,
the disturbances acting on it 35, its velocity relative to the leader 45, the
leader’s input 46, 47, 48, or the leader’s velocity 49, 50. These results
nonetheless cannot be readily generalized to the high-order MASs. Finally, the
conference version of this study in 51 considers only agents with linear high-
order dynamics. We make substantial expansion in this paper to investigate
MASs with generic-form nonlinear high-order dynamics.

The rest of the paper is organized as follows. Section 2 introduces the
notation about graph theory. Section 3 formulates the problem of leader-
follower tracking control for a linear high-order MAS, develops an observer-
based tracking control approach, and analyzes its convergence. Section 4
proceeds to study the tracking problem for a nonlinear high-order MAS. Section
5 then offers a numerical simulation example to validate the proposed design.
Finally, Section 6 gathers our conclusions.

##  2 Notation

We use a graph to delineate the information exchange topology among the leader
and followers. First, consider a network composed of $N$ independent
followers, the topology of which is modeled as an undirected graph. The
follower graph is expressed as $\mathcal{G}=(\mathcal{V},\mathcal{E})$, where
$\mathcal{V}=\\{1,2,\cdots,N\\}$ is the node set and
$\mathcal{E}\subseteq\mathcal{V}\times\mathcal{V}$ is the edge set that
contains unordered pairs of nodes. A path is a sequence of connected edges in
a graph. The follower graph is connected if there is a path between every pair
of nodes. The neighbor set of node $i$ is denoted as $\mathcal{N}_{i}$. The
adjacency matrix of $\mathcal{G}$ is $A=[a_{ij}]\in\mathbb{R}^{N\times N}$,
which has non-negative elements. The element $a_{ij}>0$ if and only if
$(i,j)\in\mathcal{E}$, and moreover, $a_{ii}=0$ for all $i\in\mathcal{V}$. For
the Laplacian matrix $L=[l_{ij}]\in\mathbb{R}^{N\times N}$, $l_{ij}=-a_{ij}$
if $i\neq j$ and $l_{ii}=\sum_{k\in\mathcal{N}_{i}}a_{ik}$. The leader is
numbered as node $0$ and can send information to its neighboring followers.
Then, we define a graph
$\bar{\mathcal{G}}=(\bar{\mathcal{V}},\bar{\mathcal{E}})$ for the entire
network, where $\bar{\mathcal{V}}=\\{0\\}\cup\mathcal{V}$, and
$\bar{\mathcal{E}}\subseteq\bar{\mathcal{V}}\times\bar{\mathcal{V}}$ is the
edge set for all nodes. The leader is globally reachable in
$\bar{\mathcal{G}}$ if there is a path in graph $\bar{\mathcal{G}}$ from every
node 0 to node $i$. We denote the leader adjacency matrix associated with
$\bar{\mathcal{G}}$ by $B=\mathrm{diag}(b_{1},\ldots,b_{N})$, where $b_{i}>0$
if the leader is a neighbor of agent $i$ and $b_{i}=0$ otherwise.

##  3 Leader-Follower Tracking with Linear High-Order Dynamics

In this section, we show the leader-follower tracking problem for linear high-
order MASs and develop an observer-based tracking control method as well as
convergence analysis.

###  3.1 Problem Formulation

Consider an MAS composed of $N+1$ agents, among which agent $0$ is the leader
and agents $1$ to $N$ are followers. Each agent has $l$th-order dynamics
($l\geq 3$) expressed as

$\displaystyle\dot{x}_{i,m}$ $\displaystyle=x_{i,m+1},\quad m=1,2,\ldots,l-1,$
(1a) $\displaystyle\dot{x}_{i,l}$ $\displaystyle=u_{i},$ (1b)

for $i=0,1,\ldots,N$, where $x_{i,m}\in\mathbb{R}$ is the $m$th state of agent
$i$, and $u_{i}$ the maneuver input. The objective is to design a distributed
control law $u_{i}$ such that follower $i$ for $i=1,2,\ldots,N$ can
convergently track the leader with
$\lim_{t\to\infty}|x_{i,m}(t)-x_{0,m}(t)|=0$ for $m=1,2,\ldots,l$.

Here, we assume that only $x_{i,1}$ for $i=0,1,\ldots,N$ is available. That
is, only the first state of an agent is measured, regardless of whether it is
the leader or a follower. This assumption considerably relaxes the usual
requirement in the literature that substantial states of an agent must be
measured. However, it also implies that the accessible information about the
agents is rather limited, which makes it more challenging to design an
effective distributed tracking controller.

###  3.2 The Proposed Algorithm

We develop an observer-based control algorithm to enable convergent tracking
in the above setting. To begin with, we consider the following controller for
follower $i$:

$\displaystyle u_{i}=$ $\displaystyle-
k_{1}\left[\sum_{j\in\mathcal{N}_{i}}a_{ij}(x_{i,1}-x_{j,1})+b_{i}(x_{i,1}-x_{0,1})\right]-\sum_{m=2}^{l}k_{m}(\hat{x}_{i,m}-\hat{x}_{0,i,m})+\hat{u}_{0,i},$
(2)

for $i=1,\ldots,N$, where $k_{m}$ for $m=1,2,\ldots,l$ are gain parameters,
$\hat{x}_{0,i,m}$ and $\hat{u}_{0,i}$ are follower $i$’s estimates of the
leader’s state $x_{0,m}$ and input $u_{0}$, respectively, and $\hat{x}_{i,m}$
is follower’s estimate of its own state $x_{i,m}$. The motivation underlying
(2) is to drive follower $i$ toward its neighbors and the leader
simultaneously. When all the followers do this, they can track the leader in a
collective manner. Next, we design the observers so as to obtain the estimates
as needed in (2).

A distributed input observer is first introduced to enable the followers to
estimate $u_{0}$, which is given by

$\displaystyle\dot{\hat{u}}_{0,i}=-\sum_{j\in\mathcal{N}_{i}}a_{ij}(\hat{u}_{0,i}-\hat{u}_{0,j})-b_{i}(\hat{u}_{0,i}-u_{0})-d_{i}\cdot\mathrm{sgn}\left[\sum_{j\in\mathcal{N}_{i}}a_{ij}(\hat{u}_{0,i}-\hat{u}_{0,j})+b_{i}(\hat{u}_{0,i}-u_{0})\right],$
(3a)
$\displaystyle\dot{d}_{i}=\tau_{i}\left|\sum_{j\in\mathcal{N}_{i}}a_{ij}(\hat{u}_{0,i}-\hat{u}_{0,j})+b_{i}(\hat{u}_{0,i}-u_{0})\right|,$
(3b)

for $i=1,\ldots,N$, where $d_{i}$ is an adaptive gain, and $\tau_{i}$ a
positive scalar. For (3a), the first two terms on its right-hand side are used
to drive the input estimation $\hat{u}_{0,i}$ to approach $u_{0}$ while
maintaining consistency with the neighboring followers; further, its third
term serves as further correction with an adaptive gain given in (3b) to
enhance the convergence.

We further propose another distributed observer for the followers to estimate
$x_{0,m}$ for $m=2,3,\ldots,l$:

$\displaystyle\dot{z}_{0,i,2}$
$\displaystyle=-b_{i}c_{0,2}z_{0,i,2}-b_{i}^{2}c_{0,2}^{2}x_{0,1}-c_{0,2}\sum_{j\in\mathcal{N}_{i}}a_{ij}(\hat{x}_{0,i,2}-\hat{x}_{0,j,2})+\hat{x}_{0,i,3},$
(4a) $\displaystyle\hat{x}_{0,i,2}$
$\displaystyle=z_{0,i,2}+b_{i}c_{0,2}x_{0,1},$ (4b)
$\displaystyle\dot{z}_{0,i,m}$
$\displaystyle=-c_{0,m}z_{0,i,m}-c_{0,m}^{2}\hat{x}_{0,i,m-1}+\hat{x}_{0,i,m+1},$
(4c) $\displaystyle\hat{x}_{0,i,m}$
$\displaystyle=z_{0,i,m}+c_{0,m}\hat{x}_{0,i,m-1},\ \ m=3,4,\ldots,l-1,$ (4d)
$\displaystyle\dot{z}_{0,i,l}$
$\displaystyle=-c_{0,l}z_{0,i,l}-c_{0,l}^{2}\hat{x}_{0,i,l-1}+\hat{u}_{0,i},$
(4e) $\displaystyle\hat{x}_{0,i,l}$
$\displaystyle=z_{0,i,l}+c_{0,l}\hat{x}_{0,i,l-1},$ (4f)

for $i=1,\ldots,N$, where $z_{0,i,m}$ and $c_{0,m}$ for $m=2,3,\ldots,l$ are
the observer’s internal states and gain parameters, respectively. The
development of (4) is inspired by 52, in which a centralized disturbance
observer is designed for a single plant. Significantly transforming the
original design, we develop the above observer, which has a distributed
structure that is uniquely suitable for the considered MAS setting. Here, the
observer in (4a)-(4b) attempts to estimate the leader’s second state
$x_{0,2}$. In (4a),
$-\sum_{j\in\mathcal{N}_{i}}a_{ij}(\hat{x}_{0,i,2}-\hat{x}_{0,j,2})$ tries to
keep the estimation consistent between follower $i$ and its neighbors; then,
(4b) performs the estimation to obtain $\hat{x}_{0,i,2}$ using the internal
variable $z_{0,i,2}$ and the leader’s first state $x_{0,1}$. Following similar
lines in (4a)-(4b), the observers in (4c)-(4f) are designed to estimate the
leader’s rest higher-order states.

Finally, we design the following observer such that follower $i$ can estimate
its own states $x_{i,m}$ for $m=2,3,\ldots,l$:

$\displaystyle\dot{z}_{i,2}$
$\displaystyle=-r_{2}z_{i,2}-r_{2}^{2}x_{i,1}+\hat{x}_{i,3},$ (5a)
$\displaystyle\hat{x}_{i,2}$ $\displaystyle=z_{i,2}+r_{2}x_{i,1},$ (5b)
$\displaystyle\dot{z}_{i,m}$
$\displaystyle=-r_{m}z_{i,m}-r_{m}^{2}\hat{x}_{i,m-1}+\hat{x}_{i,m+1},$ (5c)
$\displaystyle\hat{x}_{i,m}$ $\displaystyle=z_{i,m}+r_{m}\hat{x}_{i,m-1},\
m=3,4,\ldots,l-1,$ (5d) $\displaystyle\dot{z}_{i,l}$
$\displaystyle=-r_{l}z_{i,l}-r_{l}^{2}\hat{x}_{i,l-1}+u_{i},$ (5e)
$\displaystyle\hat{x}_{i,l}$ $\displaystyle=z_{i,l}+r_{l}\hat{x}_{i,l-1},$
(5f)

for $i=1,\ldots,N$, where $z_{i,m}$ and $r_{i,m}$ for $m=2,3,\ldots,l$ are the
internal states and gain parameters, respectively.

Putting together (2)-(5), we obtain a distributed observer-based control
algorithm to achieve high-order leader-follower tracking. Its convergence is
analyzed in Section 3.3.

######  Remark 3.1.

We highlight a comparison between the proposed approach and the study of
output-feedback leader-following tracking control in 39, 37, 38, 40, 41, 42.
These references use different observers to help a follower estimate its own
states or certain parameters. These observers are local observers as they are
designed to estimate local unknown quantities. Compared with them, the
proposed approach focuses more on distributed observer design—the observers in
(3)-(4) have a distributed structure to enable the followers to collectively
infer the leader’s input and states by information exchange. This new design
allows the followers to keep tracking the input-driven leader, setting it
apart from the references that restrictively require the leader to be input-
free or the followers to have at least certain knowledge of the leader’s
input. $\bullet$

###  3.3 Convergence Analysis

This section characterizes the convergence property for the algorithm proposed
above. Before proceeding further, we make the following assumption:

The leader’s input $u_{0}\in\mathcal{C}^{1}$ with $|\dot{u}_{0}|\leq
w<\infty$, where $w$ is unknown.

This assumption is reasonable and justifiable due to the fact that control
actuations are usually smooth and subject to ramp-up/down limits.

As the convergence depends on the estimation and tracking errors, we lay out
the definitions of these errors first. For the observer in (3), we define the
estimation error as $e_{u,i}=\hat{u}_{0,i}-u_{0}$. Its dynamics is

$\displaystyle\dot{e}_{u,i}$
$\displaystyle=-b_{i}e_{u,i}-\sum_{j\in\mathcal{N}_{i}}a_{ij}(e_{u,i}-e_{u,j})-d_{i}\cdot\mathrm{sgn}\left[\sum_{j\in\mathcal{N}_{i}}a_{ij}(e_{u,i}-e_{u,j})+b_{i}e_{u,i}\right]-\dot{u}_{0}.$
(6)

Further, define
$e_{u}=\left[\begin{matrix}e_{u,1}&e_{u,2}&\cdots&e_{u,N}\end{matrix}\right]^{\top}$.
It then follows from (6) that

$\displaystyle\dot{e}_{u}=-He_{u}-D\cdot\mathrm{sgn}(He_{u})-\dot{u}_{0}\mathbf{1},$
(7)

where $H=B+L$ and $D=\mathrm{diag}\\{d_{1},\ldots,d_{N}\\}$. It is seen that
the signum function term in (7) is discontinuous, measurable and locally
bounded. Therefore, (7) admits a Filippov solution according to 53, which is
governed by the differential inclusion $\mathcal{K}[\cdot]$:

$\displaystyle\dot{e}_{u}\in^{a.e.}\mathcal{K}\left[-He_{u}-D\cdot\mathrm{sgn}(He_{u})-\dot{u}_{0}\mathbf{1}\right].$
(8)

Now, we consider the observer in (4). Defining
$e_{0x,i,m}=\hat{x}_{0,i,m}-x_{0,m}$, we have

$\displaystyle\dot{e}_{0x,i,2}$
$\displaystyle=-c_{0,2}b_{i}e_{0x,i,2}+e_{0x,i,3}c_{0,2}\sum_{j\in\mathcal{N}_{i}}a_{ij}({e}_{0x,i,2}-{e}_{0x,j,2}),$
(9a) $\displaystyle\dot{e}_{0x,i,m}$
$\displaystyle=-c_{0,m}c_{0,2}b_{i}e_{0x,i,2}+e_{0x,i,m+1}-c_{0,m}c_{0,2}\sum_{j\in\mathcal{N}_{i}}a_{ij}({e}_{0x,i,2}-{e}_{0x,j,2}),\quad
m=3,4,\ldots,l-1,$ (9b) $\displaystyle\dot{e}_{0x,i,l}$
$\displaystyle=-c_{0,l}c_{0,2}b_{i}e_{0x,i,2}+e_{u,i}-c_{0,l}c_{0,2}\sum_{j\in\mathcal{N}_{i}}a_{ij}(e_{0x,i,2}-e_{0x,j,2}).$
(9c)

Define
$e_{0x,m}=\left[\begin{matrix}e_{0x,1,m}&e_{0x,2,m}&\cdots&e_{0x,N,m}\end{matrix}\right]^{\top}$
and
$e_{0x}=\left[\begin{matrix}e_{0x,2}^{\top}&e_{0x,3}^{\top}&\cdots&e_{0x,l}^{\top}\end{matrix}\right]^{\top}$.
Then, (9) can be written into a compact form as below:

$\displaystyle\dot{e}_{0x}=F_{1}e_{0x}+\ell_{1},$ (10)

where

$\displaystyle F_{1}$
$\displaystyle=\left[\begin{matrix}-c_{0,2}H&I&0&\cdots&0\\\
\vdots&0&\ddots&\ddots&\vdots\\\ \vdots&\vdots&\ddots&\ddots&0\\\
-c_{0,l-1}c_{0,2}H&0&\cdots&0&I\\\
-c_{0,l}c_{0,2}H&0&\cdots&\cdots&0\end{matrix}\right],\
\ell_{1}=\left[\begin{matrix}0\\\ \vdots\\\ 0\\\ e_{u}\end{matrix}\right].$

Proceeding further, we define $e_{x,i,m}=\hat{x}_{i,m}-x_{i,m}$ for the
observer in (5) and have

$\displaystyle\dot{e}_{x,i,2}$ $\displaystyle=-r_{2}e_{x,i,2}+{e}_{x,i,3},$
$\displaystyle\dot{e}_{x,i,m}$
$\displaystyle=-r_{m}r_{2}e_{x,i,2}+{e}_{x,i,m+1},\ \ m=3,4,\ldots,l-1,$
$\displaystyle\dot{e}_{x,i,l}$ $\displaystyle=-r_{l}r_{2}e_{x,i,2}.$

Define
$e_{x,m}=\left[\begin{matrix}e_{x,1,m}&e_{x,2,m}&\cdots&e_{x,N,m}\end{matrix}\right]^{\top}$
for $m=2,3,\ldots,l$ and
$e_{x}=\left[\begin{matrix}e_{x,2}^{\top}&e_{x,3}^{\top}&\cdots&e_{x,l}^{\top}\end{matrix}\right]^{\top}$.
Then,

$\displaystyle\dot{e}_{x}=F_{2}e_{x},$ (11)

where

$\displaystyle F_{2}=\left[\begin{matrix}-r_{2}I&I&0&\cdots&0\\\
\vdots&0&\ddots&\ddots&\vdots\\\ \vdots&\vdots&\ddots&\ddots&0\\\
-r_{l-1}r_{2}I&0&\cdots&0&I\\\
-r_{l}r_{2}I&0&\cdots&\cdots&0\end{matrix}\right].$

Finally, to investigate the state tracking error, we define it as
$e_{i,m}=x_{i,m}-x_{0,m}$. The dynamics of $e_{i,m}$ is

$\displaystyle\dot{e}_{i,m}$ $\displaystyle=e_{i,m+1},\ \ m=1,2,\ldots,l-1$
$\displaystyle\dot{e}_{i,l}$
$\displaystyle=-k_{1}\left[\sum_{j\in\mathcal{N}_{i}}a_{ij}(e_{i,1}-e_{j,1})+b_{i}e_{i,1}\right]-\sum_{m=2}^{l}k_{m}e_{i,m}-\sum_{m=2}^{l}k_{m}(e_{x,i,m}-e_{0x,i,m})+e_{u,i}.$

Define
$e_{m}=\left[\begin{matrix}e_{1,m}&e_{2,m}&\cdots&e_{N,m}\end{matrix}\right]^{\top}$
for $m=1,2,\ldots,l$, and
$e=\left[\begin{matrix}e_{1}^{\top}&e_{2}^{\top}&\cdots&e_{l}^{\top}\end{matrix}\right]^{\top}$.
Here, $e$ is the global tracking error with its dynamics governed by

$\displaystyle\dot{e}=F_{3}e+\ell_{3},$ (12)

where

$\displaystyle F_{3}$ $\displaystyle=\left[\begin{matrix}0&I&0&\cdots&0\\\
\vdots&0&\ddots&\ddots&\vdots\\\ \vdots&\vdots&\ddots&\ddots&0\\\
0&0&\cdots&0&I\\\
-k_{1}H&-k_{2}I&\cdots&\cdots&-k_{l}I\end{matrix}\right],\ell_{3}=\left[\begin{matrix}0\\\
\vdots\\\ 0\\\
-\sum_{m=2}^{l}k_{m}(e_{x,m}-e_{0x,m})+e_{u}\end{matrix}\right].$

As a main result, the following theorem outlines the convergence property of
the proposed algorithm.

######  Theorem 3.2.

Suppose that the controller in (2)-(5) is applied to (1). If Assumption 3.3
holds, then

$\displaystyle\lim_{t\to\infty}e_{u}(t)=0.$ (13)

Further, $\lim\limits_{t\to\infty}e_{0x}(t)=0$ if there exist
$c_{0,2},c_{0,3},\ldots,c_{0,l}>0$ such that the polynomials

$\displaystyle
h_{i}(s)=s^{l-1}+c_{0,2}s^{l-2}\lambda_{i}(H)+c_{0,2}\lambda_{i}(H)\sum_{z=0}^{l-3}c_{0,l-z}s^{z}$
(14)

for $i=1,2,\ldots,N$ are Hurwitz stable, and $\lim_{t\to\infty}e_{x}(t)=0$ if
there exist $r_{2},r_{3},\ldots,r_{l}>0$ such that the polynomial

$\displaystyle s^{l-1}+r_{2}s^{l-2}+r_{2}\sum_{z=0}^{l-3}r_{l-z}s^{z}$ (15)

is Hurwitz stable. Finally, the global tracking error
$\lim_{t\to\infty}e(t)=0$ if there exist $k_{m}$ for $m=1,2,\ldots,l$ such
that the polynomials

$\displaystyle s^{l}+k_{1}\lambda_{i}(H)+\sum_{z=2}^{l}s^{z-1}k_{z}$ (16)

for $i=1,2,\ldots,N$ are Hurwitz stable.

Proof: To prove (13), let us consider a Lyapunov functional candidate,
$V_{1}=\bar{V}_{1}(e_{u})+\tilde{V}_{1}(d_{i})$, where

$\displaystyle\bar{V}_{1}(e_{u})=\frac{1}{2}e_{u}^{\top}He_{u},\
\tilde{V}_{1}(d_{i})=\sum_{i=1}^{N}\frac{(d_{i}-\beta)^{2}}{2\tau_{i}},$

with $\beta\geq w$. The set-valued Lie derivative of $\bar{V}_{1}(e_{u})$
denoted by $\mathcal{L}\dot{\bar{V}}_{1}$ along (8) is given by

$\displaystyle\mathcal{L}\dot{\bar{V}}_{1}=\mathcal{K}\left[-e_{u}^{\top}H^{2}e_{u}-e_{u}^{\top}HD\cdot\mathrm{sgn}(He_{u})-e_{u}^{\top}H\dot{u}_{0}\mathbf{1}\right]$
$\displaystyle=\mathcal{K}\Big{[}-\sum_{i=1}^{N}d_{i}\left(\sum_{j\in\mathcal{N}_{i}}a_{ij}(\hat{u}_{0,i}-\hat{u}_{0,j})+b_{i}(\hat{u}_{0,i}-u_{0})\right)^{\top}$
$\displaystyle\quad\cdot\mathrm{sgn}\left(\sum_{j\in\mathcal{N}_{i}}a_{ij}(\hat{u}_{0,i}-\hat{u}_{0,j})+b_{i}(\hat{u}_{0,i}-u_{0})\right)-e_{u}^{\top}H^{2}e_{u}-e_{u}^{\top}H\dot{u}_{0}\mathbf{1}\Big{]}$
$\displaystyle\leq-\sum_{i=1}^{N}d_{i}\left|\sum_{j\in\mathcal{N}_{i}}a_{ij}(\hat{u}_{0,i}-\hat{u}_{0,j})+b_{i}(\hat{u}_{0,i}-u_{0})\right|-e_{u}^{\top}H^{2}e_{u}+w\|He_{u}\|_{1},$

where the fact that $\mathcal{K}[f]=\\{f\\}$ if $f$ is continuous is used. It
then is obtained that $\dot{\bar{V}}_{1}\in\mathcal{L}\dot{\bar{V}}_{1}$ 53.
Hence,

$\displaystyle\dot{V}_{1}$
$\displaystyle=\dot{\bar{V}}_{1}+\dot{\tilde{V}}_{1}=\dot{\bar{V}}_{1}+\sum_{i=1}^{N}\frac{(d_{i}-\beta)\dot{d_{i}}}{\tau_{i}}$
$\displaystyle\leq-\sum_{i=1}^{N}d_{i}\left|\sum_{j\in\mathcal{N}_{i}}a_{ij}(\hat{u}_{0,i}-\hat{u}_{0,j})+b_{i}(\hat{u}_{0,i}-u_{0})\right|$
$\displaystyle\quad+\sum_{i=1}^{N}(d_{i}-\beta)\left|\sum_{j\in\mathcal{N}_{i}}a_{ij}(\hat{u}_{0,i}-\hat{u}_{0,j})+b_{i}(\hat{u}_{0,i}-u_{0})\right|-e_{u}^{\top}H^{2}e_{u}+w\|He_{u}\|_{1}$
$\displaystyle=-e_{u}^{\top}H^{2}e_{u}-(\beta-w)\|He_{u}\|_{1}.$

Note that $H$ is positive definite 54. This, together with $\beta\geq w$,
indicates $\dot{V}_{1}\leq 0$. Thus, $V_{1}(e_{u})$ is non-increasing,
implying that $e_{u}$ and $d_{i}$ are bounded. It follows from (3b) that
$d_{i}$ is monotonically increasing. This means that $d_{i}$ should converge
to some finite value. In the meantime, $V_{1}(e_{u})$ reaches a finite limit
as it is decreasing and lower-bounded by zero. If denoting
$s(t)=\int_{0}^{t}e_{u}^{\top}(\tau)H^{2}e_{u}(\tau)d\tau$, we see that
$s(t)\leq V_{1}(0)-V_{1}(t)$ by integrating $\dot{V}_{1}(e_{u})\leq-
e_{u}^{\top}H^{2}e_{u}$. Hence, $\lim_{t\to\infty}s(t)$ exists and is finite.
Due to the boundedness of $e_{u}$ and $\dot{e}_{u}$, $\ddot{s}$ is also
bounded. This implies that $\dot{s}$ is uniformly continuous. Then,
$\lim_{t\to\infty}\dot{s}(t)=0$ by Barbalat’s Lemma 55. Therefore, we can
obtain $\lim_{t\to\infty}e_{u}(t)=0$.

To prove the asymptotic stability of $e_{0x}$, we use the Schur complement, we
can find out that the characteristic polynomial of $F_{1}$ is
$\prod_{i=1}^{N}h_{i}(s)$, where $h_{i}(s)$ is shown in (14). Note that
$\lim_{t\to\infty}\ell_{1}(t)=0$ due to (13). Hence, we have
$\lim_{t\to\infty}e_{0x}(t)=0$ based on the input-to-state stability (ISS)
theory 55. Following similar lines, we can obtain that
$\lim_{t\to\infty}e_{x}(t)=0$ given (15), and further prove that
$\lim_{t\to\infty}e(t)=0$ if (16) holds. $\bullet$

######  Remark 3.3.

The proposed controller only requires the neighboring followers to interchange
$x_{i,1}$, $\hat{u}_{0,i}$ and $\hat{x}_{0,i,2}$, because the observers can
locally estimate other quantities necessary for control. This greatly reduces
the amount of data to be exchanged between agents and makes the design more
advantageous in terms of communication costs. $\bullet$

##  4 High-Order Tracking for Nonlinear Dynamics

We have investigated leader-follower tracking for a linear high-order MAS in
the previous section. Given the importance of nonlinear leader-follower MASs,
this section moves forward to study the case when an MAS has nonlinear high-
order dynamics. We will develop an observer-based tracking control algorithm
and analyze its convergence properties.

Suppose that agent $i$’s dynamics is governed by

$\displaystyle\dot{x}_{i,m}$
$\displaystyle=x_{i,m+1}+f_{m}(\overline{x_{i,m}}),\quad m=1,2,\ldots,l-1,$
(17a) $\displaystyle\dot{x}_{i,l}$
$\displaystyle=u_{i}+f_{l}(\overline{x_{i,l}}),$ (17b)

for $i=0,1,\ldots,N$, where
$f_{m}(\overline{x_{i,m}}):\mathbb{R}^{m}\rightarrow\mathbb{R}$ for
$m=1,2,\ldots,l$ are nonlinear functions with
$\overline{x_{i,m}}=(x_{i,1},x_{i,2},\ldots,x_{i,m})$. Following Section 1, we
assume that only $x_{i,1}$ is measured and continue to hold Assumption 3.3.
The control design objective here is still to enable convergent tracking,
i.e., $\lim_{t\to\infty}|x_{i,m}(t)-x_{0,m}(t)|=0$ for $m=1,2,\ldots,l$ and
$i=1,2,\ldots,N$.

To achieve the above objective, we propose the following distributed
controller:

$\displaystyle
u_{i}=-k_{1}(x_{i,1}-\hat{x}_{0,i,1})-\sum_{m=2}^{l}k_{m}(\hat{x}_{i,m}-\hat{x}_{0,i,m})+\hat{u}_{0,i}.$
(18)

This controller must be supplemented by corresponding observers. It is noted
first that the observer in (3) can also be applied to obtain $\hat{u}_{0,i}$
here, so we continue to use it for the distributed estimation of $u_{0}$. We
then construct the following state observer to allow follower $i$ to estimate
$x_{0,m}$ for $m=1,2,\ldots,l$:

$\displaystyle\dot{\hat{x}}_{0,i,1}$
$\displaystyle=-c_{0,1}\left[\sum_{j\in\mathcal{N}_{i}}a_{ij}(\hat{x}_{0,i,1}-\hat{x}_{0,j,1})+b_{i}(\hat{x}_{0,i,1}-x_{0,1})\right]+\hat{x}_{0,i,2}+f_{1}(\overline{\hat{x}_{0,i,1}}),$
(19a) $\displaystyle\dot{z}_{0,i,2}$
$\displaystyle=-b_{i}c_{0,2}z_{0,i,2}-b_{i}^{2}c_{0,2}^{2}x_{0,1}+\hat{x}_{0,i,3}+f_{2}(\overline{\hat{x}_{0,i,2}})$
$\displaystyle\quad-
c_{0,2}\sum_{j\in\mathcal{N}_{i}}a_{ij}(\hat{x}_{0,i,2}-\hat{x}_{0,j,2})-b_{i}c_{0,2}f_{1}(x_{0,1}),$
(19b) $\displaystyle\hat{x}_{0,i,2}$
$\displaystyle=z_{0,i,2}+b_{i}c_{0,2}x_{0,1},$ (19c)
$\displaystyle\dot{z}_{0,i,m}$
$\displaystyle=-c_{0,m}z_{0,i,m}-c_{0,m}^{2}\hat{x}_{0,i,m-1}+\hat{x}_{0,i,m+1}+f_{m}(\overline{\hat{x}_{0,i,m}})-c_{0,m}f_{m-1}(\overline{\hat{x}_{0,i,m-1}}),$
(19d) $\displaystyle\hat{x}_{0,i,m}$
$\displaystyle=z_{0,i,m}+c_{0,m}\hat{x}_{0,i,m-1},\quad m=3,4,\ldots,l-1,$
(19e) $\displaystyle\dot{z}_{0,i,l}$
$\displaystyle=-c_{0,l}z_{0,i,l}-c_{0,l}^{2}\hat{x}_{0,i,l-1}+\hat{u}_{0,i}+f_{l}(\overline{\hat{x}_{0,i,l}})-c_{0,l}f_{l-1}(\overline{\hat{x}_{0,i,l-1}}),$
(19f) $\displaystyle\hat{x}_{0,i,l}$
$\displaystyle=z_{0,i,l}+c_{0,l}\hat{x}_{0,i,l-1},$ (19g)

where
$\overline{\hat{x}_{0,i,m}}=(\hat{x}_{0,i,1},\hat{x}_{0,i,2},\ldots,\hat{x}_{0,i,m})$.
To make a follower able to estimate its own states, we develop an observer as
follows:

$\displaystyle\dot{z}_{i,2}$
$\displaystyle=-r_{2}z_{i,2}-r_{2}^{2}x_{i,1}+\hat{x}_{i,3}+f_{2}(x_{i,1},\hat{x}_{i,2})-r_{2}f_{1}(x_{i,1}),$
(20a) $\displaystyle\hat{x}_{i,2}$ $\displaystyle=z_{i,2}+r_{2}x_{i,1},$ (20b)
$\displaystyle\dot{z}_{i,m}$
$\displaystyle=-r_{m}z_{i,m}-r_{m}^{2}\hat{x}_{i,m-1}+\hat{x}_{i,m+1}+f_{m}(x_{i,1},\overline{\hat{x}_{i,m}})-r_{m}f_{m-1}(x_{i,1},\overline{\hat{x}_{i,m-1}}),$
(20c) $\displaystyle\hat{x}_{i,m}$
$\displaystyle=z_{i,m}+r_{m}\hat{x}_{i,m-1},\quad m=3,4,\ldots,l-1,$ (20d)
$\displaystyle\dot{z}_{i,l}$
$\displaystyle=-r_{l}z_{i,l}-r_{l}^{2}\hat{x}_{i,l-1}+u_{i}+f_{l}(x_{i,1},\overline{\hat{x}_{i,l}})-r_{l}f_{l-1}(x_{i,1},\overline{\hat{x}_{i,l-1}}),$
(20e) $\displaystyle\hat{x}_{i,l}$
$\displaystyle=z_{i,l}+r_{l}\hat{x}_{i,l-1},$ (20f)

where
$\overline{\hat{x}_{i,m}}=(\hat{x}_{i,2},\hat{x}_{i,3},\ldots,\hat{x}_{i,m})$.

Integrating the above observers in (3), (19)-(20) into the controller in (18)
will yield a complete observer-based tracking controller. Before going further
to analyze its effectiveness, we make the following assumption:

There exist $\rho_{m}\geq 0$ such that

$\displaystyle|f_{m}(\xi)-f_{m}(\epsilon)|\leq\rho_{m}\|\xi-\epsilon\|,\
m=1,2,\ldots,l,$

where $\xi,\epsilon\in\mathbb{R}^{m}$.

Assumption 4 implies that the nonlinear functions must be of Lipschitz class.
It is commonly used in the literature on nonlinear MAS control and can be
satisfied by many practical systems.

The following theorem shows the main result about convergence of the proposed
controller.

######  Theorem 4.1.

Assume that Assumptions 3.3 and 4 hold and that the controller proposed above
is applied to (17). The state tracking error converges to zero, i.e.,
$\lim_{t\rightarrow\infty}|x_{i,m}(t)-x_{0,m}(t)|=0$ for $m=1,2,\ldots,l$ and
$i=1,2,\ldots,N$, if there exist $c_{0,m}$, $r_{n}$ and $k_{m}$ for
$m=1,2,\ldots,l$ and $n=2,\ldots,N$ such that the polynomials (14), (15) and

$\displaystyle\left(s^{l}+\sum_{z=1}^{l}s^{z-1}k_{z}\right)^{N}$ (21)

are Hurwitz stable, and if there exist matrices $Q_{i}>0$ and $\eta_{i}>0$ for
$i=1,2,3$ such that

$\displaystyle F_{4}^{\top}Q_{1}+Q_{1}F_{4}=-\eta_{1}I,$ (22a) $\displaystyle
F_{2}^{\top}Q_{2}+Q_{2}F_{2}=-\eta_{2}I,$ (22b) $\displaystyle
F_{6}^{\top}Q_{3}+Q_{3}F_{6}=-\eta_{3}I,$ (22c)
$\displaystyle\sum_{i=1}^{l}\|P_{0x,i}\|<\min\left\\{\frac{\eta_{1}}{2\|Q_{1}\|},\frac{\eta_{3}}{2\|Q_{3}\|}\right\\},$
(22d) $\displaystyle\sum_{i=2}^{l}\|P_{x,i}\|<\frac{\eta_{2}}{2\|Q_{2}\|},$
(22e)

where
$P_{0x,i}=\mathrm{diag}\\{\rho_{1}I,\rho_{2}I,\ldots,\rho_{i}I,0I,0I,\ldots,0I\\}$,
$P_{x,i}=\mathrm{diag}\\{\rho_{2}I,\rho_{3}I,\ldots,\rho_{i}I,0I,0I,\ldots,0I\\}$
for $i=1,2,\ldots,l$, and

$\displaystyle F_{4}$
$\displaystyle=\left[\begin{matrix}-c_{0,1}H&I&0&\cdots&\cdots&0\\\
0&-c_{0,2}H&I&0&\cdots&0\\\ \vdots&\vdots&0&\ddots&\ddots&\vdots\\\
\vdots&\vdots&\vdots&\ddots&\ddots&0\\\
\vdots&-c_{0,l-1}c_{0,2}H&0&\cdots&0&I\\\
0&-c_{0,l}c_{0,2}H&0&\cdots&\cdots&0\end{matrix}\right],F_{6}=\left[\begin{matrix}0&I&0&\cdots&0\\\
\vdots&0&\ddots&\ddots&\vdots\\\ \vdots&\vdots&\ddots&\ddots&0\\\
0&0&\cdots&0&I\\\ -k_{1}I&-k_{2}I&\cdots&\cdots&-k_{l}I\end{matrix}\right].$

Proof: Let us define
$e_{0x}=\left[\begin{matrix}e_{0x,1}^{\top}&e_{0x,2}^{\top}&\cdots&e_{0x,l}^{\top}\end{matrix}\right]^{\top}$,
where $e_{0x,m}$ follows the same definition as in Section 3, and define
$f_{0x,m}=\left[\begin{matrix}\cdots&f_{m}(\overline{\hat{x}_{0,i,m}})-f_{m}(\overline{x_{0,m}})&\cdots&\end{matrix}\right]^{\top}$
for $i=1,2,\cdots,N$. By (17) and (19), we have

$\displaystyle\dot{e}_{0x}=F_{4}e_{0x}+\ell_{4},$ (23)

where
$\ell_{4}=\left[\begin{matrix}f_{0x,1}^{\top}&\cdots&f_{0x,l-1}^{\top}&f_{0x,l}^{\top}+e_{u}^{\top}\end{matrix}\right]^{\top}$.
We choose a Lyapunov candidate function

$\displaystyle V_{2}(e_{0x})=\frac{1}{2}e_{0x}^{\top}Q_{1}e_{0x},$

for which there exist $\alpha_{1},\alpha_{2}>0$ such that

$\displaystyle\alpha_{1}\|e_{0x}\|^{2}\leq
V_{2}(e_{0x})\leq\alpha_{2}\|e_{0x}\|^{2}.$

By (22a), we have

$\displaystyle\dot{V_{2}}$
$\displaystyle=\frac{1}{2}e_{0x}^{\top}(Q_{1}F_{4}+F_{4}^{\top}Q_{1})e_{0x}+e_{0x}^{\top}Q_{1}\ell_{4}$
$\displaystyle\leq-\frac{1}{2}\eta_{1}\|e_{0x}\|^{2}+\|e_{0x}\|\|Q_{1}\|\|\ell_{4}\|$
$\displaystyle\leq-\frac{1}{2}\eta_{1}\|e_{0x}\|^{2}+\|e_{0x}\|\|Q_{1}\|\big{(}\|P_{0x,1}e_{0x}\|+\|P_{0x,2}e_{0x}\|$
$\displaystyle\quad+\cdots+\|P_{0x,l}e_{0x}\|+\|e_{u}\|\big{)}$
$\displaystyle=-\left(\frac{1}{2}\eta_{1}-\|Q_{1}\|\sum_{i=1}^{l}\|P_{0x,i}\|\right)\|e_{0x}\|^{2}+\|e_{0x}\|\|Q_{1}\|\|e_{u}\|,$

Define

$\displaystyle\sigma_{1}=\frac{1}{2}\eta_{1}-\|Q_{1}\|\sum_{i=1}^{l}\|P_{0x,i}\|,\
\mathcal{X}(\|e_{u}\|)=\frac{\|Q_{1}\|\|e_{u}\|}{\sigma_{1}\theta_{1}}$

for any $0<\theta_{1}<1$. By (22d), one can see that $\sigma_{1}>0$. It can be
verified that

$\displaystyle\|e_{0x}\|\geq\mathcal{X}(\|e_{u}\|)\Rightarrow\dot{V_{2}}\leq-\sigma_{1}(1-\theta_{1})\|e_{0x}\|^{2}.$

Hence, $V_{2}$ is an ISS-Lyapunov function, implying that the system (23) is
ISS 55. Then, we have $\lim_{t\to\infty}e_{0x}=0$ since
$\lim_{t\to\infty}e_{u}=0$ as indicated in (13).

Define
$f_{x,m}=\left[\begin{matrix}\cdots&f_{m}(x_{i,1},\overline{\hat{x}_{i,m}})-f_{m}(x_{i,1},\overline{x_{i,m}})&\cdots\end{matrix}\right]^{\top}$
for $i=1,2,\ldots,N$, and continue to adopt $e_{x}$ as defined in (11).
According to (20), its dynamics can be expressed as

$\displaystyle\dot{e}_{x}=F_{2}e_{x}+\ell_{5},$ (24)

where $F_{2}$ was defined in (11), and
$\ell_{5}=\left[\begin{matrix}f_{x,2}^{\top}&\cdots&f_{x,l}^{\top}\end{matrix}\right]^{\top}$.
Following similar lines to the above, we can prove that
$\lim_{t\to\infty}e_{x}=0$ if (22b) and (22e) hold.

We proceed to consider the global tracking error when the controller in (18)
is applied. We define
$f_{m}=\left[\begin{matrix}\cdots&f_{m}(\overline{x_{i,m}})-f_{m}(\overline{x_{0,m}})&\cdots\end{matrix}\right]^{\top}$
for $i=1,2,\cdots,N$. The dynamics of the tracking error
$e_{i,m}={x}_{i,m}-{x}_{0,m}$ is

$\displaystyle\dot{e}_{i,m}$
$\displaystyle=e_{i,m+1}+f_{m}(\overline{x_{i,m}})-f_{m}(\overline{x_{0,m}}),$
(25a) $\displaystyle\dot{e}_{i,l}$
$\displaystyle=-\sum_{m=1}^{l}k_{m}e_{i,m}-\sum_{m=2}^{l}k_{m}e_{x,m}+\sum_{m=1}^{l}k_{m}e_{0x,m}+e_{u,i}+f_{l}(\overline{x_{i,l}})-f_{l}(\overline{x_{0,l}}),$
(25b)

for $m=1,2,\ldots,l-1$ and $i=1,2,\ldots,N$. The notation of $e$ in (12) is
still adopted here. Now, combining (19), (20) and (25), the closed-loop system
is indicated into a compact structure as below:

$\displaystyle\dot{e}=F_{6}e+\ell_{6}+\ell_{7},$ (26)

where
$\ell_{6}=\left[\begin{matrix}f_{1}^{\top}&\cdots&f_{l}^{\top}\end{matrix}\right]^{\top}$
and

$\displaystyle\ell_{7}=\left[\begin{matrix}0\\\ \vdots\\\ 0\\\
-\sum_{m=2}^{l}k_{m}e_{x,m}+\sum_{m=1}^{l}k_{m}e_{0x,m}+e_{u}\end{matrix}\right].$

For (26), we can use the ISS theory to prove that it is asymptotically stable
if (22c)-(22d) hold. Therefore, $\lim_{t\to\infty}e(t)=0$ is established as
$\lim_{t\to\infty}\ell_{7}(t)=0$ due to (7), (23) and (24). We conclude that
$\lim_{t\to\infty}|x_{i,m}(t)-x_{0,m}(t)|=0$. Finally, we highlight that the
Hurwitz stable polynomials (14), (15) and (21) would ensure the existence of
solutions for the Lyapunov equations in (22a)-(22c). This concludes the proof.
$\bullet$

######  Remark 4.2.

Theorem 4.1 can be briefly explained as follows. The conditions (22a) and
(22d) ensure that the distributed observers in (19) are asymptotically stable;
the conditions (22b) and (22e) ensure that the local state observers in (20)
are asymptotically stable; further, the conditions (22c) and (22d) ensure that
the global closed-loop tracking errors converge to zero. When the polynomials
(14), (15) and (21) are Hurwitz stable, $F_{4}$, $F_{2}$, and $F_{6}$ will be
stable, and then the Lyapunov equations in (22a)-(22c) will admit solutions.

######  Remark 4.3.

The above design can be extended to the case when the nonlinear functions
$f_{m}(\cdot)$ is unknown but admits approximation by a known function with
bounded error. Specifically, suppose that there exist $g_{m}(\cdot)$ and
$\varphi_{m}\geq 0$ such that

$\displaystyle\left|g_{m}(\phi)-f_{m}(\phi)\right|\leq\varphi_{m},\
m=1,2,\ldots,l,$

for any $\phi\in\mathbb{R}^{m}$. We then can replace $f_{m}(\cdot)$ in
(19)-(20) by $g_{m}(\cdot)$ and obtain a tracking controller based on the
approximate nonlinearity. It can be proven that this controller will lead to
bounded-error tracking under certain mild conditions. The analysis is omitted
here for the sake of space.

Figure 1: Topology of the MAS in the simulation.

Figure 2: Third-order nonlinear MAS profiles: (a) leader’s and followers’
state trajectory profiles of $x_{i,1}$ for $i=0,1,\ldots,N$; (b) leader’s and
followers’ state trajectory profiles of $x_{i,2}$ for $i=0,1,\ldots,N$; (c)
leader’s and followers’ state trajectory profiles of $x_{i,3}$ for
$i=0,1,\ldots,N$; (d) leader’s input profile and the estimation by each
follower; (e) leader’s state trajectory profile of $x_{0,1}$ and the
estimation by each follower; (f) leader’s state trajectory profile of
$x_{0,2}$ and the estimation by each follower; (g) leader’s state trajectory
profile of $x_{0,3}$ and the estimation by each follower; (h) followers’
estimation of their own state trajectories of $x_{i,2}$ for $i=1,2,\ldots,N$;
(i) followers’ estimation of their own state trajectories of $x_{i,3}$ for
$i=1,2,\ldots,N$.

##  5 NUMERICAL STUDY

This section presents a numerical simulation example to show the effectiveness
of the proposed design. For the sake of space, we only illustrate the case of
nonlinear leader-follower tracking. Consider a third-order MAS including one
leader and five followers. The agents interchange information based on a
communication topology shown in Fig. 1. Here, node 0 is the leader, and nodes
1 to 5 are followers. The leader transmits data to only follower 1, and the
followers maintain bidirectional communication with their neighbors. The
agents’ dynamics is as described in (17), for which
$f_{m}(\overline{x_{i,m}})=\cos(\overline{x_{i,m}})^{\top}\mathbf{1}=\sum_{k=1}^{m}\cos(x_{i,k})$.
The leader’s maneuver input is set to be $u_{0}=\sin(0.2\pi t)$. When
implementing the proposed observer-based controller, we select
$c_{0,1}=c_{0,2}=c_{0,3}=5$, $r_{2}=r_{3}=4$ and $k_{1}=k_{2}=k_{3}=3$. Such a
gain set is verifiable to make the convergence conditions satisfied. The
simulation results are summarized in Fig. 2. Figs. 2-2 illustrate followers’
and the leader’s state trajectories, showing that the followers can manage to
catch up with and then keep tracking the leader, despite they differ in
initial states. Fig. 2 shows the estimation of the leader’s input by the
followers. For each follower, the estimation can quickly converge to the
actual values. Meanwhile, the followers can also effectively estimate the
leader’s states using the designed observer, with the estimation errors
approaching zero as shown in Figs. 2-2. Figs. 2 and 2 further present the
followers’ estimation of their own unmeasured states. These results validate
that the proposed design can ensure convergent tracking, despite the
nonlinearity and limited information availability.

##  6 Conclusion

We studied leader-follower tracking control for high-order MASs in this paper.
While this problem has recently attracted growing attention, the previous
studies generally require all the states of an agent to be measured, even
though the measurement information can be practically limited by the
availability of sensors. Here, we focused on the challenging but more
realistic setting where only the first state of an agent is measured. We
designed novel distributed observers, by which a follower can reconstruct
unknown or unmeasured quantities about itself and the leader, and then
performed distributed observer-based controller synthesis. We conducted the
design for both linear and nonlinear MASs and characterized the convergence
properties. A simulation result demonstrated the effectiveness of our design.
Our future work will be directed towards extending the results to directed
graphs and completely unknown nonlinearity.

## References

  * 1 Sun S, Zhang Y, Lin P, Ren W, Farrell JA. Distributed Time-Varying Optimization With State-Dependent Gains: Algorithms and Experiments. IEEE Transactions on Control Systems Technology 2021\.
  * 2 Shamma J. Cooperative Control of Distributed Multi-Agent Systems.  New York, NY, USA: Wiley-Interscience .  2008\.
  * 3 Wang Y, Garcia E, Casbeer D, Zhang F. , eds.Cooperative Control of Multi-Agent Systems: Theory and Applications.  Wiley .  2017\.
  * 4 Ren W, Cao Y. Distributed Coordination of Multi-agent Networks: Emergent Problems, Models, and Issues.  London: Springer-Verlag .  2011\.
  * 5 Qin J, Ma Q, Shi Y, Wang L. Recent advances in consensus of multi-agent systems: A brief survey. IEEE Transactions on Industrial Electronics 2017; 64(6): 4972-4983.
  * 6 Cortés J, Egerstedt M. Coordinated control of multi-robot systems: A survey. SICE Journal of Control, Measurement, and System Integration 2017; 10(6): 495-503.
  * 7 Yang T, Yi X, Wu J, et al. A survey of distributed optimization. Annual Reviews in Control 2019; 47: 278-305.
  * 8 Yan C, Fang H, Chao H. Energy-aware leader-follower tracking control for electric-powered multi-agent systems. Control Engineering Practice 2018; 79: 209-218.
  * 9 Mou S, Cao M, Stephen Morse A. Target-point formation control. Automatica 2015; 61: 113-118.
  * 10 Murray RM. Recent research in cooperative control of multivehicle systems. Journal of Dynamic Systems, Measurement, and Control 2007; 129(9): 571-583.
  * 11 Huang J, Fang H, Dou L, Chen J. An overview of distributed high-order multi-agent coordination. IEEE/CAA Journal of Automatica Sinica 2014; 1(1): 1-9.
  * 12 Ren W, Moore KL, Chen Y. High-order and model reference consensus algorithms in cooperative control of multivehicle systems. Journal of Dynamic Systems, Measurement, and Control 2007; 129(5): 678–688.
  * 13 Wen G, Hu G, Yu W, Cao J, Chen G. Consensus tracking for higher-order multi-agent systems with switching directed topologies and occasionally missing control inputs. Systems & Control Letters 2013; 62(12): 1151-1158.
  * 14 Abdessameud A, Tayebi A. Distributed consensus algorithms for a class of high-order multi-agent systems on directed graphs. IEEE Transactions on Automatic Control, 2018\.  in press.
  * 15 Jiang F, Wang L. Consensus seeking of high-order dynamic multi-agent systems with fixed and switching topologies. International Journal of Control 2010; 83(2): 404-420.
  * 16 Seo JH, Shim H, Back J. Consensus of high-order linear systems using dynamic output feedback compensator: Low gain approach. Automatica 2009; 45(11): 2659-2664.
  * 17 Xi J, Shi Z, Zhong Y. Output consensus analysis and design for high-order linear swarm systems: Partial stability method. Automatica 2012; 48(9): 2335-2343.
  * 18 Qian Y, Liu L, Feng G. Output Consensus of Heterogeneous Linear Multi-Agent Systems with Adaptive Event-Triggered Control. IEEE Transactions on Automatic Control 2019; 64(6): 2606-2613.
  * 19 You K, Xie L. Coordination of discrete-time multi-agent systems via relative output feedback. International Journal of Robust and Nonlinear Control 2011; 21(13): 1587-1605.
  * 20 Zhou B, Xu C, Duan G. Distributed and truncated reduced-order observer based output feedback consensus of multi-agent systems. IEEE Transactions on Automatic Control 2014; 59(8): 2264–2270.
  * 21 Wen G, Zhao Y, Duan Z, Yu W, Chen G. Containment of higher-order multi-leader multi-agent systems: A dynamic output approach. IEEE Transactions on Automatic Control 2015; 61(4): 1135–1140.
  * 22 Valcher ME, Misra P. On the consensus and bipartite consensus in high-order multi-agent dynamical systems with antagonistic interactions. Systems & Control Letters 2014; 66: 94-103.
  * 23 Hu J, Wu Y, Liu L, Feng G. Adaptive bipartite consensus control of high-order multiagent systems on coopetition networks. International Journal of Robust and Nonlinear Control 2018; 28(7): 2868-2886.
  * 24 Liu Y, Jia Y. Consensus problem of high-order multi-agent systems with external disturbances: An $H_{\infty}$ analysis approach. International Journal of Robust and Nonlinear Control 2010; 20(14): 1579-1593.
  * 25 Wang X, Li S, Lam J. Distributed active anti-disturbance output consensus algorithms for higher-order multi-agent systems with mismatched disturbances. Automatica 2016; 74: 30-37.
  * 26 Mu B, Zhang K, Shi Y. Integral sliding mode flight controller design for a quadrotor and the application in a heterogeneous multi-agent system. IEEE Transactions on Industrial Electronics 2017; 64(12): 9389-9398.
  * 27 Li X, Chen MZQ, Su H. Quantized consensus of multi-agent networks with sampled data and Markovian interaction links. IEEE Transactions on Cybernetics 2019; 49(5): 1816-1825.
  * 28 Yu W, Chen G, Ren W, Kurths J, Zheng WX. Distributed higher order consensus protocols in multiagent dynamical systems. IEEE Transactions on Circuits and Systems I: Regular Papers 2011; 58(8): 1924-1932.
  * 29 Zhang H, Lewis FL. Adaptive cooperative tracking control of higher-order nonlinear systems with unknown dynamics. Automatica 2012; 48(7): 1432–1439.
  * 30 Cui R, Ren B, Ge SS. Synchronised tracking control of multi-agent system with high-order dynamics. IET Control Theory & Applications 2012; 6(5): 603–614.
  * 31 Li W, Zhang JF. Distributed practical output tracking of high-order stochastic multi-agent systems with inherent nonlinear drift and diffusion terms. Automatica 2014; 50(12): 3231–3238.
  * 32 Jin X. Adaptive iterative learning control for high-order nonlinear multi-agent systems consensus tracking. Systems & Control Letters 2016; 89: 16–23.
  * 33 Sun C, Hu G, Xie L. Robust consensus tracking for a class of high-order multi-agent systems. International Journal of Robust and Nonlinear Control 2016; 26(3): 578–598.
  * 34 Mondal S, Su R. Finite time tracking control of higher order nonlinear multi agent systems with actuator saturation. Proceedings of 14th IFAC Symposium on Control in Transportation Systems 2016; 49(3): 165-170.
  * 35 Wang X, Li S, Yu X, Yang J. Distributed active anti-disturbance consensus for leader-follower higher-order multi-agent systems with mismatched disturbances. IEEE Transactions on Automatic Control 2017; 62(11): 5795-5801.
  * 36 Wen G, Yu W, Li Z, Yu X, Cao J. Neuro-Adaptive Consensus Tracking of Multiagent Systems With a High-Dimensional Leader. IEEE Transactions on Cybernetics 2017; 47(7): 1730-1742.
  * 37 Chu H. Observer-based adaptive consensus tracking for linear multi-agent systems with input saturation. IET Control Theory & Applications 2015; 9: 2124-2131(7).
  * 38 Sun J, Geng Z, Lv Y. Adaptive output feedback consensus tracking for heterogeneous multi-agent systems with unknown dynamics under directed graphs. Systems & Control Letters 2016; 87: 16 - 22.
  * 39 Sun J, Geng Z. Adaptive output feedback consensus tracking for linear multi-agent systems with unknown dynamics. International Journal of Control 2015; 88(9): 1735-1745.
  * 40 Shi P, Shen QK. Observer-based leader-following consensus of uncertain nonlinear multi-agent systems. International Journal of Robust and Nonlinear Control 2017; 27(17): 3794–3811.
  * 41 Zhang H, Lewis FL, Das A. Optimal design for synchronization of cooperative systems: state feedback, observer and output feedback. IEEE Transactions on Automatic Control 2011; 56(8): 1948-1952.
  * 42 Peng Z, Wang D, Zhang H, Sun G. Distributed neural network control for adaptive synchronization of uncertain dynamical multiagent systems. IEEE Transactions on Neural Networks and Learning Systems 2014; 25(8): 1508-1519.
  * 43 Zhang B, Jia Y, Matsuno F. Finite-time observers for multi-agent systems without velocity measurements and with input saturations. Systems & Control Letters 2014; 68: 86-94.
  * 44 Xu X, Chen S, Gao L. Observer-based consensus tracking for second-order leader-following nonlinear multi-agent systems with adaptive coupling parameter design. Neurocomputing 2015; 156: 297-305.
  * 45 Hu J, Geng J, Zhu H. An observer-based consensus tracking control and application to event-triggered tracking. Communications in Nonlinear Science and Numerical Simulation 2015; 20(2): 559–570.
  * 46 Yan C, Fang H. A new encounter between leader–follower tracking and observer-based control: Towards enhancing robustness against disturbances. Systems & Control Letters 2019; 129: 1-9.
  * 47 Yan C, Fang H. Leader-follower tracking control for multi-agent systems based on input observer design. Proceedings of the Annual American Control Conference 2018; 1(1): 478–483.
  * 48 Yan C, Fang H. Observer-based distributed leader-follower tracking control: a new perspective and results. International Journal of Control 2021; 94(1): 39-48.
  * 49 Hong Y, Chen G, Bushnell L. Distributed observers design for leader-following control of multi-agent networks. Automatica 2008; 44(3): 846–850.
  * 50 Cheng Y, Xie D. Distributed observer design for bounded tracking control of leader-follower multi-agent systems in a sampled-data setting. International Journal of Control 2014; 87(1): 41-51.
  * 51 Yan C, Fang H. Observer-based leader-follower tracking control for high-order multi-agent systems with limited measurement information. Proceedings of the IEEE 58th Conference on Decision and Control 2019; 1(1): 3904-3909.
  * 52 Yang J, Li S, Yu X. Sliding-mode control for systems with mismatched uncertainties via a disturbance observer. IEEE Transactions on Industrial Electronics 2013; 60(1): 160–169.
  * 53 Clarke F. Optimization and Nonsmooth Analysis.  John Wiley & Sons New York, NY, USA .  1983\.
  * 54 Hu J, Hong Y. Leader-following coordination of multi-agent systems with coupling time delays. Physica A: Statistical Mechanics and its Applications 2007; 374(2): 853–863.
  * 55 Khalil HK. Noninear Systems.  Prentice-Hall, New Jersey .  2002\.


# TeV and keV-MeV Excesses as Probes for Hadronic Process in BL Lacertaes

Ji-Gui Cheng  Guangxi Key Laboratory for Relativistic Astrophysics, School of
Physical Science and Technology, Guangxi University, Nanning 530004, China
Xiao-Li Huang∗ Guizhou Provincial Key Laboratory of Radio Astronomy and Data
Processing, School of Physics and Electronic Science, Guizhou Normal
Univeristy, Guiyang 550025, China  Ze-Rui Wang  College of Physics and
Electronic Engineering, Qilu Normal University, Jinan 250200, China  Jian-Kun
Huang  Guangxi Key Laboratory for Relativistic Astrophysics, School of
Physical Science and Technology, Guangxi University, Nanning 530004, China
En-Wei Liang∗ Guangxi Key Laboratory for Relativistic Astrophysics, School of
Physical Science and Technology, Guangxi University, Nanning 530004, China

###### Abstract

A hard TeV $\gamma$-ray component excess over the single-zone leptonic model
prediction (TeV excess) is observed in the spectral energy distributions
(SEDs) of some BL Lacs. Its origin is uncertain. We revisit this issue with
four BL Lacs (1ES 0229+200, 1ES 0347–121, 1ES 1101–232, and H2356–309), in
which the TeV excess is detected in their intrinsic SEDs. We represent their
SEDs with a single-zone leptohadronic model, where radiations of the electrons
and protons as well as the cascade electrons produced by the $\gamma\gamma$
and p$\gamma$ interactions within their jets are considered. We show that the
observed SEDs below the GeV gamma-ray band are attributed to the synchrotron
radiations and self-Compton process of the primary electrons, and the TeV
excess is explained with the $\gamma$-ray emission from the p$\gamma$ process
via the $\pi^{0}$ decay. The cascade emission of the electrons produced via
the $\gamma\gamma$ and p$\gamma$ interactions results in a keV-MeV excess in
the SEDs, illustrated as a bump or plateau. This extra photon field enhances
the production of TeV photons from the $p\gamma$ process, resulting in a
reduction of the proton power by about one order of magnitude. However, the
derived powers are still 3–4 orders of magnitude larger than the Eddington
limit, being challenged by the current black hole accretion physics. Applying
our model to Mrk 421, we propose that synergic observations with current and
upcoming TeV and keV-MeV telescopes for its tentative TeV and MeV excesses can
give insights to the hadronic process in its jet.

active galactic nuclei: individuals (1ES 0229+200, 1ES 0347-121, 1ES 1101-232,
H2356-309, Mrk 421) - Relativistic Jets; Non-thermal Radiation; Hadronic
process

##  1 Introduction

Most TeV $\gamma$-ray emitting active galactic nuclei (AGNs) are BL Lacertaes
(BL Lacs)111TeVCat official website: http://tevcat.uchicago.edu/. Their
observed broadband spectra energy distributions (SEDs) are characterized as
bimodal distribution with peaks at the IR-Optical-X-ray and the MeV-GeV bands
(Urry & Padovani, 1995). It is generally believed that the SEDs are attributed
to leptonic processes of relativistic electrons accelerated within the jets.
The peak at the low-energy band is explained as the synchrotron radiations of
the electrons, and the high-energy peak results from the Synchrotron-Self-
Compton (SSC) scattering process (e.g. Maraschi et al. 1992; Ghisellini &
Madau 1996). Zhang et al. (2012) presented a systematical analysis for the
SEDs of 24 TeV BL Lacs with the single-zone leptonic model. They showed that
most SEDs are indeed well represented with the model.

TeV-selected BL Lacs are the most valuable sources for studying particle
accelerations and radiation physics. The obstacle for this purpose is that the
extragalactic TeV photons are absorbed by the extragalactic background light
(EBL) via the pair production process when traveling through the universe
(Gould & Schréder, 1966). Interestingly, the intrinsic SEDs of some BL Lacs
show an apparent excess in the TeV band (the TeV excess) by correcting the EBL
absorption. The excess is an extra hard spectral component that is difficult
to be explained with the single-zone leptonic model. For instance, the High
Energy Stereoscopic System (HESS) observations show that the intrinsic photon
spectrum in an outburst of 1ES 1101-232 during March 2004 to June 2005 is very
hard with a power-law index $\sim 1.5$ in the energy range from 0.23 to 4 TeV
(Aharonian et al., 2007a). In the hadronic model, the TeV excess is suggested
to be produced through proton synchrotron radiation or the decay of neutral
pions produced via the $pp$ or $p\gamma$ interaction process (Mannheim, 1993).
Cao & Wang (2014) studied the possible hadronic origin of the hard
$\gamma$-ray spectrum and suggested that the TeV excess would originate from
the decay of neutral pion produced through p$\gamma$ interactions, in which
the soft photons are from the electron synchrotron radiation within the jet.
Sahu et al. (2013) proposed that the observed TeV orphan flare of 1ES 1959+650
is also from the photohadronic process. Mastichiadis et al. (2013)
investigated the X-ray and $\gamma$-ray variabilities of Mrk 421 in various
leptohadronic scenarios and found that the hadronic model can reproduce the
quadratic behavior between X-ray and TeV observations. Alternatively, it has
been found that the TeV excess could originate from leptonic processes of
different radiation zones. Such as inverse Compton scattering of cosmic
microwave background (IC/CMB) photons by electrons in the large scale jet
(Böttcher et al., 2008; Yan et al., 2012), or SSC process in a more distant
region away from the emitting region that is responsible for the outbursts in
the keV-GeV band (Zhang, 2009).

Hadronic models for explaining the gamma-ray emission of AGN were also
proposed. These models usually require a jet power being larger than the
Eddington luminosity by 1-2 orders of magnitude (e.g. Böttcher et al. 2013).
Although a moderate proton power could be enough for modeling the MeV-GeV bump
of the SEDs (e.g. Cerruti et al. 2015), an extremely large jet power, which is
larger than the Eddington luminosity by 4-5 orders of magnitude, is required
for attributing the TeV excess to the contributions of the $p\gamma$ process
(Cao & Wang 2014). This sharpens the issue of the breaking Eddington limit.
Despite of the Eddington luminosity is not a strict limit, it is a reasonable
approximation for the maximum jet power of a blazar (e.g. Zdziarski & Bottcher
2015).

The observed TeV excess indicates that the $e^{\pm}$-pair production and its
cascade emission within the jet may play a key role in shaping the observed
SED in the keV-GeV band (Cerruti et al., 2015; Xue et al., 2021). More
importantly, the cascade emission should offer an extra photon field for
enhancing the very high energy (VHE) gamma-ray emission via the p$\gamma$
process, leading to a reduction of the proton power that may relieve the
tension of exceeding the Eddington luminosity to some extent. This letter aims
to explore the hadronic origin of the TeV excesses in the framework of the
single-zone leptohadronic model by considering radiations of the relativistic
electrons and protons together with their cascade process in detail. We
present the selected broadband SEDs in Sec.2, and describe our model in Sec.3.
The results are reported in Sec.4. Conclusion and discussion are given in
Sec.6. Throughout the paper, cosmological parameters of $H_{\rm
0}=70{\rm\,km\,s^{-1}\,Mpc^{-1}}$, $\Omega_{\rm m}=0.3$, and
$\Omega_{\rm\Lambda}=0.7$ are adopted.

##  2 Intrinsic TeV excess in Selected-SEDs of BL Lacs

We select some BL Lacs with clear detections of a TeV excess component over
the leptonic model prediction in their SEDs. They are taken from Zhang et al.
(2012), in which a large SED sample of 24 TeV BL Lacs is compiled and
represented with the single-zone leptonic model. We derive the intrinsic SEDs
by correcting the EBL absorption with the model in Finke et al. (2010)222The
EBL model by Franceschini et al. (2008) is also adopted for double check the
EBL corrections in our analysis. The difference of the derived intrinsic
spectra induced by the two EBL models is not significant.. Based on the
leptonic model results presented in Zhang et al. (2012), we finally select the
following 4 BL Lacs that have a TeV excess component over the leptonic model
prediction with a confidence level of $3\sigma$ at the energy band of $>$0.1
TeV. The $\sigma$ value is estimated as $\sigma=\sum_{\rm i}(|F_{\rm
int,i}-F_{\rm t,i}|/\sigma_{\rm i})$, where $F_{\rm int,i}$ and $\sigma_{\rm
i}$ are the intrinsic flux and its uncertainty, respectively, and $F_{\rm
t,i}$ is the intrinsic flux predicted by the leptonic model. The selected SEDs
for the 4 BL Lacs are shown in Fig.1 and described as follows.

  * •

1ES 0229+200 ($z=0.14$; Aharonian et al. 2000). This BL Lac is known for its
hard spectrum in the TeV band. HESS observations of this source in 2005 and
2006 show a photon spectral index of $\sim 2.5$ in the energy range 500 GeV to
$\sim 15$ TeV, yielding a very hard intrinsic spectrum at very high energy
(VHE) band (Aharonian et al., 2007b). Between 2009 October and 2013 January,
the VERITAS detects an excess of 489 $\gamma$-ray events from 1ES 0229+200 in
the energy range of $0.29-7.6$ TeV, with an average integral flux of $\sim
23.3\times 10^{-9}\,{\rm photons\,m^{-2}\,s^{-1}}$(Aliu et al., 2014). The
long-term observations of the MAGIC, with a total of 265 hr good-quality data
from 2010 to 2017, suggested that the intrinsic spectral index of 1ES 0229+200
is $\sim$ 1.8 at 521 GeV (Acciari et al., 2020), for which the EBL model of
Franceschini et al. (2008) is considered. The hadronic origin of the VERITAS
observation has been discussed with various models, for example, the proton-
synchrotron and the leptohadronic models (Cerruti et al., 2015). In our work,
multi-band data include VHE and X-ray observations, which are collected from
Aharonian et al. (2007b) and Zhang et al. (2012), respectively.

  * •

1ES 0347–121 ($z=0.188$; Woo et al. 2005). It was observed by the HESS between
2006 August and December in the VHE band. The detected photon spectral index
is $\sim 3.10$ in the energy range $0.25-3$ TeV, with an integral flux $\sim
2\%$ of that of the Crab Nebula (Aharonian et al., 2007c). The hadronic origin
of this VHE observation has also been studied and discussed (Cerruti et al.,
2015). The multi-band data used in our intrinsic broadband SED analysis are
obtained from Aharonian et al. (2007c), which includes the ATOM and SWIFT
observations.

  * •

1ES 1101–232 ($z=0.186$; Wolter et al. 2000). From March 2004 to June 2005,
HESS detected 1ES 1101–232 with an excess of 649 photons. This VHE observation
shows a very hard intrinsic spectrum with a photon index $\sim 1.5$ from 0.23
to 4 TeV (Aharonian et al., 2007a). Previous studies have indicated that the
VHE emission may originate from the hadronic processes (e.g. Cao & Wang 2014;
Cerruti et al. 2015). In this work, the multi-band data of the broadband SED
including the HESS observation, as well as the semi-simultaneous X-ray
observation of the XMM-Newton and RXTE are taken from Costamante (2007).

  * •

H2356–309 ($z=0.165$; Bersanelli et al. 1992). The VHE observations of
H2356–309 have been reported twice by the HESS. The first one was from 2004
June to December with a photon index $\sim 3.09$ in the energy range from 0.2
to 1.3 TeV (Aharonian et al., 2007d). The second one was observed over a long-
term period, from 2004 to 2007, yielding the integral flux $\sim 3.06\times
10^{-12}\,{\rm photons\,cm^{-2}\,s^{-1}}$ above 240 GeV and the photon index
$\sim 3.06$ in the energy range from 200 GeV to 2 TeV (H. E. S. S.
Collaboration et al., 2010). It has been found that the H2356–309 location of
the former is coincident with the error circles of three IceCube events, and
the VHE emission may be from the hadronic processes (Sahu & Miranda, 2015).
The multi-band data used in this paper are obtained from H. E. S. S.
Collaboration et al. (2010), incorporating the simultaneous optical/UV and
X-ray observations from the XMM-Newton, NRT radio observation, and ATOM
optical observation.

##  3 Model

We employ a single-zone leptohadronic model to explain the broadband SEDs of
the selected BL Lacs. The bimodal SED from the radio to the MeV-GeV
$\gamma$-ray band are attributed to the synchrotron radiations of electrons
accelerated in the jet. The TeV excess is explained as the production of the
interaction between accelerated protons and photons (p$\gamma$ process).
Cascade emission is considered following that presented by Böttcher et al.
(2013). The injected cascade electrons include electrons produced via the
internal $\gamma\gamma$ absorption for photon fields of both the leptonic and
hadronic processes, and electrons from charged pion decay in the $p\gamma$
process (e.g. Atoyan & Dermer 2003; Romero & Vila 2008), as well as electrons
resulting from the Bethe-Heitler process (e.g. Chodorowski et al. 1992; Kelner
& Aharonian 2008). We describe the model in the following.

###  3.1 The Leptonic Radiations

Our leptonic model is adopted as Zhang et al. (2012). The emitting region is
assumed to be a single-zone sphere with a radius $R_{\rm b}$ moving with a
bulk Lorentz factor $\Gamma$ (or a velocity in unit of the light speed,
$\beta$). $R_{\rm b}$ can be estimated with the minimum variability timescale,
$R_{\rm b}=c\delta_{\rm D}\Delta t/(1+z)$, where
$\delta_{D}=1/\Gamma(1-\beta\cos\theta)$ is the Doppler factor with an viewing
angle of $\theta$ to the jet axis. For BL Lacs, the jet orientation is close
to the line of sight, and thus $\delta_{\rm D}\sim\Gamma$. The primary
electron spectrum is taken as a broken power-law function with indices $p_{1}$
and $p_{2}$ breaking at $\gamma_{\rm e,b}$ in the range of [$\gamma_{\rm
e,min}$, $\gamma_{\rm e,max}$].

The bimodal SED feature in the radio-optical-X-ray-GeV $\gamma$-ray band is
attributed to the synchrotron radiation and SSC process of the relativistic
electrons. Our calculation of the synchrotron process is based on an
approximation method (Aharonian et al., 2010), for which the synchrotron self-
absorption (Katarzyński et al., 2001) is considered. The SSC process is
calculated with analytical approximations (Khangulyan et al., 2014). The
internal $\gamma\gamma$ absorption (Finke et al., 2008) is yielded, where the
high-energy and low-energy photons are from the SSC and synchrotron
radiations, respectively. The radiative cooling of electrons can be estimated
by the energy density of the magnetic field and the synchrotron photon field,
in which the Klein-Nishina effect is incorporated (Moderski et al., 2005).

###  3.2 The Hadronic Radiation

We consider only the radiative cooling of the protons via the $p\gamma$
process. There are two main channels of this process: (1) The Bethe–Heitler
(B–H) process, ${\rm p}+\gamma\rightarrow{\rm p}+{\rm e^{-}}+{\rm e^{+}}$ for
photons with a threshold energy of $\sim 1$ MeV in the proton rest frame
(Kelner & Aharonian, 2008); (2) The photomeson process of single-pion pass
${\rm p}+\gamma\rightarrow{\rm p}+{\rm a}\pi^{0}+{\rm b}(\pi^{+}+\pi^{-})$ and
multi-pion pass ${\rm p}+\gamma\rightarrow{\rm n}+\pi^{+}+{\rm a}\pi^{0}+{\rm
b}(\pi^{+}+\pi^{-})$ (Romero & Vila, 2008) for photons with a threshold energy
of $\sim 145$ MeV in the proton rest frame (Atoyan & Dermer, 2003). A
$\pi^{0}$ particle decays to two energetic photons, which may result in the
observed TeV excess. The electrons generated through the two channels may be
cooled via the synchrotron radiation and SSC process. This electromagnetic
cascade emission may modify the SED in the radio-optical-X-ray-GeV band.

Our calculations of the $p\gamma$ process are based on the semi-analytical
method under the framework of Atoyan & Dermer (2003) and Romero & Vila (2008).
The collision rate of a proton with the Lorentz factor $\gamma_{\rm p}$
through the $p\gamma$ process with a soft photon field can be expressed as

$\omega_{\rm p\gamma,i}(\gamma_{\rm p})=\frac{c}{2\gamma_{\rm
p}^{2}}\int_{\frac{\epsilon_{\rm
th,i}^{{}^{\prime}}}{2\gamma_{p}}}^{\infty}d\epsilon\frac{n_{\rm
ph}(\epsilon)}{\epsilon^{2}}\int_{\epsilon_{\rm
th,i}^{{}^{\prime}}}^{2\epsilon\gamma_{\rm
p}}d\epsilon^{{}^{\prime}}\sigma_{\rm
p\gamma,i}(\epsilon^{{}^{\prime}})\epsilon^{{}^{\prime}},$ (1)

where the subscript ${\rm i}={\rm e},\pi$ represents for the B–H and
photomeson processes respectively, $n_{\rm ph}(\epsilon)$ is the number
density of photons with energy $\epsilon$ (the synchrotron radiation of
electrons in our case), $\epsilon^{{}^{\prime}}$ and $\epsilon_{\rm
th,i}^{{}^{\prime}}$ are the photon energy and photon threshold energy in the
proton rest frame, and $\sigma_{\rm p\gamma,i}$ is the cross section. Thus,
the cooling rate can be given by

$\tau_{\rm p\gamma,i}^{-1}(\gamma_{\rm p})=\frac{c}{2\gamma_{\rm
p}^{2}}\int_{\frac{\epsilon_{\rm
th,i}^{{}^{\prime}}}{2\gamma_{p}}}^{\infty}d\epsilon\frac{n_{\rm
ph}(\epsilon)}{\epsilon^{2}}\int_{\epsilon_{\rm
th,i}^{{}^{\prime}}}^{2\epsilon\gamma_{\rm
p}}d\epsilon^{{}^{\prime}}\sigma_{\rm p\gamma,i}(\epsilon^{{}^{\prime}})K_{\rm
p\gamma,i}(\epsilon^{{}^{\prime}})\epsilon^{{}^{\prime}},$ (2)

where $K_{\rm p\gamma,e}$ is the inelasticity factor. The cross section
$\sigma_{\rm p\gamma,e}$ and inelasticity $K_{\rm p\gamma,e}$ can be found in
Chodorowski et al. (1992). For the photomeson process, the cross section and
inelasticity can be approximated as

$\sigma_{\rm
p\gamma,\pi}(\epsilon^{{}^{\prime}})\approx\left\\{\begin{array}[]{lll}3.4\times
10^{-28}\,{\rm cm^{-2}},&&200\,{\rm
MeV}\leqslant\epsilon^{{}^{\prime}}\leqslant 500\,{\rm MeV}\\\ 1.2\times
10^{-28}\,{\rm cm^{-2}},&&\epsilon^{{}^{\prime}}>500\,{\rm
MeV},\end{array}\right.$ (3) $K_{\rm
p\gamma,\pi}(\epsilon^{{}^{\prime}})\approx\left\\{\begin{array}[]{lll}0.2,&&200\,{\rm
MeV}\leqslant\epsilon^{{}^{\prime}}\leqslant 500\,{\rm MeV}\\\
0.6,&&\epsilon^{{}^{\prime}}>500\,{\rm MeV},\end{array}\right.$ (4)

in which the two energy ranges are correspond to single-pion and multi-pion
channels, respectively.

The proton distribution in the energy range [$E_{\rm p,min}$, $E_{\rm p,max}$]
is assumed as an exponential cutoff power-law with an index $\alpha$ and an
cutoff energy $E_{\rm p,cut}$, i.e.,

$N_{\rm p}(E_{\rm p})=N_{0}\left(\frac{E_{\rm p}}{1\,{\rm
eV}}\right)^{-\alpha}{\rm exp}\left(-\frac{E_{\rm p}}{E_{\rm
p,cut}}\right),\qquad E_{\rm p,min}\leqslant E_{\rm p}\leqslant E_{\rm
p,max}.$ (5)

The electron energy produced in the B–H process depends on the proton energy
as $E_{\rm e}=\frac{1}{2}\overline{K}_{\rm p\gamma,e}E_{\rm p}$, and its
productivity is

$Q_{\rm BH,e}(E_{\rm e})=\frac{4}{\overline{K}_{\rm p\gamma,e}}N_{\rm
p}(\frac{2}{\overline{K}_{\rm p\gamma,e}}E_{\rm e})\omega_{\rm
p\gamma,e}(\frac{2}{\overline{K}_{\rm p\gamma,e}}\frac{E_{\rm e}}{m_{\rm
p}c^{2}}).$ (6)

The mean inelasticity $\overline{K}_{\rm p\gamma,e}$ can be approximated to
$2m_{\rm e}/m_{\rm p}$ (Romero & Vila, 2008). Due to that electrons produced
through a single proton-photon interaction in the B–H process has a broad
energy distribution, the productivity calculated through the mean inelasticity
and Eq.(6) is not accurate enough. Hence, we adopt an analytic method
introduced by Kelner & Aharonian (2008) to calculate the electron production.
This method is based on the differential cross section of the process and can
provide better accuracy.

Mesons ($\pi^{0}$ and $\pi^{\pm}$) produced in the p$\gamma$ interaction would
inherit about 20% of the energy from the collided proton and further pass it
to photons and electrons, $E_{\gamma}\approx 0.1E_{\rm p}$ and $E_{\rm
e}\approx 0.05E_{\rm p}$. The photon emissivity from $\pi^{0}$ decay and the
electron emissivity from $\pi^{\pm}$ decay can be calculated as

$Q_{p\gamma,\gamma}(E_{\gamma})=20N_{\rm p}(10E_{\gamma})\omega_{\rm
p\gamma,\pi}(10E_{\gamma})n_{\pi^{0}}(10E_{\gamma}),$ (7) $Q_{\rm
p\gamma,e}(E_{\rm e})=20N_{\rm p}(20E_{\rm e})\omega_{\rm
p\gamma,\pi}(20E_{\rm e})n_{\pi^{\pm}}(20E_{\rm e}),$ (8)

where $n_{\pi^{0}}\approx 0.5p_{1}+p_{2}$ and $n_{\pi^{\pm}}\approx
0.5p_{1}+2p_{2}$ are defined as the mean number of $\pi^{0}$ and $\pi^{\pm}$
created per collision, respectively. $p_{1}$ is the possibility of the
p$\gamma$ interaction through the single-pion pass and $p_{2}=1-p_{1}$ is for
the multi-pion pass. They can be derived by defining the mean inelasticity as

$\overline{K}_{\rm p\gamma,\pi}=\tau_{\rm p\gamma,\pi}^{-1}\omega_{\rm
p\gamma,\pi}^{-1}=0.2p_{1}+0.6(1-p_{1}).$ (9)

The specific luminosity of the gamma-ray photons from the $\pi^{0}$ decay in
the comoving frame is given by

$L_{\rm E}(E_{\gamma})=4\pi E_{\gamma}Q_{\rm
p\gamma,\gamma}(E_{\gamma})V\times\frac{1-e^{-\tau_{\gamma\gamma}}}{\tau_{\gamma\gamma}},$
(10)

where $V=4\pi R_{\rm b}^{3}/3$ is the volume of the radiating region and
$\tau_{\gamma\gamma}$ is the optical depth of the internal $\gamma\gamma$
interaction (Finke et al., 2008). Without considering the EBL absorption, the
intrinsic specific flux in the observed frame is $F_{\rm
E}(E_{\rm\gamma,obs})=(1+z)\delta_{\rm D}^{3}L_{\rm E}(E_{\gamma})/4\pi D_{\rm
L}^{2}$, where $E_{\rm\gamma,obs}=E_{\gamma}\delta_{\rm D}/(1+z)$, $D_{\rm L}$
is the luminosity distance.

###  3.3 The Cascade Process

Generally, the pair cascade process is initiated by the internal
$\gamma\gamma$ interaction of existing photon fields inside the radiation
region, creating a loop between electron generation, cooling, and escaping.
The system eventually reaches a temporary equilibrium and returns a stationary
secondary electron distribution. The cascade process is calculated in a time-
independent method (Böttcher et al., 2013), and the Fokker–Planck equation is
written as

$\frac{\partial}{\partial\gamma_{\rm e}}(\dot{\gamma}_{\rm e}N_{\rm e}^{\rm
cas})=Q_{\rm e}(\gamma_{\rm e})+\dot{N}_{\rm e}^{\gamma\gamma}(\gamma_{\rm
e})+\dot{N}_{\rm e}^{\rm esc}(\gamma_{\rm e}),$ (11)

where $Q_{\rm e}$ is the sum of the B–H and photomeson electrons in the
p$\gamma$ interaction, $\dot{N}_{\rm e}^{\gamma\gamma}$ is the electron
injection from the internal $\gamma\gamma$ interaction, $\dot{N}_{\rm e}^{\rm
esc}$ is the escape term, and $N_{\rm e}^{\rm cas}$ is the final electron
distribution. The particle escape can be estimated as $\dot{N}_{\rm e}^{\rm
esc}(\gamma_{\rm e})=N_{\rm e}^{\rm cas}(\gamma_{\rm e})/t_{\rm esc}$, where
$t_{\rm esc}=R_{\rm b}/c$. $\dot{N}_{\rm e}^{\gamma\gamma}(\gamma_{\rm e})$ is
given by

$\dot{N}_{\rm e}^{\gamma\gamma}(\gamma_{\rm e})=f_{\rm
abs}(\epsilon_{1})(\dot{N}_{\epsilon_{1}}^{\rm 0}+\dot{N}_{\epsilon_{1}}^{\rm
syn}+\dot{N}_{\epsilon_{1}}^{\rm SSC})+f_{\rm
abs}(\epsilon_{2})(\dot{N}_{\epsilon_{2}}^{\rm 0}+\dot{N}_{\epsilon_{2}}^{\rm
syn}+\dot{N}_{\epsilon_{2}}^{\rm SSC}),$ (12)

where $\dot{N}^{\rm 0}$ is the initial $\gamma$-rays from the $\pi^{0}$ decay,
$\dot{N}^{\rm syn}$ and $\dot{N}^{\rm SSC}$ are the synchrotron and SSC
radiations of $N_{\rm e}^{\rm cas}$, $\epsilon_{1}$ and $\epsilon_{2}$ are the
energies of photons to produce $\gamma_{\rm e}$ electrons via the pair
production, and $f_{\rm abs}$ is the internal absorption coefficient.
$\epsilon_{1}$ and $\epsilon_{2}$ are estimated as $\epsilon_{1}=\gamma_{\rm
e}/f_{\gamma}$ and $\epsilon_{2}=\gamma_{\rm e}/(1-f_{\gamma})$, where
$f_{\gamma}$ is a energy fraction of the two photons, which is taken as 0.9
following Böttcher et al. (2013).

##  4 Results

Assuming that the protons are accelerated by the relativistic shocks via the
Fermi acceleration mechanism, the energy spectral index of proton distribution
is fixed at $p_{\rm p}=2.2$ in our calculations. In the leptohadronic model,
contributions from the following spectral components shape the SEDs, i.e., (1)
the primary leptonic emission components, including the synchrotron and SSC
emission components from the primary electron population, (2) the primary VHE
$\gamma$-ray component from the p$\gamma$ process via the $\pi^{0}$ decay, in
which photons are from the primary synchrotron radiations, (3) the cascade
leptonic component, including the synchrotron and SSC emission components from
the cascade electrons produced via the internal $\gamma\gamma$ interaction and
the $p\gamma$ interaction, and (4) the cascade VHE $\gamma$-ray component from
the cascade p$\gamma$ process, in which photons are from the cascade
synchrotron radiations. Our numerical results are shown in Fig.1 and the model
parameters are reported in Tab.1. One can observe that the model can well
represent the data. The parameters of the leptonic model part are generally
consistent with those reported in the literature (e.g. Zhang et al. 2012;
Tavecchio et al. 2010; Aharonian et al. 2007c; Costamante 2007; H. E. S. S.
Collaboration et al. 2010). The cut-off energy ($E_{\rm p,cut}$) of protons is
typically $10$ TeV.

As shown in Fig.1, the synchrotron radiations of the primary electron
population overwhelmingly dominate the SEDs in the optical-soft-X-ray band
($1\sim 10^{4}$ eV), and the SSC components of the primary electron population
peak at around 10 GeV for the four BL Lacs. The cascade synchrotron radiation
bump peaks at around 0.1 MeV. The peak fluxes of both the synchrotron
radiation and the SSC emission components of the cascade electron population
are lower than that of the primary electron population by a factor of $2\sim
10$. The combinations of the primary and cascade synchrotron emission
components shape the SEDs as a distinct hump or a plateau in the keV-MeV
energy band, as shown up in the SEDs of 1ES 0347–121, 1ES 1101–232, and
H2356–309. Such a keV-MeV excess would be a distinct feature for the hadronic
process.

The TeV $\gamma$-rays are contributed by the $\pi^{0}$-decay process in the
primary (or cascade) p$\gamma$ interactions, in which the soft $\gamma$-ray
photon fields are the synchrotron radiations of the primary (or cascade)
electron population. Both the primary and cascade SSC emission decrease
rapidly at high energies because of the Klein–Nishina effect, making
negligible contributions to the TeV excess. The $\pi^{0}$-decay process
completely dominates the TeV $\gamma$-ray flux, although the internal
$e^{\pm}$-pair production is very significant for these VHE $\gamma$-ray
photons. Since we have fixed the energy spectral index of the injected proton
distribution, the shape of the hadronic spectrum depends largely on the soft
photon fields. Note that the photon threshold energy of the p$\gamma$
interaction is $\epsilon_{\rm th}>145\,{\rm MeV}/2\gamma_{\rm p}$. Thus, the
soft photons should be at least as energetic as $10$ keV for generating the
TeV $\gamma$-rays. For 1ES 0229+200 and 1ES 0347–121, the peak flux of the TeV
$\gamma$-rays from the primary p$\gamma$ interactions is larger than that from
the cascade p$\gamma$ interactions. Inversely, the TeV $\gamma$-ray peak flux
in 1ES 1101–232 and H2356–309 is mainly contributed by the cascade p$\gamma$
interactions. The intrinsic total peak flux of the TeV excess is $F_{\rm
int,p}\sim 1\times 10^{-11}\,{\rm erg\,cm^{-2}\,s^{-1}}$, comparable to the
synchrotron peak flux of the primary electrons.

##  5 Discussion

###  5.1 The Issue of Breaking Eddington Limit

The derived jet power breaks the Eddington luminosity limit is a long-standing
issue in single-zone hadronic (leptohadronic) models (e.g. Böttcher et al.
2013). This issue is sharpened while explaining the TeV excess with gamma-rays
generated through the p$\gamma$ process (Cao & Wang, 2014) due to the
inefficient radiative cooling of the protons. We compare the cooling rates of
protons and electrons in our model in Fig.2. It is found that the primary
electrons with energy $<10^{9}\sim 10^{10}$ eV are efficiently cooled via the
SSC process, and the cooling of electrons with energy $\gtrsim 10^{10}$ eV is
dominated by the synchrotron radiations due to that the Klein-Nishina effect
reduces the cross section of high energy electrons. The cooling rates of the
protons in the energy range from 1 TeV to 100 TeV via the B–H and photomeson
processes are about $10^{-15}\sim 10^{-12}\,{\rm s^{-1}}$. In most cases, the
photomeson process is less efficient than the B–H process, which has a much
lower threshold photon energy. However, the B–H process has a limited
influence in calculating the cascade process for relatively small cross-
sections and inelasticity. Electron injections from the internal pair
production and the charged pion decay are the major parts that decide the
stationary electron distribution of the cascade process.

It was also proposed that the TeV excess is resulted from the $p\gamma$
process with an extremely dense photon field (Sahu et al., 2013). However, how
to form such a photon field is a great challenge in this model. In our
analysis, the cascade synchrotron photon field is considered in the p$\gamma$
process to enhance the production of VHE photons. We estimate the proton
powers of the sources as $P_{\rm p}=\pi R_{\rm b}^{2}\Gamma^{2}cU_{\rm p}$
with the parameters reported in Tab.1. Our results are $P_{\rm p}=2.19\times
10^{50}$, $4.31\times 10^{51}$, $3.06\times 10^{51}$, and $9.85\times 10^{50}$
${\rm erg\,s^{-1}}$ for 1ES 0229+200, 1ES 0347–121, 1ES 1101–232, and
H2356–309, and their Eddington luminosities are $2.19\times 10^{47}$,
$5.63\times 10^{46}$, $1.26\times 10^{47}$, and $5.02\times 10^{46}$ ${\rm
erg\,s^{-1}}$, respectively333Since no estimated central black hole mass of
1ES 1101–232 is available, we use an average black hold mass $10^{9}M_{\odot}$
instead.. Note that the derived proton powers in our model are lower than that
without considering the cascade synchrotron photon field by one order of
magnitude, i.e., $3.06\times 10^{51}$ vs. $3.10\times 10^{52}$ ${\rm
erg\,s^{-1}}$ of 1ES 1101–232 for instance. However, they are still 3-4 orders
of magnitude larger than the Eddington luminosities.

###  5.2 Magnetic Field Strength and the Proton Acceleration

In our leptohadronic model, the magnetic field strength is about $0.1\sim 1$
G, comparable to that used in leptonic models for explaining the SEDs below
the TeV band, but 2-3 orders of magnitude smaller than that required in the
proton-synchrotron hadronic model for representing the MeV-GeV gamma-ray bump
(e.g. Böttcher et al. 2013; Cerruti et al. 2015). We examine whether the
proton can be sufficiently accelerated under the relatively small magnetic
field strength. Assuming $B=0.5$ G, the Larmor radius of a 100 TeV proton is
$\sim 6.67\times 10^{11}\,{\rm cm}$, which is still smaller than the typical
emitting region size by 5-6 orders of magnitude. Therefore, the magnetic field
strength used in our model should be sufficient for accelerating the proton to
relativistic in the radiating region.

The particle acceleration mechanism in BL Lacs is still uncertain. In our
analysis, we assume that the charged particles are accelerated by the Fermi
acceleration mechanism via the relativistic shocks, and the proton
distribution $N_{\rm p}(E_{\rm p})\propto E_{\rm p}^{-2.2}$ is adopted. Note
that particles accelerated through magnetic reconnection may result in a
harder particle distribution (e.g. Zhu et al. 2016 and reference therein). The
number density of the low energy end of the proton distribution can be
significantly reduced, leading to a decrease of the required proton power in
our model. We test this scenario with 1ES 1101–232 by taking the proton
distribution index as -1.5. The SED is still can be well represented, and the
derived $P_{p}$ is reduced by a factor of 4.5, changing from $3.06\times
10^{51}$ to $6.96\times 10^{50}$ ${\rm erg\,s^{-1}}$. Meanwhile, narrowing the
energy range, i.e., setting a larger $E_{\rm p,min}$, for the proton
distribution can also help to reduce the power. The $E_{\rm p,min}$ value is
fixed at 10 GeV in our calculations for the 4 BL Lacs. By setting $E_{\rm
p,min}=100\,{\rm GeV}$, the derived proton power changes from $3.06\times
10^{51}$ to $1.58\times 10^{51}$ ${\rm erg\,s^{-1}}$ for 1ES 1101–232.
Unfortunately, neither of the approaches we have tested can lower the jet
power to an acceptable level concerning the Eddington luminosity that the
supper-Eddington issue still exists.

###  5.3 Synergy of the TeV excess and keV-MeV excess as a probe for the
hadronic process

Our analysis shows that the flux of the cascade synchrotron emission peaks at
the keV-MeV energy band. Since it is usually lower than the primary
synchrotron flux in the keV band, it may be featured as a keV-MeV excess (a
bump or plateau) in the observed SEDs. The TeV excess and the corresponding
keV-MeV excess would be promising probes for the hadronic process. Nearby Mrk
421 ($z=0.031$; Punch et al. 1992) is the best candidate to verify this
speculation. So far, Mrk 421 has been detected in the TeV band by the
H.E.S.S., VERITAS, ARGO-YBJ, and MAGIC telescopes. It was detected during
2003-2004 with an outburst of a peak flux $\sim 135$ mCrab in the X-ray band
and $\sim 3$ Crab in the $\gamma$-ray band (Błażejowski et al., 2005). A
tentative TeV excess with large error bars is presented in its SED, as shown
in Fig.4. The possible hadronic explanation of the VHE observations has
already been discussed (e.g. Mastichiadis et al. 2013; Zech et al. 2017). Here
we represent the SED with our leptohadronic model, and the result is also
shown in Fig.4. Same as the other four BL Lacs, the TeV excess of Mrk 421 is
mainly contributed by $\pi^{0}$ decay. Interestingly, Nandikotkur et al.
(2008) reported an intriguing convex break at 235 MeV in the SED of Mrk 421
observed with CGRO/EGRET. It is possible that the keV-MeV excess moves in the
EGRET energy band during the EGRET observation campaign, shaping the convex
break.

A great opportunity for testing our speculation is available with the current
and upcoming VHE telescopes. The Large High Altitude Air Shower Observatory
(LHAASO, Bai et al. 2019) is sensitive in the 0.1 TeV to 1 PeV energy band,
and the upcoming Cherenkov Telescope Array (CTA, Actis et al. 2011) is
sensitive in the 20 GeV to 300 TeV energy band444CTA official website:
https://www.cta-observatory.org/. Their synergy observations provide an
unprecedented chance to investigate the particle acceleration and radiation
physics of the TeV-PeV $\gamma$-rays. We show the delectability of the TeV
excess in the SEDs for the four BL Lacs with the LHAASO and CTA in Fig.3,
where one-year sensitivities are used. One can find that the TeV excess is not
clearly shown up in the observed SEDs due to the strong EBL absorption effect,
which usually flattens the observed spectrum above the 0.1 TeV range. In
contrast, it should be confidently detectable with the CTA up to $\sim 10$ TeV
for the four BL Lacs. The TeV excess in 1ES 0229+200 may also be marginally
detectable with the LHAASO in the energy range of several TeV. Meanwhile, The
TeV excess of Mrk 421 predicted by our model is detectable with the CTA and
LHAASO up to $\sim 40$ TeV. Its spectrum in the $8\sim 17$ TeV energy band is
predicted as $F_{\nu}\propto\nu^{-1.8}$. We have proposed to verify this TeV
excess using the first operation year data of the LHAASO.

The keV-MeV excesses of the four BL Lacs are around $0.1\sim 10$ MeV. Although
they are marginally in the energy range of the EGRET ($0.2\sim 100$ MeV;
Hartman et al. 1999), it is not sensitive enough for the detection. However,
the keV-MeV excess of H2356–309 was tentatively observed in the X-ray band,
showing up like an X-ray plateau with a distinct feature of spectral
hardening. Missions have been proposed for observations on the MeV
$\gamma$-rays that have sensitivities higher than current/past missions by one
or two orders of magnitude. For instance, the e-ASTROGAM is designed to
improve the instrument sensitivity in the energy range of $0.3\sim 100$ MeV
with one-year sensitivity of $\sim 1.3\times 10^{-12}\,{\rm
erg\,cm^{-2}\,s^{-1}}$ at $\sim 0.5$ MeV (de Angelis et al., 2018). As shown
in Fig.3, the keV-MeV excess should be detectable with the e-ASTROGAM.

##  6 Summary

We have investigated the hadronic origin of the observed TeV excess for four
selected TeV BL Lacs (1ES 0229+200, 1ES 0347–121, 1ES 1101–232, and H2356–309)
with a single-zone leptohadronic model by considering the cascade emission
within their jets in detail. We summarize our results as follows.

  * •

Their broadband SEDs are well represented by our model. The SED in the radio-
optical-X-ray-GeV gamma-ray energy band is mainly attributed to the
synchrotron radiations and the SSC process of the primary electron population.
The model parameters are consistent with those in leptonic models. The TeV
excess is explained with the VHE $\gamma$-ray emission from the p$\gamma$
process via the $\pi^{0}$ decay, assuming that the power-law index of the
proton distribution is -2.2. The target photon fields are from the synchrotron
radiations of the primary electron population and the cascade electron
population produced via the internal $\gamma\gamma$ absorption and the
p$\gamma$ process.

  * •

The cascade synchrotron radiations result in an excess in the keV-MeV band of
the SEDs, illustrated as a distinct bump or plateau as shown in the SEDs of
1ES 0347–121, 1ES 1101–232, and H2356–309. The keV-MeV excess enhances the
production of VHE photons in the $p\gamma$ process and reduces the proton
power by about one order of magnitude. However, the derived powers are still
3-4 orders of magnitude larger than the Eddington luminosity. Further tests by
setting the power-law index of the proton energy distribution as -1.5 and
enlarging the low energy boundary of the proton energy distribution, the
derived proton power is reduced by a factor around 2-4. The breaking Eddington
limit issue is still cannot be overcome. Thus, our model is challenged by the
current accretion picture of AGNs.

  * •

Tentative TeV excess and MeV excess were observed in nearby bright TeV source
Mrk 421. As predicted by our model, its TeV excess is mainly contributed by
the cascade p$\gamma$ emission, which dominates the observed SED beyond 5 TeV
and can be detectable with the CTA and LHAASO up to $\sim 40$ TeV. The convex
break at 235 MeV in the SED of Mrk 421 observed with the EGRET might be
resulted from its MeV excess. Mrk 421 is the best candidate for testing our
model by synergic observations of the TeV excess and keV-MeV excess with
LHAASO, CTA, and future missions that are sensitive in the keV-MeV energy
band.

We thank Jin Zhang, Xiang Yu Wang, Ruo Yu Liu, and Rui Xue for help
discussion. We also thank Jin Zhang for providing us the SED data. This work
is supported by the National Natural Science Foundation of China (Grant
No.12133003 and U1731239), Guangxi Science Foundation (grant No. 2017AD22006)
and Innovation Project of Guangxi Graduate Education (YCBZ2021025).

Table 1: Leptohadronic Model Parameters

Source | z | $\delta_{\rm D}$ | B | $\Delta t$ | $\gamma_{\rm e,min}$ | $\gamma_{\rm e,b}$ | $p_{1}$ | $p_{2}$ | $N_{0,e}$ | $E_{\rm p,max}$ | $E_{\rm p,cut}$ | $N_{\rm 0,p}$
---|---|---|---|---|---|---|---|---|---|---|---|---
|  |  | $(\rm G)$ | $\rm hr$ |  |  |  |  | $(\rm cm^{-3})$ | $(\rm eV)$ | $(\rm eV)$ | $(\rm cm^{-3}\,eV^{-1})$
1ES 0229+200 | $0.139^{(1)}$ | 8.4 | 0.48 | 24 | 2 | $3.9\times 10^{5}$ | 2.08 | 3.16 | $1.3\times 10^{3}$ | $1\times 10^{15}$ | $2.8\times 10^{13}$ | $1.5\times 10^{18}$
1ES 0347–121 | $0.188^{(2)}$ | 11 | 0.65 | 12 | 100 | $1.4\times 10^{5}$ | 2.42 | 3.5 | $8.2\times 10^{4}$ | $1\times 10^{14}$ | $8\times 10^{12}$ | $4.8\times 10^{19}$
1ES 1101–232 | $0.186^{(3)}$ | 12 | 1.05 | 12 | 5 | $6.5\times 10^{4}$ | 1.8 | 4.1 | 85 | $5\times 10^{14}$ | $2.4\times 10^{13}$ | $2.2\times 10^{19}$
H2356–309 | $0.165^{(4)}$ | 7.6 | 0.5 | 24 | 2 | $6.3\times 10^{4}$ | 2.1 | 3.4 | $3.8\times 10^{3}$ | $1\times 10^{15}$ | $2.8\times 10^{13}$ | $1.1\times 10^{19}$
Mrk 421 | $0.031^{(5)}$ | 29 | 0.058 | 3 | 20 | $4.5\times 10^{5}$ | 2.36 | 4.0 | $1.9\times 10^{5}$ | $1\times 10^{15}$ | $2.8\times 10^{13}$ | $4.5\times 10^{19}$

References.(1)Aharonian et al. (2000);(2)Woo et al. (2005);(3)Wolter et al.
(2000);(4)Bersanelli et al. (1992);(5)Punch et al. (1992)

Figure 1: The EBL-absorbtion-corrected SEDs (solid black dots with error bars)
and our leptohadronic model results (solid black curves). The gray points are
the observed data without the EBL absorption correction. The purple inverted
triangles indicate the upper limit of the Fermi/LAT observations from
Tavecchio et al. (2010). Color lines represent the different radiation
processes: blue for synchrotron, green for SSC, and red for $\pi^{0}$-decay.
Solid lines are for the emission of the primary electrons and dashed lines are
for the emission of the cascade electrons. The hadronic emission without
considering the internal $\gamma\gamma$ absorption are also drawn with red
translucent solid/dashed lines for demonstration.

Figure 2: The cooling rates of the primary synchrotron and SSC as well as the
primary p$\gamma$ components. The upper panel and the lower panel are the
cooling rates of the primary electrons and protons for the 4 BL Lacs.

Figure 3: Examination of the detection capability of some instruments for the
TeV and MeV excesses in the TeV and MeV bands. The solid lines are our model
results presented in Fig.1 by considering the EBL-absorption. The dashed lines
mark the EBL-absorbed SSC component of the primary electron population. The
one-year sensitivities of the instruments are plotted with dash-dotted lines
in different colors as marked in the plot.

Figure 4: Intrinsic (black dots) and observed (grey dots) SEDs of Mrk 421 and
our theoretical modeling (showing with the same symbols as that in Fig.1)
together with examination of detection capability with some instruments for
the TeV and MeV excesses in the TeV and MeV bands (showing with the same
symbols as that in Fig.3).

## References

  * Acciari et al. (2020) Acciari, V. A., Ansoldi, S., Antonelli, L. A., et al. 2020, ApJS, 247, 16, doi: 10.3847/1538-4365/ab5b98
  * Actis et al. (2011) Actis, M., Agnetta, G., Aharonian, F., et al. 2011, Experimental Astronomy, 32, 193, doi: 10.1007/s10686-011-9247-0
  * Aharonian et al. (2000) Aharonian, F. A., Akhperjanian, A. G., Barrio, J. A., et al. 2000, A&A, 353, 847
  * Aharonian et al. (2007a) Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2007a, A&A, 470, 475, doi: 10.1051/0004-6361:20077057
  * Aharonian et al. (2007b) Aharonian, F., Akhperjanian, A. G., Barres de Almeida, U., et al. 2007b, A&A, 475, L9, doi: 10.1051/0004-6361:20078462
  * Aharonian et al. (2007c) Aharonian, F., Akhperjanian, A. G., Barres de Almeida, U., et al. 2007c, A&A, 473, L25, doi: 10.1051/0004-6361:20078412
  * Aharonian et al. (2007d) Aharonian, F., Akhperjanian, A. G., Barres de Almeida, U., et al. 2007d, A&A, 475, L9, doi: 10.1051/0004-6361:20078462
  * Aharonian et al. (2010) Aharonian, F. A., Kelner, S. R., & Prosekin, A. Y. 2010, Phys. Rev. D, 82, 043002, doi: 10.1103/PhysRevD.82.043002
  * Aliu et al. (2014) Aliu, E., Archambault, S., Arlen, T., et al. 2014, ApJ, 782, 13, doi: 10.1088/0004-637X/782/1/13
  * Atoyan & Dermer (2003) Atoyan, A. M., & Dermer, C. D. 2003, ApJ, 586, 79, doi: 10.1086/346261
  * Bai et al. (2019) Bai, X., Bi, B. Y., Bi, X. J., et al. 2019, arXiv e-prints, arXiv:1905.02773.  https://arxiv.org/abs/1905.02773
  * Bersanelli et al. (1992) Bersanelli, M., Bouchet, P., Falomo, R., Tanzi, E. G. 1992, AJ, 104, 28. doi: 10.1086/116218
  * Błażejowski et al. (2005) Błażejowski, M., Blaylock, G., Bond, I. H., et al. 2005, ApJ, 630, 130, doi: 10.1086/431925
  * Böttcher et al. (2008) Böttcher, M., Dermer, C. D., & Finke, J. D. 2008, ApJ, 679, L9, doi: 10.1086/588780
  * Böttcher et al. (2013) Böttcher, M., Reimer, A., Sweeney, K., & Prakash, A. 2013, ApJ, 768, 54, doi: 10.1088/0004-637X/768/1/54
  * Cao & Wang (2014) Cao, G., & Wang, J. 2014, ApJ, 783, 108, doi: 10.1088/0004-637X/783/2/108
  * Cerruti et al. (2015) Cerruti, M., Zech, A., Boisson, C., & Inoue, S. 2015, MNRAS, 448, 910, doi: 10.1093/mnras/stu2691
  * Chodorowski et al. (1992) Chodorowski, M. J., Zdziarski, A. A., & Sikora, M. 1992, ApJ, 400, 181, doi: 10.1086/171984
  * Costamante (2007) Costamante, L. 2007, Ap&SS, 309, 487, doi: 10.1007/s10509-007-9418-7
  * de Angelis et al. (2018) de Angelis, A., Tatischeff, V., Grenier, I. A., et al. 2018, Journal of High Energy Astrophysics, 19, 1, doi: 10.1016/j.jheap.2018.07.001
  * Finke et al. (2008) Finke, J. D., Dermer, C. D., & Böttcher, M. 2008, ApJ, 686, 181, doi: 10.1086/590900
  * Finke et al. (2010) Finke, J. D., Razzaque, S., & Dermer, C. D. 2010, ApJ, 712, 238, doi: 10.1088/0004-637X/712/1/238
  * Franceschini et al. (2008) Franceschini, A., Rodighiero, G., & Vaccari, M. 2008, A&A, 487, 837, doi: 10.1051/0004-6361:200809691
  * Ghisellini & Madau (1996) Ghisellini, G., & Madau, P. 1996, MNRAS, 280, 67, doi: 10.1093/mnras/280.1.67
  * Gould & Schréder (1966) Gould, R. J., & Schréder, G. 1966, Phys. Rev. Lett., 16, 252, doi: 10.1103/PhysRevLett.16.252
  * H. E. S. S. Collaboration et al. (2010) H. E. S. S. Collaboration, Abramowski, A., Acero, F., et al. 2010, A&A, 516, A56, doi: 10.1051/0004-6361/201014321
  * Hartman et al. (1999) Hartman, R. C., Bertsch, D. L., Bloom, S. D., et al. 1999, ApJS, 123, 79, doi: 10.1086/313231
  * Katarzyński et al. (2001) Katarzyński, K., Sol, H., & Kus, A. 2001, A&A, 367, 809, doi: 10.1051/0004-6361:20000538
  * Kelner & Aharonian (2008) Kelner, S. R., & Aharonian, F. A. 2008, Phys. Rev. D, 78, 034013, doi: 10.1103/PhysRevD.78.034013
  * Khangulyan et al. (2014) Khangulyan, D., Aharonian, F. A., & Kelner, S. R. 2014, ApJ, 783, 100, doi: 10.1088/0004-637X/783/2/100
  * Mannheim (1993) Mannheim, K. 1993, Phys. Rev. D, 48, 2408, doi: 10.1103/PhysRevD.48.2408
  * Maraschi et al. (1992) Maraschi, L., Ghisellini, G., & Celotti, A. 1992, ApJ, 397, L5, doi: 10.1086/186531
  * Mastichiadis et al. (2013) Mastichiadis, A., Petropoulou, M., & Dimitrakoudis, S. 2013, MNRAS, 434, 2684, doi: 10.1093/mnras/stt1210
  * Moderski et al. (2005) Moderski, R., Sikora, M., Coppi, P. S., & Aharonian, F. 2005, MNRAS, 363, 954, doi: 10.1111/j.1365-2966.2005.09494.x
  * Nandikotkur et al. (2008) Nandikotkur, G., Jahoda, K. M., Georganopoulos, M., et al. 2008, preprint from NASA Technical Reports Server
  * Punch et al. (1992) Punch., M., Akerlof, C. W., Cawley, M. F., Chantell, M., Fegan, D. J., Fennell, S., Gaidos, J. A., et al. 1992, Nature, 358, 477, doi: 10.1038/358477a0
  * Romero & Vila (2008) Romero, G. E., & Vila, G. S. 2008, A&A, 485, 623, doi: 10.1051/0004-6361:200809563
  * Sahu & Miranda (2015) Sahu, S., & Miranda, L. S. 2015, European Physical Journal C, 75, 273, doi: 10.1140/epjc/s10052-015-3519-1
  * Sahu et al. (2013) Sahu, S., Oliveros, A. F. O., & Sanabria, J. C. 2013, Phys. Rev. D, 87, 103015, doi: 10.1103/PhysRevD.87.103015
  * Tavecchio et al. (2010) Tavecchio, F., Ghisellini, G., Ghirlanda, G., Foschini, L., & Maraschi, L. 2010, MNRAS, 401, 1570, doi: 10.1111/j.1365-2966.2009.15784.x
  * Urry & Padovani (1995) Urry, C. M., & Padovani, P. 1995, PASP, 107, 803, doi: 10.1086/133630
  * Wolter et al. (2000) Wolter, A., Tavecchio, F., Caccianiga, A., Ghisellini, G., Tagliaferri, G. 2000, A&A, 357, 429
  * Woo et al. (2005) Woo, J.-H., Urry, C. M., van der Marel, R. P., Lira, P., Maza, J. 2005, ApJ, 631, 762, doi: 10.1086/432681
  * Xue et al. (2021) Xue, R., Liu, R.-Y., Wang, Z.-R., Ding, N., & Wang, X.-Y. 2021, ApJ, 906, 51, doi: 10.3847/1538-4357/abc886
  * Yan et al. (2012) Yan, D., Zeng, H., & Zhang, L. 2012, MNRAS, 424, 2173, doi: 10.1111/j.1365-2966.2012.21376.x
  * Zdziarski & Bottcher (2015) Zdziarski, A. A., & Bottcher, M. 2015, MNRAS, 450, L21, doi: 10.1093/mnrasl/slv039
  * Zech et al. (2017) Zech, A., Cerruti, M., & Mazin, D. 2017, A&A, 602, A25, doi: 10.1051/0004-6361/201629997
  * Zhang (2009) Zhang, J. 2009, Research in Astronomy and Astrophysics, 9, 777, doi: 10.1088/1674-4527/9/7/006
  * Zhang et al. (2012) Zhang, J., Liang, E.-W., Zhang, S.-N., & Bai, J. M. 2012, ApJ, 752, 157, doi: 10.1088/0004-637X/752/2/157
  * Zhu et al. (2016) Zhu, Y.-K., Zhang J., Zhang H.-M., Liang E.-W., Yan D.-H., Cui W., Zhang S.-N. 2016, RAA, 16, 170, doi: 10.1088/1674-4527/16/11/170