A Calculation of the coupling integrals

Tricritical points in a Vicsek model of self-propelled particles with bounded confidence

Abstract

We study the orientational ordering in systems of self-propelled particles with selective interactions. To introduce the selectivity we augment the standard Vicsek model with a bounded-confidence collision rule: a given particle only aligns to neighbors who have directions quite similar to its own. Neighbors whose directions deviate more than a fixed restriction angle are ignored. The collective dynamics of this systems is studied by agent-based simulations and kinetic mean field theory. We demonstrate that the reduction of the restriction angle leads to a critical noise amplitude decreasing monotonically with that angle, turning into a power law with exponent for small angles. Moreover, for small system sizes we show that upon decreasing the restriction angle, the kind of the transition to polar collective motion changes from continuous to discontinuous. Thus, an apparent tricritical point with different scaling laws is identified and calculated analytically. We investigate the shifting and vanishing of this point due to the formation of density bands as the system size is increased. Agent-based simulations in small systems with large particle velocities show excellent agreement with the kinetic theory predictions. We also find that at very small interaction angles the polar ordered phase becomes unstable with respect to the apolar phase. We derive analytical expressions for the dependence of the threshold noise on the restriction angle. We show that the mean-field kinetic theory also permits stationary nematic states below a restriction angle of . We calculate the critical noise, at which the disordered state bifurcates to a nematic state, and find that it is always smaller than the threshold noise for the transition from disorder to polar order. The disordered-nematic transition features two tricritical points: At low and high restriction angle the transition is discontinuous but continuous at intermediate . We generalize our results to systems that show fragmentation into more than two groups and obtain scaling laws for the transition lines and the corresponding tricritical points. A novel numerical method to evaluate the nonlinear Fredholm integral equation for the stationary distribution function is also presented. This method is shown to give excellent agreement with agent-based simulations, even in strongly ordered systems at noise values close to zero.

pacs:
05.65.+b, 64.70.qj, 87.18.Nq

I Introduction

Dynamic self-organisation and, in particular, mechanisms of flocking behavior in groups of living species remain one of the most intriguing problems at the interface of physics and biology. Numerous physical models of interacting self-propelled agents have been proposed recently to study these phenomena (see review papers Toner et al. (2005); Ramaswamy (2010); Vicsek and Zafeiris (2012)). The collective behavior much resembling the dynamics of living organisms has also been observed in a variety of synthetic systems, generally referred to as active soft matter Paxton et al. (2004); Howse et al. (2007); Golestanian et al. (2007); Thakur and Kapral (2011); Taktikos et al. (2012); Theurkauff et al. (2012); Marchetti et al. (2013). The level of consciousness of the individuals apparently plays minor role for the large scale dynamics as the same general principles that apply to groups of animals or cells seem to govern also human social phenomena, traffic, robotics, and decision making Helbing et al. (2000, 2001); Deffuant et al. (2001); Weidlich (2002). Therefore, similar modelling approaches are employed to describe their generic dynamic properties.

One of the simplest and earliest models to describe the collective motion of self-propelled agents has been introduced by Vicsek in 1995 Vicsek et al. (1995). In this original paper, a second order phase transition between the orientationally ordered, globally aligned motion state and disordered state was claimed. The continuous nature of the transition between these states was also supported by a number of publications from the same group involving original Vicsek model (VM) with angular noise Czirók et al. (1997); Czirók and Vicsek (2000). However, the nature of the transition has been questioned in a number of later studies Grégoire and H (2004); Chaté et al. (2008a). Chaté et al. have demonstrated Chaté et al. (2008a) that there exists a critical system size , beyond which a discontinuous, or first order transition is observed. The dependence of this critical system size on particle density has been theoretically estimated in Ref. Ihle (2011). It has also been reported that the kind of the transition seems to depend on particle velocities Baglietto and Albano (2009) and on the way, in which the noise is introduced into the system Aldana et al. (2007). All above mentioned factors can result in a behavior, which can be associated with instabilities leading to co-existence of the ordered and disordered phases at the transition point, or, in other words, to a first order phase transition Nagy et al. (2007). The discontinuous character of the transition has also been elucidated from analysis of Ising-like 1D models of flocking Solon and Tailleur (2013).

