Marta Menci

Roberto Natalini

Thierry Paul

MICROSCOPIC, KINETIC AND HYDRODYNAMIC HYBRID MODELS OF COLLECTIVE MOTIONS WITH CHEMOTAXIS: A NUMERICAL STUDY

A general class of hybrid models has been introduced recently, gathering the advantages of multiscale descriptions. Concerning biological applications, the particular coupled structure fits to collective cell migrations and pattern formation phenomena due to intercellular and chemotactic stimuli. In this context, cells are modelled as discrete entities and their dynamics are given by ODEs, while the chemical signal influencing the motion is considered as a continuous signal which solves a diffusive equation. From the analytical point of view, this class of model has been recently proved to have a mean-field limit in the Wasserstein distance towards a system given by the coupling of a Vlasov-type equation with the chemoattractant equation. Moreover, a pressureless nonlocal Eulertype system has been derived for these models, rigorously equivalent to the Vlasov one for monokinetic initial data. For applications, the monokinetic assumption is quite strong and far from real experimental setting. The aim of this paper is to introduce a numerical approach to the hybrid coupled structure at the different scales, investigating the case of general initial data. Several scenarios will be presented, aiming at exploring the role of the different terms on the overall dynamics. Finally, the pressureless nonlocal Euler-type system is generalized by means of an additional pressure term.

Introduction and rigorous results previously obtained

Collective motions modelling is attracting the interest of different research fields, due to the great variety of living and non-living systems exhibiting collective behaviors (see the seminal paper in [START_REF] Vicsek | Collective motion[END_REF] for an introduction the the field). Different approaches have been proposed, depending on the features models aim at reproducing and on the scale of observation. At the microscopic scale, the majority of models are based on elementary mechanical interactions among agents. Alignment, repulsion and attraction kind of interactions are usually considered to model the characteristic tendence to move towards the same direction, keeping cohesion in the group and avoiding collisions [START_REF] Vicsek | Novel Type of Phase Transition in a System of Self-Driven Particles[END_REF][START_REF] Cucker | On the mathematics of emergence[END_REF][START_REF] Cucker | Emergent behavior in flocks[END_REF][START_REF] D'orsogna | Self-Propelled Particles with Soft-Core Interactions: Patterns, Stability, and Collapse[END_REF][START_REF] Strömbom | Collective motion from local attraction[END_REF]. Concerning alignment models, the seminal model in [START_REF] Cucker | On the mathematics of emergence[END_REF] originally proposed to describe the dynamics in flocks of birds was extended in different directions [START_REF] Carrillo | Mean-field limits: from particle descriptions to macroscopic equations[END_REF][START_REF] Natalini | On The Mean Field limit for Cucker-Smale models[END_REF][START_REF] Piccoli | Control to flocking of the kinetic Cucker-Smale model[END_REF].

In the biological field, collective behaviors occur in living processes involving cells dynamics, see [START_REF] Arboleda-Estudillo | Movement Directionality in Collective Migration of Germ Layer Progenitors[END_REF][START_REF] Belmonte | Self-propelled particle model for cell-sorting phenomena[END_REF][START_REF] Hatzikirou | Cellular automata as microscopic models of cell migration in heterogeneous environments[END_REF][START_REF] Méhes | Collective motion of cells: from experiments to models[END_REF] for seminal review papers of the field. The main feature is that collective cells migration is also driven by a chemical stimuli, and not only by mechanical interactions among agents. On the one hand, the microscopic approach allows to model in fine details mechanical interactions among cells. On the other hand, when modelling the evolution in time of a chemical signal influencing the overall dynamics, the microscopic approach is actually not convenient, and a continuum approach, based on reaction-diffusion equation better fullfill the requirement. The structure of hybrid coupled models here considered gather the advantages of microscopic and macroscopic descriptions. Focusing on alignment models, in [START_REF] Di Costanzo | A hybrid mathematical model for self-organizing cell migration in the zebrafish lateral line[END_REF] a model for the morphogenesis in the zebrafish lateral line primordium was proposed. Based on the experimental data in [START_REF] Haas | Chemokine signaling mediates self-organizing tissue migration in the zebrafish lateral line[END_REF][START_REF] Lecaudey | Dynamic Fgf signaling couples morphogenesis and migration in the zebrafish lateral line primordium[END_REF], the model couple alignment and attraction-repulsion kind of interactions with chemotactic effect. The model is hybrid in the following sense: particles are considered as circular shaped entities, whereas the chemical signal is modelled as a continuum. During time evolution, the positions and velocities of each cells are described by second-order Ode, whereas the concentration of the chemical signal is described as the solution of a parabolic equation with term of source and degradation. The original structure has been extended including stochasticity and/or further cell mechanism for different biological phenomena [START_REF] Di Costanzo | A discrete in continuous mathematical model of cardiac progenitor cells formation and growth as spheroid clusters (Cardiospheres)[END_REF][START_REF] Bretti | Estimation algorithm for a hybrid pde-ode model inspired by immunocompetent cancer on-chip experiment[END_REF]. From the analytical point of view, a simplified version of the model in [START_REF] Di Costanzo | A hybrid mathematical model for self-organizing cell migration in the zebrafish lateral line[END_REF] was proposed in [START_REF] Di Costanzo | A hybrid model of collective motion of discrete particles under alignment and continuum chemotaxis[END_REF] to allow a full analytical investigation of the asymptotic behavior, together with well posedeness results in R 2 . Further analytical results on generalized version of hybrid systems can be found in [START_REF] Menci | Global solutions for a path-dependent hybrid system of differential equations under parabolic signal[END_REF][START_REF] Menci | Existence of solutions for hybrid systems of differential equations under exogenous information with discontinuous source term[END_REF][START_REF] Menci | On a coupled hybrid system of nonlinear differential equations with a nonlocal concentration[END_REF]. In all of the previous papers, the hybrid structure considered is at particle level. We begin our work summarizing the main results of the literature concerning the passage from particle to kinetic and macroscopic scales. Consider on R 2dN ∋ ((x i (t)) i=1,...,N , (v i (t)) i=1,...,N ) := (X(t), V (t)) the following vector field [START_REF] Albi | Binary interaction algorithms for the simulation of flocking and swarming dynamics[END_REF] ẋi (t) = v i , vi (t) = F i (t, X(t), V (t)) where

(2)

F i (t, X, V ) = 1 N N j=1 γ(v i -v j , x i -x j ) + η∇ x φ t (x i ) + F ext (x i ),
for any i = 1, . . . , N .

Here F ext denotes an external force, x i , v i are the position and velocity of the i -th cell and φ stands for a generic chemical signal produced by the cells themselves. In particular, φ satisfies the equation [START_REF] Aregba-Driollet | Discrete kinetic schemes for multidimensional systems of conservation laws[END_REF] ∂ s φ s (x) = D∆φ s -κφ s + f (x, X(s)), s ∈ [0, t],

for some κ, D, η ≥ 0 and function f of the form

(4) f (x, X) = 1 N N j=1 χ(x -x j ), χ ∈ C 1 c .
The choice of function χ is justified by the fact that every cell produces the chemical signal, i.e. cells are seen as regions from which the signal arises. The assumption of a compact support recovers the modelling choice of [START_REF] Bretti | Estimation algorithm for a hybrid pde-ode model inspired by immunocompetent cancer on-chip experiment[END_REF][START_REF] Di Costanzo | A discrete in continuous mathematical model of cardiac progenitor cells formation and growth as spheroid clusters (Cardiospheres)[END_REF], where each cell j = 1, .., N is modeled as a disk, centered in its position x j . The function γ : R d × R d → R d models the interactions among agents and it is supposed to be Lipschitz continuous 1 . Initial data are given by (X(0), V (0), φ(0)) = (X in , V in , φ in ).

Let us describe the main steps and results of [START_REF] Natalini | The Mean-Field limit for hybrid models of collective motions with chemotaxis[END_REF], in order to introduce the corresponding kinetic and hydrodynamic limit derived for the particle hybrid structure.

For any fixed function φ in and any t, N the mapping Φ t N = Φ t is defined as (5) Φ t N : R 2dN -→ R 2dN (X in , V in ) -→ (X(t), V (t)) solution of [START_REF] Albi | Binary interaction algorithms for the simulation of flocking and swarming dynamics[END_REF][START_REF] Albi | Pressureless Euler alignment system with control[END_REF][START_REF] Aregba-Driollet | Discrete kinetic schemes for multidimensional systems of conservation laws[END_REF][START_REF] Arboleda-Estudillo | Movement Directionality in Collective Migration of Germ Layer Progenitors[END_REF].

Note that Φ t N is not a flow. In [START_REF] Natalini | The Mean-Field limit for hybrid models of collective motions with chemotaxis[END_REF] a kinetic model was derived, corresponding to system [START_REF] Albi | Binary interaction algorithms for the simulation of flocking and swarming dynamics[END_REF][START_REF] Albi | Pressureless Euler alignment system with control[END_REF][START_REF] Aregba-Driollet | Discrete kinetic schemes for multidimensional systems of conservation laws[END_REF][START_REF] Arboleda-Estudillo | Movement Directionality in Collective Migration of Germ Layer Progenitors[END_REF]. In particular, the following non local in time Vlasov system is derived:

