Laplacian dynamics of convergent and divergent collective behaviors

Collective behaviors refer to the activities exhibited by a group of interacting units (e.g. animals or people) [1]. Numerous collective behaviors dynamically evolve with time or other factors [2–8]. The convergent and divergent evolution processes of these behaviors exist across different scales. For example, in similar environments, multiple species that interact with each other may evolve along converging pathways [9–11], while the closely related populations of a species may evolve along diverging pathways if they experience distinct selective pressures [12–14]. As another example, people may exhibit convergent social behaviors (e.g. opinions or consuming behaviors) if emergencies occur (e.g. epidemic or social turmoil), while these social behaviors diverge and become increasingly individualized after emergencies end [8, 15, 16].

Although the universality of convergent and divergent collective behaviors ensures its ubiquity in the world, it does propose critical challenges for analytic modeling due to the rich connotations and diverse details of collective behaviors. It is unclear whether the convergent and divergent dynamics of organic evolution, the convergent and divergent social behaviors of people, and any other similar phenomenon can be characterized by a unified and simple model. While specialized models [1, 2, 17–24], statistical analyses [25–29], and data-driven methods [30–33] have achieved essential progress in studying collective dynamics, advances in describing the convergent and divergent evolution processes of collective behaviors remain limited.

In this study, we propose a simple model of convergent and divergent collective behaviors and describe the dynamic switching of convergent and divergent evolution across time. This framework may serve as a minimal model to characterize diverse real collective dynamics. In addition, our model can be enriched by considering more agent heterogeneity and local interaction rules (e.g. memory and preference). Our computational experiments of this model suggest phase-transition-like phenomena between convergent and divergent evolution phases, which are controlled by the balance between convergent and divergent tendencies in local interactions. We formalize a conjecture on the latent phase transition mechanisms and summarize theoretical difficulties for verifying this conjecture.

To reach at a balance between the universality to fit in with different collective phenomena and the simpleness to support analytic derivations, we consider a n-dimensional abstract state space $S\subset\mathbb{R}^{n}$, where each element $\mathbf{s}\in S$ denotes an abstract state that can be used to model any target properties of system units (e.g. spatial coordinates, opinions, or genetic characters). A collective model updates the state of each unit $\sigma_{i}\in V$ at moment $t\in T$ according to its current state, denoted by $\mathbf{s}_{i}\left(t\right) = \mathbf{s}\left(\sigma_{i}\left(t\right)\right)$, and its interactions with its neighbor units, included by set N_i , following a local rule ψ

$$\begin{align} \mathbf{s}_{i}\left(t+\varepsilon\right)=\psi\left(\mathbf{s}_{i}\left(t\right),\mathbf{s}_{N_{i}}\left(t\right)\right). \end{align} \tag{ 1 }$$

Here $\mathbf{s}_{N_{i}}\left(t\right) = \{\mathbf{s}_{j}\left(t\right)\vert \sigma_{j}\in N_{i}\}$ stands for the states of neighbor units of unit i at moment t and ε denotes a minimum time step (see figure 1(a) for illustrations).

Zoom In Zoom Out Reset image size

The localization of interaction rules is realized by constraining N_i , the neighbor set of unit σ_i , as a sub-set of units whose states are close to σ_i

$$\begin{align} N_{i}=\{\sigma_{j}\vert d\left(\mathbf{s}_{i}\left(t\right),\mathbf{s}_{i}\left(t\right)\right)\leqslant\Delta\}\subset V, \end{align} \tag{ 2 }$$

where $d\left(\cdot,\cdot\right)$ denotes an arbitrary distance function and Δ is a parameter that determines the interaction range (see figure 1(a)). Please note that set N_i evolves across time according to equation (2) as the unit state is time-dependent.

A convergent evolution process is defined as a global pattern of the system, during which units become increasingly similar to each other in the state space. This process also corresponds to the case where the states of all units approach the global mean state. The non-triviality of this process lies in that each unit only interacts with its neighbors and, consequently, can never know the absolute coordinates of all units in the state space. Only the relative coordinate information obtained through interactions is available for units at each moment. For each unit σ_i , the relative coordinate of each of its neighbors in the state space is defined as

$$\begin{align} \mathbf{r}_{ij}\left(t\right)=\mathbf{s}_{j}\left(t\right)-\mathbf{s}_{i}\left(t\right). \end{align} \tag{ 3 }$$