The contribution of range and symmetry of aligning interactions to the instabilities remains relatively poorly understood. Recent numerical studies demonstrated essential differences between the symmetry of the ordered phase and stability of density bands depending on whether the active agents have polar or apolar alignment mechanism Peruani et al. (2011). In this paper we study the role of selectivity of interactions in the ordering of the Vicsek model (VM). This property has a number of important applications. For example, the selectivity of the interaction has previously been introduced in the study of the social models such as the bounded confidence model Deffuant et al. (2001). In a situation where the agents are prepared to align themselves only with the fellow individuals with the opinion vector not too different from their own one, the restrictive rules can become crucial for the polarization of the group and the collective decision making. The bounded confidence rules can be introduced into the Vicsek model as a restriction on the angle of interaction, which makes the alignment of the individual with too different directions of motion impossible Deffuant et al. (2001); Ben-Naim et al. (2003). Given that the Vicsek model was originally introduced to model the behavior of social agents such as birds and fish, such a selectivity rule could make their description more realistic.

Another motivation to study an interaction rule that is sensitive to orientation differences are the very recent experiments by Lu et al. Lu et al. (2013) on the collective behavior of Bacillus subtilis in the presence of a photosensitizer. To account for the cell-to-cell interaction via intercellular flagella bundling, the authors propose a Vicsek-like model where the interaction among neighbor agents becomes weaker with increasing orientation difference. Since our model explores an extreme version of this interaction, the kinetic theory from this paper might become useful for a better analytical understanding of these experiments.

In this paper, we derive the ordering behavior in the Vicsek model with bounded confidence from the microscopic kinetic equations and demonstrate how a qualitatively new behavior arises from the interaction selectivity. In particular, we find within homogeneous mean-field theory that at strong selectivity, the transition from a disordered state to a state of polar order is discontinuous. If one increases the restriction angle beyond a so-called tricritical point the transition becomes continuous. We also analyze nematic phases and fragmented states which consist of several aligned groups moving in different directions. A scaling law for the transition from disorder to an ordered state with fragments is calculated and two different tricritical points are identified. Furthermore, we present a novel numerical method to accurately evaluate the nonlinear Fredholm integral equation for stationary distribution functions. The paper is organized as follows: in Section II we introduce the kinetic theory of the Vicsek model and predict its phase behavior, describe the numerical solution methodology for the kinetic equations and the settings of the agent-based simulation. In Section III, we present the numerical and analytical results for the model: phase diagrams and data on the ordering behavior. Finally, in Sections IV and V we discuss the significance of main findings of the work and summarize the results. Details about the calculations of coupling integrals are relegated to Appendix A. In Appendix B we provide a general discussion of terms such as “phase transition” that are borrowed from equilibrium statistical mechanics but here are applied to finite nonequilibrium systems.

Ii Theory and simulation settings

ii.1 Model

We consider a two-dimensional model with point particles at number density , which move at constant speed . The particles with positions and velocities undergo discrete-time dynamics with time step . The evolution consists of two steps: streaming and collision. In the streaming step all positions are updated according to

(1)

In the subsequent collision step, the directions of the velocity vectors change. Similar to the standard VM, particles align with their neighbors within a fixed distance . However, the interaction in the present paper is selective such that the particles align only with those neighbors whose direction of motion deviates by an angle less than some fixed value from their own velocity vector (see Fig. 1). In this implementation, the Vicsek model becomes similar to so-called “Bounded confidence” models, commonly used in social sciences to study opinion dynamics Deffuant et al. (2001); Ben-Naim et al. (2003). It simulates the common social tendency to disregard opinions that appear too extreme with regard to their own perspective. Thus, the parameter can be interpreted as the degree of ignorance of a population of self-propelled agents. In particular, a circle of radius is drawn around a given particle and the average angle of motion of the particles within the circle is determined according to

(2)

where particles whose inner product is smaller than are excluded from the summation. In an extreme case, it is possible that even if particle has many neighbors in its collision circle, all are rejected due to too large differences, and thus . The regular VM is recovered in the limit . Once the average angles are known, the new particle directions follow as