(6) ∂ t ρ t + v • ∇ x ρ t = ∇ v •(ν(t, x, v)ρ t ), ρ 0 = ρ in where (7) ν(t, x, v) = γ * ρ t (x, v) + η∇ x ψ t (x) + F ext (x)
and ψ s satisfies (8) ∂ s ψ s (x) = D∆ x ψ s -κψ s + g(x, ρ s ), ψ 0 = φ in with (9) g(x, ρ s ) = R 2d χ(x -y)ρ s (y, ξ)dydξ. (Lip(f i ) 2 .

The hydrodynamic limit of Cucker-Smale models has provided up to now a large literature, whose exhaustive quotation is beyond the scope of the present paper. The approach and results in [START_REF] Natalini | The Mean-Field limit for hybrid models of collective motions with chemotaxis[END_REF] differ from previous results of the literature. In particular authors do not make use of empirical measures formalism.

Indeed one easily sees that, in the case where ρ in is monokinetic, i.e. [START_REF] Carrillo | Mean-field limits: from particle descriptions to macroscopic equations[END_REF] ρ in (x, v) = µ in (x)δ(v -u in (x)), the monokinetic form is preserved by the Vlasov equation ( 6) and the solution is furnished by the solution of the Euler type system: (11)

∂ t µ t + ∇ x •(u t µ t ) = 0 ∂ t (µ t u t ) + ∇ x •(µ t (u t ) ⊗2 ) = µ t γ(• -y, u t (•) -u t (y))µ t (y)dy + ηµ t ∇ x ψ t + µ t F ext where (12) ∂ s ψ s = D∆ x ψ s -κψ s + χ * µ s , s ∈ [0, t], Ψ 0 = φ in .
Let the Wasserstein distance of order two between two probability measures μ, ν on R m with finite second moments be defined as

(13) W 2 (μ, ν) 2 = inf γ∈Γ(μ,ν) R m ×R m |x -y| 2 γ(dx, dy)
where Γ(μ, ν) is the set of probability measures on R m × R m whose marginals on the two factors are μ and μ, ν (see [START_REF] Villani | Topics in Optimal Transportation[END_REF][START_REF] Villani | Optimal Transport: old and new[END_REF]). For the sake of completeness we here sketch the two main results of [START_REF] Natalini | The Mean-Field limit for hybrid models of collective motions with chemotaxis[END_REF]:

Theorem 1.1. Let ρ in be a compactly supported probability on R 2dN , let Φ t N be the mapping generated by the particles system (1, 2, 3, 4) as defined by [START_REF] Belmonte | Self-propelled particle model for cell-sorting phenomena[END_REF], and let τ ρ in be the function defined in [START_REF] Natalini | The Mean-Field limit for hybrid models of collective motions with chemotaxis[END_REF]Formula (41), Proposition 5.1].

Then, for any t ≥ 0,

W 2 (Φ t N #(ρ in ) ⊗N ) N ;1 , ρ t 2 ≤ τ ρ in (t)    N -1 2 d = 1 N -1 2 log N d = 2 N -1 d d > 2
where ρ t is the solution of the Vlasov equation [START_REF] Brenier | Sticky particle dynamics with interactions[END_REF][START_REF] Brenier | Sticky Particles and Scalar Conservation Laws[END_REF][START_REF] Bretti | Estimation algorithm for a hybrid pde-ode model inspired by immunocompetent cancer on-chip experiment[END_REF][START_REF] Bouchut | On zero pressure gas dynamics[END_REF] with initial condition ρ in provided by [START_REF] Natalini | The Mean-Field limit for hybrid models of collective motions with chemotaxis[END_REF]Theorem 8.1] 2 . Moreover, let us denote by φ t Z in the chemical density solution of (1, 2, 3, 4) with initial data (Z in = (X in , V in ), φ in ) and by ψ t ρ in the one solution of [START_REF] Brenier | Sticky particle dynamics with interactions[END_REF][START_REF] Brenier | Sticky Particles and Scalar Conservation Laws[END_REF][START_REF] Bretti | Estimation algorithm for a hybrid pde-ode model inspired by immunocompetent cancer on-chip experiment[END_REF][START_REF] Bouchut | On zero pressure gas dynamics[END_REF] with initial data (ρ in , φ in ).

Then

R 2dN ∥∇φ t Z in -∇ψ t ρ in ∥ 2 ∞ (ρ in ) ⊗N (dZ in ) ≤ τ c (t)    N -1 2 d = 1 N -1 2 log N d = 2 N -1 d d > 2
where τ c is defined below by [START_REF] Natalini | The Mean-Field limit for hybrid models of collective motions with chemotaxis[END_REF]Formula (54)].

2 We recall that the pushforward of a measure μ by a measurable function Φ is Φ#μ defined by hd(Φ#μ) := (h • Φ)dμ for every test function h. Moreover, (ρ in ) ⊗N stands for ρ in ⊗ ... ⊗ ρ in N times .

Finally, the functions τ (t), τ c (t) depend only on t, Lip(γ), Lip(χ), Lip(∇χ), and the supports of Φ t N #(ρ in ) ⊗N and ρ t , and satisfies the following estimate for all t ∈ R,

τ ρ in (t) ≤ e e Ct , τ c (t) ≤ e e Cct
for some constants C, C c ,depending on Lip(γ), Lip(χ), Lip(∇χ) and |supp(ρ in )|.

Theorem 1.2. Let µ t , u t , ψ t be a solution to the following system

       ∂ t µ t + ∇ x •(u t µ t ) = 0 ∂ t (µ t u t ) + ∇ x •(µ t (u t ) ⊗2 ) = µ t γ(• -y, u t (•) -u t (y))µ t (y)dy + ηµ t ∇ x ψ t + µ t F ext ∂ s ψ s = D∆ x ψ s -κψ s + χ * µ s , s ∈ [0, t], (µ 0 , u 0 , ψ 0 ) = (µ in , u in , ψ in ) ∈ H s , s > d 2 + 1.
where

µ t , u t ∈ C([0, t]; H s ) ∩ C 1 ([0, T ]; H s-1 ), ψ t ∈ C([0, t]; H s ) ∩ C 1 ([0, T ]; H s-2 ) ∩ L 2 (0, T ; H s+1 ) Then ρ t (x, v) := µ t (x)δ(v -u t (x)) solves the following system        ∂ t ρ t + v • ∇ x ρ t = ∇ v •(ν(t, x, v)ρ t ), ν(t, x, v) = γ(x, v) * ρ t + η∇ x ψ t (x) + F ext (x), ∂ s ψ s (x) = D∆ x ψ s -κψ s + g(x, ρ s ), ψ 0 = ψ in , ρ 0 (x, v) = µ in (x)δ(v -u in (x)).
We remark that the derivation of the models is obtained assuming a monokinetic structure of the initial data at the microscopic level. This infers a very strong constraint on the initial data: they have to be included in an N d submanifold of the 2N d dimensional phase space, instead of having the freedom of filling the entire phase space R 2N d . Moreover, when dealing with applications, the monokinetic structure is clearly far from reality. The following question arises naturally: neglecting monokinetic assumption, is there any relation among the different scales? Are the zero-th and first moments (in v) of the solution to the Vlasov system with general initial data ρ in still well approximated by the solution of the Euler system with initial data the zero-th and first moments of the initial condition ρ in ? One of the main objectives of the present article is to gain some insights on the answer to these questions, presenting a numerical study over different scenarios.

The paper is organized as follows: Section 2 reviews the previous numerical results established for the particle, Vlasov and Euler systems inspired by collective dynamics. In Section 3 we state the setting and motivations of our computations, whose precise results are contained in Section 4. Conclusions and perspectives conclude the paper in Section 5. Details of the numerical schemes adopted for the different simulations are exposed in Appendix A.

Let us conclude this section by recalling the three following dynamics involved in this paper, denoted by (P) for (Particles), (V) for (Vlasov) and (E) for (Euler).

(P )                              ẋi = v i , vi = F i (t, X(t), V (t)), (X(0), V (0)) ∈ R 2dN F i (t, X, V ) = 1 N N j=1 γ(v i -v j , x i -x j ) + η∇ x φ t (x i ) ∂ s φ s (x) = D∆ x φ s -κφ s + f (x, X(s)), s ∈ [0, t], φ 0 = φ in , f (x, X) = N j=1 χ(x -x j ). (V )                    ∂ t ρ t + v • ∇ x ρ t = ∇ v •(ν(t, x, v)ρ t ), ρ 0 = ρ in ∈ P(R 2d ) ν(t, x, v) = γ * ρ t (x, v) + η∇ x ψ t (x) -αv ∂ s ψ s (x) = D∆ x ψ s -κψ s + g(x, ρ s ), ψ 0 = φ in , g(x, ρ s ) = χ * ρ s (x). (E)            ∂ t µ t + ∇ x •(u t µ t ) = 0 ∂ t (µ t u t ) + ∇ x •(µ t (u t ) ⊗2 ) = µ t γ(• -y, u t (•) -u t (y))µ t (y)dy + ηµ t ∇ x ψ t -αµ t u t ∂ s ψ s = D∆ x ψ s -κψ s + χ * µ s , s ∈ [0, t].