Given this information, unit σ_i can only approach the time-dependent local mean state

$$\begin{align} \overline{\mathbf{s}}_{i}\left(t\right)=\mathbf{s}_{i}\left(t\right)+\frac{1}{\vert N_{i}\vert+1}\sum_{\sigma_{j}\in N_{i}}\mathbf{r}_{ij}\left(t\right) \end{align} \tag{ 4 }$$

during convergent evolution. Note that $\vert N_{i}\vert+1$ in the denominator of equation (4) is used to deal with the case with $\vert N_{i}\vert = 0$.

The above definitions have described the convergent evolution process of collective behaviors on a local scale. To study its non-trivial manifestation on a global scale, we need to relate equation (4) with specific global characteristics of the system.

We notice that equation (4) can be equivalently reformulated as

$$\begin{align} \overline{\mathbf{s}}_{i}\left(t\right)=\mathbf{s}_{i}\left(t\right)+\frac{1}{\vert N_{i}\vert+1}\sum_{\sigma_{j}\in V}\mathbf{A}_{ij}\left(t\right)\left(\mathbf{s}_{j}\left(t\right)-\mathbf{s}_{i}\left(t\right)\right), \end{align} \tag{ 5 }$$

where $\mathbf{A}\left(t\right)$ denotes the time-dependent adjacency matrix

$$\begin{align} \mathbf{A}_{ij}\left(t\right)=\left(1-\delta_{i,j}\right)\Theta\left[\Delta-d\left(\mathbf{s}_{i}\left(t\right),\mathbf{s}_{j}\left(t\right)\right)\right]. \end{align} \tag{ 6 }$$

In equation (6), notion $\delta_{\cdot,\cdot}$ denotes the Kronecker delta function and $\Theta\left(\cdot\right)$ is the Heaviside step function. By simple calculation, we can transform equation (4) as

$$\begin{align} &\overline{\mathbf{s}}_{i}\left(t\right)-\mathbf{s}_{i}\left(t\right)=\frac{1}{\operatorname{deg}\left(\sigma_{i}\left(t\right)\right)+1}\left(\sum_{\sigma_{j}\in V}\mathbf{A}_{ij}\left(t\right)\mathbf{s}_{j}\left(t\right)-\operatorname{deg}\left(\sigma_{i}\left(t\right)\right)\mathbf{s}_{i}\left(t\right)\right), \end{align} \tag{ 7 }$$

where $\operatorname{deg}\left(\cdot\right)$ denotes the degree in a graph defined by the adjacency matrix in equation (6).

If we consider the observable global state

$$\begin{align} &\mathbf{s}\left(t\right)=\left[\mathbf{s}_{1}\left(t\right),\ldots,\mathbf{s}_{\vert V\vert}\left(t\right)\right]^{\mathsf{T}}, \end{align} \tag{ 8 }$$

where $\mathsf{T}$ denotes the transpose, we can reformulate equation (7) into a new form

$$\begin{align} &\overline{\mathbf{s}}_{i}\left(t\right)-\mathbf{s}_{i}\left(t\right)=\frac{1}{\operatorname{deg}\left(\sigma_{i}\left(t\right)\right)+1}\left(\left[\mathbf{A}_{i1}\left(t\right),\ldots,\mathbf{A}_{i\vert V\vert}\left(t\right)\right]-\mathbf{d}_{i}\left(t\right)\right)\mathbf{s}\left(t\right) \end{align} \tag{ 9 }$$

based on vector $\mathbf{d}_{i}\left(t\right)$

$$\begin{align} &\mathbf{d}_{i}\left(t\right)=\left[\delta_{i,1}\operatorname{deg}\left(\sigma_{1}\left(t\right)\right),\ldots,\delta_{i,\vert V\vert}\operatorname{deg}\left(\sigma_{\vert V\vert}\left(t\right)\right)\right]. \end{align} \tag{ 10 }$$

We notice that vector $\mathbf{d}_{i}\left(t\right)$ in equation (9) happens to coincide with the ith row of the degree matrix in graph theory, where the ith entity is the degree of unit i at moment t while all other entities are zeros. Meanwhile, vector $\left[\mathbf{A}_{i1}\left(t\right),\ldots,\mathbf{A}_{i\vert V\vert}\left(t\right)\right]$ in equation (9) is known as the ith row of the time-dependent adjacency matrix. These observations remind us that $\left[\mathbf{A}_{i1}\left(t\right),\ldots,\mathbf{A}_{i\vert V\vert}\left(t\right)\right]-\mathbf{d}_{i}\left(t\right)$ is actually the ith row of matrix $-\mathbf{L}\left(t\right)$, where $\mathbf{L}\left(t\right)$ is the time-dependent Laplacian of the corresponding graph

