Opinion formation on social networks with algorithmic bias: dynamics and bias imbalance

We investigate opinion dynamics and information spreading on networks under the influence of content filtering technologies. The filtering mechanism, present in many online social platforms, reduces individuals’ exposure to disagreeing opinions, producing algorithmic bias. We derive evolution equations for global opinion variables in the presence of algorithmic bias, network community structure, noise (independent behavior of individuals), and pairwise or group interactions. We consider the case where the social platform shows a predilection for one opinion over its opposite, unbalancing the dynamics in favor of that opinion. We show that if the imbalance is strong enough, it may determine the final global opinion and the dynamical behavior of the population. We find a complex phase diagram including phases of coexistence, consensus, and polarization of opinions as possible final states of the model, with phase transitions of different order between them. The fixed point structure of the equations determines the dynamics to a large extent. We focus on the time needed for convergence and conclude that this quantity varies within a wide range, showing occasionally signatures of critical slowing down and meta-stability.


Introduction
The collective behavior of a system made of interacting individuals can be successfully analyzed using agent-based models, as shown in many examples across various disciplines [1][2][3].In these models, individuals (or agents) are often pictured as nodes in a network [4,5], where the links represent the possible interactions between them.Each node holds a dynamical state variable whose interpretation depends on the context of the model.In opinion dynamics this opinion variable [1] can be considered as the political party preference of an individual (e.g., liberal or conservative), her inclination towards or against some regulation or initiative, etc.The usefulness of the approach lies in the simplicity of the setting, leading, nevertheless, to complex phenomena due to collective effects.The possibility of considering various elements that are hypothesized to be relevant for opinion formation (such as structural and dynamical heterogeneities) deeply improves our understanding of the underlying social phenomena.
In recent years, human communication has changed dramatically, and moved from traditional media (face-to-face, phone, or mass media like television and the press) to online social media platforms (Google, Twitter, Facebook, etc.) [6,7].In contrast to earlier media channels, online social networks control the information that users see and send to each other by means of personalized filtering algorithms [8].These algorithms record individual information about users' preferences and then filter incoming data accordingly [9,10].Hence, people tend to be exposed to opinions they already agree with, producing so-called algorithmic bias.This reinforcement feature changes, ultimately, the global behavior of the population [13][14][15][16], promoting phenomena like "filter bubbles" or "echo chambers", where people divide in groups with opposing views that barely interact with each other.Explaining how and when the polarization of opinion groups emerges within a population is of crucial importance [17,18].It is of particular interest to understand how copying or herding processes (typical signatures of human social behavior), when coupled to algorithmic bias, may enhance or decrease polarization.These two ingredients can be implemented in the formalism of agent-based models, leading to a flexible mathematical framework open to analytical and numerical treatments.Ultimately, the results of the models can be interpreted to help devise potential strategies to mitigate the negative effects of algorithmic bias.
Previous modeling efforts [19] have been made to consider algorithmic bias in bounded confidence models [20], where the opinion variable of individuals is a real continuous variable on an interval.The filtering algorithm requires the opinions of two individuals to be similar enough to be able to interact, and bias means that similar people with similar opinions have a greater chance to meet, leading to enhanced polarization and fragmentation in opinion space.Another class of models considers opinions to be discrete (a binary variable in the simplest case) [21,22].Perra and Rocha [23] have studied algorithmic bias in such a model by considering that the opinion of an individual is influenced by its neighbors in the network, and filtered in various ways.An alternative implementation of algorithmic bias in binary-state models has been proposed in [24].In this case, the social platform records information about all the previous opinions of individuals, and then influences them to keep the opinion that has been held for the longest time, similarly to a memory or 'inertia' effect [25][26][27].All these implementations of algorithmic bias in opinion dynamics modeling suggest that polarization is a consequence of both the social behavior of individuals and the content filtering algorithms constraining their actions.
In the present work, similarly to the approach of [23], we consider that a fraction of the neighbors of an individual holding disagreeing opinions are filtered, and thus interactions with those neighbors are not possible.Recently we have proposed a general formalism within the binary-state approach that includes this implementation of algorithmic bias [28].We have extended previous theoretical tools to describe the macroscopic dynamics on networks, including mean-field and pair approximations.We have also explored modular community structures, a crucial ingredient to characterize opinion polarization, i.e. the division of the population into opinion groups.We have studied the static, asymptotic behavior of archetypal models of opinion formation in the presence of algorithmic bias and concluded that, systematically, pairwise interactions lead to polarization, while group interactions promote coexistence of opinions.Here we use the same formalism of [28], but focus on dynamical aspects of the opinion formation process.We also extend the algorithmic bias mechanism to consider situations in which the online social platform promotes one opinion over the other, thus unbalancing the dynamics in favor of the preferred opinion by the platform, in what we call bias asymmetry.
We use a prototypical model of social behavior, the language model [29,30], which takes into account both pairwise and group-based (copying) interactions depending on the value of a tunable parameter α.We also include the possibility that individuals act independently of their neighbors [31][32][33], which we denote as noise with intensity Q.This general framework enables us to consider several interaction mechanisms and leads to various opinion formation scenarios.In [28], we have considered other archetypal models of opinion formation (including voter-like and majority-vote models), and realized that the language model essentially interpolates between dynamics with either pairwise or group interactions depending on the value of α, and is thus a good candidate to explore the effects of bias asymmetry within a single model.
The paper is organized as follows.In Sec. 2 we define the opinion formation model, algorithmic bias, and social network with community structure we use.In Sec. 3 we derive a set of mean-field equations that describe the global dynamics of the model, and derive its associated fixed points (stationary states of the dynamics).In Sec. 4 we explore the local dynamics and stability of the fixed points and build the phase diagram of the model.In Sec. 5 we present a detailed study of temporal behavior of the considered opinion dynamics model using numerical simulations and some theoretical tools, with particular emphasis on the role of initial conditions, and the behavior of the time to reach the final state.Throughout this paper we will pay special attention to the effect of algorithmic bias and its asymmetry.