Previous numerical results

In this Section we present a short overview of previous seminal works of the literature concerning numerical simulations of the models in (P ), (V ) and (E), focusing in particular on the ones that will be recalled in the following sections. As a preliminary observation, note that numerical simulations of Vlasov-type equations arising in collective dynamics modelling can be found in the literature in absence of chemotaxis. To the best of our knowledge this is the first work showing numerical simulations of the coupled system in (V ). The same holds for (E). A complete lists of reference dealing with numerical approaches to Vlasov and Euler equations goes beyond the scope of the paper. In this Section we refer in particular to some seminal and/or recent papers of the literature, highlighting the difference with the structure of the models considered in our paper and with the numerical scheme adopted in the different simulations. Moreover, the numerical tests exhibited in some of the cited works will be used to test the numerical scheme adopted in this paper over well established results.

Particle dynamics Concerning the particle level, numerical simulations of (P ) can be found in [START_REF] Di Costanzo | A hybrid model of collective motion of discrete particles under alignment and continuum chemotaxis[END_REF]. In fact, the discrete-continuum hybrid system investigated corresponds to (P ) with d = 1, choosing a Cucker-Smale like interaction function [START_REF] Carrillo | Critical thresholds in 1D Euler equations with nonlocal forces[END_REF] γ

(v i -v j , x i -x j ) = 1 1 + ∥x i -x j ∥ 2 R 2 β (v j -v i )
where β > 0 is the alignment parameter and R > 0 is the radius of the disk modelling each particle. Initial data are given by initial position and velocity for each particle:

X(0) = X 0 , V (0) = V 0 , (15)
and by the initial concentration of signal, φ 0 = 0.

On a linearised version of the system, authors analytically prove that the particles' aggregate exponentially converges to a state in which all agents share the same position, and the velocity converges to zero. The behavior of the full nonlinear system is investigated with a numerical approach, showing a good agreement between numerical simulations and the theoretical results on the linearised version. For the generalized structure of (P ), introducing different cell mechanisms and stochasticity, numerical simulations oriented to applications on biological phenomena can be found in [START_REF] Di Costanzo | A hybrid model of collective motion of discrete particles under alignment and continuum chemotaxis[END_REF][START_REF] Di Costanzo | A discrete in continuous mathematical model of cardiac progenitor cells formation and growth as spheroid clusters (Cardiospheres)[END_REF][START_REF] Bretti | Estimation algorithm for a hybrid pde-ode model inspired by immunocompetent cancer on-chip experiment[END_REF].

Vlasov equation Several numerical methods for Vlasov equations have been proposed over the years (see [START_REF] Filbet | Comparison of eulerian vlasov solvers[END_REF][START_REF] Elkina | A new conservative unsplit method for the solution of the Vlasov equation[END_REF][START_REF] Albi | Binary interaction algorithms for the simulation of flocking and swarming dynamics[END_REF][START_REF] Carrillo | Particle, kinetic, and hydrodynamic models of swarming Mathematical modeling of collective behavior in socio-economic and life sciences[END_REF][START_REF] Carrillo | Nonoscillatory interpolation methods applied to Vlasov-based models[END_REF][START_REF] Marco | Numerical methods for kinetic equations[END_REF] and the references therein). Focusing on Vlasov equations arising in collective dynamics phenomena, in [START_REF] Albi | Binary interaction algorithms for the simulation of flocking and swarming dynamics[END_REF] several stochastic Monte Carlo methods, improving classic Monte Carlo approaches, are presented, whereas in [START_REF] Carrillo | Nonoscillatory interpolation methods applied to Vlasov-based models[END_REF] deterministic splitting techniques, together with semi-Lagrangian and flux balance methods, are analyzed. To the best of our knowledge, the structure of system (V ), where a Vlasov-like equation is coupled with a parabolic diffusion equation, has not been investigated from a numerical point of view.

Euler system Two are the main features of (E): the presence of the non-local integral term in the Euler equations, and its coupling with the parabolic chemotactic equation. Neglecting chemotaxis, a numerical approach to pressureless hydrodynamic alignment models has been presented in [START_REF] Carrillo | A review on attractive-repulsive hydrodynamics for consensus in collective behavior[END_REF]. Here authors adopt a Lagrangian numerical scheme to approximate the solution of one dimensional hydrodynamic system. In particular, numerical simulations with different initial data are performed, in order to investigate regularity of the solutions, comparing global solutions and solutions exhibiting finite time blow up. Adding a control term, a finite volume scheme to simulate the proposed non-local hydrodynamic system with and without control can be found in [START_REF] Albi | Pressureless Euler alignment system with control[END_REF]. Concerning the coupling aspect, the seminal paper in [START_REF] Gamba | Percolation, morphogenesis, and burgers dynamics in blood vessels formation[END_REF] and reference therein can be compared to the structure of (E), but replacing the non-local integral term with the gradient of a phenomenological pressure. Far from the numerical approach, the recent paper [START_REF] Tadmor | Swarming: hydrodynamic alignment with pressure[END_REF] presents a systematic study of the long-time swarming behavior of hydrodynamic Euler system accounting both the alignment interaction and pressure term. A detailed description of the mixed pressure-alignment structure will be the core of Section 4.4.

Setting and motivation of the present numerical study

In the present paper, we study numerically the solutions of the Vlasov and Euler systems starting with different initial data, both with and without chemotaxis.

Precisely, we study the system at particle level [START_REF] Di Costanzo | A hybrid model of collective motion of discrete particles under alignment and continuum chemotaxis[END_REF] for different values of the parameters β and η. As we have seen in Section 2, parameter β leads to conditional or unconditional flocking for the Cucker-Smale model, which can be obtained from [START_REF] Di Costanzo | A hybrid model of collective motion of discrete particles under alignment and continuum chemotaxis[END_REF] considering η = 0. Values of η ̸ = 0 introduces the coupling of the microscopic system to the chemical activity. Here we are interested in positive chemotactic phenomena (η > 0), since inspired by biological phenomena, cells move towards higher concentrations of chemicals. The path microscopic -→ V lasov -→ Euler corresponds to successive lacks of knowledge of the system: the meanfield limit, leading to the Vlasov equation, is obtained out of the microscopic system by averaging on all particles but one, and the passage from Vlasov to Euler retains only the moments of order zero and one in the momentum of the solution.

On the other side, the benefit of these succesive transformation is clear: the kinetic Vlasov equation deals with probabilities on finite dimensional phase spaces instead of the infinite dimensional ones of the original system, and the Euler setting describes the system by mean of densities and velocity fields on physical (configuration) spaces, that is quantities directly observable.

Finally, as explained in Section 1, the passage form Vlasov to Euler was rigorously obtained in [START_REF] Natalini | The Mean-Field limit for hybrid models of collective motions with chemotaxis[END_REF] by a monokinetic hypothesis on the solution to Vlasov of the form ρ t (x, ξ) = µ t (x)δ(ξ -u t (x)), shown to be rigorously equivalent for the couple (µ t , u t ) satisfying the Euler system.

Inspired by the results obtained for ρ t with monokinetic form, one might wonder if, for more general (non monokinetic) solution to Vlasov, the quantities defined by [START_REF] Cucker | Emergent behavior in flocks[END_REF] still satisfy, at least approximatively, the Euler system of equations.

From a numerical point of view, this would be even more appealing for computations based on finite difference schemes, since solving a system of two PDEs, each of them in a m-dimensional space, is much more economical than solving one, having dimension 2m.

The motivation of our work is therefore fourfold:

(i) showing that the different features of the microscopic dynamics are still present in the numerics of Vlasov and Euler situations. (ii) comparing, for ρ in non-monokinetic, the moments of the solution to Vlasov as defined by ( 16) with the solution of the Euler system with initial data ( 16)

   µ in (x) = ρ in (x, ξ)dξ, µ in (x)u in (x) = ξρ in (x, ξ)dξ.
(iii) evaluating on the settings (i) and (ii) the influence of the chemical coupling by varying the value of the parameter η from 0 to positive values. (iv) investigating the role of an additional pressure term in the non local Euler system and tuning its effect in order to recover agreement with the kinetic scale.

Main results

In the following section we present our numerical study across the different scales, simulating different scenarios for (P ), (V ) and (E). The aim of this paper is to analyze numerically some of the open problems related to generalizations of theoretical results obtained in [START_REF] Natalini | The Mean-Field limit for hybrid models of collective motions with chemotaxis[END_REF]. We present some numerical insights concerning general initial data, far from the monokinetic ones, and we investigate the role of the different terms, especially at kinetic and hydrodynamic level. We restrict ourselves to the 1-dimensional case (d = 1), leaving higher dimensional cases for future researches. In particular, for partial version of the models investigated (i.e. without chemotaxis and/or damping), we compare our results with the ones of the literature, when available, running the simulations starting with the same initial data. Details of the numerical implementation can be found in the Appendix A. We run the numerical codes on a laptop equipped with an Intel Core i7-1060NG7 processor and 16 GB RAM.