$$\begin{align} \mathbf{L}_{ij}\left(t\right)=\delta_{i,j}\sum_{k}\mathbf{A}_{ik}\left(t\right)-\mathbf{A}_{ij}\left(t\right). \end{align} \tag{ 11 }$$

Taken together, the ith row of the time-dependent Laplacian $\mathbf{L}\left(t\right)$ can be naturally discovered in the formulation of dynamics when we analyze the state of unit i.

Given the above foundation, we can analyze the observable global state of the system in equation (8) applying the time-dependent Laplacian. Specifically, after defining the observable of all local mean states

$$\begin{align} &\overline{\mathbf{s}}\left(t\right)=\left[\overline{\mathbf{s}}_{1}\left(t\right),\ldots,\overline{\mathbf{s}}_{\vert V\vert}\left(t\right)\right]^{\mathsf{T}} \end{align} \tag{ 12 }$$

and the observable of the inverse of the degrees of all units

$$\begin{align} &\mathbf{q}\left(t\right)=\operatorname{diag}\left(\left[\frac{1}{\operatorname{deg}\left(\sigma_{1}\left(t\right)\right)+1},\ldots,\frac{1}{\operatorname{deg}\left(\sigma_{\vert V\vert}\left(t\right)\right)+1}\right]\right) \end{align} \tag{ 13 }$$

using the diagonal $\operatorname{diag}\left(\cdot\right)$, we can eventually derive a global version of equation (4)

$$\begin{align} \overline{\mathbf{s}}\left(t\right)-\mathbf{s}\left(t\right)=&-\mathbf{q}\left(t\right)\mathbf{L}\left(t\right)\mathbf{s}\left(t\right). \end{align} \tag{ 14 }$$

Please see figure 1(b) for illustrations. Equation (14) describes the states of all units in a unified manner while equation (9) serves as a special case of equation (14) where only the state of unit i is analyzed.

Now, let us consider the global dynamics characterized by equation (14). We define an abstract velocity matrix to quantify the intrinsic degrees of the willingness of units to exhibit convergent evolution in the state space

$$\begin{align} &\mathbf{v}_{\textrm c}\left(t\right)=\left[\mathbf{v}_{\textrm c,1}\left(t\right),\ldots,\mathbf{v}_{\textrm c,\vert V\vert}\left(t\right)\right]\mathbf{1}^{\mathsf{T}}, \end{align} \tag{ 15 }$$

where 1 denotes an all-one row vector of an appropriate size that is applied to ensure the consistency of matrix dimensions. Under the assumption of time continuity, the dynamics for $\mathbf{s}\left(t\right)$ to approach $\overline{\mathbf{s}}\left(t\right)$ during convergent evolution is formalized as

$$\begin{align} \frac{\partial}{\partial t}\mathbf{s}\left(t\right)=-\mathbf{v}_{\textrm c}\left(t\right)\odot\mathbf{q}\left(t\right)\mathbf{L}\left(t\right)\mathbf{s}\left(t\right), \end{align} \tag{ 16 }$$

where local topology effects are described by $\mathbf{q}\left(t\right)$ (i.e. how the neighbors of a unit affect its evolution), and $\odot$ denotes the Schur product (or referred to as the element-wise product between matrices, i.e. we have $\mathbf{Z}_{ij} = \mathbf{X}_{ij}\mathbf{Y}_{ij}$ if $\mathbf{Z} = \mathbf{X}\odot\mathbf{Y}$ for arbitrary matrices X, Y, and Z of the same size). In equation (16), the time-dependent Laplacian $\mathbf{L}\left(t\right)$ serves as a bridge to make the dynamic evolution of global system states incorporate with the local interactions among units described by graph topology.

Then we turn to analyze the divergent evolution processes of collective behaviors. At the first glance, one only needs to consider the intrinsic degrees of the willingness of units to exhibit divergent evolution

$$\begin{align} &\mathbf{v}_{\textrm d}\left(t\right)=\left[\mathbf{v}_{\textrm d,1}\left(t\right),\ldots,\mathbf{v}_{\textrm d,\vert V\vert}\left(t\right)\right]\mathbf{1}^{\mathsf{T}} \end{align} \tag{ 17 }$$

