Instabilities and chaos in a kinetic equation for active nematics

We study dry active nematics at the kinetic equation level, stressing the differences with the well-known Doi theory for non-active rods near thermal equilibrium. By deriving hydrodynamic equations from the kinetic equation, we show analytically that these two description levels share the same qualitative phase diagram, as defined by the linear instability limits of spatially-homogeneous solutions. In particular, we show that the ordered, homogeneous state is unstable in a region bordering the linear onset of nematic order, and is only linearly stable deeper in the ordered phase. Direct simulations of the kinetic equation reveal that its solutions are chaotic in the region of linear instability of the ordered homogeneous state. The local mechanisms for this large-scale chaos are discussed.


Introduction
The collective properties of assemblies of active particles are currently under intense scrutiny, both at the theoretical and experimental level. Following the seminal works of Vicsek et al, and Toner and Tu [1,2], which were published in 1995, physicists have made progress towards uncovering and understanding universal, generic features of groups of interacting self-propelled particles, some of these results have started to be confronted to experimental data. This is the case, for instance, in the long-range correlations and anomalous fluctuations that are predicted by Toner and Tu [3] and by Ramaswamy et al [4][5][6] to characterize most homogeneous, orientationally-ordered, fluctuating phases of flocks. They have been observed in numerical simulations of simple models [7] and they have been argued to have been seen in vibrated granular particles [8][9][10], bacterial swarms [11,12], as well as in vitro motility assays, where grafted motor proteins displace biofilaments [13].
The Vicsek et al model [1] now serves as the reference minimal model. In this context, the collective motion is reduced to the competition between noise and the local ferromagnetic alignment of the polarities of constant-speed particles. All other possible interactions are neglected, as is the fluid in which the particles evolve (as opposed to active suspensions). Even within this somewhat restrictive, 'dry flocking', minimal context, other classes of active particles systems can be defined, varying the symmetry of the particles (polar or apolar) and of the interaction (ferromagnetic versus nematic alignment). The case of active nematics [14] was considered very early on by Ramaswamy et al [4][5][6], and has recently enjoyed a lot of attention. This may be understood here as assemblies of apolar particles that are forced to jiggle along their director and aligning nematically. Ramaswamy and collaborators, working on a linearized theory, predicted the existence of a non-equilibrium current at the origin of particularly strong number fluctuations [4][5][6].
The active nematics version of the Vicsek model was introduced in [15]. In this model, point particles possess a (uniaxial) nematic degree of freedom. They are displaced at discrete time steps along one of the two directions defined by the director and are calculated from the nematic average of local neighbors' orientation, to which some random noise was added. This model can be seen as a toy model for active, self-propelled rod-like objects interacting by steric inelastic collisions. The numerical work of [15] revealed an isotropic/nematic transition to quasi-long-range order, similar to the equilibrium Kosterlitz-Thouless transition of the XY model. The quasi-ordered phase was found to be spontaneously segregated into a large highdensity high-order band fluctuating in a sparse, disordered sea, and is shown to exhibit the strong number fluctuations predicted by Ramaswamy et al.
More recently [16], the hydrodynamic equations describing dry active nematics were derived from the Vicsek-like model of [15], and were found to coincide with those obtained by different approaches [19,20]. These equations, derived via a Boltzmann equation obtained in the dilute limit by kinetic theory techniques, were found to possess features at odds with those reported in the original Vicsek-like microscopic model. For instance, a homogeneous ordered phase is predicted to exist deep in the ordered phase whereas it was not observed in the microscopic model.
Here, we study dry active nematics directly at the level of the Boltzmann kinetic equation level in order to clarify the above issues. We consider a slightly more complicated (but more realistic) case where particles also undergo positional diffusion. Stressing the differences with the well-known Doi theory for non-active rods near thermal equilibrium, we write a simple kinetic equation (section 2) and study the linear stability of its homogeneous solutions (section 3). In the same section, we also derive the corresponding hydrodynamic equations governing the nematic order field and show analytically that these two description levels share the same qualitative phase diagram, but quantitatively diverge rather strongly as the global density is increased. We show, in particular, that the ordered, homogeneous state is unstable in a region bordering the linear onset of nematic order and is only linearly stable deeper in the ordered phase. In section 4, we perform direct simulations of the kinetic equation and reveal that its solutions are chaotic in the region of linear instability of the ordered homogeneous state. We discuss the local mechanisms at the origin of this large-scale chaos. In section 5, our results are summarized and discussed in view of the above-mentioned discrepancies between the reported behavior of the microscopic Vicsek-like model and the properties of continuous descriptions of active nematics.