4.1. Numerical Simulations of (P) and (V). We start our analysis comparing particle and kinetic levels. To the best of our knowledge, this is the first time that a Vlasov-like equation is coupled to a chemotaxis parabolic equation. Numerical simulations of (V ), neglecting chemotaxis and damping, can be found in [START_REF] Albi | Binary interaction algorithms for the simulation of flocking and swarming dynamics[END_REF][START_REF] Carrillo | Particle, kinetic, and hydrodynamic models of swarming Mathematical modeling of collective behavior in socio-economic and life sciences[END_REF]. Starting with the same initial conditions of these papers, we recover the expected behavior using a different numerical scheme. Then we add the coupling with the chemotaxis and we show the effect of the gradient term. In conclusion, we simulate the scenario of a pure chemotactic interaction. 4.1.1. alignment without chemotaxis. Neglecting chemotaxis and considering only an alignment kind of interaction, systems (P ) and (V ) reduce to a Cucker-Smale model in a particle and kinetic regime, respectively:

     ẋi = v i , vi = 1 N N j=1 1 (1 + |x i -x j | 2 ) β (v j -v i ) i = 1, ..., N (17)
       ∂ t ρ t + v∂ x ρ t = ∂ v (ν(t, x, v)ρ t ), ν(t, x, v) = R 2 v -w (1 + |x -y| 2 ) β ρ t (y, w)dydw. (18)
In the literature, numerical simulations of ( 18) are often based on Montecarlo techniques. Here, dealing with a 1-dimensional case, we implemented a first-order scheme (see Appendix A for details). The drawback of the constrain given by a CFL condition is overcome by the easier implementation of the scheme. As preliminary tests, we run the the numerical code starting with the same initial distribution of [START_REF] Albi | Binary interaction algorithms for the simulation of flocking and swarming dynamics[END_REF], i.e.

(

) ρ 0 (x, v) = 1 2πσ x σ v e -x 2 2σ 2 x e -(v+v 0 ) 2 2σ 2 v + e -(v-v 0 ) 2 2σ 2 v , with v 0 = 3.5, σ x = √ 0.1, σ v = √ 0. 19
5. We consider two different scenarios, depending on the value of β, which is the key parameter of the alignment kind of interaction, as proved in [START_REF] Ha | A simple proof of Cucker-Smale flocking dynamics and mean-field limit[END_REF]. In particular, the convergence to a uniform velocity, regardless of the initial state of the system, is ensured only for value of β ≤ 0.5. In the following Test 1 -Test 4 the phase space representation is obtained considering the space-velocity domain [-20, 20] × [-5, 5], with ∆x = ∆v = 0.01. The simulation runs over the time interval [0, T ] = [0, 5] with ∆t = 0.001.

Figure 1 shows three different screenshots of the numerical simulation of (17) (first line) and (18) (second line) performed with β = 0.05. Since β ≤ 0.5, at particle level we observe that a flocking state is reached, and particles' velocities tend to a common (null) one. In a similar way, we observe that the distribution function tends to be distributed only along the spatial dimension, concentrating around the null value in the velocity space. Figure 2 shows three screenshots of the numerical simulation performed with β = 0.95. In this case, since β > 0.5, the convergence to a uniform velocity is no longer reproduced. At particle level, the two groups start migrating with positive or negative velocities, without reaching an asymptotic flocking state. The microscopic behavior is coherently reproduced at kinetic level.

               ẋi = v i , vi = 1 N N j=1 1 (1 + |x i -x j | 2 ) β (v j -v i ) + η∂ x φ t (x i ), i = 1, ..., N ∂ s φ s = D∂ 2 x φ s -κφ s + N j=1 χ R (x -x j ), s ∈ [0, t], (20)
where R, η, κ, D > 0. In particular η > 0 tunes the effect of the chemotactic gradient, and R denotes the radius of the characteristic function centered in x j , modelling the fact that the chemical signal is produced by the particles themselves. Let now focus on the novel coupled Vlasov-chemotaxis system:

                   ∂ t ρ t + v∂ x ρ t = ∂ v (ν(t, x, v)ρ t ), ν(t, x, v) = R 2 v -w (1 + |x -y| 2 ) β ρ t (y, w)dydw + η∂ x ψ t (x), ∂ s ψ s (x) = D∂ 2 x ψ s -κψ s + R x+R x-R ρ s (y, w)dydw, s ∈ [0, t]. (21)
To the best of our knowledge, this is the first paper presenting numerical simulations of this coupled system. The obtained results, and the comparison with the corresponding particle simulations, give some interesting insights.

In Test 3 we consider the same initial scenario and parameters of Test 2, adding chemotaxis. In particular we chose R = 10∆x, D = 1, κ = 0.01, η = 1.4, and initial null concentration ψ 0 = 0. Figure 3 shows three screenshots of the performed numerical simulation. At particle level, in [START_REF] Di Costanzo | A hybrid model of collective motion of discrete particles under alignment and continuum chemotaxis[END_REF] a rigorous proof of the asymptotic behavior of the solution to [START_REF] Di Costanzo | A hybrid model of collective motion of discrete particles under alignment and continuum chemotaxis[END_REF] is given. In particular, it has been proved that the introduction of the chemotactic effect ensures convergence to a flocking state, even in the cases of conditional flocking (β > 0.5). Numerical evidence show that the behavior analytically proved on the microscopic scale, is recovered also in the kinetic regime. 4.1.3. chemotaxis without alignment. To conclude our analysis, we present a numerical simulation in the case of a pure chemotactic kind of interaction, hence neglecting the alignment term. The system reads:

               ∂ t ρ t + v∂ x ρ t = ∂ v (ν(t, x, v)ρ t ), ν(t, x, v) = η∂ x ψ t (x), ∂ s ψ s (x) = D∂ 2 x ψ s -κψ s + R x+R x-R ρ s (y, w)dydw, s ∈ [0, t]. (22)
Figure 4 shows three screenshots of the performed numerical simulation. Both at particle and kinetic level, the absence of the alignment kind of interaction results in an oscillating behavior around the center of mass. As in Test 3, due to the chemotactic gradient, particles tend to concentrate, but a consensus state is not reached.

We remark that the initial distribution of Tests 1-4 is clearly far from the monokinetic case. The comparison shown in this section shows a good level of correspondence between particle and kinetic scale in the different scenarios, and it represents the preliminary step of our analysis. In the next Sections we focus on numerical simulations of the hydrodynamic system (E), recovering previous results of the literature and exploring the novel coupling with the chemotaxis equation.

Numerical Simulations of (E).

As seen in Section 2, pressureless Euler systems have been investigated in the literature, both from a theoretical and numerical perspective. Also at this level, the novel aspect lies in the coupling with the chemotaxis equation. The model here considered writes as a hyperbolic-parabolic non local system for macroscopic quantities: density, momentum and concentration of a chemoattractant:

                       ∂ t µ t + ∂ x (u t µ t ) = 0 ∂ t (µ t u t ) + ∂ x (µ t (u t ) 2 ) = µ t γ(• -y, u t (•) -u t (y))µ t (y)dy +ηµ t ∂ x ψ t ∂ s ψ s = D∂ 2 x ψ s -κψ s + x+R x-R µ s (y)dy, s ∈ [0, t]. (23)
We investigate the role of chemotaxis, comparing the case with η = 0 and η ̸ = 0. The numerical scheme implemented is based on relaxation techniques originally proposed in [START_REF] Aregba-Driollet | Discrete kinetic schemes for multidimensional systems of conservation laws[END_REF] (see Appendix A for the details). As a preliminary test, we simulate the system in absence of chemotaxis, starting with the initial condition chosen as in [START_REF] Carrillo | A review on attractive-repulsive hydrodynamics for consensus in collective behavior[END_REF]. Then we explore the role of chemotaxis and damping in the same scenario. 4.2.1. without chemotaxis: finite time blow up of solution. Neglecting chemotaxis, system [START_REF] Filbet | Comparison of eulerian vlasov solvers[END_REF] reduces to a 1-dimensional pressureless Euler-alignment system:

   ∂ t µ t + ∂ x (u t µ t ) = 0 ∂ t (µ t u t ) + ∂ x (µ t (u t ) 2 ) = µ t γ(• -y, u t (y) -u t (•))µ t (y)dy. ( 24
)
Analytical results and numerical insights on the stability of flock solutions for this system can be found in [START_REF] Carrillo | Critical thresholds in 1D Euler equations with nonlocal forces[END_REF][START_REF] Carrillo | A review on attractive-repulsive hydrodynamics for consensus in collective behavior[END_REF]. In [START_REF] Carrillo | Critical thresholds in 1D Euler equations with nonlocal forces[END_REF], authors distinguish between supercritical and subcritical region: initial data lying in the supercritical region lead to blow up of the solution, whereas starting with data in the subrcritical region global existence of the solution is ensured. In [START_REF] Carrillo | A review on attractive-repulsive hydrodynamics for consensus in collective behavior[END_REF] different numerical simulations in the spatial domain [-0.75, 0.75] with data lying in both regions are performed. In Test 5 we approximate the solution of ( 24) starting with the same initial data defined in [START_REF] Carrillo | A review on attractive-repulsive hydrodynamics for consensus in collective behavior[END_REF]:

(25) µ 0 (x) = c 1 cos πx 1.5 , (26)
u 0 (x) = -c 2 sin πx 1.5
where c 1 is the normalization factor, and c 2 > 0 is varied in the different simulations.

The results obtained are in good agreement with the ones presented in [START_REF] Carrillo | A review on attractive-repulsive hydrodynamics for consensus in collective behavior[END_REF], where a Lagrangian numerical scheme is adopted. Figure 5 (first line) shows the density and velocity profile obtained for c = 0.2 for three different time instants. The initial data lies in the subcritical region, and the convergence of the velocity to zero corresponds to the expected global consensus. On the contrary, Figure 5 (second line) shows the density and velocity profile obtained for c = 0.5, hence initial data in the supercritical region.

Plotting the velocity at different times, we observe that the derivative is getting larger negative at the origin as time goes. This implies high and rapidly increasing values of the density solution, which almost cannot be solved numerically after t ≈ 3. The obtained results represent our starting point for numerical investigations on the effect of the introduction of a chemotactic effect and a damping term. 4.2.2. with chemotaxis and/or damping. Let now consider the coupled system in [START_REF] Filbet | Comparison of eulerian vlasov solvers[END_REF] adding a damping term:

                       ∂ t µ t + ∂ x (u t µ t ) = 0 ∂ t (µ t u t ) + ∂ x (µ t (u t ) 2 ) = µ t γ(• -y, u t (•) -u t (y))µ t (y)dy +ηµ t ∂ x ψ t -αµ t u t ∂ s ψ s = D∂ 2 x ψ s -κψ s + x+R x-R µ s (y)dy, s ∈ [0, t]. (27)
where η and α are nonnegative constants regulating the chemotactic and the damping influence, respectively. From a modelling perspective, the introduction of damping is justified and necessary for biological applications, in particular when dealing with cells moving on a substrate. Figure 6 shows the results of three different scenarios. First, we simulate the setting with data in the subcritical region already seen (c 2 = 0.2, α = 0) adding the chemical influence (η = 1). Numerical evidence shows that the global consensus is still preserved, and the presence of a chemotactic influence results in higher values of the density. In the other scenarios (Figure 6, second and third lines) we investigate the role of chemotaxis and damping in suppressing the blow up phenomena (c 2 = 0.5). Numerical evidence (Figure 6, second line) shows that a pure chemotactic effect (η = 1, α = 0) enhances blow up phenomena: density value at the origin is higher than the case without chemotaxis, and the derivative of the velocity is sharper. The presence of a damping term seems to play the crucial role, as shown in (Figure 6, third line). Here, the blow up exhibited in the same setting without damping (see Figure 5) is suppressed: velocity converges to zero and global consensus is reached.

Numerical comparison between (V) and (E).

In the following sections we present the main results of our work. One of the main aims is to get numerical insights on the behavior of the solutions of (V ) and (E), in particular starting with non monokinetic initial data, being not covered by theoretical results.

First, we introduce the following notation for the zero-th and first order moments of the solution of (V ), respectively (28)

ν t 0 (x) = ρ t (x, ξ)dξ, (29)
ν t 1 (x) = ξρ t (x, ξ)dξ,
for any t ≥ 0.

In the following tests we compare the moments of the solution to Vlasov with the solution of the Euler system with initial data (30)

   µ 0 (x) = ν 0 0 (x), Q 0 (x) = ν 0 1 (x)
, where Q t := µ t u t denotes the momentum, for any t ≥ 0. Choosing different ρ 0 we simulate different scenarios, ranging from almost monokinetic to largely non-monokinetic initial data. The aim is to compare the behavior of ν t 0 and ν t 1 with density µ t and momentum Q t solving (E) starting with [START_REF] Méhes | Collective motion of cells: from experiments to models[END_REF].

Monokinetic initial data:

We recall that analytical results obtained in [START_REF] Natalini | The Mean-Field limit for hybrid models of collective motions with chemotaxis[END_REF] ensure a complete correspondence between Vlasov moments and Euler solutions, in case of monokinetic initial data. As a preliminary test, we consider a numerical approximation of the monokinetic initial data, choosing a gaussian distribution with low value of σ v parameter. In greater detail, in Test 7 we assume (31)

ρ 0 (x, v) = 1 2πσ x σ v e -(x-x 0 ) 2 2σ 2 x + -(v-v 0 ) 2 2σ 2 v with x 0 = -2, v 0 = 1.5, σ x = √ 0.2, σ v = √ 0.001.
Clearly, even choosing a small value for σ v , we are considering an approximation of the monokinetic case, and we do not expect a full correspondence with the analytical results on this case.

Figure 7 shows the plot of ρ 0 in the phase space, and the profile of initial data for (E) (blue dotted line), which are defined in [START_REF] Méhes | Collective motion of cells: from experiments to models[END_REF] and correspond to Vlasov moments (cyan line). Figure 8 shows the behavior, at the same time instant t = 2, of ν t 0 and ν t 1 together with µ t and Q t respectively, both without (first line) and with chemotaxis (second and third line). The evolution of Vlasov dynamics preserves the monokinetic (approximated) initial structure, and Vlasov moments and Euler solutions are in good agreement. Numerical evidence shows that adding chemotaxis, the approximated monokinetic-like structure is not preserved for long time. As already observed at particle level, the presence of a chemical gradient influencing the dynamics tends to concentrate the density profile in the same point, and to asymptotically let the velocities converge to zero. This can be seen in the phase space (first column, second line). Moreover, introducing also a damping term, the momentum tends to zero faster (third line). Focusing on the Vlasov moments, we observe that as soon as the monokinetic structure is lost, the agreement between Vlasov moments and Euler solution is worse. In particular, we observe a blow up of the density profile. In conclusion, the monokinetic property seems to play a crucial role. For that reason we exploit with several tests the case of largely non monokinetic initial data. Vlasov phase space ν t 0 and µ t ν t 1 and Q t Non monokinetic initial data: As non monokinetic initial data, we consider the distribution function already introduced in (19). Figure 9 shows initial conditions for the following Test 8 and Test 9, as defined in [START_REF] Méhes | Collective motion of cells: from experiments to models[END_REF].

In Test 8 we compare the moments ν t 0 and ν t 1 of the solution of ( 18)-( 19) with the solution of system [START_REF] Gamba | Percolation, morphogenesis, and burgers dynamics in blood vessels formation[END_REF] with initial conditions [START_REF] Méhes | Collective motion of cells: from experiments to models[END_REF]. Since u 0 (x) = 0, the velocity u t remains zero, being a stationary solution for the Euler system. On the contrary, the first order moment profile shows a separation due to the effect of speeds. Being far from the monokinetic case, in Vlasov dynamics v is a variable which is not a priori related to x, whereas in Euler the velocity field depends on the position. As shown in Figure 10, there is no agreement between Vlasov moments and Euler solutions starting with that initial condition.

Let now solve [START_REF] Hatzikirou | Cellular automata as microscopic models of cell migration in heterogeneous environments[END_REF], hence considering also a chemotactic effect. Figure 11 shows three different snapshots of the numerical simulation performed. At the beginning, we observe that the spreading of the speed in the density, solution of Vlasov equation, is stronger than the chemistry, and two separate peaks appear. With respect to the previous case, the presence of the chemotactic gradient prevent the separation observed for longer time, leading to the convergence to a unique bump. Concerning the first order moment (Vlasov) and the momentum (Euler), at the beginning they are both null. Due to the presence of chemicals, u t (x) = 0 is no longer a stationary solution for Euler. In particular in presence of a damping term (α = 0.8) first order moment and momentum asymptotically converge to zero, as shown in Figure 12. We observe that, as time grows, the agreement between Vlasov moments and Euler solution is almost recovered especially for the density. Numerical evidence hence shows that, far from the monokinetic case, the dissipation of energy induced by damping term helps in recovering a good agreement. In fact, the role of the second-order moment, that we are here neglecting, is reduced and almost negligible.

