# Zero-inertia limit: From Particle Swarm Optimization to Consensus Based
Optimization

Cristina Cipriani  Department of Mathematics, Technical University Munich,
Munich, Germany<EMAIL_ADDRESS>,  Hui Huang  Department of
Mathematics and Statistics, University of Calgary, Calgary, Canada
<EMAIL_ADDRESS>and  Jinniao Qiu  Department of Mathematics and
Statistics, University of Calgary, Calgary, Canada<EMAIL_ADDRESS>
###### Abstract.

Recently a continuous description of the particle swarm optimization (PSO)
based on a system of stochastic differential equations was proposed by Grassi
and Pareschi in [grassi2020particle] where the authors formally showed the
link between PSO and the consensus based optimization (CBO) through zero-
inertia limit. This paper is devoted to solving this theoretical open problem
proposed in [grassi2020particle] by providing a rigorous derivation of CBO
from PSO through the limit of zero inertia, and a quantified convergence rate
is obtained as well. The proofs are based on a probabilistic approach by
investigating the weak convergence of the corresponding stochastic
differential equations (SDEs) of Mckean type in the continuous path space and
the results are illustrated with some numerical examples.

The authors would like to thank Prof. Lorenzo Pareschi for suggesting this
problem to them. C. C. acknowledges the partial support of the project ”Online
Firestorms And Resentment Propagation On Social Media: Dynamics,
Predictability and Mitigation” of the TUM Institute for Ethics in Artificial
Intelligence. H. H. is partially supported by the Pacific Institute for the
Mathematical Sciences (PIMS) postdoc fellowship. J. Q. is partially supported
by the National Science and Engineering Research Council of Canada (NSERC) and
by the start-up funds from the University of Calgary.

Keywords: Swarm optimization, consensus based optimization, Laplace’s
principle, tightness.

##  1\. Introduction

Over the last decades, large systems of interacting particles are widely used
in the investigation of complex systems that model collective behaviour (or
swarming), an area that has attracted a great deal of attention; see for
instance [carrillo2010asymptotic, bellomo2013microscale, ha2008particle,
motsch2014heterophilious] and references therein. Such complex systems
frequently appear in modeling phenomena such as biological swarms
[cucker2007emergent], crowd dynamics [bellomo2011modeling], self-assembly of
nanoparticles [holm2006formation], and opinion formation
[motsch2014heterophilious]. In the field of global optimization, similar
particle models are also used in metaheuristics [Aarts:1989:SAB:61990,
Back:1997:HEC:548530, Blum:2003:MCO:937503.937505, Gendreau:2010:HM:1941310],
which provide empirically robust solutions to tackle hard optimization
problems with fast algorithms. Metaheuristics are methods that orchestrate an
interaction between local improvement procedures and global/high level
strategies, and combine random and deterministic decisions, to create a
process capable of escaping from local optima and performing a robust search
of a solution space. In the sequel, we consider the following optimization
problem

$x^{*}\in\mbox{argmin}_{x\in\mathbb{R}^{d}}\mathcal{E}(x)\,,$ (1.1)

where $\mathcal{E}(x):~{}\mathbb{R}^{d}\to\mathbb{R}$ is a given continuous
cost function.

One noble example of metaheuristics is the so-called Particle Swarm
Optimization (PSO), which was initially introduced to model the intelligent
collective behavior of complex biological systems such as flocks of birds or
schools of fish [kennedy1995particle, kennedy1997particle, shi1998modified],
and it is now widely recognized as an efficient method for tackling complex
optimization problems [poli2007particle, lin2008particle]. The PSO method
solves optimization problem (1.1) by considering a group of candidate
solutions, which are represented by particles. Then the algorithm moves those
particles in the search space according to certain mathematical relationships
on the particle position and velocity. Each particle is driven to its best
known local location, which is updated once the particles find better
positions. However the mathematical understanding of PSO is still in its
infancy. Recently Grassi and Pareschi [grassi2020particle] took a significant
first step towards a mathematical theory for PSO based on a continuous
description in the form of a system of stochastic differential equations:

$\displaystyle\begin{cases}dX_{t}^{i,m}=V_{t}^{i,m}dt,\\\
dV_{t}^{i,m}=-\frac{\gamma}{m}V_{t}^{i,m}dt+\frac{\lambda}{m}(X_{t}^{\alpha}(\rho^{N,m})-X_{t}^{i,m})dt+\frac{\sigma}{m}D(X_{t}^{\alpha}(\rho^{N,m})-X_{t}^{i,m})dB_{t}^{i},\quad
i=1,\cdots,N\,,\end{cases}$ (1.2)

where the $\mathbb{R}^{d}$-valued functions $X_{t}^{i,m}$ and $V_{t}^{i,m}$
denote the position and velocity of the $i$-th particle at time $t$, $m>0$ is
the inertia weight, $\gamma=1-m\geq 0$ is the friction coefficient,
$\lambda>0$ is the acceleration coefficient, $\sigma>0$ is the diffusion
coefficient, and $\\{(B_{t}^{i})_{t\geq 0}\\}_{i=1}^{N}$ are $N$ independent
$d$-dimensional Brownian motions. We also use the notations for the diagonal
matrix

$D(X_{t}):=\mbox{diag}\\{(X_{t})_{1},\dots,(X_{t})_{d}\\}\in\mathbb{R}^{d\times
d}\,,$

where $(X_{t})_{k}$ is the $k$-th component of $X_{t}$, and the weighted
average is given by

${X}_{t}^{\alpha}(\rho^{N,m}):=\frac{\int_{\mathbb{R}^{d}}x\omega_{\alpha}^{\mathcal{E}}(x)\rho^{N,m}(t,dx)}{\int_{\mathbb{R}^{d}}\omega_{\alpha}^{\mathcal{E}}(x)\rho^{N,m}(t,dx)}$
(1.3)

with the empirical measure
$\rho^{N,m}(t,dx):=\frac{1}{N}\sum_{i=1}^{N}\delta_{X_{t}^{i,m}}(dx)$. So we
can rewrite

$D(X_{t}^{\alpha}(\rho^{N,m})-X_{t}^{i,m})dB_{t}^{i}=\sum_{k=1}^{d}(X_{t}^{\alpha}(\rho^{N,m})-X_{t}^{i,m})_{k}d(B_{t}^{i})^{k}e_{k}\,,$
(1.4)

where $e_{k}$ is the unit vector in the $k$-th dimension for $k=1,\dots,d$.
Furthermore, the initial data $(X_{0}^{i},V_{0}^{i})_{i=1}^{N}$ are
independent and identically distributed (i.i.d.) with the common distribution
$f_{0}\in\mathcal{P}_{4}(\mathbb{R}^{2d})$, where
$\mathcal{P}_{4}(\mathbb{R}^{2d})$ denotes the space of probability measures
with finite forth moment, endowed with the Wasserstein distance
[ambrosio2008gradient]. The choice of the weight function

$\omega_{\alpha}^{\mathcal{E}}(x):=\exp(-\alpha\mathcal{E}(x))$

comes from the well-known Laplace’s principle [miller2006applied, Dembo2010],
a classical asymptotic method for integrals, which states that for any
probability measure $\rho\in\mathcal{P}(\mathbb{R}^{d})$, there holds

$\lim\limits_{\alpha\to\infty}\left(-\frac{1}{\alpha}\log\left(\int_{\mathbb{R}^{d}}\omega_{\alpha}^{\mathcal{E}}(x)\rho(dx)\right)\right)=\inf\limits_{x\in\rm{supp}(\rho)}\mathcal{E}(x)\,.$
(1.5)

Thus for $\alpha$ large enough, one expects that

${X}_{t}^{\alpha}(\rho^{N,m})\approx\mbox{argmin
}\\{\mathcal{E}(X_{t}^{1,m}),\dots,\mathcal{E}(X_{t}^{N,m})\\}\,,$

which means that ${X}_{t}^{\alpha}(\rho^{N,m})$ is a global best location at
time $t$.

Figure 1. Left: the Ackley function for $d=2$ with the unique global minimum
at the point $x^{*}=(0,0)^{T}$. Right: Particles trajectories of the PSO model
(1.2) along the simulation for the 2-$d$ Ackley function with the global
minimizer $x^{*}$. The simulation parameters are: time discretization $0.01$,
number of particles $10^{3}$, $\lambda=1$, $\sigma=\frac{1}{\sqrt{3}}$,
$\alpha=30$, $m=0.1$. The initial data are sampled from a normal bi-
dimensional distribution. Figure 2. Application of the PSO dynamics (1.2) to
the 2-$d$ Ackley function $\mathcal{E}(x)$ with the global minimizer
$x^{*}=(0,0)^{T}$. Particles initially have a normal distribution around
$x^{*}$. Then all particles converge to one point, the global minimizer
$x^{*}$, and they stop moving eventually, i.e. velocity converges to zero. The
simulation parameters are the ones described below Figure 1.

Before starting our analysis of the PSO dynamics (1.2), let us illustrate
numerically the behavior of the dynamics for the benchmark Ackley function

$\mathcal{E}(x)=-20\exp\left(-\frac{0.2}{\sqrt{d}}\left|x-x^{*}\right|\right)-\exp\left(\frac{1}{d}\sum_{k=1}^{d}\cos\left(2\pi\left(x_{k}-x_{k}^{*}\right)\right)\right)+e+20$

in the case of $d=2$, and with the global minimizer $x^{*}=(0,0)^{T}$. In
Figure 1 and 2, we initialize the particles with a normal distribution around
$x^{*}$ and then apply a discretization scheme (which will be explained in
Section 4) to the system (1.2). We can see that all the particles successfully
find the global minimizer $x^{*}$, and particles’ velocity converges to zero.

As it has been shown in [grassi2020particle], in the zero-inertia limit ($m\to
0$), one may expect to obtain the recent developed Consensus Based
Optimization (CBO) dynamics [carrillo2018analytical, ha2019convergence, PTTM,
fornasier2021consensus] satisfying

$dX_{t}^{i}=\lambda(X_{t}^{\alpha}(\rho^{N})-X_{t}^{i})dt+\sigma
D(X_{t}^{\alpha}(\rho^{N})-X_{t}^{i})dB_{t}^{i},\quad i=1,\cdots,N\,,$ (1.6)

where

${X}_{t}^{\alpha}(\rho^{N}):=\frac{\int_{\mathbb{R}^{d}}x\omega_{\alpha}^{\mathcal{E}}(x)\rho^{N}(t,dx)}{\int_{\mathbb{R}^{d}}\omega_{\alpha}^{\mathcal{E}}(x)\rho^{N}(t,dx)}\mbox{
with }\rho^{N}(t,dx):=\frac{1}{N}\sum_{i=1}^{N}\delta_{X_{t}^{i}}(dx)\,.$

It has been proved that CBO is a powerful and robust method to solve many
interesting non-convex high-dimensional optimization problems in machine
learning [CJ]. By now, CBO methods have also been generalized to optimization
over manifolds [FHPS1, FHPS2, kim2020stochastic, fornasier2021anisotropic].
The objective of the present paper is to complete a theory gap suggested in
[grassi2020particle] by providing a rigorous proof of the zero-inertia limit.

On the one hand, as $N\to\infty$, the mean-field limit results (see
[bolley2011stochastic, huang2020mean, carrillo2019propagation,
sznitman1991topics, jabin2017mean, huang2021note] for instance) indicate that
our PSO dynamics (1.2) converge to the solutions of following mean-field
nonlinear Mckean systems:

$\displaystyle d\overline{X}_{t}^{m}=\overline{V}_{t}^{m}dt,$ (1.7a)
$\displaystyle
d\overline{V}_{t}^{m}=-\frac{\gamma}{m}\overline{V}_{t}^{m}dt+\frac{\lambda}{m}(X_{t}^{\alpha}(\rho^{m})-\overline{X}_{t}^{m})dt+\frac{\sigma}{m}D(X_{t}^{\alpha}(\rho^{m})-\overline{X}_{t}^{m})dB_{t}\,,$
(1.7b)

where

$X_{t}^{\alpha}(\rho^{m})=\frac{\int_{\mathbb{R}^{d}}x\omega_{\alpha}^{\mathcal{E}}(x)\rho^{m}(t,dx)}{\int_{\mathbb{R}^{d}}\omega_{\alpha}^{\mathcal{E}}(x)\rho^{m}(t,dx)},\quad\rho^{m}(t,x)=\int_{\mathbb{R}^{d}}f^{m}(t,x,dv)\,,$
(1.8)

and the initial data $(\overline{X}_{0},\overline{V}_{0})$ is the same as in
(1.2). Here $f^{m}(t,x,v)$ is the distribution of
$(\overline{X}_{t}^{m},\overline{V}_{t}^{m})$ at time $t$ , which makes the
set of equations (1.7) nonlinear. We refer to [huang2021note] for a proof the
well-posedness of PSO particle system (1.2) and its mean-field dynamic (1.7).
A direct application of the Itô-Doeblin formula yields that the law
$f_{t}^{m}:=f^{m}(t,\cdot,\cdot)$ at time $t$ is a weak solution to the
following nonlinear Vlasov-Fokker-Plank equation

$\partial_{t}f^{m}_{t}+v\cdot\nabla_{x}f^{m}_{t}=\nabla_{v}\cdot\left(\frac{\gamma}{m}vf^{m}_{t}+\frac{\lambda}{m}\left(x-{X}_{t}^{\alpha}(\rho^{m})\right)f+\frac{\sigma^{2}}{2m^{2}}D\left(x-{X}_{t}^{\alpha}(\rho^{m})\right)^{2}\nabla_{v}f^{m}_{t}\right),$
(1.9)

with the initial data
$f^{m}_{0}(x,v)=\mbox{Law}(\overline{X}_{0},\overline{V}_{0})$. On the other
hand, taking $N\to\infty$ in (1.6) leads to the mean-field CBO dynamic of the
form

$\displaystyle
d\overline{X}_{t}=\lambda(X_{t}^{\alpha}(\rho)-\overline{X}_{t})dt+\sigma
D(X_{t}^{\alpha}(\rho)-\overline{X}_{t})dB_{t}$ (1.10)

with $\rho_{t}=\mbox{Law}(\overline{X}_{t})$ satisfying the corresponding CBO
equation

$\partial_{t}\rho_{t}+\lambda\nabla_{x}\cdot(\rho_{t}({X}_{t}^{\alpha}(\rho)-x))=\frac{\sigma^{2}}{2}\sum_{j=1}^{d}\frac{\partial^{2}}{\partial
x_{j}^{2}}\left(\rho_{t}\left(x_{j}-({X}_{t}^{\alpha}(\rho))_{j}\right)^{2}\right)\,.$
(1.11)

In this paper, we prove that in the zero-inertia limit, as $m\to 0^{+}$, the
processes $\\{\overline{X}^{m}\\}$ satisfying SDEs (1.7) converge weakly to
the solution $\overline{X}$ to SDE (1.10) in the continuous path space
$\mathcal{C}([0,T];\mathbb{R}^{d})$. A convergence rate is obtained and the
generalizations to cases with memory effects are also addressed. This is
related to the study of the overdamped limit [kramers1940brownian,
choi2020quantified, duong2017variational], or large friction limit
[carrillo2020quantitative, jabin2000macroscopic, fetecau2015first] for Vlasov
type equations. However, the nolinear term $X_{t}^{\alpha}(\rho^{m})$ here
makes our model very different from theirs, which is nonstandard in the
literature. Moreover all of those results mentioned earlier are obtained
through the investigation of PDEs like (1.9) and (1.11), while in the present
paper we adopt a probabilistic approach by investigating the weak convergence
of the non-Markovian stochastic processes $\\{\overline{X}^{m}\\}$ satisfying
SDE (1.7) to the solution $\\{\overline{X}\\}$ to SDE (1.10)) in the
continuous path space.

