E. CalzolaDepartment of Computer Science, University of Verona, Strada Le Grazie 15, 37134 Verona, Italy. (elisa.calzola@univr.it) G. DimarcoDepartment of Mathematics and Computer Science & Center for Modeling, Computing and Statistics (CMCS), University of Ferrara, via Machiavelli 30, 44121 Ferrara, Italy. (giacomo.dimarco@unife.it) G. ToscaniInstitute of Applied Mathematics and Information Technologies, Via Ferrata 5, 27100 Pavia Italy.Department of Mathematics, University of Pavia, Via Ferrata 5, 27100 Pavia Italy. (giuseppe.toscani@unipv.it) M. ZanellaDepartment of Mathematics, University of Pavia, Via Ferrata 5, 27100 Pavia Italy. (mattia.zanella@unipv.it)

###### Abstract

In this work, we define a class of models to understand the impact of population size on opinion formation dynamics, a phenomenon usually related to group conformity. To this end, we introduce a new kinetic model in which the interaction frequency is weighted by the kinetic density. In the quasi-invariant regime, this model reduces to a Kaniadakis-Quarati-type equation with nonlinear drift, originally introduced for the dynamics of bosons in a spatially homogeneous setting. From the obtained PDE for the evolution of the opinion density, we determine the regime of parameters for which a critical mass exists and triggers blow-up of the solution. Therefore, the model is capable of describing strong conformity phenomena in cases where the total density of individuals holding a given opinion exceeds a fixed critical size. In the final part, several numerical experiments demonstrate the features of the introduced class of models and the related consensus effects.

Keywords: multi-agent systems, Bose-Einstein condensate, Boltzmann equation, nonlinear Fokker-Planck equation, opinion dynamics, conformity and consensus.

Mathematics Subject Classification: 35Q91, 91D30, 91B74,35Q84,82B40

## 1 Introduction

Opinion formation is a complex process influenced by several factors including social interactions, personal experiences, cultural background, and exposure to information. Understanding how opinions form and evolve is crucial in many fields in social sciences [44, 57]. Mathematical models are able to provide powerful tools for studying opinion formation, allowing researchers in other disciplines to simulate and analyze different scenarios and to gain insights into the underlying mechanisms at the basis of the opinion consensus phenomena [2, 34, 47, 52].

In this work, we concentrate on the case where opinion formation dynamics are influenced by group conformity, that is intended as the influence in opinion formation dynamics exercised by an high number of agents sharing a certain similar opinion. One possible way to understand group conformity is as a form of social pressure. Indeed, understanding the human tendency to *conform* to a larger group’s opinion has become a point of interest of many disciplines, particularly in the field of social and psychological sciences.These concepts are considered as relevant also in political choices and voting behaviors [20]. Already in 1958, Herbert Kelman [38] observed that there are different levels at which a change in attitude or actions of an individual can occur under the effect of social influence. These levels have different underlying processes leading the individual to accept the influence of the group because they want to cause a positive reaction from another person or group considered influential and/or avoid a punishment or disapproval.Social psychology works [7, 42] mainly focused their studies on two type of conformity: normative conformity, that occurs when a person behaves so to be accepted by a larger group, and informational conformity, that happens when individuals adapt their opinion or behavior to the one of people believed to have some sort of authority or accurate information. In more recent years, scientists started to investigate the biological reasons behind this, and it is worth noticing that conformity is not only a social phenomenon, it may also has a *biological* explanation. Neurosciences started investigating the brain’s reaction to conformity suggesting an emotional response of the participants to the group’s approval [11, 51], even generating a reward response [16].

To shed light on this fascinating process and to try to reproduce qualitatively the above described conformism behaviors, in this work we consider a new model in which the size of a group directly influences the dynamics of opinion formation. To this end, we will employ the methods of statistical physics, see e.g. [18].Recently, these techniques have proven to be highly effective and efficient tools for studying and analyzing opinion formation phenomena [2, 3, 5, 26, 29, 33, 49, 52] and related aspects [12, 13, 14, 22, 23, 24, 30, 39, 59]. A key area of research focuses on modeling information exchange, with particular attention to understanding and analyzing how interpersonal influence impacts processes like opinion formation and the creation of connections. This aspect is closely linked to developing graph models for large networks [15, 25, 28, 45]. In this regard, one of the key concepts from statistical physics that can be applied to the dynamics of opinion formation is the notion of emergent behavior [43], in which observable properties of large systems spontaneously arise from interaction forces defined at the microscopic level. In this direction, kinetic Boltzmann-type equations are capable to naturally link microscopic agent-based interactions with their large time collective properties, elucidating the processes involved in consensus formation.

In this work, we introduce a new microscopic mathematical description of opinion formation where each individual changes his own opinion as a consequence of the interactions with the whole population. The particular and innovative aspect of this interaction is that it is weighted with the intensity of the local density of agents sharing similar ideas. In our model, we also suppose that there is a certain amount of randomness in the interaction, modeling external factors which can be hardly controlled, such as the possibility to access information and the knowledge of every single individual. This microscopic dynamics is then inserted in a Boltzmann-type evolution equation and used to measure the variation in time of the density of individuals possessing an opinion on a given subject. Then, in order to try to get some insights on the behaviors of such dynamics we pass to the mean-field limit and we approximate the resulting Boltzmann equation through a grazing collision limit. This corresponds, from a modeling point of view, in assuming that the post-interaction opinions are very close to the pre-interaction ones, while at the same time the frequency of the interactions is suitably increased. The result obtained is a new Fokker-Planck-type model with nonlinear drift for the time evolution of the density of opinions, whose solution can present a blow-up to describe strong polarization phenomena due to the presence of group conformity. This blow-up appears only when the total density of individuals exceeds a fixed critical level and reveals strong analogies with some recent studies related to the Bose-Einstein condensate phenomenon [6, 48, 50, 53]. In this direction we mention the Fokker-Planck-type approach characterized by linear diffusion and superlinear drift introduced by Kaniadakis and Quarati [36, 37] for quantum indistinguishable particles.

In the second part of this work, we explore the new derived model from the numerical perspective point of view. Since both the microscopic as well as the Boltzmann dynamics are atypical, we first derive a new Monte Carlo method for the underlying model [8, 46] and then we numerically show that indeed the Boltzmann equation and its Fokker-Planck approximation produce analogous solutions when the scaling is chosen in such a way that the grazing collision limit is approached by the Boltzmann dynamics. We successively use this new numerical method to document the capability of the model in describing strong polarization of opinions when a singular ensemble of individuals is considered and finally we extend the model to the case of distinct ensembles of individuals. This last case is constructed in such way to mimic political parties competition and the possible migration of people from one party to another during the political competition. For this last situation, we also derive conditions which permit to forecast the formation of strong consensus among one party during the time evolution of the underlying dynamics.

The rest of the work is structured as follows. In Section 2 we present our new microscopic model of opinion formation and we derive the Boltzmann equation verified by the density of the opinions. In Section 3 we approximate the Boltzmann type equation with a Fokker-Planck model passing to the grazing collision limit. We also compute the explicit analytical expression of the steady state and recover the value of the critical mass separating the case where the steady state remains smooth and the case where the steady state exhibits a blow-up. In Section 4, we introduce a new Monte Carlo method and we use numerical simulations to validate the theoretical setting, showing in particular the resulting steady state density for various choices of the parameters and highlighting the differences in the shape of the distribution as the total mass increases or decreases, documenting blow-up situations. We also examine the case of several interacting populations of agents exchanging mass and opinions, as in the case of political factions interacting with different ideas on the same topic. We derive under which conditions one of the two populations reaches a critical mass and consequently a consensus emerges. In Section Conclusion and perspectives, we conclude the work with some future perspectives and open questions.

