Macroscopic Noisy Bounded Confidence Models with Distributed Radical Opinions

TL;DR

This paper analyzes a nonlinear Fokker-Planck equation modeling opinion dynamics with environmental noise and radical opinions, establishing mathematical properties, stability conditions, and numerical methods to understand opinion clustering and transitions.

Contribution

It provides the first rigorous analysis of the well-posedness, stability, and clustering behavior of macroscopic opinion models influenced by radicals and noise, with validated numerical schemes.

Findings

Proved well-posedness and existence of stationary solutions.

Derived explicit noise bounds for exponential convergence.

Developed numerical schemes for clustering analysis.

Abstract

In this article, we study the nonlinear Fokker-Planck (FP) equation that arises as a mean-field (macroscopic) approximation of bounded confidence opinion dynamics, where opinions are influenced by environmental noises and opinions of radicals (stubborn individuals). The distribution of radical opinions serves as an infinite-dimensional exogenous input to the FP equation, visibly influencing the steady opinion profile. We establish mathematical properties of the FP equation. In particular, we (i) show the well-posedness of the dynamic equation, (ii) provide existence result accompanied by a quantitative global estimate for the corresponding stationary solution, and (iii) establish an explicit lower bound on the noise level that guarantees exponential convergence of the dynamics to stationary state. Combining the results in (ii) and (iii) readily yields the input-output stability of the system for sufficiently large noises. Next, using Fourier analysis, the structure of opinion clusters under the uniform initial distribution is examined. Specifically, two numerical schemes for identification of order-disorder transition and characterization of initial clustering behavior are provided. The results of analysis are validated through several numerical simulations of the continuum-agent model (partial differential equation) and the corresponding discrete-agent model (interacting stochastic differential equations) for a particular distribution of radicals.

Equations

$$\dot{x}_{i}(t)=\frac{1}{N}\sum_{j=1}^{N}w\big{(}x_{j}(t)-x_{i}(t)\big{)},\quad w(\xi)=\begin{cases}\xi,\,|\xi|\leq R\\ 0,\,|\xi|>R.\end{cases}$$

$$\mathrm{d}x_{i}(t)=\frac{1}{N}\sum_{j=1}^{N}w\big{(}x_{j}(t)-x_{i}(t)\big{)}\mathrm{d}t+\sigma\mathrm{d}W_{i}(t),$$

$$\displaystyle\begin{array}[]{l}\rho_{t}=\big{[}\rho\ (w\star\rho)\big{]}_{x}+\frac{\sigma^{2}}{2}\rho_{xx},\quad t\geq 0,\ x\in\mathbb{R}.\end{array}$$

$$\displaystyle\mathrm{d}x_{i}(t)=\frac{1}{N}\sum_{j=1}^{N+N_{r}}w\big{(}x_{j}(t)-x_{i}(t)\big{)}\mathrm{d}t+\sigma\mathrm{d}W_{i}(t),\quad i=1,\ldots,N$$

$${\dot{x}}_i(t)=0,i=N+1,\dots ,N+N_r.$$

$$\displaystyle\begin{array}[]{l}\rho_{t}=\big{[}\rho\ (w\star(\rho+M\rho_{r}))\big{]}_{x}+\frac{\sigma^{2}}{2}\rho_{xx},\quad t\geq 0,\ x\in\mathbb{R}.\end{array}$$

