Scale-Free Chaos in the 2D Harmonically Confined Vicsek Model

Animal motion and flocking are ubiquitous nonequilibrium phenomena that are often studied within active matter. In examples such as insect swarms, macroscopic quantities exhibit power laws with measurable critical exponents and ideas from phase transitions and statistical mechanics have been explored to explain them. The widely used Vicsek model with periodic boundary conditions has an ordering phase transition but the corresponding homogeneous ordered or disordered phases are different from observations of natural swarms. If a harmonic potential (instead of a periodic box) is used to confine particles, then the numerical simulations of the Vicsek model display periodic, quasiperiodic, and chaotic attractors. The latter are scale-free on critical curves that produce power laws and critical exponents. Here, we investigate the scale-free chaos phase transition in two space dimensions. We show that the shape of the chaotic swarm on the critical curve reflects the split between the core and the vapor of insects observed in midge swarms and that the dynamic correlation function collapses only for a finite interval of small scaled times. We explain the algorithms used to calculate the largest Lyapunov exponents, the static and dynamic critical exponents, and compare them to those of the three-dimensional model.


I. INTRODUCTION
Experiments on the social behavior of animals in the laboratory may give rise to unreproducible results due to imposing artificial tasks on the animals and subsets of animals behaving differently [1].Experiments in natural environments may allow for observing the emergence of social behavior free from artificial laboratory constraints.Examples include flocks of birds [2] and sheep [3], fish schools [4], marching locusts [5], swarms of midges [6] and hordes of rodents [7,8].In these systems, collective behavior results from the dynamical interaction between individuals often producing power laws, which poses the question whether biological systems are in the critical region of a phase transition [9,10].
Many works have tried to apply ideas from phase transitions in statistical mechanics (scale-free behavior, finite size scaling, renormalization group, critical exponents, universality [11][12][13][14]) to the emergent collective behavior in animals and to compare them to experiments [15].Many theoretical studies have dealt with the Vicsek model (VM) of particles moving at a constant speed in a box with periodic boundary conditions and changing their velocity at discrete times by selecting the mean average velocity of all particles in their neighborhood plus an alignment noise [16][17][18].This model presents an ordering phase transition when the alignment noise falls below a critical value [18], which reminds of an equilibrium second-order phase transition between disordered and ordered homogeneous phases at the critical temperature [11].For animal collectivities, the periodic boundary conditions are artificial and it has been argued that placing the particles in a confining harmonic potential keeps the cohesion of the flock [19][20][21][22].
In recent works, we have numerically simulated the three-dimensional (3D) harmonically confined Vicsek model (HCVM) and discovered a phase transition between phases within chaos [23].
We have also analyzed the mean field HCVM near the scale-free-chaos phase transition and found that its static critical exponents are the same as in the Landau theory of equilibrium phase transitions [24].The critical exponents obtained from our numerical simulations [23,25] are close to those measured in natural midge swarms [22,26,27].In swarms, midges acoustically interact when their distances are sufficiently small [26].The distribution of speeds is peaked about some value and exhibits heavy tails for large swarms (perhaps due to the formation of clusters) [20].
The statistics of accelerations of individual midges in a swarm is consistent with postulating a linear spring force (therefore a harmonic potential) that binds insects together [20].Laboratory experiments have shown that the swarm consists of a condensed core and a vapor of insects that leave or enter it [28].Swarms of midges form over specific darker spots on the ground (wet areas, cow dung, man-made objects, etc) called markers [29].This empirical fact precludes states of the swarm that are invariant under space translations, such as the ordered and disordered homogeneous phases that are conspicuous in theories of active matter [18,27].
Other animals in a flock, such as starlings [2], define neighbors topologically not metrically.
Bird rotations propagate swiftly as linear waves [15].Metric-free models may incorporate a distributed motional bias [30], visual and auditory sensing are compared in [31], the influence of time delay is studied in [32], the influence of metric and topological interactions on flocking is studied in [33] and [34] considers a swarming model based on effective velocity-dependent gravity.
Here we consider the 2D HCVM, which may model large vertical insect swarms [26].As in the 3D case, the 2D HCVM exhibits scale-free chaos on a curve in the parameter space of alignment noise and spring constant (confinement parameter) for finitely many particles.This curve separates chaotic states in which the swarm is split into several subgroups from chaotic single-cluster states [23].On the critical curve, the swarm size is proportional to the correlation length, which is the only length scale entering power laws for macroscopic quantities.Thus, the swarm is scale free on the critical curve.Using the finite size and dynamical scaling hypotheses [13,14], we can extrapolate power laws obtained for finitely many particles to the phase transition comprising infinitely many particles [15,22,26,27], which we call scale-free-chaos phase transition [23].
For infinitely many particles, the confinement parameter and the largest Lyapunov exponent are both zero [23].In 3D, there are other critical curves that coalesce at the same rate to the line of zero confinement as the number of particles N → ∞.For finite N , these critical curves form an extended criticality region, which may describe data from natural swarms [25].We calculate the static and dynamic critical exponents and show their relation to the 3D ones.We show that the origin of the parameter space (zero confinement and zero noise) is an organizing center that helps understanding the scale-free-chaos phase transition.In particular, it is possible to calculate one static critical exponent by varying the noise at fixed N instead of the traditional approach of varying N for confinement on the critical curve at fixed noise strength [25].
The rest of the paper is organized as follows.Section II describes the HCVM and the deterministic attractors that appear at different values of the confinement strength for zero noise in numerical simulations.We find periodic, quasiperiodic and chaotic attractors.To distinguish the latter, we chiefly use the Benettin algorithm (BA) [35] to compute regions in parameter space where the largest Lyapunov exponent (LLE) is positive.The alignment noise of the HCVM changes the attractors as explained in Section III.Section IV describes the extended criticality region of the scale-free-chaos phase transition and the associated power laws with critical exponents determined from our numerical simulations.For this purpose, we use the critical curve separating single and multicluster chaos: on it, we fix the strength of the noise and vary the number of particles N .
In Section V, we determine critical exponents by fixing N and varying the noise strength on an interval of small values.Section VI discusses our results and contains our conclusions.