The rest of the paper is organized as follows: In Section 2 we verify the
tightness of the PSO model (1.7) through Aldous criteria, which allows us to
obtain the zero-inertia limit from the PSO model (1.7) towards the CBO model
(1.10) as $m\to 0$; see Theorem 2.3. Then in Section 3 we generalize the
result to the PSO model with memory effects of the local best positions.
Lastly we conclude this paper in Section 4 by reporting a few instructive
numerical experiments on validating the zero-inertia limit.

##  2\. Zero-inertia limit

Throughout this work, the letter $C$ denotes a generic constant whose value
may vary from line to line and its dependence on certain model parameters will
be specified whenever needed. We start this section with the standing
assumption on the cost function $\mathcal{E}$.

######  Assumption 1.

The given cost function $\mathcal{E}:\mathbb{R}^{d}\rightarrow\mathbb{R}$ is
locally Lipschitz continuous and satisfies the properties:

  * 1.

There exists some constant $L>0$ such $|\mathcal{E}(x)-\mathcal{E}(y)|\leq
L(|x|+|y|)|x-y|$ for all $x,y\in\mathbb{R}^{d}$;

  * 2.

$\mathcal{E}$ is uniformly bounded, i.e.
$-\infty<\underline{\mathcal{E}}:=\inf\mathcal{E}\leq\mathcal{E}\leq\sup\mathcal{E}=:\overline{\mathcal{E}}<+\infty$,
and define
$C_{\alpha,\mathcal{E}}:=e^{\alpha(\overline{\mathcal{E}}-\underline{\mathcal{E}})}$
.

The following theorem gives the well-posedness of the mean-field PSO and CBO
dynamics (1.7) and (1.10) whose proofs are analogous to [huang2021note,
Theorem 2.3] and [carrillo2018analytical, Theorem 3.1], and thus omitted.

######  Theorem 2.1.

Let Assumption 1 hold. For each $T>0$, there hold the following assertions.

(i) If
$(\overline{X}^{m}_{0},\overline{V}^{m}_{0})=(\overline{X}_{0},\overline{V}_{0})$
is distributed according to $f_{0}$ with
$f_{0}\in\mathcal{P}_{4}(\mathbb{R}^{2d})$, then for each $m\in(0,1]$, the
nonlinear SDE (1.7) admits a unique solution up to time $T$ with the initial
data $(\overline{X}^{m}_{0},\overline{V}^{m}_{0})$ and it holds further that

$\sup\limits_{t\in[0,T]}\mathbb{E}\left[|\overline{X}_{t}^{m}|^{4}+|\overline{V}_{t}^{m}|^{4}\right]\leq
e^{CT}\cdot\mathbb{E}\left[|\overline{X}_{0}|^{4}+|\overline{V}_{0}|^{4}\right]\,,$
(2.1)

where $C$ depends only on $\lambda,m,\sigma$, and $C_{\alpha,\mathcal{E}}$.

(ii) If $\overline{X}_{0}$ is distributed according to
$\rho_{0}\in\mathcal{P}_{4}(\mathbb{R}^{d})$, then SDE (1.10) admits a unique
solution up to time $T$ with the initial data $\overline{X}_{0}$ and it holds
further that

$\sup\limits_{t\in[0,T]}\mathbb{E}\left[|\overline{X}_{t}|^{4}\right]\leq
e^{CT}\cdot\mathbb{E}\left[|\overline{X}_{0}|^{4}\right]\,,$ (2.2)

where $C$ depends only on $\lambda,\sigma$ and $C_{\alpha,\mathcal{E}}$.

Solving (1.7b) for $\overline{V}_{t}^{m}$ gives

$\overline{V}_{t}^{m}=e^{-\frac{\gamma}{m}t}\overline{V}_{0}+\frac{\lambda}{m}\int_{0}^{t}e^{-\frac{\gamma}{m}(t-s)}(X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m})ds+\frac{\sigma}{m}\int_{0}^{t}e^{-\frac{\gamma}{m}(t-s)}D(X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m})dB_{s}\,,$

which implies that

$\displaystyle\overline{X}_{t}^{m}$
$\displaystyle=\overline{X}_{0}+\int_{0}^{t}\overline{V}_{\tau}d\tau=\overline{X}_{0}+\int_{0}^{t}e^{-\frac{\gamma}{m}\tau}\overline{V}_{0}d\tau+\frac{\lambda}{m}\int_{0}^{t}\int_{0}^{\tau}e^{-\frac{\gamma}{m}(\tau-s)}(X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m})dsd\tau$
$\displaystyle\quad+\frac{\sigma}{m}\int_{0}^{t}\int_{0}^{\tau}e^{-\frac{\gamma}{m}(\tau-s)}D(X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m})dB_{s}d\tau\,.$
(2.3)

Then $\overline{X}_{t}^{m}$ has the law $\rho^{m}_{t}$ for each $t\geq 0$.
Denote by $\mathcal{C}([0,T];\mathbb{R}^{d})$ the space of all
$\mathbb{R}^{d}$-valued continuous functions on $[0,T]$ equipped with the
usual uniform norm $\|\cdot\|_{0}$. Each continuous stochastic process
$\overline{X}^{m}$ may be seen as a $\mathcal{C}([0,T];\mathbb{R}^{d})$-valued
random function and it induces a probability measure (or law, denoted by
$\rho^{m}$) on $\mathcal{C}([0,T];\mathbb{R}^{d})$. We shall use the weak
convergence in the space of probability measures on
$\mathcal{C}([0,T];\mathbb{R}^{d})$. In what follows, we write
$\overline{X}^{m}\rightharpoonup\overline{X}$ or $\rho^{m}\rightharpoonup\rho$
with $\rho$ being the law of $\overline{X}$, if
$\left\\{\rho^{m}\right\\}_{m>0}$, as a sequence of probability measures on
$\mathcal{C}([0,T];\mathbb{R}^{d})$, converges weakly to $\rho$, i.e., for
each bounded continuous functional $\Phi$ on
$\mathcal{C}([0,T];\mathbb{R}^{d})$, there holds $\lim_{m\rightarrow
0^{+}}\mathbb{E}\left[\Phi(\overline{X}^{m})\right]=\mathbb{E}\left[\Phi(\overline{X})\right]$.
The weak convergence $\overline{X}^{m}\rightharpoonup\overline{X}$ is stronger
than and actually implies the convergence of $\\{\rho^{m}_{t}\\}_{m>0}$ to
$\rho_{t}$ with $\rho_{t}$ being the law of $\overline{X}_{t}$ for each $t\geq
0$, while the converse need not hold. Moreover, due to the separability and
completeness of the space $\mathcal{C}([0,T];\mathbb{R}^{d})$, Prohorov’s
theorem implies that the relative compactness is equivalent to the tightness;
see [billingsley1999convergence] for more details.

The proof of zero-inertia limit will proceed in two steps:

  * •

The tightness of the sequence of probability distributions
$\\{\rho^{m}\\}_{0<m\leq 1}$ of $\\{\overline{X}^{m}\\}_{0<m\leq 1}$ is
justified by using Aldous tightness criteria.

  * •

We will check that all the limit points of $\\{\overline{X}^{m}\\}_{0<m\leq
1}$ as $m\to 0$ satisfy mean-field CBO dynamic (1.10) which in fact admits a
unique solution.

For the sake of completeness, let us recall Aldous tightness criteria.

######  Lemma 2.1.

Let $\\{X^{n}\\}_{n\in\mathbb{N}}$ be a sequence of random variables defined
on a probability space $(\Omega,\mathcal{F},\mathbb{P})$ and valued in
$\mathcal{C}([0,T];\mathbb{R}^{d})$. The sequence of probability distributions
$\\{\mu_{X^{n}}\\}_{n\in\mathbb{N}}$ of $\\{X^{n}\\}_{n\in\mathbb{N}}$ is
tight on $\mathcal{C}([0,T];\mathbb{R}^{d})$ if the following two conditions
hold.

$(Con1)$ For all $t\geq 0$, the set of distributions of $X_{t}^{n}$, denoted
by $\\{\mu_{X_{t}^{n}}\\}_{n\in\mathbb{N}}$, is tight in $\mathbb{R}^{d}$.

$(Con2)$ For all $\varepsilon>0$, $\eta>0$, there exists $\delta_{0}>0$ and
$n_{0}\in\mathbb{N}$ such that for all $n\geq n_{0}$ and for all discrete-
valued $\sigma(X^{n}_{s};s\in[0,T])$-stopping times $\beta$ such that
$0\leq\beta+\delta_{0}\leq T$,

$\sup_{\delta\in[0,\delta_{0}]}\mathbb{P}\left(|X^{n}_{\beta+\delta}-X^{n}_{\beta}|\geq\eta\right)\leq\varepsilon\,.$
(2.4)

######  Theorem 2.2 (Tightness).

Let Assumption 1 hold and $(X_{t}^{m},V_{t}^{m})_{t\in[0,T]}$ satisfy the
system (1.7). For each countable subsequence
$\\{m_{k}\\}_{k\in\mathbb{N}}\subset[0,1]$ with
$\lim_{k\rightarrow\infty}m_{k}=0$, the sequence of probability distributions
$\\{\rho^{m_{k}}\\}_{k\in\mathbb{N}}$ of
$\\{\overline{X}^{m_{k}}\\}_{k\in\mathbb{N}}$ is tight.

###### Proof.

By Lemma 2.1, it is sufficient to justify conditions $(Con1)$ and $(Con2)$ in
Aldous tightness criteria .

$\bullet$ Step 1: Checking $(Con1)$.  First, for $0<m\leq\frac{1}{2}$,
recalling (2), we have by Fubini’s theorem (see [da2014stochastic, Theorem
4.33] for the stochastic version)

$\displaystyle\overline{X}_{t}^{m}$
$\displaystyle=\overline{X}_{0}+\int_{0}^{t}e^{-\frac{\gamma}{m}\tau}\overline{V}_{0}d\tau+\frac{\lambda}{m}\int_{0}^{t}\int_{0}^{\tau}e^{-\frac{\gamma}{m}(\tau-s)}(X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m})dsd\tau$
$\displaystyle\quad+\frac{\sigma}{m}\int_{0}^{t}\int_{0}^{\tau}e^{-\frac{\gamma}{m}(\tau-s)}D(X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m})dB_{s}d\tau$
$\displaystyle=\overline{X}_{0}+\int_{0}^{t}e^{-\frac{\gamma}{m}\tau}\overline{V}_{0}d\tau+\frac{\lambda}{m}\int_{0}^{t}\int_{s}^{t}e^{-\frac{\gamma}{m}(\tau-s)}d\tau(X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m})ds$
$\displaystyle\quad+\frac{\sigma}{m}\int_{0}^{t}\int_{s}^{t}e^{-\frac{\gamma}{m}(\tau-s)}d\tau
D(X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m})dB_{s}$
$\displaystyle=\overline{X}_{0}+\frac{m}{\gamma}(1-e^{-\frac{\gamma}{m}t})\overline{V}_{0}+\frac{\lambda}{\gamma}\int_{0}^{t}(1-e^{-\frac{\gamma}{m}(t-s)})(X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m})ds$
$\displaystyle\quad+\frac{\sigma}{\gamma}\int_{0}^{t}(1-e^{-\frac{\gamma}{m}(t-s)})D(X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m})dB_{s}\,.$
(2.5)

Here the assumption on $0<m\leq\frac{1}{2}$ ensures that
$\gamma=1-m\in[\frac{1}{2},1)$, so $\frac{1}{\gamma}$ is well defined. It
follows from Hölder’s inequality that

$|\overline{X}_{t}^{m}|^{4}\leq
64|\overline{X}_{0}|^{4}+\frac{64m^{4}}{\gamma^{4}}|\overline{V}_{0}|^{4}+\frac{64\lambda^{4}t^{3}}{\gamma^{4}}\int_{0}^{t}|X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m}|^{4}ds+\frac{64\sigma^{4}}{\gamma^{4}}\left|\int_{0}^{t}(1-e^{-\frac{\gamma}{m}(t-s)})D(X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m})dB_{s}\right|^{4},$

where we have used the fact that for any sequence $\\{a_{i}\\}_{i=1}^{n}\geq
0$ and $p\geq 2$, there holds

$\left(\sum_{i=1}^{n}a_{i}\right)^{p}\leq n^{p-1}\sum_{i=1}^{n}a_{i}^{p}\,.$

Using the moment inequality for stochastic integrals as in [mao2007stochastic,
Theorem 7.1] yields that

$\displaystyle\mathbb{E}\left[\left|\int_{0}^{t}(1-e^{-\frac{\gamma}{m}(t-s)})D(X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m})dB_{s}\right|^{4}\right]$
$\displaystyle\leq$ $\displaystyle
d^{3}\mathbb{E}\left[\sum_{k=1}^{d}\left|\int_{0}^{t}(1-e^{-\frac{\gamma}{m}(t-s)})(X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m})_{k}dB_{s}^{k}e_{k}\right|^{4}\right]$
$\displaystyle\leq$ $\displaystyle
36d^{3}t\int_{0}^{t}\mathbb{E}\left[\sum_{k=1}^{d}|(X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m})_{k}|^{4}\right]ds\leq
36d^{3}t\int_{0}^{t}\mathbb{E}\left[|X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m}|^{4}\right]ds\,.$

Thus,

$\displaystyle\mathbb{E}[|\overline{X}_{t}^{m}|^{4}]\leq
64\mathbb{E}[|\overline{X}_{0}|^{4}]+\frac{64m^{4}}{\gamma^{4}}\mathbb{E}[|\overline{V}_{0}|^{4}]+\frac{64(\lambda^{4}t^{3}+36d^{3}t\sigma^{4})}{\gamma^{4}}\int_{0}^{t}\mathbb{E}[|X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m}|^{4}]ds\,.$

Notice that

$\displaystyle\mathbb{E}[|X_{t}^{\alpha}(\rho^{m})-\overline{X}_{t}^{m}|^{4}]$
$\displaystyle=\int_{\mathbb{R}^{d}}\left|\frac{\int_{\mathbb{R}^{d}}x\omega_{\alpha}^{\mathcal{E}}(x)\rho^{m}(t,dx)}{\int_{\mathbb{R}^{d}}\omega_{\alpha}^{\mathcal{E}}(x)\rho^{m}(t,dx)}-y\right|^{4}\rho^{m}(t,dy)=\int_{\mathbb{R}^{d}}\left|\frac{\int_{\mathbb{R}^{d}}(x-y)\omega_{\alpha}^{\mathcal{E}}(x)\rho^{m}(t,dx)}{\int_{\mathbb{R}^{d}}\omega_{\alpha}^{\mathcal{E}}(x)\rho^{m}(t,dx)}\right|^{4}\rho^{m}(t,dy)$
$\displaystyle\leq\frac{\int_{\mathbb{R}^{d}}\int_{\mathbb{R}^{d}}|x-y|^{4}\omega_{\alpha}^{\mathcal{E}}(x)\rho^{m}(t,dx)\rho^{m}(t,dy)}{\int_{\mathbb{R}^{d}}\omega_{\alpha}^{\mathcal{E}}(x)\rho^{m}(t,dx)}\leq
16C_{\alpha,\mathcal{E}}\mathbb{E}[|\overline{X}_{t}^{m}|^{4}]\,,$ (2.6)