Kinetic equation for dry active nematics
One salient feature of the Vicsek-style model for active nematics that was developed by Chaté et al [15] is that the rotations of particles are governed by local nematic alignment while the streaming (translational movement) is implemented at constant speed independently of interactions. Similarly, real 'self-propelled' particles are driven along their long axis, inducing strong longitudinal diffusion [21][22][23][24][25][26].
Let us thus consider elongated objects interacting by local nematic alignment and suppose that every particle is driven along its long axis at speed v 0 , with the driving force switching the sign of the velocity randomly with rate k. One can then write, in the mean-field limit, the kinetic equations for the two one-particle number distribution functions ± f t r u ( , , ) for particles of center of mass located at r moving in ±u directions where π ∈ u [0, 2 ] is the unit vector along their long axis. Given that the two subpopulations + f t r u ( , , ) and − f t r u ( , , ) constantly exchange members at rate k, we have: t r r 0 where D r is the rotational diffusion constant, the rotational operator  is defined by  = × ∂ u u , and w r u ( , ) is a self-consistent interaction potential which has ±u-symmetry [27][28][29]. In the two-dimensional case studied here, w r u ( , ) commonly takes the form of the excluded-volumelike interaction for infinitely thin rods [27,28], 2 where l is the particle length and we have neglected the high-order spatial derivatives in the integration.
Equations (1) and (2) When the switching rate k is large, we expect the difference between + f r u ( , ) and − f r u ( , ) to be small. In particular, the damping term − kf r u 2 ( , ) M is dominant in equation (5) when ≫ k D r and f M relaxes quickly to a quasi-steady-state. In such a case, we can assume to the first order in spatial derivatives, which is assumed to be small if ≫ k v l / 0 . Substituting this result into equation (4) is the effective longitudinal diffusion coefficient. Similarly, in the presence of some transverse diffusion perpendicular to particleʼs long axis, we can write: where ⊥ D is the perpendicular component of the translational effective diffusion rate referring to particleʼs long axis and . For simplicity, in the following we rescale length scales by the particleʼs length and set l = 1. Equation (7) confirms that the excluded-volume-like potential w r u ( , ) codes for nematic alignment but has no effect on the translational motion of particles. A system of physical rods interacting by volume exclusion can only exhibit this property if some driving forces are applied to the particles to overcome the effect of mutual interaction on their spatial translocation. This is different from the well-known dynamic mean field theory of non-active rods near thermal equilibrium [27], where excluded-volume interactions also affect translational motion.

Instabilities of homogeneous solutions
We now proceed to study the linear stability of the homogeneous solutions to the kinetic equation (7). We first note that equation (7) is invariant under reversal of u, as expected for nematics. We can thus expand it in Fourier series of the angular variable θ (the argument of u) as:  (7) can be rewritten as a double hierarchy of equations [16,30] given by: t p x y n x y x y 0  The coefficients controlling the isotropic/nematic transition differ because different interaction rules are implemented. Here we have only one control parameter for the transition, the particle density. For Vicsek-style models, the noise-level is also an explicit parameter that dominates the nematic ordering transition.