$$\displaystyle\left\{\begin{array}[]{rcl}\rho_{t}=(\rho\ G_{\rho})_{x}+\frac{\sigma^{2}}{2}\rho_{xx}&\text{in}&\ \tilde{X}\times(0,T)\\ \rho(\cdot+2,t)=\rho(\cdot,t)&\text{on}&\ \partial\tilde{X}\times(0,T)\\ \rho(x,\cdot)=\rho_{0}(x)&\text{on}&\ \tilde{X}\times\{t=0\},\end{array}\right.$$

$$G_{\rho}(x,t):=w(x)\star\big{(}\rho(x,t)+M\rho_{r}(x)\big{)}.$$

$$\frac{{\sigma }^2}{2}{\rho }_{xx}+(\rho G_{\rho }{)}_x=0.$$

$${\sigma }_b^2:=\frac{4R}{\pi }\left(M+\frac{R}{\sqrt{3}}+2\right)\text{ and }{c}_b:=\frac{4R^{2}M}{\pi \sqrt{3}}.$$

$${\sigma }_s^2=\frac{4R(3+M)}{\pi }+\frac{4R^2}{\pi \sqrt{3}}\exp \left(\frac{8R(1+M)}{{\sigma }_s^2}\right).$$

$$\Vert \rho (t)-1{\Vert }_{L^2}\le \beta e^{-\lambda t}+\frac{1}{\eta }\Vert {\rho }_r{\Vert }_{L^2},$$

$${\int }_{\tilde{X}}\rho (x,t)\mathrm{d}x={\int }_{\tilde{X}}{\rho }_0(x)\mathrm{d}x=2,$$

$$\Vert \rho (t){\Vert }_{L^1}\ge {\int }_{\tilde{X}}\rho (x,t)\mathrm{d}x=2>0.$$

$$\Vert G_{\rho }{\Vert }_{L^{\infty }}\le R\left(\Vert \rho (t){\Vert }_{L^1}+2M\right).$$

$$\Vert G_{\rho }{\Vert }_{L^{\infty }}\le C\Vert \rho (t){\Vert }_{L^1}\le C\Vert \rho (t){\Vert }_{L^2}.$$

$$|G_{\rho }(x,t)|$$

$$\le \int |x-y|{\text{1}}_{|x-y|\le R}|\rho (y,t)+M{\rho }_r(y)|\mathrm{d}y$$

$$\le R{\int }_{\tilde{X}}|\rho (y,t)+M{\rho }_r(y)|\mathrm{d}y$$

$$\le R\left({\int }_{\tilde{X}}|\rho (y,t)|\mathrm{d}y+2M\right),$$

$$\Vert v{\Vert }_{L^p(\tilde{X})}\le \mu (\tilde{X}{)}^{\frac{1}{p}-\frac{1}{q}}\Vert v{\Vert }_{L^q(\tilde{X})},$$

$$\Vert (G_{\rho }{)}_x{\Vert }_{L^p}\le C_1\Vert \rho (t){\Vert }_{L^p}+C_2\Vert {\rho }_r{\Vert }_{L^p}\le C[\Vert {\rho }_r{\Vert }_{L^p}]\Vert \rho (t){\Vert }_{L^p}.$$

$$\Vert {\partial }_x^k{G}_{\rho }{\Vert }_{L^2}\le C[\Vert {\rho }_r{\Vert }_{H^{k-1}}]\Vert \rho (t){\Vert }_{H^{k-1}}.$$

$$\displaystyle\big{(}G_{\rho}(x,t)\big{)}_{x}=$$

$$=$$

$$\displaystyle+\int_{x-R}^{x+R}\big{(}\rho(y,t)+M\rho_{r}(y,t)\big{)}\ \mathrm{d}y,$$

$${\partial }_x^k{G}_{\rho }=$$

$$+{\partial }_x^{k-2}\rho (x+R,t)-{\partial }_x^{k-2}\rho (x-R,t)+M{\partial }_x^{k-2}{\rho }_r(x+R,t)-M{\partial }_x^{k-2}{\rho }_r(x-R,t).$$

$$\Vert {\partial }_x^k{G}_{\rho }{\Vert }_{L^2}$$

$$\le C\left(\Vert \rho (t){\Vert }_{H^{k-1}}+\Vert {\rho }_r{\Vert }_{H^{k-1}}\right)$$

$$\le C[\Vert {\rho }_r{\Vert }_{H^{k-1}}]\Vert \rho (t){\Vert }_{H^{k-1}},$$

$$\Vert \nu G_{\rho }{\Vert }_{H^k}\le C[\Vert {\rho }_r{\Vert }_{H^{k-1}}]\Vert \nu {\Vert }_{H^k}\Vert \rho (t){\Vert }_{H^{k-1}}.$$

$$\Vert \nu G_{\rho }{\Vert }_{H^k}\le C\left(\Vert \nu G_{\rho }{\Vert }_{L^2}+\Vert {\partial }_x^k(\nu G_{\rho }){\Vert }_{L^2}\right).$$

$$\Vert \nu G_{\rho }{\Vert }_{L^2}\le \Vert \nu {\Vert }_{L^2}\Vert G_{\rho }{\Vert }_{L^{\infty }}\le C\Vert \nu {\Vert }_{L^2}\Vert \rho (t){\Vert }_{L^2}\le C\Vert \nu {\Vert }_{H^k}\Vert \rho (t){\Vert }_{H^{k-1}},$$

Full text

Next, using Fourier analysis, the structure of opinion clusters under the uniform initial distribution is examined. Specifically, two numerical schemes for identification of order-disorder transition and characterization of initial clustering behavior are provided. The results of analysis are validated through several numerical simulations of the continuum-agent model (partial differential equation) and the corresponding discrete-agent model (interacting stochastic differential equations) for a particular distribution of radicals.

M.A.S. Kolarijani and P. Mohajerin Esfahani are with Delft Center for Systems and Control, Delft University of Technology, Delft, The Netherlands. Email: {M.A.SharifiKolarijani, P.MohajerinEsfahani}@tudelft.nl.

A.V. Proskurnikov is with Department of Electronics and Telecommunications, Politecnico di Torino, Turin, Italy and also with the Institute for Problems of Mechanical Engineering, Russian Academy of Sciences (IPME RAS). Email: anton.p.1982@ieee.org.

1. Introduction

Recent decades have witnessed enormous progress in study of complex systems and their system-theoretic properties [1, 2]. The main effort has been invested into study of “self-organization” and “spontaneous order” phenomena [3] that have inspired the development of synchronization and consensus theory [4, 5]. Paradoxically, these regular behaviors arising from local interactions between subsystems (agents, nodes) of a complex system are studied much better than various “irregular” dynamic effects such as persistent disagreement and clustering, exhibited by many real-world systems. Although some culprits of this asynchrony and dissent (e.g. symmetries and other special structures in the coupling mechanisms, exogenous forces acting on some nodes, heterogeneous dynamics of nodes, etc.) have been revealed in the literature [6, 7, 8, 9, 10], only a few mathematical models have been proposed that are sufficiently “rich” to capture the diversity of clustering behaviors in real-world networks and, at the same time, admit rigorous analysis. Long before the recent “boom” in complex systems, the lack of such models was realized in mathematical sociology. The problem of disclosing mechanisms preventing consensus and maintaining enduring disagreement between individuals [11] is nowadays referred to as the community cleavage problem or Abelson’s diversity puzzle [12, 13]. The interdisciplinary area of sociodynamical modeling [14, 15, 16, 13, 17, 18, 19, 20] has attracted enormous attention of the research community and is primarily concerned with mechanisms of opinion formation under social influence.

Only few models, proposed in the literature to describe opinion formation processes, have been secured by experimental evidence. Such models, however, play an important role and contribute, in various aspects, in comprehending complex systems’ behaviors such as birth, death and evolution of clusters in systems of interacting particles, and in developing algorithms for control of these behaviors. This explains explosion of interest in models of “opinion formation” in systems and control literature. From the control-theoretic prospect, most of these models are simply networks of interacting agents, obeying the first-order integrator model. However, the term “opinion” is now widespread and used to denote the scalar or multi-dimensional state of an agent, even if this state does not have a clear sociological interpretation111From the sociological viewpoint, opinions are cognitive orientations of individuals towards some objects or topics [13]. (belonging, e.g., to an abstract manifold [21]). The opinion is thus some value of interest, held by an agent and updated, based on displayed opinions of the other agents.

Nowadays, linear models of opinion dynamics, extending the classical French-DeGroot system in various directions (allowing, e.g., stubborn agents, asynchronous interactions and repulsion of opinions [22, 13, 17, 23]) have been thoroughly studied. These models are sufficient to explain consensus and disagreement in social groups, as well as formation of special opinion profiles (e.g., bimodal distributions, standing for opinion polarization), however, general mechanisms leading to emergence and destruction of unequal clusters are still far from being well understood. To explain them, more complicated nonlinear models have been proposed, mimicking some important features of social influence. One feature observed in social and biological systems is the homophily [24], or tendency of individuals to bond with similar ones. Homophily is related to biased assimilation [25] effects: individuals readily accept opinions consistent with their views and tend to dismiss and discount opinions contradicting to their own views. Mathematically, coupling between close opinions is stronger than that of distant opinions, which is modeled by introducing opinion-dependent influence weights. Although the possibility of such nonlinearities in opinion dynamics models was mentioned in the pioneering work [11], substantial progress has been primarily achieved in analysis of bounded confidence models proposed several decades later as extensions of the deterministic [26] and randomized gossip-based [27] consensus algorithms for multi-agent networks. Bounded confidence models stipulate that a social actor is insensitive to opinions beyond its bounded confidence set (usually, this set is an open or closed ball, centered at the actor’s own opinion), which makes the graph of interactions among the agents distance-dependent. A detailed survey of bounded confidence models and relevant mathematical results can be found in [18]. Bounded confidence models exhibit convergence of the opinions to some steady values, which can reach consensus or split into several disjoint clusters. If the state-dependent interaction graph of the system is symmetric, this follows from general properties of iterative averaging procedures, and can alternatively be proved by exploring a special Lyapunov function (“kinetic energy”) [28, 18, 29]. In the general case of asymmetric interaction graphs, such a convergence has been proved only in special situations [30, 29], but seems to be a generic behavior [31, 32, 30].

Opinions in real social groups, however, usually do not terminate at steady values yet oscillate, which is usually explained by two factors. The first reason explaining opinion fluctuation is exogenous influence, which can be interpreted as some “truth” available to some individuals [33] or a position shared by a group of close-minded opinion leaders or stubborn individulas (“radicals”) [34, 35, 36]. Important results on stability of the HK model with radicals and more general “inertial” bounded confidence models were obtained in [30]. Typically, the exogenous signal is supposed to change slowly compared to the opinion evolution and is thus replaced by a constant; the main concern is the dependence between the constant input and the resulting opinion profile. Numerical results, reported in [34, 35] demonstrate high sensitivity of the opinion clusters to the radical’s opinion and reveal some counter-intuitive effects, e.g., an increase in the number of radicals sometimes decreases the number of their followers. The second culprit of persistent opinion fluctuation is uncertainty in the opinion dynamics, usually modeled as a random drift of each opinion. The presence of a random excitation can be interpreted as “free will” and unpredictability of a human’s decision [37]; besides this, randomized opinion dynamics models are broadly adopted in statistical physics [38, 39, 40, 41] to study phase transitions in systems of interacting particles.

Even for the classical models from [26, 27], disclosing the relation between the initial and the terminal opinion profiles remains a challenging problem (including, e.g., the $2R$-conjecture [42, 43]). In presence of noise, the analysis becomes even more difficult; some progress in the study of the interplay between confidence range and noise level has been achieved in recent works [44, 45]. One of the important directions in analysis of bounded confidence models is examination of their asymptotic properties as the number of social actors becomes very large $N\to\infty$ and their individual opinions are replaced by infinitesimal “elements”. The arising macroscopic approximations of agent-based models describe the evolution of the distribution of opinion (usually supposed to have a density) and are referred to as density-based [46], continuum-agent [47, 48], Eulerian [49, 50], kinetic [51], hydrodynamical [28] or mean-field [43, 52] models of opinion formation. In the continuous-time situation, the density obeys a nonlinear Fokker-Planck (FP) equation. To study clustering behavior of the macroscopic bounded confidence models, efficient numerical methods have been proposed that are based on Fourier analysis [40, 53, 43].

From practical viewpoint, it is convenient to consider opinions staying in a predefined interval, e.g., $[0,1]$. The HK and Deffuant-Weisbuch (DW) models, as well as their continuous-time counterparts [18], imply that starting within the interval, opinions never escape from it. This property, however, is destroyed by arbitrarily small noises. To keep the opinions bounded, some “boundary conditions” are usually introduced. The absorbing boundary condition assumes that the opinions are saturated at the extreme values [math] and $1$ [40, 45]; an important result from [45] demonstrates that arbitrarily small noises in this situation destroy clusters and lead to approximate consensus (the maximal deviation of opinions is proportional to the noise level). More interesting are opinion dynamics with the periodic boundary condition, wrapping the interval $[0,1]$ into a circle. The opinion density on the circle corresponds to a $1$-periodic solution of the FP equation on the real line [53, 54, 43]. A disadvantage of the simple periodic boundary condition is the merging of two extreme opinion values [math] and $1$. To distinguish between these extreme opinions, we incorporate an “almost reflective” (precisely, an even $2$-periodic) boundary condition. Dealing with the macroscopic FP equation, the opinion density is then conveniently represented by an even $2$-periodic solution on the real line. This paper is primarily concerned with mathematical properties of such solutions.

Main Contributions. In this paper, we advance the theory of macroscopic modeling of bounded confidence dynamics. We consider a bounded confidence model with environmental noise which also includes radical opinions, which are not concentrated at a single point (as in [49, 34, 33]) but rather distributed over the interval $[0,1]$. The FP equation acquires an (infinite-dimensional) exogenous input, describing the density and total mass of the radical opinions. This setup allows us to consider the interplay between the noise and the distributed radicals concerning the behavior of the system. In particular, for the macroscopic FP equation,

• (i)

the criteria for the existence, uniqueness, and regularity of an even periodic solution are establish (Theorem 2.1);

• (ii)

the existence of stationary solution is studied and a global estimate is provided that bounds the deviation of the stationary state from the uniform distribution (Theorem 2.2);

• (iii)

a sufficient condition is presented for exponential convergence of the dynamics to stationary state (Theorem 2.3). Combining this result with the global estimate of item (ii) renders input-output stability of the system (Corollary 2.4).

Developing ideas from [40, 53, 43], we then use Fourier analysis to characterize the clustering behavior of the system under the uniform initial distribution. Specifically, two numerical schemes are presented to analyze

• (iv)

the interplay between the relative number (mass) of radical agents (with respect to normal agents) and the critical noise level for order-disorder transition (Section 5.2), and

• (v)

the impact of the noise and the radical opinions density on the number and timing of the initial clustering behavior (Section 5.3).

These schemes are then validated through several numerical simulations of the large-scale agent-based and macroscopic density-based models for a particular distributions of radical opinions.

The paper in organized as follows. Section 2 introduces the macroscopic opinion dynamics model in question. Here, we also present our main theoretical results regarding well-posedness and stability of the model. The next two sections are concerned with technical proofs of these results. Section 3 is devoted to the proofs of well-posedness of the dynamics. In Section 4, we examine the properties of the corresponding stationary equation and provide the technical proofs for theoretical results on stability of stationary state. In Section 5, two numerical schemes are presented for characterization of the clustering behavior of the model. These general schemes are then tested in Section 6 for a particular distribution of the radical opinions. These results are accompanied by numerical simulations of both the agent-based and the macroscopic models. Section 7 concludes the paper.

Notations. The convolution of two functions $f$ and $g$ is denoted by $f\star g=\int f(x)\ g(y-x)\ \mathrm{d}y$. We note that in our case one of the functions has a compact support, so the integral always exists. For a function $f(t,x)$ we use $f_{x}$ (respectively, $f_{t}$) to denote the derivatives with respect to (w.r.t.) $x$ (respectively, $t$), so that $f_{xx}$ is the second partial derivative w.r.t. $x$. We also use the notation $\partial^{i}_{x}f$ for the $i$-th order derivative w.r.t. $x$. Let $X=[0,1]$ and $\tilde{X}=[-1,1]$. We use $\mathcal{P}(X)$ to denote the the space of probability densities on $X$. That is, $\rho\in\mathcal{P}(X)$ if $\int_{X}\rho(x)\ \mathrm{d}x=1$ and $\rho(x)\geq 0$ for all $x\in X$. We also use $\mathcal{P}_{e}(\tilde{X})$ to denote the space of probability densities on $X$, extended evenly to $\tilde{X}$. That is, $\mathcal{P}_{e}(\tilde{X})$ is the space of all functions $\rho:\tilde{X}\rightarrow[0,\infty)$ such that $\int_{X}\rho(x)\ \mathrm{d}x=1$ and $\rho(x)=\rho(-x)\geq 0$ for all $x\in\tilde{X}$. $L^{p}(\tilde{X})$ denotes the Banach space of all measurable functions $f:\tilde{X}\rightarrow\mathbb{R}$ for which $\|f\|_{L^{p}(\tilde{X})}<\infty$. $H^{k}(\tilde{X})$ for $k\in\mathbb{N}$ is used to denote the Sobolev space $W^{k,2}(\tilde{X})$. We use the subscripts per (respectively, ep) to denote the closed subspace of periodic (respectively, even periodic) functions in the corresponding function space. We denote the dual space of $H^{1}_{per}(\tilde{X})$ by $H^{-1}_{per}(\tilde{X})$ and we use $\langle\cdot,\cdot\rangle$ to denote the corresponding paring of $H^{1}_{per}(\tilde{X})$ and $H^{-1}_{per}(\tilde{X})$. We use $\rightarrow$ and $\rightharpoonup$ to denote strong and weak convergences, respectively, in an appropriate Banach space. A brief overview of function spaces relevant to this study is provided in Appendix A.

2. Model Description and Main Theoretical Results

2.1. Macroscopic Model of Opinion Formation

The conventional bounded confidence model describes opinion formation process in a network of $N>1$ agents. All agents have the same confidence range $R>0$. Agent $i$’s opinion at time $t\geq 0$, denoted by $x_{i}(t)\in\mathbb{R}$, is (directly) influenced only by the opinions of agents $j$, such that $|x_{j}(t)-x_{j}(t)|\leq R$. One of the simplest continuous-time bounded confidence models is as follows [28]

[TABLE]

It can be shown [18] that the opinions obeying the model (1) always converge $x_{i}(t)\to x_{i}^{s}$ as $t\to\infty$, with $w(x_{i}^{s}-x_{j}^{s})=0$ for all $i,j$. This corresponds to either consensus ($x_{i}^{s}=x_{j}^{s}$ for all $i,j$) of the terminal opinions or their splitting into clusters, comprising one or several coincident opinions. In the latter situation, the distance between each two clusters is greater than $R$.

Dynamics of real opinions (as well as physical processes, portrayed by “opinion dynamics” models) often do not exhibit convergence to steady values, and the fluctuation of opinions persists. In order to capture this effect, random uncertainties can be introduced into the model mimicking “free will” and unpredictability of a human’s decision [37]. The simplest of these uncertainties is an additive random noise. The model (1) is then replaced by the system of nonlinear SDE

[TABLE]

where $W_{i}$ are independent standard Wiener processes and $\sigma>0$ characterizes the noise level.

Since the dynamics of a stochastic system (2) becomes quite complicated as the number of agents grows, the standard approach to examine it is the mean-field (or macroscopic) approximation, considering the opinion profile $(x_{i}(t))_{i=1}^{N}$ as a random sampling drawn from some (time-varying) probability distribution of the opinion. Precisely, it can be shown [55, 56, 57] that empirical distributions $N^{-1}\sum_{i=1}^{N}\delta_{x_{i}(t)}$ converge (in the weak sense) as $N\to\infty$ to a distribution, whose density $\rho(t,x)$ obeys the FP equation

[TABLE]

A natural extension of the bounded confidence dynamics allows the presence of $N_{r}\geq 1$ radicals (stubborn agents, zealots) that do not assimilate others’ opinions, however, influence them directly or indirectly. Typically, the radicals’ opinions are supposed to be constant (or changing very slowly compared to the opinion formation of “normal” agents). Indexing the “normal” individuals $1$ through $N$ and the radicals $(N+1)$ through $N+N_{r}$, the opinion dynamics becomes as follows

[TABLE]

Often it is supposed that the radicals share a common opinion $x_{i}\equiv T$ for $i=N+1,\ldots,N+N_{r}$, which may also be considered as some “truth” perceived by some individuals [33] or, more generally, an exogenous signal [34]. The ratio $M=N_{r}/N$ can be treated as the relative “weight” or “strength” of this external opinion. More generally, one can assume that the radicals’ opinions are spread over $\mathbb{R}$. Supposing that $N,N_{r}\to\infty$, the relative mass of the radicals $M$ remains constant, and their empirical distribution $N_{r}^{-1}\sum_{i=1}^{N_{r}}\delta_{x_{N+i}}$ converges (in the weak sense) to a distribution with sufficiently smooth density $\rho_{r}$, the density of the “normal” opinions obeys the modified FP equation as follows

[TABLE]

Note that the classical bounded confidence dynamics (1), being a special case of continuous-time consensus protocol, has an important property: the minimal and maximal opinions $\min_{i}x_{i}(t)$ and $\max_{i}x_{i}(t)$ are, respectively, non-decreasing and non-increasing. In particular, if the initial opinions are confined to some predefined interval, e.g., $x_{i}(0)\in[0,1]$, then one has $x_{i}(t)\in[0,1]$ for all $t\geq 0$. The additive noise leads to random drift of the opinion profile, thus destroying the latter important property. Since in practice bounded ranges of opinions are usually considered, the dynamics (2), (5) are usually complemented by boundary conditions [40], preventing the opinions from escaping from the predefined range.

A typical boundary condition is the periodic condition, where the opinion domain $[0,1]$ is wrapped on a circle of circumference $1$ (formally, replacing a real opinion value $x\in\mathbb{R}$ by its fractional part $\{x\}=x-\lfloor x\rfloor=x\mod 1$). A disadvantage of the periodic boundary condition is that there is no distinction between the extreme opinions [math] and $1$. In this paper, we address this issue by considering another type of boundary condition, which we call even $2$-periodic. Precisely, a real opinion $x\in\mathbb{R}$ is replaced by $f(x)$, where $f$ is an even $2$-periodic function, such that $f(x)=x$ on $[0,1]$ (and hence $f(x)=-x$ for $x\in[-1,0]$, $f(x)=2-x$ for $x\in[1,2]$ and so on). In other words, we first evenly extend the opinion domain $[0,1]$ into the interval $[-1,1]$ and then wrap it on a circle of circumference $2$ so that the extreme opinions [math] and $1$ correspond to the antipodes of this circle. We note that with this even $2$-periodic extension, the “effective” boundary condition experienced by the agents is an almost reflective one, that is, when an agent leaves the opinion domain from one end, it is reflected back into the domain from the same end. This is different from the behavior under simple periodic boundary condition where the agents leaving the domain form one end, enter the domain from the other end. However, the introduced boundary condition is “almost” reflective since the even extension causes some boundary effects: the influence of more extreme neighbors of opinion values in the $R$-neighborhood of extreme opinions [math] and $1$ is reinforced. This is due to the even extension which introduces more extreme “artificial” neighbors; see Fig. 1.

As discussed in [53, 54, 43], the FP equation (4) under the periodic conditions retains its validity, however, $\rho(t,x)$ is not a probability density on $\mathbb{R}$ but a $1$-periodic function $\rho(t,x+1)=\rho(t,x)\geq 0$, such that $\int_{0}^{1}\rho(t,x)\mathrm{d}x=1$ (that is, $\rho(t,\cdot)$ serves as a density on the interval $[0,1]$). Similarly, for the even $2$-periodic boundary condition, the equation (4) retains its validity when we replace the probability density $\rho(t,x)$ with an even $2$-periodic function, that is, $\rho(t,-x)=\rho(t,x)$ and $\rho(t,x+2)=\rho(t,x)$. On the interval $[0,1]$, the function $\rho(t,\cdot)$ again serves as a probability density: $\int_{0}^{1}\rho(t,x)\mathrm{d}x=1$. We also assume that the initial density $\rho_{0}(x)=\rho(0,x)$ and the density of radical opinions $\rho_{r}(x)$, defined on $[0,1]$, are extended (in the unique possible way) to even $2$-periodic functions on $\mathbb{R}$.

In this study, without loss of generality, we take $X=[0,1]$ and $\tilde{X}=[-1,1]$ to be the bounded opinion domain and its even extension, respectively. To summarize the discussion above, the macroscopic model for opinion dynamics considered in this study is fully described by the following PDE

[TABLE]

where

[TABLE]

Note that in (11), we are considering the dynamics over a finite time horizon $T$ for the sake of analysis, however, $T$ can be chosen arbitrarily large. We again emphasize that the initial density $\rho_{0}$ and the radical density $\rho_{r}$ are the unique even $2$-periodic extensions of the corresponding densities from $X$ to $\tilde{X}$. In essence, we are considering the same dynamics as in [54] with the extra requirement for $\rho_{0}$ (and the newly introduced density $\rho_{r}$) to be even. Finally, we note that [58] also provides a detailed treatment of this dynamics (without radicals) for a class of interaction potentials on a torus in higher dimensions.

2.2. Main Theoretical Results

To recapitulate, we are interested in even $2$-periodic solutions of PDE (11), where $\rho_{0}$ and $\rho_{r}$ are even $2$-periodic. A natural question arises whether the model is well-posed in the sense that every (sufficiently smooth) initial condition $\rho_{0}$ and input $\rho_{r}$ correspond to a unique solution. The affirmative answer is given in the following theorem.

Theorem 2.1** (Well-posedness of dynamics).**

Let the initial density of normal opinions and the radical opinions density satisfy $\rho_{0}\in H_{ep}^{3}(\tilde{X})\cap\mathcal{P}_{e}(\tilde{X})$ and $\rho_{r}\in H_{ep}^{2}(\tilde{X})\cap\mathcal{P}_{e}(\tilde{X})$, respectively. Then, PDE (11) has a unique, even, strictly positive, classical solution $\rho\in C^{1}(0,\infty;C_{ep}^{2}(\tilde{X}))$ such that $\rho(t)\in\mathcal{P}_{e}(\tilde{X})$ for all $t>0$.

This result implies that $\rho(t):=\rho(t,\cdot)$ is a (strictly positive) probability density on $X=[0,1]$ for all $t>0$, as required. For the autonomous systems (without radicals), [54, 58] provide a sufficient condition for exponential convergence of the dynamics towards uniform distribution $\rho=1$ as an equilibrium of the system. Unlike those studies, the uniform distribution is not an equilibrium of the model considered in this study. However, it is possible to extend this stability result to our model including the exogenous input, i.e., the radicals. To this end, we first consider the stationary equation corresponding to PDE (11) given by

[TABLE]

We are particularly interested in even stationary solutions $\rho^{s}\in\mathcal{P}_{e}(\tilde{X})$ of (13). Our next result characterizes the stationary state of the system.

Theorem 2.2** (Stationary behavior).**

Let $\rho_{r}\in H_{ep}^{1}(\tilde{X})\cap\mathcal{P}_{e}(\tilde{X})$ be the radical opinions density.

• •

Existence*: the stationary equation (13) has an even, strictly positive, classical solution $\rho^{s}\in C^{2}_{ep}(\tilde{X})\cap\mathcal{P}_{e}(\tilde{X})$.*

• •

Estimate*: for any $\eta>0$, if $\sigma^{2}>\sigma_{b}^{2}+\eta c_{b}$, then $\|\rho^{s}-1\|_{L^{2}}\leq\frac{1}{\eta}\|\rho_{r}\|_{L^{2}}$, where*

[TABLE]

Notice how the global estimate in Theorem 2.2 bounds the difference between the stationary solution and the uniform distribution. This result shows that, even in presence of radical opinions, the stationary solution can be made arbitrarily close to the uniform distribution by increasing the noise level beyond a minimum level $\sigma_{b}$. We note that the minimum noise level $\sigma_{b}$ is directly related to the confidence range $R$ and the relative mass $M$. Also, as the “energy” $M\|\rho_{r}\|_{L^{2}}$ of the radicals increases, in order to counteract their effect and keep the stationary profile in a somewhat uniform state, one must increase the noise level further beyond $\sigma_{b}$.

With this result in hand, we can now consider the asymptotic stability of stationary state. The next result provides a sufficient condition for exponential convergence of the dynamics to stationary state for arbitrary (and sufficiently smooth) initial density $\rho_{0}$ and radical density $\rho_{r}$.

Theorem 2.3** (Stability).**

Let $\rho_{0}\in H_{ep}^{3}(\tilde{X})\cap\mathcal{P}_{e}(\tilde{X})$ be the initial density of normal opinions and $\rho_{r}\in H_{ep}^{2}(\tilde{X})\cap\mathcal{P}_{e}(\tilde{X})$ be the radical opinions density. Also, let $\rho\in C^{1}(0,\infty;C_{ep}^{2}(\tilde{X}))$ with $\rho(t)\in\mathcal{P}_{e}(\tilde{X})$ be the solution to the dynamic equation (11). Then, $\rho(t)$ converges to a stationary state $\rho^{s}\in C^{2}_{ep}(\tilde{X})\cap\mathcal{P}_{e}(\tilde{X})$ exponentially in $L^{2}$ as $t\rightarrow\infty$ if $\sigma>\sigma_{s}$, where $\sigma_{s}>0$ uniquely solves

[TABLE]

An immediate result of Theorems 2.2 and 2.3 is that for sufficiently large noises, the dynamics will converge to a stationary state that can be made arbitrarily close to uniform distribution by increasing the noise level.

Corollary 2.4** (Input-output stability).**

For any $\eta>0$, if $\sigma^{2}>\max\{\sigma_{b}^{2}+\eta c_{b},\sigma_{s}^{2}\},$ where $\sigma_{b}$ and $c_{b}$ are defined in (14) and $\sigma_{s}>0$ uniquely solves (15), then it holds that

[TABLE]

where the constant $\beta>0$ depends on $\rho_{0}$ and $\rho_{r}$ and the convergence rate $\lambda>0$ depends on $\sigma$, $R$, and $M$.

Remark 2.5** (Connection to existing works).**

The stability result of Corollary 2.4 corresponds to the result reported in [54, Theorem 2.3] on global stability of uniform distribution $\rho=1$ for sufficiently large noises in the autonomous system without radicals. In particular, by setting $M=0$ in the estimate given in Theorem 2.2, one has $c_{b}=0$, hence $\rho^{s}=1$ is the unique stationary state of the system for $\sigma^{2}>\sigma_{b}^{2}=\frac{4R}{\pi}\left(2+R/\sqrt{3}\right)$. We note that $\sigma_{b}$ is the same minimum noise level given in [54, Theorem 2.3], taking into account a multiplicative factor of two due to the even extension considered in this study. However, direct application of Theorem 2.3 for stability of $\rho^{s}=1$ leads to a sufficient minimum noise level $\sigma_{s}>\sigma_{b}$. This is due to the fact that this result is based on conservative estimates for $\rho^{s}$. Indeed, if one incorporates the fact that $\rho^{s}=1$ and modifies some of the arguments provided in the proof of Theorem 2.3 in Section 4.3, then one can show that, in the absence of radical agents, the uniform distribution $\rho^{s}=1$ is also globally exponentially stable for $\sigma>\sigma_{b}$, reproducing the result of [54, Theorem 2.3].

Finaly, we note that, based on the results provided in [58], the input-output stability result of Corollary 2.4 can be generalized to multi-dimensional first-order stochastic interacting particle systems for a particular class of interaction potentials.

The next two sections are mainly concerned with the technical proofs of the theoretical results listed above.

3. Well-posedness of Dynamics

This section is devoted to the proof of Theorem 2.1 concerning the well-posedness of the dynamics (11). Throughout this section, all the norms are w.r.t. $\tilde{X}=[-1,1]$ (as opposed to $X=[0,1]$), unless indicated otherwise. We use $C,C_{0},C_{1},\ldots$ to represent a generic constant (depending on the model parameters) which actual values may change from line to line. In case these constants depend on a particular object of interest, say $\theta$, this dependence is explicitly indicated by $C[\theta]$.

Let us first note that because of periodicity, the mass is preserved in (11), that is,

[TABLE]

for all $t\geq 0$. In particular, we have

[TABLE]

We will be using this property in the sequel.

We start by presenting some useful estimates for the object $G_{\rho}$ defined in (12) that make it possible to extend the results provided by [54] to our model.

Lemma 3.1** (Estimates for $G_{\rho}$).**

Let $G_{\rho}$ be the function defined in (12) with $\rho_{r}\in\mathcal{P}_{e}(\tilde{X})$. If $\rho(t)\in L_{per}^{1}(\tilde{X})$, then

[TABLE]

If, moreover, $\|\rho(t)\|_{L^{1}}>0$ , then

[TABLE]

Proof.

Notice

[TABLE]

from which we can conclude the inequality (17). The first inequality in (18) then immediately follows from (17) and the assumption $\|\rho(t)\|_{L^{1}}>0$. For the second inequality in (18) notice that since $\tilde{X}$ is of finite measure $\mu(\tilde{X})=\mu([-1,1])=2$, for any measurable function $v$ we have

[TABLE]

where $1\leq p\leq q\leq\infty$. ∎

Using the estimate (18) in Lemma 3.1, one can follow similar arguments as in [54, Lemma 2.1] to show $\|\rho(t)\|_{L^{1}}=2$ and $\rho(t)\geq 0$ for all $t\geq 0$; see also [54, Corollary 2.2]. Specifically, assuming PDE (11) has a solution $\rho\in C^{1}(0,T;C_{per}^{2}(\tilde{X}))$, one can derive a priori estimate which in turn implies that the solution is non-negative so that $\rho(t)$ is a probability distribution on $X=[0,1]$ for all $t\geq 0$.

Lemma 3.2** (Estimates for $\partial_{x}^{k}G_{\rho}$).**

Let $G_{\rho}$ be the function defined in (12) with $\rho_{r}\in\mathcal{P}_{e}(\tilde{X})$.

• (i)

For $1\leq p\leq\infty$, if $\rho(t),\rho_{r}\in L_{per}^{p}(\tilde{X})$ with $\|\rho(t)\|_{L^{1}}>0$, then

[TABLE]

• (ii)

For $k\geq 2$, if $\rho(t),\rho_{r}\in H_{per}^{k-1}(\tilde{X})$ with $\|\rho(t)\|_{L^{1}}>0$, then

[TABLE]

Proof.

We have

[TABLE]

which leads to the first inequality in (20). Using the fact that $\|\rho(t)\|_{L^{2}}\geq C\ \|\rho(t)\|_{L^{1}}>0$ (see (19)), we have the second inequality in (20). Computing the higher order derivatives w.r.t. $x$, we obtain for $k\geq 2$

[TABLE]

Hence,

[TABLE]

where for the last inequality we used the fact that $\|\rho(t)\|_{H^{k-1}}\geq\|\rho(t)\|_{L^{2}}\geq C\ \|\rho(t)\|_{L^{1}}>0$. ∎

Lemma 3.3** (More estimates for $G_{\rho}$).**

Let $\nu\in H^{k}_{per}(\tilde{X})$, $\rho_{r}\in H^{k-1}_{per}(\tilde{X})\cap\mathcal{P}_{e}(X)$, and $\rho(t)\in H^{k-1}_{per}(\tilde{X})$ with $\|\rho(t)\|_{L^{1}}>0$. Then for $k\geq 2$

[TABLE]

Proof.

Notice

[TABLE]

For the first term on the right hand side (r.h.s.) of (24), we have

[TABLE]

where for the second inequality we used (18). Also, using Leibniz rule, for the second term on the r.h.s. of (24) we can write

[TABLE]

where for the last inequality we used Morrey’s inequality which implies

[TABLE]

Now, from (20) we have for $i=1$

[TABLE]

and from (21) we have for $i\geq 2$

[TABLE]

Putting all these estimates together while keeping only the highest Sobolev norms, we obtain

[TABLE]

∎

Remark 3.4** (Connection to existing works).**

Lemma 3.3 is an extension of [54, Proposition 4.1].

With these estimates in hand, we can follow the same arguments as in [54] to show well-posedness of the dynamics described by PDE (11).

Sketch of proof of Theorem 2.1.

Consider the following sequence of PDEs

[TABLE]

with smooth initial and radical distributions $\rho_{0},\rho_{r}\in C_{per}^{\infty}(\tilde{X})\cap\mathcal{P}_{e}(\tilde{X})$ for now. By standard results on linear parabolic PDEs [59, Chapter 7], there exists a sequence $\{\rho_{n}:\ n\geq 0\}$ in $C^{\infty}(0,T;C_{per}^{\infty}(\tilde{X}))$ that satisfies (28). Furthermore, using the estimate (18) in Lemma 3.1, one can follow the same procedure provided in [54, Proposition 3.1] to show $\|\rho_{n}(t)\|_{L^{1}}=\|\rho_{n}(0)\|_{L^{1}}=2$, and hence, $\rho_{n}(t)\geq 0$ for all $n\geq 1$ and $t\geq 0$; see also [54, Corollary 3.2].

Remark 3.5** (Evenness of $\rho_{n}$).**

One can use the evenness of $\rho_{0}$ and $\rho_{r}$ to show that the unique solutions $\rho_{n}$ to PDEs (28) are also even in $x$ for all $t\geq 0$. However, since this property will not be used for existence, uniqueness and regularity results provided below, we will postpone this argument to later when we deal with the evenness of the unique solution to PDE (11).

Existence with smooth data. Using Lemmas 3.1 and 3.2 and following a similar idea as in [54, Lemmas 3.5 and 3.7], we can obtain the following convergence results for a limiting object $\bar{\rho}$

[TABLE]

where $n_{k}$ denotes a subsequence. Moreover, we have the following estimate for $\{\rho_{n}:\ n\geq 1\}$ and $\bar{\rho}$

[TABLE]

We claim that $\bar{\rho}$ is the unique weak solution to (11). That is, $\bar{\rho}$ solves the weak formulation of (11) defined as

[TABLE]

for any $\eta\in L^{2}(0,T;H^{1}_{per}(\tilde{X}))$. To show this, we multiply (28) by $\eta$ with $n=n_{k}$ and integrate to obtain

[TABLE]

For the first two terms in (34), using convergence results (31) and (30), we have as $k\rightarrow\infty$

[TABLE]

and

[TABLE]

Also, the last term in (34) can be written as

[TABLE]

The limit of the first integral in (35) is zero as $k\rightarrow\infty$. Indeed, we know $G_{\rho_{n_{k}-1}}$ is bounded by the inequality (17) in Lemma 3.1, hence, $\eta_{x}G_{\rho_{n_{k}-1}}\in L^{2}(0,T;L^{2}_{per}(\tilde{X}))$. Also, (30) implies $\rho_{n_{k}}\rightharpoonup\bar{\rho}$ in $L^{2}(0,T;L^{2}_{per}(\tilde{X}))$. The limit of the second integral in (35) is also zero as $k\rightarrow\infty$. For this integral, we have

[TABLE]

where for the second inequality we used (32) and the fact that $|w\star(\rho_{n_{k}-1}-\bar{\rho})|\leq C\ \|\rho_{n_{k}-1}-\bar{\rho}\|_{L^{1}}$ by Lemma 12 (set $M=0$ in (17)). Now, notice $\|\rho_{n_{k}-1}-\bar{\rho}\|_{L^{1}(\tilde{X})}\leq\|\rho_{n_{k}-1}\|_{L^{1}(\tilde{X})}+\|\bar{\rho}\|_{L^{1}(\tilde{X})}\leq 4$. Hence,

[TABLE]

as $k\rightarrow\infty$ by the strong convergence (29). Putting all these results together, we see that $\bar{\rho}$ indeed satisfies the weak formulation (33).

To complete the existence result, we have to show $\bar{\rho}(x,0)=\rho_{0}(x)$. This condition makes sense since $\bar{\rho}\in C(0,T;L^{2}_{per}(\tilde{X}))$ by [54, Theorem 3.8] and the convergence results (30) and (31). Pick some $\eta\in C^{1}(0,T;H^{1}_{per}(\tilde{X}))$ with $\eta(T)=0$ and rewrite the weak formulation (33) as

[TABLE]

Similarly, since $\rho_{n_{k}}(x,0)=\rho_{0}(x)$, we have

[TABLE]

Let $k\rightarrow\infty$ in (37), so for arbitrary $\eta(x,0)$ we obtain from (37) and (36) that

[TABLE]

which implies $\bar{\rho}(x,0)=\rho_{0}(x)$.

Relaxed regularity on data. In order to relax regularity assumption on data to $\rho_{0},\rho_{r}\in L_{per}^{2}(\tilde{X})\cap\mathcal{P}_{e}(\tilde{X})$, we can use the mollified version of the distributions $\rho^{\epsilon}_{0}=\phi_{\epsilon}\star\rho_{0}$ and $\rho^{\epsilon}_{r}=\phi_{\epsilon}\star\rho_{r}$ with the standard positive mollifier $\phi_{\epsilon}$, follow the same procedure and take the limit $\epsilon\rightarrow 0$ at the end. See also [54, Theorem 3.12] for the details of this process.

Uniquness. Let $\xi=\bar{\rho}_{1}-\bar{\rho}_{2}$ where $\bar{\rho}_{1}$ and $\bar{\rho}_{2}$ are two weak solutions to (11) with $\rho_{0},\rho_{r}\in L_{per}^{2}(\tilde{X})\cap\mathcal{P}_{e}(\tilde{X})$. Then, for every $\eta\in L^{2}(0,T;H^{1}_{per}(\tilde{X}))$ we have

[TABLE]

We can rewrite the last integrand as

[TABLE]

to obtain

[TABLE]

Now, for the first integral on the r.h.s. of (38), we have

[TABLE]

where for the first inequality we used (17) in Lemma 3.1 and Cauchy-Schwarz inequality, and for the second inequality we used Young’s inequality. Similarly, for the the second integral on the r.h.s. of (38), we have

[TABLE]

where for the second inequality we used (32) and Lemma 3.1 (see (18) and (19)). Using (3) and (3) for (38) and setting $\eta=\xi$, we obtain

[TABLE]

By [54, Theorem 3.8], we know

[TABLE]

Thus, for all $T$, we have

[TABLE]

This implies, for a.e. $t\in[0,T]$

[TABLE]

Hence, by Grönwall’s inequality,

[TABLE]

This implies $\|\xi(t)\|_{L^{2}}=\|\bar{\rho}_{1}(t)-\bar{\rho}_{2}(t)\|_{L^{2}}=0$ since $\xi(0)=\rho_{0}-\rho_{0}=0$. Then, from continuity of $\bar{\rho}_{1}$ and $\bar{\rho}_{2}$ in time (by [54, Theorem 3.8]), we obtain uniqueness. That is, $\bar{\rho}_{1}=\bar{\rho}_{2}$ for all $t\in[0,T]$.

Regularity. Here, we first mollify the problem data $\rho_{0}$ and $\rho_{r}$ with the standard positive mollifier $\phi_{\epsilon}$ so that the solutions $\{\rho_{n}:\ n\geq 0\}$ to (28) are all smooth. This allows us to take derivatives of (28) to any order. We then take the limit $\epsilon\rightarrow\infty$ at the end. For simplicity, we omit the arguments for this last step and drop the subscript $\epsilon$.

Employing Lemma 3.3, we can extend the improved regularity results in space in [54] (Theorem 4.2). That is, for $\rho_{0}\in H_{per}^{k}(\tilde{X})\cap\mathcal{P}_{e}(\tilde{X})$ and $\rho_{r}\in H_{per}^{k-1}(\tilde{X})\cap\mathcal{P}_{e}(\tilde{X})$, we have

[TABLE]

Moreover, since $\rho_{r}$ is constant in time, we can also employ the results on improved regularity in time provided by [54, Theorem 4.3] for our model. This means, for $\rho_{0}\in H_{per}^{2k}(\tilde{X})\cap\mathcal{P}_{e}(\tilde{X})$ and $\rho_{r}\in L_{per}^{2}(\tilde{X})\cap\mathcal{P}_{e}(\tilde{X})$, we have for $i\leq k$

[TABLE]

and

[TABLE]

With these regularity results in space and time, we can derive the required regularity on the solution as stated in Theorem 2.1. Let $\rho_{0}\in H_{per}^{3}(\tilde{X})\cap\mathcal{P}_{e}(\tilde{X})$ and $\rho_{r}\in H_{per}^{2}(\tilde{X})\cap\mathcal{P}_{e}(\tilde{X})$ and $\bar{\rho}$ be the unique weak solution to PDE (11). By (41), we have $\bar{\rho}\in L^{\infty}(0,T;H^{3}_{per}(\tilde{X}))$. Hence, by Sobolev embedding theorem [60, Section 4.12], we have $\bar{\rho}(t)\in C_{per}^{2}(\bar{\tilde{X}})$ (after possibly being redefined on a set of measure zero). This gives the required regularity in space. Also, (42) and (43) imply that $\bar{\rho}_{t}\in L^{2}(0,T;H^{1}_{per}(\tilde{X}))$ and $\bar{\rho}_{tt}\in L^{2}(0,T;H^{-1}_{per}(\tilde{X}))$. Hence, by [54, Theorem 3.8], we have $\bar{\rho}_{t}\in C(0,T;L^{2}_{per}(\tilde{X}))$ (after possibly being redefined on a set of measure zero). This gives the required regularity in time. Putting these results together, we have $\bar{\rho}\in C^{1}(0,T;C_{per}^{2}(\tilde{X}))$.

Evenness. The evenness imposed on $\rho_{0}$ and $\rho_{r}$ implies that if $\rho(x,t)$ is a solution of (11), then $\rho(-x,t)$ is also a solution. Indeed, from (11) we obtain

[TABLE]

where for that last equality we used the fact that $w$ is an odd function. Then, assuming $\rho_{0}\in H_{ep}^{3}(\tilde{X})\cap\mathcal{P}_{e}(\tilde{X})$ and $\rho_{r}\in H_{ep}^{2}(\tilde{X})\cap\mathcal{P}_{e}(\tilde{X})$ (notice that $H_{ep}^{k}(\tilde{X})\subset H_{per}^{k}(\tilde{X})$), the uniqueness of the solution $\bar{\rho}\in C^{1}(0,T;C_{per}^{2}(\tilde{X}))$ to PDE (11) implies that the solution is even, that is, $\bar{\rho}\in C^{1}(0,T;C_{ep}^{2}(\tilde{X}))$.

Positivity. Using the same approach as in [58], we consider the following version of (11) in the unknown function $\rho$ with $\bar{\rho}$ being the non-negative weak solution

[TABLE]

This is a linear parabolic PDE with smooth and bounded coefficients (by Lemmas 3.1 and 3.2) for which $\bar{\rho}$ is a classic non-negative solution. Thus, by parabolic Harnack inequality [59, Section 7.1.4, Theorem 10], we have

[TABLE]

for $0<t_{1}<t_{2}<\infty$ and some positive constant $c$. Non-negativity of $\bar{\rho}(x,t)$ implies that $\inf_{x\in\tilde{X}}\bar{\rho}(x,t)$ and hence $\bar{\rho}(x,t)$ is strictly positive for all $t>0$. ∎

4. Stationary Behavior and Stability

4.1. Existence of Stationary Solution

This section mainly concerns the proof of existence result in Theorem 2.2 for stationary equation (13). All the norms in this section are w.r.t. $X=[0,1]$ (as opposed to $\tilde{X}=[-1,1]$), unless indicated otherwise. We note that norms on the even $2$-periodic spaces computed w.r.t. to $X$ and $\tilde{X}$ differ by a multiplicative constant, e.g., $\|u\|_{L^{p}(\tilde{X})}=2^{\frac{1}{p}}\|u\|_{L^{p}(X)}$. We again use $C,C_{0},C_{1},\ldots$ to represent a generic constant (depending on the model parameters) which actual values may change from line to line. In case these constants depend on a particular object of interest, say $\theta$, this dependence is explicitly indicated by $C[\theta]$.

Let us begin with providing a fixed point characterization of the solution to stationary equation (13). We note that, corresponding to the solution to dynamic equation (11), we are particularly interested in even solutions $\rho^{s}\in\mathcal{P}_{e}(\tilde{X})$ of stationary equation (13).

Lemma 4.1** (Fixed point characterization).**

$\rho^{s}\in C^{2}_{ep}(\tilde{X})\cap\mathcal{P}_{e}(\tilde{X})$* is a solution of stationary equation (13) if and only if $\rho^{s}$ is a fixed point of the operator $\mathcal{T}:\mathcal{P}_{e}(\tilde{X})\rightarrow\mathcal{P}_{e}(\tilde{X})$ defined by*

[TABLE]

where the constant $K$ is determined by the normalizing condition

[TABLE]

Proof.

The “if” part is clear since any fixed point $\rho^{s}\in C^{2}_{ep}(\tilde{X})$ satisfies the stationary equation (13). For the “only if” part, note that integrating (13) once, we have

[TABLE]

Now notice that we can set $C=0$ since we are interested in even solutions to (45). Indeed, from (45) we have

[TABLE]

Hence, for an even solution, we obtain

[TABLE]

where we used the fact that $w$ is an odd function. This implies $C=0$. Rearranging and integrating (45) once again, we have

[TABLE]

where the normalizing condition gives the constant $K$ as

[TABLE]

This completes the proof. ∎

This characterization allows us to use tools from operator theory, in particular, Schauder fixed point theorem to derive existence result for the stationary solution. Before that, we present some preliminary results for the operator $\mathcal{T}$.

Lemma 4.2** (Estimates for $\mathcal{T}$).**

Let $\mathcal{T}$ be the operator on $\mathcal{P}_{e}(\tilde{X})$ defined by (44).

• •

If $\rho,\rho_{r}\in\mathcal{P}_{e}(\tilde{X})$, then

[TABLE]

and

[TABLE]

• •

If $\rho,\rho_{r}\in L_{ep}^{2}(\tilde{X})\cap\mathcal{P}_{e}(\tilde{X})$, then

[TABLE]

• •

If $\rho,\rho_{r}\in H_{ep}^{k-2}(\tilde{X})\cap\mathcal{P}_{e}(\tilde{X})$, then for $k\geq 3$

[TABLE]

Proof.

From the definition (44) and inequality (17) in Lemma 3.1 we obtain

[TABLE]

which gives the estimate (47).

Now, observe

[TABLE]

Using (17) in Lemma 3.1 and (47), we obtain the inequality (48).

For the inequality (49), first notice

[TABLE]

where for the second inequality we used the fact that $\|\mathcal{T}\rho\|_{L^{2}}\leq C\ \|\mathcal{T}\rho\|_{L^{\infty}}$ is bounded by (47). Also, we have

[TABLE]

where for the second inequality we used (20) in Lemma 3.2 and the last inequality follows from the fact that $\|\rho\|_{L^{2}}\geq\|\rho\|_{L^{1}}>0$ (see (19)). Inserting this result in (51), we obtain inequality (49).

Similarly, for $k\geq 3$, we have (see (51))

[TABLE]

Now, notice

[TABLE]

where for the last inequality we used Lemma 3.3. Combining this result with (52), we derive a recursive inequality. Performing the recursive computations while keeping the highest Sobolev norms, we obtain

[TABLE]

Then, since $\|\rho\|^{i}_{H^{k-2}}\geq\|\rho\|_{L^{2}}\geq C\ \|\rho\|_{L^{1}}>0$, we can remove the constant $C_{0}$ and consider its effect in constants $C_{i}$. This gives the desired inequality (50). ∎

Proposition 4.3** (Lipschitz continuity of $\mathcal{T}$).**

Let $\mathcal{T}$ be the operator on $\mathcal{P}_{e}(\tilde{X})$ defined by (44) with $\rho_{r}\in\mathcal{P}_{e}(\tilde{X})$. Then $\mathcal{T}$ is Lipschitz continuous in $L^{p}$ for $1\leq p<\infty$ with Lipschitz constant

[TABLE]

Proof.

We use a similar argument to the one provided by [52]. Let $\rho_{1},\rho_{2}\in\mathcal{P}_{e}(\tilde{X})$. Using the estimate (47) in Lemma 4.2, we have for $1\leq p<\infty$

[TABLE]

where for the last inequality we used $\|\mathcal{T}\rho\|_{L^{1}(X)}=1$. Now, define

[TABLE]

and observe

[TABLE]

Similarly, we can write the normalizing constant $K_{1}$ as

[TABLE]

From (4.1), it follows

[TABLE]

and

[TABLE]

Hence,

[TABLE]

Using (4.1) and (56), we can rewrite (4.1) as

[TABLE]

Hence,

[TABLE]

Now, notice that

[TABLE]

(recall that norms are defined over $X$) and for $a>0$

[TABLE]

Thus, we have

[TABLE]

Combining (57) and (58), we obtain

[TABLE]

Finally, using (47) in Lemma 4.2 and the inequality (19) which relates norms over domains of finite measure, we have

[TABLE]

where the constant $L_{\mathcal{T}}$ is given by (53). ∎

With this preliminary results in hand, we next move on to the proof of existence result in Theorem 2.2.

Proof of Theorem 2.2 (Existence).

Following essentially the same argument as in [58, Theorem 2.3] and using Lemma 4.1, we can present the existence result for the stationary solution as the fixed point of the operator $\mathcal{T}$. First note that using the estimate (47) in Lemma 4.2, we have $\|\mathcal{T}\rho\|_{L^{2}}\leq C\ \|\mathcal{T}\rho\|_{L^{\infty}}\leq c$ for some positive constant $c$. Thus, for the purpose of finding the fixed points of $\mathcal{T}$, we can restrict $\mathcal{T}$ to act on the closed and convex set $E:=\{\rho\in L^{2}_{ep}(\tilde{X})\cap\mathcal{P}_{e}(\tilde{X}):\|\rho\|_{L^{2}}\leq c\}$. Now, notice that, using inequalities (47) and (48) in Lemma 4.2, we have for any $\rho\in E$

[TABLE]

for some constant $c^{\prime}>0$. That is, $\mathcal{T}(E)\subset E$ is uniformly bounded in $H^{1}_{ep}(\tilde{X})$. Thus, by the Rellich-Kondrachov compactness theorem [59, Section 5.7, Theorem 1], $\mathcal{T}(E)$ is precompact in $L^{2}_{ep}(\tilde{X})$. Since $E\subset L^{2}_{ep}(\tilde{X})$ is closed, this implies $\mathcal{T}(E)$ is also precompact in $E$. Also, $\mathcal{T}$ is Lipschitz continuous by Proposition 4.3. Hence, by Schauder fixed point theorem [59, Section 9.2.2, Theorem 3], it has a fixed point $\rho^{s}\in E$ which by (59) belongs to $H^{1}_{ep}(\tilde{X})$.

Regularity. The estimate (50) in Lemma (4.2) implies that if $\rho_{r}\in H^{k-2}_{ep}(\tilde{X})$, then the fixed point $\rho^{s}=\mathcal{T}\rho^{s}\in H^{k}_{ep}(\tilde{X})$. In particular, if $\rho_{r}\in H^{1}_{ep}(\tilde{X})$, then $\rho^{s}\in H^{3}_{ep}(\tilde{X})$. Hence, by Sobolev embedding theorem [60, Section 4.12], $\rho\in C^{2}_{ep}(\overline{\tilde{X}})$ (after possibly being redefined on a set of measure zero).

Positivity. The positivity of the fixed point directly follows from the representation (44). ∎

Remark 4.4** (Uniqueness).**

By Proposition 4.3, $\mathcal{T}$ is Lipschitz continuous in $L^{p}$ with Lipschitz constant $L_{\mathcal{T}}$ given by (53), and thus, is a contraction for $L_{\mathcal{T}}<1$. Hence, by Banach fixed-point theorem [59, Section 9.2.1, Theorem 1], $\mathcal{T}$ has a unique fixed point for $L_{\mathcal{T}}<1$. Setting $p=1$ in (53) gives the sufficient condition $\sigma^{2}>\frac{16R}{\ln 3}$ for uniqueness of stationary solution. This result corresponds to the sufficient condition provided in [52, Theorem 2].

Remark 4.5** (Semi-Gaussian clusters).**

For a highly concentrated radical opinion distribution with average opinion $A=\int_{X}x\ \rho_{r}(x)\ \mathrm{d}x$, we can provide an approximate solution to the stationary equation (13) as follows

[TABLE]

where $K$ is the normalizing constant (see Appendix B for the details). This result is an extension of the approximate solution provided by [43, Section 5.2]. In particular, one can reproduce the same result by setting $M=0$ and $A=0$. Equation (60) shows that for highly concentrated radicals the possible accumulation of normals around the average radical opinion $A$ in the stationary state is semi-Gaussian with variance $\frac{\sigma^{2}}{2(M+1)}$. Note that, as argued in [43], other clusters centered at opinion values other than $x=A$ may also exist. As long as these clusters are well-separated so that inter-cluster influences can be ignored, one can use the same approximation to derive a semi-Gaussian profile for the shape of these clusters (set $M=0$ and $A=x_{0}$ in (60) where $x_{0}$ denotes the center of the corresponding cluster). This analysis shows that $M$ affects the shape of the possible cluster formed at the average radical opinion $A$ in stationary state.

4.2. Global Estimate for Stationary Solution

This section is devoted to the proof of the estimate given in Theorem 2.2. In this section, all the norms are w.r.t. the domain $\tilde{X}=[-1,1]$, unless indicated otherwise.

Proof of Theorem 2.2 (Estimate).

Let $\psi=\rho^{s}-1$ so that $\int_{X}\psi(x)\ \mathrm{d}x=0$. From the stationary equation (13) we obtain

[TABLE]

where we used the fact that $w\star 1=0$. Next, we multiply this last equation by $\psi$ and integrate by part over $\tilde{X}$ to derive

[TABLE]

The extra terms are zero due to periodicity. Thus,

[TABLE]

Now, using inequality (17) in Lemma 3.1, we obtain

[TABLE]

Also, we have

[TABLE]

Hence,

[TABLE]

Using estimates (4.2) and (4.2), we can obtain form (61) (recall that uniform distribution is not an equilibrium of the system and hence $\|\psi_{x}\|_{L^{2}}\neq 0$)

[TABLE]

Now, since $\int_{X}\psi(x)\ \mathrm{d}x=0$, we can employ Poincaré inequality [59, Section 5.8.1, Theorem 1] to obtain $\|\psi\|_{L^{2}}\leq C\ \|\psi_{x}\|_{L^{2}}$. The optimal value for the Poincaré constant for $\tilde{X}=[-1,1]$ is $C=\frac{1}{\pi}$. Combining this result with inequality (4.2), we have

[TABLE]

Defining $\sigma_{b}$ and $c_{b}$ as in (14) gives the desired inequality $\|\psi\|_{L^{2}}\leq\frac{1}{\eta}\|\rho_{r}\|_{L^{2}}$, where $\eta=(\sigma^{2}-\sigma_{b}^{2})/c_{b}$. ∎

4.3. Stability of Stationary State

This section is devoted to the proof of Theorem 2.3 concerning stability of stationary state. All the norms in this subsection are w.r.t. the domain $\tilde{X}=[-1,1]$ (as opposed to $X=[0,1]$), unless indicated otherwise.

Proof of Theorem 2.3.

We follow similar arguments as the ones in [54], except we consider a general stationary state $\rho^{s}$ (instead of the uniform distribution considered in [54]). Let $\psi=\rho-\rho^{s}$ so that $\int_{X}\psi(x)\ \mathrm{d}x=0$. From the dynamic equation (11), we obtain

[TABLE]

where for the last equality we used the fact that $\rho^{s}$ is a solution to the stationary equation (13), that is,

[TABLE]

Multiplying (4.3) by $\psi$ and integrating by part over $\tilde{X}$ we obtain (the extra terms are zero due to periodicity)

[TABLE]

Now, from inequality (17) in Lemma 3.1, we have

[TABLE]

and

[TABLE]

Also, following a similar procedure as in (4.2) and (4.2) with $M=0$, we obtain $\|w\star\psi\|_{L^{2}}\leq 2R^{2}\|\psi\|_{L^{2}}/\sqrt{3}$. Finally, from (47) in Lemma 4.2, we have $\|\rho^{s}\|_{L^{\infty}}\leq\exp\left(8R(1+M)/\sigma^{2}\right)$. Using these estimates and the Young’s inequality we can rewrite (4.3) as

[TABLE]

Hence,

[TABLE]

Once again, since $\int_{X}\psi(x)\ \mathrm{d}x=0$, we can employ the Poincaré inequality [59] (Section 5.8.1, Theorem 1) $\|\psi\|_{L^{2}}\leq C\ \|\psi_{x}\|_{L^{2}}$ with optimal Poincaré constant $C=\frac{1}{\pi}$ to obtain

[TABLE]

Then, by Grönwall’s inequality, we have

[TABLE]

Now, notice that $\|\psi(0)\|_{L^{2}}\leq\|\rho_{0}\|_{L^{2}}+\|\rho^{s}\|_{L^{2}}$ is finite. Thus, if the constant factor in the exponential is negative, then $\|\psi(t)\|^{2}_{L^{2}}\rightarrow 0$ exponentially fast as $t\rightarrow\infty$. Negativity of the this constant factor corresponds to the condition $\sigma>\sigma_{s}$, where $\sigma_{s}$ solves (15). ∎

5. Characterization of Solution: Fourier Analysis

In this section, we exploit the periodic nature of the system and use Fourier analysis to study the behavior of the solution to the PDE (11) with uniform initial condition $\rho_{0}=1$. To this end, we derive a system of ordinary differential equations (ODEs) describing the evolution of Fourier coefficients of the normal opinion density $\rho$. Then, these ODEs are used for identification of the so-called order-disorder transition. In particular, a numerical scheme is presented for approximating the critical noise level at which this transition occurs. Moreover, we use these ODEs to provide another approximation scheme for characterizing the initial clustering behavior of the system including the number and the timing of possible clusters. These numerical schemes are in essence similar to the linear stability analysis previously employed by [37, 61, 40, 53, 43] for analysis of noisy bounded confidence models without radicals.

5.1. Fourier ODEs for Macroscopic Model

Notice that the set $\{\cos\left(\pi nx\right)\}^{\infty}_{n=0}$ is an orthogonal basis for the space $L^{2}_{ep}(\tilde{X})$ containing even $2$-periodic functions on $\tilde{X}=[-1,1]$. Then, the even $2$-periodic extension of the probability densities in the model allows us to consider the Fourier expansions of $\rho$ and $\rho_{r}$ in the form of

[TABLE]

By inserting the expansions (69) into (11) and setting the inner product of the residual with elements of the basis to zero (in other words, taking inverse Fourier transform), we can obtain a system of quadratic ODEs describing the evolution of Fourier coefficients $p_{n}(t)$. Considering the first frequency components $n=1,\ldots,N_{f}$ , these ODEs are expressed as

[TABLE]

where $p=(p_{1},p_{2},\ldots,p_{N_{f}})^{T}$. Note that for $n=0$, i.e., the constant term in the Fourier expansion, we obtain $\dot{p_{0}}=\ 0$. This is due to the periodic nature of the system that preserves the zeroth moment. The coefficients in (70) are given by

[TABLE]

where

[TABLE]

with $\operatorname{sinc}x=\frac{\sin x}{x}$. Also, recall that $q_{n},n\in\mathbb{N}$ are the Fourier coefficients of $\rho_{r}$.

Interestingly, one notices that the interaction between different frequency components in the quadratic terms is limited to those that are in a sense complements of each other. That is, each frequency $n$ of $\rho$ is affected by the frequency pairs $(n_{1},n_{2})$ such that either $n_{1}+n_{2}=n$ or $|n_{1}-n_{2}|=n$. This, in turn, leads to a particular structure for the matrix $Q_{n}$ in the quadratic terms. As expected, a similar behavior is seen in the linear terms: the effect of each frequency $k$ of $\rho$ on a given frequency $n$ of $\rho$ is modulated by the frequency components $n+k$ and $|n-k|$ of $\rho_{r}$.

5.2. Order-disorder Transition

A common behavior in noisy interactive particle systems is the order-disorder transition. For large values of $\sigma$, the effect of diffusion process can overcome the attracting forces among agents preventing the system from forming any cluster. This behavior has been analyzed and observed in several noisy bounded confidence models for opinion dynamics. Pineda et. al. used linear stability analysis in [37, 61] to compute the critical noise level above which the clustering behavior diappears for a modified version of Defuant model [27]. The same behavior was also reported in [39] for Defuant model. The same technique of linear stability analysis was used in [43, 53] to compute the critical noise level for noisy HK system similar to our model, except without radicals.

Here, we provide a method for approximating the critical noise level $\sigma_{c}$ at which the transition occurs. To this end, we linearize the systems at $t=0$ to obtain a system of linear ODEs expressed as

[TABLE]

The vector $c\in\mathbb{R}^{N_{f}}$ and the matrix $B\in\mathbb{R}^{N_{f}\times N_{f}}$ are defined accordingly using the objects $c_{n}$ and $b_{n}$ in (81). We emphasize that the linearization (83) is for a uniform initial condition, i.e., $p_{n}(0)=0$ for $n=1,\ldots,N_{f}$.

Looking at coefficients $c_{n}$ and $b_{n}$ in (81), we notice that the noise level $\sigma$ only appears in the diagonal entries of $B$ such that by increasing $\sigma$, these diagonal entries decrease. That is, for a large enough $\sigma$, all eigenvalues of $B$ are negative and the linearized system (83) is stable. This will be our first criterion for determining the critical noise level $\sigma_{c}$: the noise level above which all eigenvalues of $B$ are negative.

In order to consider the effect of the constant linear growth rates $c$ in (83), we further require the stationary values $\bar{p}_{n},n=1,\ldots,N_{f}$ of the linearized system (83) (i.e., the solution to the equation $c+B\bar{p}=0$) to be relatively small. In other words, taking the equilibrium of the linearized system $1+\sum_{n=0}^{N_{f}}\bar{p}_{n}\cos\left(\pi nx\right)$ as an approximation of the stationary state $\rho^{s}$, we require $\rho^{s}$ to be close to uniform distribution $\rho=1$ (representing disorder). Similar to the theoretical estimate of Theorem 2.2, we quantify this criterion by using Parseval’s identity and setting

[TABLE]

where the constant $\gamma>0$ determines the level of similarity between $\rho^{s}$ and uniform distribution. To sum up, we solve numerically for the minimum level of noise for which $B$ is Hurwitz and the inequality (84) holds.

5.3. Initial Clustering Behavior

For noises smaller than the critical noise level $\sigma_{c}$, we expect to see a clustering bahvior. In order to characterize the initial clustering behavior, we make use of the exponential growth rate $\gamma_{n}\coloneqq(b_{n})_{n}$ and linear growth rate $c_{n}$ given in (81). The proposed numerical method is as follows. We ignore the interactions between different frequencies in (70), that is, for each frequency $n=1,\ldots,N_{f}$, we consider the equation $\dot{p}_{n}=c_{n}+\gamma_{n}p_{n}$ with $p_{n}(0)=0$ (corresponding to uniform initial distribution) for initial evolution of the Fourier coefficient $p_{n}$. Then, for a given set of model parameters $(\sigma,R,M)$ and radical opinions density $\rho_{r}$, we numerically compute the dominant wave-number $n^{*}:=\operatorname*{\arg\!\max}_{n\in\mathbb{N}}\gamma_{n}$ with $\gamma_{n^{*}}>0$, that is, the unstable mode with the largest exponential growth rate. We speculate that the corresponding trigonometric term $p_{n^{*}}\cos(\pi n^{*}x)$ is the dominant component of the initial clustering behavior. The sign of $p_{n^{*}}$ depends on the linear growth rate $c_{n^{*}}$: $p_{n^{*}}>0$ if $c_{n^{*}}>0$, and $p_{n^{*}}<0$ otherwise.

Considering the even $2$-periodic extension of the model, the dominant wave-form must be interpreted on the interval $\tilde{X}=[-1,1]$. Then, the number of initial clusters $n_{\text{clu}}$ in the interval $X=[0,1]$ resulting from the wave-form $1+p_{n^{*}}\cos(\pi n^{*}x)$ is given by

[TABLE]

We also expect that the timing of this initial clustering behavior to be inversely related to the corresponding exponential growth rate $\gamma_{n^{*}}$. Indeed, by solving for the time for which the solution to the equation $\dot{p}_{n}=c_{n}+\gamma_{n}p_{n}$ is equal to $\pm 1$, we can approximate the time to initial clustering $t_{\text{clu}}$ as

[TABLE]

A similar approximation has been used in [53] in order to derive the time to the initial clustering using fluctuation theory.

6. Numerical Study

In this section, we provide a numerical study of the model at hand for a particular distribution of radical agents/opinions through simulations of the corresponding discrete- and continuum-agent models. Furthermore, we validate the result of Fourier analysis for identification of order-disorder transition (Section 5.2) and characterization of initial clustering behavior (Section 5.3).

The particular radical distribution considered in this section is a triangular distribution with average $A$ and width $2S$

[TABLE]

Although this choice may seem specific, it is rich enough for our purposes. In particular, with this choice, the zeroth, first and second moments of the radical opinions density are simply captured by the parameters $M$, $A$ and $S$, respectively. Moreover, we assume that the radicals are concentrated around their average opinion, that is, we consider small values of $S$ (w.r.t. the confidence range $R$).

For the discrete-agent model, the SDEs (5) are solved numerically using the Euler-Maruyama method for $N=500$ normal agents with time step $\Delta t=0.01$. In particular, for the radical agents, we produce a random sample of size $N_{r}=MN$ from the triangular distribution (91). The initial distribution of normal agents is taken to be uniform, that is, the initial opinions are randomly sampled from a uniform distribution on the interval $X=[0,1]$. For complete correspondence between the discrete- and continuum- agent models, we also consider the effect of even $2$-periodic extension in the simulations of the discrete-agent model. See Appendix C for the details of the numerical scheme. The details of numerical scheme for simulation of the continuum-agent model will be discussed in Section 6.1.

In the sequel, we make use of the order parameter

[TABLE]

introduced by [43] and its continuum counterpart

[TABLE]

to quantify orderedness in the clustering behavior of the model. In words, the order parameter $Q$ is the (normalized) number/mass of agents that are in $R$-neighborhood of and hence interacting with each other. In particular, in the continuum case, $Q_{c}=2R$ for a uniform distribution of opinions (absolute disorder), while $Q_{c}=1$ for a single-cluster distribution with all agents residing in an interval of width $R$ or less (complete order). In case of a clustered behavior, roughly speaking, the inverse of the order parameter is equal to the number of clusters. We also use the evolution of order parameter to characterize the timing of the clustering behavior.

In all the simulation results reported in this section the width of radicals distribution and the confidence range are fixed at $S=0.1$ and $R=0.1$, respectively.

6.1. Simulation of the Continuum-agent Model

In order to solve the continumm-agent model described by PDF (11) numerically, we use Fourier ODEs (70) to compute the coefficients of Fourier expansion of normal opinion density $\rho$ using the first $N_{f}$ terms of the expansion. However, regarding the radical opinion density, one notices that the considered triangular distribution does not satisfy the conditions of Theorem 2.1 for well-posedness of PDE (11), that is, $\rho_{r}\notin H^{2}_{ep}(\tilde{X})$. This will not be an issue since we will be working with the projection of the proposed $\rho_{r}$ in the Hilbert space $L^{2}_{ep}(\tilde{X})$. That is, we use the Fourier coefficients of $\rho_{r}$ in (70) which for the triangular distribution (91) are given by

[TABLE]

To be precise, we need the Fourier coefficients $q_{n}$ of $\rho_{r}$ for $1\leq n\leq 2N_{f}$, that is, twice the length of Fourier expansion of $\rho$; see the linear terms of (70). For the initial condition, we again consider uniform distribution $\rho_{0}=1$, which corresponds to $p_{0}=1$ and $p_{n}(0)=0$ for the Fourier coefficients.

It is also possible to employ a semi-explicit pseudo-spectral method, similar to the one provided by [43], for numerically solving (11); see Appendix D for details of this method for our model. The main difference is that the pseudo-spectral method solves the PDE for a set of discrete points in the space ($x\in X$) while solving the Fourier ODEs gives an approximation of the solution in terms of a finite basis for the corresponding Hilbert space.

These two methods (if both converge) result in the same solution. Fig. 2 compares the result of numerical simulations of the model using these two methods for a particular combination of system data. Note that, in these simulations, the number of points for the spacial discretization in pseudo-spectral method is twice the $N_{f}$ for Fourier ODEs so that the methods are compatible, i.e., both include the same set of frequency components. The left panel of Fig. 2 shows a similar result using these two methods for $N_{f}=32$ frequencies. However, as the number of frequencies considered in the simulations are decreased, we see that that the pseudo-spectral method starts to diverge while Fourier ODEs are still stable.

In the remainder of this section, we use the Fourier ODEs (70) with $N_{f}=128$ for numerical simulation of the continuum-agent model since they are computationally more efficient.

6.2. Order-disorder Transition

In this section, we numerically study the order-disorder transition in the model. In particular, we consider the effect of the relative mass of radicals $M$ on the critical noise level $\sigma_{c}$ at which this transition occurs. Furthermore, we use our simulation results to examine the approximation scheme presented in Section 5.2. In this regard, we note that the interplay between the confidence range $R$ and the critical noise level $\sigma_{c}$ have been studied in [43]. There, the authors showed that as $R$ increases, the critical noise level $\sigma_{c}$ also increases in such way that for small values of $R$, we observe a first-order transition.

6.2.1. An Illustrative Example

Our model exhibits the same order-disorder transition previously reported for similar noisy HK systems [40, 53, 43]. Fig. 3 shows this effect for a particular combination of system data in the discrete- and continuum-agent models. Notice that for $\sigma$ larger than a critical level the clustering behavior almost disappears (see the lower panel corresponding to $\sigma=0.05$ in Fig. 3(a)). To be more precise, a higher level of noise decreases the life-time of clustering behaviors with larger number of clusters.

This effect can be particularly seen in the evolution of order parameter in Fig. 3(b). In this regard, notice that for noises smaller than the critical noise level (here $\sigma<0.05$) the flat areas in the order parameter in Fig. 3(b) correspond to a clustered behavior, where the number of clusters is equal to the inverse of the order parameter. To illustrate, observe that for $\sigma=0.03$ and $\sigma=0.04$, the system reaches a single-cluster profile around the average radical opinion $A=0.7$. Notice, however, for $\sigma=0.03$ the system first goes through a 2-cluster profile corresponding to the flat area in the blue solid line at height $0.5$ in Fig. 3(b). On the other hand, for $\sigma=0.02$, we observe a 2-cluster profile at $t=10^{4}$ in Fig. 3(a). Notice, however, how the system goes through 4-cluster and 3-cluster profiles as depicted in Fig. 3(b) (the flat areas in the order parameter). Finally, for $\sigma=0.01$, we observe a very fast emergence of a 4-cluster profile (Fig. 3(b)) that has survived until $t=10^{4}$ as shown in Fig. 3(a). Here, we also notice that exact position of clusters in the discrete- and continuum-agent models differ. This particular difference between mean-field and agent-based models has been also mentioned by [37, 61]. Indeed, our numerical simulations show that even the number of clusters resulting from mean-field and agent-based models may differ; this also has been reported and explained previously in [43]. Finally, we note that for $M=0.1$, the approximation scheme explained in Section 5.2 results in $\sigma_{c}=0.043$ for $\gamma=1$ and $\sigma_{c}=0.051$ for $\gamma=0.1$ (see (84) for influence of $\gamma$).

6.2.2. Effect of $M$ on $\sigma_{c}$

Fig. 4 shows the order parameter derived numerically by simulating the continuum- and discrete-agent models. Notice how for each $M$, as noise increases, the system experiences a transition form order (clustered phase with $Q\approx 1$ in the yellow strip) to disorder (with $Q\approx 0.2$ in the dark blue area in the upper part of the plots). Also, we note that the blue strip in the lower part of plots in Fig. 4 represents clustering behaviors with larger number of clusters (similar to the behavior seen for $\sigma=0.01$ in Fig. 3).

This result shows that as the relative mass of radicals $M$ increases, the corresponding critical noise level $\sigma_{c}$, above which the system is in disordered state, also increases. The dependence of $\sigma_{c}$ on $M$ is in the form of a concave function. Furthermore, for small values of $M$, the transition seems to be discrete, signaling a first-order transition. However, for large values of $M$ the transition becomes blurry. This phenomenon was also reported in [43] for the dependence of the critical noise level on the confidence range $R$. Notice that as $M$ increases, the required noise level for disordered behavior also increases. This increase in the noise level leads to wider clusters, which, in turn, makes it difficult to differentiate order from disorder; see, e.g., the panels corresponding to $\sigma=0.04,0.05$ in Fig. 3.

Also shown in Fig. 4 (red lines) is the result of scheme provided in Section 5.2 for approximating the critical noise level. As can be seen, the scheme indeed provides a good approximation of the critical noise level. In particular, the dashed red line (for $\gamma=1$) almost perfectly separates the two phases of order and disorder.

6.3. Initial Clustering Behavior

For noises smaller than the critical noise level, agents start to form clusters; see Fig. 3. In particular, we observe a cluster of normal agents around the average radical opinion $A$ due to the force field generated by the radicals. Generally, three types of clusters may form: (1) the cluster at the average radical opinion $A$, (2) the cluster(s) at the extreme opinions $x=0$ and/or $x=1$, and (3) the cluster(s) around opinion values other than $x=0,1,A$. The third type of clusters are expected to perform a random walk with their center of mass moving like a Brownian motion (assuming clusters do not interact). The effective diffusivity of these Brownian motions is inversely related to the size of the cluster, i.e., the number of agents in the cluster. This will result in a process of consecutive merging between these clusters until complete disappearance of them. Detailed descriptions of this process are provided in [53, 43]. Notice however that this description does not apply to cluster(s) formed at $x=A$ and $x=0,1$. These clusters are affected by forces other than the normal attractions among the agents within the cluster. The cluster formed at $x=A$ is under influence of radicals and the possible clusters at the extreme opinions $x=0,1$ are reinforced due to the even $2$-periodic extension considered in our model. The behavior of these clusters (survival or dissolution) depends on their size, the exogenous force acting on them, and the effect of other clusters in their neighborhood.

In this section, we use the analysis scheme provided in Section 5.3 to investigate the effect of the zeroth and first moment of radicals ($M$ and $A$, respectively) on the initial clustering behavior of the model for noises smaller than the critical level. In particular, we investigate the effect of $M$ and $A$ on the number, position and timing of initial clusters for different values of $\sigma$. We again emphasize that we are considering a concentrated triangular distribution for radical agents and a uniform initial distribution for normal agents. Let us begin with illustrating how the objects introduced in Section 5.3, namely, exponential and linear growth rates and the dominant wave-number, can be used to characterize the initial clustering behavior.

6.3.1. An Illustrative Example

Consider the system data $(\sigma,M,A)=(0.01,0.1,0.7)$. Fig. 5 depicts the values of the exponential growth rate $\gamma_{n}$ and the linear growth rate $c_{n}$ for different frequencies. In Fig. 5(a) we observe that the unstable mode with the maximum exponential growth rate is $n^{*}=8$ with $\gamma_{n^{*}}=0.177$. Fig. 5(b) shows that the linear coefficient corresponding to this frequency is $c_{n^{*}}=0.007>0$. Then, (87) implies that the initial clustering behavior is expected to have $n_{\text{clu}}=5$ clusters. Also, using (88), we obtain $t_{\text{clu}}=18.16$ for the time to initial clustering.

Fig. 6 shows the time evolution of distribution of normal opinions/agents for the system data corresponding to Fig. 5. For the continuum-agent model, we can see a 5-cluster profile corresponding to the speculated waveform as depicted in Fig. 6(a). A similar clustering behavior is observed in the Monte Carlo simulation of the discrete-agent model in Fig. 6(b). Here, we observe three clear clusters: the cluster at average radical opinion $A=0.7$ and the two clusters at extreme opinions $x=0,1$. However, we observe an almost uniform distribution of normal agents in the opinion range $[0.1,0.5]$. This is due to the fact that the exact position of the corresponding clusters formed in the discrete-agent model varies within this range. Individual realizations of the discrete model show one, two or three clusters in this range with two clusters being the most frequent behavior as expected. This effect has been also reported by [37] in Monte Carlo simulations of a noisy Defuant model. Furthermore, we notice that the timing object $t^{*}=18.16$ also gives a good approximation for the onset of the corresponding clustering behavior for both continuum- and discrete-agent systems.

6.3.2. Effect of $M$ and $A$ on Initial Clustering

Performing a similar analysis to the one provided in the example above, we can compute the dominant wave-number ($n^{*}$), number of initial clusters ($n_{\text{clu}}$) and time to initial clustering ($t_{\text{clu}}$) for a general combination of system data. Fig. 7 shows the result of this analysis for different values of $M$ and $A$ at three different noise levels $\sigma$. Here, we only considered the values $A<1-R=0.9$ since for $1-R<A<1$ the boundary effect due to even $2$-periodic extension comes into play.

Comparing the left, middle and right panels of Fig. 7 corresponding to different levels of noise, we observe that as the level of noise increases, the number of clusters in the possible clustering behavior of the system decreases (see Fig. 7(b)), while the timing experiences a general increase (see Fig. 7(c)). This effect has been already shown in Fig. 3. In particular, with respect to the timing, we notice that as the level of noise decreases, the initial clustered profile emerges faster; see Fig. 3(b).

For low levels of noise, e.g., $\sigma=0.01$ (see the left panels in Fig. 7), the dominant wave-number does not depend on the $M$ or $A$. In this case, the most important effect of the first moment of radical opinions density $A$ is on the position of clusters. That is, the clustered profile emerges in a way that we observe a particular cluster formed at the average radical opinion $A$. The parameter $A$ also affects the timing of the clustering behavior in a periodic fashion. On the other hand, the zeroth moment of radical opinions density $M$ only affects the timing of the clustering behavior: as $M$ increases, $t_{\text{clu}}$ decreases. Fig. 8 shows the simulation results for $\sigma=0.01$ and compares the evolution of opinions for different values of $M$ and $A$. For the continuum model in the the top panels of Fig. 8 we observe that indeed a 4-cluster profile has emerged in all systems. Comparing Figs. 8(a) and 8(b) shows that $M$ only affects the timing of clustering behavior. This effect is better seen in Fig. 8(g) where we observe a faster convergence of order parameter for $\mathcal{S}_{2}$ with larger $M$. On the other hand, comparing Figs. 8(b) and 8(c) corresponding to $A=0.85$ and $A=0.7$, respectively, we observe a change in the positioning of the clusters. Monte Carlo simulations of the discrete-agent model reveals that the same general description also holds for this system. This is particularly seen in the time evolution of the order parameter in the discrete-agent model as depicted in Fig. 8(h). However, we once again note that there are differences between the behavior of the continuum- and discrete-agent models. In particular, the evolution of order parameter in Fig. 8(g) shows that the continuum-agent model has seemingly converged to steady-state with four clusters, while this is clearly not the case for the discrete-agent model as can be seen in Fig. 8(h). Indeed, in the discrete-agent model, as described in the beginning of this section, all the possible clusters formed around opinion values other than $x=0,1,A$ will necessarily disappear in the steady state profile, where the time required for their disappearance depends on the noise level and particularly the size of these clusters. Hence, unlike the discrete-agent model, for the continuum-agent model (in the limit $N\rightarrow\infty$), the system may require infinite time for this merging of the clusters to occure. This, in turn, can lead to differet behaviors in the discrete- and continuum-agent models over exponentially large times scales [43]; see also the evolution of order parameter in Fig. 3(b).

As shown in Fig. 7, for higher levels of noise, e.g., $\sigma=0.03$, we observe nonlinear effects. That is, $M$ and $A$ start to affect the dominant wave-number (see the middle and right panels of Fig. 7(a)). Nevertheless, these effects are limited as the number of clusters is still 3 or 4 for $\sigma=0.02$, and 2 or 3 for $\sigma=0.03$. Besides, we still observe a general increase in the timing of the clustering behavior as $M$ decreases. Fig. 9 shows the evolution of normal opinions/agents distribution and the corresponding order parameter for three different combinations of $M$ and $A$ at the noise level $\sigma=0.03$. Once again, in the continuum-agent model we observe a 2-cluster profile for all combinations as shown in the top panels of Fig. 9. For the discrete-agent model, we observe a 3-cluster behavior in which the cluster formed between the two clusters at $x=0$ and $x=A$ has already disappeared for $\mathcal{S}_{3}$ in Fig. 9(f) at $t=400$. Indeed, our simulations for $\sigma=0.03$ reveals a single-cluster profile around the average radical opinion $x=A$ after a large enough time; see Fig. 3.

To summarize the discussions above, for concentrated distribution of radicals, the main effect of the zeroth and first moments of radical distribution is on the timing and positioning of the possible clustering behavior, respectively. The number of clusters (to be precise, the life-time of possible transient clustered profiles) is mainly determined by the noise level of the system. This is particularly the case for lower levels of noise.

7. Final Remarks

In this paper, we considered a macroscopic model for bounded confidence opinion dynamics with environmental noise. In particular, we studied the effect of exogenous influence by adding a mass of radical (continuum) agents to the original population of the normal agents. The well-posedness of the continuum dynamics expressed as a nonlinear Fokker-Planck equation was established under some assumptions on the initial density of the normal opinions and the density of radical opinions. The long-term behavior of the model was also discussed by considering the corresponding stationary equation. In this regard, we provided a sufficient condition based on the noise level that guarantees exponential convergence of the dynamics towards the stationary state that can be made arbitrarily close to uniform distribution. In the context of opinion dynamics, we derived a theoretical bound on the minimum noise level required to counteract the effect of radical agents and keep the system in a somewhat uniform state.

Exploiting the periodicity of the considered continuum-agent model, we used Fourier analysis to provide a general framework for characterization of the clustering behavior of the system with uniform initial distribution. We then applied this framework for a particular distribution of radical opinions, namely, a relatively concentrated triangular distribution. In particular, we studied the effect of the relative mass of the radicals on the critical noise level for order-disorder transition. As expected, the analysis showed that for a larger number of radical agents, the critical noise level increases. We note that this result corresponds to the theoretical result on the global estimate for stationary state. However, comparing the theoretical lower bound on the noise level for the global estimate with its counterpart derived numerically, we find that the theoretical bound is quite conservative, which was expected considering its theoretical nature. We also considered the effect of relative mass and average opinion of radicals on the number, timing and positioning of the clusters for noises smaller than the critical noise level. Here, the noise level was shown to be the main factor in determining the number of clusters. Meanwhile, the relative mass of the radicals mainly affects the timing of the clustering behavior, that is, for larger masses of radicals, the clustering behavior is expected to emerge faster. On the other hand, the main effect of the average opinion of the radicals is on the positioning of the clusters; the clusters are positioned in a way that we see a cluster formed around the average opinions of radicals. The numerical simulations of the continuum-agent model and the corresponding discrete-agent model were in agreement with these results.

Appendix A Preliminaries on Function Spaces

The definitions provided here are mostly borrowed from [59]. Let $\{f_{k})\}_{k=1}^{\infty}$ be a sequence in a Banach space $B$ with norm $\|\cdot\|_{B}$. The strong convergence $f_{k}\rightarrow f$ implies $\|f_{k}-f\|_{B}\rightarrow 0$, while the weak convergence $f_{k}\rightharpoonup f$ implies $g(f_{k})\rightarrow g(f)$ for all bounded linear functionals $g:B\rightarrow\mathbb{R}$.