where Jensen’s inequality is applied in the first inequality. Thus we have

$\displaystyle\mathbb{E}[|\overline{X}_{t}^{m}|^{4}]\leq
64\mathbb{E}[|\overline{X}_{0}|^{4}]+\frac{64m^{4}}{\gamma^{4}}\mathbb{E}[|\overline{V}_{0}|^{4}]+\frac{1024C_{\alpha,\mathcal{E}}(\lambda^{4}t^{3}+36d^{3}t\sigma^{4})}{\gamma^{4}}\int_{0}^{t}\mathbb{E}[|\overline{X}_{s}^{m}|^{4}]ds\,.$

Using Gronwall’s inequality leads to

$\mathbb{E}[|\overline{X}_{t}^{m}|^{4}]\leq\left(64\mathbb{E}[|\overline{X}_{0}|^{4}]+\frac{64m^{4}}{\gamma^{4}}\mathbb{E}[|\overline{V}_{0}|^{4}]\right)\exp\left(\frac{1024C_{\alpha,\mathcal{E}}(\lambda^{4}T^{3}+36d^{3}T\sigma^{4})}{\gamma^{4}}T\right),\quad
t\in[0,T]\,.$ (2.7)

Recalling $0\leq m\leq\frac{1}{2}$ and $\frac{1}{\gamma}=\frac{1}{1-m}\leq 2$,
from estimate (2.7) we obtain the boundedness:

$\mathbb{E}[|\overline{X}_{t}^{m}|^{4}]\leq
C(\mathbb{E}[|\overline{X}_{0}|^{4}],\mathbb{E}[|\overline{V}_{0}|^{4}],C_{\alpha,\mathcal{E}},\lambda,d,\sigma,T)\,.$
(2.8)

Next we consider the case when $\frac{1}{2}\leq m\leq 1$. It is obvious that

$\displaystyle|\overline{X}_{t}^{m}|^{4}=|\overline{X}_{0}|^{4}+4\int_{0}^{t}|\overline{X}_{s}^{m}|^{2}\overline{X}_{s}^{m}\cdot\overline{V}_{s}^{m}ds$
$\displaystyle\leq|\overline{X}_{0}^{m}|^{4}+8\int_{0}^{t}|\overline{X}_{s}^{m}|^{2}(|\overline{X}_{s}^{m}|^{2}+|\overline{V}_{s}^{m}|^{2})ds$
$\displaystyle\leq|\overline{X}_{0}^{m}|^{4}+C\int_{0}^{t}(|\overline{X}_{s}^{m}|^{4}+|\overline{V}_{s}^{m}|^{4})ds\,.$
(2.9)

Applying Itô-Doeblin formula to (1.7b) gives

$\displaystyle|\overline{V}_{t}^{m}|^{4}$
$\displaystyle=|\overline{V}_{0}|^{4}+\frac{4\lambda}{m}\int_{0}^{t}|\overline{V}_{s}^{m}|^{2}\overline{V}_{s}^{m}\cdot(X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m})ds+\frac{4\sigma}{m}\int_{0}^{t}|\overline{V}_{s}^{m}|^{2}\overline{V}_{s}^{m}\cdot
D(X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m})dB_{s}$
$\displaystyle\quad+\frac{6\sigma^{2}}{m^{2}}\int_{0}^{t}|\overline{V}_{s}^{m}|^{2}|X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m}|^{2}ds-\int_{0}^{t}\frac{4\gamma}{m}|\overline{V}_{s}^{m}|^{4}ds$
$\displaystyle\leq|\overline{V}_{0}|^{4}+C\left(\frac{\lambda}{m}+\frac{\sigma^{2}}{m^{2}}\right)\int_{0}^{t}|\overline{V}_{s}^{m}|^{4}ds+C\left(\frac{\lambda}{m}+\frac{\sigma^{2}}{m^{2}}\right)\int_{0}^{t}|X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m}|^{4}ds$
$\displaystyle\quad+\frac{4\sigma}{m}\int_{0}^{t}|\overline{V}_{s}^{m}|^{2}\overline{V}_{s}^{m}\cdot
D(X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m})dB_{s}\,.$ (2.10)

Collecting estimates (2) and (2.9) and recalling $\frac{1}{2}\leq m\leq 1$, we
have

$\displaystyle\mathbb{E}[|\overline{X}_{t}^{m}|^{4}+|\overline{V}_{t}^{m}|^{4}]$
$\displaystyle\leq$
$\displaystyle\mathbb{E}[|\overline{X}_{0}|^{4}+|\overline{V}_{0}|^{4}]+C\int_{0}^{t}\mathbb{E}[|\overline{X}_{s}^{m}|^{4}+|\overline{V}_{s}^{m}|^{4}]ds+C\int_{0}^{t}\mathbb{E}[|X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m}|^{4}]ds$
$\displaystyle\leq$
$\displaystyle\mathbb{E}[|\overline{X}_{0}|^{4}+|\overline{V}_{0}|^{4}]+C(1+8C_{\alpha,\mathcal{E}})\int_{0}^{t}\mathbb{E}[|\overline{X}_{s}^{m}|^{4}+|\overline{V}_{s}^{m}|^{4}]ds\,,$
(2.11)

where the estimate (2) is used in the last inequality. Applying Gronwall’s
inequality yields that

$\displaystyle\mathbb{E}[|\overline{X}_{t}^{m}|^{4}+|\overline{V}_{t}^{m}|^{4}]\leq\mathbb{E}[|\overline{X}_{0}|^{4}+|\overline{V}_{0}|^{4}]\exp\left(C\cdot(1+8C_{\alpha,\mathcal{E}})t\right),\quad
t\in[0,T]\,.$ (2.12)

Finally, combining (2.8) and (2.12) yields that

$\sup_{m\in(0,1]}\sup_{t\in[0,T]}\mathbb{E}[|\overline{X}_{t}^{m}|^{4}]\leq
C(\mathbb{E}[|\overline{X}_{0}|^{4}],\mathbb{E}[|\overline{V}_{0}|^{4}],C_{\alpha,\mathcal{E}},\lambda,\sigma,d,T)=:C_{1}$
(2.13)

where the constant $C_{1}>0$ is independent of $m$. Therefore, for any
$\varepsilon>0$, there exists a compact subset
$K_{\varepsilon}:=\\{x:~{}|x|^{4}\leq\frac{C_{1}}{\varepsilon}\\}$ such that
by Markov’s inequality

$\rho_{t}^{m}((K_{\varepsilon})^{c})=\mathbb{P}(|X_{t}^{m}|^{4}>\frac{C_{1}}{\varepsilon})\leq\frac{\varepsilon\mathbb{E}[|X_{t}^{m}|^{4}]}{C_{1}}\leq\varepsilon,\quad\forall~{}0<m\leq
1\,.$ (2.14)

This means that for each $t\in[0,T]$, each countable subset of
$\\{\rho_{t}^{m}\\}_{0<m\leq 1}$ is tight, which verifies condition $(Con1)$
in Lemma 2.1.

$\bullet$ Step 2: Checking $(Con2)$.  Let $\beta$ be a
$\sigma(X^{m}_{s};s\in[0,T])$-stopping time with discrete values such that
$\beta+\delta_{0}\leq T$. Without any loss of generality, we may assume that
the concerned countable subsequence $\\{m_{k}\\}_{k\in\mathbb{N}}\subset[0,1]$
satisfies $m_{k}\leq\frac{1}{2}$ for all $k\in\mathbb{N}$; thus, we may just
consider the case of $0<m\leq\frac{1}{2}$ which indicates
$\frac{1}{2}\leq\gamma<1$. Recall (2) and compute

$\displaystyle\overline{X}_{\beta+\delta}^{m}-\overline{X}_{\beta}^{m}$
$\displaystyle=\int_{\beta}^{\beta+\delta}\overline{V}_{\tau}d\tau=\int_{\beta}^{\beta+\delta}e^{-\frac{\gamma}{m}\tau}\overline{V}_{0}d\tau+\frac{\lambda}{m}\int_{\beta}^{\beta+\delta}\int_{0}^{\tau}e^{-\frac{\gamma}{m}(\tau-s)}(X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m})dsd\tau$
$\displaystyle\quad+\frac{\sigma}{m}\int_{\beta}^{\beta+\delta}\int_{0}^{\tau}e^{-\frac{\gamma}{m}(\tau-s)}D(X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m})dB_{s}d\tau$
$\displaystyle=\int_{\beta}^{\beta+\delta}e^{-\frac{\gamma}{m}\tau}\overline{V}_{0}d\tau+\frac{\lambda}{m}\int_{0}^{\beta}\int_{\beta}^{\beta+\delta}e^{-\frac{\gamma}{m}(\tau-s)}d\tau(X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m})ds+\frac{\lambda}{m}\int_{\beta}^{\beta+\delta}\int_{s}^{\beta+\delta}e^{-\frac{\gamma}{m}(\tau-s)}d\tau(X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m})ds$
$\displaystyle\quad+\frac{\sigma}{m}\int_{0}^{\beta}\int_{\beta}^{\beta+\delta}e^{-\frac{\gamma}{m}(\tau-s)}d\tau
D(X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m})dB_{s}+\frac{\sigma}{m}\int_{\beta}^{\beta+\delta}\int_{s}^{\beta+\delta}e^{-\frac{\gamma}{m}(\tau-s)}d\tau
D(X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m})dB_{s}$
$\displaystyle=\frac{m}{\gamma}(e^{-\frac{\gamma}{m}\beta}-e^{-\frac{\gamma}{m}(\beta+\delta)})\overline{V}_{0}$
$\displaystyle\quad+\frac{\lambda}{\gamma}\int_{0}^{\beta}(e^{-\frac{\gamma}{m}(\beta-s)}-e^{-\frac{\gamma}{m}(\beta+\delta-s)})(X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m})ds+\frac{\lambda}{\gamma}\int_{\beta}^{\beta+\delta}(1-e^{-\frac{\gamma}{m}(\beta+\delta-s)})(X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m})ds$
$\displaystyle\quad+\frac{\sigma}{\gamma}\int_{0}^{\beta}(e^{-\frac{\gamma}{m}(\beta-s)}-e^{-\frac{\gamma}{m}(\beta+\delta-s)})D(X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m})dB_{s}+\frac{\sigma}{\gamma}\int_{\beta}^{\beta+\delta}(1-e^{-\frac{\gamma}{m}(\beta+\delta-s)})D(X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m})dB_{s}\,.$
(2.15)

Note that there holds $|e^{-x}-e^{-y}|\leq|x-y|\wedge 1$ for all
$x,y\in[0,\infty)$. Basic computations further indicate that for each $q\geq
1$ and $\tau\in[0,T]$,

$\displaystyle\int_{0}^{\tau}\left|e^{-\frac{\gamma(\tau-s)}{m}}-e^{-\frac{\gamma(\tau+\delta-s)}{m}}\right|^{q}\,ds\leq\int_{0}^{\tau}\left(e^{-\frac{\gamma(\tau-s)}{m}}-e^{-\frac{\gamma(\tau+\delta-s)}{m}}\right)\,ds$
$\displaystyle=\frac{m}{\gamma}\left(1-e^{-\frac{\gamma\delta}{m}}\right)-\frac{m}{\gamma}\left(e^{-\frac{\gamma\tau}{m}}-e^{-\frac{\gamma(\tau+\delta)}{m}}\right)$
$\displaystyle\leq\frac{m}{\gamma}\cdot\frac{\gamma\delta}{m}$
$\displaystyle=\delta,$

and in particular,

$\int_{\beta}^{\beta+\delta}\left(1-e^{-\frac{\gamma(\beta+\delta-s)}{m}}\right)^{q}ds\leq\int_{\beta}^{\beta+\delta}1\,ds=\delta.$

Then, it is obvious that

$\displaystyle\mathbb{E}\left[\left|\frac{m}{\gamma}(e^{-\frac{\gamma}{m}\beta}-e^{-\frac{\gamma}{m}(\beta+\delta)})\overline{V}_{0}\right|^{2}\right]\leq\frac{m^{2}}{\gamma^{2}}\cdot\frac{\gamma^{2}\delta^{2}}{m^{2}}\left(\mathbb{E}[|\overline{V}_{0}|^{4}]\right)^{\frac{1}{2}}\leq\delta^{2}\left(\mathbb{E}[|\overline{V}_{0}|^{4}]\right)^{\frac{1}{2}}.$

Next, it follows that

$\displaystyle\mathbb{E}\left[\left|\int_{0}^{\beta}(e^{-\frac{\gamma}{m}(\beta-s)}-e^{-\frac{\gamma}{m}(\beta+\delta-s)})(X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m})ds\right|^{2}\right]$
$\displaystyle\leq\mathbb{E}\left[\int_{0}^{\beta}|e^{-\frac{\gamma}{m}(\beta-s)}-e^{-\frac{\gamma}{m}(\beta+\delta-s)}|^{2}ds\cdot\int_{0}^{\beta}|X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m}|^{2}ds\right]$
$\displaystyle\leq\delta\cdot
T\sup_{s\in[0,T]}\left(\mathbb{E}\left[|X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m}|^{4}\right]\right)^{1/2},$

and analogously,

$\displaystyle\mathbb{E}\left[\left|\int_{\beta}^{\beta+\delta}(1-e^{-\frac{\gamma}{m}(\beta+\delta-s)})(X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m})ds\right|^{2}\right]$
$\displaystyle\leq\mathbb{E}\left[\int_{\beta}^{\beta+\delta}\left(1-e^{-\frac{\gamma(\beta+\delta-s)}{m}}\right)^{2}ds\cdot\int_{\beta}^{\beta+\delta}|X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m}|^{2}ds\right]$
$\displaystyle\leq\delta\cdot\mathbb{E}\left[\int_{\beta}^{\beta+\delta}|X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m}|^{2}ds\right]$
$\displaystyle\leq\delta\cdot
T\sup_{s\in[0,T]}\left(\mathbb{E}\left[|X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m}|^{4}\right]\right)^{1/2},\,.$

Further, applying Itô’s isometry gives

$\displaystyle\mathbb{E}\left[\left|\int_{0}^{\beta}(e^{-\frac{\gamma}{m}(\beta-s)}-e^{-\frac{\gamma}{m}(\beta+\delta-s)})D(X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m})dB_{s}\right|^{2}\right]$
$\displaystyle\leq
d\mathbb{E}\left[\int_{0}^{\beta}|e^{-\frac{\gamma}{m}(\beta-s)}-e^{-\frac{\gamma}{m}(\beta+\delta-s)}|^{2}|X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m}|^{2}ds\right]$
$\displaystyle\leq
d\left(\mathbb{E}\left[\int_{0}^{\beta}|e^{-\frac{\gamma}{m}(\beta-s)}-e^{-\frac{\gamma}{m}(\beta+\delta-s)}|^{4}ds\right]\right)^{1/2}\cdot\left(\mathbb{E}\left[\int_{0}^{\beta}|X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m}|^{4}ds\right]\right)^{1/2}$
$\displaystyle\leq
d\delta^{1/2}\left(T\sup_{s\in[0,T]}\mathbb{E}\left[|X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m}|^{4}\right]\right)^{1/2},$

and analogously,

$\displaystyle\mathbb{E}\left[\left|\int_{\beta}^{\beta+\delta}(1-e^{-\frac{\gamma}{m}(\beta+\delta-s)})D(X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m})dB_{s}\right|^{2}\right]\leq
d\delta^{1/2}\left(T\sup_{s\in[0,T]}\mathbb{E}\left[|X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m}|^{4}\right]\right)^{1/2}.$