Spatially homogeneous solutions
We first determine the spatially-homogeneous solutions to the kinetic equation. An obvious example is the disordered homogeneous state . To find the nematically-ordered solution, we neglect all spatial derivatives in equations (9)- (11). Assuming that the nematic director is along the x-axis, we can further set = q 0 n , so that we are left with simple ordinary differential equations for coefficients p n : The fixed point solution = * p p n n of this hierarchy can be found numerically by keeping a large number of modes and setting the others to zero.
At order 3, an explicit solution is easily found, which gives the density dependence of the magnitude of nematic order. Conventionally, the tensorial order parameter for homogeneous nematics is defined as where ∫ ρ = f u u d ( ) is the density andn is the unit vector of the nematic director. If we letn align with the x-axis, that isˆ= n 1 x andˆ= n 0 y , from equation (13), we have . Setting = p 0 n for ⩾ n 3, we obtain an equation for nematic order parameter S(t), where ρ ρ ρ = * / is the number density rescaled by the critical density ρ π = * 3 /2. For ρ < 1, the disordered solution S = 0 is stable. For ρ > 1, it is unstable and is replaced by the homogeneous solution 2 which is then stable. In other words, the transition is continuous in the mean field homogeneous limit and ρ ρ ∝˜−S ( 1)/ . Since we only include the excluded-volume-like interaction equation (3) in equation (7), equation (15) shows that the magnitude of nematic order is independent of driven activity as given by parameters ⊥ D and ∥ D . In dry active nematics, enhanced ordering effects can be accounted for if we explicitly write down how the active processes modify the rotation of particles upon collision, as done for example in [29,31].
In figure 1, we show the difference between the ordered homogeneous solutions of the full kinetic equation (7) (obtained numerically) and solutions of the hierarchy (12) keeping 300 modes and 3 modes (solution (15) above). No difference can be seen between the full numerical solution and the high-order truncation but, as expected, the cubic solution only matches those near the transition line.
The stability of the homogeneous solutions can be studied at two levels: numerically at the 'full' kinetic equation level (with a high-order truncation), and analytically at the hydrodynamic  (7), using the algorithm described in the appendix.
We have used a 3 × 3 spatial lattice, which is small enough to ensure that no long-wave length instability takes place in the system. The red curve, which coincides with the black one, was generated by solving numerically equation (12) truncated at 300 modes. The overlap of these two curves, which are generated by different approaches, testifies to the accuracy of our numerical kinetic equation solver. The difference is of order − 10 3 . The blue curve is generated by equation (15) (hydrodynamic level). This low-level truncation underestimates the value of S in the high-density regime. equations level (with an appropriate low-order truncation and closure scheme). Below we first proceed with the hydrodynamic approach, following closely the recent work of Bertin et al [16] (see also [17][18][19][20]30] for polar particle and other cases). We then treat the full kinetic equation level.

Stability at the hydrodynamic level
In order to work with a well-controlled expansion, we follow the 'Boltzmann-Ginzburg-Landau' approach used in [16,30,32]. We work in the vicinity of the main instability line: where ϵ is a small parameter. We further use the scaling ansatz δρ ϵ where δρ is the local deviation from the average density. The simplest non-trivial closure is to discard terms that are higher than ϵ 3 . In this scenario, equation (9) is unchanged, and equations (10) and (11) are truncated according to the scaling ansatz, leaving four equations for p 1 , q 1 , p 2 , q 2 : t p x y n x y r r 1 From equations (18) and (19), we see that p 2 and q 2 are enslaved to p 1 and q 1 as long as . By solving equations (18) and (19) for p 2 and q 2 , substituting their expressions into equations (16) and (17), and discarding ϵ 4 terms, we have π ρ π Equations (9), (20) and (21) constitute the closed hydrodynamic description of equation (7). They are formally identical, up to the positional diffusion terms, to those derived in [16]. They can also be rewritten in tensorial notation. Similarly to equation (13), the local nematic order parameter tensor is given by is the density andn r ( ) is the unit vector of the nematic director [19,28,29,33]. By comparing with the definition of p r ( ) The hydrodynamic equations written in terms of the reduced density ρ ρ ρ = * r r ( ) ( )/ and the alignment tensor αβ Q r ( ) then read: These hydrodynamic equations are similar to those used by Baskaran and Marchetti for analyzing the stability of fluctuations in the nematic state [28], where they not only include the enhanced diffusion by self-propelled motion but also include normal passive diffusion coefficients. One major difference is that, while we explicitly include the density dependence of the magnitude of nematic order, they focus on the deep nematic state and suppose the strength of nematics is a constant parameter. As a result, they did not observe the instability we describe below.
We now turn to the linear stability analysis of the homogeneous solutions of equations (22) and (23). Let us assume that the nematic directorn r ( ) of the homogeneous ordered solution is along the x-axis of the system, that isˆ=n r The linearized hydrodynamic equations read: t n x y r p r To keep σ > 0, which ensures that local nematic order grows with increasing local density, we restrict our analysis to the regime δ < 2. . Finding the eigenvalues of equation (27) is difficult since we need to solve a cubic equation. Here, we only determine the instability line. Let us define is an eigenvalue of the coefficient matrix of equation (27). The homogeneous nematic ordered state becomes linearly unstable if λ ⩾ Re ( ) 0; that is, For the cubic equation = g z ( ) 0, the system either has: (i) one real root and two complex conjugate roots, (ii) one real root and two identical real roots, or (iii) three different real roots.
For cases (i) and (ii) we can prove that only the single real root (not the multiple real roots) of = g z ( ) 0 can cause the linear instability of the system. Let us suppose z 1 , z 2 and z 3 are three roots of = g z ( ) 0, where z 2 and z 3 are the complex roots or the multiple roots and we have = * z z 2 3 . According to Vietaʼs formula, we have: If the linear instability is caused by the complex root or the multiple roots, then ⩾ z Re ( ) 1 2 and it is easy to see that equation (29) and (31) contradict each other. Thus, if the linear instability takes place, it must be caused by a real but non-multiple root.
For case (iii), we can prove that only the largest real root of = g z ( ) 0 can cause the linear instability. Let us suppose z 1 , z 2 and z 3 are three different roots of = g z ( ) 0, and we have > > z z z 1 2 3 . According to Vietaʼs formula, we have: , we must have ⩽ g (1) 0, and = g (1) 0 gives the instability line. By substituting it into the expression of g(z), we have Since ⩽ D 1 0 , the first term is positive. If the second and third terms together are negative, the system can always be destabilized when κ → 0. To push the system into an instability regime, we try to minimize equation (34) with respect to κ and θ, which gives The third inequality can not be satisfied since it contradicts the condition σ δ| | ⩽ D /( ) 1 0 . The instability line is given by The coefficient of κ 4 is always positive since ⩽ D 1 0 , signifying that for large enough wave numbers the fluctuations are always stable. For small wave numbers that describe large-scale fluctuations, the stability is controlled by the coefficient of κ 2 . Thus, when σ θ δ 0, 0, 0, 0, , ). Some typical dispersion relation curves are shown in figure 2(b). Above the orange line, the homogeneous ordered solution (15) is linearly stable.