Model and definitions
We consider the formalism of binary-state dynamics [21,22] as basic ground for modeling opinion formation.The model system is composed of a set of i = 1, . . ., N individuals, each one holding a binary-state (opinion) variable s i (t) = 0, 1 at time t (e.g., liberal or conservative in a political setting).We define the macroscopic state (global opinion) of the system as ρ = N −1 N i=1 s i ∈ [0, 1], i.e. the density of individuals in state 1.Individuals are represented by nodes of an (undirected) network, and links in the network correspond to some social relationship between them, such that the opinion of an individual can be influenced by its neighbors in the network.The state of node i changes according to rates depending on the specific dynamics and the network structure: with "infection" rate F k i ,m i from s i = 0 → 1, and with "recovery" rate R k i ,m i from s i = 1 → 0, where k i is the degree of the node in the network and m i ∈ [0, k i ] is the number of (nearest) neighbors of i in state s j = 1.(Here the names of the rates refer to the analogy with epidemic spreading.) In the following we specify the spreading dynamics, i.e. the functional form of the rates F k,m and R k,m , and the network structure.As for the spreading dynamics, we incorporate algorithmic bias (representing content filtering as implemented in many online social platforms) that influences and controls the way people interact.For the network structure, we include modules or communities by dividing the population into groups with tunable connectivity.

Transition rates in the language model
As mentioned above, we focus on the language model [29,30], which is able to describe both pairwise and group interactions.The transition rates are as follows: with Q ∈ [0, 1/2] and α ∈ (0, ∞) as tuning parameters.The model takes into account two mechanisms driving the dynamics: (i) noisy or idiosyncratic changes of state, with intensity Q; and (ii) herding or copying the states of neighbors with probability proportional to the fraction of neighbors in the opposite state to a power α.The rates in Eqs.(1,2) were first studied in the case without noise (Q = 0) to model the dynamics of language death [29], but the same model has been applied to other types of social human behavior, for example in the context of opinion formation in social media [34].The role of noise (Q > 0) has been extensively studied of late, as for the non-linear noisy voter model in [33] with α a real (continuous) number, and for the q-voter model [35][36][37] with α = q a positive integer.The language model encapsulates a wide variety of phenomena depending on the value of α.E.g., for the mean field (complete network) the following regions can be distinguished (i) low 0 < α < 2 (pairwise interactions); (ii) high 2 < α < 5, (group interactions); and (iii) very high α > 5.Each of these cases displays a distinct phenomenology [28,33] and represents a different archetypal way for humans to influence each other (either in pairs or in groups).
Note that the original rates in Eqs.(1, 2) fulfill the "up-down symmetry" condition R k,m = F k,k−m but, as we will show next, Eqs.(3,4) are only symmetric for b 0 = b 1 .In other words, unbalanced algorithmic bias breaks the symmetry of the system and favors one opinion over the other.

Algorithmic bias
A simple implementation of algorithmic bias has been proposed by us in a previous study [28].Here we generalize the definition of that paper by introducing two parameters characterizing the bias, instead of one.The bias intensities b 0 and b 1 (where the subscripts refer to state 0 or 1) take values in the interval [0, 1].These parameters measure the probabilities that the online platform filters out a neighbor in the opposite state, 0 or 1, (disagreeing opinion) of an individual with a given opinion, so that further interactions with that neighbor cannot take place.This mechanism of content filtering can be implemented in the formalism of any rate governed binary-state model by considering the following effective transition rates with the binomial (3,4) express the average rates of changing state, after removing with probability b 0 or b 1 a subset of neighbors in the opposite state (b 0 if neighbors are in state s = 0 and b 1 for s = 1).We define the total bias intensity b = (b 0 + b 1 )/2, and the bias asymmetry ∆b = b 1 − b 0 , with b 0 = b − ∆b/2 and b 1 = b + ∆b/2.For ∆b > 0 the social platform favors s = 0, while for ∆b < 0 it favors s = 1.In [28] we have implemented bias in a similar way to Eqs. (3,4), but with b 0 = b 1 = b (i.e.∆b = 0), such that the role of algorithmic bias is symmetric, or balanced across opinions.