(3)

where is a random number which is uniformly distributed in the interval . Note, that the updated positions (and not the old locations ) are used to determine the average angles of motion . This corresponds to the so-called forward updating rule of the standard VM, as defined in Refs. Huepe and Aldana (2008); Baglietto and Albano (2009)).

Note, that our model is qualitatively different from models with a finite sensing region such as the “blind-spot” models introduced by Couzin et al. Couzin et al. (2002), and used, for example, in Refs. Lukeman et al. (2010); Wood and Ackland (2007). In contrast to these models, particle exclusion is based on velocity differences and not spatial location. It can be easily shown that Vicsek-like models with rear blind sectors do not show tricritical points in homogeneous mean-field theory. However, the version of the VM presented here is related to the one introduced by Lu et al. Lu et al. (2013), where instead of completely excluding misaligned partners, the interactions are weighted by the orientation difference.

Figure 1: The interaction parameters in the restricted angle Vicsek model. The particle aligns itself only with those neighbors (blue arrows) whose relative angle of motion is less or equal to .

ii.2 Kinetic theory

Deriving the kinetic equation

Recently, a kinetic theory formalism for self-propelled particles has been introduced that can handle discrete time dynamics and “exotic” collision rules such as genuine multi-body and topological interactions Ihle (2011); Chou et al. (2012); Ihle (2014a). It has been shown Ihle (2013) that this approach is able to quantitatively reproduce the results of agent-based simulations, in the limit of large mean free path, , compared to the radius of interaction . In this chapter we will adapt this method to the restricted-angle model. The formalism starts with an evolution equation for a Markov chain in phase space,

(4)

where is the -particle probability density referring to an ensemble of independent Vicsek systems, which are initialized at time with some initial probability density . This initial density is assumed to be symmetric against permuting particle indices. Equation (4) describes the transition from a microscopic state to the state during one time step with transition probability . The microscopic state of the system at time is given by the -dimensional vector, , where contains the directions of motion of all particles, and describes all particle positions. The initial microscopic state at time is denoted as . The integral over the initial state translates to . Pre-collisional angles and positions are given by and , respectively. The transition probability contains the microscopic collision rules,

(5)

Here, , is the periodically continued delta function, which accounts for the -periodicity of angles. The particle velocities , are given in terms of angle variables, .

The kinetic equation for the -particle probability density, Eq. (4) is exact but intractable without simplification. Here, as done, for example, in Ihle (2011), we use Boltzmann’s molecular chaos approximation by assuming that the particles are uncorrelated prior to a collision, which amounts to a factorization of the -particle probability into a product of one-particle probabilities, . This approximation can be justified at moderate and large noise strength and when the mean free path (mfp) is large compared to the radius of interaction . Here, the mfp is defined as the distance a particle travels between collisions, , and is density-independent due to the discrete nature of the dynamics. More details on the validity of molecular chaos and a general discussion of kinetic theory approaches to Vicsek-style models can be found in Refs. Ihle (2014a, b); Peshkov et al. (2014); Chou and Ihle (2014).

Because molecular chaos neglects pre-collisional correlations, it leads to a mean-field theory. To derive this mean-field theory, we follow Refs. Ihle (2011); Malevanets and Kapral (1999); Ihle (2009) and multiply Eq. (4) by the phase space density . A subsequent integration over all particle positions and angles leads, in the large limit, to an Enskog-like kinetic equation for the one-particle distribution function, ,

(6)

where is the average number of particles in a circle of radius centered around and can be position-dependent. The subscript “” at the integral denotes integration over this circle. The local particle density is given as a moment of the distribution function, ; denotes the integration over all positions, particles can assume within the interaction circle; is the average over all pre-collisional angles of particles in the interaction circle. In Eq. (6) particle is assumed to be the focal particle. It is fixed at position and particles are supposed to be its neighbors. More details on how to interpret equations similar to Eq. (6), can be found in Ref. Ihle (2014a).

Solving the kinetic equation

We restrict ourselves to spatially homogeneous solutions of the Enskog-like equation. Then Eq. (6) becomes

(7)

with the simpler collision integral

(8)