For the sake of completeness we consider the scenario without alignment, where only a chemotactic effect rules the dynamics. At kinetic level, we have already discussed the results (see Test 4 and Figure 4). Figure 13 shows the results of the comparison. We observe that the agreement between Vlasov moments and solutions of Euler is not recovered. The zero-th order moment shows the observed oscillations, with two bumps separating and merging, as reflecting by the sign of the first-order moment. The oscillating behavior is not reproduced in the solution of Euler. Regardless of the comparison between the two scale, numerical simulations in largely non monokinetic cases shed light on an important result. As already stated, the main advantage of (E) is to clearly preserve the particle interactions in the passage from microscopic to the macroscopic scale, which are included in the nonlocal integral term on the righthand side. Neglecting the nonlocal integral term and damping, the system reduces to a pressureless Euler system, that arises in the modeling of sticky particles and has been investigated at theoretical level in several works of the literature (see [START_REF] Brenier | Sticky Particles and Scalar Conservation Laws[END_REF][START_REF] Brenier | Sticky particle dynamics with interactions[END_REF][START_REF] Bouchut | On zero pressure gas dynamics[END_REF] for seminal papers). One of the main features of pressureless gas system is the development of delta-shocks. Numerical simulations in Test 9 -Test 10 actually exhibit a blow-up in the density, which can be seen as a drawback of the pressureless structure of (E). In particular we conclude that in a general setting, the presence of the integral term and/or the chemotactic gradient is not enough dissipative to prevent a delta-shocks phenomena. In the next section we investigate the role of damping term, proposing a mixed pressure-nonlocal Euler system incorporating the system (E). Comparison between zero-th order moment of the solution of (V) and µ t , first order moment of the solution of (V) and

Q t at a)-b) t = 1, c)-d) t = 2, e)-f) t = 3.
Here η = 0 (no chemotaxis), α = 0. 4.4. Numerical comparison between (V) and (E) with pressure. We observed that, in the general case, the presence of the damping term helps in recovering agreement between Vlasov moments and the macroscopic quantities, density and momentum, solving the nonlocal pressureless Euler system. In absence of damping, or if the γ function is not enough dissipative, a discrepancy between the two scales can be expected. The structure of (E) is actually due to the monokinetic property of the solution of (V ), but in the general case, also the second order moment should be involved when taking the momenta of the Vlasov solution. And the second moment would in turn involve the third one and so on. Therefore the hierarchy of equations satisfied by the moments would not be closed. In the following we extend the Euler system (E) with an additional term approximating to a certain extent the role of the second order moment. In the literature of hyperbolic chemotaxis model, the model in [START_REF] Gamba | Percolation, morphogenesis, and burgers dynamics in blood vessels formation[END_REF], originally conceived to reproduce the early stages of the vasculogenesis process, shares a common structure with (E). The main difference is that the nonlocal integral term in (E) is replaced by a phenomenological density dependent pressure function, which does not have a clear connection with the microscopic interactions among cells, but it plays a crucial role in stabilizing the system. Here we mix the original structure of (E) adding a pressure term inspired by the model in [START_REF] Gamba | Percolation, morphogenesis, and burgers dynamics in blood vessels formation[END_REF]. Note that the nonlocal-pressure structure has been investigated in the context of hydrodynamic limit of agent-based dynamics (without chemotaxis) in a recent paper [START_REF] Tadmor | Swarming: hydrodynamic alignment with pressure[END_REF], proposing rigorous results on the long-time behavior of alignment Euler systems for a general class of entropy pressure.

Our aim is to investigate, from a numerical perspective, the effect of the additional pressure term, and to gain some insights on possible agreement between the scales, far from the monokinetic assumption on initial data. The hydrodynamic alignment Euler system we here consider reads:

                       ∂ t µ t + ∂ x (u t µ t ) = 0 ∂ t (µ t u t ) + ∂ x µ t (u t ) 2 + εP (µ t ) = µ t γ(• -y, u t (•) -u t (y))µ t (y)dy +ηµ t ∂ x ψ t -αµ t u t ∂ s ψ s = D∂ 2 x ψ s -κψ s + x+R x-R µ s (y)dy, s ∈ [0, t]. (32)
We assume a pressure law for isentropic gases P (µ) = µ p , choosing p = 2 in our simulations. Parameter ε ≥ 0 tunes the effect of the pressure term. Choosing ε = 0, we recover the pressureless system [START_REF] Hatzikirou | Cellular automata as microscopic models of cell migration in heterogeneous environments[END_REF]. In the following we present the results of numerical simulations performed for different values of ε. As in the previous case, the effect of the damping term is investigated.

As a preliminary test, we run the system starting with the (approximated) monokinetic initial data in [START_REF] Menci | Global solutions for a path-dependent hybrid system of differential equations under parabolic signal[END_REF] following the same approach already considered for Test 7.

Figure 14 shows the comparison between the first two moments of the solution of ( 21) and density-momentum solutions of (32) at three different time-steps, for different values of ε ≥ 0.

Varying the value of ε we get different solutions of system [START_REF] Menci | Existence of solutions for hybrid systems of differential equations under exogenous information with discontinuous source term[END_REF]. As expected, the presence of a pressure term introduces a dissipative behavior, which increases with the value of ε. Tuning the value of the parameter starting with ε = 0, we observe that, at least from a qualitative perspective, there exists a critical value ε * > 0 such that ε < ε * the pressure term is not enough to counterbalance the density blow up, whereas if ε > ε * the effect of the pressure term is predominant and the dynamic is completely different from the kinetic one. The results exhibited in Figure 14 show that ϵ * = 0.01 represents a good candidate as optimal value minimizing the distance between ν t 0 and µ t .

In order to quantitatively estimate the agreement between ν t 0 and µ t , and ν t 1 and Q t , we compute, at each iteration k = 1, 2, ... the value of distances

E k 0 (ε), E k 1 (ε) defined by (33) E k 0 (ε) := ∥ν k 0 -µ k (ε)∥ L 1 (Ω) E k 1 (ε) := ∥ν k 1 -Q k (ε)∥ L 1 (Ω) where µ k (ε), Q k (ε)
denote the solution of (32) obtained for a fixed value of ε at iteration k.

Table 1 exhibits the obtained values for the three time iterations considered in Figure 14. Computing the quantities at any iteration k, we observe that the values of E k 0 and E k 1 decrease as ε increases from 0 to a certain value ε * , whereas their value increase for ε > ε * , giving insights for the existence of a minimum point. To this end, we run an optimization procedure, aiming at finding, at each iteration, the value ϵ * k ∈ [0, 0.05] minimizing E k 0 . The same can be done for E k 1 (data not shown). Table 2 shows the results obtained for the selected time instants. Considering an accuracy of three decimal places, we observe that the optimal values are slightly higher than our qualitative guess. As representative value we choose ε 1 = ε 2 /5 and ε 3 = ε 2 * 5, where ε 2 represents the numerical (qualitative) candidate for the optimal value of parameter ε. Time instants correspond to the ones selected for snapshots in Figure 14.

E k 0 ε 1 = 0.002 ε 2 = 0.01 ε 3 = 0.05 t k = 0 0 0 0 t k = 1 0.
ε * k E k 0 (ε * k ) t k = 0 0 0 t k = 1 0.013 0.07 t k = 1.5 0.015 0.14 t k = 2 0.017 0.55 Table 2. Computation of ϵ *
k and corresponding value of E k 0 . Time instants correspond to the ones selected for snapshots in Figure 14.

We recall that the correspondence between Vlasov moments and Euler solutions, in the monokinetic case, is analytically proved in [START_REF] Natalini | The Mean-Field limit for hybrid models of collective motions with chemotaxis[END_REF] for the pressureless case, and there is no a priori reason to expect some sort of correspondence adding the densitydependent pressure term. The results shown in Figure 14 suggest that the moments of the solution of Vlasov system (21), starting with something that is like-monokinetic, look like solutions of Euler with pressure, at least for some values of ε.

In a second test we consider a general largely non monokinetic scenario. Initial data are shown in Figure 15 . As defined in [START_REF] Méhes | Collective motion of cells: from experiments to models[END_REF], µ 0 and Q 0 are choosen as zero and first moment of the following initial distribution (34) 16 and Figure 17 show three snapshots of the performed numerical simulations without and with damping term, respectively. In particular we plot ν 0 and ν 1 against the solution of (32) for two values of ε, corresponding to the case of pressureless Euler system (ε = 0) and choosing the value of the parameter computed as the mean value of the optimal values ε * k collected at each iteration. For the sake of completeness, in Table 3 and Table 4 we present the results obtained for the three sample time instants chosen in the corresponding Figure 16 and Figure 17.

