A Finite Volume Method for Continuum Limit Equations of Nonlocally Interacting Active Chiral Particles

The continuum description of active particle systems is an efficient instrument to analyze a finite size particle dynamics in the limit of a large number of particles. However, it is often the case that such equations appear as nonlinear integro-differential equations and purely analytical treatment becomes quite limited. We propose a general framework of finite volume methods (FVMs) to numerically solve partial differential equations (PDEs) of the continuum limit of nonlocally interacting chiral active particle systems confined to two dimensions. We demonstrate the performance of the method on spatially homogeneous problems, where the comparison to analytical results is available, and on general spatially inhomogeneous equations, where pattern formation is predicted by kinetic theory. We numerically investigate phase transitions of particular problems in both spatially homogeneous and inhomogeneous regimes and report the existence of different first and second order transitions.


Introduction
Collective motion of groups of multiple agents of various origin is a fascinating phenomenon that manifests itself in versatile systems, ranging from microscopic scale, e.g., colloidal suspensions [46,26,33], through mesoscopic scale, e.g., bacterial suspensions and microtubule bundels [56,20,61], to macroscopic scale, e.g., schools of fish and flocks of birds [42,19,21]. The study of such systems is an active field of research, both theoretically and experimentally. When a system consists of self-propelled agents, i.e., the ones that are able to move without any external forcing, it is referred to as an active matter system. The direct approach to model such systems is to describe the motion of each constituent agent with differential equations. However, when the number of agents is large, this approach becomes computationally expensive and one looks for a respective continuum description.
First approaches for the continuum description of active matter were based on symmetry arguments and conservation laws [58,8] and described the evolution of hydrodynamic variables, e.g., a marginal density function and a polar order field. This approach allowed to reproduce some of the behavior given by an agent-based description but the resulting equations were not linked to the microscopic parameters, thus, not allowing for the respective analysis of the agent-based behavior. As an alternative, the kinetic theory presents a systematic way to construct a continuum description of an agent system via density functions of agent's positions and velocities [15,16,10,13,11,40,44,6]. Mostly, one is interested in a time evolution of a density function of one agent, which is governed by a nonlinear PDE. As a result, the terms in such a PDE do depend on microscopic parameters from the agent-based model, thus, allowing one to study the continuum limit behavior of a particle system in terms of those parameters.
In this paper, we are interested in the construction of effective FVMs for numerical integration of PDEs derived as the continuum limit of nonlocally interacting self-propelled particle systems in two dimensions that takes the following form: where r ∈ R 2 is a position vector, e(ϕ) ∈ S 1 is a unit velocity vector that depends on particle's orientation ϕ ∈ R/(2πZ) =: T, t ∈ R + is time, D ϕ ∈ R + is a rotational diffusion coefficient, f : R 2 × T × R + → R + is a one-particle probability density function which quantifies the probability to find a particle with a position r, orientation ϕ at time t, w is some functional that represents nonlocal interaction between particles. Besides the high dimensionality of the problem, the presence of the nonlocal interaction requires particular attention to the performance of constructed numerical schemes. Numerical integration of such continuum limit kinetic equations for active matter systems is an ongoing research [12,57,51,32,18,7]. For a general survey on the numerical treatment of kinetic equations, we refer the reader to [23]. An important problem in active matter systems is the study of related phase transitions versus model parameters. Its knowledge allows one to analyze, predict, control, and design particle systems with experimentally desirable properties. Depending on the context, one might consider particles of different origins, which in the continuum limit frequently assume polar or nematic representation with polar or nematic interactions [9,22,48,49,43,7,47,41,14]. Phase transitions are commonly quantified in terms of a polar order parameter. It has been shown that depending on different regimes, one might observe first as well as second order transitions. In the present work, one of our goals is to establish the nature of phase transitions between different states of a model of interest by developing an accurate numerical scheme.
The paper is organized as follows. In Section 2, we introduce a particular particle model and its continuum limit PDE, which we are interested to numerically investigate. We also provide some of the analytical results that will be used later in the construction of numerical methods as well as in their performance assessment. In Section 3, we describe the construction of FVMs appropriate for numerical investigation of active matter continuum limit equations in the form of Eq. (1). In Section 4, we demonstrate how our schemes perform in recovering spatially homogeneous and inhomogeneous solutions as well as study related phase transitions. Finally, we summarize results of the present paper and outline the future work in Section 5.

Active Brownian particle flow in the continuum limit
We begin the discussion by first introducing an exemplary particle model that lies in the origins of the current research. The model is a starting point for the definition of continuum limit equations (1), which we want to numerically integrate. Moreover, it serves as a reference point for interpretation of solutions of continuum limit PDEs. Next, we will present a PDE corresponding to the particle model along with its key properties, which will be used to test the performance of finite volume schemes in Section (4).

Finite size particle model
Let U := R/(LZ) and T := R/(2πZ) be one-dimensional spaces with periodic boundaries extending from [0, L] and [0, 2π], respectively. We consider a system of N ∈ N self-propelled particles moving in a two-dimensional domain with periodic boundaries. Each particle is assigned a unique index i = 1, . . . , N . A spatial position of a particle is described with a vector with periodic components r i = (x i , y i ) ∈ U 2 . We assume that particles move in an overdamped regime with a constant speed v 0 ∈ R + . Thereby, their velocity is completely determined by their direction of motion ϕ i ∈ T. Let us denote a unit velocity vector as e(ϕ i ) = (cos ϕ i , sin ϕ i ) ∈ S 1 . As a result, we consider the state vector of a particle i to be p i = (x i , y i , ϕ i ) ∈ U 2 × T =: Ω. We consider particles as interacting Ω-valued processes that are solutions to the following coupled system of stochastic differential equations (SDEs): (2) 2 The first equation described a self-propulsion a particle in the direction of ϕ i . The right hand side of the second equation consists of a nonlocal particle interaction via alignment and the Wiener process, which plays a role of an external perturbation, respectively. The strength of nonlocal alignment is controlled with a parameter σ ∈ R + . Particles interact with their nonlocal neighborhood, which consists of all particles within the distance ∈ U and is defined as and |B i | denotes the neighborhood's cardinality. We postulate that particles' alignment mechanism is subject to a phase lag α ∈ T. If α = 0, particles align to the directions of their neighbors. Otherwise, particles perform excessive rotation reminiscent of chiral motion. In a noninteracting regime σ = 0, particles exhibit Brownian motion, modeled by a family of independent Wiener processes W i with a self-diffusion coefficient D ϕ ∈ R + . As a result, particle's dynamics are determined though an interplay of nonlocal alignment and external stochastic forces. As an initial value problem (IVP), we consider this system together with independent and identically distributed initial data p i (0) ∈ Ω for i = 1, . . . , N . Such systems of coupled Langevin equations are a common method to formalize the collective dynamics of interacting particles and they have been extensively investigated on the matter of self-organization phenomena [59] and related phase transitions. The most famous phenomenon is the spontaneous rotational symmetry breaking resulting in the emergence of orientational order, either polar or higher order one. Apart from that, recent attention has been devoted to the revealing of spatially inhomogeneous pattern formation. In regard to the aforementioned model (2), it has been shown previously [36] that it generates a large variety of intriguing spatially inhomogeneous chiral dynamics, the most prominent of which are traveling bands, dense clouds, and vortices, as well as localized chimera states [37].