II. HCVM AND ITS DETERMINISTIC ATTRACTORS
In 2D, the HCVM consists of N particles with positions x j (t) = (x j (t), y j (t)) and velocities v j (t) = v (cos θ j (t), sin θ j (t)), j = 1, . . ., N , that are updated at discrete times m∆t, m = 0, 1, . .., according to the rule, In Eq. (1b), we sum over all particles (including j) that, at time t, are inside a circle of influence with radius r 1 R 0 centered at x j (t); r 1 is the time-averaged nearest-neighbor distance within the swarm and R 0 is the dimensionless radius.At each time, ξ j (t) is a random number chosen with equal probability in the interval (−η/2, η/2).
We nondimensionalize the model using data from the observations of natural midges reported in the supplementary material of Ref. 22.We measure times in units of ∆t = 0.24 s, lengths in units of the time-averaged nearest-neighbor distance of the 20120910 A1 swarm in Table I [22], which is r 1 = 4.68 cm, and velocities in units of r 1 /∆t, whereas v = 0.195 m/s.The nondimensional version of Eq. ( 1) has ∆t = 1, r 1 = 1, speed v 0 and confinement parameter β given by For the example we have selected, v 0 = 1, whereas other entries in the same table produce orderone values of v 0 with average 0.53.For these values, the HCVM has the same behavior as described here.Thus, our HCVM describing midge swarms is far from the continuum limit v 0 ≪ 1.
Cavagna et al consider a much smaller speed for their periodic VM, v 0 = 0.05, closer to the continuum limit where derivatives replace finite differences [22].Thus, the nondimensional equations of the HCVM are x j (t+ 1) = x j (t) + v j (t + 1), v j = v 0 (cos θ j , sin θ j ), (2a) Eq. (2b) can be written equivalently as Here Θ(x) = x/|x| and ξ is a random number selected with equal probability on the specified interval of width η.
The HCVM has chaotic attractors among its solutions for appropriate values of the parameters.
We identify these attractors by calculating the largest Lyapunov exponent (LLE), which is positive for chaos.To this effect, we use the Benettin algorithm (BA) [35].We have to simultaneously solve Equations (2) and the linearized system of equations δx i (t + 1) = δx i (t) + δṽ i (t + 1), i = 1, . . ., N, (3a) in such a way that the random realizations R η are exactly the same for Equations ( 2) and (3).Here I 2 is the 2D identity matrix.The initial conditions for the disturbances, δx i (0) and δṽ i (0), can be randomly selected so that the overall length of the vector δχ = (δx 1 , . . ., δx N , δṽ 1 , . . ., δṽ   Close to these solutions, there are quasiperiodic ones.At β = 0.1, quasiperiodicity has evolved to chaos, as shown in the last panel of Fig. 1. Fig. 2 shows the bifurcation diagram of the center of mass (CM) coordinates X + Y versus At decreasing values of β, there are P2, P4, P6, quasi-periodic solutions interspersed with larger-period solutions (e.g., P14), and different chaotic attractors for smaller β.The CM trajectories illustrate the different types of attractors in the HCVM.

III. NOISY ATTRACTORS OF THE HCVM
In this section, we describe the attractors of the noisy HCVM and characterize the regions of deterministic and noisy chaos in parameter space.To this purpose, we define below scaledependent Lyapunov exponents (SDLEs) from the time traces of the CM position [23].
As the confinement parameter decreases, different attractors are shown in the left and right columns of Fig. 3 corresponding to two values of the noise strength, η = 0.1, 0.5, respectively (N = 5000).For η = 0.1 and β = 5 × 10 4 , β = 2 × 10 4 , the solutions are period-2 noisy cycles consisting of subgroups of particles forming different patterns alternating densely populated regions with sparsely populated regions.For β = 500, there are noisy quasiperiodic attractors consisting of a densely populated inner region and a number of orbits comprising different subgroups of particles.This pattern is similar to that in Fig. 1 for β = 500.For the larger noise value η = 0.5 and β = 10 9 , the right column of Fig. 3 shows that the annulus with 6 densely populated regions seen on the left column has been filled.For η = 0.5 and β = 5000, a rotating chaotic pattern with dense regions appear.These chaotic patterns change their shape for β = 1000 and β = 100.
Adding the components of X(t), we form the time series x(t) = X(t) + Y (t).To calculate the SDLE, we construct the lagged vectors: The simplest choice is m = 2 and τ = 1 (other values can be used, see below).From this dataset, we determine the maximum ε max and the minimum ε min of the distances between two vectors, ∥X α − X β ∥.Our data is confined in [ε min , ε max ].Let ε 0 , ε t and ε t+∆t be the average separation between nearby trajectories at times 0, t, and t + ∆t, respectively.The SDLE is The smallest possible ∆t is of course the time step τ = 1, but ∆t may also be chosen as an integer larger than 1.The following scheme yields the SDLE [36]: Find all the pairs of vectors in the phase space whose distances are initially within a shell of radius ϵ k and width ∆ϵ k : We calculate the Lyapunov exponent (5a) as follows: where ⟨⟩ k is the average within the shell (ε k , ε k + △ε k ).The shell dependent SDLE λ(ε) in Fig. 4 displays the dynamics at different scales for τ = 9 and m = 6 [36].In deterministic chaos, λ(ε) > 0 presents a plateau with ends ε 1 < ε 2 ≪ 1, in noisy chaos, this plateau is preceded and succeeded by regions in which λ(ε) decays as −γ ln ε, whereas it shrinks and disappears when noise swamps chaos.As η increases, λ(ε) first decays to a plateau for η = 0.1.A criterion to distinguish (deterministic or noisy) chaos from noise is to accept the largest Lyapunov exponent as the positive value at a plateau (ε 1 , ε 2 ) satisfying For η = 0.5, the region where log 10 (ε 2 /ε 1 ) = 1/2 is marked in Fig. 4 by vertical lines.Plateaus with smaller values of log 10 (ε 2 /ε 1 ) or their absence indicate noisy dynamics [36].This occurs for η = 1.
whose slope near the origin gives the LLE [38], as shown in the central panel of Figure 5.In Eq. ( 6), the brackets indicate ensemble average over all vector pairs with ∥X α − X β ∥ < r * for an appropriately selected small distance r * .The LLEs are λ 1 = 0 (P-14 and NQP) and λ 1 = 0.023 for η = 0.5.
For the SDLE to produce accurate values of the LLE, we need sufficient lagged coordinates for safely reconstructing the chaotic attractor.This is achieved if the dimension of the lagged vectors is twice the fractal dimension D 0 or larger [37].Thus, we have used m = 6 lagged coordinates and τ = 9 in Fig. 5(c) and then GZA and SDLE yield the same values of the LLE.Moreover, P-14 and NQP solutions have zero and negative slope of λ(ϵ) at the beginning of the ϵ interval, denoting deterministic and noisy solutions, respectively.

IV. SCALE-FREE CHAOS PHASE TRANSITION
In the previous sections, the numerical simulations of the 2D HCVM have shown the existence of different attractors for fixed values of the number of particles N , confinement strength β and noise η.Fig. 6 shows the phase diagram of different attractors on the plane (β, η).As in the case of the 3D HCVM [23][24][25], as β → 0 and N → ∞, there is a phase transition of scale-free chaos.
At finite N , scale free means that the correlation length is proportional to the size of the swarm, so that all other characteristic lengths are not important.There are several critical curves on the phase diagram that tend to β = 0 at the same rate as N → ∞.We describe here the critical curve β c (N ; η) separating single cluster from multicluster chaotic attractors [23].There are other critical curves that have been studied for the 3D HCVM at fixed η and increasing valued of N : (i) The critical curve β m (N ; η) of local maxima of the susceptibility as a function of β; (ii) the curve where In Eq. (7a), δ(r−r ij ) = 1 if r < r ij < r+dr and zero otherwise, and dr is the space binning factor.
The averages are over time and over five independent realizations corresponding to five different random initial conditions during 10000 iterations [23].The SCCF is the equal time connected correlation function C(r) = C(r, 0) given by Eq. (7a).Note that The correlation length ξ can be defined as the first zero of C(r), r 0 , corresponding to the first maximum of the cumulative correlation function [26]: where θ(x) is the Heaviside unit step function.For r larger than the swarm size, The susceptibility χ is the value of Q(r) at its first maximum, as in Ref. 26.Alternatively, we can use the Fourier transform of Eq. (7a), and define the critical wavenumber k c =argmax k Ĉ(k, 0), the susceptibility as χ = max k Ĉ(k, 0), and the correlation length as ξ = 1/k c [15,22,23].It turns out that k c ∝ 1/r 0 on critical curves and we can use either the real-space or the Fourier space SCCF to find correlation length and susceptibility.On the critical curves where the correlation length is proportional to the swarm size, the correlation length and the susceptibility obey power laws with critical exponents ν and γ, respectively: As N → ∞, the confinement on the critical curve β tends to zero and correlation length and susceptibility diverge [11,13].Similarly, the polarization, i.e., the average speed of the CM velocity, tends to zero as a power law with critical exponent b: For the DCCF, the dynamic scaling hypothesis implies Here z is the dynamic critical exponent and the correlation time τ k = k −z ϕ(kξ) of the normalized DCCF (NDCCF) (11) at wavenumber k obtained by solving the equation: [22,23] tmax t=0 At k c =argmax k Ĉ(k, 0), the correlation time τ kc is a function of β, η and N .For fixed N and η, there is a value of β = β c at which τ kc reaches a minimum.This minimum correlation time corresponds to the smallest time t m (β; η, N ) at which Ĉ(k c , t) = 0 [23].See Figure 5(a) of Ref. [23] for the 3D HCVM.

B. Scale free chaos and critical exponents
For finite N , the 2D HCVM exhibits a region of chaos and noisy chaos for small values of β and η as shown in Fig. 6.Deep into this region, there is a transition between two types of chaotic attractors: multicluster chaos, in which the swarm is split into several groups of particles as in The polarization order parameter goes to zero as β → 0 following the power law Eq.( 10b) with b = 0.42 ± 0.02; see Fig. 8(a).Moreover, the chaotic attractor is multifractal [37], as shown by Fig. 8(b).Then some regions of the chaotic attractor are visited more often than others, which indicates that different length and timescales coexist within the attractor.This is made manifest by calculating the multifractal dimension D q .After a long transient (30 000 time steps), a set of M values of the CM position ⃗ x i = X(t i ) + Y (t i ), i = 1, . . ., M , form a Poincaré map of the attractor.
The multifractal dimension is where θ(x) is the Heaviside unit step function, M ≈ 70000, and C q (r) is the generalized correlation function.D 0 , D 1 and D 2 are the box counting (capacity) dimension, the information dimension and the correlation dimension, respectively.As we vary q, different regions of the attractor will determine D q .D ∞ corresponds to the region where the points are mostly concentrated, while D −∞ is determined by the region where the points have the least probability to be found.If D q is a constant for all q, the CM trajectory will visit different parts of the attractor with the same probability and the point density is uniform in the Poincaré map.This type of attractor is called trivial, whereas a non constant D q characterizes a nontrivial attractor with multifractal structure.As explained above, β c =argmin β τ kc (β; η, N ), where τ kc is the minimum correlation time given by the solution of Eq. ( 12) for fixed η and N .At β c , the correlation length and the size of the swarm are proportional for different values of N as shown by Fig. 9(a).This indicates that the chaotic attractor of the HCVM is scale free for the curve β c (η; N ).The dynamical critical exponent relating correlation time and length, τ k ∼ ξ z ∼ k −z c , is z = 0.99 ± 0.03; see Fig. 11(a).For different values of N , the NDCCF g(t) of Eq. ( 11) is shown in Fig. 11(b).As shown in Fig. 11(c), the different NDCCF curves collapse into a single curve for small values of the scaled time k z c t.The same collapse of the NDCCF as a function of k z c t only for 0 < k z c t < 4 occurs using data from natural swarms, as shown in Figures 2a and 2b of Ref. 22 for z = 1.2.We ascribe the partial collapse of the NDCCF to the multifractality of the chaotic attractors shown in Fig. 8(b), which indicates that different length and time scales persist for all the critical values of β.Thus, a single rescaling of time as in Fig. 11(c) cannot collapse the full NDCCF in our simulations.Furthermore, Fig. 11(d) indicates that the LLE scales with β = β c as with φ = 0.29±0.02= zν.As for the 3D HCVM, these scaling laws illustrate the scale-free-chaos phase transition [23].

V. UNFOLDING OF THE PHASE TRANSITION AT SMALL NOISE
The scale-free-chaos phase transition in the 3D HCVM is such that all critical curves tend to β = η = 0 for finite number of particles [25].There are power laws in η that allow the calculation of critical exponents without having to increase N as we did in Section IV.The point β = η = 0 acts as an organizing center of codimension two for the phase transition in a sense that reminds of singularity theory [40].The critical curves including η = 0 and β c (η; N ) issue from the organizing center at finite N and the attractors in the regions between these lines have specific properties (nonchaotic, single-cluster chaos, multicluster chaos, deterministic chaos, etc).Two parameters, β and η, unfold these behaviors as they take on nonzero values.As N → ∞, all critical curves tend to β = 0 (scale-free-chaos phase transition).On the critical curve β c , the 2D power law is [25]  We can also obtain the dynamical critical exponent by scaling time in the graph of the ND-CCF g(t) to k z c t, with k c = 1/ξ, as shown in Fig. 13(c).Least square (LS) and reduced major axis (RMA) regressions [27] shown in this figure produce similar values of the dynamical critical exponent, z LS = 1.18 ± 0.11 and z RMA = 1.24 ± 0.12, respectively.Figs.As the number of particles goes to infinity, the LLE tends to zero according to the scaling law Eq.(14).Using noise values between 0.1 and 0.5 for N = 100, we find φ = 0.46 ± 0.01; see Fig. 13(f).This exponent is compatible with the dynamic relation φ = zν found for the 3D HCVM in Ref. [23].However, this value is different from that obtained in Section IV by considering increasing values of N at fixed η = 0.5.

VI. DISCUSSION
We have numerically simulated the 2D HCVM for different values of the noise strength η and the number of particles N .As the confinement parameter decreases, the model exhibits periodic attractors with periods 2, 4, and so on, followed by quasiperiodic attractors and, eventually, chaotic attractors with different shapes that occupy regions of finite area and comprise one or several groupings.The noise modifies these attractors.Using SDLEs, we can distinguish essentially deterministic chaos from chaos modified by noise (noisy chaos) and from predominantly noisy signals [36]; see Fig. 4 and the phase diagram in Fig. 6.Within the parameter region of deterministic and noisy chaotic solutions, there is a phase transition at β = 0, N = ∞.To study this transition using numerical simulations, we use the finite size scaling and dynamical scaling hypotheses [11, 13-15, 22, 26]: on critical curves, the characteristic timescale, and the static and dynamic connected correlation functions depend on the control parameters β and η only through the correlation length ξ.On critical curves, the swarm size and the correlation length are proportional, hence ξ ∝ N 1 d (d is the space dimension, 2 in the present work).Finite-size scaling allows us to extrapolate power laws of macroscopic variables obtained for finite N to the case of infinitely many particles, which characterize phase transitions [13].
For finite N , there is an extended critical region comprising different critical curves on the plane (η, β) that converge to β = 0 at the same rate as N → ∞.On the critical curves, the chaotic attractors are multifractal, spanning many different length scales, even as N increases.Compared to the 3D case [25], the confinement parameter on the critical curves is much smaller.In fact, the BA algorithm used to find the LLE ceases to converge before we succeed in finding a zero value of the LLE.This precludes studying the curve β 0 (N ; η) separating chaotic from non-chaotic regions, which necessarily lies below the critical curve β c (N ; η) separating single from multicluster For values on this critical curve, the chaotic attractor occupies a connected finite region of space.
It is shaped as a condensed core plus particles that enter and exit the core, as shown in Fig. 7(b) and 7(c).This is similar to observations of midge swarms in the laboratory [28] With these limitations, we have used our numerical simulations to evaluate static and dynamical critical exponents of the scale-free-chaos phase transition in 2D.We have characterized the critical curve β c (N ; η) and computed the static critical exponents ν and γ associated to correlation length and susceptibility, respectively, as well as the dynamical critical exponents z and φ.To derive the power laws and calculate critical exponents, we have used two limiting procedures.In the traditional one, we keep η fixed and calculate β c (N ; η) for increasing values of N , using to calculate γ (corresponding to the susceptibility) that the curve of local maxima of the susceptibility, β m (N ; η), collapses to β = 0 at the same rate as β c (N ; η); see Fig. 9(d).The other procedure follows from the fact that all critical curves tend to η = 0, β = 0, which produces power laws in η for fixed N [25].Using this second procedure as in Eq. (15b), we have found a value of ν (corresponding to the correlation length) compatible with that obtained by the first procedure.For the static critical exponent γ, the proportionality coefficient between χN −γ/(dν) and η is independent of η and therefore we cannot use the power laws in η to find it [25].
The exponent z characterizes the collapse of the normalized dynamical correlation function when expressed as a function of k z c t = t/ξ z .Unlike the case of critical dynamics near equilibrium, the NDCCF collapses only for an interval of small rescaled times (of width about 4), not for all rescaled times.This partial collapse is also observed in experiments [22], but not in models that try to explain experiments by an ordering phase transition between spatially homogeneous phases [22,27].We ascribe this partial collapse of the NDCCF to the multifractal nature of chaotic attractors, which contain many different length and time scales [23].As the scale-free-chaos phase transition occurs as β → 0, N → ∞, the LLE of the chaotic attractor tends to zero and it follows a new power law while doing so.This power law has a critical exponent φ, which we have also calculated.
How can we compare our results with observations of swarms of midges?On the qualitative static critical exponents.In Ref. [27] in addition to z = 1.35, ν = 0.748, γ = 1.171 are obtained, when the observed values are ν = 0.35 ± 0.10 and γ = 0.9 ± 0.2 [26].Furthermore, these theories fail to predict the limited collapse of the NDCCF [22], or the shape of the swarm [26,28].
In conclusion, we have analyzed the scale-free-chaos phase transition of the 2D HCVM based on numerical simulations of values of noise and number of particles on the critical curve β c separating single from multicluster chaotic attractors.The shape of the swarm (condensed core plus a vapor of particles entering and leaving it) and the partial collapse of the NDCCF in terms of rescaled time are the same as both 3D simulations and experimental data.The value of the static critical exponents ν and γ are close to those obtained from simulations of the 3D HCVM and from experimental data.The dynamical exponent z is different from that of the 3D HDVM and from experiments.We could not investigate the critical curve β 0 (N ; η) separating chaotic and non-chaotic attractors due to non-convergence of the Benettin algorithm used to calculate the LLE at values of the confinement parameter that are too small.It would be interesting to study the HCVM having different confinement parameters in the vertical and transversal directions.This would be closer the observations of natural swarms, which are elongated in the vertical direction, and the results might be interpolations between the 2D and 3D models.These interpolations could also be useful to study renormalization group properties of the scale-free chaos phase transition based on the quasiperiodic route to chaos [39,41].

FIG. 1 .
FIG. 1. Visual patterns for N = 5000, R 0 = 0.472, η = 0 and decreasing values of β.Period-2 solutions occur at β = 10 9 and at β = 2 × 10 5 : the particles remain inside the unit circle forming either connected patterns or patterns with large holes (reminding of the vibration mode of a plate).At β = 10 4 , there are two subgroups within the unit circle that exchange their positions in a period-4 solution.At β = 500, there are large period solutions with subgroups of particles inside the unit circle and others oscillating outside, forming orbits reminiscent of the Bohr atomic model (the different colors correspond to different times).At β = 0.1, the pattern corresponds to a chaotic attractor displaying many more subgroups (of smaller density) and also single particles.

Figure 1
Figure1shows different attractors for N = 5000 in the deterministic case η = 0. Initially, the particles are randomly placed within a sphere with unit radius and the particle velocities are pointing outwards.As β decreases from very large values, we observe, from top to bottom and from the left to the right, period two (P2) solutions, period-four (P4) solutions, quasiperiodic (QP), large period solutions and chaotic solutions.At β = 10 9 , the particles remain inside the unit circle forming a symmetric pulsating pattern that repeats itself every two iterations.P2 solutions also occur for β = 2 × 10 5 but the particles form a pattern with unoccupied spaces within the unit circle.At β = 10 4 , the period has duplicated, the particles form two subgroups and pass from one to the other, repeating their positions every four iterations.At β = 500, there are large period solutions with subgroups of particles inside the unit circle and others oscillating outside, forming orbits reminiscent of the Bohr atomic model (the different colors correspond to different times).

FIG. 2 .
FIG. 2. Bifurcation diagram of the center of mass coordinates X + Y versus β for N = 128.At different values of confinement there are P2, P4, P6, quasi-periodic solutions interspersed with larger-period solutions (e.g., P14), and different chaotic attractors for smaller β.
inflection points of the susceptibility vs β; (iii) the curve β 0 (N ; η) separating regions of zero and positive LLE.These curves satisfy the relation β 0 (N ; η) < β c (N ; η) < β i (N ; η) < β m (N ; η)[23,25].In 3D, all these critical curves converge to β = 0 at the same rate as N → ∞[23,25].In 2D, we have not found a curve β 0 (N ; η) because the algorithms used to compute the LLE do not converge at values of β that are too small: for all β < β c (N ; η) within the range of convergence of the BA, the LLE is positive.A. Correlation functions and power lawsLet us first define the static and dynamic correlation functions, the correlation length and time and the susceptibility with the corresponding power laws and critical exponents.The dynamic connected correlation function (DCCF) is[15,26]

Fig. 8 (FIG. 9 .
Fig. 8(b) shows that the box-counting dimension D 0 and D q for q > 0 undergo a downward trend with increasing N (decreasing β c ).Then the dimension of the more commonly visited sites of the attractor decreases.The chaotic attractor remains multifractal when β → 0 as N → ∞ and chaos disappears: different time scales persist [39].

Figs. 9
Figs. 9(b) and 9(c) show that, for increasing N → ∞, correlation length and susceptibility scale with β = β c (with lim N →∞ β c (η; N ) = 0) following the power laws (10a) with critical exponents ν = 0.30 ± 0.02 and γ = 0.78 ± 0.05, respectively.While the correlation length and time are calculated using numerical simulations of the HCVM on the critical curve β c (η; N ), the susceptibility does not tend to a definite value as N increases due to the very small values of β c for large N (which are much smaller for the 2D HCVM than the corresponding values for the 3D HCVM).To fix this problem of the 2D HCVM, we recall that the scale-free-chaos phase transition in the 3D HCVM has an extended critical region which collapses to β = 0 at the same rate for all values of the noise as N → ∞[23].The critical curve with larger values of β is β m (η; N ), corresponding to the local maximum of the susceptibility χ with respect to β (at fixed N and noise) calculated as in χ in Equations (8), χ = Q(r 0 ), or as χ = max k Ĉ(k, 0)[23].The curve β m (η; N ) is also scale free and it yields scaling laws as in Equations(10).Figs.10(a) and 10(b) illustrate that the power law Eq.(10a) for the susceptibility χ = Q(r 0 ) at β = β m .Moreover, Fig. 9(d) shows that the curves β c and β m tend toward β = 0 at the same rate: β c /β m → 0.035 as N → ∞.Thus, we calculate the critical exponent of the susceptibility using β m (η; N ) instead of β c (η; N ), thereby obtaining Fig. 9(c).
13(d) and 13(e) shows that the curves for different values of η ∈ [0.1, 0.5] visually collapse when time is rescaled with z ≈ 1.15.