# A note on Fokker-Planck equations and graphons

Fabio Coppini  Paris, France<EMAIL_ADDRESS>
###### Abstract.

Fokker-Planck equations represent a suitable description of the finite-time
behavior for a large class of particle systems as the size of the population
tends to infinity. Recently, the theory of graph limits have been introduced
in the classical mean-field framework to account for heterogeneous
interactions among particles. In many instances, such network heterogeneity is
preserved in the limit which turns from being a single Fokker-Planck equation
(also known as McKean-Vlasov) to an infinite system of non-linear partial
differential equations (PDE) coupled by means of a graphon. While appealing
from an applied viewpoint, few rigorous results exist on the graphon particle
system. This note addresses such limit system focusing on the relation between
initial conditions and interaction network: if the system initial datum and
the graphon degrees satisfy a suitable condition, a significantly simpler
representation of the solution is available.

Interesting conclusions can be drawn from this result: substantially different
graphons can give rise to exactly the same particle behavior, e.g., if the
initial condition is independent of the graphon, the particle behavior is
indistinguishable from the well-known mean-field limit as long as the graphon
has constant degree; step-kernels, a building block of graph limit theory, can
be used to approximate the graphon particle system with a finite system of
Fokker-Planck equations and an explicit rate of convergence.

2020 MSC: 35Q84, 05C82, 82B20, 60H20, 60K35, 35Q70.

Key words and phrases: Fokker-Planck equation, graphons, mean-field systems,
weakly interacting particles, McKean-Vlasov, graphon particle system,
interacting particle systems, step kernels, graph limits

##  1\. Introduction, model and literature

###  1.1. Introduction, aim of this note and organization

In the last years, the study of interacting particle systems with a non-
trivial dense network structure has been repeatedly addressed in the
mathematical community: see, e.g., [4, 10, 22, 9] for interacting diffusions,
[5, 6, 21] for applications in mean-field games and [15, 11] in the context of
dynamical systems. Depending on the setting, many results on interacting
particle systems are nowadays available [1, 3, 8, 7, 18] whenever the
underlying graph sequence is converging, in some sense depending case by case,
to a suitable object. More precisely, if the graph limit is a graphon (e.g.,
[17, 16]), then, as the size of the system tends to infinity, the finite-time
population behavior is suitably described by an infinite system of coupled
non-linear Fokker-Planck equations, the coupling between equations being made
by the graphon limit itself (see equation (1.1) for a first example).

This note addresses the graphon particle system obtained in the limit and
tries to better understand the relationship between the initial conditions and
the choice of the underlying network. Such a perspective is important not only
to better understand the limit system, but also in case one is interested in
reconstructing the network structure by looking at the particle dynamics: we
prove that, depending on the initial conditions, the class of suitable
graphons is potentially very large.

Graphons have been used as a model for many real-world networks, yet, to the
author’s knowledge, known results on graphon particle systems are mainly
limited to existence and uniqueness of solutions; we refer to Subsection 1.3
for the current literature.

A labeled graphon (we use the same notation of [16]) is a symmetric measurable
function $W$ defined on the unit-square

$\begin{split}W:[0,1]^{2}&\to[0,1]\\\ (x,y)&\mapsto W(x,y).\end{split}$

If the limit population is represented as a continuum of particles labeled by
the unit interval $[0,1]$, then $W(x,y)$ stands for the connection strength
between the particle labeled with $x$ and the one labeled with $y$. The
function $W$ thus describes the connection network underlying an infinite
population of particles.

Fix a finite time horizon $T>0$. For one-dimensional particles which are
interacting through a regular function $\Gamma$, the graphon particle system
is an infinite system of partial differential equations (PDE) coupled by means
of $W$, i.e.,

$\begin{split}\partial_{t}\mu_{t}(\theta,x)&=\frac{1}{2}\partial^{2}_{\theta}\mu_{t}(\theta,x)\\\
-&\partial_{\theta}\left[\int_{0}^{1}W(x,y)\mu_{t}(\theta,x)\int_{\mathbb{R}}\Gamma(\theta,\theta^{\prime})\mu_{t}(\text{\rm
d}\theta^{\prime},y)\text{\rm d}y\right],\quad x\in[0,1],\end{split}$ (1.1)

for $t\in[0,T]$ and where $\mu=\\{\mu_{t}(\cdot,x),t\in[0,T]\\}_{x\in[0,1]}$
is a collection of probability measures. The initial datum is given by a
probability measure $\mu_{0}\in{\mathcal{P}}(\mathbb{R}\times[0,1])$. From the
disintegration theorem ([13, §5]), we have that
$\mu_{0}(\cdot,x)\in{\mathcal{P}}(\mathbb{R})$ for almost every $x\in[0,1]$.
Namely, the probability measure

$\mu(x):=\mu(\cdot,x)=\\{\mu_{t}(\cdot,x),t\in[0,T]\\}\in{\mathcal{P}}(C([0,T],\mathbb{R}))$

represents the law of the trajectory associated to the $x$-labeled particle,
this last one being connected to the others by means of
$W(x,\cdot):[0,1]\to[0,1]$, see the following (1.2). Each trajectory is a
continuous function from $[0,T]$ to $\mathbb{R}$, i.e., an element of the
space $C([0,T],\mathbb{R})$.

To the system of non-linear Fokker-Planck equations (1.1), it is associated a
family of continuous processes $\\{\theta^{x}\\}_{x\in[0,1]}\subset
C([0,T],\mathbb{R})$, which solve

$\begin{cases}\theta^{x}_{t}=\theta^{x}_{0}+\int_{0}^{t}\int_{0}^{1}W(x,y)\int_{\mathbb{R}}\Gamma(\theta^{x}_{s},\theta)\mu_{s}(\text{\rm
d}\theta,y)\text{\rm d}y\,\text{\rm d}s+B^{x}_{t},\quad t\in[0,T],\\\
\mu_{t}(\cdot,y)={\mathcal{L}}(\theta^{y}_{t}),\quad\text{ for
}t\in[0,T]\text{ and }y\in[0,1],\end{cases}$ (1.2)

where the law of the initial condition $\theta^{x}_{0}$, denoted by
${\mathcal{L}}(\theta^{x}_{0})$, is thus given by $\mu_{0}(x)$. The family
$\\{B^{x}\\}_{x\in[0,1]}$ is composed of independent and identically
distributed (IID) Brownian motions, independent of the initial conditions as
well.

The link between (1.1) and (1.2) is given by the fact that
$\mu_{t}(x)={\mathcal{L}}(\theta^{x}_{t})$ for every $t\in[0,T]$ and
$x\in[0,1]$. We refer to Proposition 1.3 for a precise statement which
includes more general systems, including the degenerate case ($\sigma\equiv
0$).

System (1.1) or, equivalently, the family of non-linear processes (1.2), have
been proposed as a limit description in the literature [18, 22, 5, 1, 3], yet
very little is known on their mathematical properties. We thus aim at
mathematically addressing (1.1), firstly by showing the strong link with the
graphon theory (including unlabeled graphons, see Proposition 1.5), secondly
by proving that a simpler representation for both (1.1) and (1.2) exists,
whenever the initial conditions and the underlying graphon satisfy a suitable
condition, see Theorem 2.1 and, in particular, Corollary 2.3 and Proposition
2.6.

The note is organized as follows: the next subsection presents the general
model associated to (1.2), together with known results from the existing
literature. Related works are discussed at the end of this first section.

Section 2 provides the main result, Theorem 2.1, as well as two corollaries
and a relevant application to step kernels, see Proposition 2.6. Notably,
Corollary 2.4 addresses the classical mean-field scenario showing that, if
every particle has the same initial law, then for every graphon with constant
degree, i.e., such that $d(\cdot)=\int_{0}^{1}W(\cdot,y)\text{\rm d}y$ is
constant, the dynamics is mean-field.

Section 3 contains the mathematical set-up and the proof of the main result.

###  1.2. The model and some known results

We consider a class of models slightly more general than (1.2): fix
$\mu_{0}\in{\mathcal{P}}(\mathbb{R}\times[0,1])$ and let
$\\{\theta^{x}\\}_{x\in[0,1]}$ be the family solving the $\infty$-dimensional
coupled system:

$\begin{cases}\theta^{x}_{t}=\theta^{x}_{0}+\int_{0}^{t}F(\theta^{x}_{s})\text{\rm
d}s+\int_{0}^{t}\int_{0}^{1}W(x,y)\int_{\mathbb{R}}\Gamma(\theta^{x}_{s},\theta)\mu_{s}(\text{\rm
d}\theta,y)\text{\rm d}y\,\text{\rm
d}s+\int_{0}^{t}\sigma(\theta^{x}_{s})\text{\rm d}B^{x}_{s}\\\
\mu_{t}(\cdot,y)={\mathcal{L}}(\theta^{y}_{t}),\quad\text{ for
}t\in[0,T]\text{ and }y\in[0,1],\end{cases}$ (1.3)