Therefore, summing up the above estimates and recalling $0<m\leq
m_{0}=\frac{1}{2}$, $\frac{1}{\gamma}\leq 2$, and the relations (2) and
(2.13), we arrive at

$\displaystyle\mathbb{E}[|\overline{X}_{\beta+\delta}^{m}-\overline{X}_{\beta}^{m}|^{2}]$
$\displaystyle\leq\frac{5}{\gamma^{2}}\delta^{2}(\mathbb{E}[|\overline{V}_{0}|^{4}])^{\frac{1}{2}}+\frac{10}{\gamma^{2}}\cdot\left(\lambda^{2}\delta
T+\sigma^{2}d\left(\delta
T\right)^{1/2}\right)\sup_{s\in[0,T]}\left(\mathbb{E}\left[|X_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m}|^{4}\right]\right)^{1/2}$
$\displaystyle\leq
C\left(\mathbb{E}[|\overline{X}_{0}|^{4}],\mathbb{E}[|\overline{V}_{0}|^{4}],C_{\alpha,\mathcal{E}},\lambda,\sigma,d,T\right)\left(\delta^{\frac{1}{2}}+\delta+\delta^{2}\right).$

Hence, for any $\varepsilon>0$, $\eta>0$, there exists some $\delta_{0}>0$
such that for all $0<m\leq\frac{1}{2}$ it holds that

$\sup_{\delta\in[0,\delta_{0}]}\mathbb{P}(|\overline{X}_{\beta+\delta}^{m}-\overline{X}_{\beta}^{m}|^{2}\geq\eta)\leq\sup_{\delta\in[0,\delta_{0}]}\frac{\mathbb{E}[|\overline{X}_{\beta+\delta}^{m}-\overline{X}_{\beta}^{m}|^{2}]}{\eta}\leq\varepsilon\,.$
(2.16)

This justifies condition $Con2$ in Lemma 2.1. ∎

Next we shall identify the limit process, before which we recall a lemma on
the stability estimate of the nonlinear term $X^{\alpha}(\rho)$.

######  Lemma 2.2.

[carrillo2018analytical, Lemma 3.2] Assume that
$\rho,\widehat{\rho}\in\mathcal{P}_{4}(\mathbb{R}^{d})$. Then the following
stability estimate holds

$|X^{\alpha}(\rho)-X^{\alpha}(\widehat{\rho})|\leq
CW_{2}(\rho,\widehat{\rho})\,,$ (2.17)

where $W_{2}$ is the $2$-Wasserstein distance, and $C$ depends only on
$\alpha,L$, $\int_{\mathbb{R}^{d}}|x|^{4}\rho(dx)$, and
$\int_{\mathbb{R}^{d}}|x|^{4}\widehat{\rho}(dx)$.

######  Theorem 2.3 (Zero-inertia limit).

Let Assumption 1 hold and $(X_{t}^{m},V_{t}^{m})_{t\in[0,T]}$ satisfy the
system (1.7). Then as $m\rightarrow 0^{+}$, the sequence of stochastic
processes $\\{\overline{X}^{m}\\}_{0<m\leq 1}$ converge weakly to
$\overline{X}$, which is the unique solution to the following SDE:

$\displaystyle\overline{X}_{t}=\overline{X}_{0}+\lambda\int_{0}^{t}(X_{s}^{\alpha}(\rho)-\overline{X}_{s})ds+\sigma\int_{0}^{t}D(X_{s}^{\alpha}(\rho)-\overline{X}_{s})dB_{s}\,.$
(2.18)

Moreover it holds that

$\sup_{t\in[0,T]}\mathbb{E}[|\overline{X}^{m}_{t}-\overline{X}_{t}|^{2}]\leq
C\,m\,,$ (2.19)

where the constant $C$ depends only on
$\mathbb{E}[|\overline{X}_{0}|^{4}],\mathbb{E}[|\overline{V}_{0}|^{4}],C_{\alpha,\mathcal{E}},\lambda,\sigma,d$,
and $T$.

######  Remark 2.1.

It follows from the definition of Wasserstein distance that

$\sup_{t\in[0,T]}W_{2}^{2}(\rho^{m}_{t},\rho_{t})\leq\sup_{t\in[0,T]}\mathbb{E}[|\overline{X}^{m}_{t}-\overline{X}_{t}|^{2}]\leq
C\,m\,,$ (2.20)

which in a way is consistent with the result obtained in [choi2020quantified,
Theorem 1.3], where the authors obtained a quantified overdamped limit (with
the same rate $m$) of the singular Vlasov-Poisson-Fokker-Planck system to the
aggregation-diffusion equation. Besides, the obtained weak convergence of
$\overline{X}^{m}\rightharpoonup\overline{X}$ is in the path space
$\mathcal{C}([0,T];\mathbb{R}^{d})$, which implies and is obviously stronger
than the convergence of $\\{\rho^{m}_{t}\\}_{m>0}$ to $\rho_{t}$ for each time
$t\geq 0$.

###### Proof.

By Theorem 2.2, each subsequence $\\{\overline{X}^{m_{k}}\\}_{k\in\mathbb{N}}$
with $m_{0}\leq 1/2$ and $m_{k}$ converging to $0$ as $k\rightarrow\infty$
admits a subsequence (denoted w.l.o.g. by itself) that converges weakly. By
Skorokhod’s lemma (see [billingsley1999convergence, Theorem 6.7 on page 70]),
we may find a common probability space $(\Omega,\mathcal{F},\mathbb{P})$ on
which the processes $\\{\overline{X}^{m_{k}}\\}_{k\in\mathbb{N}}$ converge to
some process $\widehat{X}$ as random variables valued in
$\mathcal{C}([0,T];\mathbb{R}^{d})$ almost surely. In particular, we have

$\displaystyle\mathbb{P}\left(\lim_{k\rightarrow\infty}|\overline{X}^{m_{k}}_{t}-\widehat{X}_{t}|=0\right)=1,\,\,\forall\,t\in[0,T].$
(2.21)

We shall verify that the limit $\widehat{X}$ is indeed the unique solution
$\overline{X}$ to SDE (2.18).

Recall the SDE satisfied by $\overline{X}^{m_{k}}$ in (2)

$\displaystyle\overline{X}_{t}^{m_{k}}$
$\displaystyle=\overline{X}_{0}+\frac{m_{k}}{\gamma}(1-e^{-\frac{\gamma}{{m_{k}}}t})\overline{V}_{0}+\frac{\lambda}{\gamma}\int_{0}^{t}(1-e^{-\frac{\gamma}{{m_{k}}}(t-s)})(X_{s}^{\alpha}(\rho^{m_{k}})-\overline{X}_{s}^{m_{k}})ds$
$\displaystyle\quad+\frac{\sigma}{\gamma}\int_{0}^{t}(1-e^{-\frac{\gamma}{{m_{k}}}(t-s)})D(X_{s}^{\alpha}(\rho^{m_{k}})-\overline{X}_{s}^{m_{k}})dB_{s}\,.$
(2.22)

By the estimates in (2.13) and (ii) of Theorem 2.1, there exists a constant
$C_{2}$ being independent of $m_{k}$ such that

$\sup_{k\in\mathbb{N}}\sup_{t\in[0,T]}\mathbb{E}\left[|\overline{X}_{t}^{m_{k}}|^{4}\right]+\sup_{t\in[0,T]}\mathbb{E}\left[\left|\overline{X}_{t}\right|^{4}\right]\leq
C_{2}:=C(\mathbb{E}[|\overline{X}_{0}|^{4}],\mathbb{E}[|\overline{V}_{0}|^{4}],C_{\alpha,\mathcal{E}},\lambda,\sigma,d,T)<\infty.$
(2.23)

Letting $\rho(t,dx)$ be the probability distribution of $\overline{X}_{t}$ for
$t\in[0,T]$, we have

$|{X}_{t}^{\alpha}(\rho)|=\left|\frac{\int_{\mathbb{R}^{d}}x\omega_{\alpha}^{\mathcal{E}}(x)\rho(t,dx)}{\int_{\mathbb{R}^{d}}\omega_{\alpha}^{\mathcal{E}}(x)\rho(t,dx)}\right|\leq
C_{\alpha,\mathcal{E}}\int_{\mathbb{R}^{d}}|x|\rho(t,dx)\leq
C_{\alpha,\mathcal{E}}(\mathbb{E}[|\overline{X}_{t}|^{4}])^{\frac{1}{4}},$

and

$\sup_{k\in\mathbb{N}}\sup_{t\in[0,T]}|{X}_{t}^{\alpha}(\rho^{m_{k}})|\leq
C_{\alpha,\mathcal{E}}(C_{2})^{\frac{1}{4}},\quad\text{ and
}\quad\sup_{t\in[0,T]}|{X}_{t}^{\alpha}(\rho)|\leq
C_{\alpha,\mathcal{E}}(C_{2})^{\frac{1}{4}}\,.$ (2.24)

Then we compare the SDEs (2.18) and (2.22) term by term. By Lemma 2.2, we have

$|X_{t}^{\alpha}(\rho^{m_{k}})-X_{t}^{\alpha}(\rho)|^{2}\leq
CW_{2}^{2}(\rho_{t}^{m_{k}},\rho_{t})\leq
C\mathbb{E}[|\overline{X}_{t}^{m_{k}}-\overline{X}_{t}|^{2}],$

and thus,

$\displaystyle\mathbb{E}\left[\left|\frac{\lambda}{\gamma}\int_{0}^{t}(1-e^{-\frac{\gamma}{{m_{k}}}(t-s)})(X_{s}^{\alpha}(\rho^{m_{k}}_{s})-\overline{X}_{s}^{m_{k}})ds-\lambda\int_{0}^{t}(X_{s}^{\alpha}(\rho_{s})-\overline{X}_{s})ds\right|^{2}\right]$
$\displaystyle\leq
2\mathbb{E}\left[\left|\frac{\lambda}{1-m_{k}}\int_{0}^{t}(1-e^{-\frac{1-m_{k}}{{m_{k}}}(t-s)})(X_{s}^{\alpha}(\rho^{m_{k}}_{s})-X_{s}^{\alpha}(\rho_{s})+\overline{X}_{s}-\overline{X}_{s}^{m_{k}})ds\right|^{2}\right]$
$\displaystyle\quad+2\mathbb{E}\left[\left|{\lambda}\int_{0}^{t}\left(\frac{1-e^{-\frac{1-m_{k}}{{m_{k}}}(t-s)}}{1-m_{k}}-1\right)(X_{s}^{\alpha}(\rho_{s})-\overline{X}_{s})ds\right|^{2}\right]$
$\displaystyle\leq
C\mathbb{E}\left[\int_{0}^{t}\left|\overline{X}_{s}-\overline{X}_{s}^{m_{k}}\right|^{2}ds\right]+C{\lambda^{2}}\int_{0}^{t}\left|\frac{1-e^{-\frac{1-m_{k}}{{m_{k}}}(t-s)}}{1-m_{k}}-1\right|^{2}ds\cdot\mathbb{E}\left[\int_{0}^{T}\left|X_{s}^{\alpha}(\rho_{s})-\overline{X}_{s}\right|^{2}ds\right]$
$\displaystyle\leq
C\mathbb{E}\left[\int_{0}^{t}\left|\overline{X}_{s}-\overline{X}_{s}^{m_{k}}\right|^{2}ds\right]+C\int_{0}^{t}\left|\frac{1-e^{-\frac{1-m_{k}}{{m_{k}}}(t-s)}-(1-m_{k})}{1-m_{k}}\right|^{2}ds$
$\displaystyle\leq
C\mathbb{E}\left[\int_{0}^{t}\left|\overline{X}_{s}-\overline{X}_{s}^{m_{k}}\right|^{2}ds\right]+C\int_{0}^{t}\left(\left|{m_{k}}\right|^{2}+e^{-\frac{2(1-m_{k})}{{m_{k}}}(t-s)}\right)ds$
$\displaystyle\leq
C\mathbb{E}\left[\int_{0}^{t}\left|\overline{X}_{s}-\overline{X}_{s}^{m_{k}}\right|^{2}ds\right]+C\left(t\left|{m_{k}}\right|^{2}+\frac{m_{k}}{2(1-m_{k})}\right),$
(2.25)

where the constants $C$s are independent of $k$. For the stochastic integrals,
it holds analogously that

$\displaystyle\mathbb{E}\left[\left|\frac{\sigma}{\gamma}\int_{0}^{t}(1-e^{-\frac{\gamma}{{m_{k}}}(t-s)})D(X_{s}^{\alpha}(\rho^{m_{k}}_{s})-\overline{X}_{s}^{m_{k}})dB_{s}-\sigma\int_{0}^{t}D(X_{s}^{\alpha}(\rho_{s})-\overline{X}_{s})dB_{s}\right|^{2}\right]$
$\displaystyle\leq{d\sigma^{2}}\sum_{n=1}^{d}\mathbb{E}\left[\left|\frac{1}{\gamma}\int_{0}^{t}(1-e^{-\frac{\gamma}{{m_{k}}}(t-s)})(X_{s}^{\alpha}(\rho^{m_{k}}_{s})-\overline{X}_{s}^{m_{k}})_{n}dB_{s}^{n}e_{n}-\int_{0}^{t}(X_{s}^{\alpha}(\rho_{s})-\overline{X}_{s})_{n}dB_{s}^{n}e_{n}\right|^{2}\right]$
$\displaystyle=d{\sigma^{2}}\sum_{n=1}^{d}\mathbb{E}\left[\int_{0}^{t}\left|\frac{1-e^{-\frac{\gamma}{{m_{k}}}(t-s)}}{\gamma}(X_{s}^{\alpha}(\rho^{m_{k}}_{s})-\overline{X}_{s}^{m_{k}})_{n}-(X_{s}^{\alpha}(\rho_{s})-\overline{X}_{s})_{n}\right|^{2}ds\right]$
$\displaystyle\leq{2d\sigma^{2}}\sum_{n=1}^{d}\mathbb{E}\left[\int_{0}^{t}\left|\frac{1-e^{-\frac{\gamma}{{m_{k}}}(t-s)}}{\gamma}\left((X_{s}^{\alpha}(\rho^{m_{k}}_{s})-\overline{X}_{s}^{m_{k}})_{n}-(X_{s}^{\alpha}(\rho_{s})-\overline{X}_{s})_{n}\right)\right|^{2}\,ds\right]$
$\displaystyle\quad+{2d\sigma^{2}}\sum_{n=1}^{d}\mathbb{E}\left[\int_{0}^{t}\left|\left(\frac{1-e^{-\frac{\gamma}{{m_{k}}}(t-s)}}{\gamma}-1\right)(X_{s}^{\alpha}(\rho_{s})-\overline{X}_{s})_{n}\right|^{2}ds\right]$
$\displaystyle\leq
C\mathbb{E}\left[\int_{0}^{t}\left|\overline{X}_{s}^{m_{k}}-\overline{X}_{s}\right|^{2}\,ds\right]+{2d\sigma^{2}}\sup_{s\in[0,t]}\mathbb{E}\left[\left|(X_{s}^{\alpha}(\rho_{s})-\overline{X}_{s})\right|^{2}\right]\cdot\int_{0}^{t}\left|\left(\frac{1-e^{-\frac{\gamma}{{m_{k}}}(t-s)}}{\gamma}-1\right)\right|^{2}ds$
$\displaystyle\leq
C\mathbb{E}\left[\int_{0}^{t}\left|\overline{X}_{s}-\overline{X}_{s}^{m_{k}}\right|^{2}ds\right]+C\left(t\left|{m_{k}}\right|^{2}+\frac{m_{k}}{2(1-m_{k})}\right).$
(2.26)