Let $f:\tilde{X}\rightarrow\mathbb{R}$ be a measurable function on $\tilde{X}=(-1,1)$. The $L^{p}$-norm of $f$ is defined as follows

[TABLE]

Then, $L^{p}(\tilde{X})$ denotes the Banach space of all measurable functions $f:\tilde{X}\rightarrow\mathbb{R}$ for which $\|f\|_{L^{p}(\tilde{X})}<\infty$.

Let $f,g\in L^{1}_{loc}(\tilde{X})$ be locally summable functions (i.e., $f,g$ have a finite integral over every compact subset of $\tilde{X}$). We say that $g$ is the $k$-th weak (partial) derivative of $f$, if

[TABLE]

for all test functions $\phi\in C^{\infty}_{c}(\tilde{X})$ (infinitely differentiable functions $\phi:\tilde{X}\rightarrow\mathbb{R}$ with compact support in $\tilde{X}$).

$H^{k}(\tilde{X})$ for $k\in\mathbb{N}$ is used to denote the Sobolev space $W^{k,2}(\tilde{X})$ consisting of functions $f\in L^{2}(\tilde{X})$ whose weak derivatives up to order $k$ exist and belong to $L^{2}(\tilde{X})$. Note that $H^{k}(\tilde{X})$ is a Hilbert space.