where $F$, $\Gamma$ and $\sigma$ are 1-Lipschitz functions bounded by 1, and
$\\{B^{x}\\}_{x\in[0,1]}$ is a family of IID Brownian motions on $\mathbb{R}$.
For every $x\in[0,1]$, the initial condition $\theta^{x}_{0}$ is a random
variable with law given by ${\mathcal{L}}(\theta^{x}_{0})=\mu_{0}(x)$,
independent of the other initial conditions and of the Brownian motions. The
probability measure induced by the initial conditions and the family of
Brownian motions is denoted by ${\mathbf{P}}$ and the corresponding
expectation by ${\mathbf{E}}$.

We work under the following assumptions.

######  Hypothesis 1.1.

We assume that

  1. (1)

(measurability) The map $[0,1]\ni
x\mapsto\mu_{0}(x)\in{\mathcal{P}}(\mathbb{R})$ is measurable;

  2. (2)

(moment condition) ${\mathbf{E}}[|\theta^{x}_{0}(0)|^{2}]<\infty$ for every
$x\in[0,1]$.

######  Remark 1.2.

Measurability with respect to $x$ is somehow necessary since we are working
with measurable functions $W$ on the unit interval $[0,1]$; we will explicitly
use it in the following relation (2.1).

The finite second moment condition is related to the use of the 2-Wasserstein
distance between probability measures (see the Section 3 and, in particular,
(3.1)). It is obviously possible to work with initial conditions with finite
$p$-moment for every $p\;\geqslant\;1$, and the choice $p=2$ is purely
arbitrary.

We stick to the one-dimensional setting, i.e., $\theta^{x}$ taking values in
$\mathbb{R}$; however, all the proofs presented below are easily extendable to
any finite dimension.

Existence and uniqueness for system (1.3) are known.

######  Proposition 1.3.

Under the measurability assumption and the moment condition in Hypothesis 1.1,
there exists a unique pathwise solution to system (1.3).

Moreover, if we denote by $\mu(x)$ the law of $\theta^{x}$ for each
$x\in[0,1]$, then the map $[0,1]\ni
x\mapsto\mu(x)\in{\mathcal{P}}(C([0,T],\mathbb{R}))$ is measurable and
$\mu(x)$ weakly solves

$\begin{split}\partial_{t}\mu_{t}(\theta,x)\,&=\,\frac{1}{2}\partial^{2}_{\theta}\left({\sigma^{2}(\theta)}\mu_{t}(\theta,x)\right)-\partial_{\theta}\left(\mu_{t}(\theta,x)F(\theta)\right)\\\
-&\partial_{\theta}\left[\int_{0}^{1}W(x,y)\mu_{t}(\theta,x)\,\int_{\mathbb{R}}\Gamma(\theta,\theta^{\prime})\mu_{t}(\text{\rm
d}\theta^{\prime},y)\text{\rm d}y\right],\end{split}$ (1.4)

for every $x\in[0,1]$.

We refer to [1, Proposition 2.1] and [18, Proposition 2.4] for two standard
similar proofs. We point out that $W$ needs only to be bounded and not
necessarily with values in $[0,1]$, i.e, existence and uniqueness hold for
every kernel $W$, as in the notation below.

Let ${\mathcal{W}}:=\\{W:[0,1]^{2}\to\mathbb{R}\text{ bounded, symmetric and
measurable}\\}$ be the space of kernels111we always consider two kernels to be
equal if and only if they differ on a subset of Lebesgue measure zero. We
follow closely the notation in [16].. The cut-norm of $W\in{\mathcal{W}}$ is
defined as

$\left\lVert
W\right\rVert_{\square}:=\max_{S,T\subset[0,1]}\left|\int_{S\times
T}W(x,y)\text{\rm d}x\text{\rm d}y\right|$ (1.5)

where the maximum is taken over all measurable subsets $S$ and $T$ of $I$. Let
${\mathcal{W}}_{0}:=\\{W\in{\mathcal{W}}:0\;\leqslant\;W\;\leqslant\;1\\}$ be
the space of labeled graphons: for $W,\,V\in{\mathcal{W}}_{0}$ their cut-
distance is defined by

$\delta_{\square}(W,V):=\,\min_{\varphi\in S_{[0,1]}}\,\left\lVert
W-V^{\varphi}\right\rVert_{\square},$ (1.6)

where the minimum ranges over $S_{[0,1]}$ the space of invertible measure
preserving maps from $[0,1]$ into itself and where
$V^{\varphi}(x,y):=V(\varphi(x),\varphi(y))$ for $x,y\in[0,1]$.

The cut-distance $\delta_{\square}$ is a pseudometric on ${\mathcal{W}}_{0}$
since it can be zero between two different labeled graphons. If we identify
all labeled graphons with cut-distance zero, we obtain the space of graphons
$\widetilde{{\mathcal{W}}}_{0}:={\mathcal{W}}_{0}/\delta_{\square}$. A well-
known result of graph limits theory says that
$(\widetilde{{\mathcal{W}}}_{0},\delta_{\square})$ is a compact metric space
[16, Theorem 9.23].

Proposition 1.3 can thus be restated by saying that for every kernel
$W\in{\mathcal{W}}$, there exists a unique solution $\mu^{W}$ to (1.4). It is
thus natural to ask whether the application $W\mapsto\mu^{W}$ is continuous
(e.g., in the topology of the weak-convergence) with respect to
$\left\lVert\cdot\right\rVert_{\square}$. This point is discussed in the next
remark.

######  Remark 1.4.

Under suitable assumptions on the coefficients, see [3, Proposition 3.3] but
also [1, Theorem 2.1] and [14, Lemma 2.7], it is possible to show that for two
solutions $\mu^{W}$ and $\mu^{V}$ associated to $W$ and $V$ in ${\mathcal{W}}$
respectively, it holds that

$D^{2}_{T}(\mu^{W},\mu^{V})\;\leqslant\;C\left\lVert
W-V\right\rVert_{\square},\quad\text{ for some }C>0,$ (1.7)

where $D_{T}$ is some distance on probability measures metricizing the weak-
convergence. Equation (1.7) proves that the application $W\mapsto\mu^{W}$ is
(Hölder-)continuous and, thus, that similar graphons – in terms of the cut-
norm – leads to similar particle behaviors. However, it does not say anything
whether two different graphons (possibly in $\delta_{\square}$-distance) can
lead to close behaviors. As Theorem 2.1 and its corollaries show, there are
relevant examples where this happens and where it is possible to prove that
the particle dynamics coincide.

For every $W\in{\mathcal{W}}_{0}$, Proposition 1.3 shows that the infinite
system (1.4) admits a unique solution, yet it does not address the unlabeled
class of $W$ in $\widetilde{{\mathcal{W}}}_{0}$. We explicit such relation in
the next proposition.

######  Proposition 1.5.

Let $U$ be a uniform random variable on $[0,1]$ and $\theta^{U}$ the non-
linear process in $\\{\theta^{x}\\}_{x\in[0,1]}$ with random label $U$. Then,
the law of $\theta^{U}$ is given by
$\widetilde{\mu}=\int_{0}^{1}\mu(x)\text{\rm d}x$ and it is independent of the
class of $W$ in $\widetilde{{\mathcal{W}}}_{0}$.

###### Proof.

This result was firstly stated in [3, Proposition 2.1] in the case of
interacting oscillators. In case the particles are living in $\mathbb{R}$ the
proof does not change. Indeed, let $\varphi$ be a measure preserving map from
$[0,1]$ to itself, we observe that $\mu(\varphi(x))$ solves:

$\begin{split}\partial_{t}\mu_{t}(\theta,\varphi(x))\,&=\,\frac{1}{2}\partial^{2}_{\theta}\left({\sigma^{2}(\theta)}\mu_{t}(\theta,\varphi(x))\right)-\partial_{\theta}\left(\mu_{t}(\theta,\varphi(x))F(\theta)\right)\\\
-&\partial_{\theta}\left[\int_{0}^{1}W(\varphi(x),y)\mu_{t}(\theta,\varphi(x))\,\int_{\mathbb{R}}\Gamma(\theta,\theta^{\prime})\mu_{t}(\text{\rm
d}\theta^{\prime},y)\text{\rm d}y\right],\end{split}$ (1.8)

where the last term is equal to

$\partial_{\theta}\left[\int_{0}^{1}V(x,y)\mu_{t}(\theta,\varphi(x))\,\int_{\mathbb{R}}\Gamma(\theta,\theta^{\prime})\mu_{t}(\text{\rm
d}\theta^{\prime},\varphi(y))\text{\rm d}y\right]$

with $V(x,y)=W(\varphi(x),\varphi(y))$. Thus
$\\{\mu(\varphi(x))\\}_{x\in[0,1]}$ solves the same system of
$\\{\nu(x)\\}_{x\in[0,1]}$ solution to (1.4) with $V\in{\mathcal{W}}_{0}$ and
initial condition $\\{\mu_{0}(\varphi(x))\\}_{x\in[0,1]}$. Clearly
$\widetilde{\mu}=\int_{0}^{1}\mu(x)\text{\rm
d}x=\int_{0}^{1}\mu(\varphi(x))\text{\rm d}x=\int_{0}^{1}\nu(x)\text{\rm
d}x=\widetilde{\nu}$. ∎

Proposition 1.5 makes clear that, as for a graphon the relevant information
does not depend on the chosen labeling, the same holds true for the marginal
law of the graphon particle system: it does not change under relabeling of the
particles. With the exception of [3] for interacting oscillators, this result
has never been stated in the literature.