## 2 Microscopic interactions for group conformity and consensus

In this section, we propose a new model for opinion interactions within a population of homogeneous agents. This takes into account group conformity and leads to clusterization effects as discussed in the previous section. In recent years, a variety of heterogeneous social forces have been analyzed to understand, from a mathematical perspective, the process of consensus formation [2]. In this direction, several models for opinions, based on kinetic theory [18] have been introduced to describe the emergence of structures and patterns in real observed opinion dynamics, see e.g. [27, 52] and the references therein. Those models have been often based on the assumption that the agents’ opinion is modified by binary interactions with other individuals and possibly influenced by exogenous factors. In this direction, we mention recent works on the control of opinion formation dynamics [2], the case of interactions mediated by network topologies [1, 4, 40, 54] and multigroup dynamics, where agents may interact in a different fashion depending on the fact that they belong to a group or another, see e.g. [3, 26]. Apart from the above recalled mechanisms, in a realistic context, there are more complex phenomena involved like stubbornness, competence, exposure to information, confirmation bias, heterogeneous influences or group dynamics, among others. The latter aspect is in particular taken into account in this work.

In the following, we comply with an alternative path and we are specifically interested by modeling how individuals modify their opinions in relation to the number of people sharing a certain common idea. This process can be thought as a consequence of psychological factors influencing decision making processes towards behavioral and group conformity patterns [19, 35]. To that aim, let us start with considering a population of indistinguishable agents characterized by their opinion $x\in\mathcal{I}=[-1,1]$, where the extreme values $-1$ and $+1$ represent, respectively, a strongly negative and a highly positive opinion on a specific subject. For simplicity we suppose that each agent modifies his opinion through interaction with a fixed background distribution $g(z):\mathcal{I}\to\mathbb{R}_{+}$ with

$\int_{\mathcal{I}}zg(z)dz=m\in{(-1,1)},$ |

where $m$ is the mean opinion of the background and where this is supposed to mimic a distribution of opinion among a given population of individuals. We further introduce the distribution function $f=f(x,t):\mathcal{I}\times\mathbb{R}_{+}\to\mathbb{R}_{+}$ such that $f(x,t)dx$ represents the fraction of agents with opinion in $[x,x+dx]$ at time $t\geq 0$. Hence, if $(x,z)\in\mathcal{I}\times\mathcal{I}$ is a pair of pre-interaction opinions, the post-interaction opinions $(x^{\prime},z^{\prime})$ is obtainedby the following interaction scheme

$\begin{split}x^{\prime}&=x+\epsilon P(x)(z-x)+D[f](x,t)\xi,\\z^{\prime}&=z,\end{split}$ | (2.1) |

where $\epsilon\in(0,1/2]$, $P(\cdot)\in[0,1]$ is an interaction function and $D[\cdot]\geq 0$ represents the local relevance of the diffusion and it is such that $D[f](\pm 1,t)\equiv 0$ for all $t\geq 0$. In addition, $\xi$ is a random variable with moments

$\left\langle\xi\right\rangle=0,\qquad\textrm{and}\qquad\left\langle\xi^{2}%\right\rangle=\frac{\epsilon}{\lambda}>0,$ |

with $\lambda>0$, having denoted with $\left\langle\cdot\right\rangle$ the expectation with respect to the random variable $\xi$.The physical admissibility of the proposed interaction rule (2.1) is guaranteed if $|x^{\prime}|\leq 1$ for $|x|\leq 1$, meaning that, starting from an initial opinion $x\in\mathcal{I}$, the post-interaction opinion belongs to the $\mathcal{I}$. To this end, we can determine suitable bounds on the support of the random variable $\eta$, see also [21, 55, 56]. We start by observing that

$\begin{split}|x^{\prime}|&\leq(1-\epsilon P(x))x+\epsilon P(x)z+D[f](x,t)|\xi|%\\&\leq(1-\epsilon P(x))|x|+\epsilon P(x)+D[f](x,t)|\xi|,\end{split}$ |

since $|x|\leq 1$, from which we get that the sufficient condition to guarantee $|x^{\prime}|\leq 1$ is provided by

$D[f](x,t)|\xi|\leq(1-\epsilon P(x))(1-|x|).$ |

This condition is satisfied if a constant $c>0$ exists and is such that

$\begin{cases}|\xi|\leq c(1-\epsilon P(x))\\cD[f]\leq 1-|x|,\end{cases}$ | (2.2) |

for all $x,z\in\mathcal{I}$ and $f$. Hence, being $P\in[0,1]$ the first condition in (2.2) can be enforced by requiring

$|\xi|\leq c(1-\epsilon),$ |

and it is sufficient to control the support of the random variable $\xi$ such that $|\xi|\leq c(1-\epsilon)$. The second condition in (2.2) enforces $D[f](\pm 1,\cdot)=0$ and it is needed to avoid opinions to exit from their support. In the linear case, where $D[f](x,t)=D(x)$, several choices for such local diffusion function have been analyzed in the past [52]. One possible form is $D(x)=1-|x|$, which enforces the choice $c=1$ and therefore the bound on the support of $|\xi|\leq 1-\epsilon$. Another possibility is $D(x)=1-x^{2}$, which enforces the choice $c=\frac{1}{2}$, and therefore $|\xi|\leq\frac{1}{2}(1-\epsilon)$.

Afterward, we will consider a different situation and we will deal with the following non linear case

$D[f](x,t)=\sqrt{\dfrac{H(x)}{1+\beta(H(x)f(x,t))^{\alpha}}}$ | (2.3) |

where $H(\cdot)\geq 0$ is a concave function such that $H(\pm 1)\equiv 0$, and where the parameters $\alpha,\beta\geq 0$ are introduced to weight the impact of the group conformity in the self-thinking process. In more details, within the introduced choice, around the state $x\in\mathcal{I}$ self-thinking forces are mitigated by condensation effects which enforce compromise forces. In the following, we will concentrate on the following choices for the function $H(\cdot)$, weighting the local relevance of the diffusion: $H(x)={1-x^{2}}$ and $H(x)=(1-x^{2})^{2}$, see [52] for further details. In (2.3), the parameter $\beta>0$ tunes the influence of nonlinear effects. In the following, we will determine the regime in which an opinion condensation occurs for sufficiently large values of $\alpha>0$. Following now the analyses presented in [47] we can observe that

$\left\langle x^{\prime}+z^{\prime}\right\rangle=x+z+\epsilon P(x)(z-x),$ |

and consequently in such a model the mean opinion is not conserved on average unless $x\equiv z$ as expected since the distribution of $z$ is supposed not to vary with $x$. Furthermore, for $\epsilon\ll 1$ we have

$\left\langle(x^{\prime})^{2}+(z^{\prime})^{2}\right\rangle=x^{2}+z^{2}+2%\epsilon xP(x)(z-x)+\dfrac{D^{2}[f](x,t)}{\lambda}+o(\epsilon),$ |