Stability at the kinetic equation level
The homogeneous ordered solution to the kinetic equation (7) (22) and (23) is linearly stable. The green curve marks the same limit, but for the many-mode truncation of the kinetic equation as given in equations (41) , δ = 0.01. The black curves in polar plots are generated by equation (34) and the red curves are generated by equation (38). These branches in the polar plots show that the instabilities are developed perpendicular and parallel to the nematic director for positive and negative D 0 , respectively. The red dots and blue circles near the yellow curve are located where numerical integration of equation (7)  space: We calculated numerically the maximum eigenvalue of the above system, truncated at the same high-order as for the solution, for wavenumbers = k k close to zero. We find that this eigenvalue changes sign and becomes negative above the green line in figure 2(a). (These results were obtained by keeping 300 modes; we checked that no signifiant change occurs when keeping more modes.) For large D 0 , the region of linear instability of the homogeneous ordered solution is thus much larger for the kinetic equation (7) than for the hydrodynamic equations (22) and (23). The hydrodynamic description is quantitatively valid only in the near threshold region of isotropicnematic transition.

Beyond linear instability: chaotic solutions of the kinetic equation
We now investigate numerically the solutions of the kinetic equation (7) in the parameter region where no homogeneous solution is linearly stable. We use an alternating-direction implicit algorithm, which is unconditionally stable and accurate to second-order for the linear part (see the appendix for details).
We first test the phase boundary predicted by linear instability analysis (green line in figure 2(a)). We integrate the kinetic equation for parameter values around the predicted line, starting from the disordered homogeneous state with small initial noise, in a domain of size 120 with periodic boundary conditions.
As shown in figure 2(a), all runs above the phase boundary (blue circles) end up in the homogeneous nematic state, while those below (red dots) never settle in a stationary state. This confirms the results of the linear stability analysis. Figure 3 shows the patterns found by numerical integration of the kinetic equation, starting from a slightly perturbed homogeneous In the regime of linear instability of the homogeneous nematic state, the system tends to evolve first to a high density nematic band, which later breaks. The two snapshots shown for each run presented are taken in these two stages. The width of this transient nematic band increases with the global density δ. Increasing the anisotropy of diffusion , the nematic band becomes highly concentrated and narrower. All these bands are unstable. The dynamical behavior is analyzed in detail in the next section using a larger periodic domain.