where is the area of the collision circle and the average particle number in this circle, , is proportional to the constant particle number density . For stationary solutions where , Eq. (7) constitutes a nonlinear Fredholm integral equation of the second kind with a singular kernel. Further below, we will present a novel numerical method to solve this equation with high precision. In principle, spatially inhomogeneous solutions, including steep soliton-like waves, could be obtained by the Lattice-Boltzmann-like method of Ref. Ihle (2013) but this is beyond the scope of this paper. Instead, we will handle inhomogeneous states by agent-based simulations only.

A convenient starting point for analytical and numerical studies of Eq. (7) is the angular Fourier expansion of both the distribution function and the collision integral,

(9)

The emergence of a globally ordered state breaks the rotational symmetry. Particles start to flow preferentially in a common but arbitrary direction . Since we are only interested in steady states of homogeneous systems in the thermodynamic limit, it suffices to set and to only keep the cosine terms in the Fourier expansion FOO (a). The zeroth mode is proportional to the average density,

(10)

because . The first mode, , serves as polar order parameter and is determined by the average -component of the momentum density ,

(11)

with . By definition, we have .

Formally, Eqs. (6), (7) look identical to the kinetic equation derived previously in Refs. Ihle (2011, 2014a) for the regular VM. The difference hides in the definition of the average angle of the focal particle: In the regular VM all particles within interaction range are accepted, irrespective of their orientation. This leads to a straightforward evaluation of the integrals over the pre-collisional angles . For the restricted angle model, the integration domain has to be split into subdomains because the number of angular arguments in depends on the values of the pre-collisional angles.

For simplicity we assume low densities, , and only include self-interactions and binary collisions. This amounts to truncating the sum in Eq. (6) after . As pointed out in the supplemental material of Ref. Ihle (2013), in order to enforce exact mass conservation, any such truncation must be accompanied by a consistent rescaling of the Poissonian weight factor. In our case, must be replaced by . The self-interaction or self-diffusion term with is the same as in the regular VM. However, in the evaluation of the binary collision term, two cases have to be distinguished, (i) the direction of particle deviates too strongly from that of particle , that is , and (ii) the “opinion” of particle is accepted by particle . In the first case, whereas for case (ii) which can be reformulated by means of trigonometric identities as

(12)

Multiplying Eq. (6) by and integrating over will lead to an infinite set of algebraic equations for the Fourier coefficients ,

(13)

where

(14)

The coupling constants in Eq. (13) are given by angular integrals,

(15)

where the binary collision couplings depend on the angle and have been split into two types,

(16)

Here, corresponds to the situation where the focal particle rejects its neighbor’s “opinion”, whereas in the directions of both particles and contribute to the average direction . The first set of integrals in Eq. (16) is easily evaluated. For one finds,

(17)

where is Kronecker’s delta. For we have,

(18)

In order to transform into simple trigonometric integrals, only the first identity of Eq. (12) is needed because the integrals for are set up such that and by definition . This gives

(19)

The discussion of this integral is relegated to Appendix A. One major difference of the hierarchy equations, Eq. (13) to those of the regular VM is, that the coupling matrices are asymmetric with respect to interchanging and . As discussed in Appendix A, this is due to the social bias of an agent to favor its own “opinion”.

To verify the consistency of Eq. (13), we evaluate the first hierarchy equation for . Mass conservation requires that the homogeneous density stays invariant in every iteration, . Therefore, must be equal to . From Eq. (15) we see that and that these are the only nonzero coefficients entering the equation for . Substituting these coefficients into Eq. (13) leads to

(20)

Using Eq. (10) and the definition of we obtain the correct result, .

At any noise the hierarchy equations, Eq. (13), have a trivial solution where all coefficients for . This solution describes the disordered state; there is no preference for any direction. Similar to the regular VM, for low noise we expect an ordered solution with , that bifurcates off the disordered solution at the threshold noise . On this branch of the solution, at a noise slightly below , the first mode dominates all higher modes, that is the ratios go to zero if . Hence, to find the branching point we assume stationarity, , and neglect all modes for in Eq. (13). Then, the second hierarchy equation yields,

(21)