and also the energy, playing the role of the variance with respect to the mean opinion, is not conserved on average in a single transition not even for $x\equiv z$. Looking now at the intensity of the transition we get

$\left\langle x^{\prime}-x\right\rangle=\epsilon P(x)(z-x)$ | (2.4) |

and

$\left\langle(x^{\prime}-x)^{2}\right\rangle=\epsilon^{2}P^{2}(x)(z-x)^{2}+%\epsilon\dfrac{D^{2}[f](x,t)}{\lambda}.$ | (2.5) |

The result of such process is that the compromise tendency becomes dominating in (2.1). The new rule for the change in opinions defined in (2.1) can be now encoded in a Boltzmann-type kinetic equation which, in weak form, writes

$\begin{split}&\dfrac{d}{dt}\int_{\mathcal{I}}\varphi(x)f(x,t)dx=\\&=\int_{\mathcal{I}\times\mathcal{I}}B_{\epsilon}[f](x,t)\left\langle\varphi(x%^{\prime})-\varphi(x)\right\rangle f(x,t)g(z)dx\,dz,\end{split}$ | (2.6) |

where $B_{\epsilon}[f](x,t)>0$ is a collision kernel. This equation describes the time evolution of the density of individuals sharing a given opinion where the ideas change in time with the law (2.1). The scope of the term $B_{\epsilon}[f](x,t)$ is to introduce a variable frequency enforcing larger interactions to highly populated opinions. The specific form chosen is the following

$B_{\epsilon}[f](x,t)=1+\beta(H(x)\min\left\{f(x,t),\epsilon^{-\alpha}\right\})%^{\alpha}.$ |

This is introduced with the aim of modeling the fact that opinions with larger consensus find typically more space in discussions in particular on social platforms and media channels having consequently more possibilities to influence the population. In other words, the model rewards discussions where the leading general opinion is present, while comparatively disregards interactions in which the leading opinion is not present.

###### Remark 2.1

One possible alternative way to model the effects of a strong exposure to leading information and social pressure in shaping the ideas of individuals could consist in modifying the interaction dynamics proposed in (2.1) in the following way