A chaotic, disordered phase?
Next, we investigate the dynamics of the system in the linear instability regime using larger domains. Starting from an isotropic and spatially-homogeneous initial condition, local ordered nematic domains form at the beginning, accompanied by the quick development of density inhomogeneities. Further coarsening of these structures leads to the coexistence of particleenriched nematic domains (where ρ ρ > * ) in a particle-poor isotropic sea (where ρ ρ < * ). The system exhibits large-scale chaos with statistically stationary properties. In figure 4, we show snapshots of the reconstructed density and nematic order fields in this asymptotic chaotic regime. The rescaled density field δ ρ =˜− r r ( ) ( ) 1 represented by the color scale was obtained by integrating the distribution function f r u ( , ) over the angular variable. The magnitude Vicsek-like active nematics model [15]. It is not clear at this stage whether or not this chaotic dynamics eventually gives rise to a bona fide disordered phase in the asymptotic large-size limit. Unfortunately, we are currently unable to simulate much larger systems.
Similarly, the two-point orientational correlation functions display, somewhat surprisingly, algebraic tails typical of quasi-long-range order ( figure 5(c)). This may only be a finite-size effect: if the chaotic dynamics possesses a large but finite correlation length (something we are not able to probe yet), then we expect these correlation functions to have exponential tails.

Fragmentation of nematic bands
We now take a close look at how a nematic band along the x-axis breaks (see figure 6). For > D 0 0 , from our previous stability analysis, the density inhomogeneity develops along θ = − π π , 2 2 ; that is, perpendicular to nematic director n. The term δ − ∂ Initially, the nematic director in the high-density region is almost parallel to the nematic ordered stripe boundary (see figure 6(a), where only a fraction of the space domain is shown). To illustrate the mechanism for the deviation of nematic director, we investigate whether this alignment is stable with respect to small fluctuations of the director , where we have assumed that there is no spatial variation of ρS along x-axis. The spatial variation of ρS along the y-axis is significant since the main density inhomogeneity is developed in that direction. Near the stripe boundaries, we always have a region with ρ ∂˜> S 0 y 2 , which makes the fluctuations δ ⊥ n y unstable. This instability induces changes of the order parameterʼs orientation; thus, the nematic directors in figure 6(b) become oblique to the density profile boundaries after the band forms. This also explains why the nematic directors in figure 4(a), (b) are often found to be oblique to the domain boundaries. In figure 6(b), together with the obliqueness of directors, we observe a leakage of particles from the high-density region via the relatively faster diffusion along the directorʼs alignment. With this decrease the density of ordered nematic stripe enters into the unstable region, and the spatial instability takes place again. Figure 6(c)  show the fluxes around the twisted-spindle shaped region after its formation. The density currents show that it absorbs particles from the medium on both sides and generates outward flux on both tips. In this way, it can grow and extend itself quickly into the low-density medium.