We use the subscript per to denote the closed subspace of periodic functions in the corresponding function space, e.g.,

[TABLE]

Similarly, we use the subscript ep to denote the closed subspace of even periodic functions in the corresponding function space, e.g.,

[TABLE]

We denote the dual space of $H^{1}_{per}(\tilde{X})$ by $H^{-1}_{per}(\tilde{X})$, that is, the space of bounded linear functionals on $H^{1}_{per}(\tilde{X})$. Moreover, we use $\langle\cdot,\cdot\rangle$ to denote the corresponding paring of $H^{1}_{per}(\tilde{X})$ and $H^{-1}_{per}(\tilde{X})$. That is, for $f\in H^{1}_{per}(\tilde{X})$ and $g\in H^{-1}_{per}(\tilde{X})$, we use $\langle g,f\rangle$ to denote the real number $g(f)$. Since periodic boundary condition allows for integration by parts without extra terms, $H^{-1}_{per}(\tilde{X})$ has most of the properties of the space $H^{-1}(\tilde{X})$, the dual space of $H^{-1}_{0}(\tilde{X})$; see [59, Section 5.9.1] for a detailed description of the space $H^{-1}(\tilde{X})$. In particular, one can extend the result in [59, Section 5.9, Theorem 3] to derive  [54, Theorem 3.8]. For reader’s convenience, the corresponding theorem is presented below.