Continuum limit
When the number of particles becomes large, we look for a continuum description of a particle system. We expect that such a description is more efficient compared to the finite size particle system. The dynamical density functional theory provides a number of ways to derive density functions of particle state variables starting from the Langevin dynamics like (2). Even though it is possible to derive continuum limit equations in terms of a joint many-particle density functions, one usually restricts oneself to consider a one-particle density function. For (2), we have shown [36] using the framework of the Fokker-Planck equation [6] that the continuum description of (2) with macroscopic scaling [34] is given by a one-particle probability density function evolving according to the Vlasov-Fokker-Planck [53,24,44] equation subject to an initial condition: The density function f (r, ϕ, t) quantifies the probability to find a particle at a given position r ∈ U 2 with a given orientation ϕ ∈ T at time t. As before, e(ϕ) = (cos ϕ, sin ϕ) ∈ S 1 is a unit vector in the direction of self-propulsion ϕ. Note that ∇ r = (∂ x , ∂ y ) denotes a spatial gradient. The rotational velocity or torque exerted by a nonlocal neighborhood is found to be The region of nonlocal interaction is now defined as a cylinder The size of the nonlocal neighborhood, which is quantified via |B i | in (2), is now measured as One should keep in mind that in general this neighborhood mass is time-dependent but we usually omit this dependence for notational simplicity. Eq. 5 determines the polarization of a particle flow around a given point (up to a phase shift α). This can be formulated by introducing an interaction kernel K(r, ϕ) = H (r) sin(ϕ + α), which is a product of a Heaviside step function H (r) = H( − r ), which ensures that only a particle flow within the distance is accounted for, and a shifted alignment function. We can therefore express the rotational velocity functional in a more general form as Whereas this functional form of the angular velocity w[f ](r, ϕ, t) does not change the dynamics of the one-particle density function, the convolutional form of the alignment interaction allows for a substantial decrease of temporal complexity of numerical algorithms by means of the discrete Fourier transform [52]. We observe that by writing the diffusion term as ∂ ϕϕ f = ∂ ϕ (f ∂ ϕ log f ), we can combine the last two terms in (3) in the form of a gradient flow as The new functional ξ denotes the potential function of the flow in the angular direction. However, we cannot extend this form to all the right hand side of (3) unless we design a system to be spatially homogeneous.
For the latter case, we use structure preserving numerical strategies developed for gradient flow structures in the construction of a numerical scheme.
One of the transitions we are interested to investigate is the one in terms of the polarization of a particle flow. This is commonly measured using the global polar order parameter defined as Here, the absolute value R gives the aforementioned measure while the phase Θ can be interpreted as a mean direction of the particle flow. If the flow is completely synchronized so that a ϕ-marginal of f is a point mass density, the magnitude equals its maximal value R = 1. If the flow is uniformly distributed so that f = const, the magnitude equals its minimal value R = 0. For any partially synchronized solution with respect to the angular variable ϕ, the order parameter magnitude assumes intermediate values R ∈ (0, 1). Note that in general the right hand side of (8) must be normalized but since we consider f as a probability density, the normalization term equals one. In the following, when referring to the order parameter, we will often refer to its magnitude R since it provides the main structural information about the particle flow. The global polar order parameter (8) provides a global information about the momentum field, which is not enough in a spatially inhomogeneous context. From the continuum PDE (3), we see that it is worthwhile to consider a nonlocalized version of (8) as One can show that in terms of such a nonlocal polar order field, the PDE (3) becomes The presence of the magnitude R in front of the angular flux emphasizes that the rotational rate of change of a particle flow is proportional to the polarization at that point.