and assume that units depart from associated time-dependent local mean states by moving in a direction opposite to equation (16). Under this condition, the joint dynamics of convergent and divergent evolution is nothing more than a Newtonian motion in the state space, $\frac{\partial}{\partial t}\mathbf{s}\left(t\right) = -\left[\mathbf{v}_{\textrm c}\left(t\right)-\mathbf{v}_{\textrm d}\left(t\right)\right]\odot\mathbf{q}\left(t\right)\mathbf{L}\left(t\right)\mathbf{s}\left(t\right)$.

However, it is unreasonable for our framework to be so trivial because Newtonian mechanics may be impractical in the abstract state space. The information of the mentioned opposite direction can be unavailable in real cases. For instance, if the state space describes the preference of drinks, it is impossible to precisely define an opposite direction of preferring coffee, i.e. it is difficult for units to tell preferring which drink is opposite to preferring coffee. This is because the definition of an opposite drink of coffee is ill-posed, at least in most situations. To describe these non-negligible cases where Newtonian mechanics is unacceptable, we need to consider a weaker assumption.

Our idea arises from the relation between $\mathbf{L}\left(t\right)$ and $-\mathbf{L}^{\dagger}\left(t\right)$, where $\dagger$ denotes the Moore–Penrose pseudoinverse [34]. To understand their relation, let us consider an instance where the covariance matrix units is defined as $\mathbf{C}\left(t\right) = \mathbf{L}\left(t\right)$. Two units are expected to evolve inversely (i.e. they have a negative covariance) if they are connected by an edge in the graph defined by equation (6). In this case, adjusting edges is equivalent to adjusting the differences (i.e. gradients) between units. Therefore, equation (16) can utilize $-\mathbf{L}\left(t\right)$ in the derivative to minimize those differences. If we further consider the associate precision matrix $\mathbf{Q}\left(t\right): = \mathbf{C}^{\dagger}\left(t\right) = \mathbf{L}^{\dagger}\left(t\right)$, we can relate the partial correlation [35, 36] between units σ_i and σ_j with $\mathbf{L}^{\dagger}\left(t\right)$

$$\begin{align} \operatorname{corr}\left(\sigma_{i},\sigma_{j}\vert V\setminus\{\sigma_{i},\sigma_{j}\}\right) & = -\frac{\mathbf{Q}_{ij}\left(t\right)}{\sqrt{\mathbf{Q}_{ii}\left(t\right)\mathbf{Q}_{jj}\left(t\right)}}, \end{align} \tag{ 18 }$$

$$\begin{align} & = -\frac{\mathbf{L}^{\dagger}_{ij}\left(t\right)}{\sqrt{\mathbf{L}^{\dagger}_{ii}\left(t\right)\mathbf{L}^{\dagger}_{jj}\left(t\right)}}. \end{align} \tag{ 19 }$$

Based on equation (19), two units are expected to evolve consistently (i.e. have a positive partial correlation) if they correspond to a negative value in $\mathbf{L}^{\dagger}\left(t\right)$ (i.e. a positive value in $-\mathbf{L}^{\dagger}\left(t\right)$). In this case, adjusting edges is equivalent to controlling similarities between units. During divergent evolution, units are expected to minimize their similarities, which can be realized by the term $\mathbf{L}^{\dagger}\left(t\right)$ in the derivative of the following dynamics

$$\begin{align} \frac{\partial}{\partial t}\mathbf{s}\left(t\right)=\mathbf{v}_{\textrm d}\left(t\right)\odot\mathbf{q}^{-1}\left(t\right)\mathbf{L}^{\dagger}\left(t\right)\mathbf{s}\left(t\right). \end{align} \tag{ 20 }$$

In equation (20), we use $\mathbf{q}^{-1}\left(t\right)$, the inverse matrix of $\mathbf{q}\left(t\right)$, to characterize local topology effects.

Now, we have the opportunity to define the joint global dynamics of convergent and divergent evolution based on equations (16) and (20). Specifically, the joint dynamics is given as

$$\begin{align} \frac{\partial}{\partial t}\mathbf{s}\left(t\right)=&\big[\mathbf{v}_{\textrm d}\left(t\right)\odot\mathbf{q}^{-1}\left(t\right)\mathbf{L}^{\dagger}\left(t\right)-\mathbf{v}_{\textrm c}\left(t\right)\odot\mathbf{q}\left(t\right)\mathbf{L}\left(t\right)\big]\mathbf{s}\left(t\right), \end{align} \tag{ 21 }$$