Theorem**.**

[54, Theorem 3.8]** Let the function $f:\tilde{X}\times[0,T]\rightarrow\mathbb{R}$ be such that

[TABLE]

Then, $f\in C(0,T;L^{2}_{per}(\tilde{X}))$ after possibly being redefined on a set of measure zero. Moreover, the mapping $t\mapsto\|f(t)\|_{L^{2}(\tilde{X})}^{2}$ is absolutely continuous, with

[TABLE]

for almost every $t\in[0,T]$.

Appendix B Approximate Solution to Stationary Equation

In order to provide an approximate solution to the stationary equation (13), we assume radicals are highly concentrated around a particular opinion value $x=A$. To be precise, we assume that the average opinion of radicals is $A=\int_{X}x\ \rho_{r}(x)\ \mathrm{d}x$ and the variance of radicals $\sigma_{r}^{2}=\int_{X}(x-A)^{2}\ \rho_{r}(x)\ \mathrm{d}x$ is much smaller than the confidence range $R$. It helps to think of the limit being a point mass of radicals located at opinion value $x=A$. We further assume that the noise level $\sigma$ is also much smaller than $R$ so that the inter-cluster influences (from other possible clusters) can be ignored. Using these assumptions, we can expect this particular cluster of normal agents to be concentrated around $A$. This implies that in order to evaluate the integral in (46), we only need to consider values of $y$ near $A$. Under these assumptions, for $R<A<L-R$, we can write