######  Remark 1.6.

Observe that two different initial conditions
$\mu^{1}_{0},\mu^{2}_{0}\in{\mathcal{P}}(\mathbb{R}\times[0,1])$ but such that
$\int_{0}^{1}\mu^{1}_{0}(x)\text{\rm d}x=\int_{0}^{1}\mu^{2}_{0}(x)\text{\rm
d}x$, lead to two different solutions $\mu^{1}$ and $\mu^{2}$ to equation
(1.4). Following the notation of Proposition 1.5, the corresponding marginals
$\widetilde{\mu}^{i}=\int_{0}^{1}\mu^{i}(x)\text{\rm d}x$, for $i=1,2$, are in
general not equal.

In particular, by choosing
$\mu^{2}_{0}\equiv\int_{0}^{1}\mu^{1}_{0}(x)\text{\rm d}x$, one can derive the
following statement: the behavior of a uniform random particle can be
described by two different solutions, i.e., $\widetilde{\mu}^{2}$ in case only
the marginal law of the initial conditions is known, or $\widetilde{\mu}^{1}$
in case the joint law is known. Differently stated, one cannot interchange
expectation related to $U$ and the dynamics semigroup.

The way in which the dynamics mixes the initial conditions is an interesting
challenge for further studies.

###  1.3. Related works

The graphon framework [17, 16] has been introduced in the theory of particle
systems [7, 19, 20, 22] as an useful ingredient to construct inhomogeneous
random graph sequences with nice statistical properties (edge independence,
graph homomorphism, etc.). Most of the known literature on particle systems
has been focused on the converging properties of particle systems on such
random graph sequences, see, e.g., [1, 3, 18, 22]. They show that some of the
classical mean-field arguments (as propagation of chaos [23]) can be extended
to deal with random graphs and to include labeled graphons in the limit
description. Nevertheless, only little attention has been put into the study
of the limit object (1.4) which remains, to the author’s knowledge, rather
unknown.

The recent work [3] shows a direct connection between particle systems and the
graphon theory: under suitable hypothesis on (1.2), there exists a Hölder-
continuous mapping between the space of graphons
$(\widetilde{{\mathcal{W}}}_{0},\delta_{\square})$ and $\widetilde{\mu}$ as in
Proposition 1.5. This allows to consider general graph sequences, as
exchangeable random graphs [12], which leads to a (possibly) random graphon
$W$ in the limit. However, despite the rather understood convergence
properties, only a few insights are available on the limit particle system, we
refer to the examples in Subsection 2.3 of [3].

Although interesting on itself, the study of the limit Fokker-Planck equation
does not necessarily provide a suitable description of particle systems on
diverging time scales, as already raised in [9, 8] for the mean-field case.
With the exception of dissipative dynamics [2], a substantial understanding of
the phase space of (1.4) is needed to study the long-time behavior of finite
particle systems on graphs as shown in [8]. It is thus important to better
understand system (1.4) in order to address the finite particle system
counterpart.

##  2\. Main result, examples and applications

We suppose that the initial condition
$\mu_{0}\in{\mathcal{P}}(\mathbb{R}\times[0,1])$ and the labeled graphon
$W\in{\mathcal{W}}_{0}$ are fixed.

Before stating the main result, we introduce an equivalence relation on the
unit interval and give the definition of degree in $W$. For $x$ and $y$ in
$[0,1]$, define the relation $\sim$ by

$x\sim y\quad\text{if and only if}\quad\mu_{0}(x)=\mu_{0}(y),$ (2.1)

i.e., we identifies two labels whenever the initial conditions of the
corresponding particles have the same law. Denote by $J$ the quotient space
$[0,1]/$$\scriptstyle\sim$ and, for $\bar{x}\in J$, denote its orbit by

$[\bar{x}]:=\left\\{x\in[0,1]:x\sim\bar{x}\right\\}.$

For some $x\in[0,1]$, we denote the degree of $x$ by
$d(x)=\int_{0}^{1}W(x,y)\text{\rm d}y$. More generally, we define the degree
of $x$ with respect to a (measurable) subset $A\subset[0,1]$ by

$d_{A}(x)=\int_{A}W(x,y)\text{\rm d}y.$ (2.2)

Observe that under the measurability assumption (1.1), $[\bar{y}]$ is a
measurable subset of $[0,1]$ for every $\bar{y}\in J$; with a slight abuse of
notation, we denote $d_{\bar{y}}(\cdot):=d_{[\bar{y}]}(\cdot)$.

###  2.1. Main result

We have the following theorem.

######  Theorem 2.1.

Assume Hypothesis 1.1. Suppose that, for every $\bar{x}$ and $\bar{y}$ in $J$
it holds that

$d_{\bar{y}}(x)=\int_{[\bar{y}]}W(x,y)\text{\rm
d}y=\int_{[\bar{y}]}W(\bar{x},y)\text{\rm d}y=d_{\bar{y}}(\bar{x})\quad\text{
for all }x\in[\bar{x}].$ (2.3)

Then for every $\bar{x}\in J$, $\mu(x)=\mu(\bar{x})$, for all $x\in[\bar{x}]$.

Condition (2.3) is requiring that, if we partition the interval $[0,1]$ in $J$
subsets $\\{[\bar{x}]\\}_{\bar{x}\in J}$, then the density of neighbors with
respect to any of these subsets, i.e., $d_{\bar{y}}(\cdot)$ for some
$\bar{y}\in J$, is piecewise constant on $[0,1]$ and completely characterized
by the values on $J$. See Figure 1 (B) and Figure 2 (A) for graphons that
satisfy (2.3).

Theorem 2.1 basically states that, if one can group the particles so that
within each group they have the same initial condition (in law) and are equi-
connected to the other groups in the sense of (2.3), then the solution $\mu$
to (1.4) is constant on these groups.

We observe that this result does not depend on the class of $W$ as unlabeled
graphon. As a consequence, one can relabel the particles such that $\mu$ is
piecewise constant as a function of $x\in[0,1]$.

######  Remark 2.2.

It is possible to weaken (2.3) by considering a larger $J$ based on vertices
which have not only the same initial law, but also the same degree in $W$. We
refrain from increasing the complexity in the construction of $J$ to keep the
main ideas as clear as possible. This aspect is briefly discussed in the
application to step kernels below.

A direct consequence of Theorem 2.1 is that the information in
$\\{\mu(x)\\}_{x\in[0,1]}$ is contained in the eventually smaller object
$\\{\mu(\bar{x})\\}_{\bar{x}\in J}$.

######  Corollary 2.3.

Under the assumptions of Theorem 2.1. The infinite system (1.4) is suitably
described by the (possible finite) system of coupled partial differential
equations

$\begin{split}\partial_{t}\nu_{t}(\theta,\bar{x})\,&=\,\frac{1}{2}\partial^{2}_{\theta}\left({\sigma^{2}(\theta)}\nu_{t}(\theta,\bar{x})\right)-\partial_{\theta}\left(\nu_{t}(\theta,\bar{x})F(\theta)\right)\\\
-&\partial_{\theta}\left[\int_{J}W(\bar{x},\bar{y})\nu_{t}(\theta,\bar{x})\,\int_{\mathbb{R}}\Gamma(\theta,\theta^{\prime})\nu_{t}(\text{\rm
d}\theta^{\prime},\bar{y})\text{\rm d}\bar{y}\right],\quad\bar{x}\in
J,\end{split}$ (2.4)

with $\nu_{0}=\mu_{0}$. Notably, $\mu(x)=\nu(\bar{x})$ for every $\bar{x}\in
J$ and $x\in[\bar{x}]$.

###### Proof.

Existence and uniqueness for system (2.4) directly follow from Proposition
1.3. Moreover, system (2.4) is written in closed form and can be thus solved
independently of (1.4). ∎

Graphons have proven to be an important tool for establishing the convergence
of particle systems on graph sequences. Notably, it is now possible to show
that interacting particle systems on apriori different graph sequences have
the same asymptotic behavior whenever the limit of these sequences coincides
as unlabeled graphon. However, different unlabeled graphons can lead to the
same particle behavior: in these cases, the limit (1.4) tends to be a formal
object instead of giving a precise description. Theorem 2.1 and Corollary 2.3
precise that both network structure and initial conditions are necessary to
understand the system evolution.

###  2.2. The case of uniform initial conditions

We now turn to an interesting application when $|J|=1$, i.e., when the law of
the initial condition is constant with respect to $x\in[0,1]$. As already
expressed in Remark 1.6, this covers the case of the behavior of a uniform
random particle in (1.4) when only the marginal law of the initial condition,
i.e., $\int_{0}^{1}\mu_{0}(x)\text{\rm d}x$, is available.

It turns out that condition (2.3) boils down to a constant degree assumption
on $W$.

(a) Constant graphon

(b) Disconnected graphon

(c) Cayley graphon

Figure 1. Examples of graphons that satisfy (2.5) (constant degree). Observe
that the graphon in (B) is composed of two connected components, the smaller
one being more densely connected (darker) with respect to the other one so to
satisfy (2.5). Finally, observe that (A), (B) and (C) are different both as
labeled and unlabeled graphons, i.e., different in
$\delta_{\square}$-distance.