$x^{\prime}=\left\{\begin{array}[]{ll}\min\{x+(z-x)(1+\beta(H(x)f)^{\alpha})+H(%x)\xi,1\},&z>x\\\max\{x+(z-x)(1+\beta(H(x)f)^{\alpha})+H(x)\xi,{-}1\},&z<x.\\\end{array}\right.$ | (2.7) |

This different microscopic dynamics condenses the effects of all complex mechanisms involved in the formation of opinions through a mean field type of action where an agent align himself with the mean opinion expressed by the distribution of $g(z)$, i.e. the background fixed distribution of the society’s idea. This alignment, as before, can be due to the fact that a person finds out about a topic over the net, reading newspapers, watching news, through social media connections or as a result of face-to-face discussions. To this effect, we add, as opposite to the previous interaction, a linear random variable which models the unpredictable role played by the environment and/or personal situation of the single individual such as beliefs, knowledge or influence. Finally, one weights the interactions with the size of the local density of people sharing a similar idea to characterize the propensity to conformity and the social pressure which leads individuals to align. Such type of microscopic dynamics leads to similar mesoscopic models as the one detailed in (2.6).

## 3 Nonlinear Fokker-Planck kinetic models for group conformity

To get some insights about the model introduced in (2.6) with microscopic opinion updates defined in (2.1), here we develop an expansion techniques to reduce the complexity of the kinetic integro-differential equation (2.6). The idea consists in obtaining a Fokker-Planck-type partial differential equations in suitable regimes of parameters identified by the modeling choices which permits to mimic this slow gradual process of the formation of ideas and beliefs. This derivation has roots in the grazing collision limit [31, 58] and has been further explored in a variety of applications in socio-economic and life sciences in the recent past, see e.g. [32, 47].

Let $\varphi(x)$ be a test smooth function. From (2.4)-(2.5) we get by Taylor expansion

$\begin{split}&\left\langle\varphi(x^{\prime})-\varphi(x)\right\rangle=\dfrac{d%}{dx}\varphi(x)\left\langle x^{\prime}-x\right\rangle\\&\qquad+\frac{1}{2}\dfrac{d^{2}}{dx^{2}}\varphi(x)\left\langle(x^{\prime}-x)^{%2}\right\rangle+\dfrac{1}{6}\dfrac{d^{3}}{dx^{3}}\varphi(\tilde{x})\left%\langle(x^{\prime}-x)^{3}\right\rangle,\end{split}$ | (3.8) |

where $\tilde{x}\in(\min\{x^{\prime},x\},\max\{x^{\prime},x\})$. Plugging the above expansion in the model (2.6) we get

$\begin{split}&\dfrac{d}{dt}\int_{\mathcal{I}}\varphi(x)f(x,t)dx=\\&\qquad\epsilon\int_{\mathcal{I}\times\mathcal{I}}B_{\epsilon}[f](x,t)\varphi^%{\prime}(x)P(x)(z-x)f(x,t)g(z,t)dz\,dx+\\&\qquad\dfrac{\epsilon}{2\lambda}\int_{\mathcal{I}\times\mathcal{I}}B_{%\epsilon}[f](x,t)\varphi^{\prime\prime}(x)D^{2}[f](x,t)f(x,t)g(z)dz\,dx+\\&\qquad R(f)(x,t),\end{split}$ |

where $R(\cdot)$ is the so-called reminder term

$\begin{split}&R(f)(x,t)\\&=\dfrac{1}{2}\int_{\mathcal{I}\times\mathcal{I}}B_{\epsilon}[f](x,t)\dfrac{d^%{2}}{dx^{2}}\varphi(x)\epsilon^{2}P^{2}(x)(z-x)^{2}f(x,t)dx\\&+\dfrac{1}{6}\int_{\mathcal{I}\times\mathcal{I}}B_{\epsilon}[f](x,t)\left%\langle(x^{\prime}-x)^{3}\right\rangle f(x,t)g(z)dx\,dz,\end{split}$ | (3.9) |

with

$\begin{split}&\left\langle(x^{\prime}-x)^{3}\right\rangle=\epsilon^{3}P^{3}(x)%(z-x)^{3}\\&\qquad+\dfrac{3\epsilon^{2}}{\lambda}P(x)(z-x)D^{2}[f](x,t)+D^{3}[f](x,t)%\left\langle\xi^{3}\right\rangle\end{split}$ |

The modeling justification of the above expansion is now clear: we look for a regime where the single interactions between an individual and the background society are able only to slightly influence the personal believes. We look indeed to this process through a Taylor expansion and so to situations where $x^{\prime}$ is close to $x$ as a result of a single exchange of ideas. Hence, by introducing the time scaling $\tau=t\epsilon$ and the scaled distribution $f_{\epsilon}(x,\tau)=f(x,t/\epsilon)$, we look now to the long time behavior of this process. This gives

$\begin{split}&\dfrac{d}{d\tau}\int_{\mathcal{I}}\varphi(x)f_{\epsilon}(x,\tau)%dx\\&=\int_{\mathcal{I}\times\mathcal{I}}B_{\epsilon}[f_{\epsilon}](x,\tau)\varphi%^{\prime}(x)P(x)(z-x)f_{\epsilon}(x,\tau)g(z,t)dz\,dx\\&+\dfrac{1}{2\lambda}\int_{\mathcal{I}\times\mathcal{I}}B_{\epsilon}[f_{%\epsilon}](x,\tau)\varphi^{\prime\prime}(x)D^{2}[f_{\epsilon}](x,t)f_{\epsilon%}(x,t)g(z)dz\,dx\\&+\dfrac{1}{\epsilon}R(f_{\epsilon})(x,t).\end{split}$ | (3.10) |

Moreover, from (3.9) we also get for $\epsilon\to 0^{+}$

$|R(f_{\epsilon})|\lesssim(\epsilon^{2}+\epsilon^{3})\min\{f_{\epsilon},%\epsilon^{-\alpha}\}^{\alpha}\to 0.$ |

Therefore, from (3.10), the limit $\epsilon\to 0^{+}$ gives

$\begin{split}&\dfrac{d}{d\tau}\int_{\mathcal{I}}\varphi(x)f_{0}(x,\tau)dx\\&\qquad=\int_{\mathcal{I}\times\mathcal{I}}B_{0}[f_{0}](x,\tau)\varphi^{\prime%}(x)P(x)(z-x)f_{0}(x,\tau)g(z)dz\,dx\\&\qquad+\dfrac{1}{2\lambda}\int_{\mathcal{I}\times\mathcal{I}}B_{0}[f_{0}](x,%\tau)\varphi^{\prime\prime}(x)D^{2}[f_{0}](x,t)f_{0}(x,t)g(z)dz\,dx.\end{split}$ |

We can now integrating back by parts obtaining

$\begin{split}&\partial_{\tau}f_{0}(x,\tau)\\&=\partial_{x}\left[P(x)(x-m)f_{0}(x,\tau)B_{0}[f_{0}](x,t)+\dfrac{1}{2\lambda%}\partial_{x}\left(H(x)f_{0}(x,\tau)\right)\right]\end{split}$ |

thanks to an appropriate choice of the boundary conditions

$\begin{split}P(x)(x-m)f_{0}(x,\tau)B_{0}[f_{0}](x,t)+\dfrac{1}{2\lambda}%\partial_{x}\left(H(x)f_{0}(x,\tau)\right)\Big{|}_{x=\pm 1}=0\\H(x)f_{0}(x,\tau)\Big{|}_{x=\pm 1}=0\end{split}$ |

affirming that there is no mass exiting from the boundaries of the domain, i.e. opinions can be at maximum extremely negative ($x=-1$) or extremely positive ($x=1$).

Coming back now for simplicity to the original variables, we observe that starting from (2.6) we reduced ourselves to study the evolution of the distribution solution of a new nonlinear Fokker-Planck model reading

$\begin{split}&\partial_{t}f(x,t)=\partial_{x}\Big{[}P(x)(x-m)(1+\beta(H(x)f(x,%t))^{\alpha})f(x,t)\\&+\dfrac{1}{2\lambda}\partial_{x}(H(x)f(x,t))\Big{]}.\end{split}$ | (3.11) |

Interesting enough, the model obtained from this expansion shares many similarities with the well studied Kaniadakis-Quarati equation [37, 53] with the main difference that in our case the solution lives in the bounded interval $\mathcal{I}$. The original Kaniadakis-Quarati model was indeed a Fokker-Planck equation with quadratic drift describing the dynamics of bosons exhibiting, in the spatially homogeneous setting, condensation behaviors. We may then expect that our model to be eventually able to generate in the long time similar behaviors.

We now notice that, under the assumption $P\equiv 1$ and $H(x)=1-x^{2}$, we get

$\begin{split}&\partial_{t}f(x,t)=\partial_{x}\Big{[}(x-m)(1+\beta((1-x^{2})f(x%,t))^{\alpha})f(x,t)+\\&+\dfrac{1}{2\lambda}\partial_{x}((1-x^{2})f(x,t))\Big{]}.\end{split}$ | (3.12) |

This Fokker-Planck model, in the absence of conformity forces, $\beta\equiv 0$, reduces to the usual Fokker-Planck model for opinion consensus in a bounded interval introduced in [52] and for which the large time distribution reads

$f_{\infty}(x)=c_{m,\lambda}\left(1+x\right)^{-1+\lambda(1+m)}\left(1-x\right)^%{-1+\lambda(1-m)},$ | (3.13) |

with $c_{m,\lambda}>0$ a normalization constant. This equilibrium distribution corresponds to a beta-type distribution in the interval $\mathcal{I}$. Another possible choice, enforcing faster decay at the boundaries is $H(x)=(1-x^{2})^{2}$, which, under the assumption $P\equiv 1$ gives

$\begin{split}\partial_{t}f(x,t)=&\partial_{x}\Big{[}(x-m)(1+\beta((1-x^{2})^{2%}f(x,t))^{\alpha})f(x,t)+\\&+\dfrac{1}{2\lambda}\partial_{x}((1-x^{2})^{2}f(x,t))\Big{]}.\end{split}$ | (3.14) |

Therefore, in the absence of conformity forces, $\beta=0$, (3.14) triggers the emergence of the equilibrium distribution

$\begin{split}&f_{\infty}(x)=c_{m,\lambda}(1+x)^{-2+\lambda m/2}(1+x)^{-2-%\lambda m/2}\exp\left\{-\dfrac{\lambda(1-mx)}{1-x^{2}}\right\},\end{split}$ | (3.15) |

where, as before, $c_{m,\lambda}>0$ is a normalization constant. This second case depicts with respect to the case $H(x)=1-x^{2}$, a dynamics where strong positive or negative opinions can be handled with more difficulties by individuals in a population.

In contrast to the smooth cases (3.13) and (3.15) characterizing dynamics of opinion formation in absence of mechanisms related to conformity, in this work, we are interested to the case of emergence of consensus when group pressures are active and may possibly modify the structure of the solution or in other words the ideas of individuals. For this reason, in the next section we will explore from the analytical point of view the specific cases $\alpha,\beta>0$ where the conformity forces are active leading to the new nonlinear Fokker-Planck model for opinion formation introduced in (3.11). This study will highlight the different nature of our new approach introducing a new way to describe consensus phenomena not necessarily based on bounded confidence interactions [34].

### 3.1 Blow up, existence of finite critical masses, steady states and entropy

First of all, we want to prove the existence of a finite critical mass for the nonlinear Fokker-Planck model (3.11). This quantity identifies the passage from a smooth steady state to a steady state presenting blow up. From a modeling point of view, the existence of a critical mass means that the size of the population plays a critical role in the consensus process. When, indeed the number of persons sharing a similar idea overtakes a threshold, consensus becomes unavoidable due to social pressure. We successively want to show that the steady state for this model corresponds to the unique minimizer of a suitable entropy functional. From the technical point of view, we follow in the sequel the approach proposed in [10] in the case of condensates of bosons.

Let observe now that equation (3.11) belongs to a more general class of Fokker-Planck equations having the form

$\partial_{t}f(x,t)=\partial_{x}\left[(x-m)\mathcal{D}[f]+\sigma^{2}\partial_{x%}\mathcal{R}[f]\right],$ | (3.16) |

being $\mathcal{D}[f]>0$ the drift functional and $\mathcal{R}[f]$ a diffusion functional. This general form will be then used in the following to characterize the solution. We start now by rewriting equation (3.11) for $x\in(-1,1)$ with respect to $g(x,t)=H(x)f(x,t)$ to get

$\begin{split}&\dfrac{1}{H(x)}\partial_{t}g(x,t)\\&\quad=\frac{\partial}{\partial x}\left[\frac{x-m}{H(x)}g(x,t)(1+\beta g^{%\alpha}(x,t))+\dfrac{1}{2\lambda}\partial_{x}g(x,t)\right],\end{split}$ | (3.17) |

from which we get that the unique steady state $g_{\infty}(x)$, if it exists, is such that

$2\lambda\frac{x-m}{H(x)}g_{\infty}(x)(1+\beta g_{\infty}^{\alpha}(x))+\partial%_{x}g_{\infty}(x)=0,$ |

which gives

$\partial_{x}\log\dfrac{g_{\infty}(x)}{(1+\beta g_{\infty}^{\alpha}(x))^{1/%\alpha}}=-2\lambda\dfrac{x-m}{H(x)}.$ | (3.18) |

It is therefore of interest to understand the condition under which $g^{\infty}(x)$ solution to (3.18) is a density function, i.e. $g^{\infty}(x)\geq 0$, additionally weighted by a given mass $M$ defined later. To this end, we observe that (3.11) recast in the form (3.16) gives as drift functional

$\mathcal{D}[g]=g(1+\beta g^{\alpha}),$ |

where $\mathcal{D}[g]>0$ for all $g>0$. We observe also that $\mathcal{D}[g]$ is superlinear, meaning that

$\int_{1}^{\infty}\frac{1}{\mathcal{D}(\rho)}{d}\rho\leq K<\infty.$ | (3.19) |

The steady state of (3.17) of given mass $\mu=\int_{\mathcal{I}}f(x,t)dx>0$ solves then

$\dfrac{1}{2\lambda}\partial_{x}g+\dfrac{x-m}{H(x)}\mathcal{D}[g]=0,$ |

which can be rewritten as

$\dfrac{1}{\mathcal{D}[g]}\partial_{x}g+2\lambda\partial_{x}V(x)=0,$ | (3.20) |

supposing that a $V(x)\geq 0$ such that

$\dfrac{d}{dx}V(x)=\dfrac{x-m}{H(x)},$ |

exists. To continue the computation we need now to specify the shape of the function $V(x)$. To that aim, the two leading examples for the shape of the concave function $H(x)$, which scope is to make the diffusion disappear close to the boundaries, are analyzed. In the case $H(x)=1-x^{2}$, we get

$V(x)=\log\left[(1-x)^{\frac{m-1}{2}}(1+x)^{-\frac{m+1}{2}}\right],$ |

whereas for $H(x)=(1-x^{2})^{2}$ we get

$V(x)=\dfrac{1}{4}\left(\dfrac{2(mx-1)}{x^{2}-1}+\log\left(\dfrac{1-x}{1+x}%\right)^{m}\right).$ |

We define now $\Phi^{\prime}(g)$ as

$\Phi^{\prime}(g)=-\int_{g}^{\infty}\frac{1}{\mathcal{D}(\rho)}d\rho,$ | (3.21) |

which gives, in the case $\mathcal{D}[g]=g(1+\beta g^{\alpha})$

$\Phi^{\prime}(g)=\dfrac{1}{\alpha}\log\dfrac{g^{\alpha}}{1+\beta g^{\alpha}}.$ | (3.22) |

Substituting now $\Phi(\cdot)$ into (3.20), we obtain that the steady state is solution to

$\partial_{x}\left[\Phi^{\prime}(g)+2\lambda V(x)\right]=0,$ | (3.23) |

or equivalently

$\Phi^{\prime}(g)=-2\lambda V(x)-c,$ | (3.24) |

with $c$ a constant such that $c>0$. Since $\Phi^{\prime}(\cdot)$ is strictly increasing, denoting by $[\Phi^{\prime}]^{-1}$ its inverse, we can alternatively express the steady state of (3.17) as follows

$g_{\infty,c}=[\Phi^{\prime}]^{-1}\left(-2\lambda V(x)-c\right).$ | (3.25) |

As a consequence of (3.21), we have that

$\Phi(g)=\int_{0}^{g}\Phi^{\prime}(\rho)d\rho,$ |

where we observe that the function $\Phi(\rho)$, which has a derivative that is non-decreasing, is then convex for $\rho\geq 0$. This latter observation permits therefore to consider an entropy functional associated to the steady state solution(3.23)

$\mathcal{H}[g](t)=\int_{\mathcal{I}}\left(2\lambda V(x)g+\Phi(g)\right)\dfrac{%1}{H(x)}dx,$ |

which has the property to decrease in time along the solution to equation (3.17). In fact we have that

$\begin{split}\frac{d}{dt}\mathcal{H}[g](t)=&\int_{-1}^{1}\left(2\lambda V(x)+%\Phi^{\prime}(g)\right)\dfrac{1}{H(x)}\partial_{t}g\;dx\\=&\int_{-1}^{1}\left(2\lambda V(x)+\Phi^{\prime}(g)\right)\frac{\partial}{%\partial x}\left\{\mathcal{D}[g]\frac{\partial}{\partial x}\left[2\lambda V(x)%+\Phi^{\prime}(g)\right]\right\}\;dx\\=&-\int_{-1}^{1}\mathcal{D}[g]\left(\frac{\partial}{\partial x}\left[2\lambda V%(x)+\Phi^{\prime}(g)\right]\right)^{2}\;dx\leq 0.\end{split}$ |

Thus, we have obtained that, starting from an initial density with a given mass $\mu$, the solution is such that $\mathcal{H}[\cdot]$ is dissipated monotonically. Formally, we have obtained that the steady state is the one given by (3.25) which, by construction, preserves the initial mass $\mu$ over time.

We can distinguish now two different situations depending on the initial value of $\mu$. The first case refers to a smooth equilibrium distribution, while the second case refers to a solution presenting a blow-up where concentration of a part of the total mass around the mean opinion $m=\int_{\mathcal{I}}xf(x,t)dx$ can be observed. To that aim, we focus now on the specific case $m=0$. The steady state solution reads in this situation

$g_{\infty,c}(x)=\dfrac{e^{-2\lambda V(x)-c}}{\left(1-\beta e^{-2\lambda V(x)-c%}\right)^{1/\alpha}},$ |

which may eventually have a critical state and thus presents a blow-up, only if $1-\beta e^{-c}e^{-2\lambda V(x)}=0$. In the case $H(x)=1-x^{2}$, this implies that the constant in equation (3.24) has to take the value $c=\log\beta$, being also $x=0$ the unique minimum of $V(x)\geq 0$ with $V(0)=0$. On the other hand, in the case $H(x)=(1-x^{2})^{2}$ and $m=0$, with similar computations one gets $c=\log\beta-\lambda$. We want now to find the values of the parameters in (3.25) such that the integral of $g_{\infty,c}$ is finite, in other words the solution found has a physical sense.To that aim, let us go back to (3.25) with $H(x)=1-x^{2}$ and $m=0$. In this situation, we have that the critical mass $\mu_{c}(g_{\infty})$ is

$\begin{split}\mu(g_{\infty})&=\int_{-1}^{1}g_{\infty,\log\beta}\;{d}x\\&=\int_{-1}^{1}\left[\Phi^{\prime}\right]^{-1}\left(\lambda\log\left(1-x^{2}%\right)-\log\beta\right)dx\\&=2\int_{0}^{1}\left[\Phi^{\prime}\right]^{-1}\left(\lambda\log(1-x^{2})-\log%\beta\right)dx.\end{split}$ |

Following [10] we change the variables in the above integral with $s$ so that $\Phi^{\prime}(s)=\lambda\log(1-x^{2})-\log\beta$ and we get

$\begin{split}\mu(g_{\infty})=&-\dfrac{1}{\lambda}\int_{+\infty}^{0}\dfrac{s%\left(\beta e^{\Phi^{\prime}(s)}\right)^{1/\lambda}\Phi^{\prime\prime}(s)}{%\sqrt{1-\left(\beta e^{\Phi^{\prime}(s)}\right)^{1/\lambda}}}ds\\=&\dfrac{1}{\lambda}\int_{0}^{+\infty}\dfrac{\beta^{1/\lambda}\left(\dfrac{s^{%\alpha}}{1+\beta s^{\alpha}}\right)^{\frac{1}{\lambda\alpha}}}{(1+\beta s^{%\alpha})\sqrt{1-\beta^{1/\lambda}\left(\dfrac{s^{\alpha}}{1+\beta s^{\alpha}}%\right)^{\frac{1}{\lambda\alpha}}}}ds.\end{split}$ | (3.26) |

The integrand in the last integral appearing in (3.26) as $s\to\infty$ behaves like $s^{-\alpha/2}$, so there exists a finite critical mass only if $\alpha>2$.

A similar computation can be done in the case $H(x)=(1-x^{2})^{2}$ and $m=0$. In this case, the critical mass $\mu(g_{\infty})$ related to the drift $\mathcal{D}(g)=g(1+\beta g^{\alpha})$ is obtained as

$\begin{split}\mu(g_{\infty})&=\int_{-1}^{1}g_{\infty,\log\beta-\lambda}dx=\int%_{-1}^{1}[\Phi^{\prime}]^{-1}\left(-\dfrac{\lambda}{x^{2}-1}-\log\beta+\lambda%\right)dx\\&=2\int_{0}^{1}[\Phi^{\prime}]^{-1}\left(-\dfrac{\lambda}{x^{2}-1}-\log\beta+%\lambda\right)dx.\end{split}$ |

Hence, we perform the change of variables from $\Phi^{\prime}(s)=-\dfrac{\lambda}{x^{2}-1}-\log\beta+\lambda$ and we get

$\begin{split}\mu(g_{\infty})&=\int_{0}^{+\infty}\dfrac{\lambda\Phi^{\prime%\prime}(s)s}{\left(\Phi^{\prime}(s)+\log\beta-\lambda\right)^{3/2}}\dfrac{1}{%\sqrt{\Phi^{\prime}(s)+\log\beta}}ds\\&=\int_{0}^{+\infty}\dfrac{\lambda\left(\log\dfrac{\beta s}{(1+\beta s^{\alpha%})^{1/\alpha}}\right)^{-1/2}}{(1+\beta s^{\alpha})\left(\log\left(\dfrac{s^{%\alpha}}{1+\beta s^{\alpha}}\right)^{1/\alpha}+\log\beta-\lambda\right)^{3/2}}%ds,\end{split}$ | (3.27) |

where the integrand behaves as $s^{-\alpha/2}$ for $s\to+\infty$ and a critical mass exists, as before, only if $\alpha>2$.To conclude, as stated in the introduction, our new model in the case $m=0$ and for values of the function $H(x)=1-x^{2}$ or $H(x)=(1-x^{2})^{2}$ presents blow-up when the total density of individuals exceeds a fixed critical level and $\alpha>2$. This result will be used in the numerical section to describe strong polarization phenomena due to the presence of social pressure.

### 3.2 Critical value for non centered steady states

In the case $m\neq 0$, the computations done in the previous section which lead to the existence of a finite critical mass cannot be explicitly done. Thus, in what follows we make use of a different argument to prove an analogous result. In particular we start from the explicit expression of the steady state obtained from (3.17).

Let us first consider the case $H(x)=1-x^{2}$ . To this extent, consider that (3.20) implies that the distribution $g$ of the stationary state satisfies the differential equation

$\frac{1}{g(1+\beta g^{\alpha})}\partial_{x}g=-2\lambda\frac{x-m}{1-x^{2}}.$ | (3.28) |

Since

$\frac{1}{g(1+\beta g^{\alpha})}=\frac{1}{g}-\frac{1}{\alpha}\frac{\beta\alpha g%^{\alpha-1}}{1+\beta g^{\alpha}},$ |

and

$\frac{1}{g(1+\beta g^{\alpha})}\partial_{x}g=\partial_{x}\log\left(\frac{g}{(1%+\beta g^{\alpha})^{1/\alpha}}\right),$ |

we get

$\partial_{x}\log\left(\frac{g^{\alpha}}{1+\beta g^{\alpha}}\right)=-2\alpha%\lambda\partial_{x}\log\left[(1+x)^{-\frac{m+1}{2}}(1-x)^{\frac{m-1}{2}}\right].$ |

Therefore

$\frac{g^{\alpha}}{(1+\beta g^{\alpha})}=C^{\alpha}(1+x)^{\lambda\alpha(1+m)}(1%-x)^{\lambda\alpha(1-m)},$ | (3.29) |

that implies

$g_{\infty,C}=\frac{C(1+x)^{\lambda(1+m)}(1-x)^{\lambda(1-m)}}{(1-\beta C^{%\alpha}(1+x)^{\lambda\alpha(1+m)}(1-x)^{\lambda\alpha(1-m)})^{1/\alpha}},$ |

and, finally

$f_{\infty,C}=\frac{C(1+x)^{\lambda(1+m)-1}(1-x)^{\lambda(1-m)-1}}{(1-\beta C^{%\alpha}(1+x)^{\lambda\alpha(1+m)}(1-x)^{\lambda\alpha(1-m)})^{1/\alpha}}.$ | (3.30) |

Now, it is important to remark that the recovered steady state must be a non negative function. This means that, for a given positive constant $C$ the denominator of the fraction in (3.30) must be non negative. This implies an upper bound on the possible values taken by the constant $C$. In this regard it is sufficient to note that for $g\geq 0$, $g^{\alpha}/(1+\beta g^{\alpha})\leq 1/\beta$. Consequently the constant $C$ appearing in (3.29) has to satisfy the condition

$C(1+x)^{\lambda(1+m)}(1-x)^{\lambda(1-m)}\leq\left(\frac{1}{\beta}\right)^{{1}%/{\alpha}}.$ |

On the other hand, the function $(1+x)^{r}(1-x)^{s}$, for $\left|x\right|\leq 1$, takes the maximum in $\overline{x}=(r-s)/(r+s)$, so that

$(1+x)^{r}(1-x)^{s}\leq\left(\frac{2r}{r+s}\right)^{r}\left(\frac{2s}{r+s}%\right)^{s}.$ |

In our case, since $r=\lambda(1+m)$, $s=\lambda(1-m)$, the maximum value is attained at $\overline{x}=m$ which implies the condition

$C(1+m)^{\lambda(1+m)}(1-m)^{\lambda(1-m)}\leq\left(\frac{1}{\beta}\right)^{%\frac{1}{\alpha}}.$ |

Finally, the constant $C$ in (3.30) can not be greater than $C_{M}$, given by

$C_{M}=\left(\frac{1}{\beta}\right)^{\frac{1}{\alpha}}(1+m)^{-\lambda(1+m)}(1-m%)^{-\lambda(1-m)}.$ |

In correspondence to this value, the steady state $f_{\infty,{C_{M}}}(x)$ is such that

$\lim_{x\to m}f_{\infty,{C_{M}}}(x)=+\infty,$ |

i.e. we observe a blow-up. Moreover, a direct computation shows that

$\lim_{x\to m}\frac{1-\beta C_{M}^{\alpha}(1+x)^{\lambda\alpha(1+m)}(1-x)^{%\lambda\alpha(1-m)}}{(x-m)^{2}}=\lambda\alpha(1-m^{2})^{-1},$ |

which implies that the steady state $f_{\infty,{C_{M}}}$ tends to infinity at the order $2/\alpha$ as $x\to m$ as for the case $m=0$. Consequently,the integral

$\mu_{c}=\int_{-1}^{1}f_{\infty,{C_{M}}}(x)\mathrm{d}x$ | (3.31) |

is bounded as soon as $\alpha>2$. Therefore, as $\alpha>2$ we conclude with the existence of a finite critical mass $\mu_{c}$, coherently with what proven in (3.26).Next, we consider the case $H(x)=(1-x^{2})^{2}$. From (3.20) we get that the distribution $g$ of the stationary state is such that

$\dfrac{\partial_{x}g}{g(1+\beta g^{\alpha})}=-2\lambda\dfrac{x-m}{(1-x^{2})^{2%}}.$ |

Hence, proceeding as before, we have

$\partial_{x}\log\left(\dfrac{g^{\alpha}}{1+\beta g^{\alpha}}\right)=-\alpha%\lambda\partial_{x}\left[\dfrac{1-mx}{1-x^{2}}+\dfrac{1}{2}\log\left(\dfrac{1-%x}{1+x}\right)^{m}\right],$ |

and therefore

$\dfrac{g^{\alpha}}{1+\beta g^{\alpha}}=C^{\alpha}e^{-\alpha\lambda\frac{1-mx}{%1-x^{2}}}\left(\dfrac{1-x}{1+x}\right)^{-\frac{\alpha\lambda m}{2}},$ |

which implies

$g_{\infty,C}=\dfrac{Ce^{-\lambda\frac{1-mx}{1-x^{2}}}\left(\frac{1-x}{1+x}%\right)^{-\frac{\lambda m}{2}}}{1-\beta C^{\alpha}e^{-\alpha\lambda\frac{1-mx}%{1-x^{2}}}\left(\frac{1-x}{1+x}\right)^{-\frac{\alpha\lambda m}{2}}},$ | (3.32) |

and, finally

$f_{\infty,C}=(1-x^{2})^{-2}\dfrac{Ce^{-\lambda\frac{1-mx}{1-x^{2}}}\left(\frac%{1-x}{1+x}\right)^{-\frac{\lambda m}{2}}}{\left(1-\beta C^{\alpha}e^{-\alpha%\lambda\frac{1-mx}{1-x^{2}}}\left(\frac{1-x}{1+x}\right)^{-\frac{\alpha\lambdam%}{2}}\right)^{1/\alpha}}.$ | (3.33) |

To determine the critical value of the constant $C$ in (3.32) we get as before the condition

$Ce^{-\lambda\frac{1-mx}{1-x^{2}}}\left(\dfrac{1-x}{1+x}\right)^{-\frac{\lambdam%}{2}}\leq\left(\frac{1}{\beta}\right)^{1/\alpha},$ |

whose maximum is attained at $x=m$, which implies the condition

$Ce^{-\lambda}\left(\dfrac{1-m}{1+m}\right)^{-\frac{\lambda m}{2}}\leq\left(%\frac{1}{\beta}\right)^{1/\alpha}.$ |

Therefore, the constant $C>0$ cannot be greater than $C_{M}$ given by

$C_{M}=\left(\frac{1}{\beta}\right)^{1/\alpha}e^{\lambda}\left(\dfrac{1-m}{1+m}%\right)^{\frac{\lambda m}{2}},$ |

for which the steady state $f_{\infty,C_{M}}(x)$ is such that

$\lim_{x\to m}f_{\infty,C_{M}}(x)=+\infty.$ |

Furthermore, we get

$\lim_{x\to m}\dfrac{1-\beta C_{M}^{\alpha}e^{-\alpha\lambda\frac{1-mx}{1-x^{2}%}}\left(\frac{1-x}{1+x}\right)^{-\frac{\alpha\lambda m}{2}}}{(x-m)^{2}}=\alpha%\lambda(1-m^{2})^{-2},$ |

implying that also the steady state $f_{\infty,C_{M}}$ defined in (3.33) tends to infinity at the order $2/\alpha$ as soon as $x\to m$ and the integral

$\mu_{c}=\int_{-1}^{1}f_{\infty,C_{M}}dx$ |

is bounded for $\alpha>2$. Therefore, also in the case $H(x)=(1-x^{2})^{2}$ we obtain the existence of the critical mass $\mu_{c}$ coherently with what proven in (3.27) in the case $m=0$.

## 4 Numerical results

In this section, we explore the capability of the proposed model (2.6) and of its grazing collision limit (3.11) to describe certain features related to opinion formation. In particular, as already stated, we focus on the study of conformism within a population and how this mechanism plays a role in shaping ideas in a community. We distinguish in the sequel the case in which opinion is formed and regularly distributed among its possible spectrum, i.e. $x\in[-1,1]$, and we refer to this situation to as the sub-critical case and the condition in which a predominant idea possessed by a given size of population forces other individuals to align in order to conform and avoid strong criticisms from the leading part of the society. We refer to the latter situation as the super-critical case. To highlight the link between the introduced microscopic dynamics and the emerging equilibrium distribution we stick to a Monte Carlo approach [46]. For the seek of completeness, we mention also the deterministic numerical methods for similar equations developed in [9, 17, 41] and the references therein.In the first part, we validate the proposed theoretical approach through a body of numerical results. In the second part, we consider several tests highlighting the possibility to reach dynamically a supercritical state in a scenario involving more than one population.

### 4.1 Validating the theoretical setting

In this first set of simulations, we show how the Monte Carlo method applied to the Boltzmann model (2.6) behaves for different values of the modeling and scaling parameters. In particular, we are interested in numerically showing that our method is able to describe the steady state solutions both in the smooth as well as in the blow-up situation and that this steady state corresponds to the one obtained in the previous section for the Fokker-Planck limit. In the following, we first compute the exact stationary solution for various sets of parameters $\alpha,\lambda$ and $m$ and then we compare it with MC. We start from a uniform distribution in $[-1,1]$ and we use $N_{s}=10^{5}$ samples with $H(x)=1-x^{2}$ and $\epsilon=10^{-3}$. For the simulations shown in Figure 4.1 we choose $\alpha=3$, $\lambda=4$ and $m=0$. The critical mass $\mu_{c}$ for this choice of parameters reads as in (3.31). Given this value, we compute the stationary states for respectively $\mu=\mu_{c}/2$ and $\mu=\mu_{c}+0.5$. In the first case, as expected, the stationary solution is smooth since its mass is below the critical one, while in the second case instead we can clearly see the formation of a singularity in $m=0$, which is very well caught by our Monte Carlo scheme.

To underline the importance of the relation existing between the parameters and the possible blow-up around the mean position $m$, we perform another numerical test, now fixing the mass $\mu=1.5$ and showing the different stationary states as we change the values of $\alpha$ and $\lambda$ still keeping $H(x)=1-x^{2}$. The initial condition is as before the uniform distribution in $[-1,1]$ and, again, we used $N_{s}=10^{5}$ particles with scaling parameter $\epsilon=10^{-3}$. We report the results in Figure 4.2. On the left, we can see the stationary state for the choice $\alpha=2.9$, $\lambda=7$, and $m=0.1$, where the fixed mass $\mu=1$ is below the critical one, $\mu_{c}=2.112$, so that the stationary solution is smooth. On the contrary, on the right, we chose $\alpha=4$, $\lambda=11$, and $m=0.1$ and for such parameters we observe the super-critical situation: the mass starts to accumulate around $m$. In Figure 4.3 we show the case of convergence to equilibrium in the sub and super critical cases where now $H(x)=(1-x^{2})^{2}$. Again, the Monte Carlo method is able to capture the smooth and non smooth equilibrium states with high accuracy.

In order to show the sensitivity of the Monte Carlo simulations to the mesh size used to reconstruct the density $f$ that is needed to compute the collision frequency and the diffusion term as detailed in algorithm 1, we perform several simulations with different grid points. To that aim, in Figure 4.4, there are shown these results obtained when trying to catch the stationary states with different mesh sizes indexed by $n$, the number of bins necessary to reconstruct $f$, as $\alpha=3$ and $\lambda=4$ are fixed, both for the sub-critical case on the left as well as for the super-critical case on the right. We can see how lower values of $n$ may lead to instabilities near the mean value $m$. At the same time, we observe that even higher values of mesh points may not be sufficient to obtain a good representation of the steady states and that the numerical method can fail in representing the solution.

Figure 4.5 reports instead the numerical convergence of the Boltzmann model (2.6) to the Fokker-Planck one (3.11) for different values of the scaling parameter $\epsilon$. We clearly see that, as soon as $\epsilon$ becomes sufficiently small, the Boltzmann model collapses towards the Fokker-Planck solution both for $H(x)=1-x^{2}$ as well as for $H(x)=(1-x^{2})^{2}$.On the other hand, when $\epsilon$ is large enough the two dynamics lead to quite different results especially in the smooth case. Figure 4.6 shows for completeness the time evolution in the sub and super critical cases of the Boltzmann model in the Fokker-Planck setting starting from an uniform distribution for a mean opinion $m=0.3$ and $H(x)=1-x^{2}$.

Finally we consider, to document the qualitative behavior of the model, a situation in which the mean opinion $m$ is not fixed in time while instead it is driven by some external forces. In Figure 4.7 we show on the top the case $m(t)=0.4\sin(0.5t)$, while on the bottom we set $m(t)=0.5\mbox{tanh}(t-50)$. In both simulations we set $\alpha=4$, $\lambda=5$, $\mu=1.5$ and $H(x)=(1-x^{2})$.

### 4.2 Interacting populations

We now concentrate on the time evolution of opinions between two competing populations of individuals. We suppose that personal ideas may change only as a result of interaction with people belonging to the same group. However, ideas are also shared with all individuals, independently on the group to which the agents belong, and, consequently, one individual may decide to move from one group to the other one when interacting with individuals of the other group with similar ideas. These dynamics are constructed in such a way to mimic the competition between political parties where one party tries to increase its supporters by moving closer to the opinion of individuals belonging to another party in order to capture their interest and tempt them to change faction.We suppose then that these two populations are characterized by the time evolution of their internal opinions. The process of opinion formation is regulated by equation (3.11), where the variable $x\in[-1,1]$ represents the opinion on a certain topic, with $x=-1$ corresponding to a completely negative opinion on the subject, while $x=1$, on the other hand, signifies a totally positive opinion on the same topic. The function $f_{1}(t,x)$ now represents the density of opinions at time $t>0$ of the first population while $f_{2}(t,x)$ of the second one. We suppose that, at time $t=0$, their densities $f_{1,0}$ and $f_{2,0}$ are given by

$f_{1}(x,0)=\begin{cases}C_{1}\exp\left\{-\frac{1}{2}\frac{(x+0.1)^{2}}{0.05^{2%}}\right\}&x\in[-1,1]\\0&\textrm{elsewhere},\end{cases}$ |

and

$f_{2}(x,0)=\begin{cases}C_{2}\exp\left\{-\frac{1}{2}\frac{(x-0.2)^{2}}{0.1^{2}%}\right\}&x\in[-1,1]\\0&\textrm{elsewhere},\end{cases}$ |

with $C_{1},\,C_{2}$ constants such that

$\int_{\mathcal{I}}f_{1}(x,0)\mathrm{d}x=\rho_{1}(t=0)\quad\mbox{and}\quad\int_%{\mathcal{I}}f_{2}(x,0)\mathrm{d}x=\rho_{2}(t=0),$ |

where $\rho_{i}(t)$, $i=1,2$, indicates the total mass of population $i$ at time $t$ and for each $t\geq 0$, $\rho_{1}(t)+\rho_{2}(t)=m_{tot}$, i.e. the total number of people, remains constant over time.

The competing populations model is then driven by the following system of equation

$\begin{cases}\partial_{t}f_{1}(t,x)=-K_{1}(f_{1},f_{2})+K_{2}(f_{2},f_{1})+%\frac{1}{\tau}Q^{1}_{BE}(f_{1},f_{1}),\\\partial_{t}f_{2}(t,x)=K_{1}(f_{1},f_{2})-K_{2}(f_{2},f_{1})+\frac{1}{\tau}Q^{%2}_{BE}(f_{2},f_{2}),\end{cases}$ | (4.34) |

where the operators $K_{i}$, $i=1,2$, regulate the exchange between the two groups and $Q^{i}_{BE}$ is the Fokker-Planck operator modeling the shape of political positions. The interaction between the two groups is defined as

$\displaystyle K_{1}(f_{1},f_{2})=$ | $\displaystyle\frac{f_{1}(t,x)}{\rho_{1}(t)m_{tot}}\int_{-1}^{1}\kappa_{1}(x,x_%{*})f_{2}(t,x_{*})\mathrm{d}x_{*},$ | (4.35) | ||

$\displaystyle K_{2}(f_{2},f_{1})=$ | $\displaystyle\frac{f_{2}(t,x)}{\rho_{2}(t)m_{tot}}\int_{-1}^{1}\kappa_{2}(x,x_%{*})f_{1}(t,x_{*})\mathrm{d}x_{*},$ |

with

$\kappa_{1,2}(x,x_{*})=\chi\left(\left|x-x_{*}\right|\leq\Delta_{1,2}\right)%\gamma_{1,2}.$ | (4.36) |

The constants $0\leq\gamma_{1,2}\leq 1$ represent the transition rates between populations and the presence of the indicator function $\chi\left(\left|x-x_{*}\right|\leq\Delta_{1,2}\right)$ models the fact that the transition probability is str