[TABLE]

Inserting this result in (46), we have

[TABLE]

which can also be expressed as (by modefyying the normalizing constant $K$)

[TABLE]

In Fig. 3(a), this approximate solution is shown for $\sigma=0.03$ and $\sigma=0.04$, where the system has converged to a single cluster profile in both continuum- and discrete-agent models corresponding to the assumptions for derivation of the approximate solution. This result shows that the approximate solution is indeed a good approximation as it almost perfectly matches the numerical solution of the continuum-agent model.

Appendix C Euler-Maruyama Method for Discrete-agent Model

The interacting SDEs considered for the simulation of the discrete-agent model in this study is

[TABLE]

where $x_{i}^{ext},\ i=1,\ldots N$ are the opinions of normal agents and $x_{r_{i}}^{ext},\ i=1,\ldots N_{r}$ are the opinions of radical agents with $N_{r}=MN$. The superscript $ext$ corresponds to the even $2$-periodic extension as explained below. In order to solve (96) numerically, we employ the Euler-Maruyama method. Algorithm 1 summarizes the numerical scheme. As described in Section 6, we assume that the radicals have a triangular distribution centered at $A$ with width $2S$. That is, we produce a random sample of radicals with size $N_{r}$ from the triangular distribution (91) (Step 0). In particular, for complete correspondence between the discrete- and continuum- agent models, we also consider the effect of even $2$-periodic extension in our simulations. To this end, we use even $2$-periodic extensions of x and $\textbf{x}_{r}$ for calculating the sum on the r.h.s. of (96) (vectors denoted by $\textbf{x}^{ext}$ and $\textbf{x}_{r}^{ext}$ in Steps 0, 1 and 2). Also, due to periodicity, in each iteration, the opinion values outside the support $X=[0,1]$ are reflected back to $X$ (Step 5).

Appendix D Pseudo-spectral Method for Continuum-agent Model

This is a modification of the algorithm given in [43] for our model. Let us first recall the PDE (11) for the continuum-agent model

[TABLE]

where

[TABLE]

Using the first $N_{f}$ terms of Fourier expansions of $\rho$ and $\rho_{r}$, we can write

[TABLE]

Inserting this into (98), we obtain

[TABLE]

where $f_{k}$ is given by (82). Hence,

[TABLE]

With Fourier coefficients of $G$ in terms of Fourier coefficients of $\rho$ in hand, we can apply pseudo-spectral method for solving (97) as described in Algorithm 2. As shown, the multiplication $h=\rho\ G$ on the r.h.s. of (97) is performed in the time domain (Step 4), while the differentiations w.r.t. $x$ are performed in frequency domain (Step 5). Note that the symmetric nature of solution is preserved in the algorithm (Step 1). Also, preservation of mass is satisfied by setting $\hat{\rho}_{0}(t+1)=\hat{\rho}_{0}(t)$ (Step 5). We also note that the algorithm is semi-explicit (see the first equation in Step 5).