which is essentially a type of generalized Laplacian dynamics with more complicated details (see [37–39] for the elementary Laplacian dynamics). See figure 1(c) for illustrations. Although the idea underlying equation (21) is natural and simple, analyzing the dynamic properties of equation (21) is intricate because the competition between convergent and divergent evolution is characterized by the non-trivial interaction between $\mathbf{L}\left(t\right)$ and $\mathbf{L}^{\dagger}\left(t\right)$ under the effects of local topology. To offer a clear vision, we primarily consider a special case of equation (21) in our subsequent analysis.

The considered special case of equation (21) is a homogeneous system of units that exhibits uniform motions in the state space. To create uniform motions (i.e. with constant velocities) and make the system homogeneous (here, being homogeneous refers to the case where all units share the same intrinsic properties), we need to define

$$\begin{align} \mathbf{v}_{\textrm c}\left(t\right)\equiv x\mathbf{1}^{\mathsf{T}}\mathbf{1},\;\forall t\in T, \end{align} \tag{ 22 }$$

$$\begin{align} \mathbf{v}_{\textrm d}\left(t\right)\equiv y\mathbf{1}^{\mathsf{T}}\mathbf{1},\;\forall t\in T, \end{align} \tag{ 23 }$$

where $x,y\in\mathbb{R}$ are specific constants that we can manipulate (e.g. used as controlling factors). Equations (22) and (23) make all units share the same and uniform intrinsic willingness degrees, x or y, to exhibit convergent or divergent evolution.

Therefore, in a homogeneous system with uniform motions, equation (21) reduces to a simple form

$$\begin{align} \frac{\partial}{\partial t}\mathbf{s}\left(t\right)& =\left[y\mathbf{q}^{-1}\left(t\right)\mathbf{L}^{\dagger}\left(t\right)-x\mathbf{q}\left(t\right)\mathbf{L}\left(t\right)\right]\mathbf{s}\left(t\right), \end{align} \tag{ 24 }$$

$$\begin{align} & :=\mathbf{H}\left(t\right)\mathbf{s}\left(t\right) \end{align} \tag{ 25 }$$

which is easier to analyze. In computational experiments, we apply the following numerical iteration

$$\begin{align} \mathbf{s}\left(t+\varepsilon\right)&=\varepsilon\mathbf{H}\left(t\right)\mathbf{s}\left(t\right)+\mathbf{s}\left(t\right) \end{align} \tag{ 26 }$$

to derive the results of equation (25). Note that equation (26) only serves as a simple illustration. One can also consider more complicated numerical approaches.

For the sake of visualization, we mainly consider a 1-dimensional state space in our work, yet a state space with higher dimensions can be directly analyzed without any modification. In figures 1(d)–(l), we present diverse instances of this homogeneous system with different parameter settings, where a regular initialization condition of the system is applied to offer a clear vision (i.e. initialized units are uniformly distributed and form clear cluster structures in the state space). As shown in figures 1(d)–(i), different ratios between x and y in equation (24), i.e. different ratios between the intrinsic willingness degrees of convergent and divergent evolution, create distinct trends of collective behaviors. A larger x can lead to convergent evolution while a larger y makes units divergent. As suggested by figures 1(j)–(l), the change of local interaction range controlled by a varying Δ in equation (2) leads to distinct final patterns in the state space after evolution. A larger local interaction range can amplify the trends of convergent and divergent evolution because more units are involved. These phenomena can also be observed in figure 2, where a random initialization condition is applied to the system (i.e. the initial states of units are randomly sampled following a specific distribution). These results provide qualitative illustrations of the global dynamics of collective behaviors and suggest how the dynamics is shaped by the competition between convergent and divergent evolution trends or by the local interaction range.

Zoom In Zoom Out Reset image size

To quantitatively analyze the convergent and divergent evolution of collective behaviors, we study the properties of the eigenvalue spectrum of $\mathbf{H}\left(t\right)$. Specifically, we analyze the standardized moments of the eigenvalue spectrum to reflect the latent dynamics of the evolution process (note that the first and the second standardized moment are not considered because they are constant [40]). Meanwhile, we consider $d_{k}\left(\mathbf{s}_{i}\left(t\right)\right)$, the distance between state $\mathbf{s}_{i}\left(t\right)$ and its k-nearest state at moment t given a distance function $d\left(\cdot,\cdot\right)$ in equation (2). We define the average k-nearest distance at moment t as