Summary and discussion
To summarize, we have studied the kinetic equation for active nematics in the dilute regime, both analytically and by direct numerical simulations. In deriving the kinetic equation, we deliberately stayed in the fast reversal limit, that is, the reversal rate k in equations (1) and (2) is kept much larger than the self-diffusion rate D r , which leads to a kinetic equation with pure nematic symmetry. When k is comparable to, or even smaller than, D r , such a simple kinetic equation is not expected to be valid. Since self-diffusion itself gives rise to an effective reversal rate, it contributes a term v D / r 2 to the renormalized spatial diffusion coefficient. With small k, we can expect enhanced spatial diffusion, and the drift terms in equations (1) and (2) will induce local polarities enslaved to the nematic order. Furthermore, the implicit numerical algorithm we use to solve the kinetic equation (7) is not suitable for equations with drift terms.
At the level of the linear instability of spatially-homogeneous solutions of the kinetic equation (7), we have found that the homogeneous nematically-ordered solution is linearly unstable to long wavelength transversal modes in a region bordering the transition line defined by the instability of the homogeneous disordered solution (between the blue and green lines in figure 2). This scenario is in agreement with recent studies at the hydrodynamic level [16,19]. Here, deriving and performing analytically the linear stability of the hydrodynamic equations, we have quantified the consequences of the reduction to the hydrodynamic level (the difference between the green and orange lines in figure 2). To our knowledge, this is achieved for the first time in the context of active matter (see, however, the recent preprint by T Ihle [34]).
At the nonlinear level, we have found that in the region of instability of the homogeneous ordered solution, the kinetic equation exhibits chaotic behavior. We have described this chaos as the fragmentation of the larger high-density high-order structures formed by the instability mechanism and the subsequent renewal of such 'spindle-shaped' structures (figures 4 and 6). Clearly, more work will be needed to characterize the nature and the asymptotic properties of this chaos and, in particular, to determine whether or not it constitutes a disordered phase.
Indeed, the chaotic dynamics described above are different from the chaos found at the hydrodynamic level, but our findings are nevertheless reminiscent of the dynamics reported in [15] for the Vicsek-style microscopic model for active nematics, where intermittent destabilization events of a large high-density, high-order macroscopic band were described. As already suggested in the recent work of Bertin et al [16], this Vicsek-style microscopic model, which should be described at the kinetic level by a Boltzmann equation similar to the kinetic equation studied here, will have to be revisited. Indeed, no homogeneous ordered phase was reported for it and the segregated phase with its intermittent destruction of a single large, dense, ordered band was found to be quasi-long-range ordered. While this is not, strictly speaking, at odds with the results presented here (see the relatively-high order parameter values in figure 5(a), and the correlation function in figure 5(c)), this calls for larger-scale simulations of both the microscopic model and the kinetic equation. The hydrodynamic equations (22) and (23) have the same mathematical structure as those derived in [16] and used in [37], except for the difference in dynamic parameters. While these hydrodynamic equations also give rise to large scale chaotic behavior, as shown very recently in [37], the ordered structures found there are different from those arising in the full kinetic equations studied here. In the kinetic equation, the system generates twisted spindle shaped local ordering structures that are unlike the bowlike traveling structures found in hydrodynamic simulations. By keeping the higher modes in equations (9)-(11), we can recover similar division processes of ordered domains. Moreover, the chaos found in the nonlinear regime also differs from that found in the hydrodynamic description of active nematics suspensions, where interesting low Reynolds number turbulence has been reported [38][39][40][41]. The mechanism for linear instability, however, shares similarity, where the density dependence of nematic strength is required for the initiation of the destabilization of ordered state [40].
At this stage, we can nevertheless already come back to the question of 'giant number fluctuations' in active nematics systems. These anomalous fluctuations have become a landmark of active matter studies following the seminal works of Toner and Tu, as well as Ramaswamy and co-workers [2,4,5]. In their study of the Vicsek-style model for active nematics, Chaté et al reported that the standard deviation Δn of the number of particles present in a sub-system containing on average n particles scales anomalously: Δ ∼ n n, in agreement with predictions of Ramaswamy et al based on a simple linearized theory. We have performed the same measurement in the chaotic regime of the kinetic equation studied here, and found the same 'anomalous' scaling ( figure 7). This may appear to be a surprising result since the kinetic equation is not supposed to incorporate noise terms at the origin of giant number fluctuations. In fact, the results obtained here are simply explained by the strong density segregation observed: having sharp interfaces delimiting dense-ordered regions moving in a sparse-disordered medium trivially yields the above scaling without recursing to the Ramaswamy et al calculation. Also note that this calculation, as well as that of Toner and Tu for polar flocks, assumes that one is in the homogeneous ordered phase.
Again, the question of determining whether the chaotic phase reported here is asymptotically ordered or not, and/or whether the quasi-long-range ordered phase reported in [15] is in fact disordered, appears as the next crucial question down the road towards a comprehensive understanding of dilute active nematics.  figure 4 (chaotic phase). As described in [35,36], both spatial-averaging-first (SAF) and temporalaveraging-first (TAF) procedures were used, yielding the same results. To calculate the fluctuations, the domain is divided into a number of small subsystems of equal size. The number of particles in each subsystem is calculated by numerical integration. For the SAF procedure, the mean square deviation of the number of particles in subsystems from the global average was calculated for a single snapshot. Then, an average over all snapshots is performed. For TAF procedure, we first calculate the temporal average of number of particles in a given subsystem. Then the mean square deviation from the average number 〈 〉 N for that subsystem is calculated and averaged over all subsystems. For both SAF (square) and TAF (circle) measurements, here the root mean square deviation ΔN grows proportionally with 〈 〉 N , showing particularly strong number fluctuations. The solid line has unit slope. The curvature at large N is due to finite-size effects.
Appendix. Numerical method Equation (7) can be solved by using the alternating-direction implicit method [42]. Using a uniform spatial and angular grid with = + ▵