# Persistent Coverage Control for Teams of Heterogeneous Agents*
– Extended Version –

Alberto Mellone1,2, Giovanni Franzini1,3, Lorenzo Pollini1 and Mario
Innocenti1 *This paper is an extended version of: “A. Mellone, G. Franzini, L.
Pollini and M. Innocenti, Persistent coverage control for teams of
heterogeneous agents, to appear in Proc. 57th IEEE Conference on Decision and
Control, Miami Beach, FL, USA, Dec. 2018”. *This work was partially supported
by the University of Pisa PRA-2017-15 and PRA-2017-41 research grants.
1University of Pisa, Department of Information Engineering, Largo L. Lazzarino
1, 56122 Pisa, Italy. 2A. Mellone is currently with Imperial College London,
Department of Electrical and Electronic Engineering, SW7 2AZ, London, UK. The
work presented herein was carried out when he was Graduate Student at
University of Pisa. 3G. Franzini is currently with United Technologies
Research Centre Ireland, 4th Floor, Penrose Business Center, Penrose Wharf,
T23 XN53, Cork, Republic of Ireland. The work presented herein was carried out
when he was Ph.D. Candidate at University of Pisa. E-mail addresses:
<EMAIL_ADDRESS>(A. Mellone<EMAIL_ADDRESS>(G. Franzini),
<EMAIL_ADDRESS>(L. Pollini<EMAIL_ADDRESS>(M.
Innocenti).

###### Abstract

A distributed cooperative control law for persistent coverage tasks is
proposed, capable of coordinating a team of heterogeneous agents in a
structured environment. Team heterogeneity is considered both at vehicles’
dynamics and at coverage capabilities levels. More specifically, the general
dynamics of nonholonomic vehicles are considered. Agent heterogeneous sensing
capabilities are addressed by means of the descriptor function framework, a
set of analytical tools for controlling agents involved in generic coverage
tasks. By means of formal arguments, we prove that the team performs the task
and no collision occurs between agents nor with obstacles. A numerical
simulation validates the proposed strategy.

##  I Introduction

Distributed robotic systems are currently one of the most extensively studied
topics worldwide. Besides the theoretical implications, the practical
advantages they yield make them particularly attractive. Indeed, a team of
autonomous agents, each one characterized by a reduced computational and
implementation complexity, can enforce local, basic actions to pursue a global
desired behavior in a cost-efficient way. Coverage problems are one of the
branches of multi-agent systems study, featuring autonomous vehicles exploring
the environment to gather data through on-board sensors or to provide
resources to different regions. Applications are extremely differentiated,
from surveillance to geographical surveys, lawn mowing, floor cleaning, and
environmental monitoring [1, 2, 3].

The coverage type considered in this paper is known as _persistent coverage_
and consists in agents recursively exploring the environment, thus modeling
scenarios in which information that has been collected in the past should be
updated over time. Alternately, resources previously delivered are supposed to
undergo a degradation process and, in turn, require the team to sweep again
already visited spots.

In the literature, different types of techniques have been used to address
this problem. Solutions based on waypoint guidance have been proposed for
agents with single-integrator [4, 5] and unicycle kinematics [6]. Path
planning techniques are proposed in [7, 8]. In the former reference, the
authors discuss the geometric design of paths aimed at performing the task,
that the agents must follow, whereas in [8] the single-integrator kinematics
of the agents is explicitly taken into account in the planning. A gradient
based law is instead proposed in [9] for a team of single-integrator agents.
Event-driven strategies for single-integrator agents are proposed in [10]. The
majority of the proposed solutions consider teams composed by agents with the
same kinematics, and with sensing capabilities characterized by the same
mathematical model, although the parameters may vary from agent to agent. The
control of a team of heterogeneous agents seems to be still un-addressed in
the literature, at the best of the authors’ knowledge.

The main contribution of the paper is a distributed cooperative control law
for persistent coverage tasks, designed for coordinating the motion of a team
composed by agents with general sensing capabilities and different dynamics,
extending the results presented in [11].

More specifically, we consider agents with general nonholonomic dynamics, thus
providing a control law that better complies with a rather wide class of
vehicles. Sensing heterogeneity is addressed by means of the analytical tools
offered by the _descriptor function framework_. Introduced in [12], the
framework addresses the need for a unified method to approach different kinds
of coverage tasks. It provides a mathematical abstraction, the _descriptor
function_ , which models both the task requirements and the agents’
contributions. This framework allows for an abstraction from the actual
coverage capabilities of each agent, as well as from the actual scenario. A
control law is then designed to make the agents cooperatively carry out the
assigned task by solely resorting to information retrieved through local
sensors or inter-agent communication. In addition, a unified approach to avoid
collisions between the agents and with the obstacles is provided, ensuring the
safe execution of the task. The validity and efficacy of the proposed
technique is assessed by means of formal arguments and is shown through a
simulation.

#### Notation

The field of reals is denoted with $\mathbb{R}$. The following subsets of
$\mathbb{R}$ are introduced: $\mathbb{R}_{0}^{+}=\left\\{x\in\mathbb{R}\colon
x\geq 0\right\\}$, and $\mathbb{R}^{+}=\mathbb{R}_{0}^{+}\backslash\\{0\\}$.
The $n\times n$ identity matrix is denoted with $\mathbf{I}_{n}$. Given a
positive semi-definite matrix $\mathbf{A}\in\mathbb{R}^{n\times n}$, the
weighted Euclidean norm of $\mathbf{v}\in\mathbb{R}^{n}$ is denoted with
$\|\mathbf{v}\|_{\mathbf{A}}=\sqrt{\mathbf{v}^{T}\mathbf{A}\mathbf{v}}$. The
support operator $\operatorname{supp}\left\\{\cdot\right\\}$ of a real- or
vector-valued function $\mathbf{f}(\mathbf{x})$ is defined as
$\operatorname{supp}\left\\{\mathbf{f}\right\\}=\left\\{\mathbf{x}\colon\mathbf{f}(\mathbf{x})\neq\mathbf{0}\right\\}$,
i.e. it is the set of points of its domain where $\mathbf{f}$ is not
identically zero. The operators $\operatorname{ver}\left\\{\cdot\right\\}$ and
$\operatorname{diag}\left\\{\cdot\right\\}$ represent the vertical and block
diagonal concatenations.

##  II Team Modeling

###  II-A Preliminaries

Consider a _team_ $\mathcal{T}$ composed by $N$ agents, allowed to operate in
a closed and bounded topological set $Q\subset\mathbb{R}^{n}$, with
$n\in\\{2,3\\}$. We denote with $\mathbf{p}_{i}\in\mathcal{C}(Q)$ the $i$-th
agent _pose_ , i.e. its position and orientation in $Q$, where
$\mathcal{C}(Q)$ represents the agent configuration space on $Q$.

The operational area may be populated by _obstacles_. Without loss of
generality, obstacles will be considered convex polygons ($n=2$) or polyhedra
($n=3$). The $k$-th obstacle occupies the region $Q_{\text{obs}}^{k}\subseteq
Q$. Given the obstacle set $\mathcal{O}$, we define
$Q_{\text{obs}}=\bigcup_{k\in\mathcal{O}}Q_{\text{obs}}^{k}$.

###  II-B Descriptor Function Framework Overview

A brief overview of the descriptor function framework is now provided,
reviewing the elements essential for the paper. The interested reader shall
refer to [13] for further details. The main idea behind the framework is the
use of a common abstraction, the _descriptor function_ , for modeling both the
agent capabilities and the deployment requested to the team in order to
accomplish the task.

The _agent descriptor function_ (ADF) $d_{i}\colon\mathcal{C}(Q)\times
Q\to\mathbb{R}_{0}^{+}$ represents the agent $i\in\mathcal{T}$ capability of
performing the assigned task or, equivalently for a coverage task, the amount
of sensing that it instantly provides. The ADF is assumed to be continuous and
differentiable over $Q$, its support being a connected set. For ADFs with
unbounded support it is reasonable to define:

$\operatorname{supp}_{\epsilon}\left\\{d_{i}(\mathbf{p}_{i},\cdot)\right\\}=\left\\{\mathbf{q}\in
Q\colon
d_{i}(\mathbf{p}_{i},\mathbf{q})>\epsilon,\,\epsilon\in\mathbb{R}_{0}^{+}\right\\}$
(1)

In the remainder, the distinction between spatially bounded and unbounded ADFs
will not be explicitly addressed, and, with a slight abuse of notation,
$\operatorname{supp}\left\\{\cdot\right\\}$ will also denote
$\operatorname{supp}_{\epsilon}\left\\{\cdot\right\\}$ in the latter case.

######  Remark 1.

Given the positive definiteness, continuity and differentiability of the ADFs,
the following result holds:

$\operatorname{supp}\left\\{\frac{\partial
d_{i}}{\partial\mathbf{p}_{i}}\right\\}\subseteq\operatorname{supp}\left\\{d_{i}\right\\}$
(2)

The sum of all the ADFs is called the _current task descriptor function_
(CTDF) and represents the cumulative amount of sensing that the team is
instantly achieving. It is denoted with $d\colon\mathcal{C}(Q)^{N}\times
Q\to\mathbb{R}_{0}^{+}$, and is defined as:

$d(\mathbf{p},\mathbf{q})=\sum_{i\in\mathcal{T}}d_{i}(\mathbf{p}_{i},\mathbf{q}),\quad\mathbf{p}=\operatorname{ver}\left\\{\mathbf{p}_{i}\right\\}\in\mathcal{C}(\mathbb{R}^{n})^{N}$
(3)

The _desired task descriptor function_ (TDF)
$d_{*}\colon\mathbb{R}_{0}^{+}\times Q\to\mathbb{R}_{0}^{+}$ describes how the
agents should be distributed in the operational area in order to maximize the
goal achievement. Therefore, $d_{*}(t,\mathbf{q})$ defines how much of the
available sensing capability is needed at time $t$ at point $\mathbf{q}$. With
the definition of the CTDF and of the TDF, it is natural to introduce an error
between the amount of sensing that the task requires and that is actually
provided by the team at each $\mathbf{q}\in Q$. This error is quantified by
the _task error function_ (TEF), defined as:

$e(t,\mathbf{p},\mathbf{q})=d_{*}(t,\mathbf{q})-d(\mathbf{p},\mathbf{q})$ (4)

The TEF models the excess or the lack of sensing over time at each point of
the environment.

###  II-C Agents Dynamics

We consider the following dynamics for the $i$-th agent:

$\displaystyle\mathbf{B}_{i}(\mathbf{x}_{i})\ddot{\mathbf{x}}_{i}+\mathbf{c}_{i}(\mathbf{x}_{i},\dot{\mathbf{x}}_{i})=\mathbf{D}_{i}(\mathbf{x}_{i})\boldsymbol{\tau}_{i}+\mathbf{A}_{i}(\mathbf{x}_{i})\boldsymbol{\lambda}_{i}$
(5a)
$\displaystyle\mathbf{A}_{i}(\mathbf{x}_{i})^{T}\dot{\mathbf{x}}_{i}=\mathbf{0}$
(5b)

where $\mathbf{x}_{i}\in\mathbb{R}^{q_{i}}$ is the state vector,
$\mathbf{B}_{i}\colon\mathbb{R}^{q_{i}}\to\mathbb{R}^{q_{i}\times q_{i}}$ is
the inertia matrix,
$\mathbf{c}_{i}\colon\mathbb{R}^{q_{i}}\times\mathbb{R}^{q_{i}}\to\mathbb{R}^{q_{i}}$
is the vector gathering the centripetal, Coriolis and potential terms,
$\mathbf{D}_{i}\colon\mathbb{R}^{q_{i}}\to\mathbb{R}^{q_{i}\times m_{i}}$
transforms the inputs $\boldsymbol{\tau}_{i}\in\mathbb{R}^{m_{i}}$ to
generalized forces,
$\mathbf{A}_{i}\colon\mathbb{R}^{q_{i}}\to\mathbb{R}^{q_{i}\times
q_{i}-m_{i}}$ is the matrix associated to the Pfaffian representation of the
$q_{i}-m_{i}$ kinematic constraints and
$\boldsymbol{\lambda}_{i}\in\mathbb{R}^{q_{i}-m_{i}}$ is the vector of
Lagrange multipliers [14].

The agent’s pose vector is extracted from the state using a continuous and
differentiable map
$\mathbf{h}_{i}:\mathbb{R}^{q_{i}}\rightarrow\mathcal{C}(\mathbb{R}^{n})$:

$\mathbf{p}_{i}=\mathbf{h}_{i}(\mathbf{x}_{i})$ (6)

In the rest of the paper we will consider the following equivalent
representation of the dynamics in (5):

$\dot{\mathbf{z}}_{i}=\begin{bmatrix}\dot{\mathbf{x}}_{i}\\\
\dot{\mathbf{v}}_{i}\end{bmatrix}=\begin{bmatrix}\mathbf{G}_{i}(\mathbf{x}_{i})\mathbf{v}_{i}\\\
\mathbf{0}\end{bmatrix}+\begin{bmatrix}\mathbf{0}\\\
\mathbf{I}_{m}\end{bmatrix}\mathbf{u}_{i}$ (7)

where $\mathbf{v}_{i}\in\mathbb{R}^{m_{i}}$ is known as _pseudo-velocity
vector_ , $\mathbf{u}_{i}\in\mathbb{R}^{m_{i}}$ is the new system input, known
as _pseudo-acceleration vector_ , and
$\mathbf{G}_{i}:\mathbb{R}^{q_{i}}\rightarrow\mathbb{R}^{q_{i}\times m_{i}}$
has columns that span the null space of $\mathbf{A}_{i}$, i.e.
$\mathbf{A}_{i}^{T}\mathbf{G}_{i}=\mathbf{0}$.

######  Proposition 1 ([14], Chap. 11.4).

Assume that $\mathbf{x}_{i}$ and $\mathbf{v}_{i}$ are measurable, and that
$\mathbf{G}_{i}(\mathbf{x}_{i})^{T}\mathbf{D}_{i}(\mathbf{x}_{i})$ is
invertible for all $\mathbf{x}_{i}$. Then the dynamics (7) is equivalent to
(5) under the nonlinear state-feedback law:

$\boldsymbol{\tau}_{i}=\left(\mathbf{G}_{i}(\mathbf{x}_{i})^{T}\mathbf{D}_{i}(\mathbf{x}_{i})\right)^{-1}\\!\left(\widetilde{\mathbf{B}}_{i}(\mathbf{x}_{i})\mathbf{u}_{i}+\widetilde{\mathbf{c}}_{i}(\mathbf{x}_{i},\mathbf{v}_{i})\right)$
(8)

where

$\displaystyle\\!\\!\\!\\!\\!\widetilde{\mathbf{c}}_{i}(\mathbf{x}_{i},\\!\mathbf{v}_{i})\\!$
$\displaystyle=\\!\mathbf{G}_{i}(\mathbf{x}_{i})^{T}\\!\\!\left(\\!\mathbf{c}_{i}(\mathbf{x}_{i},\\!\mathbf{G}_{i}(\mathbf{x}_{i})\mathbf{v}_{i})\\!+\\!\\!\mathbf{B}_{i}(\mathbf{x}_{i})\dot{\mathbf{G}}_{i}(\mathbf{x}_{i})\mathbf{v}_{i}\\!\right)$
(9) $\displaystyle\widetilde{\mathbf{B}}_{i}(\mathbf{x}_{i})$
$\displaystyle=\mathbf{G}_{i}(\mathbf{x}_{i})^{T}\mathbf{B}_{i}(\mathbf{x}_{i})\mathbf{G}_{i}(\mathbf{x}_{i})$
(10)

Therefore, if the vector $\mathbf{u}_{i}$ is obtained to control the system
with dynamics (7), relation (8) allows for the computation of the input
$\boldsymbol{\tau}_{i}$ for the system (5).

Note that the dynamics (7) along with (6) allow for the description of a
sufficiently wide variety of vehicle dynamic models, thus enhancing team
heterogeneity.

The following definitions will be used in the remainder of the paper:
$\mathbf{x}=\operatorname{ver}\left\\{\mathbf{x}_{i}\right\\}$,
$\mathbf{v}=\operatorname{ver}\left\\{\mathbf{v}_{i}\right\\}$,
$\mathbf{z}=[\mathbf{x}^{T},\mathbf{v}^{T}]^{T}$.

##  III Persistent Coverage Control

###  III-A Problem Statement

The _persistent coverage_ task deals with the recursive exploration of the
operational area: the information gathered by the agents through their sensors
becomes obsolete as time passes, thus requiring the team to visit regions of
$Q$ where the amount of actual information has faded.

To model this phenomenon, first we define the function
$I\colon\mathbb{R}_{0}^{+}\times Q\to\mathbb{R}_{0}^{+}$, that quantifies the
amount of useful information available at time $t$ in $\mathbf{q}\in Q$. The
following equation models the _information decay process_ :

$\dot{I}(t,\mathbf{q})=\delta
I(t,\mathbf{q})+d(\mathbf{p},\mathbf{q}),\quad\delta<0$ (11)

i.e. the agents through the CTDF give a positive contribution to information
gathering, while the information decay rate $\delta$ yields its degradation.
Denoting with $C^{*}>0$ the desired level of information that should be
maintained over $Q$, the TDF describing the persistent coverage task is:

$d_{*}(t,\mathbf{q})=\max\left\\{0,C^{*}-I(t,\mathbf{q})\right\\}$ (12)

while the TEF models the difference between regions insufficiently covered and
the current coverage provided by the team. Agents are expected to move towards
its minimization.

To quantify how effectively the task is being fulfilled, the error index
$\xi\colon\mathbb{R}_{0}^{+}\times\mathcal{C}(Q)^{N}\to\mathbb{R}_{0}^{+}$ is
introduced:

$\xi(t,\mathbf{p})=\int_{Q}f(e(t,\mathbf{p},\mathbf{q}))\sigma(\mathbf{q})\text{d}\mathbf{q}$
(13)

where $f\colon\mathbb{R}\to\mathbb{R}_{0}^{+}$ is a _penalty function_ defined
as $f(e)=\max\\{0,e\\}^{p}$ with $p=\\{2,3,...\\}$, and $\sigma\colon
Q\to\mathbb{R}_{0}^{+}$ is a weight that specifies the point importance in the
environment.

The task is properly accomplished if the error function is kept as low as
possible.

######  Remark 2.

The penalty function $f(\cdot)$ is continuous in $\mathbb{R}$ and strictly
convex in $\mathbb{R}^{+}$, along with $\partial^{n}f/\partial e^{n}$ for
$n<p$.

######  Remark 3.

Note that due to the definition of the TDF in (12), the error index
$\xi(\cdot)$ is always bounded:

$0\leq\xi(t,\mathbf{p})\leq\int_{Q}f(C^{*})\sigma\text{d}\mathbf{q}$ (14)

for $t\geq t_{0}$.

###  III-B Obstacles and Collision Avoidance

Besides covering, agents are expected to avoid collisions with others, as well
as with obstacles in the environment. Hence, the control action must guarantee
that the distance between the agent and either another agent or an obstacle
does not drop below a _safety threshold_ $r>0$. We make the following
assumption on the agents detection capabilities:

######  Assumption 1.

Each agent is able to detect other agents or obstacles when the relative
distance is under a _detection range_ $R>r$.

A team deployment is _safe_ if $\mathbf{p}\in\mathcal{P}(r)$, where:

$\mathcal{P}(x)=\left\\{\mathbf{p}\in\mathcal{C}(\mathbb{R}^{n})^{N}\colon\rho_{i,j}>x\;\land\;\rho_{i,k}^{o}>x\,\right.\\\
\left.\forall i,j\in\mathcal{T},\,i\neq j,\,\forall k\in\mathcal{O}\right\\}$
(15)

with $\rho_{i,j}$ and $\rho_{i,k}^{o}$ denoting the distances of agent $i$
from agent $j$ and from the $k$-th obstacle point closest to it, the latter
given by $\mathbf{c}_{k}\colon\mathcal{C}(Q)\to Q_{obs}$:

$\displaystyle\rho_{i,j}=\left\|\mathbf{S}(\mathbf{p}_{i}-\mathbf{p}_{j})\right\|,\quad\rho_{i,k}^{o}=\left\|\mathbf{S}\mathbf{p}_{i}-\mathbf{c}_{k}(\mathbf{p}_{i})\right\|$
(16)
$\displaystyle\mathbf{c}_{k}(\mathbf{p}_{i})=\operatorname*{arg\,min}_{\mathbf{q}\in
Q_{\text{obs}}}\|\mathbf{S}\mathbf{p}_{i}-\mathbf{q}\|$ (17)

The matrix $\mathbf{S}\in\mathbb{R}^{n\times\dim\\{\mathcal{C}(Q)\\}}$
extracts the position components from the agent’s pose vector.

To guarantee the team safety, the following _collision/obstacle avoidance
function_ is defined:

$\displaystyle v(\mathbf{p})$
$\displaystyle=v^{c}(\mathbf{p})+v^{o}(\mathbf{p})$ (18)
$\displaystyle=\sum_{i\in\mathcal{T}}v_{i}^{c}(\mathbf{p})+\sum_{i\in\mathcal{T}}v_{i}^{o}(\mathbf{p}_{i})$
(19)
$\displaystyle=\sum_{i\in\mathcal{T}}\sum_{j\in\mathcal{T}\backslash\\{i\\}}\\!\\!l(\rho_{i,j})+\sum_{i\in\mathcal{T}}\sum_{k\in\mathcal{O}}l(\rho_{i,k}^{o})$
(20)

where $l\colon\mathbb{R}\rightarrow\mathbb{R}_{0}^{+}$ is adapted from [15]:

$l(x)=\left(\min\left\\{0,\frac{x^{2}-R^{2}}{x^{2}-r^{2}}\right\\}\right)^{2}$
(21)

######  Remark 4.

The function $l(\cdot)$ is such that if $x\geq R$, then $l(x)=0$, while for
$x\in\left(r,R\right)$ the function is strictly decreasing, and
$\lim_{x\rightarrow r^{+}}l(x)=+\infty$. Given the avoidance function
definition in (20), collisions do not occur and
$\mathbf{p}(t)\in\mathcal{P}(r)$ is always verified, as long as $v(\cdot)$
attains finite values.

###  III-C Control Law Definition

The following control law is chosen for the agents:

$\mathbf{u}_{i}(t,\mathbf{x}_{i},\mathbf{p})=\mathbf{u}_{\xi,i}(t,\mathbf{x}_{i},\mathbf{p})+\mathbf{u}_{v,i}(\mathbf{x}_{i},\mathbf{p})-\mu\mathbf{v}_{i},\,\,\mu>0$
(22)

where:

$\displaystyle\mathbf{u}_{\xi,i}(t,\mathbf{x}_{i},\mathbf{p})$
$\displaystyle=-\beta\mathbf{G}_{i}(\mathbf{x}_{i})^{T}\frac{\partial\mathbf{h}_{i}(\mathbf{x}_{i})}{\partial\mathbf{x}_{i}}^{T}\\!\\!\frac{\partial\xi(t,\mathbf{p})}{\partial\mathbf{p}_{i}}^{T}\\!\\!\\!,\,$
$\displaystyle\beta>0$ (23)
$\displaystyle\mathbf{u}_{v,i}(\mathbf{x}_{i},\mathbf{p})$
$\displaystyle=-\gamma\mathbf{G}_{i}(\mathbf{x}_{i})^{T}\frac{\partial\mathbf{h}_{i}(\mathbf{x}_{i})}{\partial\mathbf{x}_{i}}^{T}\\!\\!\frac{\partial
v(\mathbf{p})}{\partial\mathbf{p}_{i}}^{T}\\!\\!\\!,\,$
$\displaystyle\gamma>0$ (24)

While $\mathbf{u}_{\xi,i}(\cdot)$ accelerates agent $i$ towards regions where
$I(\cdot)$ is lower, thus inducing a minimization of the error function,
$\mathbf{u}_{v,i}(\cdot)$ guarantees that it does not collide neither with
other agents nor with obstacles. Finally, the term $-\mu\mathbf{v}_{i}$ acts
as a damping term for the control law.

######  Theorem 1.

If the initial team deployment is safe, i.e.
$\mathbf{p}(t_{0})\in\mathcal{P}(r)$, then, under the control law (22), the
team performs the persistent coverage task and no collision occurs, i.e.
$\mathbf{p}(t)\in\mathcal{P}(r)$ for all $t\geq t_{0}$.

###### Proof.

Consider the following function:

$V(t,\mathbf{z})=\xi(t,\mathbf{p}(\mathbf{x}))+\frac{\gamma}{\beta}v(\mathbf{p}(\mathbf{x}))+\frac{1}{2\beta}\sum_{i\in\mathcal{T}}\left\|\mathbf{v}_{i}\right\|^{2}$
(25)

If $V(\cdot)$ attains finite values, and the agents act in order to decrease
it, then the persistent coverage task is safely executed. Note that if the
agents starts from a safe deployment, i.e.
$\mathbf{p}(t_{0})\in\mathcal{P}(r)$, then $V(\cdot)$ is finite at $t_{0}$.

Differentiation of (25) with respect to time yields:

$\dot{V}=\frac{\partial\xi}{\partial
t}+\sum_{i\in\mathcal{T}}\left[\left(\frac{\partial\xi}{\partial\mathbf{p}_{i}}+\frac{\gamma}{\beta}\frac{\partial
v}{\partial\mathbf{p}_{i}}\right)\frac{\partial\mathbf{h}_{i}}{\partial\mathbf{x}_{i}}\mathbf{G}_{i}+\frac{1}{\beta}\dot{\mathbf{v}}_{i}^{T}\right]\mathbf{v}_{i}$
(26)

Since $\dot{\mathbf{v}}_{i}=\mathbf{u}_{i}$ (see (7)), substitution of (22) in
(26) gives:

$\dot{V}=\frac{\partial\xi}{\partial
t}-\frac{\mu}{\beta}\sum_{i\in\mathcal{T}}\left\|\mathbf{v}_{i}\right\|^{2}$
(27)

The term $\partial\xi/\partial t$ can be further expanded obtaining:

$\displaystyle\dot{V}$ $\displaystyle=-\int_{Q}\frac{\partial f}{\partial
e}\,\dot{I}\,\sigma\,\text{d}\mathbf{q}-\frac{\mu}{\beta}\sum_{i\in\mathcal{T}}\left\|\mathbf{v}_{i}\right\|^{2}$
(28) $\displaystyle=-\int_{Q}\frac{\partial f}{\partial e}\,\delta
I\,\sigma\,\text{d}\mathbf{q}-\int_{Q}\frac{\partial f}{\partial
e}\,d\,\sigma\,\text{d}\mathbf{q}-\frac{\mu}{\beta}\sum_{i\in\mathcal{T}}\left\|\mathbf{v}_{i}\right\|^{2}$
(29) $\displaystyle=E_{0}+E_{1}+E_{2}$ (30)

Note that $E_{0}\geq 0$, since $\delta<0$, whereas $E_{1}\leq 0$ and
$E_{2}\leq 0$. The sign of $E_{0}$ is the consequence of the information
degradation that characterizes the persistent coverage task. This produces an
increment of $V(\cdot)$ with time. However, note that the only term in (25)
which is explicitly time-dependent is the error index $\xi(\cdot)$. Since
$\xi(\cdot)$ is bounded (see Remark 3), then the sign of $E_{0}$ cannot
produce an unlimited growth in $V(\cdot)$. The agents contribution is
expressed by the terms $E_{1}$ and $E_{2}$, which always give a negative
contribution. Hence, the persistent coverage task is properly executed by the
agents, and $V(\cdot)$ attains finite values, proving the safe execution of
the task.

∎

##  IV Control Law Decentralization

Apart from the terms
$\mathbf{G}_{i}(\mathbf{x}_{i})^{T}(\partial\mathbf{h}_{i}(\mathbf{x}_{i})/\partial\mathbf{x}_{i})^{T}$
and $-\mu\mathbf{v}_{i}$, which are fully known as long as the state
$\mathbf{z}_{i}$ is observable and the vehicle model is known, agent $i$ needs
to compute the gradients of the error $\xi(\cdot)$ and the avoidance functions
terms $v^{c}(\cdot)$ and $v^{o}(\cdot)$ in order to obtain
$\mathbf{u}_{i}(\cdot)$. It will be shown that such quantities can be obtained
in a distributed way without compromising the task fulfillment.

To this end, we introduce the proximity graph $\mathcal{G}(t)$, with the
agents as nodes. An edge between two agents exists if they are _neighbors_ ,
that is $\rho_{i,j}\leq R_{\text{com}}$, with $R_{\text{com}}$ denoting agents
_communication range_. The set of neighbors of agent $i$ at time $t$ will be
denoted with $\mathcal{N}_{i}(t)\subseteq\mathcal{T}$. In the following, we
will consider the next two assumptions holding:

######  Assumption 2.

The graph $\mathcal{G}(t)$ is always connected.

######  Assumption 3.

Let $r^{\text{cov}}_{i}=\max\\{\|\mathbf{q}_{i}-\mathbf{q}\|\colon
d_{i}(\mathbf{p}_{i},\mathbf{q})>0\\}$ and
$\overline{r}^{\text{cov}}=\max_{i\in\mathcal{T}}\\{r^{\text{cov}}_{i}\\}$. We
assume that $R_{\text{com}}>2\overline{r}^{cov}$. As a result, agents with
overlapping ADFs are neighbors.

###  IV-A Collision and Obstacle Avoidance

To enforce the collision and obstacle avoidance policy agent $i$ must compute
the gradient of $v(\cdot)$, that is:

$\frac{\partial
v(\mathbf{p})}{\partial\mathbf{p}_{i}}=2\\!\\!\\!\sum_{j\in\mathcal{T}/\\{i\\}}\frac{\partial
l(\rho_{i,j})}{\partial\rho_{i,j}}\frac{\partial\rho_{i,j}}{\partial\mathbf{p}_{i}}+\sum_{k\in\mathcal{O}}\frac{\partial
l(\rho_{i,k}^{o})}{\partial\rho_{i,k}^{o}}\frac{\partial\rho_{i,k}^{o}}{\partial\mathbf{p}_{i}}$
(31)

where:

$\frac{\partial l(x)}{\partial
x}=\left\\{\begin{array}[]{ll}0,&\,\,x>R\,\,|\,\,x<r\\\
4\frac{\left(R^{2}-r^{2}\right)\left(x^{2}-R^{2}\right)}{\left(x^{2}-r^{2}\right)^{3}}x,&\,\,R\geq
x>r\\\ \text{undefined},&\,\,x=r\end{array}\right.$ (32)
$\frac{\partial\rho_{i,j}}{\partial\mathbf{p}_{i}}=\frac{\left(\mathbf{p}_{i}-\mathbf{p}_{j}\right)^{T}\mathbf{S}^{T}\mathbf{S}}{\rho_{i,j}},\quad\frac{\partial\rho_{i,k}^{o}}{\partial\mathbf{p}_{i}}=\frac{\left(\mathbf{S}\mathbf{p}_{i}-\mathbf{c}_{k}(\mathbf{p}_{i})\right)^{T}\mathbf{S}}{\rho_{i,k}^{o}}$
(33)

Therefore, the collision and obstacle avoidance function gradients must be
computed only when another agent $j$ or an obstacle $k$ are within the agent
$i$ detection radius, that is $\rho_{i,j}\leq R$ and $\rho_{i,k}^{o}\leq R$,
respectively. Thus, the collision/obstacle avoidance is implicitly
distributed.

###  IV-B Coverage Control and Information Level Estimation

To accomplish the coverage task, agents need to compute the quantity
$\partial\xi/\partial\mathbf{p}_{i}$. Using the chain rule and observing that
the ADFs are dependent on the respective agent’s pose only, the following
holds:

$\frac{\partial\xi}{\partial\mathbf{p}_{i}}=-\int_{Q}\frac{\partial
f}{\partial e}\frac{\partial
d}{\partial\mathbf{p}_{i}}\sigma\text{d}\mathbf{q}=-\int_{Q_{i}}\frac{\partial
f}{\partial e}\frac{\partial
d_{i}}{\partial\mathbf{p}_{i}}\sigma\text{d}\mathbf{q}$ (34)

where $Q_{i}(\mathbf{p}_{i})=\operatorname{supp}\left\\{\partial
d_{i}(\mathbf{p}_{i},\cdot)/\partial\mathbf{p}_{i}\right\\}$. The weight
$\sigma$ is assumed known to each agent prior to the deployment. The gradient
$\partial d_{i}/\partial\mathbf{p}_{i}$ can be computed autonomously by each
agent. To obtain the term $\partial f/\partial e$ the knowledge of the TEF
$e(t,\mathbf{p},\mathbf{q})$ on $\mathbf{q}\in Q_{i}(\mathbf{p}_{i})$ is
needed. We recall that
$e(t,\mathbf{p},\mathbf{q})=d_{*}(t,\mathbf{q})-d(\mathbf{p},\mathbf{q})$, see
(4).

Because of Remark 1 and Assumption 3, if all agent $i$’s neighbors share with
it their poses and their ADF parameters, agent $i$ is able to compute the CTDF
$d(\mathbf{p},\mathbf{q})$ on $\mathbf{q}\in Q(\mathbf{p}_{i})$.

The computation of $d_{*}(t,\mathbf{q})$ requires the knowledge of the
attained level of information $I(t,\mathbf{q})$. We now introduce an algorithm
for the decentralized estimation of $I(t,\mathbf{q})$, adapted from [16].
Whereas it was originally formulated in a discrete time context, the algorithm
is here rearranged to fit in a continuous time scenario. The conventions of
the descriptor function framework will be employed. For the sake of
compactness, the quantities’ dependence on $\mathbf{p}$ and $\mathbf{q}$ will
be dropped when this will not compromise clarity.

Each agent autonomously computes a continuous time estimation
$\hat{I}_{i}(t)$, which is shared with the neighbors periodically during the
_update instants_ $t_{k}=kT$, with $k=\\{1,2,\dots\\}$. Given (11), it follows
that for $t\in\left[t_{k-1},t_{k}\right)$:

$I(t)=e^{\delta(t-t_{k-1})}I(t_{k-1})+\sum_{i\in\mathcal{T}}\int_{t_{k-1}}^{t}e^{\delta(t-\tau)}d_{i}(\tau)\text{d}\tau$
(35)

Since no communication occurs between update instants, agent $i$ estimates
$I(t)$ considering only its contribution to the task:

$\hat{I}_{i}(t)=e^{\delta(t-t_{k-1})}\hat{I}(t_{k-1})+\int_{t_{k-1}}^{t}\\!\\!e^{\delta(t-\tau)}d_{i}(\tau)\text{d}\tau,$
(36)

Then, at each update instant $t_{k}$ the agent communicates to the neighbors
its present estimation $\hat{I}_{i}(t_{k})$:

$\hat{I}_{i}(t_{k})=e^{T}\hat{I}_{i}(t_{k-1})+\tilde{d}_{i}(t_{k}),\quad\tilde{d}_{i}(t_{k})=\int_{t_{k-1}}^{t_{k}}\\!\\!e^{\delta(t-\tau)}d_{i}(\tau)\text{d}\tau$
(37)

After the agents exchange their estimates, the first correction is performed:

  * •

for all
$\mathbf{q}\in\operatorname{supp}\left\\{\tilde{d}_{i}(t_{k})\right\\}$

$\hat{I}_{i}^{-}(t_{k})=\hat{I}_{i}(t_{k})+\\!\\!\sum_{j\in\mathcal{N}_{i}(t_{k})}\\!\\!\max\left\\{0,\,\,\hat{I}_{j}(t_{k})-e^{T}\hat{I}_{i}(t_{k-1})\right\\}$
(38)

  * •

whereas for all
$\mathbf{q}\notin\operatorname{supp}\left\\{\tilde{d}_{i}(t_{k})\right\\}$

$\hat{I}_{i}^{-}(t_{k})=\max_{j\in\mathcal{N}_{i}(t_{k})}\left\\{\hat{I}_{j}(t_{k}),\,\,e^{T}\hat{I}_{i}(t_{k-1})\right\\}$
(39)

At this point, since for all
$\mathbf{q}\notin\operatorname{supp}\left\\{\tilde{d}_{i}(t_{k})\right\\}$
only the highest contribution between agent $i$’s coverage (which, by
definition, results from the time degradation of the one attained at time
$t_{k-1}$) and those of neighbors is considered, agent $i$ obtains a raw
coverage underestimation. In fact, potential overlappings between
$\tilde{d}_{j}(t_{k})$ and $\tilde{d}_{l}(t_{k})$, with
$j,l\in\mathcal{N}_{i}(t_{k})$, which yield a higher coverage level than the
one currently estimated for all
$\mathbf{q}\notin\operatorname{supp}\left\\{\tilde{d}_{i}(t_{k})\right\\}$,
would be ignored. However, such information can be easily retrieved through a
second correction, which improves each agent’s estimation. In fact, having
received $\hat{I}_{j}(t_{k})$ from neighbors, agent $i$ is able to compute the
portion of $\tilde{d}_{i}(t_{k})$ overlapped with those of neighbors. Let

$O_{i}(t_{k})=\bigcup_{j\in\mathcal{N}_{i}(t_{k})}\left\\{\mathbf{q}\in
Q\colon\hat{I}_{j}(t_{k})-e^{T}\hat{I}_{i}(t_{k-1})>0\right\\}$ (40)

be such overlapped area and define

$\tilde{d}_{i}^{o}(t_{k})=\left\\{\begin{array}[]{ll}\tilde{d}_{i}(t_{k}),&\mathbf{q}\in
O_{i}(t_{k})\\\ 0,&\mathbf{q}\notin O_{i}(t_{k})\end{array}\right.$ (41)

Neighbors then exchange their $\tilde{d}_{i}^{o}(t_{k})$ and perform the last
correction:

  * •

for all
$\mathbf{q}\in\operatorname{supp}\left\\{\tilde{d}_{i}(t_{k})\right\\}$

$\hat{I}_{i}(t_{k})=\hat{I}_{i}^{-}(t_{k}),$ (42)

  * •

whereas for all
$\mathbf{q}\notin\operatorname{supp}\left\\{\tilde{d}_{i}(t_{k})\right\\}$

$\hat{I}_{i}(t_{k})=\hat{I}_{i}^{-}(t_{k})+\\!\\!\sum_{j\in\mathcal{N}_{i}(t_{k})}\\!\\!\\!\tilde{d}_{j}^{o}(t_{k})-\max_{j\in\mathcal{N}_{i}(t_{k})}\\!\\!\tilde{d}_{j}^{o}(t_{k})$
(43)

The following result ensures the existence of regions where the estimation
error is zero:

######  Theorem 2.

Assume that $\hat{I}_{i}(0)=I(0)$ for all $\mathbf{q}\in Q$, and for all
$i\in\mathcal{T}$. In addition, assume that each agent can travel for a
maximum distance $\ell^{\text{max}}_{T}$ during a period $T$. Then, there
exists a region centered in the agent position at time $t_{k}$,
$\mathbf{S}\mathbf{p}_{i}(t_{k})$, and with radius
$r^{*}=R_{\text{com}}-\overline{r}^{\text{cov}}-(N-1)\ell^{\text{max}}_{T}$
denoted with:

$Z_{i}(t_{k})=\left\\{\mathbf{q}\in
Q\colon\|\mathbf{S}\mathbf{p}_{i}(t_{k})-\mathbf{q}\|\leq r^{*}\right\\}$ (44)

such that $\hat{I}_{i}(t_{k})=I(t_{k})$ for all $\mathbf{q}\in Z_{i}(t_{k})$,
and $i\in\mathcal{T}$.

###### Proof.

The proof consist in a transposition of Theorem IV.6 proof in [16], which is
proved in a discrete time scenario, to the continuous time case considered in
this paper. Theorem IV.6 in [16] states that each agent’s estimation is
correct within a distance
$R_{\text{com}}-\overline{r}^{\text{cov}}-(N-1)u^{\text{max}}$ at each
discrete time instant, with $u^{\text{max}}$ being the maximum allowable
control input norm for agents characterized by single-integrator kinematics
($\mathbf{x}(t+1)=\mathbf{x}(t)+\mathbf{u}(t)$), assuming that Assumptions 2
and 3 hold and that $\hat{I}_{i}(0)=I(0)$ for all $i\in\mathcal{T}$. In our
case, the maximum allowed travel distance $\ell^{\text{max}}_{T}$ during a
period $T$ is the continuous time analogous of the discrete time maximum
control input norm $u^{\text{max}}$ in [16] for single integrators. ∎

###  IV-C Decentralized Control Law

The decentralized control law is then defined as follows:

$\hat{\mathbf{u}}_{i}(t,\mathbf{x}_{i},\mathbf{p})=\hat{\mathbf{u}}_{\xi,i}(t,\mathbf{x}_{i},\mathbf{p}_{i})+\mathbf{u}_{v,i}(\mathbf{x}_{i},\mathbf{p})-\mu\mathbf{v}_{i},\,\,\mu>0$
(45)

where the coverage term is computed as follows:

$\hat{\mathbf{u}}_{\xi,i}(t,\mathbf{x}_{i},\mathbf{p}_{i})=-\beta\mathbf{G}_{i}(\mathbf{x}_{i})^{T}\frac{\partial\mathbf{h}_{i}(\mathbf{x}_{i})}{\partial\mathbf{x}_{i}}^{T}\\!\\!\frac{\partial\hat{\xi}_{i}(t,\mathbf{p}_{i})}{\partial\mathbf{p}_{i}}^{T}\\!\\!,\,\,\beta>0$
(46)

The task error index is estimated by each agents using the distributed
estimation discussed in Section IV-B:

$\frac{\partial\hat{\xi}_{i}(t,\mathbf{p}_{i})}{\partial\mathbf{p}_{i}}=-\int_{Q_{i}(\mathbf{p}_{i})}\\!\\!\frac{\partial
f(\hat{e}(t,\mathbf{q}))}{\partial\hat{e}(t,\mathbf{q})}\frac{\partial
d_{i}(\mathbf{p}_{i},\mathbf{q})}{\partial\mathbf{p}_{i}}\sigma\text{d}\mathbf{q}$
(47)

where the TEF estimation is given by:

$\displaystyle\hat{e}_{i}(t,\mathbf{q})$
$\displaystyle=\hat{d}_{*,i}(t,\mathbf{q})-\\!\\!\sum_{j\in\mathcal{N}_{i}(t)}\\!\\!d_{j}(\mathbf{p}_{j}(t),\mathbf{q})-d_{i}(\mathbf{p}_{i}(t),\mathbf{q})$
(48) $\displaystyle\hat{d}_{*,i}(t,\mathbf{q})$
$\displaystyle=\max\left\\{0,C^{*}-\hat{I}_{i}(t,\mathbf{q})\right\\}$ (49)

Observe that in the definition of $\hat{e}_{i}(\cdot)$ we assumed that agent
$i$ has the instantaneous knowledge of the ADFs of its neighbors. This is
possible as long as agents are able to exchange their poses and ADFs
parameters at a sufficiently high rate. Therefore, it should be pointed out
that two types of communication protocols need to be enforced. The first one,
allowing estimation updates as described in Section IV-B, must be guaranteed
with period $T$ and requires the exchange of high amounts of data, since each
agent should, at least, send its own estimated coverage map and receive those
of neighbors. The second one requires communication of quantities, such as
positions, orientations, and ADF-related parameters that are fast and easy to
handle and to exchange.

###  IV-D Comments on the Information Level Estimation

With reference to Theorem 2, the expression of $r^{*}$ implies that, if a
guarantee on the quality of the coverage estimation is to be sought,
communication updates should happen more frequently if vehicles are allowed to
move faster (and thus have larger $\ell^{\text{max}}_{T}$). Agents speed can
be reduced by choosing higher values of the damping $\mu$ in (22).

Moreover, in a team with a high number of agents, the region around each agent
where, at time $t_{k}$, its estimation is correct is sensibly reduced. The
reason is twofold: first, the cardinality $N$ of the team directly affects the
radius $r^{*}$; furthermore, the estimator definition in (36), which considers
only the agent contribution between the update instants, does not take into
account the contributions of the other agents, thus having higher chances of
errors along with a higher number of team members.

In fact, the level of information provided by (36) is an underestimation of
the real attained level:

$\hat{I}_{i}(t,\mathbf{q})\leq I(t,\mathbf{q}),\quad t\geq t_{0},\quad\forall
i\in\mathcal{T},\,\,\forall\mathbf{q}\in Q$ (50)

for the aforementioned reasons.

Having shown that the centralized control law guarantees the task fulfillment,
it should be pointed out that the decentralized scheme only induces a slight
performance degradation, without compromising the mission. Indeed, since (50)
holds, agent $i$ could be attracted by regions with an already high level of
information $I$, with a remarkable disparity $\hat{I}_{i}\leq I$. This would
be a reasonable trade-off due to the employment of a distributed estimation.
Conversely, agent $i$ would never be repulsed from regions with a very low
level of true attained coverage, since $\hat{I}_{i}$ would be low in those
regions as well. In addition, because of its definition,
$\hat{\xi}_{i}(t,\mathbf{p})\geq\xi(t,\mathbf{p})$ for all $i\in\mathcal{T}$,
and for $t\geq t_{0}$, thus remarking that with only locally achievable
information agent $i$ overestimates the error function, assuming that the task
is farther from being fulfilled than it truly is.

##  V Simulation Results

To validate the proposed distributed control law a persistent coverage
scenario in a $100\times 100$ square environment involving $6$ agents was
simulated. The information decay rate was set to $\delta=-0.1$ and
$C^{*}=\sigma(\mathbf{q})=1$ for all $\mathbf{q}\in Q$. Seven obstacles were
placed in the environment: four of them are virtual walls surrounding the
operational area, preventing agents from exiting it. The remaining ones are
polygons placed within $Q$. Three of the agents have double integrator
dynamics while the other three are dynamic unicycles with unitary mass and
moment of inertia (for the equations see [14, Chap. 11.4]). Having set $r=0.5$
and $R=3$, the initial poses were chosen such that
$\mathbf{p}(t_{0})\in\mathcal{P}(r)$. The double integrators carry isotropic
_Gaussian ADFs_:

$d_{G}(\mathbf{p}_{i},\\!\mathbf{q})\\!=\\!A\exp\left(-\frac{1}{2}\left\|\\!\begin{bmatrix}\cos\theta&\sin\theta\\\
-\sin\theta&\cos\theta\end{bmatrix}\\!(\mathbf{q}-\mathbf{S}\mathbf{p}_{i})\right\|^{2}_{\boldsymbol{\Sigma}^{-1}}\right)$
(51)

with $A=3$ and $\Sigma=\operatorname{diag}\left\\{3,3\right\\}$. The unicycles
are characterized by a _Gaussian ADF with limited field of view_:

$d_{G,\text{FOV}}(\mathbf{p}_{i},\mathbf{q})=d_{G}(\mathbf{p}_{i},\mathbf{q})f_{\text{FOV}}(\mathbf{p}_{i},\mathbf{q})$
(52)

with

$\displaystyle
f_{\text{FOV}}(\mathbf{p}_{i},\mathbf{q})=f_{\text{FOV},r}(\mathbf{p}_{i},\mathbf{q})f_{\text{FOV},l}(\mathbf{p}_{i},\mathbf{q})$
(53) $\displaystyle
f_{\text{FOV},r}(\mathbf{p}_{i},\mathbf{q})=\left(1+e^{-k\left(r_{1}\cos\left(\phi/2\right)+r_{2}\sin\left(\phi/2\right)\right)}\right)^{-1}$
(54) $\displaystyle
f_{\text{FOV},l}(\mathbf{p}_{i},\mathbf{q})=\left(1+e^{-k\left(-r_{1}\cos\left(\phi/2\right)+r_{2}\sin\left(\phi/2\right)\right)}\right)^{-1}$
(55) $\displaystyle
r_{1}=(q_{x}-x_{i})\cos\left(\theta_{i}-\pi/2\right)+(q_{y}-y_{i})\sin\left(\theta_{i}-\pi/2\right)$
(56) $\displaystyle
r_{2}=-(q_{x}-x_{i})\sin\left(\theta_{i}-\pi/2\right)+(q_{y}-y_{i})\cos\left(\theta_{i}-\pi/2\right)\\!\\!$
(57)

The function
$f_{\text{FOV}}:\mathcal{C}(Q)\times\mathbb{R}^{2}\rightarrow[0,1]$ shapes the
field of view and is equal to $1$ inside it, while smoothly decreasing to $0$
as its boundaries are approached. The parameter $k$ models the slope and was
set to $2$, while $\phi$ is the field of view angle and was set to $\pi/2$. As
regards the control gains, the following values were assigned:
$\alpha=\beta=30$, $\mu=20$. The communication range is $R_{\text{com}}=150$,
thus ensuring the connectivity of the graph $\mathcal{G}(t)$. The estimate
update period $T$ has been set to $1$.

Figure 1: Agents trajectories.

(a) Normalized true $\xi(t)$.

(b) Normalized agents estimation error, for $t\in[100,110]$.

Figure 2: Evolution of $\xi(t)$ and agents estimation error.

(a) Inter-agent distances.

(b) Agent-obstacle distances.

Figure 3: Inter-agent and agent-obstacle distances evolution.

Agents trajectories are shown in Figure 1. It can be seen that the persistent
task is efficiently performed and all the area is continuously covered by the
agents. The efficiency is confirmed by the evolution of the error index
$\xi(\cdot)$ shown in Figure 2a. The agents estimation error is shown in
Figure 2b. Note that $\hat{\xi}(t)\geq\xi(t)$ during the analyzed period, that
is agents overestimate the error index confirming what pointed out in Section
IV-D. The error grows between the updated instants, and becomes approximately
zero when updates are exchanged, proving the effectiveness of the distributed
estimation algorithm. Figures 3a and 3b confirm that the task is safely
executed, and no collisions occur.

##  VI Conclusions

We presented a distributed cooperative control strategy for persistent
coverage tasks execution. The proposed control law is designed for coping with
agents having different types of dynamics and heterogeneous sensing
capabilities. The latter was modeled and handled by means of the tools
provided by the descriptor function framework. Task execution safety was
guaranteed by means of obstacles and inter-agent collision avoidance
functions. An algorithm for the estimation of the attained level of coverage
was derived in order to decentralize the control law. The effectiveness of the
control law and of the estimation algorithm was formally proved and showed by
means of a numerical simulation.

## References

  * [1] H. Sahin and L. Guvenc, “Household robotics: autonomous devices for vacuuming and lawn mowing,” _IEEE Control Syst. Mag._ , vol. 27, no. 2, pp. 20–96, 2007.
  * [2] M. Dunbabin and L. Marques, “Robots for environmental monitoring: significant advancements and applications,” _IEEE Robot. Automat. Mag._ , vol. 19, no. 1, pp. 24–39, 2012.
  * [3] G. S. C. Avellar, G. A. S. Pereira, L. C. A. Pimenta, and P. Iscold, “Multi-UAV routing for area coverage and remote sensing with minimum time,” _Sensors_ , vol. 15, no. 11, pp. 27 783–27 803, 2015.
  * [4] P. F. Hokayem, D. M. Stipanović, and M. W. Spong, “On persistent coverage control,” in _Proc. 46th IEEE Conference on Decision and Control_ , New Orleans, LA, USA, Dec. 2007, pp. 6130–6135.
  * [5] C. Franco, G. Lopez-Nicolas, C. Sagues, and S. Llorente, “Persistent coverage control with variable coverage action in multi-robot environment,” in _Proc. 52nd IEEE Conference on Decision and Control_ , Florence, Italy, Dec. 2013, pp. 6055–6060.
  * [6] C. Franco, D. M. Stipanović, G. Lopez-Nicolas, C. Sagues, and S. Llorente, “Persistent coverage control for a team of agents with collision avoidance,” _European Journal of Control_ , vol. 22, pp. 30–45, 2015.
  * [7] C. Song, L. Liu, G. Feng, Y. Wang, and Q. Gao, “Persistent awareness coverage control for mobile sensor networks,” _Automatica_ , vol. 49, no. 6, pp. 1867–1873, 2013.
  * [8] J. M. Palacios-Gasós, Z. Talebpour, E. Montijano, C. Sagüés, and A. Martinoli, “Optimal path planning and coverage control for multi-robot persistent coverage in environments with obstacles,” in _Proc. 2017 IEEE International Conference on Robotics and Automation_ , Singapore, May 2017, pp. 1321–1327.
  * [9] Y. Wang and I. I. Hussein, “Awareness coverage control over large-scale domains with intermittent communications,” _IEEE Trans. Automat. Contr._ , vol. 55, no. 8, pp. 1850–1859, 2010.
  * [10] N. Zhou, X. Yu, S. B. Andersson, and C. G. Cassandras, “Optimal event-driven multi-agent persistent monitoring of a finite set of targets,” in _Proc. 55th IEEE Conference on Decision and Control_ , Las Vegas, NV, USA, Dec. 2016, pp. 1814–1819.
  * [11] G. Franzini and M. Innocenti, “Effective coverage control for teams of heterogeneous agents,” in _Proc. 56th IEEE Conference on Decision and Control_ , Melbourne, Australia, Dec. 2017, pp. 2836–2841.
  * [12] M. Niccolini, M. Innocenti, and L. Pollini, “Near optimal swarm deployment using descriptor functions,” in _Proc. 2010 IEEE International Conference on Robotics and Automation_ , Anchorage, AK, USA, 2010, pp. 4952–4957.
  * [13] M. Niccolini, “Swarm abstractions for distributed estimation and control,” Ph.D. dissertation, Univ. of Pisa, Pisa, Italy, 2011.
  * [14] B. Siciliano, L. Sciavicco, L. Villani, and G. Oriolo, _Robotics - Modelling, Planning and Control_. Springer-Verlag London, 2009.
  * [15] D. M. Stipanović, P. F. Hokayem, M. W. Spong, and D. D. S̆iljak, “Cooperative avoidance control for multiagent systems,” _J. Dyn. Sys., Meas., Control_ , vol. 129, pp. 699–707, 2007.
  * [16] J. M. Palacios-Gasós, E. Montijano, C. Sagüés, and S. Llorente, “Distributed coverage estimation and control for multirobot persistent tasks,” _IEEE Trans. Robot._ , vol. 32, no. 6, pp. 1444–1460, 2016.

  *[ADF]: agent descriptor function
  *[CTDF]: current task descriptor function
  *[TDF]: desired task descriptor function
  *[TEF]: task error function
  *[ ADF]: agent descriptor function


# Nonlinear dynamics of topological Dirac fermions in 2D spin-orbit coupled
materials

Rajesh K. Malla  Theoretical Division, Los Alamos National Laboratory, MS
B262, Los Alamos, New Mexico 87545, USA  Center for Nonlinear Studies, MS
B258, Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA
Wilton J. M. Kort-Kamp∗ Theoretical Division, Los Alamos National Laboratory,
MS B262, Los Alamos, New Mexico 87545, USA

###### Abstract

The graphene family materials are two-dimensional staggered monolayers with a
gapped energy band structure due to intrinsic spin-orbit coupling. The mass
gaps in these materials can be manipulated on-demand via biasing with a static
electric field, an off-resonance circularly polarized laser, or an exchange
interaction field, allowing the monolayer to be driven through a multitude of
topological phase transitions. We investigate the dynamics of spin-orbit
coupled graphene family materials to unveil topological phase transition
fingerprints embedded in the nonlinear regime and show how these signatures
manifest in the nonlinear Kerr effect and in third-harmonic generation
processes. We show that the resonant nonlinear spectral response of
topological fermions can be traced to specific Dirac cones in these materials,
enabling characterization of topological invariants in any phase by detecting
the cross-polarized component of the electromagnetic field. By shedding light
on the unique processes involved in harmonic generation via topological
phenomena our findings open an encouraging path towards the development of
novel nonlinear systems based on two-dimensional semiconductors of the
graphene family.

Graphene is the typical go-to material to investigate the optoelectronic
response of two-dimensional (2D) systems Neto2009 ; Peres2010  because of its
extraordinary electron mobility Graphenereview , tunable linear electronic
conductivity Tunable2011 , and potential for strong-light matter interactions
at sub-wavelength scales Grapheneplasmons . Although it supports quantum Hall
states in the presence of strong magnetic fields QHall1 ; QHall2 ; QHall3 ,
the intrinsically weak spin-orbit coupling WeakSOC  severely limits graphene’s
suitability to study topological phase transitions TI1 ; TI2  in low-
dimensional materials. The search for alternative topologically non-trivial 2D
structures with characteristics similar to graphene recently ended with the
synthesis of silicene Silicene , germanene Germanene , stanene Stanene , and
plumbene Plumbene . The newer members of the graphene family materials (GFM)
have a buckled honeycomb lattice with silicon, germanium, tin, or lead atoms
occupying the lattice sites Gomez2016 ; Molle:2017aa ; Mannix:2017aa . Unlike
graphene, the strong spin-orbit coupling in these systems results in a gapped
bulk energy band-structure and protected one-way edge states characteristic of
topological insulators. Strikingly, they can be driven through a variety of
phase transitions (Fig. 1a,b) by actively controlling the mass gap via
external interactions EzawaJapanese , which could enable unprecedented all-in-
one material multi-optoelectronic functionalities with potential applications
to spintronics and valleytronics.

Various recent works have investigated ultrafast and nonlinear effects in
topologically trivial 2D materials Gnonlinear1 ; Mikhailov2008 ; Kerr ;
Gnonlinear2 ; Transition ; Exp2 ; Exp3 ; Exp4 ; Exp5 ; Sipe2014 ; Sipe2015 ;
Pedersen2016 ; Mikhailov2016 ; HHG1 ; Lysne2020 ; Linear1  and in topological
systems HHG4 ; HHG5 ; HHG6 ; Chacon2020 ; RSeq . Unravelling the interplay
between topology and nonlinear effects in spin-orbit coupled monolayer
semiconductors of the graphene family is a natural step at the materials
science forefront, which could aid in developing next-generation technologies
that meet the urgent demands for, among others, higher performance radio-
frequency modulators, optically gated transistors, and practical spintronic-
based devices. Nevertheless, studies on the optoelectronic properties of these
materials have largely focused in their linear responseTunablebandgap1 ;
Tunablebandgap2 ; WiltonPRB ; WiltonPRL ; WiltonNonlocal ; Circular3 ;
Wiltonnature ; Tabert2014 ; Farias2018 ; Wu2018 ; Wu2020 , which include
investigations of signatures of the tunable band gap Tunablebandgap1 ;
Circular3 , spatial dispersion Tabert2014 ; WiltonNonlocal , the interplay
between the quantum Hall effect Tunablebandgap2  and photo-induced topology
WiltonPRB , as well as of topological phase transitions in quantum forces
Wiltonnature ; Farias2018 , spin-orbit photonic interactions WiltonPRL , and
light beam shifts Wu2018 ; Wu2020 . The crossroads between nonlinear dynamics
and topology in the extended graphene family, potentially allowing access to
topological phase transition signatures, material symmetries, selection rules,
and relaxation mechanisms otherwise screened by spurious effects in the linear
response, remains uncharted.

We explore the interplay between topological Dirac fermions in the GFM and
optical fields beyond the linear regime. We develop a comprehensive and
unified description of their nonlinear dynamical conductivity tensor,
including non-trivial topological phases due to an externally applied
electrostatic field, a circularly polarized laser, and exchange interaction.
Our results show a rich structure of nonlinear optoelectronic responses across
the phase diagram, and we demonstrate that the third-order optoelectronic
conductivity contribution to the fundamental and third-harmonics encodes
fingerprints of topological phase transitions. We show that characterization
of the Chern number is possible near bandgap resonances via polarization
sensitive detection of the cross-polarized component of the field harmonics.
Our work sets the cornerstones for investigating topological phase transitions
beyond the linear response in the family of graphene-like elemental 2D
materials.

Results

Figure 1: Topological phases in the graphene family. (a) Schematic
representation of a graphene family monolayer under the influence of a static
electric field ($E_{z}$), an off-resonance circularly polarized laser with
intensity $I_{C}$, and the staggered exchange field $\lambda_{M}$. The in-
plane lattice constant is $a$ and the buckling distance between sites $A$ and
$B$ is $2l$. (b) 3D topological phase diagram and corresponding Chern numbers
(${\cal C}$, ${\cal C}_{s}$, ${\cal C}_{\eta}$, ${\cal C}_{s\eta}$) in spin-
orbit coupled graphene family materials. The momentum dependent meron
structures for a given spin are shown in (c) and (d) for valleys $K$ and (d)
$K^{\prime}$, respectively.

Topological properties of the graphene family -  Let us consider the following
Dirac-like Hamiltonian describing staggered GFM in the low-energy regime
EzawaJapanese ,

${\hat{H}}_{s}^{\eta}={\bm{\tau}}\cdot{\bf d}^{\eta}_{s}=\hbar v_{F}(\eta
k_{x}\tau_{x}+k_{y}\tau_{y})+\Delta_{s}^{\eta}\tau_{z},$ (1)

where ${\bf d}^{\eta}_{s}=\\{\eta\hbar v_{F}k_{x},\hbar
v_{F}k_{y},\Delta_{s}^{\eta}\\}$ describes a meron structure in momentum space
(Fig. 1c,d), $\hat{\textbf{p}}=\hbar(k_{x},k_{y})$ is the particle momentum,
$\tau_{i}$ are the Pauli matrices, $\eta,s=\pm 1$ are valley and spin indexes,
$v_{F}$ is the Fermi velocity, and the mass term $\Delta_{s}^{\eta}=\eta
s\lambda_{SO}-\lambda_{E}-\eta\lambda_{L}+s\lambda_{M}$ corresponds to half
energy band gap. Here, $\lambda_{SO}$ represents the intrinsic spin-orbit
coupling energy and the last three terms in $\Delta_{s}^{\eta}$ account for
interactions with external fields, which allow for tailoring the Dirac mass
gap. Indeed, the second term $\lambda_{E}$ corresponds to the potential
difference between sub-lattices in the buckled structure in the presence of a
static electric field $E_{z}$ applied perpendicularly to the plane of the
monolayer Tunablebandgap1 . The third component $\lambda_{L}$ describes the
anomalous quantum Hall effect and arises due to the coupling between an off-
resonant circularly polarized laser and the GFM Circular2 ; Tunablebandgap2 .
The final term $\lambda_{M}$ depicts the staggered antiferromagnetic exchange
interaction Antiferromagnetic . For simplicity, we have neglected Rashba
couplings due to their small magnitude compared to the other contributions. It
is worth mentioning that the Hamiltonian in Eq. (1) applies more generally to
other 2D systems, including antiferromagnetic manganese chalcogenophosphates
(MnPX3, X = S, Se) Li  and perovskites Liang . The formalism developed here is
also valid for these systems.

Within the Dirac picture the electronic properties of the GFM can be fully
characterized through a set of four topological invariants EzawaJapanese ,
namely the Chern ${\cal C}=\sum_{\eta,s}{\cal C}_{s}^{\eta}$, spin Chern
${\cal C}_{s}=\sum_{\eta,s}s\leavevmode\nobreak\ {\cal C}_{s}^{\eta}/2$,
valley Chern ${\cal C}_{\eta}=\sum_{\eta,s}\eta\leavevmode\nobreak\ {\cal
C}_{s}^{\eta}$, and spin-valley Chern $C_{s\eta}=\sum_{\eta,s}\eta
s\leavevmode\nobreak\ {\cal C}_{s}^{\eta}/2$ numbers. Here ${\cal
C}_{s}^{\eta}=\eta\leavevmode\nobreak\ \text{sign}[\Delta_{s}^{\eta}]/2$ is
the Pontryagin number TI1 ; TI2 , a topological quantity that counts how many
times the vector ${\bf d}^{\eta}_{s}/|{\bf d}^{\eta}_{s}|$ wraps a sphere in
momentum space (Fig. 1c,d). By varying the external parameters $\lambda_{E},\
\lambda_{L}$, and $\lambda_{M}$ the monolayer can be driven through various
phase transitions as depicted in Fig. 1b, where we show the topological phases
in the planes $\lambda_{E}=0$, $\lambda_{L}=0$, and $\lambda_{M}=0$. The
electronic states quantum spin Hall insulator (QSHI), anomalous quantum Hall
insulator (AQHI), band insulator (BI), antiferromagnetic insulator (AFI), and
polarized- spin quantum Hall insulator (PS-QHI) have non-zero mass gaps for
all four Dirac cones. The system undergoes a topological phase transition when
the band gap of at least one Dirac cone changes sign. Hence, the boundaries
between insulating states are determined by the condition
$\Delta^{\eta}_{s}=0$, which defines the so called single Dirac cone phases
(diagonal solid lines in Fig. 1b). At the intersection between boundary lines
there are two Dirac gaps closed, and the system is either in a spin, valley,
or spin-valley polarized semimetal state. Finally, note that there are four
points in the 3D phase diagram,
$\\{\lambda_{E}/\lambda_{SO},\lambda_{L}/\lambda_{SO},\lambda_{M}/\lambda_{SO}\\}=\\{-1,-1,1\\},\
\\{1,1,1\\},\ \\{1,-1,-1\\},\ \\{-1,1,-1\\}$, where three Dirac cones
simultaneously close, corresponding to unique electronic states that remain
unexplored to date.

Nonlinear light-matter interactions in the GFM - The nonlinear dynamics of
charge carriers in a given Dirac cone interacting with a linearly polarized
electromagnetic plane wave impinging normally on the monolayer can be
described via the density matrix ${\hat{\rho}}(t)$, which satisfies the
equation of motion $i\hbar\ \partial\hat{\rho}(t)/\partial
t=\left[{\hat{H}}_{s}^{\eta}+\hat{H}_{i}(t),\hat{\rho}(t)\right]-\Gamma({\hat{\rho}(t)}-{\hat{\rho}}^{(0)}).$
Here, $\hat{H}_{i}(t)=e{\hat{\bm{p}}}\cdot{{\bm{A}}}(t)/mc$ is the field-
monolayer interaction Hamiltonian in the velocity gauge Sipevelocity , where
${\bm{A}}(t)=A(t)\hat{\textbf{x}}$ is the electromagnetic vector potential. We
assume that at large times the system relaxes to the equilibrium density
matrix
$\hat{\rho}^{(0)}=\sum_{l,{\bm{k}}}f_{l{\bm{k}}}|l{\bm{k}}\rangle\langle
l{\bm{k}}|$ with a phenomenological decay rate $\Gamma$, where $f_{l{\bm{k}}}$
is the Fermi-Dirac distribution for fermions in the valence or conduction
bands ($l=\pm 1$) with momentum $\hbar\textbf{k}$, and $|l{\bm{k}}\rangle$ are
the corresponding eigenstates of the unperturbed Hamiltonian $H^{\eta}_{s}$
(see Methods). In general, the equation of motion cannot be solved
analytically due to the time dependence of the interaction Hamiltonian, and
numerical methods are often employed. Here, instead, we use perturbation
theory by expanding the density matrix
${\hat{\rho}}(t)=\sum_{n}{\hat{\rho}}^{(n)}(t)$ in powers of the vector
potential amplitude and iteratively solve for $\rho^{(n)}(t)\propto A(t)^{n}$.
We use the standard expression for total current
${\bm{j}}(t)=-e{\text{T}r}[{\hat{\rho}}(t)({\hat{\bm{v}}}+\frac{e}{mc}{\hat{\bm{A}}}(t))]$,
which contains contributions arising from all orders in the perturbative
expansionPedersen2017 ; Aversa1995 . Here, $\hat{v}=\hbar^{-1}\nabla_{{\bf
k}}\hat{H}_{s}^{\eta}$ is the velocity operator. The $n$-th order ac electric
current is given by

${\bm{j}}^{(n)}(t)=-e\textrm{Tr}[{\hat{\rho}}^{(n)}(t){\hat{\bm{v}}}]-\frac{e^{2}}{mc}\textrm{Tr}[{\hat{\rho}}^{(n-1)}(t){\bm{A}}(t)].$
(2)

Although the last term in the previous equation does not vanish in general,
one can demonstrate that it is exactly canceled by another contribution
arising from the first term and, therefore, its contributions do not appear in
the final expressions for the conductivities. The cross-canceling of
contributions involving the first and second terms in Eq. (2) is actually
critical to reobtain the well know expression for the conductivity in the
linear regime, given by Kubo’s formalism (Methods). The $n$-th order optical
conductivity is a tensor of rank $n+1$, accounts for both intra and interband
transition contributions as well as topological effects emerging from the
Berry connection EzawaJapanese , and can be directly computed after Fourier
transforming the above equation and expressing the current as a product
between the electric field and conductivity.

Figure 2:  Identifying topological properties of individual Dirac cones in the
linear regime. Real (solid) and imaginary (dashed) components of the charge
$\sigma_{yx}^{(1)}$ (blue), spin $\sigma_{yx}^{(1),s}$ (red), valley
$\sigma_{yx}^{(1),\eta}$ (black), and spin-valley $\sigma_{yx}^{(1),s\eta}$
(green) transverse linear conductivities are shown for the four points in
phase space where only when Dirac cone is open, namely
$\\{\lambda_{E}/\lambda_{SO},\lambda_{L}/\lambda_{SO},\lambda_{M}/\lambda_{SO}\\}$
equal to (a) $\\{-1,-1,1\\}$ (b) $\\{1,1,1\\}$ (c) $\\{1,-1,-1\\}$ (d)
$\\{-1,1,-1\\}$. The insets at the top of each panel show the corresponding
structure of Dirac cones. In panel (a) the charge and spin conductivities are
identical and in panel (b) the same occurs for the valley and spin-valley
conductivities. All conductivities are expressed in units of
$2\sigma_{0}/\pi$, where $\sigma_{0}=\alpha c/4$ and $\alpha$ is the fine-
structure constant, we used $\hbar\Gamma=0.05\lambda_{SO}$, and assumed the
monolayer to be neutral.

The longitudinal $\sigma_{xx}^{(1)\eta,s}(\omega)$ and transverse
$\sigma_{yx}^{(1)\eta,s}(\omega)$ linear optical conductivities per Dirac cone
obtained from Eq. (2) are in full accordance with previous studies, and have
been investigated in details in the $\lambda_{M}=0$ plane in Refs. 33, 41, 49,
47. For completeness, we explore in Fig. 2 the charge, spin, valley, and spin-
valley transverse linear conductivities (summed over valley and spin indexes
as shown in the Methods) in the four phases with only one Dirac cone open,
since these have been previously overlooked in the literature. By measuring
the dynamic linear conductivities at these unique points in phase space one
can investigate topological properties emerging from each Dirac cone
individually. Indeed, the nonzero Dirac gap in each of the cases considered
leads to Chern numbers ${\cal C}$, ${\cal C}_{\eta}$ $=\pm 1/2$, and ${\cal
C}_{s}$, ${\cal C}_{s\eta}$ $=\pm 1/4$, which clearly manifest in the zero
frequency limit of the transverse conductivities for neutral monolayers, as
shown in Fig. 2. Note that in panels 2a,b the charge (spin) and the valley
(spin-valley) currents are identical, which implies that the nonzero Dirac gap
belongs to the $K$ valley, while in panels 2c,d the charge (spin) and valley
(spin-valley) conductivities are the negative of each other, hence the nonzero
Dirac gap belongs to $K^{\prime}$ valley. Finally, the sign of the spin or
spin-valley linear conductivities in the quasi-static regime allows to
identify the spin of the charge carriers.

The second order response of the system described by the Hamiltonian in Eq.
(1) vanishes due to the centrosymmetry of the GFM. The first nonlinear
contribution comes from the third-order response and can be described by a
rank four conductivity tensor
$\tilde{\sigma}_{\alpha\beta\gamma\delta}^{(3)\eta,s}(\omega_{1},\omega_{2},\omega_{3})$,
which is invariant under simultaneous permutations of the indexes
$\beta,\gamma,\delta$ and frequencies $\omega_{1},\omega_{2},\omega_{3}$,
where $\alpha,\beta,\gamma,\delta=x,y$ (see Methods). For simplicity we will
concentrate in cases where $\beta=\gamma=\delta=x$ and we will assume that the
incident electromagnetic radiation is monochromatic, i.e.,
${\bm{A}}(t)=A_{0}\cos(\omega t)\leavevmode\nobreak\ \hat{\textbf{x}}$,
allowing us to separated the third-order conductivity into terms oscillating
with frequencies $\omega$ and $3\omega$. The first term adds to the linear
conductivity response and gives a correction to the fundamental harmonic that
is quadratic in the vector potential amplitude (Kerr effect), while the second
one describes the nonlinear process of third harmonic generation. Therefore,
up to third-order in perturbation theory, the optoelectronic conductivity of
the system at the fundamental and third-harmonics due to each Dirac mass gap
can be cast as

$\displaystyle\tilde{\sigma}_{\alpha x}^{\eta,s}(\omega,I_{0})$
$\displaystyle=$ $\displaystyle\sigma_{\alpha
x}^{(1)\eta,s}(\omega)+\frac{6\pi}{c}I_{0}\sigma_{\alpha
x}^{(3)\eta,s}(\omega),$ $\displaystyle\tilde{\sigma}_{\alpha
x}^{\eta,s}(3\omega,I_{0})$ $\displaystyle=$
$\displaystyle\frac{2\pi}{c}I_{0}\sigma_{\alpha x}^{(3)\eta,s}(3\omega)\,,$
(3)

where $I_{0}=\omega^{2}A_{0}^{2}/8\pi c^{2}$ is the incident field intensity
and $\tilde{\sigma}_{\alpha
x}^{(3)\eta,s}(\omega)\\!\\!=\\!\\!\tilde{\sigma}_{\alpha
xxx}^{(3)\eta,s}(\omega,\omega,-\omega)$, $\tilde{\sigma}_{\alpha
x}^{(3)\eta,s}(3\omega)\\!\\!=\\!\\!\tilde{\sigma}_{\alpha
xxx}^{(3)\eta,s}(\omega,\omega,\omega)$ is a shorthand notation. In the next
section we discuss topological fingerprints embedded in these nonlinear
conductivities.

Topological signatures in the GFM nonlinear response - Let us first look into
topological signatures buried in the phase of the optical field by
investigating the polarization state (helicity) of light resulting from third-
harmonic generation processes. We consider the low temperature regime for a
neutral monolayer, i.e., the Fermi level lies in the middle of the band gap,
and assume $\hbar\omega=0.2\lambda_{SO}$ and $\hbar\Gamma=0.05\lambda_{SO}$.
In Fig. 3a we plot the phase diagram associated with the difference
$h_{3\omega}=(I_{L}^{(3)}-I_{R}^{(3)})/I_{0}\propto\text{Im}\left[\sigma_{xx}^{(3)}(3\omega)^{*}\sigma_{yx}^{(3)}(3\omega)\right]$
between intensities $I_{L}^{(3)}$ and $I_{R}^{(3)}$ of light emitted with left
and right circular polarization at frequency $3\omega$, respectively, where
here the conductivities $\sigma_{xx}^{(3)}(3\omega)$,
$\sigma_{yx}^{(3)}(3\omega)$ account for all Dirac cones (Methods).

Figure 3: Topological phase transitions in the third-harmonic polarization
state. (a) Phase diagram of the difference between intensities of third-
harmonic emission with left and right circular polarization $h_{3\omega}$
(a.u.). (b) 2D phase diagrams of $h_{3\omega}$ at various horizontal planes
with fixed value of $\lambda_{M}$. Here, $\hbar\omega=0.2\lambda_{SO}$ and
$\hbar\Gamma=0.05\lambda_{SO}$. Figure 4:  Nonlinear charge, spin, and valley
dynamics in the graphene family materials. Panels (a)-(d) show the third-
harmonic generation for charge (blue), spin (black), and valley (red)
transverse conductivities along various topological phase transitions paths as
highlighted in the (e) 2D phase diagram of $\sigma_{yx}^{(3)}(3\omega)$ for
$\lambda_{M}=0$ plane. While (a)-(c) paths lie on the plane $\lambda_{M}=0$,
(d) is along the diagonal of the 3D phase diagram connecting points
($\lambda_{E}/\lambda_{SO},\lambda_{L}/\lambda_{SO},\lambda_{M}/\lambda_{SO})=(0,0,0)\rightarrow(2,2,2)$.
The paths are parametrized in terms of $\lambda_{E}/\lambda_{SO}$, the
conductivities are expressed in the units of
$\alpha^{2}\hbar^{3}c^{2}v_{F}^{2}/\lambda_{SO}^{4}$, and the frequency and
relaxation rates are the same as in Fig. 3.

A direct comparison with the phase diagram in Fig. 1 clearly reveals the role
of topology in the nonlinear optical response of the GFM. Far from the phase
boundaries the contribution from each Dirac gap to $h_{3\omega}$ decreases as
$1/|\Delta_{s}^{\eta}|^{4}$, resulting in weak variations of helicity in the
middle of the QSHI, BI, AQHI, PS-QHI, and AFI phases.

The effects of nonlinearity are most relevant as one approaches the semi-
metallic states (white dashed lines) marking the boundary between different
Chern numbers, where $h_{3\omega}$ is dominated by the contribution from the
Dirac cone with the smallest gap and the phase transition can be easily
identified. These conclusions are not restricted to the planes
$\lambda_{E}=0$, $\lambda_{L}=0$, and $\lambda_{M}=0$ as can be seen in Fig.
3b, where we present slices of the phase diagram corresponding to fixed values
of $\lambda_{M}$. Note in these panels that as $\lambda_{M}$ increases the
phase boundaries come closer to each other and at $\lambda_{M}/\lambda_{SO}=1$
the three phase boundaries cross at a single point,
$\\{\lambda_{E}/\lambda_{SO},\lambda_{L}/\lambda_{SO},\lambda_{M}/\lambda_{SO}\\}=\\{1,1,1\\}$,
corresponding to the closing of three Dirac gaps. Finally, we mention that the
variations of $h_{3\omega}$ across a phase transition depend on whether the
initial and final phases have trivial charge transport properties (Chern
number ${\cal{C}}$ equal to zero). For instance, in the $\lambda_{M}=0$ we
observe that $h_{3\omega}$ is a symmetric function around
$\lambda_{E}/\lambda_{SO}=1$ when we consider the topological phase transition
between the QSHI and BI phases, both with ${\cal{C}}=0$. On the other hand, a
remarkably asymmetric behavior takes place around $\lambda_{L}/\lambda_{SO}=1$
(which is energetically equivalent to the previous case) when one considers a
topological phase transition involving the QSHI and the AQHI (${\cal{C}}=-2$)
phases.

The phase diagram of $h_{3\omega}$ includes contributions from both third-
harmonic longitudinal and transverse conductivities. However, only the latter
encodes topological features of the GFM and deserves to be independently
studied. To this end one could, for instance, use a linear polarizer filter to
block the co-polarized scattered field and detect only its cross-polarized
component, which is proportional to $\sigma_{yx}^{(3)}$. In Fig. 4 we
investigate the charge conductivity $\sigma_{yx}^{(3)}(3\omega)$ along various
paths in the phase diagram (parametrized in terms of
$\lambda_{E}/\lambda_{SO}$). For completeness we also show the associated spin
$\sigma_{yx}^{(3),\textrm{s}}(3\omega)$ and valley
$\sigma_{yx}^{(3),\eta}(3\omega)$ third harmonic transverse conductivities.

In panel 4a we plot the conductivities along the $\lambda_{E}$ axes as we move
from the QSHI to the BI phase by closing two Dirac gaps when
$\lambda_{E}/\lambda_{SO}=1$. In this case the Chern numbers
(${\cal{C}},{\cal{C}}_{s},{\cal{C}}_{\eta}$) change from ($0,1,0$) to
($0,0,-2$) and we observe that the charge conductivity is always zero. The
resonant peaks and dips near $\lambda_{E}/\lambda_{SO}=1$ in the spin and
valley conductivities divulge information about the change in ${\cal{C}}_{s}$
and ${\cal{C}}_{\eta}$ across the phase transition. For instance, as the
system is driven from the QSHI to the BI phase the resonance at
$3\hbar\omega=2\Delta_{s}^{\eta}$ (dashed vertical lines) in
$\sigma_{yx}^{(3),\textrm{s}}(3\omega)$ and $\sigma_{yx}^{(3),\eta}(3\omega)$
flips from negative to positive values, which is a signature indicating a
decrease in the corresponding Chern numbers across the phase transition. Also,
the magnitude of the resonance peaks and dips in the valley conductivity is
twice as those of the spin conductivity, which coincides with the fact that
the change in ${\cal{C}}_{\eta}$ is two times that in ${\cal{C}}_{s}$. These
features are general and can be applied to any paths in the phase diagram
regardless of the nature of the transition involved.

In panel 4b we plot the conductivities as the system moves from the AQHI to
the BI phase via the PS-QHI phase. Along this path the monolayer crosses two
phase boundaries, with Chern numbers
(${\cal{C}},{\cal{C}}_{s},{\cal{C}}_{\eta}$) changing as $(-2,0,0)$
$\rightarrow$ ($-1,1/2,-1$)$\rightarrow$ ($0,0,-2$). Note that the signs of
the resonances at $3\hbar\omega=2\Delta_{s}^{\eta}$ for the valley
conductivity are opposite to those for the charge conductivity, while their
strengths are the same, meaning that the variations in ${\cal C}$ and ${\cal
C}_{\eta}$ across the phase transition are equal in magnitude. The sign of the
third harmonic resonance for the charge (valley) conductivity changes from
negative (positive) in the AQHI phase to positive (negative) in the PS-QHI
phase, which agrees with an increase (decrease) in ${\cal{C}}$
(${\cal{C}}_{\eta}$). The same feature repeats near the phase transition
between the PS-QHI and BI phases, signaling another increase (decrease) in
${\cal{C}}$ (${\cal{C}}_{\eta}$). Moreover, the behavior of the spin
conductivity in Fig. 4b reflects the fact that ${\cal C}_{s}$ increases along
the first phase boundary and then decreases along the second transition,
further confirming that the nonlinear resonant behavior of the conductivities
encode topology fingerprints.

In the third panel of Fig. 4 the system is driven along the phase boundaries
highlighted in path (c), where at least one of the Dirac gaps is closed all
the time and the Chern numbers change according to $(-3/2,1/4,-1/2)$
$\rightarrow$ ($-1/2,3/4,-1/2$) $\rightarrow$ ($-1/2,1/4,-3/2$). It is
important to notice that ${\cal C}$ does not change near
$\lambda_{E}/\lambda_{SO}=1$ and ${\cal C}_{\eta}$ remains the same across the
transition at $\lambda_{L}/\lambda_{SO}=1$ ($\lambda_{E}=0$ in the
parametrized curve). This fact is clearly reflected in Fig. 4c, where the
third harmonic resonances associated with charge (valley) conductivity have
the same signs as the system is driven through $\lambda_{E}/\lambda_{SO}=1$
($\lambda_{L}/\lambda_{SO}=1$). Finally, in Fig. 4d we plot the conductivities
for a path (not shown in Fig. 4e) going from the origin of the phase diagram
to the point
$\\{\lambda_{E}/\lambda_{SO},\lambda_{L}/\lambda_{SO},\lambda_{M}/\lambda_{SO}\\}=\\{2,2,2\\}$
across the diagonal $\lambda_{E}=\lambda_{L}=\lambda_{M}$. Once again, the
behavior of the conductivities and the change in Chern numbers agrees well
with signatures we explained above.

Next we show how one can identify signatures of the Chern numbers in each
topological phase by investigating the nonlinear spectral response of the GFM.
In Fig. 5 we present the frequency dispersion of the real part of the third-
order optoelectronic longitudinal and transverse conductivities for both
fundamental and third harmonics at four different points in phase space. These
correspond to cases where all gaps are open (Fig. 5a) or one (Fig. 5b), two
(Fig. 5c), or three (Fig. 5d) Dirac gaps are closed. The corresponding curves
for the imaginary part (not shown) of the conductivities can be obtained
directly from our formalism or via generalized Kramers-Kronig relations for
nonlinear systems Hutchings1992 . We expect that the nonlinear conductivities
should have non-negligible contributions when at least one Dirac gap matches
one or three times the energy of the incident photons. Indeed, Fig. 5 shows
that Kerr $\sigma_{\alpha x}^{(3)}(\omega)$ and third harmonic $\sigma_{\alpha
x}^{(3)}(3\omega)$ conductivities present localized spectral resonances at
$2\Delta_{s}^{\eta}=\hbar\omega$ (dashed vertical gridlines) and
$2\Delta_{s}^{\eta}=3\hbar\omega$ (dotted vertical gridlines), respectively,
with richer and more involved spectral responses in the cases where more Dirac
cones are open. We mention that although the third-harmonic generation
conductivity also includes resonances near $\omega$ and $2\omega$, they are
difficult to identify in the spectrum because their magnitude are
significantly smaller than the resonances near $3\omega$.

Note that near resonances the longitudinal Kerr and third harmonic
conductivities always shows normal and anomalous dispersion, respectively,
regardless of the resonant Dirac cone. This is due to the fact that both of
these conductivities depend only on the absolute value of the resonant mass
gap, hence not allowing to directly distinguish between energetically
equivalent but topologically different points in the phase space. On the other
hand, the transverse Hall conductivity is sensitive to the sign of the
resonant $\Delta_{s}^{\eta}$ and, therefore, encodes information of the
corresponding Pontryagin number ${\cal C}_{s}^{\eta}$. Note in Fig. 5, for
instance, that $\sigma_{yx}^{(3)}(\omega)$ and $\sigma_{yx}^{(3)}(3\omega)$
have positive and negative resonances depending on the resonant Dirac gap.
Therefore, just as the linear Hall conductivity allows to probe topological
properties of the GFM at finite frequencies WiltonPRB , so it does at the
nonlinear regime. Indeed, the Chern number associated to the open Dirac cones
at any point in phase space can be computed by summing the sign of
$\sigma_{yx}^{(3)}(\omega)$ or $\sigma_{yx}^{(3)}(3\omega)$ precisely at the
resonant frequencies, accounting appropriately for degeneracy among mass gaps,
and multiplying the result by $+1/2$ or $-1/2$, respectively. Consider, for
example, the case in Fig. 5a. Note that in this case the third-order
conductivity at the fundamental (third) harmonic evaluates to negative
(positive) values in three of the four resonances. Thus, the sum of their
signs equals $-2$ ($+2$) and results in a Chern number ${\cal C}=-1$, as
expected for the PS-QHI phase in the $\lambda_{M}=0$ plane. Note that this
argument can be extended to the nonlinear spin, valley, and spin-valley Hall
currents, allowing for obtaining the full set of Chern numbers characterizing
the topological phases in buckled two-dimensional semiconductors. In addition,
unlike the linear case where to obtain the Chern number one needs to measure
the linear Hall conductivity at various points around a resonance to obtain
the sign of its derivative WiltonPRB , here we need to evaluate
$\sigma_{yx}^{(3)}(3\omega)$ at most at the four (in the non-degenerate case)
resonant frequencies. Finally, we mention that the dependence of the
transverse (longitudinal) conductivities on $\textrm{sign}[\Delta_{s}^{\eta}]$
($|\Delta_{s}^{\eta}|$) are not exclusive to the third order nonlinear case
and similar conclusions should hold for higher order contributions as well.

Figure 5:  Nonlinear spectral response of spin-orbit coupled graphene family
materials. Real part of the third-order optical conductivities corresponding
to Kerr effect (black, frequency $\omega$) and third harmonic generation (red,
frequency $3\omega$) in both longitudinal (solid, $\alpha=x$) and and
transverse (dashed, $\alpha=y$) directions for
$\\{\lambda_{E}/\lambda_{SO},\lambda_{L}/\lambda_{SO},\lambda_{M}/\lambda_{SO}\\}$
equal to (a) $\\{1.5,1.2,0\\}$, (b) $\\{0.7,0.3,0\\}$, (c) $\\{0,1,0\\}$, and
(d) $\\{1,1,1\\}$. The dashed (dotted) vertical gridlines indicate the
resonance conditions for the third (fundamental) harmonic. Conductivities are
expressed in units of
$\alpha^{2}\hbar^{5}\omega^{2}c^{2}v_{F}^{2}/\lambda_{SO}^{6}$ and we assume
$\hbar\Gamma/\lambda_{SO}=0.1$.

Conclusions

In summary, we developed a comprehensive investigation of nonlinear light-
matter interactions in spin-orbit coupled graphene materials, with focus on
third harmonic generation and nonlinear corrections to the fundamental
harmonic (Kerr effect). We have shown that by controlling various external
parameters the monolayer can be driven through different phase transitions to
explore the interplay between nonlinear effects and topological properties of
two-dimensional materials. Specifically, we demonstrated that the helicity of
the nonlinear scattered field encodes information about the location of the
phase transitions boundaries, the dependence of the cross-polarized component
of the field (proportional to the transverse conductivity) on the external
parameters near resonances divulges details about the topological phases
involved in the transition, and the dispersive nonlinear response of the
system enables characterization of the Chern number in each phase. The recent
progress in synthesis of topological semiconductors of the graphene family
materials Silicene ; Germanene ; Stanene ; Plumbene  and well established
nonlinear characterization photonic techniques Kerr ; Gnonlinear2 ; Exp2 ;
Exp3 ; Exp4 ; Exp5  suggest that our findings can be accessed experimentally
with current technologies. We envision that the the effects predicted here
will greatly impact research at the crossroads between nonlinear optics,
topological materials, spintronics, and valleytronics.

Methods

In order to compute the nonlinear conductivity we need the unperturbed
eigenvectors and the velocity matrix elements associated to the Hamiltonian in
Eq. (1). In the following we will omit the valley and spin indexes $\eta$ and
$s$ whenever possible, and the reader should keep in mind that, unless
otherwise stated, all results presented are valid for a single Dirac cone. The
total conductivities are given by summing the contributions of all mass gaps
in the band structure. The eigenvector associated to the Hamiltonian in Eq.
(1) has the form

$|l\bm{k}\rangle=\frac{\hbar v_{F}|{\bf
k}|}{\left[2\epsilon\left(\epsilon+l\Delta\right)\right]^{1/2}}\left(\begin{matrix}\frac{\epsilon+l\Delta}{\eta\hbar
v_{F}|{\bf k}|}e^{-i\eta\theta}\\\ l\end{matrix}\right),$ (4)

where $|{\bf k}|=\sqrt{k_{x}^{2}+k_{y}^{2}}$,
$\epsilon_{l}=l\epsilon=l\sqrt{\hbar^{2}v_{F}^{2}|{\bf k}|^{2}+\Delta^{2}}$,
$l$ is the conduction ($l=+1$) and valence ($l=-1$) band index, and
$\theta=\tan^{-1}(k_{y}/k_{x})$. The matrix elements $\langle l{\bf
k}|\hat{v}_{\alpha}|l^{\prime}{\bf k}\rangle$ of the velocity operator
$\hat{v}_{\alpha}=\hbar^{-1}\partial\hat{H}^{\eta}_{s}/\partial k_{\alpha}$
are given by

$\displaystyle\langle 1\bm{k}|\hat{v}_{x}|-1\bm{k}\rangle$ $\displaystyle=$
$\displaystyle\langle-1\bm{k}|\hat{v}_{x}|1\bm{k}\rangle^{*}=-v_{F}\frac{\Delta\cos\theta+i\epsilon\sin\theta}{\epsilon},$
(5) $\displaystyle\langle 1\bm{k}|\hat{v}_{x}|1\bm{k}\rangle$ $\displaystyle=$
$\displaystyle-\langle-1\bm{k}|\hat{v}_{x}|-1\bm{k}\rangle=\frac{\hbar
v_{F}^{2}k\cos\theta}{\epsilon},$ (6) $\displaystyle\langle
1\bm{k}|\hat{v}_{y}|-1\bm{k}\rangle$ $\displaystyle=$
$\displaystyle\langle-1\bm{k}|\hat{v}_{y}|1\bm{k}\rangle^{*}=-v_{F}\frac{i\epsilon\cos\theta-\Delta\sin\theta}{\epsilon},$
(7) $\displaystyle\langle 1\bm{k}|\hat{v}_{y}|1\bm{k}\rangle$ $\displaystyle=$
$\displaystyle-\langle-1\bm{k}|\hat{v}_{y}|-1\bm{k}\rangle=\frac{\hbar
v_{F}^{2}k\sin\theta}{\epsilon}.$ (8)

We give a brief derivation for both the linear and third order optical
conductivity terms for a particular Dirac cone. To this end, the equation of
motion $i\hbar\partial\hat{\rho}(t)/\partial
t=\left[{\hat{H}}_{s}^{\eta}+\hat{H}_{i}(t),\hat{\rho}(t)\right]$ for the
density matrix operator is solved perturbatively by expanding the density
matrix in powers of the potential vector amplitude
${\hat{\rho}}(t)=\sum_{n}{\hat{\rho}}^{(n)}(t)$ with
${\hat{\rho}}^{(n)}(t)\propto A^{n}$. To leading order in the potential
vector, the first order density matrix upon Fourier transformation to the
frequency domain can be cast as

$\langle l{\bf k}|\rho^{(1)}(\omega)|l^{\prime}{\bf
k}\rangle=\frac{-ieE_{\alpha}(\omega)}{\omega}\frac{f_{l^{\prime}{\bf
k}}-f_{l{\bf k}}}{\epsilon_{l^{\prime}}-\epsilon_{l}+\hbar\omega}\langle l{\bf
k}|\hat{v}_{\alpha}|l^{\prime}{\bf k}\rangle,$ (9)

where $E_{\alpha}(\omega)=i\omega A_{\alpha}(\omega)/c$. Substituting
$\hat{\rho}^{(1)}(\omega)$ in equation (2) and making use of the identity
$[\hbar\omega\left(\epsilon_{l^{\prime}}-\epsilon_{l}+\hbar\omega\right)]^{-1}=[\epsilon_{l^{\prime}}-\epsilon_{l}]^{-1}[(\hbar\omega)^{-1}-(\epsilon_{l^{\prime}}-\epsilon_{l}+\hbar\omega)^{-1}]$,
we eliminate the first term in the current $\propto Tr(\rho^{(0)})$ and re-
obtain the well know Kubo formula after expressing the current as the product
of conductivity and electric field

$\sigma_{\alpha\beta}^{(1),\eta,s}(\omega)=-4i\hbar^{2}\sigma_{0}\sum\limits_{{\bf
k},l,l^{\prime}}\frac{f_{l^{\prime}{\bf k}}-f_{l{\bf
k}}}{\epsilon_{l^{\prime}}-\epsilon_{l}}\frac{\langle l{\bf
k}|\hat{v}_{\beta}|l^{\prime}{\bf k}\rangle\langle l^{\prime}{\bf
k}|\hat{v}_{\alpha}|l{\bf
k}\rangle}{\epsilon_{l^{\prime}}-\epsilon_{l}+\hbar\omega},$ (10)

where $\sigma_{0}=\alpha c/4$ and the phenomenological relaxation rate is
accounted for by replacing $\omega\rightarrow\omega+i\Gamma$. Equation (10)
accounts for both interband and intraband transitions. The interband
contribution is straightforward and follows directly from (10) by setting
$l^{\prime}\neq l$. The intraband contribution can be derived by equating
enforcing $l^{\prime}=l$, in which case
$(f_{l^{\prime}}-f_{l})/(\epsilon_{l^{\prime}}-\epsilon_{l})\rightarrow\partial
f_{lk}/\partial\epsilon_{l}$. The total charge, spin, and valley linear
conductivities are then computed as $\sigma_{\alpha
x}^{(1)}(\omega)=\sum_{\eta,s}\sigma_{\alpha x}^{(1)\eta,s}(\omega)$,
$\sigma_{\alpha x}^{(1),\textrm{s}}(\omega)=\sum_{\eta,s}s\sigma_{\alpha
x}^{(1)\eta,s}(\omega)/2$, $\sigma_{\alpha
x}^{(1),\eta}(\omega)=\sum_{\eta,s}\eta\sigma_{\alpha x}^{(1)\eta,s}(\omega)$,
and $\sigma_{\alpha x}^{(1),s\eta}(\omega)=\sum_{\eta,s}s\eta\sigma_{\alpha
x}^{(1)\eta,s}(\omega)/2$, respectively.

For the third order conductivities we first solve the equation of motion for
third order correction $\rho^{(3)}(t)$, which is expressed in the Fourier
space as

$\displaystyle\langle l{\bf k}|\rho^{(3)}(\omega)|l^{\prime}{\bf
k}\rangle=\frac{e^{3}E_{\alpha}(\omega_{1})E_{\beta}(\omega_{2})E_{\gamma}(\omega_{3})\\!/\omega_{1}\omega_{2}\omega_{3}}{i\left(\epsilon_{l^{\prime}}-\epsilon_{l}+\hbar\omega_{1}+\hbar\omega_{2}+\hbar\omega_{3}\right)}\\!\\!\\!\sum\limits_{l^{\prime\prime},l^{\prime\prime\prime}}$
$\displaystyle\\!\\!\frac{\langle l{\bf
k}|\hat{v}_{\delta}\\!|l^{\prime\prime\prime}{\bf k}\rangle\\!\langle
l^{\prime\prime\prime}{\bf k}|\hat{v}_{\gamma}\\!|l^{\prime\prime}{\bf
k}\rangle\\!\langle l^{\prime\prime}{\bf k}|\hat{v}_{\beta}\\!|l^{\prime}{\bf
k}\rangle}{\epsilon_{l^{\prime}}-\epsilon_{l^{\prime\prime\prime}}+\hbar\omega_{1}+\hbar\omega_{2}}\\!\\!\left(\\!\\!\frac{f_{l^{\prime}{\bf
k}}-f_{l^{\prime\prime}{\bf
k}}}{\epsilon_{l^{\prime}}\\!-\\!\epsilon_{l^{\prime\prime}}\\!+\\!\hbar\omega_{1}\\!}\\!-\\!\frac{f_{l^{\prime\prime}{\bf
k}}-f_{l^{\prime\prime\prime}{\bf
k}}}{\epsilon_{l^{\prime\prime}}\\!-\\!\epsilon_{l^{\prime\prime\prime}}\\!+\\!\hbar\omega_{2}}\\!\right)$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!-\frac{\langle l{\bf
k}|\hat{v}_{\delta}\\!|l^{\prime\prime\prime}{\bf k}\rangle\\!\langle
l^{\prime\prime\prime}{\bf k}|\hat{v}_{\gamma}\\!|l^{\prime\prime}{\bf
k}\rangle\\!\langle l^{\prime\prime}{\bf k}|\hat{v}_{\beta}\\!|l^{\prime}{\bf
k}\rangle}{\epsilon_{l^{\prime\prime}}-\epsilon_{l}+\hbar\omega_{2}+\hbar\omega_{3}}\\!\\!\left(\\!\\!\frac{f_{l^{\prime\prime}{\bf
k}}-f_{l^{\prime\prime\prime}{\bf
k}}}{\epsilon_{l^{\prime\prime}}\\!-\\!\epsilon_{l^{\prime\prime\prime}}\\!+\\!\hbar\omega_{2}}\\!-\\!\frac{f_{l^{\prime\prime\prime}{\bf
k}}-f_{l{\bf
k}}}{\epsilon_{l^{\prime\prime\prime}}\\!-\\!\epsilon_{l}\\!+\\!\hbar\omega_{3}}\\!\right).$

Substituting the above expression in Eq. (2) and using the same identity as
before with the replacement
$\omega\rightarrow\omega_{1}+\omega_{2}+\omega_{3}$, we obtain the third order
optical conductivity per Dirac cone as

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\sigma_{\alpha\beta\gamma\delta}^{(3)}(\omega_{1},\omega_{2},\omega_{3})$
$\displaystyle=$
$\displaystyle\dfrac{-e\hbar(\omega_{1}+\omega_{2}+\omega_{3})}{E_{\alpha}(\omega_{1})E_{\beta}(\omega_{2})E_{\gamma}(\omega_{3})}\times$
(13) $\displaystyle\times$ $\displaystyle\sum\limits_{{\bf
k},l,l^{\prime}}\frac{\langle l{\bf k}|\hat{v}_{\alpha}|l^{\prime}{\bf
k}\rangle}{\epsilon_{l^{\prime}}-\epsilon_{l}}\langle l^{\prime}{\bf
k}|\rho^{(3)}(\omega)|l{\bf k}\rangle.$ (14)

The phenomenological relaxation rate can be accounted in a similar manner as
in the linear case. Note that, Eq. (13) includes contributions from Berry
connection, as well as intra and interband transitions Sipevelocity ;
Aversa1995 . The nonlinear contributions to the fundamental harmonic
associated to the Kerr effect follow by setting
$\omega_{1}=\omega_{2}=-\omega_{3}$, while the third harmonic conductivity is
can be derived for $\omega_{1}=\omega_{2}=\omega_{3}$. Similarly to the linear
case the total charge, spin, valley, and spin-valley conductivities due to all
Dirac cones are given by $\sigma_{\alpha
x}^{(3)}(n\omega)=\sum_{\eta,s}\sigma_{\alpha x}^{(3)\eta,s}(n\omega)$,
$\sigma_{\alpha x}^{(3),\textrm{s}}(n\omega)=\sum_{\eta,s}s\sigma_{\alpha
x}^{(3)\eta,s}(n\omega)/2$, $\sigma_{\alpha
x}^{(3),\eta}(n\omega)=\sum_{\eta,s}\eta\sigma_{\alpha
x}^{(3)\eta,s}(n\omega)$, and $\sigma_{\alpha
x}^{(3),s\eta}(n\omega)=\sum_{\eta,s}s\eta\sigma_{\alpha
x}^{(3)\eta,s}(n\omega)/2$, respectively, where $n=1,3$.

Data availability

The data that support the plots within this paper and other findings of this
study are available from the corresponding author on reasonable request.

## References

  * [1] A. H. C. Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, A. K. Geim, The electronic properties of graphene. Rev. Mod. Phys. 81, 109 (2009).
  * [2] N. M. R. Peres, Colloquium: The transport properties of graphene: An introduction. Rev. Mod. Phys. 82, 2673 (2010).
  * [3] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, Electric Field Effect in Atomically Thin Carbon Films, Science 306, 666 (2004).
  * [4] M. Liu, X. Yin, E. Ulin-Avila, B. Geng, T. Zentgraf, L. Ju, F. Wang, and X. Zhang, A graphene-based broadband optical modulator, Nature 474, 64 (2011).
  * [5] F. H. L. Koppens, D. E. Chang, and F. J. G. de Abajo, Graphene plasmonics: a platform for strong light-matter interactions, Nano Lett. 11, 3370 (2011).
  * [6] K. S. Noveselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos, and A. A. Firsov, Two-dimensional gas of massless Dirac fermions in graphene, Nature (London) 438, 197 (2005).
  * [7] Y. Zhang, Y. -W. Tan, H. L. Stormer, and P. Kim, Experimental observation of the quantum Hall effect and Berry’s phase in graphene, Nature 438, 201 (2005).
  * [8] N. M. R. Peres, F. Guinea, and A. H. Castro Neto, Electronic properties of disordered two-dimensional carbon. Phys. Rev. B 73, 125411 (2006).
  * [9] S. Konschuh, M. Gmitra, and J. Fabian, Tight-binding theory of the spin-orbit coupling in graphene, Phys. Rev. B 82, 245412 (2010).
  * [10] M. Z. Hasan and C. L. Kane, Colloquium: Topological insulators, Rev. Mod. Phys. 82, 3045 (2010).
  * [11] X.-L. Qi and S.-C. Zhang, opological insulators and superconductors, Rev. Mod. Phys. 83, 1057 (2011).
  * [12] P. Vogt, P. DePadova, C. Quaresima, J. Avila, E. Frantzeskakis, M. C. Asensio, A. Resta, B. Ealet, and G. Le Lay, Silicene: Compelling Experimental Evidence for Graphenelike Two-Dimensional Silicon, Phys. Rev. Lett. 108, 155501 (2012).
  * [13] M. E. Dávila, L. Xian, S. Cahangirov, A. Rubio, and G. L.Lay, Germanene: a novel two-dimensional germanium allotrope akin to graphene and silicene, New J. Phys.16, 095002 (2014).
  * [14] F. -f. Zhu, W. -j. Chen, Y. Xu, C. -l. Gao, D. -d. Guan, C. -h. Liu, D. Qian, S. -C. Zhang, and J. -f. Jia, Epitaxial growth of two-dimensional stanene, Nat. Mater. 14, 1020 (2015).
  * [15] J. Yuhara, B. He, N. Matsunami, M. Nakatake, and G. Le Lay, Graphene’s latest cousin: Plumbene epitaxial growth on a ”Nano WaterCube”, Adv. Mater. 31, 1901017 (2019).
  * [16] A. Castellanos-Gomez, Why all the fuss about 2D semiconductors?, Nat. Photonics 10, 202 (2016).
  * [17] A. Molle, J. Goldberger, M. Houssa, Y. Xu, S.-C. Zhang,and D. Akinwande, Buckled two-dimensional Xene sheets, Nat. Mater. 16, 163 (2017).
  * [18] A. J. Mannix, B. Kiraly, M. C. Hersam, and N. P. Guisinger, Synthesis and chemistry of elemental 2D materials, Nat. Rev. Chem. 1, 0014 (2017).
  * [19] M. Ezawa, Monolayer Topological Insulators: Silicene, Germanene, and Stanene, J. Phys. Soc. Jpn. 84, 121003 (2015).
  * [20] S. A. Mikhailov, Non-linear electromagnetic response of graphene, Europhys. Lett. 79, 27002 (2007);
  * [21] S. A. Mikhailov and K. Ziegler, Nonlinear electromagnetic response of graphene: frequency multiplication and the self-consistent-field effects, J. Phys. Condens. Matter, 20, 384204 (2008).
  * [22] K. F. Mak, C. Lee, J. Hone, J. Shan, and Tony F. Heinz, Atomically Thin MoS2: A New Direct-Gap Semiconductor, Phys. Rev. Lett. 105, 136805 (2010).
  * [23] H. Zhang, S. Virally, Q. L. Bao, L. K. Ping, S. Massar, N. Godbout, and P. Kockaert, Z-scan measurement of the nonlinear refractive index of graphene , Opt. Lett. 37,1856 (2012).
  * [24] E. Hendry, P. J. Hale, J. Moger, A. K. Savchenko, and S. A. Mikhailov, Coherent Nonlinear Optical Response of Graphene, Phys. Rev. Lett. 105, 097401 (2010).
  * [25] J. J. Dean and H. M. van Driel, Graphene and few-layer graphite probed by second-harmonic generation: Theory and experiment, Phys. Rev. B 82, 125411 (2010).
  * [26] R. Wu, Y. Zhang, S. Yan, F. Bian, W. Wang, X. Bai, X. Lu, J. Zhao, and E. Wang, Purely Coherent Nonlinear Optical Response in Solution Dispersions of Graphene Sheets, Nano Lett. 11, 5159 (2011)
  * [27] N. Kumar, J. Kumar, C. Gerstenkorn, R. Wang, H. -Y. Chiu, A. L. Smirl, and H. Zhao, Third harmonic generation in graphene and few-layer graphite films, Phys. Rev. B 87, 121406(R) (2013).
  * [28] S.-Y. Hong, J. I. Dadap, N. Petrone, P.-C. Yeh, J. Hone, and R. M. Osgood, Jr., Optical Third-Harmonic Generation in Graphene, Phys. Rev. X 3, 021014 (2013).
  * [29] J. L. Cheng, N. Vermeulen, and J. E. Sipe, Third order optical nonlinearity of graphene, New J. Phys. 16, 053014 (2014).
  * [30] J. L. Cheng, N. Vermeulen, and J. E. Sipe, Third-order nonlinearity of graphene: Effects of phenomenological relaxation and finite temperature, Phys. Rev. B 91, 235320 (2015).
  * [31] F. Hipolito, T. G. Pedersen, and V. M. Pereira, Nonlinear photocurrents in two-dimensional systems based on graphene and boron nitride, Phys. Rev. B 94, 045434 (2016).
  * [32] S. A. Mikhailov, Quantum theory of the third-order nonlinear electrodynamic effects of graphene, Phys. Rev. B 93, 085403 (2016).
  * [33] A. Singh, K. I. Bolotin, S. Ghosh, and A. Agarwal, Nonlinear optical conductivity of a generic two-band system with application to doped and gapped graphene, Phys. Rev.B 95, 155421(2017).
  * [34] M. Lysne, Y. Murakami, M. Schl̈er, and P. Werner. High-harmonic generation in spin-orbit coupled systems, Phys. Rev. B 102, 081121(R) (2020).
  * [35] O. Zurrón, A. Picón, and L. Plaja. Theory of high-order harmonic generation for gapless graphene, New J. Phys. 20, 053033 (2018).
  * [36] D. Bauer and K. K. Hansen, High-Harmonic Generation in Solids with and without Topological Edge States, Phys. Rev. Lett. 120, 177401 (2018).
  * [37] C. Jürß and D. Bauer, High-harmonic generation in Su-Schrieffer-Heeger chains, Phys. Rev. B 99,195428 (2019).
  * [38] R. Silva, A. Jiménez-Galán, B. Amorim, O. Smirnova and M. Ivanov. Topological strong-field physics on sub-laser-cycle timescale, Nat. Photon. 13, 849 (2019)
  * [39] A. Chacón, W. Zhu, S. P. Kelly, A. Dauphin, E. Pisanty, A. Picón, C. Ticknor, M. F. Ciappina, A. Saxena, and M. Lewenstein, Circular dichroism in higher-order harmonic generation: Heralding topological phases and transitions in Chern insulators, Phys. Rev. B 102, 134115 (2020).
  * [40] D. Baykusheva, A. Chacón, D. Kim, D. E. Kim, D. A. Reis, S. Ghimire, Strong-field physics in three-dimensional topological insulators, arXiv:2008.01265.
  * [41] L. Stille, C. J. Tabert, and E. J. Nicol, Optical signatures of the tunable band gap and valley-spin coupling in silicene, Phys. Rev. B 86, 195405 (2012).
  * [42] M. Ezawa, Photoinduced Topological Phase Transition and a Single Dirac-Cone State in Silicene, Phys. Rev. Lett. 110, 026603 (2013).
  * [43] C. J. Tabert and E. J. Nicol, Dynamical polarization function, plasmons, and screening in silicene and other buckled honeycomb lattices. Phys. Rev. B 89, 195410 (2014).
  * [44] P. Rodriguez-Lopez, W. J. M. Kort-Kamp, D. A. R. Dalvit, and L. M. Woods, Nonlocal optical response in topological phase transitions in the graphene family, Phys. Rev. Materials 2, 014003 (2018).
  * [45] C. J. Tabert and E. J. Nicol, Valley-Spin Polarization in the Magneto-Optical Response of Silicene and Other Similar 2D Crystals, Phys. Rev. Lett. 110, 197402(2013).
  * [46] P. Ledwith, and W. J. M. Kort-Kamp, and D. A. R Dalvit, Topological phase transitions and quantum Hall effect in the graphene family, Phys. Rev. B, 97, 165426 (2018).
  * [47] P. Rodriguez-Lopez, W. J. M. Kort-Kamp, D. Dalvit, and L. M. Woods, Casimir force phase transitions in the graphene family, Nat. Commun.8, 14699 (2017).
  * [48] M. B. Farias, W. J. M. Kort-Kamp, and D. A. R. Dalvit, Quantum friction in two-dimensional topological materials. Phys. Rev. B, 97, 161407(R) (2018).
  * [49] W. J. M. Kort-Kamp, Topological phase transitions in the photonic spin Hall effect, Phys. Rev. Lett. 119,147401 (2017).
  * [50] W. Wu, W. Zhang, S. Chen, X. Ling, W. Shu, H. Luo, S. Wen, and X. Yin, Optics Express 26, 23705 (2018).
  * [51] W. Wu, S. Chen, W. Xu, Z. Liu, R. Lou, L. Shen, H. Luo, S. Wen, and X. Yin, Photonics Research 8, B47 (2020).
  * [52] T. Kitagawa, T. Oka, A. Brataas, L. Fu, and E. Demler, Transport properties of nonequilibrium systems under the application of light: Photoinduced quantum Hall insulators without Landau levels, Phys. Rev. B 84, 235108 (2011).
  * [53] M. Ezawa, Spin valleytronics in silicene: Quantum spin Hall-quantum anomalous Hall insulators and single-valley semimetals, Phys. Rev. B 87, 155415 (2013).
  * [54] X. Li, T. Cao, Q. Niu, J. Shi, and J. Feng, Coupling the valley degree of freedom to antiferromagnetic order, Proc. Natl. Acad. Sci. USA 110, 3738 (2013).
  * [55] Q.-F. Liang, L.-H. Wu, and X. Hu, Electrically tunable topological state in [111] perovskite materials with an antiferromagnetic exchange field, New J. Phys. 15, 063031(2013).
  * [56] J. E. Sipe and E. Ghahramani, Nonlinear optical response of semiconductors in the independent-particle approximation, Phys. Rev. B 48, 11705 (1993).
  * [57] C. Aversa and J. E. Sipe, Nonlinear optical susceptibilities of semiconductors: Results with a length-gauge analysis, Phys. Rev. B 52 , 14636 (1995).
  * [58] A. Taghizadeh, F. Hipolito, and T. G. Pedersen, Linear and nonlinear optical response of crystals using length and velocity gauges: Effect of basis truncation, Phys. Rev. B 96, 195413 (2017).
  * [59] D. C. Hutchings, M. Sheik-Bahae, D. J. Hagan, E. W. Van Stryland, Kramers-Krönig relations in nonlinear optics, Opt. Quant. Electron. 24, 1 (1992).

Acknowledgements

Research presented in this article was supported by the Laboratory Directed
Research and Development program of Los Alamos National Laboratory under
project number 20190574ECR. The authors also thank the Center for Nonlinear
Studies at LANL for financial support under project 20190495CR.

Author contributions

W.K.-K. proposed the study and developed the theory together with R.K.M.
R.K.M. performed the numerical calculations. Both authors discussed the
results and contributed to the final version of the manuscript. W.K.-K.
supervised the entire study.

∗ Corresponding author<EMAIL_ADDRESS>