Global stability for McKean–Vlasov equations on large networks

Abstract

We investigate the mean-field dynamics of stochastic McKean differential equations with heterogeneous particle interactions described by large network structures. To express a wide range of graphs, from dense to sparse structures, we incorporate the recently developed graph limit theory of graphops into the limiting McKean–Vlasov equations. Global stability of the splay steady state is proven via a generalised entropy method, leading to explicit graph structure-dependent decay rates. We highlight the robustness of the entropy approach by extending the results to the closely related Sakaguchi–Kuramoto model with intrinsic frequency distributions. We also present central examples of random graphs, such as power law graphs and the spherical graphop, and analyse the limitations of the applied methodology.

1. Introduction

Our analysis is motivated by stochastic McKean (or Kuramoto-type) differential equations [1, 2], which describe the behaviour of $N$ particles which diffuse and interact with each other. We are specifically interested in heterogeneous interaction patterns that go beyond the case that all particles interact with all others. General interactions are then described by a graph/network structure with the corresponding equation system given as:

(1) \begin{equation} \mathrm{d} X_t^i = -\frac{\kappa }{N r_N} \sum _{j\neq i}^N A^{ij} \nabla D(X_t^i - X_t^j)\, \mathrm{d} t + \sqrt{2}\, \mathrm{d} B_t^i, \quad i=1,\ldots, N, \end{equation}

where $X_t^i$ is the $i$ -th graph node with values on the $d$ -dimensional flat torus of length $L\gt 0$ , denoted as $U\,:\!=\,[-\frac{L}{2}, \frac{L}{2}]^d$ . The entries $A^{ij}\in \{0,1\}$ of the adjacency matrix $A^{(N)}\,:\!=\,(A^{ij})_{i,j=1,\ldots, N}$ represent the (undirected) graph edges and $1\geq r_N\gt 0$ is a rescaling factor relevant for sparse graphs. The nodes are coupled along edges according to a (periodic) interaction potential $D\,:\,U\to{\mathbb{R}}$ with relative coupling strength $\kappa \gt 0$ . The stochasticity in the system is modelled by independent Brownian motions $B_t^i$ on $U$ . Here, we describe the dynamics of (1) in the mean-field limit when the number of graph nodes tends towards infinity.

In the classical case of homogeneous interaction patterns, that is, $A^{ij}= 1$ for all $i,j=1,\ldots N$ in (1), the mean-field limit (see [3, 4]) is given as the (homogeneous) McKean–Vlasov equation:

(2) \begin{align} \begin{aligned} \partial _t \rho &= \kappa{\mathrm{div}}_x[\rho (\nabla _x D \star \rho )]+ \Delta _x \rho, \quad t\gt 0,\\ \rho (0) &= \rho _0, \end{aligned} \end{align}

where the solution $\rho (t,x)$ describes the average particle density at time $t\geq 0$ and position $x\in U$ . Such non-linear and non-local Fokker–Planck equations [5] have been analysed from a wide variety of perspectives. Most prominently, synchronisation phenomena for the Kuramoto model [6–8] and its (noisy) mean-field formulations [9–13] correspond to the potential $D(x) = -\cos \!(\frac{2\pi x}{L})$ in equation (2). Other important examples are the Hegselmann–Krause model [14, 15] in opinion dynamics and the Keller–Segel model for bacterial chemotaxis [16]. We refer to [14, 17] for a recent treatment of a large range of interaction potentials (for homogeneous interactions), proving global stability via entropy methods and providing a bifurcation analysis.

To incorporate heterogeneous interactions in the mean-field limit of (1), we require a theory of graph limits, which has made significant advances in recent years. For dense graphs, the infinite node limits can be expressed as graph functions, called graphons [18, 19]. A major shortcoming of this theory is that graphs with a subquadratic number of edge connections cannot be handled or the nodes are forced to have unbounded degree in the limit. Such graph cases require a refined treatment which has been developed for different subcases in a fragmented fashion. For sparse graphs of bounded degree, there is the local convergence theory of Benjamini and Schramm [20], as well as the stronger local–global convergence theory for measure-based graphings as limiting objects [21, 22]. Sparse graphs of unbounded degree, such as power law graphs, which are most challenging but crucial for applications, have been treated via rescaled graphon convergence theory for $L^p$ -graphons [23].

These different graph limit frameworks have slowly started to be incorporated into dynamical mean-field settings. For various limiting equations solution theory and finite-graph approximation has been rigorously introduced and first graph limit-dependent stability results have been shown [24–30]. For non-dense graphs, there exists an appropriate choice of the rescaling factor $r_N$ in (1) that results expressive limit. We refer to [23, 27, 31] for the detailed analysis.

In this work, we provide first entropy-based global stability results for solutions to the mean-field limit of (1). To include a wide range of interaction structures – from dense, to bounded and unbounded sparse graphs – we formulate the mean-field problem in the very general graph limit theory of graphops (or graph operators). Backhausz and Szegedy introduced this theory in [32], which unifies the above-mentioned approaches by lifting them to the more abstract level of probability space-based action convergence and operator theory. As we will see, this approach also has the advantage that it integrates rather naturally with the tools used in a PDE context. The mean-field formulation formally leads to the non-linear graphop McKean–Vlasov equation (see [33] for a discussion on the derivation):

(3) \begin{align} \begin{aligned} \partial _t \rho &= \kappa{\mathrm{div}}_x\Big (\rho V[A](\rho )\Big )+ \Delta _x \rho, \quad t\geq 0,\\ \rho (0) &= \rho _0, \end{aligned} \end{align}

where the solution $\rho (t,x,\xi )$ additionally depends on a graph limit variable $\xi \in \Omega$ (where $(\Omega, \mathcal{A},\mu )$ is a Borel probability space specified below). In this abstract formulation, the heterogeneous interactions are expressed by a bounded linear operator $A$ (with the precise properties stated in Section 3), resulting in

(4) \begin{equation} V[A](\rho )(t,x,\xi )\,:\!=\, \int _{U} \nabla _x D(x-\tilde{x})(A\rho )(t, \tilde{x},\xi ) \,\mathrm{d} \tilde{x}, \end{equation}

where we used the shorthand notation $(A\rho )(t, \tilde{x},\xi )\,:\!=\,[A\rho ](t, \tilde{x}, \cdot )(\xi )$ , that is, $A$ is an operator acting on the graph variable function space. Existence and uniqueness have been established in the recent years for prototypical Vlasov-type equations and their approximations of dynamics on finite networks [27, 31, 33–35]. Yet, as far as the authors are aware, entropy-based approaches have not yet been applied for any graphop PDE models. This work aims to provide a first application of the entropy method in this unifying graph operator setting. We emphasise the robustness of our method by also obtaining explicit global stability results for the closely related variant of (3) in 1D: the mean-field Sakaguchi–Kuramoto models with intrinsic frequency distributions combined with heterogeneous interactions. To our best knowledge, this is also the first entropy approach for mean-field Kuramoto models with heterogeneous interactions, including those involving graphons.

General assumption: Throughout the paper, we always assume the interaction potential satisfies $D\in \mathcal{W}^{2,\infty }(U)$ .

1.1 Main results of the paper

Our main result is formulated in Theorem15. There we assume reasonable regularity on a normalised initial datum and a square-integrable graphop $A$ . Then, it states that classical solutions of the graphop McKean–Vlasov equation (3) converge exponentially fast to the splay (or ‘incoherent’) steady state $\rho _\infty \,:\!=\, \frac{1}{L^d}$ in entropy, provided the interaction strength $\kappa$ is small enough. The bound on $\kappa$ as well as the exponential decay rate are explicit, and they depend on the operator norm of the graphop.

1.2 Structure of the paper

In Section 2, we review the case of global stability for homogeneous interactions and introduce the entropy method and necessary functional inequalities. Section 3 presents our main result of global stability for heterogeneous interaction patterns described by graphops. We further discuss the approaches limitations and the main theorem’s formulation for graphons. In Section 4, we apply our method to the Sakaguchi–Kuramoto model with added frequency distribution. Section 5 investigates explicit graphop examples such as spherical graphops and power law graphons.

2. Homogeneous interaction patterns

We first recall that for all-to-all coupled interactions – assuming $\kappa \gt 0$ small enough – the splay steady state $\rho _\infty \,:\!=\, \tfrac{1}{L^d}$ is a global steady state to which every initial configuration converges exponentially fast. This section we review the results developed in [14, 17] on which the following sections are built.

2.1 Global stability via entropy methods

We start by revisiting the homogeneous interaction case by assuming $A\rho = \rho$ in (3). This is the foundation for our strategy of proving decay estimates for more complex heterogeneous interactions. In the homogeneous case, the network variable $\xi$ can be omitted, $\rho (t,x,\xi )\equiv \rho (t,x)$ and $V[A](\rho ) = \nabla _x D \star \rho$ , that is, we recover equation (2). We will introduce the entropy method and related functional inequalities that we will later extend and use for the heterogeneous interaction case.

Definition 1. Let the space of absolutely continuous Borel probability measures on $U$ be denoted as $\mathcal{P}_{\text{ac}}(U)$ . Then, for each $\rho \in \mathcal{P}_{\text{ac}}(U)$ , we define the relative entropy functional:

(5) \begin{equation} H(\rho |\rho _\infty )\,:\!=\, \int _{U} \rho \log \left(\frac{\rho }{\rho _\infty }\right) \,\mathrm{d} x \end{equation}

with $\rho _\infty (x) \,:\!=\, \frac{1}{L^d}$ .

We use this functional to show exponential decay of solutions. The relative entropy is an adequate measure of “distance” between a function $\rho \in \mathcal{P}_{\text{ac}}(U)$ and the steady state $\rho _\infty$ . Indeed, $H(\rho |\rho _\infty )\geq 0$ for all $\rho \in \mathcal{P}_{\text{ac}}(U)$ , which can be seen by applying Jensen’s inequality with the convex function $x\log x$ . It further holds that $H(\rho _\infty |\rho _\infty ) = 0$ and the functional is linked to the $L^1(U)$ distance via the Csiszár–Kullback–Pinsker (CKP) inequality (see [36]):