In addition, it is obvious that

$\displaystyle\left|\frac{m_{k}}{\gamma}(1-e^{-\frac{\gamma}{m_{k}}t})\overline{V}_{0}\right|\leq
Cm_{k}\left|\overline{V}_{0}\right|.$ (2.27)

Therefore, recalling $m_{k}\leq\frac{1}{2}$, combining the estimates
(2.25)-(2.27) and subtracting both sides of SDEs (2.18) from those of (2.22),
we have

$\displaystyle\mathbb{E}[|\overline{X}^{m_{k}}_{t}-\overline{X}_{t}|^{2}]\leq
C\int_{0}^{t}\mathbb{E}[|\overline{X}^{m_{k}}_{s}-\overline{X}_{s}|^{2}]ds+C\,m_{k},\quad
t\in[0,T].$

By Gronwall’s inequality it implies that

$\sup_{t\in[0,T]}\mathbb{E}[|\overline{X}^{m_{k}}_{t}-\overline{X}_{t}|^{2}]\leq
Cm_{k}\rightarrow 0,\text{ as }k\rightarrow\infty\,,$ (2.28)

where $C$ depends only on
$\mathbb{E}[|\overline{X}_{0}|^{4}],\mathbb{E}[|\overline{V}_{0}|^{4}],C_{\alpha,\mathcal{E}},\lambda,\sigma,d$,
and $T$. In view of both the convergences (2.21) and (2.28), we must have
$\widehat{X}=\overline{X}$. Finally, due to the arbitrariness of the
subsequence $\\{\overline{X}^{m_{k}}\\}_{k\in\mathbb{N}}$ and the uniqueness
of $\overline{X}$, we conclude that as $m\rightarrow 0^{+}$, the sequence of
stochastic processes $\\{\overline{X}^{m}\\}_{0<m\leq 1}$ converge weakly to
the unique solution $\overline{X}$ to SDE (2.18), with the estimate (2.19)
following in the same way as (2.28).

∎

######  Remark 2.2.

When proving the convergence of $\\{\overline{X}^{m}\\}$ satisfying SDEs (1.7)
to the solution $\overline{X}$ of (1.10), we cannot expect the convergence of
the associated velocity processes $\\{\overline{V}^{m}\\}$ due to the
indifferentiability of the limit $\\{\overline{X}_{t}\\}_{t\geq 0}$ with
respect to time $t$ if $\sigma\neq 0$. Therefore, we do not investigate
convergence of the joint Markovian process
$\\{(\overline{X}^{m},\overline{V}^{m})\\}$ and consider instead solely the
process $\\{\overline{X}^{m}\\}$ which satisfies a stochastic integral
equation (2.22) of Volterra type, being path-dependent and thus non-Markovian.
This non-Markovianity prevents us from using the usual techniques for weak
convergence with martingale problems but prompts us to identify the limit by
measuring directly the distance between $\overline{X}^{m_{k}}$ and
$\overline{X}$ in the above proof.

##  3\. Generalization to the case with memory effects

In [grassi2020particle], the authors considered a PSO model which involves the
memory of the local best positions, and it is of the form

$\displaystyle d\overline{X}_{t}^{m}=\overline{V}_{t}^{m}dt,$ (3.1)
$\displaystyle
d\overline{Y}_{t}^{m}=\nu\left(\overline{X}_{t}^{m}-\overline{Y}_{t}^{m}\right)S^{\beta}\left(\overline{X}_{t}^{m},\overline{Y}_{t}^{m}\right)dt,$
(3.2) $\displaystyle
d\overline{V}_{t}^{m}=-\frac{\gamma}{m}\overline{V}_{t}^{m}dt+\frac{\lambda_{1}}{m}\left(\overline{Y}_{t}^{m}-\overline{X}_{t}^{m}\right)dt+\frac{\lambda_{2}}{m}\left({Y}_{t}^{\alpha}(\overline{\rho}^{m})-\overline{X}_{t}^{m}\right)dt$
$\displaystyle\quad\qquad+\frac{\sigma_{1}}{m}D\left(\overline{Y}_{t}^{m}-\overline{X}_{t}^{m}\right)dB_{t}^{1}+\frac{\sigma_{2}}{m}D\left({Y}_{t}^{\alpha}(\overline{\rho}^{m})-\overline{X}_{t}^{m}\right)dB_{t}^{2}\,,$
(3.3)

where $B^{1}$ and $B^{2}$ are two mutually independent d-dimensional Wiener
processes and similarly to the previous section, we introduce the following
regularization of the global best position

${Y}_{t}^{\alpha}(\overline{\rho}^{m})=\frac{\int_{\mathbb{R}^{d}}y\omega_{\alpha}(y){\overline{\rho}}^{m}(t,dy)}{\int_{\mathbb{R}^{d}}\omega_{\alpha}(y){\overline{\rho}}^{m}(t,dy)},\quad{\overline{\rho}}^{m}(t,y)=\iint_{\mathbb{R}^{d}\times\mathbb{R}^{d}}f^{m}(t,dx,y,dv)\,.$
(3.4)

Here the equation (3.2) of $\overline{Y}^{m}$ is the time continuous
approximation to the evolution of the local best position, and $S^{\beta}$
with $\beta\gg 1$ is hyperbolic tangent
$S^{\beta}(x,y)=\tanh(\beta(\mathcal{E}(x)-\mathcal{E}(y))$. The corresponding
mean-field PSO equation is

$\displaystyle\partial_{t}f_{t}^{m}+v\cdot\nabla_{x}f_{t}^{m}$
$\displaystyle+\nabla_{y}\cdot(\nu(x-y)S^{\beta}(x,y)f_{t}^{m})=\nabla_{v}\cdot\bigg{(}\frac{\gamma}{m}vf_{t}^{m}+\frac{\lambda_{1}}{m}(x-y)f_{t}^{m}+\frac{\lambda_{2}}{m}(x-{Y}_{t}^{\alpha}({\overline{\rho}}^{m}))f_{t}^{m}$
$\displaystyle+\Big{(}\frac{\sigma_{2}^{2}}{2m^{2}}D(x-{Y}_{t}^{\alpha}({\overline{\rho}}^{m}))^{2}+\frac{\sigma_{1}^{2}}{2m^{2}}D(x-y)^{2}\Big{)}\nabla_{v}f_{t}^{m}\bigg{)}\,.$
(3.5)

We want to prove that the zero-inertia limit $(m\to 0)$ leads to the following
mean-field CBO dynamic

$\displaystyle\begin{cases}\overline{X}_{t}=\overline{X}_{0}+\lambda_{1}\int_{0}^{t}(\overline{Y}_{s}-\overline{X}_{s})ds+\sigma_{1}\int_{0}^{t}D(\overline{Y}_{s}-\overline{X}_{s})dB_{s}^{1}+\lambda_{2}\int_{0}^{t}(Y_{s}^{\alpha}(\overline{\rho})-\overline{X}_{s})ds+\sigma_{2}\int_{0}^{t}D(Y_{s}^{\alpha}(\overline{\rho})-\overline{X}_{s})dB_{s}^{2},\\\
\overline{Y}_{t}=\overline{Y}_{0}+\nu\int_{0}^{t}\left(\overline{X}_{s}-\overline{Y}_{s}\right)S^{\beta}\left(\overline{X}_{s},\overline{Y}_{s}\right)ds\,,\end{cases}$

and its corresponding partial differential equation is

$\displaystyle\partial_{t}\rho_{t}+\nabla_{x}\cdot\left(\lambda_{1}(y-x)+\lambda_{2}({Y}_{t}^{\alpha}({\overline{\rho}})-x)\rho_{t}\right)+\nabla_{y}\cdot\left(\nu(x-y)S^{\beta}(x,y))\rho_{t}\right)$
$\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{j=1}^{d}\frac{\partial^{2}}{\partial
x_{j}^{2}}\left(\rho_{t}(\sigma_{1}^{2}(x-y)_{j}^{2}+\sigma_{2}^{2}(x-{Y}_{t}^{\alpha}({\overline{\rho}}))_{j}^{2})\right)\,,$
(3.6)

where $\overline{\rho}(t,y)=\int_{\mathbb{R}^{d}}\rho(t,dx,y)$.

Since the proof of the zero-inertia limit for the PSO dynamics with memory
effects follows similar arguments as developed in the previous section and no
essential innovation is needed to be explained, we only sketch the proof for
the tightness.

######  Theorem 3.1 (Tightness).

Let Assumption 1 hold and
$(\overline{X}_{t}^{m},\overline{Y}_{t}^{m},\overline{V}_{t}^{m})_{t\in[0,T]}$
satisfy the system (3.1)-(3.3). For each countable subsequence
$\\{m_{k}\\}_{k\in\mathbb{N}}\subset[0,1]$ with
$\lim_{k\rightarrow\infty}m_{k}=0$, the sequence of probability distributions
$\\{\rho^{m_{k}}\\}_{k\in\mathbb{N}}$ of
$\\{(\overline{X}^{m_{k}},\overline{Y}^{m_{k}})\\}_{k\in\mathbb{N}}$ is tight.

###### Proof.

The proof is similar to Theorem 2.2.

$\bullet$ Step 1: Checking $(Con1)$.  For $0<m\leq\frac{1}{2}$, we first solve
(3.3) for $\overline{V}^{m}$ and obtain

$\displaystyle\overline{V}_{t}^{m}$
$\displaystyle=e^{-\frac{\gamma}{m}t}\overline{V}_{0}+\frac{\lambda_{1}}{m}\int_{0}^{t}e^{-\frac{\gamma}{m}(t-s)}(\overline{Y}_{s}^{m}-\overline{X}_{s}^{m})ds+\frac{\sigma_{1}}{m}\int_{0}^{t}e^{-\frac{\gamma}{m}(t-s)}D(\overline{Y}_{s}^{m}-\overline{X}_{s}^{m})dB_{s}^{1}$
$\displaystyle\quad+\frac{\lambda_{2}}{m}\int_{0}^{t}e^{-\frac{\gamma}{m}(t-s)}(Y_{s}^{\alpha}(\overline{\rho}^{m})-\overline{X}_{s}^{m})ds+\frac{\sigma_{2}}{m}\int_{0}^{t}e^{-\frac{\gamma}{m}(t-s)}D(Y_{s}^{\alpha}(\overline{\rho}^{m})-\overline{X}_{s}^{m})dB_{s}^{2}\,.$

Here $\overline{\rho}^{m}(t,y)=\int_{\mathbb{R}^{d}}\rho^{m}(t,dx,y)$. By
Fubini’s theorem, similar arguments as in (2) yield that

$\displaystyle\overline{X}_{t}^{m}$
$\displaystyle=\overline{X}_{0}+\frac{m}{\gamma}(1-e^{-\frac{\gamma}{m}t})\overline{V}_{0}+\frac{\lambda_{1}}{\gamma}\int_{0}^{t}(1-e^{-\frac{\gamma}{m}(t-s)})(\overline{Y}_{s}^{m}-\overline{X}_{s}^{m})ds+\frac{\sigma_{1}}{\gamma}\int_{0}^{t}(1-e^{-\frac{\gamma}{m}(t-s)})D(\overline{Y}_{s}^{m}-\overline{X}_{s}^{m})dB_{s}^{1}$
$\displaystyle\quad+\frac{\lambda_{2}}{\gamma}\int_{0}^{t}(1-e^{-\frac{\gamma}{m}(t-s)})(Y_{s}^{\alpha}(\overline{\rho}^{m})-\overline{X}_{s}^{m})ds+\frac{\sigma_{2}}{\gamma}\int_{0}^{t}(1-e^{-\frac{\gamma}{m}(t-s)})D(Y_{s}^{\alpha}(\overline{\rho}^{m})-\overline{X}_{s}^{m})dB_{s}^{2}\,.$
(3.7)

Following the same computations as in Theorem 2.2 gives

$\displaystyle\mathbb{E}[|\overline{X}_{t}^{m}|^{4}]\leq
C\mathbb{E}[|\overline{X}_{0}|^{4}+|\overline{V}_{0}|^{4}]+C\int_{0}^{t}\mathbb{E}[|Y_{s}^{\alpha}(\rho^{m})-\overline{X}_{s}^{m}|^{4}]ds+C\int_{0}^{t}\mathbb{E}[|\overline{Y}_{s}^{m}-\overline{X}_{s}^{m}|^{4}]ds\,,$

where $C$ depends only on $\lambda_{1},\sigma_{2},\lambda_{2},\sigma_{2},d$,
and $T$. Put
$\widetilde{\rho}^{m}(t,x)=\int_{\mathbb{R}^{d}}\rho^{m}(t,x,dy)$. In a
similar way to (2) we have

$\displaystyle\mathbb{E}[|Y_{t}^{\alpha}(\overline{\rho}^{m})-\overline{X}_{t}^{m}|^{4}]$
$\displaystyle=\int_{\mathbb{R}^{d}}\left|\frac{\int_{\mathbb{R}^{d}}y\omega_{\alpha}^{\mathcal{E}}(y)\overline{\rho}^{m}(t,dy)}{\int_{\mathbb{R}^{d}}\omega_{\alpha}^{\mathcal{E}}(y)\overline{\rho}^{m}(t,dy)}-x\right|^{4}\widetilde{\rho}^{m}(t,dx)=\int_{\mathbb{R}^{d}}\left|\frac{\int_{\mathbb{R}^{d}}(y-x)\omega_{\alpha}^{\mathcal{E}}(y)\overline{\rho}^{m}(t,dy)}{\int_{\mathbb{R}^{d}}\omega_{\alpha}^{\mathcal{E}}(y)\overline{\rho}^{m}(t,dy)}\right|^{4}\widetilde{\rho}^{m}(t,dx)$
$\displaystyle\leq\frac{\int_{\mathbb{R}^{d}}\int_{\mathbb{R}^{d}}|x-y|^{4}\omega_{\alpha}^{\mathcal{E}}(y)\overline{\rho}^{m}(t,dy)\widetilde{\rho}^{m}(t,dx)}{\int_{\mathbb{R}^{d}}\omega_{\alpha}^{\mathcal{E}}(y)\overline{\rho}^{m}(t,dy)}\leq
8C_{\alpha,\mathcal{E}}\mathbb{E}[|\overline{X}_{t}^{m}|^{4}+|\overline{Y}_{t}^{m}|^{4}]\,.$

Thus it yields that

$\displaystyle\mathbb{E}[|\overline{X}_{t}^{m}|^{4}]\leq
C\mathbb{E}[|\overline{X}_{0}|^{4}+|\overline{V}_{0}|^{4}]+C\int_{0}^{t}\mathbb{E}[|\overline{Y}_{s}^{m}|^{4}+|\overline{X}_{s}^{m}|^{4}]ds\,,$
(3.8)

where $C$ depends only on $\lambda_{1},\sigma_{2},\lambda_{2},\sigma_{2},d,T$
and $C_{\alpha,\mathcal{E}}$.

Recall that

$\displaystyle\overline{Y}_{t}^{m}=\overline{Y}_{0}+\nu\int_{0}^{t}\left(\overline{X}_{s}^{m}-\overline{Y}_{s}^{m}\right)S^{\beta}\left(\overline{X}_{s}^{m},\overline{Y}_{s}^{m}\right)ds$