######  Corollary 2.4.

Under the assumptions of Theorem 2.1. Suppose that
$\mu_{0}(x)=\mu_{0}\in{\mathcal{P}}(\mathbb{R})$ for every $x\in[0,1]$. Then
there exists $p\in[0,1]$ such that condition (2.3) is equivalent to

$p=\int_{0}^{1}W(x,y)\text{\rm d}y,\quad\text{ for all }x\in[0,1].$ (2.5)

In particular, $\\{\mu(x)\\}_{x\in[0,1]}$ is label independent, i.e.,
$\mu(\cdot)\equiv\mu\in{\mathcal{P}}(C([0,T],\mathbb{R}))$, and $\mu$ solves
the classical McKean-Vlasov equation

$\partial_{t}\nu_{t}(\theta)\,=\,\frac{1}{2}\partial^{2}_{\theta}\left({\sigma^{2}(\theta)}\nu_{t}(\theta)\right)-\partial_{\theta}\left(\nu_{t}(\theta)F(\theta)\right)-\partial_{\theta}\left[\nu_{t}(\theta)\,p\int_{\mathbb{R}}\Gamma(\theta,\theta^{\prime})\nu_{t}(\text{\rm
d}\theta^{\prime})\right],$ (2.6)

with initial condition $\mu_{0}$.

###### Proof.

The hypothesis on $\mu_{0}$ forces $|J|=1$ and thus condition (2.3) becomes
(2.5). Corollary 2.3 yields the result. ∎

The constant degree assumption on $W$ is satisfied by a large class of non-
trivial graphons as the examples shown in Figure 1.

######  Remark 2.5.

Observe that Corollary 2.4 can be derived from the results in [10] combined
with the ones on particle systems on graphons, e.g., [18, 1]. Indeed, in [10]
it is shown that, for a particle system defined on a graph sequence and with
IID initial conditions, a sufficient condition to obtain the mean-field limit
(2.6) is that each vertex in the (renormalized) graph sequence has the same
asymptotic degree density. If one chooses such sequence to converge to a
graphon with constant degree, then it satisfies both the hypothesis in [10,
Theorem 1.1] and in, e.g., [1, Theorem 4.1]. Combining [10, Theorem 1.1] and
[1, Theorem 4.1], we obtain that the same finite particle system converges to
(2.6) but also to (1.4). As a consequence, the solutions of (1.4) and (2.6)
must be the same, giving an undirect proof of Corollary 2.4.

We step to a representative example whenever the graphon is a piecewise
constant function, as the one in Figure 2.

###  2.3. Finite approximation through step-kernels

Step kernels ([16, §7.1]) represent one of the fundamental blocks of graphon
theory since they provide a way to approximate any $W\in{\mathcal{W}}_{0}$
with constant-wise functions. A function $W\in{\mathcal{W}}$ is a step kernel
if there is a partition ${\mathcal{P}}=\\{S_{i}\\}_{i=1,\dots,k}$ of $[0,1]$
into measurable sets such that $W$ is constant on every product set
$S_{i}\times S_{j}$. We use the following notation

$W_{{\mathcal{P}}}(x,y)=\sum_{i,j=1}^{k}w_{ij}\,\mathbf{1}_{S_{i}\times
S_{j}}(x,y),\quad\text{ for }x,y\in[0,1],$ (2.7)

where $\\{w_{ij}\\}_{i,j=1,\dots,k}$ are bounded real numbers. See Figure 2
for an example with equidistant partition.

For simplicity, we suppose that $\mu_{0}$ is label independent, i.e.,
$\mu_{0}(x)=\mu_{0}\in{\mathcal{P}}(\mathbb{R})$. However, observe that the
strategy presented here can be applied by suitably approximating the initial
condition with a constant-wise initial datum: indeed, the continuity of the
solution $\mu$ to (1.4) with respect to the initial condition is a classical
result [23].

For a step kernel, condition (2.3) is clearly satisfied by taking
$J=\\{x_{i}\\}_{i=1,\dots,k}$ with $x_{i}\in S_{i}$ for every $i=1,\dots,k$,
thus refining the strategy presented before Theorem 2.1, as anticipated in
Remark 2.2.

Applying Theorem 2.1 and Corollary 2.3, we have the following representation
for the graphon particle system (1.4) on the step graphon $W_{\mathcal{P}}$.

######  Proposition 2.6.

Suppose that $\mu_{0}\in{\mathcal{P}}(\mathbb{R})$ and let $W$ be the step
kernel in (2.7). Then the graphon particle system (1.4) is suitably described
by the finite collection of probability measures $\\{\mu^{i}\\}_{i=1,\dots,k}$
which solve

$\partial_{t}\mu^{i}_{t}(\theta)\,=\,\frac{1}{2}\partial^{2}_{\theta}\left({\sigma^{2}(\theta)}\mu^{i}_{t}(\theta)\right)-\partial_{\theta}\left(\mu^{i}_{t}(\theta)F(\theta)\right)-\partial_{\theta}\left[\sum_{j=1}^{k}w_{ij}\mu^{i}_{t}(\theta)\,\int_{\mathbb{R}}\Gamma(\theta,\theta^{\prime})\mu^{j}_{t}(\text{\rm
d}\theta^{\prime})\right],$ (2.8)

and where $\mu^{i}_{0}=\mu_{0}$, for $i=1,\dots,k$.

As shown in Proposition 2.6, the representation given in Corollary 2.3 becomes
natural in the case of step kernels. While Proposition 2.6 could be derived
with direct computations222A step kernel can be seen as a Stochastic Block
Model, for which the representation (2.8) is somehow known., Theorem 2.1 goes
a step further: if one replaces every single constant block $S_{i}\times
S_{j}$ in (2.7) by a (suitable scaled) graphon with constant degree, equation
(2.8) does not change.

We also observe that combining Proposition 2.6 with the continuity estimates
as in Remark 1.4, allows to approximate the graphon particle system (1.4) with
a finite number of coupled Fokker-Planck equations. Since it is beyond the
scope of this note, we do not purse such analysis, yet we make this point a
bit more precise in the next remark.

######  Remark 2.7.

By using the Weak Regularity Lemma [16, Corollary 9.13], one can approximate
every graphon by a step kernel with an explicit control on the error. Assuming
the continuity of (1.4) with respect to $W$, recall Remark 1.4, Proposition
2.6 opens the way to approximate the infinite graphon particle system (1.4) by
using a finite (thus apriori numerically solvable) system of coupled Fokker-
Planck equations, with precise bounds on the error.

(a) Step kernel

(b) Scale free graphon

Figure 2. Suppose that $J=\\{x_{0},x_{1},x_{2}\\}$ and that $[x_{0}]=[0,1/3)$,
$[x_{1}]=[1/3,2/3)$ and $[x_{2}]=[2/3,1]$. Then, the step kernel in the figure
(A) satisfies (2.3). It also represents a way to approximate the scale-free
graphon (B). The grayscale stands for different values in $[0,1]$.

##  3\. Proof of the main result

###  3.1. Distance between probability measures

For two probability measures
$\bar{\mu},\,\bar{\nu}\in{\mathcal{P}}(C([0,T],\mathbb{R}))$, define their
2-Wasserstein distance as

$D_{T}(\bar{\mu},\bar{\nu})=\inf_{X,Y}\left\\{{\mathbf{E}}\left[\sup_{t\in[0,T]}\left|X_{t}-Y_{t}\right|^{2}\right]:{\mathcal{L}}(X)=\bar{\mu},\,{\mathcal{L}}(Y)=\bar{\nu}\right\\}^{1/2}$
(3.1)

where the infimum is taken on all random variables $X$ and $Y$ with values in
$C([0,T],\mathbb{R})$ and law ${\mathcal{L}}$ equal to $\bar{\mu}$ and
$\bar{\nu}$ respectively. From (3.1) we obtain that for every $s\in[0,T]$ and
for every bounded 1-Lipschitz function $f$

$\left|\int_{\mathbb{R}}f(\theta)\,\bar{\mu}_{s}(\text{\rm
d}\theta)-\int_{\mathbb{R}}f(\theta)\,\bar{\nu}_{s}(\text{\rm
d}\theta)\right|=\left|\int_{\mathbb{R}}f(\theta)\,\left[\bar{\mu}_{s}(\text{\rm
d}\theta)-\bar{\nu}_{s}(\text{\rm
d}\theta)\right]\right|\;\leqslant\;D_{s}(\bar{\mu},\bar{\nu}).$ (3.2)

Observe that (3.1) also makes sense with $T=0$ and $C([0,T],\mathbb{R})$
replaced by $\mathbb{R}$.

###  3.2. Proof of Theorem 2.1

Recall the definition of the 2-Wasserstein distance $D_{T}$ in (3.1), we aim
at showing that

$\max_{\begin{subarray}{c}\bar{x}\in J\\\
x\in[\bar{x}]\end{subarray}}D_{T}(\mu(\bar{x}),\mu(x))=0.$ (3.3)

We first assume $F\equiv 0$ and $\sigma\equiv 1$.