(6) \begin{equation} \|\rho - \rho _\infty \|_{L^1(U)} \leq \sqrt{2 H(\rho |\rho _\infty )}. \end{equation}

For functions $\rho \in \mathcal{P}_{\text{ac}}(U)$ with $\sqrt{\rho } \in H^1(U)$ , it is bounded from above via the log-Sobolev inequality (see e.g. [37] for a direct proof):

(7) \begin{equation} H(\rho |\rho _\infty ) \leq \frac{L^2}{4\pi ^2} \int _{U} |\nabla _x \log (\rho )|^2 \rho \,\mathrm{d} x. \end{equation}

For the existence of classical solutions to the homogeneous equation, we refer to following result.

Theorem 2 ([17, Theorem 2.2] adapted from [14, Theorem 3.12]). For $\rho _0 \in H^{3 + d}(U) \cap \mathcal{P}_{\text{ac}}(U)$ , there exists a unique classical solution $\rho \in C^1(0,\infty ;\, C^2(U))$ of the homogeneous McKean–Vlasov equation (2) such that $\rho (t,\cdot ) \in \mathcal{P}_{\text{ac}}(U) \cap C^2(U)$ for all $t\gt 0$ . Additionally, it holds that $\rho (t,\cdot )\gt 0$ and $H(\rho (t)|\rho _\infty )\lt \infty$ for all $t\gt 0$ .

Proposition 3 ([17, Proposition 3.1]). Let the interaction potential of the homogeneous McKean–Vlasov equation (2) satisfy $D\in \mathcal{W}^{2,\infty }(U)$ . Let $\rho$ be the classical solution with initial data $\rho _0 \in H^{3 + d}(U) \cap \mathcal{P}_{\text{ac}}(U)$ and $H(\rho _0|\rho _\infty )\lt \infty$ . If further the coupling coefficient satisfies

(8) \begin{equation} \kappa \lt \frac{2 \pi ^2}{L^2 \| \Delta _x D \|_{L^\infty (U)} }, \end{equation}

then the solution $\rho$ is exponentially stable in entropy with the decay estimate

(9) \begin{equation} H(\rho (t)|\rho _\infty ) \leq \mathrm{e}^{-\alpha t} H(\rho _0|\rho _\infty ), \quad t\geq 0, \end{equation}

where

\begin{equation*} \alpha \,:\!=\, \frac {4 \pi ^2}{L^2} - 2\kappa \| \Delta _x D\|_{L^\infty (U)} \gt 0. \end{equation*}

The detailed proof can be found in [17, Proposition 3.1]. For later reference, we discuss the main steps in order to introduce the reader to the method and necessary inequalities that are generalised to the more involved heterogeneous interactions in Section 3.

Proof of Proposition 3. Our goal is to bound the time derivative of a solution in relative entropy $H(\rho (t)| \rho _\infty )$ by a negative multiple of the relative entropy itself:

(10) \begin{equation} \frac{\mathrm{d}}{\mathrm{d} t}H(\rho (t)|\rho _\infty ) \leq -\alpha H(\rho (t)|\rho _\infty ), \quad t\geq 0, \end{equation}

with some $\alpha \gt 0$ . Then Gronwall’s Lemma provides us with an exponential decay rate $\alpha$ of solutions in relative entropy.

To this end, we differentiate the relative entropy of solution $\rho$ with respect to time, which is possible as $\rho$ is a classical solution. Using equation (2), the general assumption $D\in \mathcal{W}^{2,\infty }(U)$ and integrating by parts leads to following two terms:

(11) \begin{align} \frac{\mathrm{d}}{\mathrm{d} t} H(\rho (t)|\rho _\infty ) &= \int _{U} \Big (\Delta _x \rho + \kappa{\mathrm{div}}_x (\rho \nabla _x D \star \rho ) )\Big ) \log \left(\frac{\rho }{\rho _\infty }\right)\,\mathrm{d} x\nonumber \\ & \qquad + \underbrace{\int _{U} \rho \frac{\rho _\infty }{\rho } \frac{1}{\rho _\infty } \Big (\Delta _x \rho + \kappa{\mathrm{div}}_x (\rho \nabla _x D \star \rho )\Big )\,\mathrm{d} x}_{=0} \nonumber \\ &= - \int _{U} |\nabla _x \rho |^2\frac{1}{\rho }\, \mathrm{d} x - \kappa \int _{U} \rho \nabla _x D \star \rho \left(\frac{\rho _\infty }{\rho }\right) \nabla _x \left(\frac{\rho }{\rho _\infty }\right)\,\mathrm{d} x\nonumber \\ &= - \int _{U} |\nabla _x\log (\rho )|^2 \rho dx + \kappa \int _{U} \rho \left(\Delta _x D \star \rho \right)\,\mathrm{d} x. \end{align}

The first term in (11) can be estimated via the log-Sobolev inequality (7). Note that if the second term vanishes, this estimate would provide us directly with an exponential decay result via Gronwall’s lemma.

To estimate the second term of interactions in (11), we replace $\rho$ by $\rho -\rho _\infty$ (possible due to $\int _{U} \Delta _x D(x) dx = 0$ following from the periodicity of $D$ ) and apply the CKP inequality (6) to estimate:

(12) \begin{align} \kappa \int _{U} \rho (\Delta _x D\star \rho ) dx &\leq \kappa \|\Delta _x D \|_{L^\infty (U)} \|\rho - \rho _\infty \|^2_{L^1(U)} \nonumber \\ &\leq 2\kappa \|\Delta _x D \|_{L^\infty (U)} H(\rho |\rho _\infty ). \end{align}

In total, identity (11) in combination with (12) yields

(13) \begin{equation} \frac{d}{dt} H(\rho |\rho _\infty ) \leq (-\frac{4 \pi ^2}{L^2} + 2\kappa \|\Delta _x D\|_\infty ) H(\rho |\rho _\infty ), \end{equation}

which proves the claimed result.

Remark 4. The condition on the coupling coefficient (8) can be further relaxed for coordinate-wise even interaction potentials, that is, $D(\ldots, x_i,\ldots ) = D(\ldots, -x_i, \ldots )$ for $i=1,\ldots, d$ . Then the term $\|\Delta _x D\|_{L^\infty (U)}$ can be replaced with $\|\Delta _x D_u\|_{L^\infty (U)}$ , where $D_u$ is the interaction contribution of the Fourier coefficients with negative sign (see details in [17]).

Remark 5. For $d=1$ , that is $U = [-\frac L2, \frac L2]$ , the noisy Kuramoto model ¹ assumes an interaction of form $D(x) = -\cos \!((\frac{2\pi }{L}) x)$ . For this model, the above condition (8) simplifies to $\kappa \lt \frac 12$ .

3. Heterogeneous interaction patterns

Before proving global stability for heterogeneous couplings, we first gather the necessary precise definitions for graphop McKean–Vlasov equations (3) and derived facts on graphops that are relevant to our setting. For an extensive introduction to graphops, we refer to [32].

Definition 6.

• • Let $(\Omega, \mathcal{A},\mu )$ be a Borel probability space with ² $\mathrm{supp}\,{\mu } = \Omega$ . Let further $L^p(\Omega )\,:\!=\,L^p(\Omega, \mu )$ , $p\in [1,\infty ]$ be the corresponding real-valued Lebesgue spaces.

• • We call $P$ -operators linear operators $A\,:\,L^\infty (\Omega ) \to L^1(\Omega )$ with bounded norm:

  \begin{equation*} \|A\|_{\infty \to 1}\,:\!=\, \sup _{v\in L^\infty (\Omega )}\frac {\|A v\|_{L^1(\Omega )}}{\|v\|_{L^\infty (\Omega )}}. \end{equation*}

• • We call $P$ -operators $A$ self-adjoint if for all $f,g\in L^\infty (\Omega )$ it follows that

  \begin{equation*} (Af,g)_{L^2(\Omega )}\,:\!=\, \int _\Omega gAf\, \mathrm {d}\mu (\xi ) = \int _\Omega f Ag\,\mathrm {d}\mu (\xi )= (Ag,f)_{L^2(\Omega )}. \end{equation*}
  Note this definition of self-adjointness includes operators not necessarily acting on Hilbert spaces.

• • A $P$ -operator $A$ is $c$ -regular for some $c\in{\mathbb{R}}$ if $A1_\Omega = c 1_\Omega$ .

• • An operator $A$ is said to be positivity-preserving if for any function $f\geq 0$ it follows that also $Af\geq 0$ .

• • Graphops are positivity-preserving and self-adjoint $P$ -operators.

• • The graphop $A$ is associated with a graphon kernel $W\,:\, \Omega \times \Omega \mapsto{\mathbb{R}}$ , if

  (14) \begin{equation} (A\rho )(\xi )\,:\!=\, \int _{\Omega } W(\xi, \tilde{\xi }) \rho (\tilde{\xi })\,\mathrm{d}\mu (\tilde{\xi }). \end{equation}
  As $A$ is self-adjoint, it follows that the graphon is symmetric, that is, $W(\xi, \tilde{\xi }) = W(\tilde{\xi },\xi )$ for all $\xi, \tilde{\xi } \in \Omega$ .

• • The $P$ -operator $A$ is of finite $(p,q)$ -norm, where $1\leq p,q\leq \infty$ , if the expression

  (15) \begin{equation} \|A\|_{p \to q}\,:\!=\, \sup _{v\in L^\infty (\Omega )}\frac{\|A v\|_{L^{q}(\Omega )}}{\|v\|_{L^p (\Omega )}} \end{equation}
  is finite. By standard arguments, such an operator can be uniquely extended to a bounded operator $A\,:\,L^p(\Omega ) \to L^q(\Omega )$ , which again is a $P$ -operator (or graphop if the original $A$ was). The set of such $P$ -operators is denoted as $B_{p,q}(\Omega )$ .