and with the help of Eqs. (10), (15), (17) a transcendental equation for follows:

(22)

where is the amplification factor for the mode , which is proportional to the -component of the momentum density. In the limit of no restriction, , we find and recover the threshold equation for the regular VM (Eq. (12) in Ihle (2014a) where terms with are truncated). Expanding Eq. (22) in the low density limit, , leads to

(23)

For , this expression agrees with the small density expansion for the regular VM, Eq. (13) in Ihle (2014a),

(24)

Furthermore, the function is non-negative and increases monotonically with for , as anticipated. Investigating the additional limit of strong restriction, Eq. (22) gives

(25)

Eq. (II.2.2) predicts that goes to zero at infinite “ignorance” where all particles perform independent random walks and never align with anybody. Thus, as one would intuitively guess, in this case the theory claims that no ordered state exists. Fig. 2 shows the numerical solution of Eq. (22) for two different normalized densities, and and compares it with the asymptotic formulas, Eqs. (23) and (II.2.2). One sees that for small , the low- expansion agrees quite well with the exact result for all angles . The asymptotic power law for is superlinear, . However, due to the change of curvature of the -curve from positive to negative, at intermediate angles it appears as if there is linear scaling, . As discussed further below, this is what we observed in agent-based simulation which were not performed for very small angles, , due to numerical limitations.

Figure 2: Prediction of the mean-field theory for the critical noise amplitude at different densities (a) . (b) . The restriction angle is given in units of .

ii.3 Analytical solution and mean-field tricritical point for polar order

Near the flocking threshold, higher order Fourier modes are suppressed, and Eq. (7) can be straightforwardly solved by setting higher modes to zero, truncating the infinite hierarchy (13) after the first few equations. The Fourier coefficients are normalized by means of the density ,

(26)

which amounts to the choice . Then, the first three nontrivial hierarchy equations from Eq. (13) become,

(27)
(28)

Here, we kept only the modes to and neglected all others because we are mainly interested in understanding the nature of the bifurcation to a (homogeneous) ordered state. The coupling constants are obtained by symmetrization of the coefficients defined in Eq. (15),

(29)

The depend on the restriction angle and are given in Appendix A. Starting from the third line, equations (27) can be solved successively leading to expressions of the higher modes in terms of the first mode ,

(30)

with the abbreviations

(31)

Substituting Eqs. (30) into the first line of Eq. (27) yields a closed expression for the first mode ,

(32)

with defined in Eq. (22) and

(33)
(34)

where at the threshold one has (see Eq. (22)), and . The character of the bifurcation to the ordered state, that is to non-zero , depends on the sign of the coefficient , evaluated at the threshold. Therefore, the condition for a tricritical point, where the character changes from subcritical to supercritical is

(35)

This is only possible if at least one of the following quantities vanish, , , or . In the low density limit that our approach is based on, the threshold noise is small, , and thus and are of order one and cannot be zero. Furthermore, it is easy to see from Eq. (69) that can only be zero in the limit where no ordered state exists. However, the equation has a solution at the angle . Thus, assuming spatially homogeneous states, we found an apparent tricritical point: for all angles smaller than the flocking transition is discontinuous, in the mean-field limit of large mean free path. For angles larger than and in a small system, the transition is predicted to be continuous. However, one has to keep in mind that in the regular VM, the homogeneous flocking state has a long-wave instability right next to the flocking threshold Bertin et al. (2009); Ihle (2011). This leads to inhomogeneous, soliton-like states that change the order of the phase transition to discontinuous in systems larger than a critical size Ihle (2013), even if the transition of a homogeneous system is predicted to be continuous. Something similar is expected in our model and will be investigated further below. For a justification of the term “phase transition” in finite systems, see Appendix B.

Note that the inhomogeneous reorganization of the system due to emerging soliton-like waves at large system sizes is qualitatively different from the well-known finite-size effects of equilibrium systems in the (grand)canonical ensemble. In these systems, the trivial fact that the correlation length cannot exceed the system size is used to set up finite-size scaling and to extract the behavior at infinite system size. As shown in Ref. Baglietto and Albano (2008) it is possible to perform such a finite-size scaling analysis of the VM at system sizes below and to obtain critical exponents as well as consistent hyper-scaling relations for a second-order phase transition. However, unless is infinite for the particular model and the parameters used, such an analysis would not extract the correct behavior in the thermodynamic limit because it would miss possible density instabilities.