Let $\bar{x}\in J$ and $x\in[\bar{x}]$, without loss of generality, we can
suppose that the associated realizations of the Brownian motion are the same.
Using the fact that $x\in[\bar{x}]$, the initial conditions cancel and we have
that

$\begin{split}\theta^{x}_{t}-\theta^{\bar{x}}_{t}=&\int_{0}^{t}\int_{0}^{1}W(x,y)\int_{\mathbb{R}}\Gamma(\theta^{x}_{s},\theta)\mu(\text{\rm
d}\theta,y)\text{\rm d}y\,\text{\rm d}s\\\
&-\int_{0}^{t}\int_{0}^{1}W(\bar{x},y)\int_{\mathbb{R}}\Gamma(\theta^{\bar{x}}_{s},\theta)\mu(\text{\rm
d}\theta,y)\text{\rm d}y\,\text{\rm d}s,\end{split}$ (3.4)

In particular, this can be rewritten as

$\begin{split}&\int_{0}^{t}\int_{0}^{1}W(x,y)\int_{\mathbb{R}}\left[\Gamma(\theta^{x}_{s},\theta)-\Gamma(\theta^{\bar{x}}_{s},\theta)\right]\mu(\text{\rm
d}\theta,y)\text{\rm d}y\,\text{\rm d}s\\\
&+\int_{0}^{t}\int_{0}^{1}\left[W(x,y)-W(\bar{x},y)\right]\int_{\mathbb{R}}\Gamma(\theta^{\bar{x}}_{s},\theta)\mu(\text{\rm
d}\theta,y)\text{\rm d}y\,\text{\rm d}s.\end{split}$ (3.5)

We now use hypothesis (2.3) and add the following term in the previous
equation

$0=\int_{J}\left[\int_{[\bar{y}]}\left[W(x,y)-W(\bar{x},y)\right]\int_{\mathbb{R}}\Gamma(\theta^{\bar{x}}_{s},\theta)\mu_{s}(\text{\rm
d}\theta,\bar{y})\text{\rm d}y\right]\text{\rm d}\bar{y},$ (3.6)

where the integrals are sums whenever $J$ or $[\bar{y}]$ are countable (recall
that $[\bar{y}]$ is a measurable subset of $[0,1]$ because of Hypothesis
(1.1)).

By taking the squares and using $(a+b)^{2}\;\leqslant\;2(a^{2}+b^{2})$, this
leads to

$\begin{split}&\left|\theta^{x}_{t}-\theta^{\bar{x}}_{t}\right|^{2}\;\leqslant\;2T^{2}\int_{0}^{t}\left|\theta^{x}_{s}-\theta^{\bar{x}}_{s}\right|^{2}\text{\rm
d}s+\\\
&+2T^{2}\int_{0}^{t}\int_{J}\left|\int_{[\bar{y}]}[W(x,y)-W(\bar{x},y)]\int_{\mathbb{R}}\Gamma(\theta^{\bar{x}}_{s},\theta)[\mu_{s}(\text{\rm
d}\theta,y)-\mu_{s}(\text{\rm d}\theta,\bar{y})]\text{\rm
d}y\right|^{2}\text{\rm d}\bar{y}\,\text{\rm d}s,\end{split}$ (3.7)

where we have used Cauchy-Schwartz inequality as well as the fact that
$\Gamma$ is 1-Lipschitz. Applying Cauchy-Schwartz again, and using (3.2), we
are left with

$\begin{split}\left|\theta^{x}_{t}-\theta^{\bar{x}}_{t}\right|^{2}\;\leqslant\;2T^{2}\int_{0}^{t}\left|\theta^{x}_{s}-\theta^{\bar{x}}_{s}\right|^{2}\text{\rm
d}s+4T^{2}\int_{0}^{t}\int_{J}\left[\int_{[\bar{y}]}D^{2}_{s}(\mu(y),\mu(\bar{y}))\text{\rm
d}y\right]\,\text{\rm d}\bar{y}\,\text{\rm d}s,\end{split}$ (3.8)

Thus, taking the supremum with respect to the time and the expected value
${\mathbf{E}}$ this leads to

$\begin{split}{\mathbf{E}}\left[\sup_{t\in[0,T]}\left|\theta^{x}_{t}-\theta^{\bar{x}}_{t}\right|^{2}\right]\;\leqslant\;&2T^{2}\int_{0}^{T}{\mathbf{E}}\left[\sup_{u\in[0,s]}\left|\theta^{x}_{u}-\theta^{\bar{x}}_{u}\right|^{2}\right]\text{\rm
d}s+\\\
&+4T^{2}\int_{0}^{T}\int_{J}\left[\int_{[\bar{y}]}D^{2}_{s}(\mu(y),\mu(\bar{y}))\text{\rm
d}y\right]\,\text{\rm d}\bar{y}\,\text{\rm d}s.\end{split}$ (3.9)

Finally, we can take the maximum with respect to $\bar{x}\in J$ and
$x\in[\bar{x}]$ and use the characterization of (3.1), to obtain

$\max_{\begin{subarray}{c}\bar{x}\in J\\\
x\in[\bar{x}]\end{subarray}}D^{2}_{T}(\mu(\bar{x}),\mu(x))\;\leqslant\;6T^{2}\int_{0}^{T}\max_{\begin{subarray}{c}\bar{x}\in
J\\\ x\in[\bar{x}]\end{subarray}}D^{2}_{s}(\mu(\bar{x}),\mu(x))\text{\rm d}s,$
(3.10)

which implies (3.3).

Whenever $F$ is not zero, the proof is basically the same: using the Lipschitz
properties of $F$, a term equal to
$\left|\theta^{x}_{s}-\theta^{\bar{x}}_{s}\right|^{2}$ appears, thus adding
the constant factor of $2T^{2}$ to the final equation (3.9).

If $\sigma$ is not constant, then we may use the Burkholder-Davis-Gundy
inequality (and the Lipschitz property of $\sigma$) to bound:

${\mathbf{E}}\left[\sup_{t\in[0,T]}\left|\int_{0}^{T}\left(\sigma(\theta^{x}_{s})-\sigma(\theta^{\bar{x}}_{s})\right)\text{\rm
d}B_{s}\right|^{2}\right]\;\leqslant\;CT^{2}\int_{0}^{T}{\mathbf{E}}\left[\sup_{u\in[0,s]}\left|\theta^{x}_{u}-\theta^{\bar{x}}_{u}\right|^{2}\right]\text{\rm
d}s,$

with $C$ a universal positive constant. The rest of the proof remains
unchanged. ∎

## Acknowledgments

The author is thankful to Gianmarco Bet, Giambattista Giacomin and Francesca
Nardi for discussions while writing this note.

## References

  * [1] E. Bayraktar, S. Chakraborty, and R. Wu.  Graphon mean field systems.  arXiv:2003.13180 [math], Mar. 2020.  arXiv: 2003.13180.
  * [2] E. Bayraktar and R. Wu.  Stationarity and uniform in time convergence for the graphon particle system.  arXiv:2008.10173 [math], 2020.
  * [3] G. Bet, F. Coppini, and F. R. Nardi.  Weakly interacting oscillators on dense random graphs.  arXiv:2006.07670 [math], June 2020.  arXiv: 2006.07670.
  * [4] S. Bhamidi, A. Budhiraja, and R. Wu.  Weakly interacting particle systems on inhomogeneous random graphs.  Stochastic Processes and their Applications, 129(6):2174–2206, June 2019.
  * [5] P. E. Caines and M. Huang.  Graphon Mean Field Games and the GMFG Equations.  2018 IEEE Conference on Decision and Control (CDC), 2018.
  * [6] R. Carmona, D. Cooney, C. Graves, and M. Lauriere.  Stochastic Graphon Games: I. The Static Case.  arXiv:1911.10664 [math], Nov. 2019.  arXiv: 1911.10664.
  * [7] H. Chiba and G. S. Medvedev.  The mean field analysis of the Kuramoto model on graphs i. The mean field equation and transition point formulas.  Discrete & Continuous Dynamical Systems - A, 39(1):131, 2019.
  * [8] F. Coppini.  Long time dynamics for interacting oscillators on graphs.  arXiv:1908.01520 [math], Dec. 2019.  arXiv: 1908.01520.
  * [9] F. Coppini, H. Dietert, and G. Giacomin.  A law of large numbers and large deviations for interacting diffusions on Erdős–Rényi graphs.  Stochastics and Dynamics, 20(02):2050010, Apr. 2020.
  * [10] S. Delattre, G. Giacomin, and E. Luçon.  A Note on Dynamical Models on Random Graphs and Fokker-Planck Equations.  Journal of Statistical Physics, 165(4):785–798, Nov. 2016.
  * [11] J.-F. Delmas, D. Dronnier, and P.-A. Zitt.  An Infinite-Dimensional SIS Model.  arXiv:2006.08241 [math], June 2020.  arXiv: 2006.08241.
  * [12] P. Diaconis and S. Janson.  Graph limits and exchangeable random graphs.  Rendiconti di Matematica, 28:33–61, 2008.
  * [13] R. Dudley.  Real Analysis and Probability.  Cambridge University Press, 2 edition, 2002.
  * [14] D. S. Kaliuzhnyi-Verbovetskyi and G. S. Medvedev.  The Mean Field Equation for the Kuramoto Model on Graph Sequences with Non-Lipschitz Limit.  SIAM J. Math. Analysis, 2017.
  * [15] D. Keliger, I. Horvath, and B. Takacs.  Local-density dependent Markov processes on graphons with epidemiological applications.  arXiv:2008.08109 [math], 2020.
  * [16] L. Lovász.  Large Networks and Graph Limits, volume 60 of Colloquium Publications.  American mathematical society edition, 2012.
  * [17] L. Lovász and B. Szegedy.  Limits of dense graph sequences.  Journal of Combinatorial Theory, Series B, 96(6):933–957, Nov. 2006\.
  * [18] E. Luçon.  Quenched asymptotics for interacting diffusions on inhomogeneous random graphs.  Stochastic Processes and their Applications, 2020.
  * [19] E. Luçon and W. Stannat.  Mean field limit for disordered diffusions with singular interactions.  The Annals of Applied Probability, 24(5):1946–1993, Oct. 2014.
  * [20] G. S. Medvedev.  The continuum limit of the Kuramoto model on sparse random graphs.  arXiv:1802.03787 [nlin], Feb. 2018.  arXiv: 1802.03787.
  * [21] F. Parise and A. Ozdaglar.  Graphon games: A statistical framework for network games and interventions.  arXiv:1802.00080 [cs], June 2020.  arXiv: 1802.00080.
  * [22] G. H. Reis and R. I. Oliveira.  Interacting diffusions on random graphs with diverging degrees: hydrodynamics and large deviations.  Journal of Statistical Physics, 2019.  arXiv: 1807.06898.
  * [23] A.-S. Sznitman.  Topics in propagation of chaos.  In P.-L. Hennequin, editor, Ecole d’Eté de Probabilités de Saint-Flour XIX - 1989, volume 1464, pages 165–251. Springer Berlin Heidelberg, 1991.