ρ 0 (x, v) = 1 √ 2πσ x σ v e -(x 2 -x 1 ) 2σ 2 x e -(v-v 1 ) 2 2σ 2 v + e -(x 2 -x 2 ) 2σ 2 x e -(v-v 2 ) 2 2σ 2 v with x 1 = -2, v 1 = 1.5, x 2 = 2, v 2 = -2.5, σ x = √ 0.2, σ v = √ 0.5. Figure
Starting with non monokinetic initial data, the agreement between Vlasov moments and Euler solutions, even adding pressure, is not so evident as in the approximated monokinetic case. Focusing on the density profile, as shown in Figure 16, first column, choosing an optimal value for ε helps in suppressing the blow-up occurring with ε = 0. At the beginning, the agreement obtained by adding the pressure term seems promising. As time grows, the blow-up is still reduced but the density and zero-moment profiles are no longer overlapped. In fact, due to nonzero initial momentum, the zero-th order moment of Vlasov solution travels, moving leftwards. The comparison for longer time improves introducing also a damping term, considering a nonnull value of α both in [START_REF] Di Costanzo | A discrete in continuous mathematical model of cardiac progenitor cells formation and growth as spheroid clusters (Cardiospheres)[END_REF] and in [START_REF] Menci | Existence of solutions for hybrid systems of differential equations under exogenous information with discontinuous source term[END_REF]. Actually, the pressure term affects only the Euler system, whereas the damping term involves a coherent dissipation both in the Vlasov and in the Euler dynamics. The obtained results are shown in Figure 17. The combined effect of the pressure term, reducing the density blow up, and of the damping term both at kinetic and macroscopic level, results in an improved agreement. Moreover, for longer time, also a good agreement between first moment and momentum is reached, as shown in Figure 17 f).

ε * k E k 0 (ε * k ) t k = 0 0 0 t k = 1

Conclusions

In this paper we investigated numerically three different levels of descriptions of hybrid models for dynamics of interacting cells subject to chemotaxis: the microscopic position-velocity description (P ), the phase space density Vlasov type description (V ) and the physical space density-velocity field Euler type description (E). At particle level, in the last years several analytical and numerical results have been presented [START_REF] Di Costanzo | A hybrid model of collective motion of discrete particles under alignment and continuum chemotaxis[END_REF][START_REF] Menci | Global solutions for a path-dependent hybrid system of differential equations under parabolic signal[END_REF][START_REF] Menci | Existence of solutions for hybrid systems of differential equations under exogenous information with discontinuous source term[END_REF][START_REF] Menci | On a coupled hybrid system of nonlinear differential equations with a nonlocal concentration[END_REF], whereas kinetic and corresponding hydrodynamic limit have been analytically studied in [START_REF] Natalini | On The Mean Field limit for Cucker-Smale models[END_REF][START_REF] Natalini | The Mean-Field limit for hybrid models of collective motions with chemotaxis[END_REF]. In particular, a complete correspondence between Vlasov moments and Euler solutions have been proved in [START_REF] Natalini | The Mean-Field limit for hybrid models of collective motions with chemotaxis[END_REF] under monokinetic assumption, which is a strong requirement in applications. In our work we start investigating the behavior of the solutions far from the monokinetic case, obtaining numerical insights that can be seen as a first step for a rigorous study in general contexts. Although each passage between the three levels consists in loosing some information -averaging on every but one particle and letting the number of particles diverge for (P ) → (V ) and taking only zero-th and first moments in velocity for (V ) → (E) -we show that in many situations the most striking features of the particle level (P ) are preserved by passing from the paradigm of N particles to the one of Vlasov at N → ∞. Concerning the passage (V ) → (E), numerical evidence shows that the fidelity of Euler versus Vlasov is increased by the presence of the chemical interaction and damping term, leading to positive results even for situations which are not covered by rigorous analytical results. With respect to previous Euler systems coupled with chemotactic effect, starting with the seminal work in [START_REF] Gamba | Percolation, morphogenesis, and burgers dynamics in blood vessels formation[END_REF], the Euler system (E) keeps memory of the interactions at the microscopic level through the nonlocal integral term. This is an important feature in biological applications, especially when dealing with experimental data. However, two main drawbacks need to be considered. First, it is well-known that pressureless Euler systems exhibit finite time blow-up of the solution, depending on the initial data ([13], Theorem 3.1). Moreover, from a numerical point of view, the approximation of the integral term is computationally costly and need to be improved, aiming at 2D and 3D simulations. In order to control blow up phenomena, and to reduce the effect of neglected momenta, we investigated the effect of an additional pressure term, still keeping the nonlocal integral term. In particular, we let the effect of the pressure term variable through the parameter ε. Tuning its value we investigate the existence of an optimal value ε * realizing an improved correspondence between Vlasov moments and Euler solutions. From a numerical point of view, the transition (V ) → (E) realizes an economical numerical gain, passing from a 2d PDE to a system of two 1d PDEs. Establishing correlations among the different scales of descriptions allows to improve the computational effort, and it is one of the main perspective for applications to biological systems, aiming at developing the most appropriate macroscopic model as their digital twins.

Appendix A. Numerical schemes

In the following we detail the numerical schemes used for the numerical simulations of the different equations.

Particle level: numerical scheme for (P) For the ODE part of the coupled system, we implement the well-known implicit Euler method, whereas the chemotaxis equation is discretized computing at each time the value of the solution on grid point discretizing the spatial domain. We consider the simulation time interval [0, T ], for some T > 0, divided into subintervals of same size ∆t, assuming that T is a multiple of ∆t to avoid rounding. We define n T := T ∆t ∈ N, and denote the time grid points as {t 0 , . . . , t n T }. Finally, (x k i , v k i ) denote approximations of position and velocity of agent i at time t k , for any i = 1, ..., N , k = 1, ..., n T starting with initial condition (x 0 i , v 0 i ).

Focusing on the chemotaxis equation, we eliminate the stiff degradation term performing the classical exponential transformation φ t (x) = e -κs u t (x), where u solves [START_REF] Natalini | The Mean-Field limit for hybrid models of collective motions with chemotaxis[END_REF] ∂

t u = D∂ 2 x u + e κt N j=1 χ(x -x j ),
We approximate the solution by means of a second order centered finite differences method in space and a classical explicit-implicit Crank-Nicholson integration in time, assuming periodic boundary conditions. The computational domain Ω × [0, T ], for some Ω ⊂ R, is divided in cells of side ∆x × ∆t, where ∆x denotes the space step, and the time step is the same used for the ODE part. We denote with n Ω the total number of the spatial grid points {ω 0 , . . . , ω n Ω }. In greater details, denoting with u k l the approximation of u at the grid point (y l , t k ), the numerical scheme reads (36)

u k+1 l -u k l ∆t = D u k l-1 -2u k l + u k l+1 (∆x) 2 + 1 2 e κ(k+1)∆t N j=1 χ(y l -x k+1 j ) + 1 2 e κk∆t N j=1 χ(y l -x k j )
The coupling between the particle dynamics and the chemotaxis equation is realized by means of the gradient (the spatial derivative) of φ, evaluated, at each time instant t k , in the position of agent i at that time. The values of the gradient at grid points is approximated with centered difference, i.e.

∂ l φ k ≈ φ k l+1 -φ k l-1 2∆x
. Since x k i does not necessarely correspond to one of the grid point y l , we approximate ∂ x φ at time t k with a linear interpolation of the values of the gradient at the nearest grid points, in the following denoted as Φ x k i (∂ l φ k ). Recalling the choice of γ performed in [START_REF] Carrillo | Critical thresholds in 1D Euler equations with nonlocal forces[END_REF], choosing R = 1, the numerical scheme implemented for (P ) reads

                 v k+1 i = v k i + ∆t N N j=1 v k+1 j -v k+1 i 1+ ∥ x k i -x k j ∥ 2 β + η∆tΦ x k i (∂ l φ k )
x k+1 i = x k i + ∆tv k+1 i φ k l = e -κk∆t u k l with u k l approximated as in [START_REF] Piccoli | Control to flocking of the kinetic Cucker-Smale model[END_REF] (37) Kinetic level: numerical scheme for (V) Let now focus on system (V). The computational domain Ω × [0, T ] is divided in cells of side ∆x × ∆v × ∆t where ∆x, ∆v, are the space and velocity steps and ∆t is the time step. Let us assume that T is a multiple of ∆t to avoid rounding.

We choose Ω as a rectangular domain of size n x ∆x × n v ∆v, with n x , n v ∈ N the number of inner point discretizing the space and velocity domain. The generic point of the 2-dimensional grid is hence denoted as (x i , v j ) with i = 0, ..., n x , j = 0, ..., n v . We consider periodic boundary condition with respect to the x variable, and null condition with respect to the v variable. For any time instant t k we denote the approximation of ρ and ν at (x i , v j ) with ρ k i,j and V k i,j , respectively. The first-order scheme implemented reads: |V k i,j | ≥ 0 must be satisfied.

ρ k+1 i,j -ρ k i,j ∆t = - v j 2∆x ρ k i+1,j -ρ k i-1,j + |v j | 2∆x ρ k i+1,j -2ρ k i,j + ρ k i-1,j - 1 2∆v ρ k i+1,j V k i+1,j -ρ k i-1,j V k i-1,j + 1 2∆v ρ k i+1,j V k i+1,j -2ρ k i,j V k i,j + ρ k i-1,j V k i-1,j (38)
As at particle level, we approximate the solution of the chemotaxis equation by means of a second order centered finite differences method in space and a classical explicitimplicit Crank-Nicholson integration in time. The only difference is given by the source term. We approximate the two-dimensional integral in (21) using a quadrature formula. In greater details, it holds: Concerning the coupling through the gradient, the approximation of ∂ x ψ is required at spatial grid point. We approximate its value at each time instant using centered finite difference