It is interesting to see that (i) the tricritical angle does not depend on density, at least in the low density limit considered here, and (ii) its mathematical cause is the vanishing of the coupling between the modes and in the hierarchy equation for . Analyzing Eq. (32) at the tricritical point and leads to the mean-field exponent of for the order parameter scaling,

(36)

The structure of the denominator of the coefficients and also provides an estimate on the validity of the three-mode expansion. If the distance to the threshold noise is increased, the order parameter grows. Once it is so large that goes to zero, the approximation is expected to break down violently. This already happens not too far from the threshold and can be seen in Fig. 3, where the blue dashed curve describes the analytical solution of Eq. (27). This approximation neglects all modes with larger than three. Near the threshold, as anticipated, it agrees perfectly with the numerical solution of the Fredholm equation which is explained in a later chapter. However, deeper in the ordered phase, at smaller noise, the 3-mode-approximation starts to show unphysical behavior.

Figure 3: Comparison of the scaling of the order parameter with noise obtained analytically using Eq. (27) (all modes higher than are set to zero) and with the high- approximation method described in section II.6. Lines only serve as guides to the eye. Here, , and .

The physical reason for this deviation is the negligence of the higher modes , which in reality would start to grow when increases. In fact, using the hierarchy equations (13) it can be shown that for and all modes for become equal to the same value . This is a simple consequence of the fact that at zero noise all particles take on the same orientation and the distribution function becomes equal to the periodically-continued delta-function, .

ii.4 Nematic solutions and tricritical points

Figure 4: The critical noises for the transition from disorder to polar order (black curves), , and to nematic order (red curves), , are plotted as a function of restriction angle . The curves were calculated by means of Eqs. (22) and (38), respectively. The circles denote tricritical points at which the transition changes from discontinuous (dashed lines) to continuous (solid lines).

As described in section III, in agent-based simulations we sometimes encounter groups of particles which move in opposite directions. These states are characterized by a large nematic order parameter but only a small polar order parameter . In social science’s models of opinion dynamics, this is called polarization Hegselmann and Krause (2002); Kurmyshev et al. (2011); Baronchelli et al. (2007). To analytically explore such a possibility we reanalyze the infinite hierarchy for the Fourier coefficients, Eq. (13) with respect to nematic order. Compared to a state of polar order, a perfect nematic state has the additional symmetry of . This requires all odd Fourier coefficients in the series, Eq. (9), and therefore also the polar order parameter , to vanish exactly. Setting in Eq. (13), one realizes that the odd and even coefficients are decoupled: while all odd coefficients , become zero, the even coefficients can be non-zero and depend only on . This is because the coefficients and always vanish if neither nor is true. Setting and normalizing the coefficients as prescribed by Eq. (26), we obtain a new hierarchy for stationary nematic states,

(37)

where modes and higher are neglected. The coupling coefficients are given in Appendix A. Similar to the procedure for finding polar order, we check whether a nematic state can bifurcate from a disordered state at some critical noise value . Near such a (so far hypothetical) bifurcation, all modes higher than are negligible and the first equation of the hierarchy, Eq. (37), gives the following consistency condition,

(38)

where is the amplification factor of the mode . Interestingly, while this transcendental equation has no nontrivial solutions for the regular Vicsek model at any density, it does have a solution if the restriction value is smaller than a cut-off angle . Fig. 4 shows a plot of the nematic threshold noise as a function of together with the polar threshold noise at density . The cut-off does not depend on particle density, at least in the low density approximation applied here. Fig. 4 tells us that at small enough and , stationary nematic states do exist. We also see that at small , the threshold noise for nematic states is only slightly below the threshold for polar states, whereas the relative difference goes to one if the cut-off is approached from below. While our analytical discovery of stationary nematic states does not guarantee that these states are stable and relevant for the long-time behavior of agent-based simulations, our results are consistent with the apparent longevity of states with large nematic but small polar order in microscopic simulations. Of course, to fully explore the competition of polar, nematic and other apolarly ordered states, a comprehensive stability analysis similar to Refs. Chou et al. (2012); Menzel (2012), and more simulations are needed. This is beyond the scope of this paper. However, similar to the case of polar order, it is straightforward to analyze whether the bifurcation from a disordered to a (homogeneous) nematic state is continuous or discontinuous. For this purpose, we successively solve the hierarchy equations, Eq. (37), by first expressing the higher modes and in terms of :