Remark 7. As $(\Omega, \mathcal{A}, \mu )$ is a probability space, we have $L^p(\Omega ) \subseteq L^{p'}(\Omega )$ for $ 1\leq p' \leq p\leq \infty$ . As a result, it holds true that $B_{p,q}(\Omega ) \subseteq B_{p',q'}(\Omega ) \subseteq B_{\infty, 1}(\Omega )$ for $p'\geq p$ , $q'\leq q$ .

3.1 Existence of classical solutions

The splay steady state of (3) is given by $\rho _\infty (x,\xi )\,:\!=\, \frac{1}{L^d}$ . This means $\rho _\infty$ satisfies the stationary equation:

\begin{equation*}\kappa {\mathrm {div}}_x\Big (\rho V[A](\rho )\Big )+ \Delta _x \rho = 0.\end{equation*}

This can be seen as $\rho _\infty$ is constant in both $x$ and $\xi$ and thus

(16) \begin{equation} \Delta _x D \star (A\rho _\infty )= A \rho _\infty \int _{U} \Delta _x D\,\mathrm{d} x = 0 \end{equation}

where we again used the periodicity of the potential $D$ .

Definition 8. For the probability space $(\Omega, \mathcal{A}, \mu )$ , the relative entropy (for heterogeneous coupling) is chosen as:

(17) \begin{equation} \hat{H}(\rho |\rho _\infty )\,:\!=\, \int _{\Omega }\int _{U} \rho \log \left(\frac{\rho }{\rho _\infty }\right) \, \mathrm{d} x \,\mathrm{d}\mu (\xi ). \end{equation}

Observe that (17) can be viewed as the average entropy over the heterogeneous node space. As for the homogeneous entropy $H(\rho |\rho _\infty )$ considered in (5), the CKP inequality (12) (now applied to the measure $\mathrm{d} x\times \mathrm{d} \mu$ on $U \times \Omega$ ) provides the lower bound:

\begin{equation*} \|\rho - \rho _\infty \|_{L^1(U\times \Omega )} \leq \sqrt {2\hat {H}(\rho |\rho _\infty )}, \quad \rho \in \mathcal {P}_{\text {ac}}(U\times \Omega ). \end{equation*}

We assumed the measure of $L^1(\Omega, \mu )$ satisfies $\mathrm{supp}\, \mu = \Omega$ , thus $\hat{H}(\rho |\rho _\infty ) = 0$ precisely when $\rho = \rho _\infty$ .

Definition 9. We say an initial datum $\rho _0$ of (3) is admissible if, for almost every $\xi \in \Omega$ , it holds that $\rho _0(\cdot, \xi )\in H^{3+d}(U)\cap \mathcal{P}_{\text{ac}}(U)$ . Additionally, for each $x\in U$ , we have $\rho _0(x,\cdot ) \in L^\infty (\Omega, \mu )$ .

Proposition 10. Let $\rho _0$ be an admissible initial datum for the graphop McKean–Vlasov equation (3) with arbitrary graphop $A$ . Then it follows that $\rho (\cdot, \cdot, \xi )$ is a unique classical solution (of equation (3) with $\xi$ -fixed interaction term) for a.e. $\xi \in \Omega$ . Furthermore, $\rho (t,\cdot, \xi ) \in H^{3 + d}(U) \cap \mathcal{P}_{\text{ac}}(U)$ , $\rho (t,\cdot, \cdot ) \in \mathcal{P}_{\text{ac}}(U\times \Omega )$ and $\rho (t,\cdot, \cdot )\gt 0$ for a.e. $\xi \in \Omega$ and all $t\gt 0$ .

Proof. We show that for each fixed $\xi \in \Omega$ and time $T\gt 0$ , the existence proof for a classical non-negative solution of [14] can be applied to our case. To see this, we show that the graphop $A$ still allows the necessary bounds for compactness arguments in norms w.r.t. $t$ and $x$ . Consider the solution $\rho _n(t,x,\xi )$ for $n\in \mathbb{N}$ to the frozen linear equation with smooth initial data $\rho _0(x,\xi ) \in \mathcal{P}_{\text{ac}}(U)\cap C^\infty (\overline{U})$ :

(18) \begin{equation} \begin{cases} \partial _t\rho _n(t,x,\xi ) = \Delta _x \rho _n(t,x,\xi ) \\ \qquad + \kappa\,{\mathrm{div}}_x[\rho _{n}(t,x,\xi ) \nabla _x D \star (A \rho _{n-1})(t,x,\xi ) ], & t\in [0,T], x\in U,\\ \rho _n(t,x,\xi ) = \rho _n(t,x+Le_i,\xi ), \quad t\in [0,T],& x\in \partial U,\\ \rho _n(0,x,\xi ) = \rho _0(x,\xi ), &x\in U. \end{cases} \end{equation}

For a.e. fixed $\xi \in \Omega$ , the frozen equation can be solved with classical results of linear parabolic equations with bounded coefficients. This leads to a sequence of unique solutions $(\rho _n(\xi ))_{n\in \mathbb{N}} \in \mathcal{P}_{\text{ac}}(U)\cap C^\infty (\overline{U} \times [0,T])$ . We show that the sequence $(\rho _n)_{n\in \mathbb{N}}$ is bounded in the desired norms, by multiplying the linearised equation (18) with $\rho _n(t,\xi )$ and integrating the $x$ -variable. This leads to

(19) \begin{align} \frac{\mathrm{d}}{\mathrm{d} t} \|\rho _n(t,\xi ) \|^2_{L^2(U)} &+ 2\|\nabla _x \rho _n(t,\xi ) \|^2_{L^2(U)} \nonumber \\ &\leq 2\kappa \int _U | \rho _n(t,\xi ) \nabla _x D \star A\rho _{n-1}(\xi ) \nabla _x \rho _n(t,\xi ) | dx\nonumber \\ &\leq \frac{\kappa \epsilon ^2}{2} \|\nabla _x \rho _n(t,\xi ) \|^2_{L^2(U)} \nonumber \\ &\qquad + \frac{2\kappa }{\epsilon ^2} \|\rho _n(t,\xi ) \|^2_{L^2(U)} \|\nabla _x D \star A \rho _{n-1}(t,\xi ) \|^2_{L^\infty (U)}. \end{align}

To estimate the last term above, we use the fact that $A$ , as a linear bounded operator acting solely on the network variable, commutes with the $x$ -integral and that the frozen equation (18) is mass-preserving, that is, $\|\rho _n(t,\xi )\|_{L^1(U)} = 1$ for $a.e.\,\xi \in \Omega$ and all $t\gt 0$ :

(20) \begin{align} \|\nabla _x D \star A \rho _{n-1}(t,\xi ) \|^2_{L^\infty (U)} &\leq \|\nabla _x D \|^2_{L^\infty (U)} \| A \rho _{n-1}(t,\xi ) \|^2_{L^1(U)} \nonumber \\ &= \|\nabla _x D\|^2_{L^\infty (U)} (A \| \rho _{n-1}(t,\xi ) \|_{L^1(U)})^2 \nonumber \\ &= \|\nabla _x D\|^2_{L^\infty (U)} (A 1_\Omega (\xi ))^2. \end{align}

Inserting this estimate in (19) and choosing $\epsilon = (2\kappa )^{-\frac 12}$ , we obtain

(21) \begin{align} \begin{aligned} \frac{d}{dt} \|\rho _n(t,\xi ) \|^2_{L^2(U)} +& \|\nabla _x \rho _n(t,\xi ) \|^2_{L^2(U)} \\ &\leq \kappa ^2 \|\nabla _x D\|^2_{L^\infty (U)} (A 1_\Omega (\xi ))^2 \|\rho _n(t,\xi ) \|^2_{L^2(U)}. \end{aligned} \end{align}

Gronwall’s Lemma yields the upper bound

(22) \begin{equation} \|\rho _n(t,\xi ) \|^2_{L^2(U)} \leq C(\xi, T) \|\rho _0(\xi ) \|^2_{L^2(U)},\quad n\in \mathbb{N}, \end{equation}

with

\begin{equation*} C(\xi, T)\,:\!=\, \exp \left (\kappa ^2 \|\nabla _x D \|_{L^\infty (U)}^2 (A1_\Omega (\xi ))^2 T\right ), \end{equation*}

which is finite for $a.e.\ \xi \in \Omega$ . Integrating (21) w.r.t. the time variable and using the bound (21) yields a uniform-in- $n$ bound for $(\rho _n(\xi ))_{n\in \mathbb{N}}$ in $L^2(0,T,H^1(U))$ . The further bootstrapping steps for initial datum of general regularity, uniqueness, desired solution regularity, and positivity of solutions follow as described in [14, 17]. Let us comment on the specific requirement $\rho _0 \in H^{3+d}(U)$ . It is needed in order to show, using the equation’s specific structure, that $\partial _{t}\rho \in L^2(0,\infty ;\,H^{2+d}(U))$ . As the embedding $H^{2+d}(U) \hookrightarrow H^{1+d}(U)$ is compact, it follows by Aubin–Lions Lemma that $\partial _t \rho \in C(0,\infty ;\,H^{1+d}(U))$ . Using the compactness of the Sobolev embedding one more time yields the desired $\partial _t \rho \in C(0,\infty ;\, C(U))$ .

For each $t\gt 0$ and a.e. $\xi \in \Omega$ , the solution is mass-preserving; hence, $\rho (t,\cdot, \xi )\in \mathcal{P}_{\text{ac}}(U)$ . As

\begin{equation*} \int _\Omega \int _U \rho (t,x,\xi ) \,\mathrm {d} x \,\mathrm {d}\mu (\xi ) = \int _\Omega 1\,\mathrm {d}\mu (\xi ) = 1, \end{equation*}

Fubini–Tonelli’s theorem implies that $\rho (t,\cdot, \cdot ) \in \mathcal{P}_{\text{ac}}(U\times \Omega )$ .

If we add further assumptions on the graphop, we obtain some regularity in the $\xi$ -variable.

Corollary 11. Let $\rho _0$ be an admissible initial datum for the graphop McKean–Vlasov equation (3) with graphop $A$ . If $A\in B_{2,2}(\Omega )$ holds and $\hat{H}(\rho _0|\rho _\infty ) \lt \infty$ , then $\hat{H}(\rho (t)|\rho _\infty )\lt \infty$ for all $t\gt 0$ . If $A1_\Omega (\xi ) \leq c$ for almost every $\xi \in \Omega$ with some constant $c\geq 0$ , then for all $t\gt 0$ it holds that $\rho (t, \cdot, \xi ) \in H^{3+d}(U)\cap \mathcal{P}_{\text{ac}}(U)$ for a.e. $\xi \in \Omega$ and $\rho (t, x,\cdot )\in L^\infty (\Omega )$ for all $x\in U$ .

Proof. The proof of $\hat{H}(\rho (t)|\rho _\infty )\lt \infty$ for all $t\gt 0$ , assuming $A\in B_{2,2}(\Omega )$ , is deferred to the proof of Theorem15 as the steps are identical but finiteness works for arbitrary $\kappa \geq 0$ .

If $A1_\Omega (\xi )\leq c$ , then the estimate (20) can be refined further by:

\begin{align*} \|\nabla _x D \star A \rho _{n-1}(\xi ) \|^2_{L^\infty (U)} &\leq \|\nabla _x D\|^2_{L^\infty (U)} c^2. \end{align*}

Plugging this into (22) results in

\begin{equation*} \sup _{\xi \in \Omega } \|\rho _n(t,\xi ) \|^2_{L^2(U)} \leq \exp \left [c^2 \|\nabla _x D \|^2_{L^\infty (U)} t\right ] \sup _{\xi \in \Omega } \|\rho _0(\xi ) \|_{L^2(U)} \end{equation*}

and proves the claim.

Remark 12. The assumption $A1_\Omega (\xi ) \leq c$ of Corollary 11 clearly includes $c$ -regular graphops with any $c\geq 0$ .

Remark 13. The regularity results of Proposition 10 and Corollary 11 with regard to the network variable $\xi$ are expected to be improvable. One would hope that for the weaker assumption $\|A\|_{2\to 2}\lt \infty$ that the solutions remain $L^2(\Omega )$ -integrable for all times. However, the above estimates (21) do not provide such a bound. Part of the difficulty is rooted in the fact that the evolution equation (3) provides no explicit regularising term over time in the $\xi$ -variable and there is no mixing (and not even a geometrical link) of the $x$ and $\xi$ directions. This is reminiscent of a parabolic equation with degenerate diffusion (in $\xi$ ) for which no hope of a coupling mechanism is present [38].

Similar existence analysis for the non-linear heat equation on sparse graphs [26 , Section 3] and the Kuramoto model [31] show $L^2(\Omega )$ regularity of solutions. However, also these existence results are restricted to the cases of $c$ -regular graphops.

3.2 Global stability

Before establishing the global convergence result for solutions to (3) in entropy (17), we recall the definition of the numerical radius of an operator. This quantity allows the sharpest formulation of the convergence rate which is possible with our method.

Definition 14. The numerical radius of a graphop $A\in B_{2,2}(\Omega )$ is given as:

\begin{equation*} n(A)\,:\!=\,\sup \{(Af,f) \mid f\in L^2(\Omega ), \|f\|_{L^2(\Omega )} = 1 \}. \end{equation*}

For graphops in $B_{2,2}(\Omega )$ , the numerical radius is an equivalent norm to the operator norm ³ with bounds:

(23) \begin{equation} n(A) \leq \|A\|_{2\to 2} \leq 2 n(A). \end{equation}

Theorem 15. Consider the graphop McKean–Vlasov equation (3) with any graphop $A$ that satisfies $n(A) \lt \infty$ . Let $\rho _0$ be any admissible initial datum with $\hat{H}(\rho _0|\rho _\infty )\lt \infty$ . Let further the coupling coefficient satisfy

(24) \begin{equation} \kappa \lt \frac{2 \pi ^2}{L^2 \| \Delta _x D \|_{L^\infty (U)} n(A) }, \end{equation}

Then, the classical solution $\rho$ is exponentially stable with the decay estimate

(25) \begin{equation} \hat{H}(\rho (t)|\rho _\infty ) \leq \mathrm{e}^{-\hat{\alpha }(A) t} \hat{H}(\rho _0|\rho _\infty ), \quad t\geq 0, \end{equation}

where

\begin{equation*} \hat {\alpha }(A)\,:\!=\, \frac {4 \pi ^2}{L^2} - 2\kappa \| \Delta _x D\|_{L^\infty (U)}n(A) \gt 0. \end{equation*}

Proof. Our strategy is to generalise the proof of Proposition 3.

Let us first note that, as $\|\rho (t,\xi )\|_{L^1(U)} =1$ for a.e. $\xi \in \Omega$ and $t\gt 0$ , the mapping $\xi \mapsto \|\rho (t,\xi ) - \rho _\infty \|_{L^1(U)}$ has finite $L^p(\Omega )$ -norm, for any $p\in [1,\infty ]$ . Specifically, $\xi \mapsto \|\rho (t,\xi ) - \rho _\infty \|_{L^1(U)} \in L^2(\Omega )$ .

As $\rho (\xi )$ is a classical solution for a.e. $\xi$ , it follows that the mapping $(t\mapsto H(\rho (t)|\rho _\infty )) \in C^1(0,\infty )$ . Then, $(t\mapsto \hat{H}(\rho (t)|\rho _\infty ))\in C^1(0,\infty )$ follows from the calculation below together with the $L^\infty (\Omega )$ bound noted above. The same calculations as for (11), and again using the log-Sobolev inequality (7) (for the measure $\mathrm{d} x\times \mathrm{d}\mu$ ), yields for $t\geq 0$

(26) \begin{equation} \frac{d}{dt} \hat{H}(\rho |\rho _\infty ) \leq - \frac{4\pi ^2}{L^2}\hat{H}(\rho |\rho _\infty ) + \kappa \int _{\Omega }\int _{U} \rho [\Delta _x D \star (A\rho )] \,\mathrm{d} x \,\mathrm{d}\mu (\xi ). \end{equation}

To estimate the second term in (26), we use the fact that $A$ , as a linear bounded operator acting solely on the network variable, commutes with the $x$ -integral:

(27) \begin{equation} \Delta _x D \star (A\rho )= A\,\left(\Delta _x D \star \rho \right). \end{equation}

For the special case $\rho \equiv \rho _\infty$ , we even have

(28) \begin{equation} \Delta _x D \star (A\rho _\infty )= A \rho _\infty \int _{U} \Delta _x D \,\mathrm{d} x = 0, \end{equation}

where the last equality follows again from the periodicity of $D$ . Hence, we can replace both occurrences of $\rho$ by $\rho -\rho _\infty$ and use Hölder’s inequality in $U$ with $p = 1$ and $p^* = \infty$ to estimate

(29) \begin{align} \kappa \int _{\Omega }&\int _{U} \rho (\Delta _x D \star (A\rho ) )\, \mathrm{d} x \,\mathrm{d}\mu (\xi ) \nonumber \\ &\leq \kappa \|\Delta _x D \|_{L^\infty (U)} \int _{\Omega } \|A[\rho - \rho _\infty ]\|_{L^1(U)} \|\rho - \rho _\infty \|_{L^1(U)}\,\mathrm{d}\mu (\xi ). \end{align}

Given that a graphop $A$ preserves positivity and denoting $\nu _\pm = \mp \max \{0,\pm \nu \}$ , it follows that

\begin{equation*}|A\nu | = |A\nu _+ - A\nu _-| \leq |A\nu _+| + |A\nu _-| = A\nu _+ + A\nu _- = A|\nu |.\end{equation*}

With this, we can estimate (29), use the definition of the numerical radius of $A$ , and then apply the CKP inequality (6) for a.e. $\xi \in \Omega$ . This leads to

\begin{align*} \kappa \int _{\Omega }&\int _{U} \rho (\Delta _x D \star (A\rho ) ) \,\mathrm{d} x\,\mathrm{d}\mu (\xi )\\ &\leq \kappa \|\Delta _x D \|_{L^\infty (U)} \int _{\Omega } (A\|\rho - \rho _\infty \|_{L^1(U)}) \|\rho - \rho _\infty \|_{L^1(U)} \,\mathrm{d}\mu (\xi )\\ &\leq \kappa \|\Delta _x D \|_{L^\infty (U)} n(A)\int _{\Omega } \|\rho - \rho _\infty \|_{L^1(U)}^2 \,\mathrm{d}\mu (\xi )\\ &\leq 2\kappa \|\Delta _x D \|_{L^\infty (U)} n(A)\hat{H}(\rho |\rho _\infty ). \end{align*}

Finally, applying Gronwall’s lemma to the estimation of (26) leads to the desired result.

Remark 16. Theorem 15 naturally includes graphops that satisfy the (stronger) condition $\|A\|_{p\to q}\lt \infty$ with $p\lt 2$ , $q\gt 2$ . In such cases, the conditions of Theorem 15 are still met since $n(A)\leq \|A\|_{2\to 2} \leq \|A\|_{p\to q}\lt \infty$ . But compared to the $(p,q)$ -norms for $p\leq 2, q\geq 2$ , the numerical radius always gives the sharpest estimates with regard to our line of estimations.

3.3 Limitations of the entropy method approach

Generalising the exponential stability of Theorem15 to graphops that do not have finite $\|A\|_{2\to 2}$ norm seems not feasible with our method of proof. This is due to the necessity of relating the time derivative of the defined entropy back to the entropy itself (10). Nonetheless, below we show that in these cases, solutions at least remain bounded for all times and are therefore not exponentially unstable.

Let us look at the considerations in detail:

1. (i) For exponential decay, the proof of Theorem15 requires the existence of a constant $C(A)\lt \infty$ such that

   (30) \begin{align} \begin{aligned} \int _{\Omega } ( A\|\rho - \rho _\infty \|_{L^1(U)})& \|\rho - \rho _\infty \|_{L^1(U)} \,\mathrm{d}\mu (\xi )\\ &\leq C(A) \int _{\Omega } \|\rho - \rho _\infty \|_{L^1(U)}^2 \,\mathrm{d}\mu (\xi ). \end{aligned} \end{align}
   This is necessary ⁴ to subsequently apply the CKP inequality (6) for a.e. $\xi \in \Omega$ and close the differential inequality. But if $A$ is unbounded in $L^2(\Omega )$ , no such constant can exist.

   Let us first assume $\|A\|_{2\to 2}\lt \infty$ . Then the numerical radius $n(A)$ is finite as well due to the norm equivalence (23). From the definition of the numerical radius follows that $C(A) = n(A)$ is optimal in (30).

   Now if $A$ is an unbounded operator in $L^2(\Omega )$ with $\|A\|_{2\to 2} = \infty$ , this implies that the numerical radius $n(A)$ is also unbounded and hence no constant $C(A)\lt \infty$ exists to bound (30). This follows again from (23), as we can approximate $A$ by a sequence of bounded operators $A_n$ with growing operator norm such that $\lim _{n\to \infty } \|A\rho -A_n\rho \|_{L^2(\Omega )} = 0$ for all $\rho \in \mathcal{D}(A)\subseteq L^2(\Omega )$ .

2. (ii) If $\|A\|_{2\to 2}$ is unbounded, one could still consider a Hölder inequality estimate. For the case that only the weaker condition $\|A\|_{p\to 2}\lt \infty$ for some $p\gt 2$ is satisfied, the Cauchy–Schwarz inequality and $\|A\rho \|_{L^2(\Omega )} \leq \|A\|_{p\to 2} \|\rho \|_{L^p(\Omega )}$ yields

   \begin{align*} \begin{aligned} &\int _{\Omega } ( A\|\rho - \rho _\infty \|_{L^1(U)}) \|\rho - \rho _\infty \|_{L^1(U)} \,\mathrm{d}\mu (\xi )\leq \\ &\quad \|A\|_{p\to 2} \left (\int _{\Omega } \|\rho - \rho _\infty \|_{L^1(U)}^{p} \,\mathrm{d}\xi \right )^{\frac 1{p}}\left (\int _{\Omega } \|\rho - \rho _\infty \|_{L^1(U)}^{2}\,\mathrm{d}\mu (\xi )\right )^{\frac 1{2}}. \end{aligned} \end{align*}
   However, as $p\gt 2$ the $L^p(\Omega )$ -term cannot be bounded by an $L^2(\Omega )$ -term and hence the differential inequality cannot be closed by the CKP inequality.

   Similarly, if only the condition $\|A\|_{2\to q}\lt \infty$ for some $q\lt 2$ holds true, then Hölder’s inequality with $\frac 1q + \frac 1{q^*} = 1$ leads to a $L^{q^*(\Omega )}$ -norm factor and as $q^* \gt 2$ again it cannot be bounded by $L^2(\Omega )$ .

3. (iii) Let us consider to modify the entropy functional (17). Instead of integrating $\xi$ with respect to the underlying measure of the probability space $(\Omega, \mathcal{A}, \mu )$ , we could consider an additional probability Borel measure $\tilde{\mu }$ on $\Omega$ :

   (31) \begin{equation} \hat{H}_{\tilde{\mu }}(\rho |\rho _\infty )\,:\!=\,\int _\Omega \int _U \rho \log \left(\frac{\rho }{\rho _\infty }\right)\,\mathrm{d} x\,\mathrm{d}\tilde{\mu }(\xi ). \end{equation}
   Then closing of the differential inequality for $H_{\tilde{\mu }}(\rho |\rho _\infty )$ along the proof of Theorem15 requires a constant $\tilde{C}(A)\lt \infty$ such that
   (32) \begin{align} \begin{aligned} \int _{\Omega } ( A\|\rho - \rho _\infty \|_{L^1(U)}) &\|\rho - \rho _\infty \|_{L^1(U)}\,\mathrm{d}\tilde{\mu }(\xi ) \\ &\leq \tilde{C}(A) \int _{\Omega } \|\rho - \rho _\infty \|_{L^1(U)}^2 \,\mathrm{d}\tilde{\mu }(\xi ). \end{aligned} \end{align}
   This is possible if and only if $A$ has a bounded numerical radius in the $L^2(\Omega, \tilde{\mu })$ sense. For the counter example of power law graphons without finite $L^2(\Omega )$ operator norm, we show in §5.2 that no reasonable measure $\tilde{\mu }$ can help.

   Note that in general $A$ is not necessarily a graphop in $L^2(\Omega, \tilde{\mu })$ as it might not even be symmetric in $L^2(\Omega, \tilde{\mu })$ . Thus, the numerical radius equivalence (23) is not guaranteed; however, the lower bound $r_{\tilde{\mu }}(A)\leq \|A\|_{2\to 2,\tilde{\mu }}$ still holds.

4. (iv) The above does not exclude the possibility that completely different approaches can lead to stability results for more general graphops, for example, a different entropy functional or a different structuring of the graph and position dependent interaction term in specific cases.

If we drop the goal of exponential decay, we can still infer that the steady state $\rho _\infty$ is not exponentially unstable.

Corollary 17. Let the assumptions and notation of Theorem 15 be given but let $A$ be an arbitrary graphop. Then, classical solutions $\rho$ to (3) are bounded for all times, that is,

(33) \begin{equation} \hat{H}(\rho (t)|\rho _\infty ) \le c,\quad t\geq 0 \end{equation}

with a constant $c\geq 0$ which depends on the equation coefficients and $\hat{H}(\rho _0|\rho _\infty ).$

Proof. The proof of Theorem15 yields the estimate:

(34) \begin{align} \begin{aligned} \frac{d}{dt} & \hat{H}(\rho |\rho _\infty ) \leq -\frac{4\pi ^2}{L^2} \hat{H}(\rho |\rho _\infty )\\ & + \kappa \|\Delta _x D \|_{L^\infty (U)} \int _{\Omega } (A\|\rho - \rho _\infty \|_{L^1(U)}) \|\rho - \rho _\infty \|_{L^1(U)} \,\mathrm{d}\mu (\xi ). \end{aligned} \end{align}

In the following estimation, we use the fact that the graphop $A$ satisfies $\|A\|_{\infty \to 1}\lt \infty$ , that $\|\rho (\xi ) - \rho _\infty \|_{L^1(U)} \leq 2$ for all $\xi \in \Omega$ due to $\rho (t,\cdot, \xi ),\rho _0\in \mathcal{P}_{\text{ac}}(U)$ and the CKP inequality (6):

(35) \begin{align} &\int _{\Omega } ( A\|\rho - \rho _\infty \|_{L^1(U)}) \|\rho - \rho _\infty \|_{L^1(U)} \,\mathrm{d}\mu (\xi )\nonumber \\ &\quad \leq \|A\|_{\infty \to 1} \sup _{\xi \in \Omega } \|\rho - \rho _\infty \|_{L^1(U)}\left (\int _{\Omega } \|\rho - \rho _\infty \|_{L^1(U)}\,\mathrm{d}\mu (\xi )\right ) \nonumber \\ &\quad \leq 2\|A\|_{\infty \to 1} \left (\int _{\Omega } \|\rho - \rho _\infty \|_{L^1(U)}\,\mathrm{d}\mu (\xi )\right ) \nonumber \\ &\quad \leq \sqrt{8}\|A\|_{\infty \to 1} \sqrt{\hat{H}(\rho |\rho _\infty )}. \end{align}

Plugging (35) into (34) yields

(36) \begin{equation} \frac{d}{dt} \hat{H}(\rho |\rho _\infty ) \leq -\frac{4\pi ^2}{L^2} \hat{H}(\rho |\rho _\infty ) + \sqrt{8}\kappa \|\Delta _x D \|_{L^\infty (U)}\|A\|_{\infty \to 1} \sqrt{\hat{H}(\rho |\rho _\infty )}. \end{equation}

Excluding the trivial case $\rho (t) = \rho _\infty$ , we denote $v(t)\,:\!=\, \sqrt{H(\rho (t)|\rho _\infty )} \gt 0$ and the non-negative constants $a\,:\!=\,\frac{4\pi ^2}{L^2}$ and $b\,:\!=\, \sqrt{8}\kappa \|\Delta _x D \|_{L^\infty (U)}\|A\|_{\infty \to 1}$ , we obtain the differential inequality:

\begin{equation*} \dot {v}(t) \leq -\frac {a}{2} v(t) + b, \quad t\geq 0. \end{equation*}

If for any $t\geq 0$ , the right-hand side is positive, and this is equivalent to the bound $v(t) \lt \frac{2b}{a}$ . For any other $t\geq 0$ , $v$ is non-increasing; thus, $v(t) \leq \max \{v(0), \frac{2b}{a}\}$ for all $t\geq 0$ . As $v(t)^2 = \hat{H}(\rho (t)|\rho _\infty )$ , the claimed result follows.

The result of Corollary 17 tells us that we are not too far from a global stability result for general graphops. This gives hope for an extension with an appropriate method.

3.4 Global stability for graphons

Let us now discuss the special cases of graphops in the form of integral operators with associated graphons (cf. Definition 6).

Definition 18. The graphon norm is defined as:

(37) \begin{equation} \|W\|_p\,:\!=\, \begin{cases} \left (\int _{\Omega }\int _{\Omega } W(\xi, \tilde{\xi })^{p}\,\mathrm{d}\tilde{\xi } \,\mathrm{d}\mu (\xi ) \right )^{\frac 1{p}},& p\in [1,\infty ),\\ \sup _{(\xi, \tilde{\xi }) \in \Omega ^2} |W(\xi, \tilde{\xi })|,& p=\infty . \end{cases} \end{equation}

Given a graphon $W$ that satisfies $\|W\|_p\lt \infty$ , it follows that for the associated graphop $A_W$ we have $A_W\in B_{p,p^*}(\Omega )$ with $\|A\|_{p\to p^*} \leq \|W\|_p$ . Additionally, for $p\gt 1$ the operator $A_W$ is compact. Indeed, the boundedness follows by applying Hölder’s inequality for $\frac 1p + \frac 1{p^*} = 1$ :

(38) \begin{align} \begin{aligned} \|A_W f \|^{p^*}_{p^*} &= \int _\Omega \left (\int _\Omega W(\xi, \tilde{\xi }) f(\tilde{\xi }) \,\mathrm{d}\mu (\tilde{\xi })\right )^{p^*} \mathrm{d}\mu (\xi )\\ &\leq \|W\|_{p^*}^{p^*} \left (\int _\Omega f(\tilde{\xi })^{p} \,\mathrm{d}\mu (\tilde{\xi }) \right )^{\frac{p^*}{p}}. \end{aligned} \end{align}

The compactness of $A_W$ for the case $\|W\|_{p^*}\lt \infty$ for $p^*\gt 1$ follows by finite rank approximation, for example, an approximation of $W$ with polynomials which are dense in $L^{p^*}(\Omega \times \Omega )$ .

Graphops that have a graphon density satisfying $\|W\|_{2}\lt \infty$ are specifically convenient to treat due to the underlying Hilbert space structure. As (positivity-preserving) Hilbert–Schmidt operators, their spectrum consists exclusively of the (positive) point spectrum, $\sigma (A)= \sigma _p(A)$ , with the only possible accumulation point at $0$ . Specifically, the spectral radius $\operatorname{rad}(A)$ of $A$ coincides with its numerical radius and bounds the operator norm:

\begin{equation*} \operatorname {rad}(A_W)\,:\!=\,\sup \{|\lambda | \mid \lambda \in \sigma (A_W)\} = \lambda _{\max }^{A_W} = n(A_W) \end{equation*}

where $\lambda _{\max }^{A_W}$ denotes the largest eigenvalue of $A_W$ .

Theorem15 directly implies the following:

Corollary 19. Let the assumptions and notations of Theorem 15 be given. Furthermore, let the graphop $A_W$ in (3) be associated with a graphon density $W\,:\,\Omega \times \Omega \to{\mathbb{R}}$ , such that $\|W\|_{2}\lt \infty$ . If the coupling coefficient satisfies

(39) \begin{equation} \kappa \lt \frac{2 \pi ^2}{L^2 \| \Delta _x D \|_{L^\infty (U)} \lambda _{\max }^{A_W} }, \end{equation}

then the solution $\rho$ to equation (3) is exponentially stable with the decay estimate:

(40) \begin{equation} \hat{H}(\rho (t)|\rho _\infty ) \leq \mathrm{e}^{-\hat{\alpha }(W) t} \hat{H}(\rho _0|\rho _\infty ), \quad t\geq 0, \end{equation}

where

\begin{equation*} \hat {\alpha }(W)\,:\!=\, \frac {4 \pi ^2}{L^2} - 2\kappa \| \Delta _x D\|_{L^\infty (U)}\lambda _{\max }^{A_W} \gt 0. \end{equation*}

Proof. The result follows directly from Theorem15 and (41).

Remark 20. As $A_W$ is symmetric in $L^2(\Omega )$ , it follows that

(41) \begin{equation} \lambda _{\max }^{A_W} = n(A_W) = \|A_W\|_{2\to 2} \leq \|W\|_2. \end{equation}

Hence, we can use $\lambda _{\max }^{A_W} \leq \|W\|_2$ in the estimates of Corollary 19, which leads to less optimal but more accessible coupling conditions and decay estimates.

4. Sakaguchi–Kuramoto model with frequency distribution

Let us consider global stability of the splay state for a well-studied variant of the McKean–Vlasov equation (2): The Sakaguchi–Kuramoto mean-field equation [4] with intrinsic frequency distribution $g$ . It is given as:

\begin{equation*} \begin {aligned} \partial _t \rho &= \partial _x(-\omega \rho + \kappa \rho V[A,g](\rho )) + \beta ^{-1}\partial _{xx} \rho, \quad t\geq 0,\\ \rho (0) &= \rho _0, \end {aligned} \end{equation*}

where a frequency-dependent transport term is added. The Vlasov term is dependent on the frequency density function $g$ via

(42) \begin{equation} V[A,g](\rho )\,:\!=\, \int _{\mathbb{R}} (\nabla D * A\rho ) \,g \mathrm{d}\omega . \end{equation}

We also explicitly include an inverse temperature parameter $\beta \gt 0$ to the diffusion term in order to discuss the limit $\beta \to 0+$ . In this case, solutions depend on $\rho (t,x,\xi, \omega )$ with $x\in [-\pi, \pi ]$ , $\xi \in \Omega$ , and where $\omega \in{\mathbb{R}}$ is the frequency variable. We consider frequency distributions according to a probability space $({\mathbb{R}},{\mathcal{B}(\mathbb{R})}, g(\omega )\mathrm{d}\omega )$ with an arbitrary density function $\|g\|_{L^1({\mathbb{R}})} = 1$ . Let us point out that replacing $g\mathrm{d}\omega$ with the Dirac distribution $\delta _0$ and choosing $\beta =1$ reduces the model again to (3).

4.1 Homogeneous case

For simplicity and in order to compare the result to the existing literature, let us first omit the additional network variable, that is, we set $A={\mathrm{id}}$ and consider initial data $\rho _0$ independent of $\xi$ . We assume that $\int _U \rho _0(x,\omega ) \,\mathrm{d} x = 1$ for all $\omega \in{\mathbb{R}}$ and as a result it holds that $\int _U \rho (t,x,\omega ) \,\mathrm{d} x = 1$ for all $\omega \in{\mathbb{R}}$ , $t\gt 0$ . Analogous to the extension to the entropy for graphop interactions (17), we can define the frequency-averaged relative entropy:

\begin{equation*} \overline {H}(\rho |\rho _\infty ) \,:\!=\, \int _{\mathbb {R}} \int _U \rho \log (\frac {\rho }{\rho _\infty }) \,\mathrm {d} x \,g\,\mathrm {d}\omega . \end{equation*}

For each $\omega \in{\mathbb{R}}$ , the transport term $\partial _t \rho = -\omega \partial _x(\rho )$ conserves the relative (spatial-)entropy $H(\rho ( \omega )|\rho _\infty )$ (as defined in (5)); hence, the computation is completely analogous to the proof of Theorem15. Specifically, the Vlasov term (42) corresponds to (4) with the all-to-all coupling graphop $A_g\rho \,:\!=\, \int _{\mathbb{R}} \rho (\omega ) g(\omega )\, \mathrm{d}\omega$ on the probability space $({\mathbb{R}},{\mathcal{B}(\mathbb{R})}, g\mathrm{d}\omega )$ . This amounts to a straightforward relabelling of the frequency variable $\omega$ as the network variable $\xi$ . With the only difference being that it is not required for $\mathrm{supp}\, g$ to be $\mathbb{R}$ , thus solutions and the steady state $\rho _\infty$ are naturally only relevant on the support of $g$ .

The graphop $A_g$ can be expressed with the graphon $W(\omega, \tilde{\omega }) = 1$ with $\|W\|_{L^2(g\mathrm{d}\omega )}=1$ . Thus, Corollary 19 (together with the spectral estimate (41)) yields $g$ -independent global stability for

(43) \begin{equation} \kappa \lt \frac{2\pi ^2}{L^2\beta \| \partial _{xx} D \|_{L^\infty (U)} } \,=\!:\, \kappa _0(\beta ), \end{equation}

then the solution $\rho$ to equation (3) is exponentially stable with the decay estimate:

\begin{equation*} \overline {H}(\rho (t)|\rho _\infty ) \leq \mathrm {e}^{-\alpha t} \overline {H}(\rho _0|\rho _\infty ), \quad t\geq 0, \end{equation*}

where the $g$ -independent decay rate is given as:

\begin{equation*} \alpha \,:\!=\, \frac {4 \pi ^2}{L^2\beta } - 2\kappa \| \partial _{xx} D\|_{L^\infty (U)} \gt 0. \end{equation*}

Using the CKP inequality (6) on the product space $L^1({\mathbb{R}}\times U, g\mathrm{d}\omega \times \mathrm{d} x ) \,=\!:\, L^1(g\mathrm{d}\omega \mathrm{d} x)$ , this means we have the following type of estimate:

\begin{equation*} \|\rho -\rho _\infty \|_{L^1(g\mathrm {d}\omega \mathrm {d} x )} \leq e^{-\frac {\alpha }{2} t} c,\quad t\geq 0, \end{equation*}

for some constant $c \gt 0$ .

Remark 21. Let us compare the stability estimates of (43) with the established stability results in the literature with the standard setting $L=2\pi$ and $D(x)= -\cos \!(x)$ . Sakaguchi [4] as well as Strogatz and Mirollo [13] have proven that the critical coupling strength – which marks the onset of synchronisation phenomena and the loss of the stability of the incoherent steady state – is given as:

\begin{equation*} \kappa _c(\beta ) \,:\!=\, 2\left [ \int _{\mathbb {R}} \frac {\beta ^{-1}}{\beta ^{-2} + \omega ^2} g(\omega )\,\mathrm {d}\omega \right ]^{-1}\stackrel {\beta \to \infty }{\longrightarrow } \frac {2}{\pi g(0)}. \end{equation*}

Using Hölder’s inequality, we see that

\begin{equation*} \kappa _c(\beta ) \geq 2\left [\int _{\mathbb {R}} g(\omega ) \,\mathrm {d}\omega \sup _{\omega \in {\mathbb {R}}} \frac {\beta ^{-1}}{\beta ^{-2} + \omega ^2} \right ]^{-1} = \frac {2}{\beta } \gt \frac {1}{2\beta } = \kappa _0(\beta ), \end{equation*}

where $\kappa _0(\beta )$ was defined in (43). Hence, while our results are independent of any specific choice of frequency distribution $g(\omega )\mathrm{d}\omega$ , the global stability results are not sharp. It stays below the critical value for all $\beta \gt 0$ and converges to $0$ as $\beta \to \infty$ .

4.2 Heterogeneous case

In order to consider the presence of both network structure and frequency distribution, we prove the following lemma formulated for general “combined graphops”.

Lemma 22. Let graphop $A_i\in \mathcal{B}_{2,2}(\Omega _i)$ and underlying probability space $(\Omega _i,\mathcal{A}_i,\mu _i)$ be given for $i=1,2$ with numerical radius $n_{\mu _i}(A_i)$ . Let $\rho _0(x,\cdot, \cdot ) \in L^\infty (\Omega _1\times \Omega _2)$ for each $x\in U$ . Consider McKean–Vlasov equations (3) with Vlasov terms of form:

(44) \begin{equation} V[A_1,A_2](\rho )(x,\xi _1,\xi _2)\,:\!=\, (\nabla D * A_1A_2\rho )(x,\xi _1,\xi _2). \end{equation}

Then solutions $\rho (t,x,\xi _1,\xi _2)$ with network variables $\xi _i\in \Omega _i$ for $i=1,2$ fulfil the global stability results of Theorem 15 , setting $\mu =\mu _1\times \mu _2$ and replacing the quantity $n(A)$ with $n_{\mu _1}(A_1)n_{\mu _2}(A_2)$ .

Proof. In this setting, the relative entropy (17) contains the product measure $\mu = \mu _1\times \mu _2$ . Following the proof of Theorem15, the only difference in estimating the time derivative of the relative entropy is in estimating the interaction term in (26). Denoting $f(t,\xi _1,\xi _2) \,:\!=\, \|\rho (t,\xi _1,\xi _2)-\rho _\infty \|_{L^1(U)}$ , we can estimate the second term as:

\begin{align*} \kappa & \int _{\Omega _2}\int _{\Omega _1}\int _{U} \rho (\Delta _x D \star (A_1A_2\rho ) ) \,\mathrm{d} x\,\mathrm{d}(\mu _1\times \mu _2)(\xi _1,\xi _2)\\ &\leq \kappa \|\Delta _x D \|_{L^\infty (U)} n_{\mu _1\times \mu _2}(A_1A_2) \|f\|^2_{L^2(\mu _1\times \mu _2)}\\ &= \kappa \|\Delta _x D \|_{L^\infty (U)} n_{\mu _1}(A_1)n_{\mu _2}(A_2) \|f\|^2_{L^2(\mu _1\times \mu _2)}. \end{align*}

The last equality can be validated via Fubini’s theorem, the self-adjointness of $A_i$ in $L^2(\Omega _i)$ and the fact that $A_1A_2f = A_2A_1f$ as each operator is linear, bounded in $L^2(\Omega _i)$ and solely acts on its distinct variable.

Remark 23. Note that “combined graphop” interactions are different from multiplex networks which would correspond to $A_1 + A_2$ where both graphops act on the same network variable. In the Sakaguchi–Kuramoto model (4), we have one arbitrary graph structure and one all-to-all frequency coupling. This can be interpreted as a graph structure on the average frequency.

With the realisation of §4.1 that intrinsic frequencies can be treated analogously to all-to-all coupling and with Lemma 22, we are now able to obtain a global stability result for Sakaguchi–Kuramoto models with heterogeneous network interactions.

Proposition 24. Consider the Sakaguchi–Kuramoto model (4) with arbitrary graphop $A\in B_{2,2}(\Omega )$ , probability space $(\Omega, \mathcal{A}, \mu )$ , and arbitrary frequency distribution $\|g\|_{L^1({\mathbb{R}})} = 1$ . Then the global stability result of Theorem 15 is fulfilled for the relative entropy:

(45) \begin{equation} \hat{H}_{g\times \mu }(\rho |\rho _\infty )\,:\!=\, \int _{\mathbb{R}}\int _{\Omega }\int _{U} \rho \log \left(\frac{\rho }{\rho _\infty }\right) \, \mathrm{d} x \,\mathrm{d} \mu (\xi )\,g\mathrm{d}\omega . \end{equation}

Proof. We apply Lemma 22 with $A_1\rho = A_g\rho = \int _{\mathbb{R}} \rho (\omega ) g(\omega )\,\mathrm{d}\omega$ , $\Omega _1 ={\mathbb{R}}$ , $\mathrm{d}\mu _1 = g\mathrm{d} \omega$ as defined in (42) and $A_2=A$ . Then, we obtain the global stability result of Theorem15 for the relative entropy $\hat{H}_{g\times \mu }(\rho |\rho _\infty )$ .

Remark 25. The attentive reader might have noticed that for the here presented results, we do not require the frequency distribution to be absolutely continuous with respect to the Lebesgue measure. All results work for general Borel probability measures. Nonetheless, we choose to adhere to the literature established notation with a density function $g$ for direct comparison.

The Sakaguchi–Kuramoto example has shown that the developed entropy method to prove global stability for the splay steady state is robust even when combining heterogeneous interactions with arbitrary intrinsic frequency distributions. Let us add that the versatility of entropy functionals has also recently been showcased by providing explicit stability estimates in the large coupling strength regime [39].

5. Graph examples

In this section, we consider solutions to (3) for explicit graph interaction structures and apply the established global stability results.

5.1 Spherical graphop

Let us consider the spherical graphop [32]. It is an operator defined as:

\begin{equation*} A\,:\, L^2({\mathbb {S}}^2, \mu ) \to L^2({\mathbb {S}}^2, \mu ), \quad (A\rho )(\xi )\,:\!=\, \int _{\xi ^\perp } \rho \,\mathrm {d}\nu _\xi (\tilde {\xi }), \end{equation*}

where ${\mathbb{S}}^2\,:\!=\, \{ \xi \in{\mathbb{R}}^3 : |\xi |_2 = 1\}$ , $\mu$ is the uniform probablity measure on ${\mathbb{S}}^2$ , and the integration takes place along the $\xi$ -equator, defined as $\xi ^ \perp \,:\!=\, \{\tilde{\xi } \in \mathbb{S}^2 \mid \xi ^T\tilde{\xi } = 0 \}$ . For each $\xi \in{\mathbb{S}}^2$ , the measure $\nu _\xi$ denotes the uniform probability measure on the (1-dim) submanifold $\xi ^\perp$ .

As no density function exists with respect to a $\xi$ -independent measure, this is a graphop that has no graphon representation. Furthermore, as the degree of each $\xi$ is not finite, it is also not representable as a graphing, but rather a more general graphop located “in-between” graphons and graphings [32]. For discussion on its induced finite network structure and resulting numerical stability estimates, we refer to [33].

Proposition 26. Consider the McKean–Vlasov equation (3) with a spherical graphop in the Vlasov term. Then solutions are globally stable, provided

(46) \begin{equation} \kappa \lt \frac{2 \pi ^2}{L^2 \| \Delta _x D \|_{L^\infty (U)} }, \end{equation}

with the decay rate estimate $\hat{\alpha }(A)\,:\!=\, \frac{4 \pi ^2}{L^2} - 2\kappa \| \Delta _x D\|_{L^\infty (U)}.$

Proof. As a preparation, consider $f\in L^2({\mathbb{S}}^2, \mu )$ , then, for each $\xi \in{\mathbb{S}}^2$ , we can find a (non-unique) parameterisation of $f|_{\xi ^\perp }$ as:

(47) \begin{equation} f_\xi (\tau ) \,:\!=\, f(\cos \!(\tau ) v_1^\xi + \sin\! (\tau ) v_2^\xi ),\quad \tau \in [0,2\pi ), \end{equation}

where the vectors $\{\xi, v_1^\xi, v_2^\xi \}$ form an orthonormal basis of ${\mathbb{R}}^3$ for each $\xi \in{\mathbb{S}}^2$ . Further, we can transform each element of ${\mathbb{S}}^2$ into spherical coordinates. We denote $\xi = \xi (\phi, \theta )$ with $\phi \in [0,\pi )$ and $\theta \in [0,2\pi )$ and transformation factor $|\sin\! (\theta )|$ . In order to apply Theorem15, we show $\|A\|_{2\to 2} = 1$ :

With the considerations from above, we have

\begin{align*} \|Af\|^2_{L^2({\mathbb{S}}^2)} &= \int _{{\mathbb{S}}^2} \left (\int _{\xi ^\perp } f\, \mathrm{d}\nu _\xi (\tilde{\xi })\right )^2 \mathrm{d}\mu (\xi ) = \int _{{\mathbb{S}}^2} \left (\frac{1}{2\pi }\int _0^{2\pi } f_\xi (\tau )\,\mathrm{d}\tau \right )^2 \mathrm{d}\mu (\xi )\\ &\leq \int _{{\mathbb{S}}^2} \frac{1}{2\pi }\int _0^{2\pi } f_\xi (\tau )^2\,\mathrm{d}\tau \,\mathrm{d}\mu (\xi ) \\ &= \frac{1}{2\pi } \int _0^{2\pi } \frac{1}{4\pi }\int _0^{\pi } \int _0^{2\pi } f_{\xi (\phi, \theta )}^2(\tau ) |\sin\! (\theta )|\,\mathrm{d}\tau \,\mathrm{d}\theta \,\mathrm{d}\phi \\ &= \frac{1}{2\pi } \int _0^{2\pi } \|f\|^2_{L^2({\mathbb{S}}^2)}\,\mathrm{d}\phi = \|f\|^2_{L^2({\mathbb{S}}^2)}. \end{align*}

For the inequality above, we used Cauchy–Schwarz. The second to last equality holds, as for any fixed $\phi \in [0,2\pi )$ the inner two integrations exactly integrate $f^2$ once over ${\mathbb{S}}^2$ . Self-adjointness can be shown in a similar way using the spherical coordinates and resulting symmetries. As one can validate that $A$ is a Markov graphop, that is, satisfying $A1_{{\mathbb{S}}^2} = 1_{{\mathbb{S}}^2}$ , it holds that $\|A\|_{2\to 2} = 1$ .

To estimate the convergence of solutions to the IVP (3) with a spherical graphop coupling, we can now directly apply Theorem15 with $n(A)\leq \|A\|_{2\to 2} \leq 1$ .

Remark 27. The stability result of Proposition 26 is identical to the case of all-to-all coupling given in Proposition 3. We point out that the spherical graphop describes rather sparse interactions compared to an all-to-all coupling. Thus, it is worth investigating whether incorporating additional graphop observables into the analysis could clarify or sharpen the stability results. However, this is beyond the scope of the current study. We refer to future research and first numerical considerations of [33] which show improved convergence compared to the all-to-all coupled case. It is worth noting that there is significant potential to develop reliable numerics by leveraging the particular graph structure of the dynamics [40].

5.2 Graphon examples

Erdös–Rényi random graph

For the start, let us mention Erdös–Rényi random graphs which take the simple graphon form $W(x,y) = p$ with $p\in (0,1)$ . Then, we can apply Corollary 19 for solutions to (3) with the Erdös–Rényi random graphop in the Vlasov interaction term. In particular, (40) provides the decay rate

\begin{equation*} \hat {\alpha }(W)\,:\!=\, \frac {4 \pi ^2}{L^2} - 2\kappa \| \Delta _x D\|_{L^\infty (U)} p \gt 0, \end{equation*}

for solutions in relative entropy, given that the interaction coefficient fulfils $\kappa \lt \frac{2 \pi ^2}{L^2 \| \Delta _x D \|_{L^\infty (U)} p }$ . When compared to the homogeneous interaction case this decay aligns with the intuition, as the mean-field interaction strength $\kappa$ is simply reduced to $\kappa p$ .

Power law random graphs

Let us now consider power law random graphs, as constructed in [23]. They are important examples of intermediately sparse graphs that correspond to unbounded $L^p$ functions which the standard $L^\infty$ graphon convergence theory for dense graphs cannot handle.

To introduce power law graphs, let a set of $[N]$ vertices be given with $N\in \mathbb{N}$ . For distinct indices $i,j\in [N]$ , $i\neq j$ , the vertices are connected, that is, $A^{ij}=1$ , with the probability:

\begin{equation*} p(i,j) = \min \{1, N^\beta (ij)^{-\alpha }\},\quad \alpha \in (0,1), \beta \in (2\alpha -1, 2\alpha ). \end{equation*}

This results in a superlinear expected number of edges and an expected edge density $N^{\beta - 2\alpha }$ . To construct the empirical graphon $W^{(N)}_{\alpha, \beta }(\xi, \tilde{\xi })$ and the finite particle interaction term (1), one has to include the rescaling factor $r_N = N^{\beta - 2\alpha }$ , see [23, 31]. Consequently, such graphs converge to the power law graphon for $N\to \infty$ (in the cut metric [23] or graphop action sense [32] which are equivalent in this case), denoted as:

(48) \begin{equation} W_\alpha (\xi, \tilde{\xi })\,:\!=\, (1-\alpha )^2 (\xi \tilde{\xi })^{-\alpha },\quad \xi, \tilde{\xi }\in [0,1],\alpha \in (0,1). \end{equation}

With the choice $\Omega = [0,1]$ and the Lebesgue measure $\mu = \lambda$ , we represent the associated graphop for $\alpha \in (0,1)$ as:

(49) \begin{align} \begin{aligned} A_\alpha \,:\!=\,A_{W_\alpha }\,:\, &\quad L^\infty ([0,1], \lambda ) \to L^1([0,1], \lambda )\\ &\quad f \mapsto (A_\alpha f)(\xi )\,:\!=\, \int _{[0,1]} W_\alpha (\xi, \tilde{\xi }) f(\tilde{\xi }) \,\mathrm{d}\tilde{\xi }. \end{aligned} \end{align}

An increase in $\alpha \in (0,1)$ leads to a stronger localisation of the power law graph around the origin $(0,0)$ . While $W_\alpha$ is unbounded, for each fixed $\alpha \in (0,1)$ it is an $L^p([0,1]^2)$ graphon for $p\in [1,\frac 1\alpha )$ . Due to the bound (38), the associated graphop can be extended to an operator:

(50) \begin{equation} A_\alpha \in B_{p^*,p}(\Omega ) \text{ for each } p \in [1,\textstyle \frac 1\alpha ) \end{equation}

with $\|A_\alpha \|_{p^*\to p} \leq \|W_\alpha \|_{p}$ .

In the case $\alpha \in (0,\frac 12)$ , it follows that the power law graphon satisfies $n(A) = \|A_\alpha \|_{2\to 2} \leq \|W_\alpha \|_{2} = \frac{(1-\alpha )^2}{(1-2 \alpha )}\lt \infty$ and Corollary 19 provides the decay:

\begin{equation*} \hat {H}(\rho (t)|\rho _\infty ) \leq \mathrm {e}^{-\hat {\alpha }(W_\alpha ) t} \hat {H}(\rho _0|\rho _\infty ),\quad t\geq 0, \end{equation*}

with $\hat{\alpha }(W_\alpha )\,:\!=\, \frac{4 \pi ^2}{L^2} - 2\kappa \| \Delta _x D\|_{L^\infty (U)} \frac{(1-\alpha )^2}{(1-2 \alpha )} \gt 0$ as long as the condition

\begin{equation*} \kappa \lt \frac {2 \pi ^2(1-2 \alpha )}{L^2 \| \Delta _x D \|_{L^\infty (U)} (1-\alpha )^2} \end{equation*}

is satisfied.

For fixed $\alpha \in [\frac 12,1)$ , it only follows that $\|A_\alpha \|_{p^*\to p} \leq \|A_\alpha \|_{2\to p}\lt \infty$ with $p\lt \frac 1\alpha \leq 2$ is guaranteed. In fact, as discussed in Section 3.3, this is not sufficient to prove exponential decay with our method and entropy $\hat{H}(\rho |\rho _\infty )$ .

Finally, even modifying the measure $\tilde{\mu }$ in the entropy, as defined in (31), cannot resolve the issues for power law graphops with $\alpha \gt \frac 12$ . We prove this in the following lemma by showing that the necessary estimate (30) (discussed in § 3.3) does not hold true for any $\tilde{\mu }$ with $\mathrm{supp}\ \tilde{\mu } = [0,1]$ .

Lemma 28. Let the probability space $([0,1], \mathcal{B}([0,1]), \lambda )$ with the power law graphop $A_\alpha$ , as defined in (49) for $\alpha \in (\frac 12,1)$ be given. Then no probability Borel measure $\tilde{\mu }$ with $\mathrm{supp} \tilde{\mu } = [0,1]$ exists such that

\begin{equation*} n_{\tilde {\mu }}(A_\alpha ) = \sup \{(A_\alpha f,f)_{L^2([0,1],\tilde {\mu })} \mid f \in L^\infty ([0,1],\tilde {\mu }), \|f\|_{L^2([0,1],\tilde {\mu })} = 1 \} \lt \infty . \end{equation*}

Proof. For any such measure $\tilde{\mu }$ and $\rho \in L^\infty ([0,1])$ , we have the following identity:

(51) \begin{equation} \int _{[0,1]} \rho (\xi ) (A_\alpha \rho )(\xi ) \,\mathrm{d}\tilde{\mu }(\xi ) = (1-\alpha )^2 \left (\int _{[0,1]} \xi ^{-\alpha } \rho (\xi )\, \mathrm{d}\tilde{\mu }(\xi )\right ) \left (\int _{[0,1]} \xi ^{-\alpha } \rho (\xi ) \,\mathrm{d}\xi \right ). \end{equation}

Now, for $n\in \mathbb{N}$ define the sequence:

\begin{equation*} \rho _n(\xi )\,:\!=\, \begin {cases} \tilde {\mu }([0,\frac {1}{n}])^{-\frac 12},& \xi \in [0,\frac {1}{n}],\\ 0,& \xi \in (\frac 1n,1]. \end {cases} \end{equation*}

Then $\rho _n \in L^\infty (\tilde{\mu })$ and $\|\rho _n\|_{L^2(\tilde{\mu })} = 1$ for all $n\in \mathbb{N}$ . In order for $n_{\tilde{\mu }}(A_\alpha )\lt \infty$ to hold, it is necessary that:

(52) \begin{equation} \int _{[0,1]} \rho _n(\xi ) (A \rho _n)(\xi )\,\mathrm{d}\tilde{\mu }(\xi ) = \tilde{\mu }([0,\frac 1n])^{-1} \int _{[0,\frac 1n]} \xi ^{-\alpha }\,\mathrm{d}\tilde{\mu }(\xi ) \,n^{\alpha - 1} \leq C \end{equation}

for all $n\in \mathbb{N}$ . But as $\xi ^{-\alpha } \geq n^{\alpha }$ for all $\xi \in [0,\frac 1n]$ , it follows that

\begin{equation*} \tilde {\mu }([0,\frac 1n])^{-1} \int _{[0,\frac 1n]} \xi ^{-\alpha } \,\mathrm {d}\tilde {\mu }(\xi ) \,n^{\alpha - 1} \geq n^{2\alpha -1}. \end{equation*}

Plugging this estimate back in (52) leads to the condition $n^{2\alpha -1} \leq C$ , which needs to be satisfied for all $n\in \mathbb{N}$ . But as $2\alpha -1 \gt 0$ by assumption, this cannot be satisfied.

Remark 29. We have discussed the limitations of the entropy method in § 3.3 . Here, we computed the limitations explicitly on the examples of power law graphops. In summary, the main culprit stems from the lack of regularisation mechanisms in the equation in the network variable. Solutions evolve only according to a spatial diffusion term and a spatial drift, which do not directly improve the regularity of solutions in the network variable over time. To achieve improved regularity through the equation’s dynamics, one could consider modifying the interaction term such that the two variables are coupled more directly.

Conclusion

In this work, we investigated graphop McKean–Vlasov equations as a mean-field formulation of a system of interacting stochastic McKean (or noisy Kuramoto-type) differential equations with diverse heterogeneous interaction patterns. The resulting solutions depend not only on space and time but also on an additional network variable $\xi$ . We have shown existence of classical solutions for each fixed $\xi$ and finite $\xi$ -averaged entropy. We further applied the entropy method to show global stability of the chaotic steady state, provided the graphop is bounded in $L^2$ and the dynamic’s interaction strength does not exceed a graphop-dependent threshold. Crucially, the method achieves explicit decay rates and can deal with graphs of dense, intermediately dense and sparse structures of unbounded degree. This has been demonstrated on various prototypical examples. We have also extended the results to the closely related Sakaguchi–Kuramoto model with frequency distribution and heterogeneous interactions, highlighting the method’s robustness with respect to model variations. We have discussed the limitations and demonstrated them explicitly on the examples of power law graphops. With the provided results, we hope to showcase a rather general graph limit theory of graphops which is well suited to be incorporated in established PDE methods.