∂ i ψ k ≈ ψ k i+1 -ψ k i-1
2∆x . Macroscopic level: numerical scheme for (E) We describe the numerical scheme for the most general Euler system simulated in this paper [START_REF] Menci | Existence of solutions for hybrid systems of differential equations under exogenous information with discontinuous source term[END_REF]. The scheme used for system [START_REF] Gamba | Percolation, morphogenesis, and burgers dynamics in blood vessels formation[END_REF] and ( 27) can be obtained simply considering ε = 0, η = 0, α = 0 and ε = 0, respectively.

Denoting W = W (x, t) = (µ t (x), µ t (x)u t (x)) T the vector of unknowns, density and momentum, we rewrite the hyperbolic part as [START_REF] Villani | Topics in Optimal Transportation[END_REF] W t + A(W ) x = S(W ),

where A and S are defined as follow: We approximate the solution at each time step starting from the relaxation tecniques originally proposed in [START_REF] Aregba-Driollet | Discrete kinetic schemes for multidimensional systems of conservation laws[END_REF]. In few words, the idea is to approximate hyperbolic system of equations through a BGK relaxation, which leads to a linear advection system with a relaxation source term. The advantage lies in the fact that at each time step t k , in order to find the solution at t k+1 , it is required to solve a system of linear, transport problems, which reduces the computational complexity.

We denote with I ⊂ R the computational domain. We approximate the solution of Euler system W k i = (µ k i , µ k i u k i ) T and of the parabolic equation ψ k i at each point of the grid x i = i∆x, i = 1, ..., n x -1. At boundary, we assume W 0 = W n x = (0, 0) T . The parabolic equation is approximated using the classical explicit-implicit Crank-Nicholson method, as previously explained. The algorithm requires to first solve a homogeneous system of equations and then a time integration of the source part is performed. Let now focus on the homogeneous part of [START_REF] Villani | Topics in Optimal Transportation[END_REF]. Using BGK approximation with two velocities of equal speed and opposite sign, i.e. λ 1 = λ = -λ 2 we denote with f r = (f µ r , f µu r ), r = 1, 2, the components of BGK model approximating W . For a complete description of the relaxation technique see [START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimensions[END_REF]. At each time iteration we solve Denoting with (f r ) k i , the approximated value of f r , at time t k at grid point x i , we compute (f r ) k+1/2 i running the following numerical upwind scheme:

(f r ) k+1/2 i = (f r ) k i - ∆t 2∆x λ k r (f r ) k i+1 -(f r ) k i-1 + ∆t 2∆x |λ k r | (f r ) k i+1 -2(f r ) k i + (f r ) k i-1
for any i = 1, ..., n x , r = 1, 2, with λ k 1 = max i |u k i | + P ′ (µ k i ) , λ k 2 = -λ k 1 . The stability of the scheme is guaranteed choosing, at each time step, ∆t k = 0.9 ∆x

λ k 1 .
The solution W at time t k+1 is obtained adding the discretization of the source term S(W ). The integral term and the chemotaxis gradient are treated explicitly, whereas the damping term in implicit, due to stiffness problem. The approximated value of the integral term, in the following denoted as I k i , is obtained by a quadrature formula, which has to be solved at each iteration. Finally we get (46)

W k+1 i = (f 1 ) k+1/2 i + (f 2 ) k+1/2 i + ∆t 0, µ k i I k i + ηµ k i ψ k i+1 -ψ k i-1 2∆x
.

1

  through this paper we define Lip(f ) for f : R n → R m , m, n ∈ N, as Lip(f ) := m i=1

Figure 1 .

 1 Figure 1. Test 1: Numerical simulation of Cucker-Smale model with β = 0.05, at particle level (first line) and kinetic (second line) level.

Figure 2 .

 2 Figure 2. Test 2: Numerical simulation of Cucker-Smale model with β = 0.95, at particle level (first line) and kinetic (second line) level.

Figure 3 .

 3 Figure 3. Test 3: Numerical simulation of Cucker-Smale model with chemotaxis at particle level (first line) and kinetic (second line) level, with parameters as in Test 2 and η = 1.4.

Figure 4 .

 4 Figure 4. Test 4: Numerical simulation of (22) at particle level (first line) and kinetic (second line) level, with parameters as in Test 3.

Figure 5 .

 5 Figure 5. Test 5: Numerical simulation of (24)-(26) for c 2 = 0.2 (first line) and c 2 = 0.5 (second line).

Figure 6 .

 6 Figure 6. Test 6: Numerical simulation of (27) with initial condition (25)-(26). First line: c 2 = 0.2, η = 1, α = 0; second line: c 2 = 0.5, η = 1, α = 0; third line c 2 = 0.5, η = 0, α = 1.

Figure 7 .

 7 Figure 7. Test 7, Test 11: Initial condition a) ρ 0 in the phase space, b) initial density µ 0 and c) initial momentum Q 0 .

Figure 8 .

 8 Figure 8. Test 7: approximation of monokinetic initial data.Comparison between zero-th order moment of the solution of (V) and µ t (second column), first order moment of the solution of (V) and Q t (third column): without chemotaxis (first line: η = 0, α = 0 at t = 2), with chemotaxis (second line: η = 0.2, α = 0 at t = 2), with chemotaxis and damping (third line: η = 0.2, α = 2 at t = 2).

Figure 9 .

 9 Figure 9. Test8, Test9, Test10: a) initial density µ 0 and b) initial momentum Q 0 (here null).

Figure 10 .

 10 Figure 10. Test 8: Comparison between zero-th order moment of the solution of (V) and µ t , first order moment of the solution of (V) andQ t at a)-b) t = 1, c)-d) t = 2, e)-f) t = 3.Here η = 0 (no chemotaxis), α = 0.

Figure 11 .Figure 12 .

 1112 Figure 11. Test 9: Comparison between zero-th order moment of the solution of (V) and µ t , first order moment of the solution of (V) and Q t at a)-b) t = 0.5, c)-d) t = 1, e)-f) t = 7. Here η = 3, α = 0.

Figure 13 .

 13 Figure 13. Test 10: Comparison between zero-th order moment of the solution of (V) and µ t , first order moment of the solution of (V) and Q t a)-b) t = 0.5, c)-d) t = 2, e)-f) t = 3. Here η = 3.4, α = 0 and the integral term in (E) is negletted.

Figure 14 .

 14 Figure 14. Test 11: Screenshots of numerical simulation of[START_REF] Menci | Existence of solutions for hybrid systems of differential equations under exogenous information with discontinuous source term[END_REF] for different values of ε. Comparison between zero-th order moment of the solution of (V) and µ t , first order moment of the solution of (V) and Q t at a)-b) t = 1, c)-d) t = 1.5, e)-f) t = 2. Initial condition and parameters are the same of Test 7.

Figure 15 .Figure 16 .Figure 17 .

 151617 Figure 15. Test 12, Test 13: a) initial density µ 0 and b) initial momentum Q 0 .

  ρ s (y, ξ)dydξ ≈ m,n: x m ∈[x i -R∆x,x i +R∆x], n=0,...,n v ρ k m,n ∆x∆v

  • -y, u(•) -u(y))µ(y)dy + ηµ∂ x ψ -αµu

2 W 1 2W

 21 (x, 0) + A(W (x,0)) λ (x, 0) -A(W (x,0)) λ  

Table 1 .

 1 Value of E k 0 and E k 1 for different values of ε.

		09	0.08	0.16
	t k = 1.5	0.23	0.18	0.32
	t k = 2	0.75	0.67	0.74
	E k 1	ε 1 = 0.002 ε 2 = 0.01 ε 3 = 0.05
	t k = 0	0	0	0
	t k = 1	0.06	0.05	0.23
	t k = 1.5	0.32	0.20	0.43
	t k = 2	0.99	0.85	0.95

Table 3 .

 3 4.146 0.034 t k = 2 4.235 0.533 t k = 4 4.236 0.588 Computation of ϵ * k and corresponding value of E k 0 . Time instants correspond to the ones selected for snapshots in Figure 16.

		ε * k	E k 0 (ε * k )
	t k = 0	0	0
	t k = 1 3.891 0.154
	t k = 2 3.082 0.094
	t k = 4 2.639 0.178

Table 4 .

 4 Computation of ϵ * k and corresponding value of E k 0 . Time instants correspond to the ones selected for snapshots in Figure17.

(M.M) Università Campus Bio-Medico di Roma Via Àlvaro del Portillo 00128 Roma, Italy Email address: m.menci@unicampus.it (R.N) IAC, Via dei Taurini, 19, 00185 Roma RM Italy Email address: roberto.natalini@cnr.it (T.P.) CNRS Laboratoire Ypatia des Sciences Mathématiques (LYSM), Roma, Italia & Sorbonne Université, CNRS, Université Paris Cité, Laboratoire Jacques-Louis Lions (LJLL), F-75005 Paris, France

Email address: thierry.paul@upmc.fr