Spatially homogeneous formulation
It is straightforward to check that any constant function satisfies (3). But since we are interested in probability density functions in Ω as its solutions, we find In terms of a particle system, this uniform density function corresponds to the chaotic behavior of the system with particles uniformly distributed in U 2 having orientations uniformly distributed in T. One also says that this solution represents a globally disordered or incoherent state. In order to find solutions except for the trivial one, we note that the model (3) admits a major simplification if we assume that the solutions are spatially homogeneous, i.e., f (r, ϕ, t) = f (ϕ, t). Under that assumption, equation (3) simplifies to a (1+1)-dimensional PDE, which we can also consider as the continuum Kuramoto-Sakaguchi model [54] with diffusion: where the nonlocal interaction term (5) becomes a global one w dϕ with the neighborhood mass omitted since |C| = T f (ϕ, t) dϕ = 1. We shall consider (11) with an initial condition: Note that this equation is of a gradient form In the absence of rotations, i.e., α = 0, the free energy associated to equation (11) is given by [60] We can therefore represent (11) in a general gradient flow structure . Moreover, one can show that this energy functional decays along solutions of (11) according to However, for chiral interactions with |α| > 0, the interaction potential is not symmetric and we cannot write down a respective Liapunov functional. Therefore, we will not consider the free energy dissipation of the constructed numerical schemes as in [14]. The Kuramoto model is well-known nowadays and it is often used to model synchronization phenomena in various systems [38,50,5,17]. It was discovered in [39] that the addition of a phase lag parameter α to the model allows one to obtain a new type of solutions, termed chimera states [4], where both synchronized and disordered populations of oscillators coexist. We showed in [37] that by extending an oscillator model to a self-propelled particle model, we can obtain such chimera states with and without spatial homogeneity. It appears that for spatially homogeneous states in the continuum limit with noise, we can find a closed form expression for a corresponding density function.
We note that for spatially homogeneous systems nonlocal (9) and global (8) polar order parameters become equal and one can write the rotational velocity functional (5) in terms of this polar order parameter as First, we look for stationary solutions of (11). Noting that this equation is invariant under phase translations f (ϕ, t) → f (ϕ + ϕ 0 , t) ∀ϕ 0 ∈ T, we may put Θ = 0 without loss of generality. Therefore, we find the following condition for stationary solutions To solve this equation, we integrate it once and look for a solution in the form where the function c(ϕ) is to be determined from the ODE and γ = σR/D ϕ . We find that the solution is where c 1 , c 2 are constants to be determined. The first constant is given by the periodicity f (0) = f (2π), which can be shown to hold if and only if c 2 = 0. Next, from the normalization condition T f (ϕ, t) dϕ = 1, we find that c 1 = (2πI 0 (γ)) −1 , where I 0 (γ) = 1/(2π) T exp(γ cos ϕ)dϕ denotes the modified Bessel function of the first kind [45]. As a result, we have found a nontrivial stationary solution of the form In such a form, the solution is not particularly useful since the density function is recursively contained in the definition of the order parameter (13). However, we are able to determine the latter the other way. If we multiply (15) by cos ϕ and integrate over the domain T, we find that R = 1 2πI 0 (γ) T e γ cos(ϕ+α) e iϕ dϕ = e −iα I 1 (γ) I 0 (γ) .
Since we are looking for a nontrivial solution with R = 0, we conclude that α = 0 and the stationary solution is for which the order parameter magnitude satisfies R = I 1 [γ(R)]/I 0 [γ(R)]. This result tells us that stationary states are possible only in the absence of the phase lag, i.e., in the form of irrotational motion. We refer to [14] for the detailed analysis of such states. Our next step is to investigate solutions to (4) in the presence of a nonzero phase lag α. By noting that in this case the particle flow moves uniformly either to the left or to the right depending on the sign of α, we are looking for solutions in the form of a traveling wave. We introduce an ansatz f (ϕ, where v is the speed of the traveling wave, which is unknown. After this substitution, the compatibility condition (14) for stationary solutions becomes Integrating it with respect to ω yields where c 1 ∈ R is some constant. Solving this equation, we find the solution to be where c 2 ∈ R is some constant. One of the constants is fixed due to the periodicity constraint, i.e., g(0) = g(2π). Namely, this implies This subsequently implies Next, due to the normalization condition T g(ω) dω = 1, we put c 2 = c 0 as a normalization constant and find [29] This defines a profile of the traveling wave solution f (ϕ, t) = g(ϕ − vt) = g(ω). One may apply the ansatz backwards in order to obtain the complete form f (ϕ, t) of the solution of (11) but for the subsequent numerical analysis, we will use its profile solely. We refer to [36] for other forms of this traveling wave solution.
Proof. From the definition of the spatially homogeneous polar order parameter (13), we can write By differentiating the first equation, we finḋ

Phase space discretization
Our main goal is to develop a finite volume scheme that reproduces a correct behavior of nonlocally interacting particle flow governed by a 3+1 dimensional integro-differential PDE (3). Given that we have much more background on a spatially homogeneous system (11), we first introduce a one-dimensional finite volume scheme, i.e., for density functions of an angular variable. We next proceed to the description of a complete three-dimensional scheme that is applied to density functions of spatial and angular variables. For both proposed methods, we prove mass and positivity preservation as well as derive CFL conditions on their stability.
The approach we are going to pursue is the following. First, we perform a phase-space discretization of a PDE of interest, thereby deriving a semidiscrete system of ODE for finite volume cells. We derive the set of equations on uniform meshes but the generalization to nonuniform ones is straightforward. Afterwards, by noting that the dynamics of a velocity field in spatial and angular directions qualitatively differs, we attempt a dimensionality splitting technique in order to effectively cope with the dynamics changes due to spatial and angular fluxes. As a result, we obtain a FVM that is second order accurate both in time and in phase-space variables.

One-dimensional scheme for spatially homogeneous PDEs
In this section, we develop a finite volume scheme for continuum limit PDEs under the assumption of spatial homogeneity, i.e., for equation of the form (11). Let T L = {0, . . . , L − 1} denote a discreet onedimensional torus with L points. We divide a domain T into finite volume cells C k = [ϕ k− 1 2 , ϕ k+ 1 2 ], k ∈ T L of a uniform length ∆ϕ = 2π/L with the center of a cell ϕ k = k∆ϕ, which correspond to a site k in the torus T L . Note that since the space is periodic, we have ϕ k = ϕ k+L , k ∈ T L .
We define the cell averages [12] f k : The cell averages f k are functions of time but for the sake of compactness, we will henceforth omit the explicit time dependence of the computed quantities. The semidiscrete finite volume scheme is obtained by integrating the PDE (11) over each cell C k , k ∈ T L . It is consequently formulated as the following system of ODEs for the cell averages: where F ϕ k± 1 2 denote angular fluxes. Note that the right hand side of this expression is a second order centered difference of an original flux. In order to find the numerical approximations of the above fluxes at cell interfaces, we need to be able to compute the corresponding values of a solution itself as well as a velocity field. In this paper, we adopt a piecewise linear reconstruction of the numerical solution f (ϕ, t) at each time point. Saying that, we represent a density function in each cell as a first order polynomial as where (∂ ϕ f ) k is a cell average of a partial derivative with respect to ϕ, which is to be determined for this reconstruction to work. The knowledge of values of a solution at neighboring cell centers allows us to approximate the slopes (∂ ϕ f ) k using a second order centered difference method: Unfortunately, it might occur that this slope approximation might lead to negative values of a reconstructed numerical solution (18), which we intent to circumvent. For such cases, we recalculate the slope by imposing a slope limiter that keeps reconstructed values nonnegative. In this paper, we chose to use a generalized minmod limiter Note that the values, which are corrected with this slope limiter, are generally of first order. However, in all numerical tests we present in this paper, it is practically not imposed and the numerical scheme stays effectively of second order in ∆ϕ. At this point, the piecewise linear reconstruction is defined and we can apply (18) in the calculation of numerical fluxes, required in (17). First, we need to know the values of a solution at each cell interface.
They are computed as where f (ϕ k+ 1 2 −0) and f (ϕ k− 1 2 +0) denote function values at cell interfaces ϕ k+ 1 2 and ϕ k− 1 2 from inside a cell C k , respectively. We use the cell interface values to define the numerical fluxes in (17) as upwind fluxes as where positive and negative parts of angular velocities are denoted by The exact velocities themselves are given in the PDE (11) but need to be numerically approximated at cell interfaces. We note that the PDE is of a gradient flow structure, i.e., we can write for some potential function ξ, as we showed in (7), and we use this fact to define velocities at cell interfaces using a second order centered difference method as The velocity potential for a homogeneous system is defined as therefore, its numerical approximation proceeds as follows As a result, we have obtained the following representation of an approximated velocity potential, which is to be used in (23), as where we have used the fact that since f is a probability density function, its piecewise linear reconstruction yields Tf (ϕ, t)dϕ = n∈T L f n ∆ϕ. We note that the above approximation of the velocity potential is exact in ∆ϕ given a piecewise linear reconstruction of a density function.
Theorem 3.1. Consider the IVP (11)- (12) with periodic boundaries and the semidiscrete FVM (17) with a positivity-preserving piecewise linear reconstruction (18). Assume that the system of ODEs (17) is discretized by the forward Euler method or by a higher-order strong stability preserving (SSP) ODE solver, whose time step can be expressed as a convex combination of several forward Euler steps. Then, computed cell averages remain nonnegative f k (t) ≥ 0 ∀k ∈ T L ∀t > 0, provided that the following CFL condition is satisfied: and the velocities at cell interfaces are defined in (23).
Proof. According to the forward Euler method, we discretize (17) as We note that we can express cell averages of a solution as a linear combination of corresponding values at cell interfaces, defined in (20), as Using this fact and expressing numerical fluxes as upwind fluxes introduced in (21), we have

Now we group the terms according to cell interface values and find
The last two terms are always nonnegative. To guarantee positivity preservation of f k (t + ∆t), we must require that the values in parentheses of the first two terms remain nonnegative. This yields the following conditions: the combination of which gives the desired CFL condition.
Remark 3.1.1. We indicate that it was proved in [12] that the one-dimensional numerical scheme considered so far preserves the conservation of mass and provides the decay of discrete free energy for systems with symmetric interaction potential. It also has a very good property: numerical steady states are characterized by ξ k = constant for all k, as in the continuum setting, due to (23).
Numerical studies of spatially homogeneous PDEs (11) using the finite volume scheme of this section are conducted in Sections 4.2-4.5. In the following, we generalize the scheme for general spatially inhomogeneous three-dimensional PDEs (3).

Three-dimensional scheme for spatially inhomogeneous PDEs
We start from the discretization of a phase space Ω into finite volume cells. Dimensions corresponding to x, y, and ϕ variables are divided into N , M , and L cells, respectively. Linear sizes of cells are ∆x = 1/N , ∆y = 1/M , and ∆ϕ = 2π/L. Let U N , U M , and T L denote discreet one-dimensional tori with N, M, L ∈ N points, respectively, i.e., We define a uniform grid consisting of cells The discretization of Ω consists of three-dimensional cells which can be enumerated with a threedimensional torus as We define cell averages [12] As before, for the sake of compactness, we will omit the dependence of most computed quantities on time t henceforth.
The semidiscrete finite volume scheme for a three-dimensional system is obtained by integrating the PDE (3) over each cell C i,j,k , (i, j, k) ∈ Ω N,M,L and is formulated by the following system of ODEs for for i ∈ U N , j ∈ U M , and k ∈ T L . In order to define the above fluxes, we extend the same piecewise linear reconstruction method from the previous section. The numerical solution in each cell C i,j,k is thus approximated as a polynomial To be able to use this representation, we need to find each slope (∂ x f ) i,j,k , (∂ y f ) i,j,k , and (∂ ϕ ) i,j,k . To ensure that the solution is second-order accurate, we define the slopes using the centered difference approximations It might occur that a reconstructed solution becomes negative in some cell C i,j,k . In such cases, we recalculate a corresponding slope using a slope limiter. In this paper, we use a generalized minmod limiter (19), whose application yields We note again that the values, which are corrected with such a limiter, become of first order in a corresponding dimension. But since the number of such points is usually small compared to the total number of grid points, the overall order of the scheme is effectively not reduced. In the numerical tests of this paper, these slope limiters have practically not been triggered. Now that the piecewise linear reconstruction (25) is completely determined, we are able to compute solution values at each cell interface the following way: We compute all fluxes in the semidiscrete system of ODEs (24) as upwind fluxes as where u ± 2 ,k , and w ± i,j,k+ 1 2 denote positive and negative parts of velocities defined as before according to (22). The definition of upwind fluxes required the knowledge of velocities at cell interfaces. For advection in spatial directions, they are directly obtained from the PDE (3) by direct substitution of grid points: For the angular direction, we again use the gradient flow structure of the angular subflow and consider the last two terms in the PDE (3) as defined by a velocity potential ξ[f ](r, ϕ, t) (7). One can show, that for a spatially nonhomogenous system, it reads We thus determine the angular velocity at cell interfaces from its potential using a second order centered difference scheme Now, we need to perform the discretization of the potential (29) using the piecewise linear reconstruction of the solution (25). First, the numerator for ξ i,j,k is calculated as f l,m,n ∆x∆y∆ϕ.
As a result, we find the discretized velocity potential ξ i,j,k , which is to be used in (30), to read which is exact in ∆x, ∆y, and ∆ϕ.
Theorem 3.2. Consider the IVP (3)-(4) with periodic boundaries and the semidiscrete FVM (24) with a positivity-preserving piecewise linear reconstruction (25). Assume that the system of ODEs (24) is discretized by the forward Euler method or by a higher-order SSP ODE solver, whose time step can be expressed as a convex combination of several forward Euler steps. Then, computed cell averages remain nonnegative f i,j,k ≥ 0 ∀i ∈ U N , ∀j ∈ U M , ∀k ∈ T L , ∀t > 0, provided that the following CFL condition is satisfied: with the coefficients a = max and the velocities at cell interfaces are defined in (28), (30).
Moreover, the mass of the discretized system is conserved, i.e., Proof. The proof of this theorem follows the same lines as in Theorem 3.1. We only comment that in this case, one should express cell averages of a solution f i,j,k as a linear combination of corresponding values at cell interfaces, defined in (26), as The proof of the conservation of mass follows the same reasoning as for the one-dimensional case, which can be found in [12].
Numerical studies of spatially inhomogeneous PDEs (3) using the finite volume scheme of this section are conducted in Sections 4.7-4.8. Up to now, we have presented the second order discretization of the phase space of the problem keeping the time domain continuous. This way, we have formulated the problems of solving PDEs (3) and (11) as the problems of solving systems of ODEs (24) and (17), respectively. In the next section, we consider approaches to perform time discretization so as to keep the FVM of second order in time as well.

Dimensionality splitting
This section considers further discretization approaches for three-dimensional PDEs only. As a result of the phase space discretization from the previous section, we reformulate the IVP for PDEs (3)-(4) as the IVP for the system of ODEs (24), which we state here for convenience: for (i, j, k) ∈ Ω N,M,L and the fluxes are as defined in the previous section. From the PDE (3) itself and from the derivation of the system of ODEs (31), we have seen that the fluxes qualitatively differ for spatial and angular dimensions. Therefore, it might become unreasonable to tackle all of them at once. We will employ this fact in the further construction of our FVM. The class of methods that allow us to separate dynamics of an autonomous system of ODEs into several subsystems is known as splitting methods [30]. In this paper, we are interested in such splitting methods that lead to the second order accuracy in time. These splitting methods have been classically used in simulations for kinetic equations in plasma physics, see [35,25] for instance. First, we note the following. We can reenumerate the three-dimensional grid of the problem, which is enumerated with three-dimensional indices from Ω N,M,L , with one-dimensional indices ranging from 0 to N M L. In other words, the three-dimensional set of indices Ω N,M,L can be bijectively mapped into a one-dimensional set with N M L indices. This allows us to express the above model shortly as We note that we can separate the vector field on the right hand side as d dtf where F [1] and F [2] are defined by spatial and angular fluxes from (31), respectively. Now, according to the splitting procedure, instead of the problem (32), we consider two subproblems: Let φ  and build the composition of Φ and Φ * with halved step sizes, we obtain by the theorem of the composition of methods [30], a new method ∆t/2 , which is known as the Strang (Marchuk) splitting. The Strang splitting formula is of order 2. However, for that technique to work, we must be able to integrate two subproblems (33) exactly, which is usually not the case. Therefore, our aim is to obtain a second order method which consists of the composition of approximate direct methods only.
Proposition 3.1. Consider the composition of numerical methods If the methods Φ [1] and Φ [2] are of second order at least, the resulting composition method Ψ is of second order.
Proof. Let us denote approximate solutions to subproblems (33) using numerical flows asf 1 (∆t/2) = Φ ∆t/2 (f 3 (0)) if f 3 (0) =f 2 (∆t). By performing Taylor expansion of each solution up to second order around ∆t = 0 and expressing them in terms of a respective previous solution, we obtain the second order Taylor expansion of the exact flow generated by (32).
Note that one may start by considering the approximate flows Φ ∆t from the very beginning and combine them as Φ ∆t = Φ ∆t to obtain a first order approximation to (32). By performing the composition of this method with its adjoint with halved step sizes [30] and replacing all adjoint flows with the direct ones, one obtains Ψ ∆t = Φ ∆t/2 , which can also be shown to be second order accurate in time, provided each of the submethods is of second order too. However, since (34) requires fewer time steppers, we choose it for all the subsequently reported results. For such an FVM, we state similar results about prositivity preservation and mass conservation, as in previous sections.
where the numerical fluxes F x i± 1 2 ,j,k , F y i,j± 1 2 ,k , F ϕ i,j,k± 1 2 are defined as in (27) using a positivity-preserving piecewise linear reconstruction (25). Assume that each system of ODEs (33) is discretized by a second order method whose time step can be expressed as a convex combination of several forward Euler steps. Then, the computed cell averages remain nonnegative f i,j,k ≥ 0, (i, j, k) ∈ Ω N,M,L for all t ≥ 0, provided that the following CFL condition is satisfied: and the velocities at cell interfaces are defined in (28), (30). Moreover, the mass of the discretized system is conserved, i.e., d dt Ω f (r, ϕ, t) drdϕ = 0.
Proof. The proof of this theorem follows the same lines as in Theorem 3.1 but it should be applied to each subsystem (33) separately. By expressing cell averages of a solution f i,j,k of the first subsystem as a linear combination of corresponding values at cell interfaces in angular direction, defined in (26), as and a solution of the second subsytem as a linear combination of corresponding values at cell interfaces in spatial directions as the result follows. The proof of the conservation of mass follows the same reasoning as for the one-dimensional case, which can be found in [12].
We implemented the presented finite volume schemes in C++. To be able to perform numerical analysis of the schemes with meaningful phase space discretizations, we parallelized algorithms using the message passing interface (MPI) standard [27]. Our implementation can be found under [2].

Numerical tests
In this section, we demonstrate the performance of the developed numerical scheme. As it has been mentioned, the assumption of spatial homogeneity of solutions allows us to greatly simplify the analysis of the system of interest as well as gain theoretical insight on the behavior of its solutions. Therefore, we first conduct numerical experiments under the assumption of spatial homogeneity in Sections 4.2-4.6, including the study of phase transitions of traveling wave solutions in particular. The phase transitions are quantified with respect to the global polar order parameter that measures the degree of polarization in a particle flow. Next, we test the numerical scheme in a general setup, where nonstationary spatially inhomogeneous solutions are expected to exist, including the study of related phase transitions as well. These studies are presented in Sections 4.7-4.8.

Error norms
In the following, we will examine the accuracy of numerical schemes in both spatially homogeneous and inhomogeneous setups. In the former, analytic solutions are known while in the latter, they are not. Therefore, we introduce different norms for different cases. If we know an exact solution, we will use the following norms to quantify convergence errors [55]: wheref h is a numerical solution with the reconstruction defined by (18) or (25) and phase space discretization h, f is an exact solution, C i,j,k = [x i− 1 2 , x i+ 1 2 ]×[y j− 1 2 , y j+ 1 2 ]×[ϕ k− 1 2 , ϕ k+ 1 2 ], i = 0, . . . , N −1, j = 0, . . . , M − 1, k = 0, . . . , L − 1 is a cell on a uniform grid with N M L points as defined previously. The integrals in the above expressions are computed using Gauss-Legendre quadrature [52]. Note that in case we use quasiuniform initial conditions, we need to align both solutions in a proper way. In the reported results, we shift an exact solution such that its first moment coincides with the first moment of the reconstructed solution.
In cases where we do not have an exact solution, which is the case when we retrieve spatially inhomogeneous patterns, we first compute a reference solution with the finest discretization h 1 and compare the rest of the solutions with cruder discretizations h 2 to it using where the grid size is chosen as a least common denominator in each dimension, i.e. C i,j,k , i = 0, . . . , lcd(N 1 , N 2 )− 1, j = 0, . . . , lcd(M 1 , M 2 ) − 1, k = 0, . . . , lcd(L 1 , L 2 ) − 1. If one needs to compute an error for one and two dimensional domains, it is done straightforwardly by omitting two or one dimensions, respectively, in the above definitions.

Stationary phase synchronization (1D)
We start the inspection of performance of constructed numerical schemes by first analyzing the simpler spatially homogeneous systems, whose time evolution is governed by PDEs (11). It is known that for sufficiently high diffusion levels D ϕ (or equivalently for small coupling coefficients σ), the asymptotic solution consists of chaotically moving particles, whose distribution is given by a uniform density function (10). For sufficiently low diffusion levels (or large coupling coefficients), particles self-organize into spatially homogeneous polar groups, which are stationary in the absence of a phase lag, i.e., for α = 0, or rotate with constant frequency for α = 0. In this section, we investigate how the numerical scheme performs in the former case. Namely, we consider a spatially homogeneous version of the continuum limit equation with α = 0, also known as the continuum Kuramoto model [14] with diffusion: where the angular velocity induced by particles' interactions is w The initial condition f (ϕ, 0) in the form of a trigonometric series is used to model an irregular but sufficiently smooth function, which is required by the numerical scheme. We shall refer to such initial conditions as quasirandom initial conditions in subsequent discussions. The series coefficients are chosen in such a way that f (ϕ, 0) is nonnegative and properly normalized, i.e. a 0 = 1 2π , a k , b k ∼ U(−ε, ε), k = 1, . . . , K, K ∈ N. Note that we cannot use a uniform probability density function (10) as an initial condition since it is already a solution to the problem (37). We also remark that the normalization of the density function is not generally required by the scheme but continuum limit PDEs, we consider in this paper, describe the behavior of probability density functions.
It is well known that the problem (37), i.e. the continuum Kuramoto model for identical oscillators with diffusion, exhibits a second-order phase transition with respect to either coupling strength σ or diffusion level D ϕ . The phase transition if of second order and occurs at D ϕ = σ 2 . Its numerical investigation was already described in detail in [14], therefore, we do not consider it here. Instead, we only test how our finite volume scheme performs on the solutions of (37) for parameters from the region of stability of a polar order solution. In this case, this solution is a von Misés density function where Θ ∈ T is the average direction of a particle flow, whose value depends on initial conditions, and I 0 is the modified Bessel function of the first kind. Fig. 1(a) illustrates that a numerical solution approximates the exact one very well. Note that since the numerical solution was obtained from quasiuniform initial conditions according to (37), it is manually centered so that Θ = π, for comparison reasons. The numerical solutions are taken at t = 200, when all of them has converged to steady states. However, such a long time is not required for all presented solutions. The time the system takes to converge to a steady state depends on the value of a diffusion coefficient D ϕ (or reversely the coupling strength σ). The closer the value to the order-disorder transition point D ϕ = σ 2 , the longer the time is. This is a well known bottleneck effect near bifurcation points. For the solution with D ϕ = 0.4σ, it takes around t = 200 simulation time units to converge.
The results of error convergence are presented in Fig. 1(b). One can see that in the current setup, the scheme is second order accurate as is guaranteed by its construction. Numerical solutions were compared to the aforementioned von Misés density function (38) using error norms defined by (35). Initial conditions were again quasiuniform, the time step was ∆t = 10 −5 , and the errors were computed at t = 100, when the steady state has been reached.

Nonstationary phase synchronization (1D)
Next, we keep the assumption of spatial homogeneity but assume α = 0. If the phase lag is added to the particle alignment interaction, a particle flow starts to rotate. Its continuum description in terms of a probability density function becomes skewed and assumes a traveling wave form. A slight generalization to the previous example leads to the following continuum Kuramoto-Sakaguchi model [54] with diffusion: where w[f ](ϕ, t) = σ T f (ϕ , t) sin(ϕ − ϕ − α)dϕ and the choice of an initial condition follows the same considerations as in the previous example. This problem has again two solutions, a uniform density function (10) and a skewed unimodal density function (16) where c 0 ∈ R is a normalization constant. The stability of solutions to the problem (39) now depends on the values of a phase lag parameter α, diffusion level D ϕ , and coupling strength σ. With respect to these parameters, the system exhibits a second order phase transition, which we will investigate later (see Section 4.6). The phase transition occurs at D ϕ = σ 2 cos α. Here, we illustrate the performance of the scheme for parameter values from the region of stability of the skewed density function (40), i.e. for D ϕ < σ 2 cos α. Fig. 2(a) illustrates that the numerical solution approximates the exact one well. Note that Fig. 2(a) shows profiles of actual solutions which have been manually centered for comparison reasons. The numerical solutions were taken at t = 1000, which is the time the solution with α = 1, D ϕ = 0.25, σ = 1 takes to converge to the traveling wave form (40). This is again because of the critical slowing-down close to the phase transition line (see the discussion in the previous section). In contrast to the zero phase lag case, the phase transition occurs for smaller diffusion levels for fixed α. Moreover, the time the system takes to converge increases exponentially with α even when diffusion is absent [37].
We illustrate the error convergence for parameters from the same stability region (cf. Fig. 2(b)). One can see that the errors are again of second order as expected. The norms were computed according to Eqs. (35) with Eq. (40) as a reference solution. The time step was chosen ∆t = 10 −5 . The norms were computed at t = 100, when the system had converged to a traveling wave solution.

Stationary phase synchronization (3D)
Spatially homogeneous PDE (11) was obtained from the original one (3) under the assumption of spatial homogeneity of the system. On one hand, it allowed us to derive analytical results to understand model's behavior but on the other hand, this does not correspond to the original particle flow. We next consider the complete (3+1)-dimensional equation (3) but start from the analysis in a parameter region where spatially homogeneous solutions are stable. Again, we first assume a simpler case where the phase lag is not taken into account, i.e. α = 0 and consider the following IVP where e(ϕ) = (cos ϕ, sin ϕ) ∈ S 1 is a unit velocity vector and the angular torque is defined as w[f ](r, ϕ, t) = σ |C(r)| C(r) f (r , ϕ , t) sin(ϕ − ϕ)dr dϕ . Coefficients c 0 , c nml ∈ R in the initial condition are chosen such that the density function f (r, ϕ, 0) is nonnegative and normalized, and the normalization term |C(r)| is defined by Eq. (6). The shifts of the arguments are chosen at random α nml , β nml , γ nml ∼ U(0, 2π) in order to approximate a fluctuating density field. One can consider such an initial condition as a generalization of initial conditions from spatially homogeneous examples (37), (39) to the three-dimensional case.
From the linear stability analysis [36], we know that for α = 0 the problem (41) has the same set of solutions as the spatially homogeneous problem (37), i.e. a uniform density function (10) and the von Misés density function (38). It appears that these solutions are stable against spatially inhomogeneous perturbations in the absence of phase lag. The reason why spatial patterns, such as traveling bands, do not develop for parameter values close to the order-disorder transition line is because the alignment interaction is normalized in w[f ](r, ϕ, t). We remark that in case the normalization term |C| is removed, one actually observes the emergence of such spatial patterns, as is well known for active matter systems with polar interactions [41,43]. An exemplary solution to the problem (41) is presented in Fig. 3(a). One can see that the solution is indeed homogeneous with respect to x and y but has a unimodal symmetric profile in ϕ. Its projection onto the ϕ-axis is qualitatively similar to the one in Fig. 1(a).
The analysis of solutions for different grid sizes shows the second order convergence in terms of ∆ϕ (cf. Fig. 3(b)). Since the steady state is eventually spatially homogeneous, we cannot test the error convergence in terms of ∆x and ∆y. The errors were computed using Eqs. (35) with Eq. (38) as a spatially homogeneous reference solution. The time step was chosen ∆t = 5 · 10 −4 . The errors were computed at t = 100, when the numerical solution had converged to the steady state. The model parameters were taken from the region of

Nonstationary phase synchronization (3D)
We now proceed by allowing nonzero phase lag α = 0 in a general three-dimensional system (3). We know that in addition to a uniform disordered (10) and spatially homogeneous ordered motion (40), spatially inhomogeneous solutions emerge for sufficiently high values of α. Before we proceed to the analysis of such solutions, we would look into performance of our numerical scheme in a parameter region, where a spatially homogeneous skewed unimodal density function (40) is stable against spatially inhomogeneous perturbations in order to be consistent with the previous development. We consider the following IVP For α sufficiently small so as to guarantee the stability of spatially homogeneous densities (40), an exemplary numerical solution is illustrated in Fig. 4(a). One can see that it is indeed homogeneous in x and y but has a characteristic skewed shape in ϕ (the blue region is more pronounced above the plane than below). Its projection onto the ϕ-axis gives a qualitatively similar solitary wave as the one in Fig. 2(a). The wave moves transversal to the ϕ axis with some constant velocity v, the sign of which depends conversely on α. Thus, it has the form of a plane wave in Ω. Such a form of the solution corresponds to a partially synchronized steadily rotating particle flow, which was investigated and termed a nonlocalized self-propelled chimera state in [37].
As in the previous example, because of spatial homogeneity of the solution, we are able to determine the order of error convergence versus ∆ϕ only (cf. Fig. 4(b)). The errors were computed with respect to a spatially homogeneous plane wave solution whose ϕ-profile is given by Eq. (40) using Eqs. (35). The time step was chosen ∆t = 5 · 10 −4 . The errors were computed at t = 200, when numerical solutions converged to the plane wave form.

Phase transitions of spatially homogeneous solutions
In this section, we analyze how phase transitions of spatially homogeneous solutions are captured by the finite volume scheme. Moreover, we concentrate on the nonstationary case with nonzero phase lag and refer the reader to [14] for the numerical studies of phase transitions of steady states in the continuum Kuramoto model. The phase transitions between spatially homogeneous polar order and disordered motion are commonly quantified with respect to the polar order parameter defined in (13), i.e., it is an average orientation on the unit circle with respect to a given density function. Since we know that the continuum Kuramoto-Sakaguchi equation (11) has a solution of the traveling wave form, it is enough to study its profile g(ω) = f (ϕ − vt, 0) = f (ϕ, t) in the analysis of related phase transitions. The polar order parameter can thus be expressed as where the average direction Θ can be shifted to the origin without loss of generality due to the translation invariance of Eq. (11), i.e. Θ ≡ 0. By expanding the above system into the real and imaginary parts, we see that the order parameter must satisfy the following set of self-consistent equations (SCEs) This system does not have an analytical solution but we can solve it numerically for polarization R and velocity v. The numerical solution of this system in terms of the model parameters α and D ϕ is presented in Figs. 5(a) and (b). Given the above result, we are able to investigate phase transitions related to a spatially homogeneous rotating system (42). For each parameter set, we solve it starting from quasiuniform initial conditions, described earlier. We discretize the domain T into L = 256 points and perform all computations of this section with the time step ∆t = 10 −4 until t = 10 4 . The results are presented in Figs. 5(c) and (e) for the order parameter magnitude R, and in Figs. 5(d) and (f) for the group velocity v. Because either the coupling strength σ or the diffusion coefficient D ϕ may be eliminated by an appropriate rescaling of time in Eq. (42), we fix σ = 1 without loss of generality and investigate the behavior in terms of α and D ϕ . Since a phase transition may happen by changing either the diffusion level D ϕ or the phase lag α, we look into both possibilities. First, we fix α = 0.5, 1.0, 1.5 and investigate the behavior of R = R(D ϕ ) (cf. Fig. 5(c)) and v = v(D ϕ ) (cf. Fig. 5(d)). As we expect, the phase transition occurs at D ϕ = σ 2 cos α. The black line additionally shows how the order parameter behaves next to the transition line, according to which is known from the hydrodynamic description of the particle model (2) and derived in [36]. The grid size for D ϕ was chosen 0.0025. For comparison reasons, we show the results obtained with the FVM and those of SCEs (43). Next, we fix D ϕ = 0.1, 0.25, 0.4 and investigate the behavior of R = R(α) (cf. Fig. 5

(e))
and v = v(α) (cf. Fig. 5(f)). Again, as expected, the phase transition occurs at α = arccos A few remarks on phase transitions between spatially homogeneous states are in order. We quantify them in terms of a polar order parameter magnitude R. Therefore, under the assumption of spatial homogeneity, we observe a second order transition from partially synchronized to disordered state (cf. Figs. 5(c) and (e)). However, it is known that for polar active matter systems, the transition from uniform disordered to ordered motion is separated by a region with density-segregated patterns such as traveling bands. It has been established that for the Boltzmann equation for self-propelled particle systems, the respective transitions may be of first and second order depending on a system size [57]. As we mentioned earlier, we do not observe such a spatially inhomogeneous regime close to the order-disorder transition line because nonlocal alignment interactions in (3) are normalized, which makes spatially homogeneous solutions more stable against spatially inhomogeneous perturbations [36].

Spatially inhomogeneous solutions
We have seen that the presented finite volume scheme is capable of reproducing correct behavior under the assumption of spatial homogeneity. Therefore, we proceed to the general PDE (3) in a parameter region where spatially inhomogeneous patterns occur. We consider the same IVP (42) from the previous section. From the linear stability analysis of this PDE from the point of view of kinetic theory [36], we know that there exists a region in the parameter space of /v 0 , α, and D ϕ where the spatially homogeneous skewed unimodal density function (16) becomes unstable against spatially dependent perturbations and one observes numerous spatially inhomogeneous patterns.
In this paper, we would like to concentrate on the analysis of one of such patterns, which we refer to as a localized chimera state [37]. From the point of view of the particle model (2), such a state simultaneously consists of two rather distinct groups of particles. The first group forms a subset of particles that gather into a circular polarized cloud, which rotates in the background of the rest of chaotically moving particles. In the continuum limit, this state is characterized by the formation of a high density skewed ellipsoidal region in Ω (cf. Fig. 6(a)) that follows a helical path. We can obtain the aforementioned high density circular cloud as a projection of such a solution into spatial coordinates (cf. Fig. 6 The spatial profile has a form of a bivariate solitary wave which possesses a characteristic front-end asymmetry. This can be observed if we look at the solution profile along the line (cf. Fig. 7(b), a white dashed line), centered at the point of maximal density and directed according to the velocity field at this point. We look for the point of maximal spatial density as The velocity field can be retrieved from the momentum field, in turn, obtained as where e(ϕ) = (cos ϕ, sin ϕ) ∈ S 1 . Note that the momentum field such defied is isomorphic to the global polar order parameter, we used earlier (8). As a result, the line can be parameterized as where ϕ max = arg(u(r max )). Using the piecewise linear reconstruction of the density function (25), we can find the approximate values of the density function at any point on the line (cf. Fig. 7(b)). Moreover, for solutions, where the radius of rotation of the localized cluster is sufficiently large so that the cluster does not rotate around a fixed point, the transverse profile of the spatial projection f (r, t) has a symmetric form. We remark that Eqs. (44), (45) constitute a hydrodynamic description of the kinetic PDE (3) where polar order is expected to emerge, and is often used to get analytical insights into the dynamics. However, the known drawback of this approach is that it is limited to the regimes close to equilibrium. As opposed to that, our FVMs are applicable to any region in the parameter space, which we will employ in the following.
In a general spatially inhomogeneous setup, we are able to calculate error convergence in both angular and spatial variables. Because we do not know any exact solution with spatial dependence, we use Eqs. (36) to compute the norms. The reference solution was the one with N × M × L = 40 × 40 × 1024 grid points. First, we fix N = 40, M = 40 and vary the angular grid size ∆ϕ. We observe the second order convergence for ∆ϕ sufficiently small whereas it approaches the first order for largest grid sizes (cf. Fig. 6(c)). The reason for the first order behavior is that such discretizations cannot capture high density gradients in ϕ so that the numerical error is accumulated rather fast. Subsequently, since the dynamics in angular and spatial dimensions are coupled, this results in solutions, diffused away from the correct dynamics. The time step was chosen ∆t = 5 · 10 −4 . Next, we fix L = 128 and vary the spatial grid size ∆x, ∆y. We again observe the second order error convergence versus ∆r = (∆x) 2 + (∆y) 2 even for quite small grid sizes with N × M = 20 × 20 points (cf. Fig. 6(d)). We explain such robustness of the results by the fact that the solutions were computed for quite large = 0.3, resulting in long range interactions. The time step was chosen ∆t = 5 · 10 −4 . The reference solution was the one with N × M × L = 160 × 160 × 128 grid points.

Phase transitions of spatially inhomogeneous solutions
In Section 4.6, we found that in spatially homogeneous systems, the transition between polar order and disorder is of second order. Now that we have an established protocol to generate spatially inhomogeneous solutions, we are interested to learn about their related phase transitions. By choosing an appropriate scale for the microscopic particle velocity v 0 and interparticle interaction radius , we find a phase diagram in the (α, D ϕ )-parameter space where three distinct solutions are observed, i.e., a spatially homogeneous disordered motion (SHDM), a spatially homogeneous ordered motion (SHOM), and a spatially inhomogeneous motion (SNM) in the form of a localized chimera state (cf. Fig. 7(a)). The figure was obtained from the kinetic linear stability analysis of Eq. (16), performed in detail in [36]. We fix v 0 = 0.25, σ = 1, = 0.3. Phase transitions discussed earlier correspond to the transition between SHDM and SHOM, with the bifurcation occurring at D ϕ = σ 2 cos α (a black solid line). In the gray region, starting from quasirandom initial conditions, one observes the formation of spatial structures. The formation happens in two stages. First, the system smooths out initial spatial perturbations but gradually polarizes until the occurrence of a skewed angular profile similar to the one in Fig. 2(a). Second, the remaining spatial perturbations start to act on the solution and accumulate eventually into a bivariate unimodal shape. We remark that due to finite numerical precision, the second step might not be triggered for any quasirandom initial conditions given that spatial variations become of order of round-off error O(10 −18 ). To circumvent that, one might either look for such initial conditions that preserve enough spatial perturbations at the time point of maximal synchronization or use multiprecision arithmetic libraries, like the one we used in our implementation.
As one can see in Fig. 7(a), the new phase transitions should occur between SHOM and SNM by varying either the diffusion coefficient D ϕ or the phase lag parameter α. We inspect each route separately. Before we do that, we need to establish an appropriate order parameter to measure the level of spatial localization induced by a PDE solution as well as be able to detect changes in spatial variation of solutions upon varying model parameters. First, in a similar way as we might consider the global polar order parameter R(t) (8) as a measure of angular localization of orientation vectors e(ϕ) belonging to S 1 , manifested in a momentum field definition (45), we define an order parameter that measures the level of spatial localization of elements belonging to S 1 × S 1 in the following way: One observes a second order phase transition between SNM and SHOM with bifurcation at Dϕ ≈ 0.0125. Black dashed lines denote unstable branches and have been computed from the system of SCEs (43). (e,f) Localization P and polar R order parameters versus α, respectively. One observes two types of phase transitions, namely, the first order one, accompanied with a hysteresis loop, on the path OA 1 with bifurcation points α ≈ 1.33, 1.44 and the second order one on the path OA 2 with bifurcation at α ≈ 1.515. Black dashed lines denote unstable branches and have been obtained from (43). Black dotted lines are drawn "by hand" in place of unknown unstable branches. Colored lines in (c,d,e,f) indicate directions of bifurcation paths.
This parameter provides the following information. For systems with all the probability mass compressed in one point, i.e., for point measures, the spatial localization is most pronounced and the magnitude of the order parameter attains its maximal value P = 1. In the opposite case, for systems with uniform distribution of matter, no spatial localization is observed and the order parameter magnitude attains its minimal value P = 0. For partial localization inside a particle flow, we therefore have P ∈ (0, 1). The phase Ψ is irrelevant to our purposes. Second, to detect changes in spatial structure of solutions while changing model parameters, we introduce the following maximum absolute spatial deviation measure [57]: where spatially averaged solutions are computed as For SHOM, this measure attains values of order O(10 −14 ), when the magnitude of spatial variations is of a round-off error for double precision floating point values.
We now describe the transitions between SNM and SHOM. We start with a parameter point well inside a region where SNM is a stable solution, i.e. O = (α, D ϕ ) = (1.45, 0.0075) (cf. Fig. 8(b)), and proceed in a continuation-like manner. First, we fix the phase lag parameter α and increase the diffusion level D ϕ ∈ [0.0075, 0.02] with a parameter step size ∆D ϕ = 0.000625 (cf. Fig. 7(a), a vertical path OB 2 ). Starting from quasirandom initial conditions (42), we let the system to converge to a solitary wave form of a localized chimera state and take this solution as an initial condition for a subsequent computation. Then, we change the diffusion level, take as a new initial condition the final solution from a previous parameter, and let the system equilibrate for T = 100 time units with ∆t = 5 · 10 −3 . Afterwards, we accumulate values of δ r (t) (48) and continue integration until spatial deviations cease to fluctuate with dδ r (t)/dt < 5 · 10 −5 . We quantified the rate of change dδ r (t)/dt as a linear fit to a set of values δ r (t) over the last 50 time units. The result of this procedure is shown in Fig. 7(c). As we see, the spatial order parameter P assumes a continuous path versus the diffusion constant D ϕ . We then start from a point B 2 = (1.45, 0.02) and go in the reverse direction gradually decreasing D ϕ with the same step size as before. For each new parameter, we take as an initial condition the final state of the system from a previous parameter and impose small spatial perturbations of the same form as in (42) in order to allow spatial perturbations to grow provided that SHOM is unstable. We note that for spatially homogeneous systems, spatial variations are of order of a round-off error, and without such initial spatial perturbations, spatial deviations never grow even for parameter values where SHOM is indeed unstable. During this reverse round of continuation simulations, we integrate the system until dδ r (t)/dt < 10 −6 . As a result, by varying the diffusion constant D ϕ , we observe a supercritical transition between SNM and SHOM on the path OB 2 at D ϕ ≈ 0.0125. Additionally, we provide the results of the continuation procedure in terms of the polar order parameter R (8) (cf. Fig. 7(d)) to make the comparison with SHOM transitions.
As the next step, we study phase transitions between SNM and SHOM versus the phase lag α. There are two ways, they can occur. To begin with, let us follow the right path OA 2 in Fig. 7(a). As an initial parameter point, we again set O = (1.45, 0.0075), keep D ϕ constant, and vary α ∈ [1.45, 1.54] with a parameter step size ∆α = 0.005. We follow the same continuation protocol as before and report the results in Fig. 7(e,f). We see that the transition between SNM and SHOM versus α on this path follows the similar scenario as the previously described transition along OB 2 . That is, it is of second order with the bifurcation point α ≈ 1.515. This comes as no surprise as both bifurcations occur close to the order-disorder transition line D ϕ = σ 2 cos α, where the effect of diffusion is substantial. Here, the increase of D ϕ is qualitatively similar to the decrease of α. By doing so, a spatially localized region gradually smooths around.
The second way, the transition between SNM and SHOM occurs, is along the left path OA 1 in Fig. 7(a). By performing the same continuation procedure for the parameter values α ∈ [1.3, 1.45], we observe a hysteresis loop (cf. Figs. 7(e,f)), characteristic to first order transitions. Along the SNM→SHOM path, we come across a bifurcation point α ≈ 1.33 of a saddle-node type. Along the SHOM→SNM path, we find a bifurcation point α ≈ 1.44 of a subcritical type. Apart from the results of the linear stability analysis, which showed us where SNM is observable starting from any initial conditions (except for unstable solutions), we discover the existence of a bistability region where both SHOM and SNM are stable solutions, i.e., α ∈ (1.33, 1.44) with D ϕ = 0.0075. Moreover, we observe some discrepancy between the results of the linear stability analysis and the continuation method. According to the stability analysis, starting from α ≈ 1.38, SHOM should become unstable against spatially dependent perturbations while the continuation method provides α ≈ 1. 44. This is because the stability analysis was performed under the assumption of small microscopic particle velocities v 0 in a region of small diffusion, which is not the case here. Therefore, the numerical analysis of the PDE provides us with a better understanding of how solutions behave far from the order-disorder transition line.
Along bifurcation paths with respect to both parameters, SNM undergoes qualitatively similar transformations. Starting from second order transition points (D ϕ ≈ 0.0125 in Fig. 7(c) and α ≈ 1.515 in Fig. 7(e)), an ellipsoidal shape forms inside a high density layer (cf. Fig. 8(c)) but the layer itself does not disappear completely. In Fig. 8(b), one can observe coexistence of a localized cluster with such a layer for a parameter point even in the middle of a bifurcation path. By decreasing parameters to minimal values with SNM being stable, the localized cluster is most clearly pronounced and the secondary layer dissolves (cf. Fig. 8(a)). The exemplary videos demonstrating temporal evolution of SNM solutions for different parameter values can be found in [3, 1]. Fig. 8 demonstrates qualitative changes in SNM with respect to the phase lag α but one obtains similar results by decreasing the diffusion level D ϕ towards zero.

Conclusion
In this paper, we have discussed the problem of modeling systems of infinitely large populations of nonlocally interacting active Brownian particles. We have developed finite volume schemes to solve a class of nonlinear Vlasov-Fokker-Planck equations obtained in the continuum limit of such systems. According to the continuum limit methodology, these PDEs govern a temporal evolution of nonnegative probability density functions. Taking that into account, we considered the application of positivity-preserving slope limiters and the SSP-RK time discretization in our schemes in order to preserve the probabilistic nature of solutions. Because the problems of interest describe motion of isolated particle flows, the schemes additionally guarantee the conservation of probability mass. Given the theoretical insights on the system properties, we have considered separately one-dimensional problems for spatially homogeneous systems and general three-dimensional problems for spatially inhomogeneous ones. For each case, we have demonstrated that our finite volume schemes yield the second order with respect to discretization procedures. We have demonstrated that the schemes correctly reproduce known stationary and traveling wave solutions. In addition to the presentation of various continuum limit dynamics, we have performed the analysis on phase transitions between three classes of solutions, i.e., SHDM, SHOM, and SNM. This analysis has revealed the existence of both first and second order transitions with respect to changes in either the diffusion level or the phase lag.