Modular structure
The social network of interactions between individuals is fully specified by the adjacency matrix A ij , with elements equal to 1 if i and j are connected and 0 otherwise.For simplicity we consider the degree distribution P k as the only relevant structural feature of the network, with average degree z = k P k k.We use the standard configuration model [38] to produce synthetic networks in the corresponding numerical simulations.
When the network displays modular (community) structure, we can divide the population into groups with higher connectivity to nodes inside the group than to those outside.Assuming for simplicity two modules, nodes i = 1, . . ., N 1 are in group 1 of size N 1 , and nodes i = N 1 + 1, . . ., N 1 + N 2 = N in group 2 of size N 2 .The two groups have different connectivity depending on whether links join nodes of the same group or between groups.In this way, two nodes in the same group are more likely to be connected than if they belong to different groups.This group homophily [39,40] is a property that many real-world social networks show and is strongly correlated with the opinion of individuals.In other words, the opinion of an individual is ultimately related to the group it belongs to, in what we call a polarization state, or homophily for those who share the same ideas as one.In order to characterize the macroscopic state accordingly, we need at least two variables: , the density of nodes in state 1 in groups 1 and 2, respectively, with total ρ , which measures the degree of opinion dissimilarity between groups.
We consider four average degrees, z 1 , z 12 , z 21 , and z 2 , defining the connectivity inside and between groups.Parameter z 1 (z 2 ) is the average degree only considering links that join nodes within group 1 (2), while z 12 (z 21 ) is the average degree only considering links that depart from group 1 and end up in group 2 (from 2 to 1).The total number of links that go from group 1 to 2 is the same as those that go from 2 to 1, so we have the constraint N 1 z 12 = N 2 z 21 .

Mean-field description
We derive a set of mean-field evolution equations for the two macroscopic variables ρ 1 (t) and ρ 2 (t), of general validity in the thermodynamic (N → ∞) and highly connected (z 1 , z 12 , z 21 , z 2 → ∞) limit, with constant ratios of the average degrees.Even in a finite system with high connectivity, the mean field description is a good approximation of the dynamics, and it captures the phenomenology of the model well.For large values of N , the opinion variables fluctuate slightly around their average values, i.e., ρ 1 (t) ≈ ρ 1 (t) , ρ 2 (t) ≈ ρ 2 (t) .Thus, throughout the following mean-field description, ρ 1 (t) and ρ 2 (t) refer to average values over realizations.
In order to derive the mean-field equations [28] we first define the average rate of changing state [22] as where x is the probability of finding a neighbor in state 1, and P k k/z is the probability that a link connects to a node with degree k.In order to obtain a closed description of the dynamics, we must relate the probability x to the description variables ρ 1 and ρ 2 .The probability x depends on the group to which the node belongs.In the case of group 1 we have with p 1 = N 2 z 21 /N 1 z 1 = z 12 /z 1 , and, similarly, for group 2 exchanging the labels 1 ↔ 2, with p 2 = N 1 z 12 /N 2 z 2 = z 21 /z 2 .Eq. ( 6) is the ratio of the number of links coming out of nodes in state 1 that connect to nodes in group 1, and the number of links coming out of nodes in group 1.
If we consider algorithmic bias (b 1 > 0), we must use the effective F * k,m of Eq. ( 3) instead of F k,m to calculate the average rate f * [x] in the network using Eq. ( 5).Applying an argument based on the highly connected limit (z → ∞) [28], the effective average rate with bias reduces to When An analogous procedure can be applied to the effective recovery rate R * k,m of Eq. ( 4), leading to r * [x], which in the presence of up-down After defining the effective average rates f * [x], r * [x] and the probabilities x 1,2 , we obtain a system of two differential (mean-field) equations for the dynamics of the state variables ρ(t) with components ρ 1 (t) and ρ 2 (t): This mean-field description of the opinion formation model has the social behavioral parameters (Q, α), the algorithmic bias parameters (b, ∆b), and the group connectivity parameters (p 1 , p 2 ), together with the initial condition ρ 1 (0), ρ 2 (0).

Fixed point structure and stationary solutions
The stationary results of Eqs.(9, 10), i.e. the infinite time limit ρ 1 (∞) and ρ 2 (∞), correspond to the stable fixed points.The fixed points ρ st 1 , ρ st 2 are obtained from the condition The study of the solutions of Eqs.(11,12) as a function of the parameters is a first step in understanding the dynamics and general behavior of the model.In Fig. 1 we show the positions in (ρ 1 , ρ 2 )-phase space of all possible fixed points in the model, together with color and name coding to identify and refer to them easily in the following sections and figures.Note the presence of the collective opinion states most relevant to our discussion:

Local dynamics and stability
A basic step in understanding the dynamics of the model is to explore the vector fields of Eqs.(9,10) and the associated trajectories close to the fixed points This can be done by means of a linearization of the dynamical equations.The linearization process leads to an exponential solution, where v 1 = (v 11 , v 12 ) and v 2 = (v 21 , v 22 ) are the eigenvectors with associated eigenvalues λ 1 and λ 2 of the Jacobian matrix J ij = ∂µ i ∂ρ j , evaluated at the fixed point ρ st 1 , ρ st 2 .C 1 and C 2 can be calculated from the initial condition as The stability of a fixed point is determined by the sign of (the real part of) the associated eigenvalues: for λ 1,2 < 0 it is stable, for λ 1,2 > 0 it is unstable, and for λ 1 < 0, λ 2 > 0 or λ 1 > 0, λ 2 < 0 it is a saddle point.Only when the fixed point is stable, we expect trajectories to converge to the fixed point in the long time limit [ρ . The eigenvalues can be calculated as where τ = J 11 + J 22 is the trace and ∆ = J 11 J 22 − J 12 J 21 the determinant of the Jacobian matrix evaluated at the fixed point ρ st 1 , ρ st 2 .If τ 2 > 4∆ the eigenvalues only have a real part (the case for all fixed points in Fig. 1).The condition for a transition (bifurcation) is that one of the eigenvalues becomes zero (a so-called marginal stability), or equivalently ∆[ρ st 1 , ρ st 2 ] = 0.This condition together with Eqs.(11,12) determines the transition lines and the phase diagrams.At a transition we expect some fixed points to appear or disappear (usually in couples).
In Fig. 2 we show the phase diagram and vectors fields in the prototypical scenario of pairwise interactions (low α) as a function of bias asymmetry ∆b.The phase diagrams in Fig. 2a and Fig. 2b correspond to the well-known cusp catastrophe [41], where transitions are saddle node bifurcations in which two fixed points (stable and saddle point or saddle point and unstable) merge and disappear for high enough bias asymmetry.In the case of two equal groups (p 1 = p 2 , Fig. 2a), when tuning bias asymmetry, polarization is destroyed favoring the consensus states.After that, for a specific high value of the asymmetry, one of the two consensus states disappears (Fig. 2c) and the only remaining state is ρ st 1 = ρ st 2 ≈ 0 for ∆b > 0, or ρ st 1 = ρ st 2 ≈ 1 for ∆b < 0. Every time a pair of fixed points merge, there is a region in phase space where the dynamics ρ 1 (t), ρ 2 (t) becomes very slow and meta-stable states appear (in Section 5 we discuss this in more detail).Note that in the symmetric version of the model (∆b = 0), when more than one stable fixed point exists, initial conditions determine the final state.However, when the exogenous ingredient of bias asymmetry is introduced by the social media (∆b = 0), it is possible to 'select and control' the final opinion of the system.
Another relevant scenario is that of asymmetric groups p 1 = p 2 , where one of them is either better connected and/or bigger in size, i.e. for p 1 < p 2 group 1 has more nodes or links than group 2, and the other way around for p 1 > p 2 .The two polarized states (ρ st 1 ≈ 0, ρ st 2 ≈ 1 and ρ st 1 ≈ 1, ρ st 2 ≈ 0) are not symmetric, depending on which is the opinion of the majority and minority groups.For this reason, there are two transition lines in the phase diagram of Fig. 2b with cusps at different positions, one for ∆b > 0 and the other for ∆b < 0. Thus, bias asymmetry promotes (instead of destroying) polarization in the region between the two cusps.This result has a clear interpretation: if the social platform favors the opinion of the minority group, polarization will become stronger as it will be harder for the majority group to convince the other, leading the system towards consensus.The value of ∆b at the cusp is the 'optimal' one if we wish to balance such majority-minority scenario (i.e., asymmetry in group sizes and connections) by using an exogenous algorithmic bias.
In Fig. 3 we plot the eigenvalues of all fixed points as a function of bias asymmetry in the same case specified in Fig. 2. The eigenvalues provide us with a lot of information about the nature of the fixed points and the dynamics close to them.The sign of the eigenvalues in Fig. 3 agree with the schematic representation of the vector fields in Fig. 2c, and it determines the stability analysis of the fixed points.Every time an eigenvalue becomes zero, a pair of fixed points disappears, defining a transition or bifurcation.In Fig. 4 we show the phase diagram and vector field in the case of group interactions (high α) as a function of bias asymmetry ∆b.Figs.4a and 4b correspond to the so-called butterfly catastrophe [41].We observe some differences with respect to the pair interaction case, besides the rich phenomenology of additional fixed points.The first difference is that the coexistence and consensus states can be both stable for some parameter values, at odds with results in Fig. 2. With respect to the dependence of the consensus states on bias asymmetry, for low noise we have a similar behavior as for pair interactions, while for high noise it is possible that a consensus state, which is not observed for ∆b = 0, appears for some value ∆b = 0. Standard polarization is also destroyed for a critical value of the asymmetry in the group interaction case.The difference is that new stable polarized states are possible in the group case (pol-coex 2 and pol-coex 4 in Fig. 1), whose behavior with respect to bias asymmetry is non-trivial.Similarly to the consensus states, these new polarized states may appear for a particular value of the asymmetry, even though they are not present in the symmetric case.
In Fig. 5 we show the eigenvalues of the fixed points of Fig. 1 and Fig. 4c as a function of ∆b for different values of the noise Q.These figures provide us with information about the stability, dynamics and transitions that are possible in this case.Note that, depending on the value of Q, we may find some of the fixed points or not, and that different transitions happen for specific values of ∆b, in accordance with the phase diagram in Fig. 4b.The vector field scheme in Fig. 4c is a low-noise scenario where all fixed points are present (Fig. 5c).For other values of the noise Q (in Fig. 5a and Fig. 5b), not all the fixed of Fig. 4c are possible, and the transitions may occur in different orders as we increase ∆b.

Global dynamics, convergence times and meta-stable states
The global dynamics ρ 1 (t), ρ 2 (t) of the system has a non-trivial dependence on the initial condition ρ 1 (0), ρ 2 (0) and the model parameters (Q, α, b, ∆b, p 1 , p 2 ), besides time t.The determination of the fixed points, the stability, and the local (linearized) dynamics are a good guideline to predict and understand, at least qualitatively, the dynamical behavior of the system.There are other aspects of the dynamics that cannot be fully explained by the fixed points and the linear dynamics approach.Among these, we study with particular attention what we call meta-stable states, where the dynamics slows down strongly and stays for a long time around a determined value of the state variables.This phenomenon is observed in the model, especially when two fixed points merge and disappear (the elliptical striped zones in Fig. 2c and Fig. 4c).As there are no fixed  points around these zones, it is not possible to evaluate the eigenvalues and explore the local dynamics.Thus, we need to use a different theoretical method to characterize the meta-stable states.We also explore the time needed to reach the final states and the dependence on the initial conditions, theoretically and by means of Monte Carlo simulations.

Numerical simulations
Before introducing the theoretical description of the meta-stable states, we analyze the results coming from Monte Carlo simulations.Implementing the rules of the model (Section 2), we obtain stochastic trajectories ρ 1 (t), ρ 2 (t), from which we calculate the average global state of the system ρ(t) and polarization P (t) that characterize the dynamics.In Fig. 6 and Fig. 7 we show numerical results for pair interactions (α = 1) on a large (N = 20000) z−regular type of network with modular structure.In order to compare the simulations with the phenomenology coming from the theory, we use the same parameter values as in Fig. 2a and Fig. 3, i.e.Q = 0.01, b = 0.8 and p 1 = p 2 = 0.1.In Fig. 6 we show the dynamics for various homogeneous initial conditions and bias asymmetries.In the top panels, from left to right, one of the consensus states becomes unstable for a determined value of the bias asymmetry and then, independently on the initial condition, all trajectories evolve towards the remaining stable consensus state.This corresponds to crossing (horizontally) the green transition line in the phase diagram of Fig. 2a.Note that before the transition, when the two consensus states are possible, and depending on the initial condition ρ 1 (0) = ρ 2 (0) = ρ(0), the dynamics evolves towards one state or the other.We can thus define a threshold initial condition ρ 0 that separates the basin of attraction of the consensus states.This threshold depends on the bias asymmetry ρ 0 (∆b): for no asymmetry (∆b = 0) it is ρ 0 = 0.5, it increases for ∆b > 0 (decreases for ∆b < 0), and it is not defined above the transition point as only one consensus state is stable.In Fig. 7, we show the dynamics for different polarized initial conditions and bias asymmetries.In the top panels, from left to right, the polarized state becomes unstable for a determined value of the bias asymmetry and then all trajectories evolve towards a non-polarized (P = 0) consensus state.This corresponds to crossing (horizontally) the blue transition line in the phase diagram of Fig. 2a.Note that the dynamical results shown in Fig. 6 and Fig. 7 are in good agreement with the qualitative description of the vector fields in Fig. 2c.
A significant dynamic phenomenon observed in Figs. 6 and 7 (specially in panels 6b, 7b, 7c and 7d) is the presence of trajectories that get trapped for a long period of time in some state, but eventually get released and end in one of the possible final (stable) states.In Fig. 7d we see that such a meta-stable state appears above the transition point ∆b * , and that its duration decreases as we increase the asymmetry.We represent this meta-stable state as an elliptical stripped zone in the vector fields Fig. 2c, and we characterize them theoretically in what follows (Section 5.2).

Saddle node bifurcations and meta-stability
Close to a critical (bifurcation) point, it is possible to obtain the normal form of the dynamics [42], which goes beyond the linearization of Eqs.(13,14).Assume that we have a fixed point ρ st 1 and ρ st 2 with an eigenvalue equal to zero λ 1 = 0, for a determined value of a tuning parameter, e.g.∆b = ∆b * (bifurcation or critical point).We perform a change of variables to the eigenvector basis u(t) = P −1 ρ(t), where the columns of the matrix P are the eigenvectors v 1,2 , that is, The evolution of the transformed variables u(t) becomes and the size of the network is N = 20000, divided in two groups (communities) of equal size (N 1 = N 2 = N/2), averaged over 1000 realizations of the dynamics.In the top panels (a-c), trajectories correspond to various initial conditions ρ(0) for a fixed value of the bias asymmetry ∆b (specified in the title).In the bottom panels (d-f), a color gradient is used for bias asymmetry ∆b in the range [0, 0.5] for a fixed initial condition ρ 1 (0), ρ 2 (0) (specified in the title).
The transformed vector field fulfills U (∆b * , u st ) = 0 (here u st refers only to the fixed point at the bifurcation point ∆b = ∆b * ) and ∂U i /∂u j = −λ i δ ij (again at ∆b * and u st ), i.e., the linear part is uncoupled.In this case, when λ 1 = 0, there is a center manifold u(t) = h(∆b, u 1 (t)), i.e. a special trajectory where the time dependence of all variables is governed by the slow u 1 (t).This satisfies u st = h(∆b * , u st 1 ) and the orthogonality condition ∂h i /∂u 1 = 0.The center manifold can be obtained from Eq. ( 19) as a series expansion.Once the functions h(∆b, u 1 ) have been determined, we obtain a single equation for u 1 (t), = β (10) (∆b − ∆b * ) + β (11) where β (10) , β (11) , β (02) , β (03) are coefficients whose expressions can be derived ‡ from the series expansion of the center manifold [43].‡ For example, for the polarization transition in Fig. 3 we obtain β (10) = 0.226376, β (02) = 0.672844, β (11) = 0.789372, and β (03) = −1.55356.We consider the most common bifurcation found as a function of the bias asymmetry ∆b (see Fig. 2c and Fig. 4c), i.e. the saddle node bifurcation with β (10) = 0 and β (02) = 0.In its simplified form we have where the higher order terms can be disregarded.Assuming positive coefficients β (10) > 0, β (02) > 0, and for ∆b < ∆b * , the solution of Eq. ( 21) is β (10) (∆b * − ∆b) β (02) tanh β (10) while for ∆b > ∆b * it is β (10) (∆b − ∆b * ) β (02) tan β (10) where C is determined from the initial condition u 1 (0).Note how, for ∆b < ∆b * in Eq. ( 22), the system goes to a stable fixed point in the infinite time limit as tanh(∞) = 1 and u 1 (∞) takes a finite value, while for ∆b > ∆b * in Eq. ( 23) the system slows down close to u 1 (t) ≈ u st 1 and then diverges (corresponding to a meta-stable state).According to Eq. ( 23), close to the bifurcation point ∆b ∆b * , the solution scales as u 1 (t) − u st 1 ∼ (∆b − ∆b * ) 1/2 and time as t ∼ (∆b From the previous derivation we infer that, as we increase bias asymmetry ∆b, we find a critical value ∆b * (see Fig. 2c and Fig. 4c) where two fixed points, a saddle and a stable point, merge and disappear at a saddle node bifurcation.At the critical point we observe a meta-stable (slow) region in both variables and time.The size of this region scales with the distance to the critical point as (∆b − ∆b * ) 1/2 in the variables and as (∆b − ∆b * ) −1/2 in time.The zone is centered at the position where the two fixed points merge, i.e. the fixed point at the critical point ∆b = ∆b * , and it is elongated along the slow eigendirection v 1 .

Convergence times
An important quantity to analyze the global dynamical behavior of the model is the convergence time needed to reach the final (stationary) state.Often, when modeling the opinion dynamics of a population, we are not only interested in the possible final states of the system, but also in this convergence time and its dependence on the parameters of the model.For reason of convenience, we write the time dependence of the global quantities as follows: such that g 1,2 (0) = 1 and g 1,2 (∞) = 0.In the linear regime of Eqs.(13,14), g 1,2 (t) are a sum of two exponential functions (e −λ 1 t and e −λ 2 t ) with amplitudes depending on the eigenvectors and initial conditions.The quantities λ −1 1 and λ −1 2 estimate well the two time scales involved in the time evolution close to the fixed point.For the general global dynamics g 1,2 (t), the exponential form is not valid and we have a complicated time-dependence.We measure the time scale t f of the global dynamics as the smallest of the solutions t f fulfilling the following conditions: where %V AL is an arbitrary percentage measuring how close the system is to the final state when the convergence time t f is reached.In Fig. 8 we show the dependence of the convergence time t f , extracted from numerical simulations, on bias asymmetry and initial conditions in the case of pair interactions, for the same parameters as in Figs. 6 and 7. Note that the system reaches different final states depending on the values of ∆b, ρ 1 (0) and ρ 2 (0).This can be clearly seen in Figs.8c and 8d, where a dark blue line separates two zones where the system reaches different final states, and where the behavior of the convergence times changes (in Figs.8a and 8b it corresponds to the minimum of the curves).The line that separates the two behaviors is what we call a threshold initial condition ρ 0 (∆b) in Section 5.1, i.e. the limit of the basin of attraction of the possible (consensus and polarization) final states.The dependence of the inverse of the convergence time (t −1 f ) with bias asymmetry ∆b, increasing or decreasing, shows a clear correspondence with the eigenvalues of the final (stable) state (polarization or consensus) in Fig. 3.When the final state is polarization (∆b < 0.02 in Figs.8b and 8d) it seems that the trend, increasing or decreasing, as a function of ∆b is not that clear.This can be also understood from Fig. 3, where there is one eigenvalue increasing and another decreasing with ∆b.Generally, the smallest eigenvalue dominates the dynamics, unless the initial condition is aligned with the fast eigendirection.That is the reason why we mainly observe a decreasing behavior in the simulations, while an increasing trend is also possible in some situations.
In Fig. 8a and 8b we can identify the characteristic scaling behavior t −1 f ∼ |∆b − ∆b * | 1/2 , of the saddle-node bifurcation.Before the transition point (∆b ∆b * ), this is directly related to the eigenvalue of the final state λ ∼ (∆b * − ∆b) 1/2 and can be considered as critical slowing down, while after the transition (∆b ∆b * ), the metastable state dominates the dynamics (see Fig. 7d) with an equivalent time scaling as in Eq. (23).Note that there are two saddle-node bifurcations, Figs.8b (consensus) and 8a (polarization), with different transition points ∆b * .

Summary and conclusions
In this paper we have studied the role of algorithmic bias and community structure in the potential rise of polarization of opinions in online social networks.We have devoted special attention to the temporal behavior of an archetypal two-state opinion-formation model, the language model, as well as to the role of the bias asymmetry ∆b, i.e. the possibility that the online platform favors one opinion over the other.We have derived a pair of mean-field differential equations for the relevant variables of the dynamics, the density ρ 1 (t) [ρ 2 (t)] of nodes in group 1 (2) holding opinion 1.This theoretical description accurately captures the phenomenology of the model and shows a good fit with numerical simulations.
The possible final opinion states reproduced by the model are: . All states are found in the whole spectrum between pair and group interactions displayed by the language model as the α parameter is changed.For some parameter values in the group interaction case, we also find additional polarized (polarization-coexistence) states with ρ 1 ≈ 0 and ρ 2 ≈ 1/2 (and the three equivalent states exchanging the groups 1 ↔ 2 and states 0 ↔ 1).Using linear stability analysis, we have determined the phase diagrams for pair and group interactions.In general, we find that sufficiently strong asymmetry in the bias is capable to destroy first the stability of the polarized states, and then one of the consensus states via saddle-node bifurcations.The phase diagram of the consensus states corresponds to the well-known cusp catastrophe for pair interactions and butterfly catastrophe for group interactions.Thus, bias asymmetry is a means to 'select' final states of the dynamics, controlling the global behavior of the system.
When the population is divided in two asymmetric groups in terms of size or connectivity (which can be thought of as a majority and minority groups), a somewhat different situation relating to the polarized states is produced by the model.In a range of values of the bias asymmetry (−∆b, ∆b), polarization is not necessarily suppressed but also favored, while above that range it is only destroyed.The two polarized states are not equivalent, depending on which is the opinion of the majority and minority groups.If the social platform benefits the opinion of the minority group then polarization is promoted, while in the opposite case polarization is suppressed in favor of consensus.The values of the bias asymmetry that delimit this behavior can be understood as the case where the structural asymmetry of the network (in group size and connectivity) is compensated by the 'favoritism' of algorithmic bias.
Results of numerical simulations confirm the possible final states predicted by the mean-field theory and the role of bias asymmetry.We have found that the convergence time (time to reach a stationary final state) depends on bias asymmetry, initial conditions, and the other parameters of the model in non-trivial ways.By means of the eigenvalues of the linearized dynamics, we were able to characterize this dependence close to the final states, gaining a better understanding of the type of transitions we have, and what happens to the dynamics when one of the final states loses its stability.An unavoidable phenomenon we observe when one of the stable solutions disappears is the presence of meta-stable states.This means that the system becomes trapped for a long period of time in a region of the phase space.Using the normal form of the bifurcations we have derived the scaling relations of the meta-stable zone as a function of the bias asymmetry.This shows a double square root law: ρ 1,2 (t) ∼ (t f ) −1 ∼ (∆b − ∆b * ) 1/2 .
In conclusion, we have explored a simplified opinion formation model including some of the arguably most relevant features driving opinion dynamics on online social networks: spreading processes (pair or group copying together with independent behavior of individuals), algorithmic bias, and an underlying networked community structure.The joint effect of these ingredients produces a complex phase space of collective social behavior including coexistence, consensus, and polarization of opinions.We used this formalism as testing ground to study the influence of algorithmic bias on online communication dynamics.We showed that bias imbalance can have a crucial effect on the final opinion state and the dynamics in general.Polarization and consensus states are destroyed for high enough bias asymmetry at different transitions, while just after the transition these destroyed final states become meta-stable.We characterized all possible transitions via phase diagrams.For the local dynamics (close to the final states) we used a linearization of the dynamical equations, while for the global dynamics the normal form of the bifurcation allowed us to detect meta-stable states.Finally, we calculated convergence times both from the theoretical description and by means of numerical simulations.
The aim of such simplified modelling cannot be the quantitative reproduction of some empirical observations, e.g., related to elections.Still, we think that the richness of the behavior of the model, the relatively large number of relevant parameters, and the non-triviality of the temporal evolution of opinions are all features with strong relation to real-world systems.For example, the effect of the parameters on polarization is sometimes counter-intuitive: How algorithmic bias influences the outcome of the dynamics depends strongly on the type of interaction.In our model we can tune interactions from pairwise to group -in reality both are present and highly heterogeneous.Another relevant observation is the frequent appearance of meta-stable states.As real-world phenomena evolve over finite times (e.g., opinions matter on election day), meta-stable states may be crucial in forecasting efforts.In the future we would like to combine modeling with empirical data analysis, partly from observational data and partly from controlled social experiments, in order to further understand the interplay between human collective action and online algorithms.
these states are possible in the low (pair) and high (group) α regimes for certain values of the model parameters (Q, α, b, ∆b, p 1 , p 2 ).

Figure 1 .
Figure 1.Schematic representation of fixed points of opinion dynamics in (ρ 1 , ρ 2 )phase space.The left (right) side corresponds to the low (high) α regime of the model.Joined colored circles inside a square represent the opinion states s = 0 (blue) or s = 1 (red) of each group (e.g., two circles of the same color represent consensus, while circles of different colors represent polarization).The legend displays the color coding and names of the most relevant fixed points.Pairs of fixed points [circle (stable) and square (unstable)] have the same color if they disappear/appear together for specific values of the parameters.

Figure 2 .Figure 3 .
Figure 2. Phase diagrams for (a) p 1 = p 2 = p = 0.1, and (b) p 1 = 0.05, p 2 = 0.1, and vector fields (c) for fixed values of the model parameters α = 1 (linear or pairwise regime) and b = 0.8.In the phase diagrams (a,b) the varying parameters are (Q, ∆b), i.e. the noise and bias asymmetry.The transition lines (green and blue) delimit the parameter regions where the different possible fixed points are stable.Circles inside a square indicate the corresponding fixed point following the scheme of Fig. 1.The phase diagram (a) corresponds to the case of two equal groups, while in (b) there is a majority (big circle) and minority (small circle) group.In the region below the green line both consensus states are stable, and above only one of them remains, while below the blue line polarization is stable.In panel (c), the left vector field is a typical situation below the blue line of the phase diagram (a) for ∆b = 0.The other vector fields (from left to right) show how this changes as we increase ∆b and cross the various transition lines.The elliptical striped zones are regions where the dynamics is very slow and meta-stable states are possible (see Section 5).

Figure 4 .
Figure 4. Phase diagrams for (a) α = 4, b = 0.8, and (b) α = 6, b = 0.75, and vector fields (c) for fixed connectivity p 1 = p 2 = p = 0.1.In the phase diagrams (a,b) the varying parameters are (Q, ∆b), i.e., the noise and bias asymmetry.The transition lines (green, red, dark and light blue, and yellow) delimit the parameter regions where fixed points are stable.Circles inside a square indicate the corresponding fixed points following the scheme of Fig.1.In panel (c), the left vector field is a typical situation below the blue line (and above the red, light blue and yellow lines) of the phase diagram (b) for ∆b = 0.The other vector fields from left to right show how this changes as we increase ∆b and cross the various transition lines.The elliptical striped zones are regions where the dynamics is very slow and meta-stable states are possible (see Section 5).

Figure 5 .
Figure 5. Eigenvalues of the various fixed points of Fig.1and Fig.4cfor parameter values α = 6 (group regime), b = 0.75, p 1 = p 2 = p = 0.1 and Q = 0.01, 0.005, 0.002, as a function of bias asymmetry ∆b.The color and name coding in the legends is equivalent to that of Fig.1.Lines of the same color (light and dark blue, red and yellow) correspond to a pair of fixed points, dashed (saddle point or unstable) and solid (stable), that merge together and disappear for a particular value of ∆b.

Figure 7 .
Figure 7. Average polarizationP = |ρ 1 − ρ 2 | , starting from a polarized initial condition ρ 2 (0) = 1 − ρ 1 (0).All model, connectivity and network parameters are the same as in Fig.6.In the top panels (a-c) a fixed bias asymmetry is used for different initial conditions, while in the bottom panels (d-f) a fixed initial condition is used for different bias asymmetries (colors).