$$\begin{align} \mathbb{E}\left(d_{k}\left(t\right)\right)=\frac{1}{\vert V\vert}\sum_{\sigma_{i}\in V}d_{k}\left(\mathbf{s}_{i}\left(t\right)\right) \end{align} \tag{ 27 }$$

to characterize the global state of the system at moment t. In figure 3, these two kinds of metrics are calculated on three representative instances of the homogeneous system in equation (25). As shown in figure 3, distance $\mathbb{E}\left(d_{k}\left(t\right)\right)$ exhibits distinct changes during convergent and divergent evolution (we define $d\left(\cdot,\cdot\right)$ as the Euclidean distance for convenience). Before the dynamics converges to its limit, the standardized moments of the eigenvalue spectrum dramatically increase during divergent evolution while decreasing during divergent evolution.

Zoom In Zoom Out Reset image size

The results shown in figure 3 inspire us to distinguish between convergent and divergent phases using $\frac{\mathbb{E}\left(d_{k}\left(t\right)\right)}{\mathbb{E}\left(d_{k}\left(0\right)\right)}$, the relative ratio between the average k-nearest distances at moment 0 and moment t. This ratio is larger than one in the divergent phase while it is smaller than one in the convergent phase. In figure 4, we analyze the phase space of the homogeneous system in terms of $\frac{\mathbb{E}\left(d_{k}\left(t\right)\right)}{\mathbb{E}\left(d_{k}\left(0\right)\right)}$ under different conditions of x, y, and Δ. As shown in figures 4(a)–(f), convergent and divergent phases share a non-trivial and relatively blurry boundary on the plane of $\left(x,y\right)$, which changes across different local interaction range conditions. Phase-transition-like phenomena can be observed if we consider $\frac{\mathbb{E}\left(d_{k}\left(t\right)\right)}{\mathbb{E}\left(d_{k}\left(0\right)\right)}$ as an order parameter and treat the ratio between the intrinsic willingness degrees of convergent and divergent evolution in equation (24), $\frac{x}{y}$, as a control parameter in figures 4(g)–(l). A kind of phase-transition-like processes from the divergent phase to the convergent phase can be observed as $\frac{x}{y}$ increases, during which the collective behaviors of units experiences drastic changes. Specifically, before $\frac{x}{y}$ reaches a certain value, all units generally exhibit divergent evolution. In this phase, units tend to be further away from their neighbors as time goes by (i.e. become increasingly different). After $\frac{x}{y}$ exceeds a specific threshold, units change to exhibit convergent evolution. In this phase, units have stronger tendencies to approach to their neighbors in the state space (i.e. become increasingly similar).

Zoom In Zoom Out Reset image size

Here we emphasize that these processes should be treated as phase-transition-like phenomena rather than strict phase transitions. It remains elusive whether these phenomena really arise from specific underlying phase transitions or from other unknown mechanisms with similar appearances. Diverse theoretical difficulties are met when we attempt to confirm the detailed properties of these phase-transition-like phenomena (e.g. the order of transition and scaling relation), which impede us from verifying the potential existence of the phase transition. Therefore, the existence of potential phase transition in the Laplacian dynamics of convergent and divergent collective behaviors remains as a hypothesis. Below, we summarize the main difficulties of studying this conjecture.

The main challenges for studying the observed phase-transition-like phenomena arise from the non-triviality of predicting the limiting behaviors of the homogeneous system (i.e. the classification of convergent and divergent phases depends on the limit state of the system after the evolution process converges).

Although the motions of units in the state space have been constrained as uniform, the homogeneous system in equation (25) is shaped by complicated local topology effects, which is in contrast to a more theoretically favorable case where local topology effects are completely excluded

$$\begin{align} \frac{\partial}{\partial t}\mathbf{s}\left(t\right)&=\left[y\mathbf{L}^{\dagger}\left(t\right)-x\mathbf{L}\left(t\right)\right]\mathbf{s}\left(t\right):=\mathbf{M}\left(t\right)\mathbf{s}\left(t\right). \end{align} \tag{ 28 }$$