(39)
(40)

with the abbreviations

(41)

Inserting these expressions into the first hierarchy equation, Eq. (37), leads to a closed-form expression for ,

(42)

with

(43)
(44)

One has , and similar to the polar case, see Eq. (35), the zeros of the quadratic coefficient define tricritical points,

(45)

where the bifurcation changes from subcritical to supercritical or vice versa. Surprisingly, below the cut-off angle we find two tricritical points, and . Further analyzing and at the critical line , we find that and for both and . Thus, the transition would be discontinuous in these two regions. Note, that is always negative for . In Fig. 4, discontinuous bifurcations are described by dashed lines. Inbetween the two points, e.g. at , is positive and, depending on , can be either positive or negative. This means that the transition to a homogeneous nematic state from the disordered state is continuous in this middle section, denoted by a solid red line in Fig. 4.

ii.5 Fragmentation into ordered groups

Figure 5: The scaled critical noises for mean-field transitions from the disordered state to symmetric states with fragments versus the scaled angle for a variety of densities and fragment numbers . The filled circles represent tricritical points where the character of the transition changes from discontinuous to continuous or vice versa.

Models of bounded confidence in social science can have states that consist of several clusters of opinions in which agents achieve local consensus Hegselmann and Krause (2002); Kurmyshev et al. (2011); Baronchelli et al. (2007). This is called opinion fragmentation. Translated into the language of our system, fragmentation would correspond to several groups of particles whose members interact with each other in every iteration because their angular differences are smaller than , but there is hardly any interaction between members of different groups. The state of polar order would then correspond to global consensus where only one group exists. The nematic state would be an example for fragmentation into two groups. A natural question to ask is whether our model also allows states with three, four and more distinct groups. It turns out that the analysis of the hierarchy equation, Eq. (13), for polar and nematic order can be straightforwardly generalized to fragmentation into groups. In this paper, we will restrict ourselves to highly symmetrical arrangements where the mean angles of the participating groups point into the directions with . In these arrangements, the distribution function has the “mirror”-symmetry . This restriction allows us to still use the Fourier cosine expansion, Eq. (9). The more general case requires the inclusion of sine terms and will be left for the future. Symmetric fragmentation into groups requires the vanishing of all coefficients whose mode numbers are not multiples of due to the -fold symmetry, . For example, at , three groups move into the three main directions and , and are described by the Fourier coefficients . Similar to the nematic case with , these coefficients only couple to themselves, e.g. they do not generate coefficients such as or which are supposed to remain zero. All expressions and coefficients for the nematic state can be easily generalized by formal replacements such as , , , , and so on. In particular, the threshold noise for the bifurcation of the disordered state to a state with symmetric fragments follows from the condition that the amplification rate for the mode is equal to one,

(46)

where

(47)

Setting in Eq. (46), the cut-off selectivity angle above which stationary fragmented states are impossible, follows from the transcendental equation

(48)

with . The solution is , and therefore the cut-off for a fragmented state of main directions is . The scaling has a simple physical interpretation: The difference between the mean angles of adjacent groups is equal to . As long as is smaller than this difference , members of the two groups have only a negligible chance of interaction. It is plausible to assume that the maximum possible restriction angle scales with the available angular range , , which is indeed what we find analytically. The ratio means that if is smaller than of the “opinion” difference , fragmentation becomes possible. In this notation, the critical noises for the transition from disorder to polar or nematic order are just special cases with and .

Analyzing Eqs. (46) and (47), we find that the critical noises for all possible symmetric fragmentations with follow an approximate scaling law

(49)

with the universal scaling function . According to relation (49), plotting the scaled critical noise