with $S^{\beta}(x,y)=\tanh(\beta(\mathcal{E}(x)-\mathcal{E}(y))$. Using the
fact that $|S^{\beta}|\leq 1$ then it follows

$\mathbb{E}[|\overline{Y}_{t}^{m}|^{4}]\leq
C\mathbb{E}[|\overline{Y}_{0}|^{4}]+C\int_{0}^{t}\mathbb{E}[|\overline{Y}_{s}^{m}|^{4}+|\overline{X}_{s}^{m}|^{4}]ds\,,$

where $C$ depends only on $\nu$ and $T$. This together with (3.8) implies

$\mathbb{E}[|\overline{X}_{t}^{m}|^{4}+|\overline{Y}_{t}^{m}|^{4}]\leq
C\mathbb{E}[|\overline{X}_{0}|^{4}+|\overline{Y}_{0}|^{4}+|\overline{V}_{0}|^{4}]+C\int_{0}^{t}\mathbb{E}[|\overline{Y}_{s}^{m}|^{4}+|\overline{X}_{s}^{m}|^{4}]ds\,.$

By Gronwall’s inequality it yields that

$\sup_{t\in[0,T]}\mathbb{E}[|\overline{X}_{t}^{m}|^{4}+|\overline{Y}_{t}^{m}|^{4}]\leq
C(\mathbb{E}[|\overline{X}_{0}|^{4}+|\overline{Y}_{0}|^{4}+|\overline{V}_{0}|^{4}],\lambda_{1},\sigma_{2},\lambda_{2},\sigma_{2},d,T,C_{\alpha,\mathcal{E}},\nu)\,,$
(3.9)

which verifies $(Con1)$ for the case of $0<m\leq\frac{1}{2}$. We omit the
discussions for the case of $\frac{1}{2}<m\leq 1$.

$\bullet$ Step 2: Checking $(Con2)$.  Let $\beta$ be a
$\sigma(X^{m}_{s};s\in[0,T])$-stopping time with discrete values such that
$\beta+\delta_{0}\leq T$. Set $m_{0}=\frac{1}{2}$ w.l.o.g.. Then for all
$0<m\leq m_{0}$, one has $\frac{1}{2}\leq\gamma<1$. Similar to (2), one has

$\displaystyle\overline{X}_{\beta+\delta}^{m}-\overline{X}_{\beta}^{m}$
$\displaystyle=\frac{m}{\gamma}(e^{-\frac{\gamma}{m}\beta}-e^{-\frac{\gamma}{m}(\beta+\delta)})\overline{V}_{0}$
$\displaystyle\quad+\frac{\lambda_{2}}{\gamma}\int_{0}^{\beta}(e^{-\frac{\gamma}{m}(\beta-s)}-e^{-\frac{\gamma}{m}(\beta+\delta-s)})(Y_{s}^{\alpha}(\overline{\rho}^{m})-\overline{X}_{s}^{m})ds+\frac{\lambda_{2}}{\gamma}\int_{\beta}^{\beta+\delta}(1-e^{-\frac{\gamma}{m}(\beta+\delta-s)})(Y_{s}^{\alpha}(\overline{\rho}^{m})-\overline{X}_{s}^{m})ds$
$\displaystyle\quad+\frac{\sigma_{2}}{\gamma}\int_{0}^{\beta}(e^{-\frac{\gamma}{m}(\beta-s)}-e^{-\frac{\gamma}{m}(\beta+\delta-s)})D(Y_{s}^{\alpha}(\overline{\rho}^{m})-\overline{X}_{s}^{m})dB_{s}^{2}+\frac{\sigma_{2}}{\gamma}\int_{\beta}^{\beta+\delta}(1-e^{-\frac{\gamma}{m}(\beta+\delta-s)})D(Y_{s}^{\alpha}(\overline{\rho}^{m})-\overline{X}_{s}^{m})dB_{s}^{2}$
$\displaystyle\quad+\frac{\lambda_{1}}{\gamma}\int_{0}^{\beta}(e^{-\frac{\gamma}{m}(\beta-s)}-e^{-\frac{\gamma}{m}(\beta+\delta-s)})(\overline{Y}_{s}^{m}-\overline{X}_{s}^{m})ds+\frac{\lambda_{1}}{\gamma}\int_{\beta}^{\beta+\delta}(1-e^{-\frac{\gamma}{m}(\beta+\delta-s)})(\overline{Y}_{s}^{m}-\overline{X}_{s}^{m})ds$
$\displaystyle\quad+\frac{\sigma_{1}}{\gamma}\int_{0}^{\beta}(e^{-\frac{\gamma}{m}(\beta-s)}-e^{-\frac{\gamma}{m}(\beta+\delta-s)})D(\overline{Y}_{s}^{m}-\overline{X}_{s}^{m})dB_{s}^{1}+\frac{\sigma_{1}}{\gamma}\int_{\beta}^{\beta+\delta}(1-e^{-\frac{\gamma}{m}(\beta+\delta-s)})D(\overline{Y}_{s}^{m}-\overline{X}_{s}^{m})dB_{s}^{1}\,.$
(3.10)

Using the estimate (3.9) it follows from the same computations in Theorem 2.2
that

$\displaystyle\mathbb{E}[|\overline{X}_{\beta+\delta}^{m}-\overline{X}_{\beta}^{m}|^{2}]\leq
C\left(\delta^{\frac{1}{2}}+\delta+\delta^{2}\right)\,,$ (3.11)

where $C$ depends only on
$\mathbb{E}[|\overline{X}_{0}|^{4}+|\overline{Y}_{0}|^{4}+|\overline{V}_{0}|^{4}],\lambda_{1},\sigma_{2},\lambda_{2},\sigma_{2},d,T,C_{\alpha,\mathcal{E}}$,
and $\nu$.

Having a look at

$\displaystyle\overline{Y}_{\beta+\delta}^{m}-\overline{Y}_{\beta}^{m}=\nu\int_{\beta}^{\beta+\delta}\left(\overline{X}_{s}^{m}-\overline{Y}_{s}^{m}\right)S^{\beta}\left(\overline{X}_{s}^{m},\overline{Y}_{s}^{m}\right)ds,$

we have

$|\overline{Y}_{\beta+\delta}^{m}-\overline{Y}_{\beta}^{m}|^{2}\leq\nu^{2}\delta\int_{0}^{T}|\overline{X}_{s}^{m}-\overline{Y}_{s}^{m}|^{2}dt\,.$

By estimate (3.9), we have

$\mathbb{E}[|\overline{Y}_{\beta+\delta}^{m}-\overline{Y}_{\beta}^{m}|^{2}]\leq\nu^{2}\delta\int_{0}^{T}\mathbb{E}[|\overline{X}_{s}^{m}-\overline{Y}_{s}^{m}|^{4}]^{\frac{1}{2}}dt\leq
C\delta\,,$ (3.12)

where $C$ depends on
$\mathbb{E}[|\overline{X}_{0}|^{4}+|\overline{Y}_{0}|^{4}+|\overline{V}_{0}|^{4}],\lambda_{1},\sigma_{2},\lambda_{2},\sigma_{2},d,T,C_{\alpha,\mathcal{E}}$
and $\nu$. This together with (3.11) justifies $(Con2)$. ∎

Let us recall

$\displaystyle\overline{X}_{t}^{m}$
$\displaystyle=\overline{X}_{0}+\frac{m}{\gamma}(1-e^{-\frac{\gamma}{m}t})\overline{V}_{0}+\frac{\lambda_{1}}{\gamma}\int_{0}^{t}(1-e^{-\frac{\gamma}{m}(t-s)})(\overline{Y}_{s}^{m}-\overline{X}_{s}^{m})ds+\frac{\sigma_{1}}{\gamma}\int_{0}^{t}(1-e^{-\frac{\gamma}{m}(t-s)})D(\overline{Y}_{s}^{m}-\overline{X}_{s}^{m})dB_{s}^{1}$
$\displaystyle\quad+\frac{\lambda_{2}}{\gamma}\int_{0}^{t}(1-e^{-\frac{\gamma}{m}(t-s)})(Y_{s}^{\alpha}(\overline{\rho}^{m})-\overline{X}_{s}^{m})ds+\frac{\sigma_{2}}{\gamma}\int_{0}^{t}(1-e^{-\frac{\gamma}{m}(t-s)})D(Y_{s}^{\alpha}(\overline{\rho}^{m})-\overline{X}_{s}^{m})dB_{s}^{2}$
(3.13)

and

$\displaystyle\overline{Y}_{t}^{m}=\overline{Y}_{0}+\nu\int_{0}^{t}\left(\overline{X}_{s}^{m}-\overline{Y}_{s}^{m}\right)S^{\beta}\left(\overline{X}_{s}^{m},\overline{Y}_{s}^{m}\right)ds\,.$
(3.14)

Then following the lines of the proof in Theorem 3.2, one can easily obtain

######  Theorem 3.2 (Zero-inertia limit).

Let Assumption 1 hold and
$(\overline{X}_{t}^{m},\overline{Y}_{t}^{m})_{t\in[0,T]}$ satisfy the system
(3)–(3.14). Then as $m\rightarrow 0^{+}$, the sequence of stochastic processes
$\\{(\overline{X}^{m},\overline{Y}^{m})\\}_{0<m\leq 1}$ converge weakly to
$(\overline{X},\overline{Y})$ which is the unique solution to the following
coupled SDE:

$\displaystyle\overline{X}_{t}$
$\displaystyle=\overline{X}_{0}+\lambda_{1}\int_{0}^{t}(\overline{Y}_{s}-\overline{X}_{s})ds+\sigma_{1}\int_{0}^{t}D(\overline{Y}_{s}-\overline{X}_{s})dB_{s}^{1}+\lambda_{2}\int_{0}^{t}(Y_{s}^{\alpha}(\overline{\rho})-\overline{X}_{s})ds+\sigma_{2}\int_{0}^{t}D(Y_{s}^{\alpha}(\overline{\rho})-\overline{X}_{s})dB_{s}^{2},$
$\displaystyle\overline{Y}_{t}$
$\displaystyle=\overline{Y}_{0}+\nu\int_{0}^{t}\left(\overline{X}_{s}-\overline{Y}_{s}\right)S^{\beta}\left(\overline{X}_{s},\overline{Y}_{s}\right)ds\,.$

Moreover it holds that

$\sup_{t\in[0,T]}\mathbb{E}\left[\big{|}\overline{X}^{m}_{t}-\overline{X}_{t}\big{|}^{2}+\big{|}\overline{Y}^{m}_{t}-\overline{Y}_{t}\big{|}^{2}\right]\leq
C\,m\,,$ (3.15)

where the constant $C$ depends only on
$\mathbb{E}[|\overline{X}_{0}|^{4}+|\overline{Y}_{0}|^{4}+|\overline{V}_{0}|^{4}],\lambda_{1},\sigma_{2},\lambda_{2},\sigma_{2},d,\beta,T,C_{\alpha,\mathcal{E}}$,
and $\nu$.

##  4\. Numerical examples on the zero-inertia limit

We conclude this paper with a few instructive numerical experiments on
validating the zero-inertia limit. We will focus on the mono-dimensional case
since it allows us to see more clearly how the distribution of particles
evolves in time depending on the inertia parameter $m$, and hence show the
zero-inertia limit. Different benchmark functions have been used and tested,
but we will report here the case of the Ackley function shown in Figure 3.
Following the same structure of the paper, we will first analyze the case
without memory effect and then we will generalize as in Section 3 to the case
with memory. Extensive discussions on other numerical implementations and
experiments are presented in [grassi2020particle].

Figure 3. Ackley function in a big (left) and small (right) domain with its
many local minima.

###  4.1. Small inertia limit without memory

Given the system of stochastic differential equations in (1.2), the particle
system can be solved by using a semi-implicit discretization scheme

$\displaystyle\begin{cases}X_{n+1}^{i,m}=X_{n}^{i,m}+\Delta tV_{n+1}^{i,m},\\\
V_{n+1}^{i,m}=\frac{m}{m+\gamma\Delta t}V_{n}^{i,m}+\frac{\lambda\Delta
t}{m+\gamma\Delta t}(X_{n}^{\alpha,m}-X_{n}^{i,m})+\frac{\sigma\sqrt{\Delta
t}}{m+\gamma\Delta t}D(X_{n}^{\alpha,m}-X_{n}^{i,m})\theta_{n}^{i},\quad
i=1,\cdots,N\,,\end{cases}$ (4.1)

where $X_{n}^{i,m}$ and $V_{n}^{i,m}$ are, respectively, the position and
velocity of the i-th particle at the discrete time $n\Delta t$ with $\Delta t$
being the time discretization, and the diagonal matrix
$D(X_{n}^{\alpha,m}-X_{n}^{i,m})$ simply coincides with
$X_{n}^{\alpha,m}-X_{n}^{i,m}$ as we are considering the mono-dimensional
case. Moreover, $X_{n}^{\alpha,m}$ is defined as in (1.3) and
$\theta^{i}_{n}\sim\mathcal{N}(0,1)$ $\forall i,n$. We compare this particle
system with the CBO dynamic of the form (1.10), which can be solved using the
Euler-Maruyama scheme

$X^{i}_{n+1}=X_{n}^{i}+\Delta t\lambda(X_{n}^{\alpha}-X_{n}^{i})+\sqrt{\Delta
t}\sigma(X_{n}^{\alpha}-X_{n}^{i})\theta_{i}^{n}\,.$ (4.2)

As already mentioned, we consider the minimization of the Ackley function with
minimum at $x=0$ and, starting from the same initial distribution of
particles, we solve the PSO system (4.1) for different inertia values. Then,
we compare the evolution of the distribution of particles with the one of the
particles moving according to CBO system (4.2). In order to be able to compare
the results, we fix the parameters $\lambda=1$, $\sigma=\frac{1}{\sqrt{3}}$
and $\alpha=30$, while $\theta_{n}^{i}$ are sampled from $\mathcal{N}(0,1)$
and fixed for each $i=1,...,N$ and $n\in[0,T/\Delta t]$. Moreover, $T$ is set
to 1 and the time discretization is $\Delta t=0.01$, with a total number of
particles $N=10^{4}$.

Figure 4. Comparison of the CBO (4.2) and PSO (4.1) dynamics for different
inertia values (which are changing over the rows) and at many time steps
(changing over columns), starting from a normal distribution.

Figure 4 shows in each row the evolution of the CBO distribution and the one
of the PSO system with $m$ fixed that is decreasing over rows. The initial
particles are always sampled from the same distribution, which is in this case
a Gaussian centered in $0$ with variance $1$. Clearly, the PSO system with
$m=0.8$ leads to the correct minimum in 0 at the final time step $t=1$, but
the distribution of the particles is different from the one of the CBO. While,
for any $t\in[0,T]$, if the inertia value is decreased to $0.1$, or even to
$0.001$, the two distributions, namely the one obtained via CBO and the PSO
one, are indistinguishable, as the last two rows of Figure 4 show.
These considerations are confirmed in Figure 5 where we compare the
distributions obtained in Figure 4 using the Wasserstein 2 distance between
the CBO distribution and the PSO distribution. On the left of Figure 5, the
Wasserstein distance is plotted for each time step. Moreover, since we want to
show the influence of the inertia parameter, we take the mean value of the
Wasserstein distance over all time steps and plot it as a function of the
inertia values. This is shown on the right of Figure 5 where we also add the
mean value of the Kullback-Leibler divergence since the latter is a well-known
measure used to compare distributions, especially in statistics.

Figure 5. Left: Wasserstein distance over time ; Right: Wasserstein distance
and Kullback-Leibler divergence (in mean) over the inertia values.

Moreover, since it is necessary to start with an initial distribution that is
close to the global minimizer, we also try to see what happens when the
initial distribution is a uniform distribution between $-3$ and $3$ and
compare the evolution of its particles according to the CBO and PSO dynamics,
with varying inertia parameters. The result is shown in Figure 6. In this
case, the difference between the CBO distribution and the one of the PSO
dynamics is way higher in the case of big inertia value, but, as before, goes
to zero as soon as $m$ converges to $0$.

Figure 6. Evolution of an initial uniform distribution according to CBO (4.2)
and PSO (4.1) dynamics and their comparison for different time steps (on the
columns) and different inertia values (on the rows).

###  4.2. Small inertia limit with memory effect

The PSO model which involves the memory of the local and global best
positions, underlying (3.1)–(3.3), can similarly be solved via

$\displaystyle\begin{cases}X_{n+1}^{i,m}=X_{n}^{i,m}+\Delta
tV_{n+1}^{i,m},\quad i=1,\cdots,N\\\ Y_{n+1}^{i,m}=\nu\Delta
t(X_{n+1}^{i,m}-Y_{n}^{i,m})S^{\beta}(X_{n+1}^{i,m},Y_{n}^{i,m}),\\\
V_{n+1}^{i,m}=\frac{m}{m+\gamma\Delta t}V_{n}^{i,m}+\frac{\lambda_{1}\Delta
t}{m+\gamma\Delta t}(Y_{n}^{i,m}-X_{n}^{i,m})+\frac{\lambda_{2}\Delta
t}{m+\gamma\Delta t}(Y_{n}^{\alpha,m}-X_{n}^{i,m})\\\
\qquad\quad+\frac{\sigma_{1}\sqrt{\Delta t}}{m+\gamma\Delta
t}D(Y_{n}^{i,m}-X_{n}^{i,m})\theta_{n}^{1,i}+\frac{\sigma_{2}\sqrt{\Delta
t}}{m+\gamma\Delta
t}D(Y_{n}^{\alpha,m}-X_{n}^{i,m})\theta_{n}^{2,i}\,,\end{cases}$ (4.3)

where $Y_{n}^{i,m}$ is the local best that the i-th particle has memory of,
and $Y_{n}^{\alpha,m}$ is the regularized global best, defined as in (3.4).
Clearly, the corresponding CBO dynamics is the following

$\displaystyle\begin{cases}X_{n}^{i}=X_{n}^{i}+\lambda_{1}\Delta
t(Y_{n}^{i}-X_{n}^{i})+\lambda_{2}\Delta
t(Y_{n}^{\alpha}-X_{n}^{i})+\sigma_{1}\sqrt{\Delta
t}D(Y_{n}^{i}-X_{n}^{i})\theta_{n}^{1,i}+\sigma_{2}\sqrt{\Delta
t}(Y_{n}^{\alpha}-X_{n}^{i})\theta_{n}^{2,i}\\\ Y_{n}^{i}=Y_{n}^{i}+\nu\Delta
t(X_{n}^{i}-Y_{n}^{i})S^{\beta}(X_{n}^{i},Y_{n}^{i})\end{cases}$ (4.4)

Once again, since we want to show the convergence of the PSO distribution
obtained from (4.3) with a small inertia value to the one attained via the CBO
system (4.4), we need to set some of the parameters to the same values in
order to be able to compare the results. Their values are the following

$\lambda_{1}=\lambda_{2}=1\quad\sigma_{1}=\sigma_{2}=\frac{1}{\sqrt{3}}\quad\alpha=30\quad\beta=30\quad\nu=\frac{1}{2}$
(4.5)

and, as before, the effect of the Brownian motion leads to
$\theta_{n}^{1,i},\theta_{n}^{2,i}$ which are sampled from a normal
distribution and set to a fixed value $\forall i=1,..,N$ and $\forall
n\in[0,T/\Delta t]$. The time discretization and number of particles are set
to the same values as in the case without memory, namely $T=1$, $\Delta
t=0.01$, and $N=10^{4}$. The difference with the previous case is that now the
distribution of which we want to show convergence, is actually a function of
both the particles’ position and their local best and, as such, it is bi-
dimensional. The result of the evolution of the particles according to the CBO
dynamics (4.4) is shown in the first row of Figure 7, while the second row
represents the evolution according to the PSO dynamics (4.3) with a big
inertia parameter, e.g. $m=0.8$. This is compared with the evolution presented
in the third row in which the inertia is set to a very low value, e.g.
$m=0.001$.

Figure 7. First Row: evolution in time of the initial gaussian distribution
according to CBO dynamics (4.4); Second Row: evolution in time of the initial
Gaussian distribution according to PSO dynamics (4.3) with $m=0.8$; Third Row:
evolution produced by the PSO dynamics (4.3) with $m=0.001$.

In the case without memory, it is easy to see the convergence at each time
step of the PSO system with small inertia to the CBO, but it’s not as clear
now that the distribution we are interested in is bi-dimensional. To be able
to compare the evolution and to check how similar the distributions are at
every time step, we show in Figure 8 different plots: we look at the
particles’ mean positions (left) and their local best (center) for each time
step and, in both cases, show their standard deviation as a colored area
around the mean. It is interesting to see how the particles are moving and
where they are attracted to, especially because the pattern of the PSO
dynamics with small inertia is the same as the one of CBO. Finally, on the
right plot of Figure 8, we show how both our similarity measures, namely the
Wasserstein distance and the Kullback-Leibler divergence, are decreasing along
with the inertia parameter, validating the small inertia limit also in the
general case with memory.

Figure 8. Left: mean particles’ position with their standard deviation;
Center: mean particles’ local best with the standard deviation too; Right:
plot of the Wasserstein distance and Kullback–Leibler divergence between the
distribution obteined via CBO dynamics (4.4) and the one obtained through PSO
system (4.3) with small inertia value, i.e. $0.001$.

## References


# Detection of $\sim\text{100 days}$ periodicity in the gamma-ray light curve
of the BL Lac 4FGL 2022.7+4216

Anuvab Banerjee,1 Ajay Sharma,1 Avijit Mandal,1 Avik Kumar Das,2 Gopal
Bhatta,3 and Debanjan Bose4
1S. N. Bose National Centre for Basic Sciences, Block JD, Sector III, Salt
Lake, Kolkata 700 106, India
2Department of Physical Sciences, Indian Institute of Science Education and
Research Mohali, Knowledge City, Sector 81, SAS Nagar, Punjab 140306, India
3Institute of Nuclear Physics, Polish Academy of Sciences, PL-31342 Kraków,
Poland
4School of Astrophysics, Presidency University, 86/1 College Street, Kolkata
700073, West Bengal, India E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Study of quasi-periodic oscillations (QPO) in blazars is one of the crucial
methods for gaining insights into the workings of the central engines of
active galactic nuclei. QPOs with various characteristic time scales have been
observed in the multi-wavelength emission of blazars, ranging from the radio
to gamma-ray frequency bands. In this study, we carry out a comprehensive
variability analysis of the BL Lac object 4FGL 2022.7+4216 detected by the
Fermi-LAT, over a period of more than three years, from April 27, 2019 to
August 09, 2022. By utilizing multiple widely-used methods of time-series
analyses, we detect the presence of quasi-periodic fluctuations with a period
of $\sim$100 days with a confidence level exceeding $4\sigma$. This is the
first time such a variability feature pertaining to this source is being
reported. We propose that the observed QPO may be related to the precession of
the blazar jet with a high Lorentz factor or to the motion of a plasma blob
through the helical structure of the jet. However, for a decisive conclusion
on the physical origin of such fluctuation, further multi-wavelength
complementary observations, especially Very Long Baseline Interferometric
observations, would be required.

###### keywords:

BL Lacertae objects: individual (4FGL 2022.7+4216) – galaxies: active –
galaxies: jet – methods: observational

††pubyear: 2023††pagerange: Detection of $\sim\text{100 days}$ periodicity in
the gamma-ray light curve of the BL Lac 4FGL 2022.7+4216–Detection of
$\sim\text{100 days}$ periodicity in the gamma-ray light curve of the BL Lac
4FGL 2022.7+4216

##  1 Introduction

Blazars are a subclass of active galactic nuclei (AGN) with their relativistic
jets pointed toward our line of sight. Blazar continuum emission exhibits
rapid flux variability over a wide range of time scales, such that
multifrequency variability studies offer significant insights into the
location, size, and underlying physical processes of the emission regions.
Moreover, as the high energy $\gamma$-ray emission of blazars originate from
their jets, the exploration of $\gamma$-ray variability provides clues
regarding the jet dynamics of such systems.

Even though the blazar variability is often of stochastic nature such that the
power spectral density can be fairly represented by a power-law over a wide
range of temporal frequencies(see Bhatta & Dhital, 2020; Sobolewska et al.,
2014, and the references therein), there have been claims of the detection of
quasi-periodic oscillations (QPO) in different wavebands with a time scale
ranging from a few days to years (e.g. Rani et al., 2009; Bhatta, 2019; Gupta
et al., 2019; Bhatta, 2021). In particular, since the $\gamma$-rays in blazars
are widely believed to be originated from the highly collimated relativistic
jets, QPO detection in $\gamma$-ray is crucial to infer the jet dynamics and
particle acceleration mechanisms. It was, however, pointed out that many of
such detection last for only 2-4 cycles and the significance of those
detection is likely to be overestimated (Gupta, 2014).

Owing to the continuous monitoring of Fermi-LAT to explore the $\gamma$-ray
sky, a few significant detection of $\gamma$-ray QPOs have been, however,
reported; namely, 34.5-day QPO in PKS 2247-131 (Zhou et al., 2018a),
$\sim$47-day QPO in 3C454.3 (Sarkar et al., 2021), and $\sim$7.6-day QPO in
CTA 102 (Sarkar et al., 2020). Using a systematic search for periodicities
using a comprehensive likelihood estimation on a large sample of $\gamma$-ray
detected blazars, Peñil et al. (2020) found only 11 sources showing strong
signatures of periodicity of $>4\sigma$ significance.

The BL Lacerate object 4FGL J2202.7+4216 at redshift $z=0.069$ was detected in
a high flux state on May 1, 2019, with a daily averaged flux reaching $(1.5\pm
0.2)\times 10^{-6}~{}\text{photons cm}^{-2}\text{s}^{-1}$ by Fermi-LAT
(Garrappa & Buson, 2019). Subsequently, on August 19, 2020, Fermi-LAT has
captured the source with daily averaged flux $(2.8\pm 0.2)\times
10^{-6}~{}\text{photons cm}^{-2}\text{s}^{-1}$ (Ojha & Valverd, 2020). The
$\gamma$-ray flaring activity has coincided with the optical brightening where
R magnitude $<12$ has been registered by several optical monitoring campaigns
(Grishina & Larionov, 2020; Jankowsky & Wagner, 2020; Steineke et al., 2020).
Strong intra-night optical variability has been detected by Automatic
Telescope for Optical Monitoring (ATOM) monitoring program (Jankowsky &
Wagner, 2020). The detection of a 356 GeV photon with a strong possibility to
be associated with the source has also been reported (Garrappa & Buson, 2019).
Therefore, it can be regarded as a plausible candidate for the very high
energy (VHE) emitting source.

In this work, we present a study the Fermi/LAT gamma-ray light curve of the BL
Lac object 4FGL 2022.7+4216 detected by the Fermi-LAT spanning more than three
years. By employing multiple methods of time-series analyses, we report
presence of QPO with a characteristic timescale of $\sim$100 days. In Section
2, we provide an outline of the data acquisition and processing method used
for the Fermi/LAT telescope, including relevant details. In Section 3,
analysis and the results of time series methods, that is, LSP and WWZ and
auto-regressive method, are presented. In Section 4, the result and its
implications are discussed in the light of standard model of AGN and
conclusions are summarized.

##  2 OBSERVATIONS AND DATA REDUCTION

A study of the source $4FGL2022.7+4216$ has been done using the Fermi-LAT
data. Fermi is a space-based gamma-ray observatory onboard two instruments
(One is the Large Area Telescope (LAT) and another is Gamma-ray Burst Monitor
(GBM)). The LAT has a large effective area of > 8000 $cm^{2}$ at $\sim 1GeV$.
It is a pair conversion detector, that mainly operates in the energy range 20
MeV - 300 GeV and has a FOV of 2.4 sr Atwood et al. (2009), covering about
20$\%$ of the sky at any time, and the whole sky in every three hours.
In this study, the source was analyzed in the time domain MJD 58552-59847 and
was selected with 15∘ circular region of interest(ROI) centered at RA: 330.68
and Dec: 42.2778.
The recommended Fermi science tools ‘FERMITOOLS’ package was used to analyze
the Fermi-LAT data of the source. The data was obtained in the energy range 20
MeV - 300GeV from the fermi PASS 8 database and it was filtered using gtselect
FERMITOOLS tool with constraints evclass=128 and evtype=3 and to prevent the
source contamination with from earth limb, we set a criterion on the zenith
angle, should be less than 90∘. We filtered the data with the standard filter
${}^{\prime}\text{(DATA\\_QUAL > 0) \&\& (LAT\\_CONFIG == 1)}^{\prime}$ using
the gtmktime tool to obtain high-quality data in the good time intervals
(GTIs). The integrated livetime as a function of sky position and off-axis
angle and exposure were computed using gtltcube and gtexposure tasks
respectively. An unbinned likelihood analysis was performed using the gtlike
tool (Cash, 1979; Mattox et al., 1996) which provided the significance of each
source within the region of interest (ROI) in the form of test statistics (TS
$\sim\sigma^{2}$) (Mattox et al., 1996). The 7-day binned light curve was
obtained by integrating the source fluxes for the intervals where
$\text{TS}>25$ (above $\sim 5\sigma$ significance; Figure-1). The diffuse
$\gamma$-ray emissions of the galactic and extra-galactic were modeled using
two files : gll_iem_v07.fits and iso_P8R2_SOURCE_V3_v1.txt and adopted
instrument response functions were P8R2_SOURCE_V3  to obtain the spectrum of
the source. The iterative fitting of the light curve was done using ENRICO
software (Sanchez & Deil, 2013).

Figure 1: The Fermi-LAT light curve of the blazar 4FGL2022.7+4216 with 7 day
binning (TS > 25) in MJD 58552-59847.  Figure 2: Top: The Lomb-Scargle
periodogram corresponding to the lightcurve and the 4$\sigma$ significance
line. The period around $0.01~{}\text{days}^{-1}$ happens to be $>4\sigma$
significant. Bottom left: WWZ map of the source 4FGL2022.7+4216 light curve in
MJD 58552-59847, colors in color bar are indicating the wavelet power levels.
The decipherable power concentration around $0.01~{}\text{days}^{-1}$
indicates the plausibility of the presence of a periodicity. The time stamps
are relative to MJD = 58552. Lower right panel: Blue solid line is
representing the average wavelet power across time and the cyan solid line
marks the $4\sigma$ significance level of average WWZ power. The statistical
significance of $0.01~{}\text{days}^{-1}$ is corroborated by this method as
well.

##  3 DATA ANALYSIS & RESULTS

We adopted several quantitative tests to search for QPO in the gamma-ray light
curve of the source, namely the ‘Weighted Wavelet Z-transform (WWZ)’ method,
‘Lomb-Scargle Periodogram (LSP)’ method, and ‘REDFIT’ method. Below we
describe the results obtained by the application of these methods.

###  3.1 Lomb-Scargle Peridogram

The LSP, introduced by Lomb (Lomb, 1976) and extended later by Scargle
(Scargle, 1982), is actually a variant of the traditional discrete Fourier
transform (DFT). However, it has a distinct advantage over DFT, which is that
for uneven sampling, it reduces the effect of sampling irregularities by
iteratively fitting sinusoidal waves to the data. We compute the LSP using the
astropy LombScargle class by considering the minimum and maximum temporal
frequencies to be $1/T$ and $1/2\Delta{t}$ respectively, where $T$ is the
total observation period and $\Delta{t}$ is the time-binning. The resulting
LSP shows a prominent peak at $0.01~{}\text{days}^{-1}$.

However, it is quite well known that variability observed in blazar
lightcurves are also associated with underlying red noise, which can lead to
apparent periodic behavior for a few cycles within the low-frequency regime
(Vaughan, 2005). We, therefore, attempt to determine the statistical
significance of the periodic feature using the method proposed by
Emmanoulopoulos et al. (2013). We approximate the observed power spectra by a
power-law model and simulate 1000 light curves with the best-fit power
spectral slope and flux distribution as that of the original light curve. The
periodicity feature at $0.01~{}\text{days}^{-1}$ is found to be $>4\sigma$
significant (Top panel of Figure 2).

Furthermore, it is conventional to test the presence of periodicity using
Generalized Lomb-Scargle Periodoram
(GLSP)111https://pyastronomy.readthedocs.io/en/latest/pyTimingDoc/pyPeriodDoc/gls.html
as well which accounts for the measurement errors in the analysis. We obtain
the periodicity with the same $\sim$100 day period using GLSP as well, which
further strengthens our claim.

Statistically speaking, however, since we do not have a priori expectations of
having a particular periodicity, it is more rigorous to estimate the ‘global
significance’, which is the fraction of surrogate light curves showing a LSP
peak larger than the observed confidence level at any frequency value. In this
way, it is ensured that the peaked component is searched over a larger
frequency interval, and consequently, we are sampling from a larger population
of ‘false positive’ signals. If we do not restrict the peak to reside within a
particular frequency range, then we find that the peaked component at 100 days
is associated with $>99\%$ global confidence.

###  3.2 Weighted Wavelet Z-Transform

The standard LSP method has the limitation that it attempts to fit the
sinusoidal profile across the entire domain of observation and does not
account for the fact that the features coming from the real astrophysical
observation could be time-dependent, i.e. the amplitude and frequency could
evolve over time. Therefore, in order to characterize the periodicity features
and their evolution, the wavelet transform method turns out to be a more
suitable tool, which convolves the light curve with the time and frequency-
dependent kernel and attempts to localize the periodicity feature in time and
frequency. For the purpose of our analysis, we use the Morlet kernel Grossmann
& Morlet (1984) with the functional form:

$f[\omega(t-\tau)]=\exp[I\omega(t-\tau)-c\omega^{2}(t-\tau)^{2}]$ (1)

and the WWZ map is given by,

$W[\omega,\tau:x(t)]=\omega^{1/2}\int x(t)f^{*}[\omega(t-\tau)]dt$ (2)

We use publicly available software222https://github.com/eaydin/WWZ to estimate
the WWZ power as a function of frequency and evolution of time. In the time-
frequency plane, the color-scaled WWZ map demonstrates a decipherable
concentration of power around $0.01~{}\text{days}^{-1}$, which is apparent in
the average WWZ power as well (Bottom left panel of Figure 2). We estimate the
significance of this peak WWZ power in the same way by simulating 1000 light
curves according to the best-fit PSD model and resampling the simulated light
curves according to the source light curves. The peak significance was found
to be $>4\sigma$ here as well (Bottom right panel of Figure 2).

Figure 3: Power spectrum and the significance levels using REDFIT method. The
blue and magenta curves denote the theoretical AR(1) spectrum and the mean
simulated AR(1) spectrum respectively. The green, orange and red curves
represent 90, 95, and 99% significance levels respectively.

###  3.3 REDFIT

In the REDFIT method, the unevenly spaced time series data are fitted with a
first-order autoregressive process (AR1), avoiding interpolation and its
inevitable bias in the time domain (Schulz & Mudelsee, 2002). This method is
also used to test the significance of the flux peaks in time series against
the background of red-noise in the first-order autoregressive process. The
usage of AR1 process is justified by the autoregressive nature of the emission
flux of blazars (Schulz & Mudelsee, 2002; Kushwaha et al., 2020), where the
current emission flux is dependent on its previous flux state. We use a REDFIT
program to estimate the spectrum using LSP (Lomb-scargle periodogram) and WOSA
(Welch-overlapped-segment-averaging) procedures with the number of WOSA
segments ($n_{50}=1$); a Welch window was chosen to reduce spectral leakage.
The bias-corrected power spectrum alongside the theoretical and simulation-
generated AR(1) processes are provided in Figure-3. It shows a prominent peak
of about 100 days with a confidence level of 99$\%$ (red curve), which is the
maximum significance provided by the REDFIT program.

##  4 DISCUSSION AND CONCLUSIONS

In this work, we report the detection of 100 days periodicity pertaining to
the $\gamma$-ray detection BL Lac object 4FGL 2022.7+4216 using three
different methods of time series analysis. Below we illustrate a few physical
scenarios which can give rise to the periodic behaviour of $\gamma$-ray
lightcurve and thereby infer the plausible mechanisms operative in our context
on the basis of time scale.

  * •

In the case of binary supermassive black hole (SMBH) systems, if the secondary
black hole pierces the accretion disk of the primary black hole, QPO may be
observed as a consequence of the impact flashes (Valtonen et al., 2008). Such
a QPO of a period of $\sim$12 years has been reported earlier in the context
of OJ 287 (Valtonen et al., 2008). However, the time scale corresponding to
this process is $\sim$ years, and, therefore, the present detection of 100-day
periodicity is difficult to explain under the purview of this scenario.

  * •

Since blazar emission is dominated by jets, it is quite likely that the QPO
signature will have some bearing on the jet emission features. The precession
of the blazar jets could be induced by the interaction of the secondary source
and could manifest itself as quasi-periodicity owing to the varying Lorentz
factor (Begelman et al., 1980). Furthermore, the orientation of the jet could
also be influenced by the Lense-Thirring precession of the inner edge of the
disk. However, such mechanisms typically result in periodicity $\sim$ 1-2
years time-scale (Rieger, 2007) and fall outside of the ballpark of the
present detection. However, for blazars with jets closely aligned with our
line of sight, the detected periodicity could be significantly shorter because
of the Doppler boosting effect (Rieger, 2004). Such a mechanism has been
inferred to explain the 34.5-day QPO in the case of BL Lac PKS 2247–131 (Zhou
et al., 2018b). In our case, therefore, such a mechanism could be operative if
the jet is associated with large Lorentz factors.

  * •

The variable Doppler factor arising out of the movement of the plasma blob
along the internal helical structure of the jet could be another plausible
physical mechanism of jet-driven quasi-periodicity (Camenzind & Krockenberger,
1992). Such helical structures could come from the interaction of the jet with
the surrounding medium (Godfrey et al., 2012) or the hydrodynamic instability
effects (Hardee & Rosen, 1999). Depending upon the parameters like the pitch
angle, the viewing angle, and the Doppler boosting factor, the variability
time-scale can range from $\sim$ a few days to $\sim$ months (Rani et al.,
2009). The observed 100-day periodicity could be a result of such structural
effects of the jet.

Given the time scale of the QPO, we infer it is most likely originating from a
precessing jet with a high Lorentz factor or because of the motion of the
plasma blob along a curved jet. In the case of the one-zone leptonic model
where the plasma blob moves along the postulated helical trajectory, the time
dependence of the viewing angle will cause varying Doppler factor and
consequent intensity variation even without the intrinsic changes in jet
emission. The time dependence of the viewing reads as

$\cos{\theta_{\text{obs}}(t)}=\sin{\phi}\sin{\psi}\cos{2\pi{t}/P_{\text{obs}}}+\cos{\phi}\cos{\psi}$
(3)

where $\psi$ is the jet angle relative to our line of sight and $\phi$ is the
pitch angle of the blob. $P_{\text{obs}}$ is the observed periodicity
(Sobacchi et al., 2016). The Doppler factor varies with time as
$\delta(t)=1/(\Gamma(1-\beta\cos{\theta(t)}))$, where
$\Gamma=1/\sqrt{1-\beta^{2}}$ is the bulk Lorentz factor of the jet motion
where $\beta=\frac{v_{\text{jet}}}{t}$. The periodicity in the rest frame of
the blob is estimated by

$P_{\text{rf}}=\frac{P_{\text{obs}}}{1-\beta\cos{\psi}\cos{\phi}}.$ (4)

For typical values of $\phi=2^{\circ}$, $\psi=5^{\circ}$ and $\Gamma=8.5$, the
rest frame periodicity becomes $\sim$24 years for $P_{\text{obs}}\sim 100$
days as we have obtained in our case. During one period, the blob traverses a
distance $D=c\beta P_{\text{rf}}\cos{\phi}\sim 7$ parsec. Since the prominent
QPO signature of $\sim$100 days is observed throughout the entire domain of
$\sim$1200 days in the WWZ map, we expect roughly $\sim$12 cycles of
oscillations during this period. Therefore, the projected distance during
these 12 cycles would be estimated as $D_{p}\sim 12D\sin{\psi}=7.2$ pc. Such
parsec scale helical jets have been found in the case of several other blazar
sources as well (Vicente et al., 1996; Tateyama et al., 2002).

Roy et al. (2021) considered that the inclination angle of the jet axis
relative to the line of sight ($\psi$) could be time-dependent, which can
explain the time variation of the amplitude of oscillation. Therefore, the
amplitude variation would then be a direct offshoot of the geometric bending
structure of the jet. However, a more comprehensive investigation, including
detailed very long baseline interferometry (VLBI) monitoring needs to be
undertaken to confirm the presence of such curvature within a length scale of
$7$pc.

## Acknowledgements

This study made use of the Fermi–LAT data, obtained from the Fermi Science
Support Center (FSSC), distributed by NASA’s Goddard Space Flight Center
(GSFC). D. Bose acknowledges the support of Ramanujan Fellowship-
SB/S2/RJN-038/2017. A. Sharma and A. Mandal are grateful to S. N. Bose
National Centre for Basic Sciences under the Department of Science and
Technology (DST), Government of India, for providing the necessary support to
conduct this research.

## Data Availability

The work uses publicly available data from Fermi-LAT.

## References

  * Atwood et al. (2009) Atwood W., et al., 2009, The Astrophysical Journal, 697, 1071
  * Begelman et al. (1980) Begelman M. C., Blandford R. D., Rees M. J., 1980, Nature, 287, 307
  * Bhatta (2019) Bhatta G., 2019, MNRAS, 487, 3990
  * Bhatta (2021) Bhatta G., 2021, ApJ, 923, 7
  * Bhatta & Dhital (2020) Bhatta G., Dhital N., 2020, ApJ, 891, 120
  * Camenzind & Krockenberger (1992) Camenzind M., Krockenberger M., 1992, A&A, 255, 59
  * Cash (1979) Cash W., 1979, The Astrophysical Journal, 228, 939
  * Emmanoulopoulos et al. (2013) Emmanoulopoulos D., McHardy I. M., Papadakis I. E., 2013, MNRAS, 433, 907
  * Garrappa & Buson (2019) Garrappa S., Buson S., 2019, The Astronomer’s Telegram, 12718, 1
  * Godfrey et al. (2012) Godfrey L. E. H., et al., 2012, ApJ, 758, L27
  * Grishina & Larionov (2020) Grishina T. S., Larionov V. M., 2020, The Astronomer’s Telegram, 13930, 1
  * Grossmann & Morlet (1984) Grossmann A., Morlet J., 1984, SIAM journal on mathematical analysis, 15, 723
  * Gupta (2014) Gupta A. C., 2014, Journal of Astrophysics and Astronomy, 35, 307
  * Gupta et al. (2019) Gupta A. C., Tripathi A., Wiita P. J., Kushwaha P., Zhang Z., Bambi C., 2019, Monthly Notices of the Royal Astronomical Society, 484, 5785
  * Hardee & Rosen (1999) Hardee P. E., Rosen A., 1999, ApJ, 524, 650
  * Jankowsky & Wagner (2020) Jankowsky F., Wagner S., 2020, The Astronomer’s Telegram, 13956, 1
  * Kushwaha et al. (2020) Kushwaha P., Sarkar A., Gupta A. C., Tripathi A., Wiita P. J., 2020, Monthly Notices of the Royal Astronomical Society, 499, 653
  * Lomb (1976) Lomb N. R., 1976, Astrophysics and space science, 39, 447
  * Mattox et al. (1996) Mattox J. R., et al., 1996, The Astrophysical Journal, 461, 396
  * Ojha & Valverd (2020) Ojha R., Valverd J., 2020, The Astronomer’s Telegram, 13964, 1
  * Peñil et al. (2020) Peñil P., et al., 2020, ApJ, 896, 134
  * Rani et al. (2009) Rani B., Wiita P. J., Gupta A. C., 2009, ApJ, 696, 2170
  * Rieger (2004) Rieger F. M., 2004, ApJ, 615, L5
  * Rieger (2007) Rieger F. M., 2007, in Aschenbach B., Burwitz V., Hasinger G., Leibundgut B., eds, Relativistic Astrophysics Legacy and Cosmology - Einstein’s Legacy. p. 276, doi:10.1007/978-3-540-74713-0_63
  * Roy et al. (2021) Roy A., Sarkar A., Chatterjee A., Gupta A. C., Chitnis V., Wiita P. J., 2021, MNRAS, 510, 3641
  * Sanchez & Deil (2013) Sanchez D., Deil C., 2013, arXiv preprint arXiv:1307.4534
  * Sarkar et al. (2020) Sarkar A., Kushwaha P., Gupta A. C., Chitnis V. R., Wiita P. J., 2020, A&A, 642, A129
  * Sarkar et al. (2021) Sarkar A., Gupta A. C., Chitnis V. R., Wiita P. J., 2021, MNRAS, 501, 50
  * Scargle (1982) Scargle J. D., 1982, Astrophysical Journal, Part 1, vol. 263, Dec. 15, 1982, p. 835-853., 263, 835
  * Schulz & Mudelsee (2002) Schulz M., Mudelsee M., 2002, Computers & Geosciences, 28, 421
  * Sobacchi et al. (2016) Sobacchi E., Sormani M. C., Stamerra A., 2016, Monthly Notices of the Royal Astronomical Society, p. stw2684
  * Sobolewska et al. (2014) Sobolewska M. A., Siemiginowska A., Kelly B. C., Nalewajko K., 2014, The Astrophysical Journal, 786, 143
  * Steineke et al. (2020) Steineke R., et al., 2020, The Astronomer’s Telegram, 13958, 1
  * Tateyama et al. (2002) Tateyama C. E., Kingham K. A., Kaufmann P., de Lucena A. M. P., 2002, ApJ, 573, 496
  * Valtonen et al. (2008) Valtonen M. J., et al., 2008, Nature, 452, 851
  * Vaughan (2005) Vaughan S., 2005, A&A, 431, 391
  * Vicente et al. (1996) Vicente L., Charlot P., Sol H., 1996, A&A, 312, 727
  * Zhou et al. (2018a) Zhou J., Wang Z., Chen L., Wiita P. J., Vadakkumthani J., Morrell N., Zhang P., Zhang J., 2018a, Nature Communications, 9, 4599
  * Zhou et al. (2018b) Zhou J., Wang Z., Chen L., Wiita P. J., Vadakkumthani J., Morrell N., Zhang P., Zhang J., 2018b, Nature communications, 9, 1