In such an ideal case, we have the opportunity to express $0,\omega_{2}\left(t\right)\leqslant\cdots\leqslant\omega_{\vert V\vert}\left(t\right)$, the eigenvalues of matrix $\mathbf{M}\left(t\right)$, based on the eigenvalues $0\leqslant \lambda_{2}\left(t\right)\leqslant\cdots\leqslant\lambda_{\vert V\vert}\left(t\right)$ of matrix $\mathbf{L}\left(t\right)$

$$\begin{align} \omega_{i}\left(t\right)=\frac{y}{\lambda_{i}\left(t\right)}-x\lambda_{i}\left(t\right),\;\forall i\geqslant 2. \end{align} \tag{ 29 }$$

Meanwhile, matrix $\mathbf{M}\left(t\right)$ shares the same group of eigenvectors, $\mathbf{1}^{\mathsf{T}},\mathbf{u}_{2}\left(t\right),\ldots,\mathbf{u}_{\vert V\vert}\left(t\right)$, as matrix $\mathbf{L}\left(t\right)$. Equation (29) is derived from the fact that the eigenvalues and eigenvectors of matrix $\mathbf{M}\left(t\right)$ are $0,\frac{1}{\lambda_{2}\left(t\right)},\ldots,\frac{1}{\lambda_{\vert V\vert}\left(t\right)}$ and $\mathbf{1}^{\mathsf{T}},\mathbf{u}_{2}\left(t\right),\ldots,\mathbf{u}_{\vert V\vert}\left(t\right)$, respectively [41, 42]. The explicit eigen-decomposition of $\mathbf{M}\left(t\right)$ provides rich information about the evolution dynamics, making it possible to analyze equation (28) following the dynamic system theory [42]. However, this possibility vanishes due to the local topology effects in equation (25). To quantify these non-trivial effects, we can measure the cosine distance between the vector of eigenvalues of $y\mathbf{L}^{\dagger}\left(t\right)$ and that of $y\mathbf{q}^{-1}\left(t\right)\mathbf{L}^{\dagger}\left(t\right)$ (similarly, between the vector of eigenvalues of $x\mathbf{L}\left(t\right)$ and that of $x\mathbf{q}\left(t\right)\mathbf{L}\left(t\right)$). Meanwhile, we can calculate the average cosine distance between the eigenvectors of $y\mathbf{L}^{\dagger}\left(t\right)$ and those of $y\mathbf{q}^{-1}\left(t\right)\mathbf{L}^{\dagger}\left(t\right)$ (similarly, between the eigenvectors of $x\mathbf{L}\left(t\right)$ and those of $x\mathbf{q}\left(t\right)\mathbf{L}\left(t\right)$). These cosine distances reflect the local topology effects on the eigenvalues and eigenvectors of $x\mathbf{L}\left(t\right)$ and $y\mathbf{L}^{\dagger}\left(t\right)$. In figures 5(a) and (b), we present an instance calculated on a representative homogeneous system, suggesting that these non-negligible local topology effects are volatile until the evolution dynamics converges. Therefore, equation (29) departs from the actual properties of equation (25) in most cases.

Zoom In Zoom Out Reset image size

Moreover, another challenge arises from that $\mathbf{H}\left(t\right)$ in equation (25) is not commutative until the evolution dynamics converges to its limit. We can measure the deviation of $\mathbf{H}\left(t\right)$ from being commutative based on

$$\begin{align} \chi\left(t,t^{^{\prime}}\right)=\sum_{i,j}\vert\left[\mathbf{H}\left(t\right),\mathbf{H}\left(t^{^{\prime}}\right)\right]_{ij}\vert, \end{align} \tag{ 30 }$$

where $\left[\cdot,\cdot\right]$ denotes the commutator and $\left[\cdot,\cdot\right]_{ij}$ is the $\left(i,j\right)$th element of the commutator. We refer to $\chi\left(\cdot,\cdot\right)$ as the commutator error. In figure 5(c), we calculate an instance of $\chi\left(\cdot,\cdot\right)$ by setting $t^{^{\prime}} = t+\varepsilon$. The derived results suggest that $\mathbf{H}\left(t\right)$ becomes trivially commutative only after it arrives at its limit. The absence of a commutative $\mathbf{H}\left(t\right)$ in most cases makes it invalid to express the solution of equation (25) as $\mathbf{s}\left(t\right) = \exp\left(\int_{0}^{t}\mathbf{H}\left(\tau\right)\mathsf{d}\tau\right)\mathbf{s}\left(0\right)$, rejecting the possibility of simplifying equation (25).