###### Abstract

We present new multi-epoch Very Long Baseline Array (VLBA) observations of a
set of TeV blazars drawn from our VLBA program to monitor all TeV-detected
high-frequency peaked BL Lac objects (HBLs) at parsec scales. Most of these
sources are faint in the radio, so they have not been well observed with VLBI
by other surveys. Our previous measurements of apparent jet speeds in TeV HBLs
showed apparent jet speeds that were subluminal or barely superluminal,
suggesting jets with velocity structures at the parsec-scale. Here we present
apparent jet speed measurements for eight new TeV HBLs, which for the first
time show a superluminal tail to the apparent speed distribution for the TeV
HBLs.

###### keywords:

TeV blazars; VLBI; active galactic nuclei

x 10.3390/—— xx Academic Editors: Jose L. Gómez, Alan P. Marscher and Svetlana
G. Jorstad Received: date; Accepted: date; Published: date PARSEC-SCALE
STRUCTURE AND KINEMATICS OF FAINT TEV HBLS B. Glenn Piner 1 and Philip G.
Edwards 2,∗ Glenn Piner and Philip Edwards Correspondence:
<EMAIL_ADDRESS>Tel.: +61-2-9372-4717

##  1 Introduction

At TeV energies (1012eV), three orders of magnitude higher than those studied
by satellite-based detectors, gamma-ray astronomy is conducted with ground
based telescopes such as H.E.S.S., VERITAS and MAGIC. Over 175 TeV gamma-ray
sources have now been catalogued (http://tevcat.uchicago.edu/), with over one
third of these being extragalactic objects. The majority of these (46 of 69)
are classified as HBL (High-frequency–peaked BL Lac) objects, for which the
synchrotron peak of the Spectral Energy Distribution (SED) lies at frequencies
above 1016.5 Hz.

Many of the well-studied HBL TeV sources have shown dramatic variability in
their gamma-ray emission Krawczynski et al. (2004); Aharonian et al. (2007).
The most rapid variations suggest extremely small emitting volumes and/or time
compression by large relativistic Doppler factors of up to $\sim$100 and
challenge our understanding of relativistic jets Aharonian et al. (2007);
Begelman et al. (2008); Ghisellini and Tavecchio (2008).

The only way to directly obtain information on the parsec-scale structure of
these blazar jets is by imaging the radio sources using the technique of Very
Long Baseline Interferometry (VLBI). However, many HBLs are fainter at radio
wavelengths (typically tens of milli-janskys — see Table 1) than the more
powerful quasars and BL Lac objects, and are not included in VLBI monitoring
programs such as MOJAVE Lister et al. (2016) and TANAMI Ojha et al. (2010).
Properties that can be measured from VLBI images — the apparent jet speed,
radio core brightness temperature, core dominance, and jet–to–counter-jet
brightness ratio — provide information on fundamental properties of the jet,
such as the bulk Lorentz factor and viewing angle.

Yet despite the high Doppler factors inferred from TeV observations, our
previous VLBA observations have established that TeV sources have only modest
brightness temperatures, and jet component motions that are sub-luminal or
only slightly superluminal Piner et al. (1999); Edwards and Piner (2002);
Giroletti et al. (2004); Piner and Edwards (2004, 2005); Piner et al. (2008,
2009); Tiet et al. (2012); Piner and Edwards (2014). This has been called
the“doppler crisis” Tavecchio (2006) or “bulk Lorentz factor crisis” Henri and
Saugé (2006).

Table 1: Current status of our VLBA monitoring program. TeV source names are those used by TeVCat, redshifts are those given in Piner and Edwards (2014), where an asterisk denotes a tentative value or limit. The NVSS flux density is measured at 1.4 GHz. The number of epochs refers to the number of VLBA images made in our monitoring program. References are to papers presenting these image. TeV source name | Redshift | NVSS counterpart | NVSS flux | Number | Refs
---|---|---|---|---|---
|  |  | density (mJy) | of epochs |
SHBL J001355.9$-$185406 | 0.094 | NVSS J001356$-$185406 | 29.2 | 4 | Piner and Edwards (2014)
KUV 00311$-$1938 | 0.506 * | NVSS J003334$-$192133 | 18.5 | 4 | Piner and Edwards (2014)
1ES 0033+595 | 0.240 * | NVSS J003552+595005 | 147.3 | 5 | Piner and Edwards (2014)
RGB J0136+391 | 0.400 * | NVSS J013632+390559 | 60.0 | 4 | Piner and Edwards (2014)
RGB J0152+017 | 0.080 | NVSS J015239+014717 | 61.4 | 5 | Piner and Edwards (2014)
1ES 0229+200 | 0.140 | NVSS J023248+201716 | 82.4 | 5 | Piner and Edwards (2013, 2014)
PKS 0301$-$243 | 0.266 | NVSS J030326$-$240710 | 700.2 | 2 | —
IC 310 | 0.019 | NVSS J031642+411928 | 168.1 | 2 | —
RBS 0413 | 0.190 | NVSS J031951+184536 | 20.9 | 5 | Piner and Edwards (2013, 2014)
1ES 0347$-$121 | 0.188 | NVSS J034922$-$115914 | 23.9 | 5 | Piner and Edwards (2013, 2014)
1ES 0414+009 | 0.287 | NVSS J041652+010526 | 119.6 | 5 | Piner and Edwards (2013, 2014)
1ES 0502+675 | 0.314 | NVSS J050755+673724 | 25.4 | 5 | Piner and Edwards (2013, 2014)
PKS 0548$-$322 | 0.069 | NVSS J055040$-$321620 | 343.7 | 5 | Piner and Edwards (2013, 2014)
RX J0648.7+1516 | 0.179 | NVSS J064847+151625 | 64.2 | 4 | Piner and Edwards (2014)
1ES 0647+250 | 0.450 | NVSS J065046+250259 | 96.2 | 5 | Piner and Edwards (2014)
RGB J0710+591 | 0.125 | NVSS J071030+590817 | 158.4 | 5 | Piner and Edwards (2013, 2014)
1ES 0806+524 | 0.138 | NVSS J080949+521858 | 182.4 | 3 | Piner and Edwards (2013)
RBS 0723 | 0.198 | NVSS J084712+113350 | 32.8 | 1 | —
1RXS J101015.9$-$311909 | 0.143 | NVSS J101015$-$311906 | 73.5 | 4 | Piner and Edwards (2014)
1ES 1011+496 | 0.212 | NVSS J101504+492601 | 377.7 | 3 | Piner and Edwards (2013)
1ES 1101$-$232 | 0.186 | NVSS J110337$-$232924 | 120.3 | 5 | Tiet et al. (2012)
Markarian 421 | 0.031 | NVSS J110427+381232 | 767.4 | 17 | Piner et al. (1999); Piner and Edwards (2005); Piner et al. (2010)
Markarian 180 | 0.045 | NVSS J113626+700925 | 327.1 | 7 | Tiet et al. (2012)
RX J1136.5+6737 | 0.134 | NVSS J113629+673706 | 45.3 | 1 | —
1ES 1215+303 | 0.130 | NVSS J121752+300700 | 571.6 | 2 | —
1ES 1218+304 | 0.184 | NVSS J122121+301036 | 71.0 | 5 | Tiet et al. (2012)
MS 1221.8+2452 | 0.218 | NVSS J122424+243623 | 24.5 | 4 | Piner and Edwards (2014)
H 1426+428 | 0.129 | NVSS J142832+424022 | 58.0 | 5 | Piner et al. (2008, 2010)
1ES 1440+122 | 0.163 | NVSS J144248+120040 | 68.8 | 4 | Piner and Edwards (2014)
PG 1553+113 | 0.500 * | NVSS J155543+111124 | 312.0 | 7 | Tiet et al. (2012)
Markarian 501 | 0.034 | NVSS J165352+394536 | 1558.0 | 9 | Edwards and Piner (2002); Piner et al. (2009, 2010)
H 1722+119 | 0.170 * | NVSS J172504+115215 | 120.4 | 4 | Piner and Edwards (2014)
1ES 1727+502 | 0.055 | NVSS J172818+501311 | 200.7 | 2 | —
1ES 1741+196 | 0.083 | NVSS J174357+193508 | 301.2 | 4 | Piner and Edwards (2014)
1ES 1959+650 | 0.047 | NVSS J195959+650854 | 249.6 | 13 | Piner and Edwards (2004, 2005); Piner et al. (2010)
PKS 2155$-$304 | 0.116 | NVSS J215852$-$301330 | 489.3 | 7 | Piner and Edwards (2004, 2005); Piner et al. (2010)
B3 2247+381 | 0.119 | NVSS J225005+382437 | 103.4 | 4 | Piner and Edwards (2014)
1ES 2344+514 | 0.044 | NVSS J234705+514217 | 250.4 | 8 | Piner and Edwards (2004); Piner et al. (2010)
H 2356$-$309 | 0.165 | NVSS J235907$-$303740 | 62.1 | 5 | Tiet et al. (2012)

Of the 46 HBLs detected at TeV energies:

  * •

11 have jet kinematics published previously by us. (Some of these are also in
MOJAVE.)

  * •

7 have, or soon will have, speeds determined by the MOJAVE program Lister et
al. (2016).

  * •

20 are included in the current phase of our program, with first epoch VLBA
results for all 20 published Piner and Edwards (2014). We present here the
kinematic results for the first portion of these sources.

  * •

4 are too far south to be studied with the VLBA, with several of these are
part of the TANAMI monitoring program Ojha et al. (2010).

  * •

4 are recent detections which are yet to be monitored with VLBI.

The 39 TeV sources for which we have at least one VLBA image are listed in
Table 1, together with their redshift, association in the NRAO VLA Sky Survey
(NVSS) catalog Condon et al. (1998), 1.4 GHz flux density, and details of our
VLBA observations. Images and data available at the project website:
www2.whittier.edu/facultypages/gpiner/research/archive/archive.html .

Figure 1: VLBA images at 8 GHz of four Tev blazars in out VLBA monitoring
program: (a) RGB J0152+017, (b) 1ES 0229+200, (c) RBS 0413 (0317+185), (d) 1ES
0347$-$121\.

TeV photons are attenuated by the infrared background Aharonian et al. (2006)
and, as is apparent in Table 1, the majority of extragalactic TeV sources are
at relatively low redshift ($z<$ 0.2). Conversely, studies of TeV gamma-ray
spectra have offered a means of constraining the infrared background Aharonian
et al. (2006); Stecker et al. (2016).

##  2 Results

Images of four of the sources currently being monitored are shown in Figure 1.
All show parsec-scale morphologies typical of this class: a compact core
(which hosts the supermassive black hole powering the source), and a weaker,
one-sided jet that transitions to a decollimated structure with larger opening
angle at a few tens of milli-arcseconds from the core (see image of 0229+200).
Multi-epoch studies of these jets over the course of several years allow the
apparent speeds of the jet components to be determined.

Our previous VLBA studies indicated the absence of rapidly moving features in
the jets of TeV HBLs; jet components were either nearly stationary or slowly
moving ($<\sim 1c$) Tiet et al. (2012). With the addition of multi-epoch data
from eight previously unpublished sources, including the four sources in
Figure 1, the revised distribution of apparent jet speeds is shown in Figure
2. This Figure incorporates the results of our previously published data and
jet speeds for four sources that have been monitored as part of the MOJAVE
project Lister et al. (2016).

With the addition of new data, the tail of the distribution now extends to
mildy superluminal apparent speeds for the first time; however, the majority
of the TeV HBLs have peak apparent speeds of only about 1$c$. Combining these
slow apparent speeds with the high Doppler factors ($\delta$) implied by the
TeV data to solve for the Lorentz factor ($\Gamma$) and viewing angle
($\theta$) results in unreasonbly small viewing angles ($\theta<<$ 1∘). This
would imply tiny jet opening angles, enormous linear sizes, and huge numbers
of parent objects, and indicates that the combination of both high Doppler
factor and slow apparent speed in the same jet region is unphysical. If more
realistic viewing angles of a few degrees are assumed, then the observed
apparent speeds imply more modest Lorentz and Doppler factors for the radio
jet. The lack of detection of counterjets in the VLBI images for any TeV HBL
Giroletti et al. (2008); Piner and Edwards (2014) requires that the Doppler
factor cannot be arbitrarily low, and values of $\delta$ and $\Gamma$ of a few
are most consistent with the combined VLBI data.

##  3 Discussion

A variety of mechanisms have been proposed to try and reconcile the Doppler
crisis Georganopoulos and Kazanas (2003); Gopal-Krishna et al. (2004);
Ghisellini et al. (2005); Gopal-Krishna et al. (2006, 2007); Marscher (2014).
The most natural explanation is for a range of Doppler factors to coexist in
the same jet on parsec scales through jet stratification. One example is a jet
that decelerates along its length Georganopoulos and Kazanas (2003). In such a
jet, the fast inner part sees blueshifted photons from the slower outer part,
reducing the high Lorentz factor required in the fast portion. This is a
general feature of models with velocity structures; radiative interaction
among the different regions allows the SED to be reproduced without the
extremely high Lorentz factors and Doppler factors characterizing single-zone
models. Another alternative is a transverse velocity structure with a fast
central spine and a slower outer sheath. Radiative interaction between the
spine and sheath naturally decelerates the spine, producing both radial and
transverse velocity structures in the same jet Ghisellini et al. (2005).

Figure 2: A histogram of the peak jet speeds in TeV blazars. Blue denotes
sources from our previously published data, red are from this work, and yellow
are from MOJAVE. The highest apparent speed are observed in RBS 0413,
(6.0$\pm$1.2)$c$, and RGB J0710+591, (5.8$\pm$1.5)$c$.

If such spine–sheath jets are present in TeV HBLs, then the outer layer is
expected to dominate the radio emission due to its SED shape, even with a
lower Doppler factor than the spine Ghisellini et al. (2005). An observational
signature of this would be a limb-brightened transverse profile for the jet in
VLBI images. There is evidence for limb-brightening close to the core in a
number of HBLs, e.g., Mkn 501 Piner and Edwards (2005), 1ES 0502+675 Piner and
Edwards (2014), and H 1722+119 Piner and Edwards (2014).

Other possible jet velocity structures have also been proposed, including
multiple blobs Tavecchio et al. (2011) fast moving “needles” within the main
jet Ghisellini and Tavecchio (2008), “minijets” powered by magnetic
reconnection events Giannios (2013), and turbulent subregions within the jet
Marscher (2014).

##  4 Conclusions

Our on-going VLBA monitoring of the growing number of TeV gamma-ray emitting
HBLs has revealed that the distribution of peak apparent jet speeds in these
sources extends to moderate superluminal speeds, $\sim$6$c$, but the majority
display subluminal speeds, in contrast with the distribution for other classes
of active galactic nuclei Lister et al. (2016). It has recently been proposed
that jet kinematics may offer a better classification for blazars than the SED
peak frequency Hervet et al. (2016), with HBLs tending to display quasi-
stationary knots arising from recollimation shocks.

A possible physical explanation for this is based on TeV blazars having
intrinsically weak jets that interact with the external medium forming a slow
surrounding layer. Radiative interaction between the spine and the sheath
decelerates the spine, and eventually disrupts the jet. Such jets are
prominent in TeV-selected samples because selection favors rare high-
synchrotron peak sources, which are drawn from the low end of the luminosity
function where the source density is largest Giommi et al. (2012).

###### Acknowledgements.

The National Radio Astronomy Observatory is a facility of the National Science
Foundation operated under cooperative agreement by Associated Universities,
Inc. This research has made use the TeVCat online source catalog
(http://tevcat.uchicago.edu). This research has made use of NASA’s
Astrophysics Data System. This research has made use of the NASA/IPAC
Extragalactic Database (NED) which is operated by the Jet Propulsion
Laboratory, California Institute of Technology, under contract with the
National Aeronautics and Space Administration.

## References

  * Krawczynski et al. (2004) Krawczynski, H.; Hughes, S.B.; Horan, D.; Aharonian, F.; Aller, M.F.; Aller, H.; Boltwood, P.; Buckley, J.; Coppi, P.; Fossati, G.; et al..  Multiwavelength Observations of Strong Flares from the TeV Blazar 1ES 1959+650.  Astrophys. J. 2004, 601, 151–164.
  * Aharonian et al. (2007) Aharonian, F.; Akhperjanian, A.G.; Bazer-Bachi, A.R.; Behera, B.; Beilicke, M.; Benbow, W.; Berge, D.; Bernlöhr, K.; Boisson, C.; Bolz, O.; et al..  An Exceptional Very High Energy Gamma-Ray Flare of PKS 2155-304.  Astrophys. J. Lett. 2007, 664, L71–L74.
  * Begelman et al. (2008) Begelman, M.C.; Fabian, A.C.; Rees, M.J.  Implications of very rapid TeV variability in blazars.  Mon. Not. R. Astron. Soc. 2008, 384, L19–L23.
  * Ghisellini and Tavecchio (2008) Ghisellini, G.; Tavecchio, F.  Rapid variability in TeV blazars: the case of PKS2155-304.  Mon. Not. R. Astron. Soc. 2008, 386, L28–L32.
  * Lister et al. (2016) Lister, M.L.; Aller, M.F.; Aller, H.D.; Homan, D.C.; Kellermann, K.I.; Kovalev, Y.Y.; Pushkarev, A.B.; Richards, J.L.; Ros, E.; Savolainen, T.  MOJAVE: XIII. Parsec-scale AGN Jet Kinematics Analysis Based on 19 years of VLBA Observations at 15 GHz.  Astron. J. 2016, 152, 12.
  * Ojha et al. (2010) Ojha, R.; Kadler, M.; Böck, M.; Booth, R.; Dutka, M.S.; Edwards, P.G.; Fey, A.L.; Fuhrmann, L.; Gaume, R.A.; Hase, H.; et al..  TANAMI: tracking active galactic nuclei with austral milliarcsecond interferometry. I. First-epoch 8.4 GHz images.  Astron. Astrophys. 2010, 519, A45.
  * Piner et al. (1999) Piner, B.G.; Unwin, S.C.; Wehrle, A.E.; Edwards, P.G.; Fey, A.L.; Kingham, K.A.  VSOP and Ground-based VLBI Imaging of the TeV Blazar Markarian 421 at Multiple Epochs.  Astrophys. J. 1999, 525, 176–190.
  * Edwards and Piner (2002) Edwards, P.G.; Piner, B.G.  The Subluminal Parsec-Scale Jet of Markarian 501.  Astrophys. J. Lett. 2002, 579, L67–L70.
  * Giroletti et al. (2004) Giroletti, M.; Giovannini, G.; Feretti, L.; Cotton, W.D.; Edwards, P.G.; Lara, L.; Marscher, A.P.; Mattox, J.R.; Piner, B.G.; Venturi, T.  Parsec-Scale Properties of Markarian 501.  Astrophys. J. 2004, 600, 127–140.
  * Piner and Edwards (2004) Piner, B.G.; Edwards, P.G.  The Parsec-Scale Structure and Jet Motions of the TeV Blazars 1ES 1959+650, PKS 2155-304, and 1ES 2344+514.  Astrophys. J. 2004, 600, 115–126.
  * Piner and Edwards (2005) Piner, B.G.; Edwards, P.G.  VLBA Polarization Observations of Markarian 421 after a Gamma-Ray High State.  Astrophys. J. 2005, 622, 168–177.
  * Piner et al. (2008) Piner, B.G.; Pant, N.; Edwards, P.G.  The Parsec-Scale Jets of the TeV Blazars H1426+428, 1ES 1959+650, and PKS 2155-304: 2001-2004.  Astrophys. J. 2008, 678, 64–77.
  * Piner et al. (2009) Piner, B.G.; Pant, N.; Edwards, P.G.; Wiik, K.  Significant Limb-Brightening in the Inner Parsec of Markarian 501.  Astrophys. J. Lett. 2009, 690, L31–L34.
  * Tiet et al. (2012) Tiet, V.C.; Piner, B.G.; Edwards, P.G.  Increasing the Number of TeV Blazars with Parsec-Scale Kinematics.  Fermi & Jansky Proceedings – eConf C1111101, 2012, [arXiv:astro-ph.HE/1205.2399].
  * Piner and Edwards (2014) Piner, B.G.; Edwards, P.G.  First-epoch VLBA Imaging of 20 New TeV Blazars.  Astrophys. J. 2014, 797, 25.
  * Tavecchio (2006) Tavecchio, F.  Extragalactic Jets: The High Energy View.  The Tenth Marcel Grossmann Meeting; Novello, M.; Perez Bergliaffa, S.; Ruffini, R., Eds., 2006, p. 512, [astro-ph/0401590].
  * Henri and Saugé (2006) Henri, G.; Saugé, L.  The Bulk Lorentz Factor Crisis of TeV Blazars: Evidence for an Inhomogeneous Pileup Energy Distribution?  Astrophys. J. 2006, 640, 185–195.
  * Piner and Edwards (2013) Piner, B.G.; Edwards, P.G.  The Parsec-Scale Structure of the Newer TeV Blazars.  European Physical Journal Web of Conferences, 2013, Vol. 61, European Physical Journal Web of Conferences, p. 04021.
  * Piner et al. (2010) Piner, B.G.; Pant, N.; Edwards, P.G.  The Jets of TeV Blazars at Higher Resolution: 43 GHz and Polarimetric VLBA Observations from 2005 to 2009.  Astrophys. J. 2010, 723, 1150–1167.
  * Condon et al. (1998) Condon, J.J.; Cotton, W.D.; Greisen, E.W.; Yin, Q.F.; Perley, R.A.; Taylor, G.B.; Broderick, J.J.  The NRAO VLA Sky Survey.  Astron. J. 1998, 115, 1693–1716.
  * Aharonian et al. (2006) Aharonian, F.; Akhperjanian, A.G.; Bazer-Bachi, A.R.; Beilicke, M.; Benbow, W.; Berge, D.; Bernlöhr, K.; Boisson, C.; Bolz, O.; Borrel, V.; et al..  A low level of extragalactic background light as revealed by $\gamma$-rays from blazars.  Nature 2006, 440, 1018–1021.
  * Stecker et al. (2016) Stecker, F.W.; Scully, S.T.; Malkan, M.A.  An Empirical Determination of the Intergalactic Background Light from UV to FIR Wavelengths Using FIR Deep Galaxy Surveys and the Gamma-Ray Opacity of the Universe.  Astrophys. J. 2016, 827, 6.
  * Giroletti et al. (2008) Giroletti, M.; Giovannini, G.; Cotton, W.D.; Taylor, G.B.; Pérez-Torres, M.A.; Chiaberge, M.; Edwards, P.G.  The jet of Markarian 501 from millions of Schwarzschild radii down to a few hundreds.  Astron. Astroph. 2008, 488, 905–914.
  * Georganopoulos and Kazanas (2003) Georganopoulos, M.; Kazanas, D.  Decelerating Flows in TeV Blazars: A Resolution to the BL Lacertae-FR I Unification Problem.  Astrophys. J. Lett. 2003, 594, L27–L30.
  * Gopal-Krishna et al. (2004) Gopal-Krishna.; Dhurde, S.; Wiita, P.J.  Do the Mildly Superluminal VLBI Knots Exclude Ultrarelativistic Blazar Jets?  Astrophys. J. Lett. 2004, 615, L81–L84.
  * Ghisellini et al. (2005) Ghisellini, G.; Tavecchio, F.; Chiaberge, M.  Structured jets in TeV BL Lac objects and radiogalaxies. Implications for the observed properties.  Astron. Astrophys. 2005, 432, 401–410.
  * Gopal-Krishna et al. (2006) Gopal-Krishna.; Wiita, P.J.; Dhurde, S.  Bulk motion of ultrarelativistic conical blazar jets.  Mon. Not. R. Astron. Soc. 2006, 369, 1287–1292.
  * Gopal-Krishna et al. (2007) Gopal-Krishna.; Dhurde, S.; Sircar, P.; Wiita, P.J.  Influence of the jet opening angle on the derived kinematical parameters of blazar jets having uniform and stratified bulk motion.  Mon. Not. R. Astron. Soc. 2007, 377, 446–452.
  * Marscher (2014) Marscher, A.P.  Turbulent, Extreme Multi-zone Model for Simulating Flux and Polarization Variability in Blazars.  Astrophys. J. 2014, 780, 87.
  * Tavecchio et al. (2011) Tavecchio, F.; Becerra-Gonzalez, J.; Ghisellini, G.; Stamerra, A.; Bonnoli, G.; Foschini, L.; Maraschi, L.  On the origin of the $\gamma$-ray emission from the flaring blazar PKS 1222+216.  Astron. Astrophys. 2011, 534, A86.
  * Giannios (2013) Giannios, D.  Reconnection-driven plasmoids in blazars: fast flares on a slow envelope.  Mon. Not. R. Astron. Soc. 2013, 431, 355–363.
  * Hervet et al. (2016) Hervet, O.; Boisson, C.; Sol, H.  An innovative blazar classification based on radio jet kinematics.  Astron. Astrophys. 2016, 592, A22.
  * Giommi et al. (2012) Giommi, P.; Padovani, P.; Polenta, G.; Turriziani, S.; D’Elia, V.; Piranomonte, S.  A simplified view of blazars: clearing the fog around long-standing selection effects.  Mon. Not. R. Astron. Soc. 2012, 420, 2899–2911.