In sum, the above properties create obstacles for theoretically analyzing the homogeneous system in equation (25). The situation becomes even more non-trivial if we relax the constraints on $\mathbf{v}_{\textrm c}$ and $\mathbf{v}_{\textrm d}$ to consider heterogeneous cases. These challenges may serve as valuable directions for future theoretical studies.

In this study, we propose a minimal model to describe the convergent and divergent evolution of collective behaviors. Despite its mathematical simplicity, the model exhibits abundant non-trivial behaviors and suggests the potential existences of phase-transition-like phenomena between convergent and divergent evolution controlled by local interactions. The non-triviality of the proposed model may inspire new mathematical techniques in analyzing the physics of collective dynamics. Moreover, the simplicity of the proposed model in computational aspects enables it to serve as a foundation for developing large-scale numerical simulations of collective behaviors.

While the current work primarily focuses on the theoretical analysis of the proposed model, more explorations in model applications can be considered in future studies. The proposed model may find its place in studying the social dynamics of people and animals [1, 2, 8] (note that numerous animals, such as chimpanzees, have complex social behaviors). For instance, the state space can be used to describe opinions, where the Newtonian mechanics may be not applicable (i.e. opinions are usually vague and do not support a precise definition of the opposite direction) [8]. The Laplacian dynamics of convergent and divergent evolution of opinions can be used to model the emergence and collapse of consensus in social networks. In the divergent phase, units tend to differentiate their opinions to improve the uniqueness of their attitudes. In the convergent phase, units attempt to forge consensus and form a certain number of clusters where each cluster of units share a representative opinion. In somewhere between these two phases, units may seek common points while reserving differences in their opinions. In general cases where the intrinsic willingness degrees of convergent or divergent evolution are not constrained as homogeneous and constant but are allowed to vary (e.g. let $\mathbf{v}_{\textrm c}\left(t\right)$ and $\mathbf{v}_{\textrm d}\left(t\right)$ change across time or be correlated with unit states), the proposed model may be able to reproduce the intricate evolution of public opinions. A possible way to relate our model with existing models in opinion dynamics is to make $\mathbf{v}_{\textrm c}\left(t\right)$ and $\mathbf{v}_{\textrm d}\left(t\right)$ time-dependent and correlated with the size of the neighbor set. As an illustration, one can let $\mathbf{v}_{\textrm c,i}\left(t\right)$ or $\mathbf{v}_{\textrm d,i}\left(t\right)$ be positively correlated with $\operatorname{deg}\left(\sigma_{i}\left(t\right)\right)$ for every unit σ_i . Under this condition, the actual velocity of convergent or divergent evolution of a unit critically depends on the size of cluster where this unit belongs to, making it possible to describe social impacts (e.g. how convincing power depends on the size of social groups [8]) by our model. To this aim, one can apply our framework to study the opinion dynamics characterized by diverse classic models. In the case where opinions are discrete, one can discretize the state space and apply our model to study the cellular automata developed by Latané [43] and Nowak et al [44], whose statistical mechanics has been previously analyzed at the mean-field limit [45] but remains elusive in more general cases. Given the extensive applications of this cellular automata and its variants in studying social opinions [43, 44], learning [46], the effects of leaders on opinion formation [47–49], and the mitigation of social impacts due to population diversity [50] (see [51] for a review), our model may serve as an alternative choice for formalizing these collective behaviors apart from the mean-filed approximation. Except for the cellular automata mentioned above, the Sznajd dynamics [52–57] can also be analyzed applying our model, which suggests the possibility for our model to be used in studying opinion spreading under social impacts [52–54], voting [58], and the interaction of economics and opinions [54]. In the case where opinions are continuous, our model can also be defined on the dynamics generated by the Deffuant model [59], the Hegselmann–Krause model [60], and their density-based variants (see [61] for a review). In these models, an assumption referred to as the bounded confidence is applied, implying that units only have effective interactions with peers with similar opinions [61]. When the differences between opinions are not tolerable, interactions become non-effective [61]. This property can be naturally and equivalently described by the local interaction rule in our model. Apart from opinion dynamics, other topics in social dynamics (see a summary in [8]) may also be analyzed by our model, which remains for future explorations.

This project is supported by the Artificial and General Intelligence Research Program of Guo Qiang Research Institute at Tsinghua University (2020GQG1017) as well as the Tsinghua University Initiative Scientific Research Program.

All data that support the findings of this study are included within the article (and any supplementary files).