Ghais Zriki
email: ghais.zriki@cefe.cnrs.fr

Rumsaïs Blatrix
email: rumsais.blatrix@cefe.cnrs.fr

Lise Roy
email: lise.roy@univ-montp3.fr

Predation interactions among henhouse-dwelling arthropods, with a focus on the poultry red mite Dermanyssus gallinae

Predation interactions among henhouse-dwelling arthropods, with a focus on the poultry red mite Dermanyssus gallinae Running title: Predation interactions involving Dermanyssus gallinae in poultry farms

Ghais Zriki, Rumsais Blatrix, Lise Roy

INTRODUCTION

The poultry red mite (PRM) Dermanyssus gallinae (De Geer) is the most important pest of laying hens worldwide. [START_REF] Chauve | The poultry red mite Dermanyssus gallinae (De Geer, 1778): current situation and future prospects for control[END_REF][START_REF] Sparagano | Prevalence and key figures for the poultry red mite Dermanyssus gallinae infections in poultry farm systems[END_REF] PRM is a nidicolous ectoparasite attacking resting hens at night. After relatively rapid blood meal, it retires to hide in different microhabitats like cracks, crevices, and dry droppings in the farm building. [START_REF] Nakamae | The new parasite ecology of chicken mites Dermanyssus gallinae parasitizing and propagating on chickens even in the daytime[END_REF]4 PRM has significant effect on the health and welfare of hens, it can cause anemia, decreased egg production and increased hen mortality. 5,6 The typical conventional control of PRM by means of synthetic acaricides is often not sufficient. In addition, the use of synthetic products has become increasingly reduced by stricter legislation regarding active ingredients. 7,8 Therefore, different alternative methods of control have been developed such as plant-derived product, 9 inert substances such as diatomaceous earth and silica, 10 electronic perches, biological control by means of natural enemies like entomopathogenic fungi 11,12 and predatory mites, 13,14 and research on vaccines is making progress. 15,16 The particular lifestyle of PRM as a nidicolous parasite living in a diversity of habitats in poultry farms makes this mite less likely to be reached by chemical treatment and more accessible for arthropod predators than ectoparasites living on the host. This suggests that PRM is an ideal target of biological control. Although biological control is well developed and has been successfully adopted to control pest arthropods in crop farming, this method has more recently begun to be developed against pests in livestock production and is still in its infancy. [START_REF] Chauve | The poultry red mite Dermanyssus gallinae (De Geer, 1778): current situation and future prospects for control[END_REF]7 To date, five predatory mites were experimentally shown to have potential in PRM biological control: Cheyletus eruditus (Schrank) and C. malaccensis (Oudemans), Androlaelaps casalis (Berlese), Hypoaspis aculeifer (Canestrini), and Stratiolaelaps scimitus (Womersley). 13,14,[17][18][START_REF] Toldi | Life cycle of the predatory mite Cheyletus malaccensis (Acari: Cheyletidae) fed on Poultry Red Mite Dermanyssus gallinae (Acari: Dermanyssidae)[END_REF] The two predatory mites A. casalis and C. eruditus are currently mass-reared and used as biocontrol agents of PRM in laying poultry farms, 7,[START_REF] Knapp | Use of predatory mites in commercial biocontrol: current status and future prospects[END_REF] but their efficiency in the field should be improved by complementary research. 18,[START_REF] Knapp | Use of predatory mites in commercial biocontrol: current status and future prospects[END_REF] Arthropods associated to poultry production [START_REF] Geden | Predation by predators of the house fly in poultry manure: Effects of predator density, feeding history, interspecific interference, and field conditions[END_REF][START_REF] Brady | The Mites of Poultry Litter: Observations on the Bionomics of Common Species, with a Species List for England and Wales[END_REF][START_REF] Proctor | Feather mites (Acari: Astigmata): ecology, behavior, and evolution[END_REF][START_REF] Walter | Mites in Soil and Litter Systems[END_REF][25] include several predatory taxa some of which are known to feed on various prey mite species. 13,25 However, the ability of the majority of these predators to feed on PRM has not been investigated. Amongst the naturally-occurring arthropod predators in poultry farms, various taxa are known to dwell in poultry manure. [START_REF] Geden | Predation by predators of the house fly in poultry manure: Effects of predator density, feeding history, interspecific interference, and field conditions[END_REF][START_REF] Walter | Mites in Soil and Litter Systems[END_REF][25][26] Several are typically active hunters that have been recurrently observed into microhabitats other than this substrate, including in PRM traps. 13 This makes them promising agents for PRM control in layer farm buildings.

The effect of predator communities on their prey's population may depend on direct and indirect interactions between these predators like competition and intraguild predation. [START_REF] Radloff | Large predators and their prey in a southern African savanna: a predator's size determines its prey size range[END_REF][START_REF]freshly fed preys, white + grey for unfed preys)[END_REF] Intraguild predation is very common among generalist predators that exploit common food resources [START_REF] Radloff | Large predators and their prey in a southern African savanna: a predator's size determines its prey size range[END_REF]29 which could significantly affect the dynamics of their shared preys. [START_REF]freshly fed preys, white + grey for unfed preys)[END_REF]30,31 The effect of intraguild predation on prey regulation can be antagonistic, though many case reports do not confirm it, and it may even be synergistic. [32][33][34] When focal and alternative preys share the same predatory species, [START_REF] Holt | The ecological consequences of shared natural enemies[END_REF] the availability of alternative prey can also influence the behavior of the shared predator and can lead to either increased or decreased predation on focal prey. [START_REF] Desneux | Potential of an alternative prey to disrupt predation of the generalist predator, Orius insidiosus , on the pest aphid, Aphis glycines , via short-term indirect interactions[END_REF][START_REF] Chailleux | Natural enemy-mediated indirect interactions among prey species: potential for enhancing biocontrol services in agroecosystems: NEmediated indirect interactions among prey species for biocontrol in agroecosystems[END_REF] This depends on many factors including the relative size of prey populations and predator preferences.

Astigmatic mites are microbivores/detritivores and they are the most frequent and abundant taxa in poultry manure. 25,[START_REF] Brady | Litter Mites and their Effects on Poultry[END_REF] They may serve as main or alternative prey for many generalist predatory mites. [START_REF] Barker | Bionomics of Androlaelaps casalis (Berlese) (Acarina: Laelapidae) a predator of mite pests of stored cereals[END_REF][START_REF] Žďárková | Compatibility of Cheyletus eruditus (Schrank)(Acari: Cheyletidae) and Cephalonomia tarsalis (Ashmead)(Hymenoptera: Bethylidae) in biological control of stored grain pests[END_REF] Evaluating the preferences of potential predators of PRM between astigmatic mites (as possible alternative or competing prey) and PRM (as focal prey) is essential to predict the effect of these predators on pest regulation. 41 The physiological status of PRM in terms of duration since the last blood meal may have a substantial effect on predation. One can expect from the fresh blood meal either a facilitating effect on predation owing to the weakening of the highly extended cuticle (making it possibly easier to be penetrated by the predator's chelicera) and the substantial slowing of PRMs' movements (akinesis is observed quickly after feeding [START_REF] Entrekin | Aggregation of the Chicken Mite, Dermanyssus gallinae (Acari: Dermanyssidae)[END_REF], or an antagonistic effect owing to the oxidative stress produced by the ingestion of fresh blood and/or the toxic products of its digestion (see adaptation mechanisms in hematophagous arthropods to dealing with feeding on fresh blood [START_REF] Graca-Souza | Adaptations against heme toxicity in blood-feeding arthropods[END_REF][START_REF] Lehane | The biology of blood-sucking in insects[END_REF] ). The physiological status of prey in interaction with other factors (species, sex, and size of the proposed preys) was shown to be an important driver of prey selection by a spider predator that feeds on mosquitoes. [START_REF] Jackson | A spider that feeds indirectly on vertebrate blood by choosing female mosquitoes as prey[END_REF][START_REF] Jackson | Convergence between a mosquito-eating predator's natural diet and its prey-choice behaviour[END_REF] Lastly, adult individuals of certain prey species have a greater ability to escape and better defense responses against predators than juveniles. [START_REF] Pastorok | Prey Vulnerability and Size Selection by Chaoborus Larvae[END_REF][START_REF] Sabelis | Predatory arthropods, Natural Enemies: The population biology of predators, parasites and diseases[END_REF] Prey consumption by several phytoseiid mite predators was inversely related to prey size. [START_REF] Fadamiro | Functional responses and prey-stage preferences of three species of predacious mites (Acari: Phytoseiidae) on citrus red mite, Panonychus citri (Acari: Tetranychidae)[END_REF][START_REF] Blackwood | Prey-Stage Preference in Generalist and Specialist Phytoseiid Mites (Acari: Phytoseiidae) When Offered Tetranychus urticae (Acari: Tetranychidae) Eggs and Larvae[END_REF][START_REF] Ganjisaffar | Prey stage preference and functional response of the predatory mite Galendromus flumenis to Oligonychus pratensis[END_REF] Hence, predators may be more or less prone to feed on freshly fed PRM and/or on different stages of PRM in farms.

This study aimed to characterize the potential trophic interactions involving PRM, alternative preys like microbivorous mites and arthropod predators that usually share the same microhabitats with PRM or those prone to hunt this mite in such microhabitats. Our specific objectives were to (1) identify potential predators of PRM based on their ability to consume various forms of PRM (different stages and physiological statuses) when no physical barriers hinder the access of theses predators to their preys, (2) evaluate the effect of the presence of alternative preys like astigmatic mites on the predation on PRM, and (3) assess whether predation interactions can occur between PRM predators.

MATERIALS AND METHODS

Arthropod sources

Arthropod predators and microbivore mites

To maximize the diversity of arthropod predators to be tested, we sampled arthropods from several barn layer farms located in the Drôme department (Rhône-Alpes-Auvergne region, France). These farms were selected based on previous in farm surveys. 25 Barn layer farms have slatted flooring under which manure is allowed to accumulate over long periods (flock duration = ca. 1 year) which permits an important development and establishment of manure-dwelling arthropods. [START_REF] Lay | Hen welfare in different housing systems[END_REF] We focused on 13 taxa of arthropods including 12 manure-dwelling taxa and one taxon of spiders. Certain mite taxa were multispecific and others were monospecific (Table 1). These arthropods were not reared or maintained on any transitional diet in the lab. Manure samples were kept into plastic containers covered with nylon-filter lids (mesh size 80 µm) before the extraction of arthropods.

Manure-dwelling arthropod individuals were extracted by dry sieving of manure samples using a series of stacked sieves with decreasing mesh size (from 1000 µm to 180 µm). In order to identify and differentiate mite taxa, we used the definition of morphospecies in Roy et al. 2017. 25 After performing predation tests on arthropods, the following taxa were identified at the species level: Pseudoscorpionida, Heteroptera (Insecta), Macrocheles spp. (Mesostigmata). For other morphoespecies, the taxonomic level was simply the level discernible under the stereomicroscope according to Roy et al. 25 (species level for monospecific morphoespecies, higher levels for others; see Table 1).

PRM

PRM aggregates were collected from two farm buildings in sealable plastic bags and kept fasting for one to three weeks in an incubator at 15 °C before tests. To produce freshly fed PRMs, fasted individuals were introduced into a PVC cylindrical container (60 L, 40 cm in diameter) with a chick for 2 h at 25 ± 5 °C in complete darkness. The top of the container was sealed with mite-proof nylon mesh (100 × 100 µm, PE171.6, Diatex, France). Pieces of folded paper were put in the container to provide shelters for PRMs to aggregate after the blood meal and facilitate the collection of fed individuals.

Experimental arenas for predation tests

Predation tests were conducted following the protocol by El Adouzi et al. [START_REF] El Adouzi | Detecting pyrethroid resistance in predatory mites inhabiting soil and litter: an in vitro test[END_REF] In short, each predator was tested in an experimental arena constituted by a well (flat-bottomed, 7 mm diameter) of a transparent polystyrene microplate (Nunc™ 167008 F96 MicroWell 96-well × 400 µL, external dimensions: 128 × 86 mm Cell Culture Microplate, Fyn, Denmark). We added 2-4 µl of 1.5 % agarose gel into each well to prevent dehydration of the arthropods, a technique validated by El Adouzi et al. [START_REF] El Adouzi | Detecting pyrethroid resistance in predatory mites inhabiting soil and litter: an in vitro test[END_REF] Microplates were covered with stretched plastic paraffin film (Parafilm®, Bemis Co., USA).

Predation test experiment

A test consisted of confining one single predator and one or two preys (depending on the modality, see below) into a well of a 96-well microplate, wells of which were used as replicated experimental arenas. After a fixed contact duration, prey mortality was recorded in each arena. Controls consisted in the same single or paired prey(s) isolated into wells on the same microplate without any predator (same number of wells for controls as for tests with predators). Several different modalities of predator/preys(s) combination were tested on the same microplate and at the same time (= a series). To minimize the effect of random factors, (1) each predator modality was tested on different microplates successively, (2) two to four different modalities of predator/prey(s) combination were tested together in each series, (3) the set of modalities to be tested together was randomly rearranged for each series, so as to be free of dependencies between modalities. Each modality was replicated dozens of times on a microplate (one predatory individual per taxon and per well on 2 to 4 columns, ie 16 to 32 individual tests, and a similar number of control wells), and this was repeated two to three time in different series.

Predators were fasted for 24 h before being tested. Microplates were maintained in a climatic chamber at 25.0 ± 0.5 °C in complete darkness.

Prey condition in each experimental arena was examined under a stereomicroscope after 24 hours of test duration. A prey was deemed to be dead if no movement was triggered by contact with a thin paintbrush.

Arachnid predators do not swallow prey but only suck internal body liquids, [START_REF] Koehler | Predatory mites (Gamasina, Mesostigmata), Invertebrate Biodiversity as Bioindicators of Sustainable Landscapes[END_REF] making it difficult to differentiate predation events from prey natural mortality. Thus, in each series, prey mortality in predation tests was corrected by prey natural mortality recorded in the corresponding controls to obtain the predation rate.

Predation test types

Prey-choice tests

Prey choice tests were used to evaluate (1) the predation rate of each predator on each prey species and (2) the preference of predators between these offered preys. The predation rate (1) was estimated as the frequency of predation on a given kind of prey. Two prey individuals with (test) or without (control) one predator individual were introduced per well.

The preference (2) was tested between two PRM developmental stages (unfed juveniles and unfed adult females), or two physiological statuses of PRM (freshly fed and unfed), or between unfed PRM and another prey mite, with a focus on Astigmatic mites. The preference according to the physiological status and between PRM and Astigmatic mites was tested both with PRM juvenile and adult females. All in all, twelve putative predatory taxa were tested in 31 modalities of prey-choice test (Table 2).

Predator-to-predator interactions

We evaluated whether six native predators engaged in predator-prey interaction when no shared prey was offered, following five combinations:

(1) A. casalis vs Cheyletus spp., (2) A. casalis vs M. penicilliger, (3) A. casalis vs Lam. nodosus, (4) Lam. nodosus vs M. penicilliger, and (5) Dendrolaelaps spp. vs Pro. parascolyti. The choice of these combinations was driven by the availability of predators in the successive sampling campaigns. For each predator pair, one individual of each of the two species (test) or a single individual (control) were introduced per well.

Note that since each of the two predators confronted in the test wells is also a prey, two controls (one for each taxon) are associated with each test. Combinations one to five were tested 39, 24, 24, 22 and 24 times respectively.

Complementary tests

We performed a few complementary tests with the same protocol as preychoice tests (2.4.1) to document trophic relationships between tested arthropods excluding PRM. Five predatory taxa and microbivore mites were subjected to predation tests according to different combinations described in Table 2. (Modalities 32 to 37).

Data analysis

A predation test replicate was discarded from analyses if any of predator(s) or prey(s) had escaped from the experimental arena, or if the predator was dead at the moment of rating. In the case of predator-topredator interaction tests, the replicate was discarded when the two predators were dead at the time of rating. In addition, In addition, a test modality was considered invalid in a series where natural mortality of prey(s) (or predators in the case of predator-to-predator interactions), i. e. in controls, exceeded 15% or differed significantly between the two prey species (or the two predators in predator-to-predator interactions; we used the Chi-squared test to check for independence of natural morality rates between the two prey species).

To represent the potential trophic relationships between tested arthropod taxa, we considered the frequency of predation per putative predator on each prey species from all tests. The predation rate for each predator-prey pair in each test series was estimated by correcting prey mortality with natural mortality determined in the corresponding control (on the same series) using the Abbott formula [START_REF] Abbott | A method of computing the effectiveness of an insecticide[END_REF] . The predation rate on each prey species and for each predator was calculated as the average corrected prey mortality over all test series. Predation rates on PRM were calculated by considering all PRM stages and physiological statuses tested indiscriminately.

To assess preferences of predators between prey species, Chi-squared tests were applied to test for dependence between prey mortality and prey species in the presence of each predator. In prey-choice tests, to consider the three possible configurations (both prey dead, one prey dead, no prey dead), 2 x 2 contingency tables were constructed as follows: prey A dead or live x prey B dead or live. In tests of predator-to-predator interactions, only tests where one and only one predator had died were included in the analysis and 2 x 2 contingency tables were constructed as follows: predator A or B x live or dead. All chi-squared tests were conducted within the R environment, [START_REF] Team | R: A language and environment for statistical computing[END_REF] P values were calculated by Monte Carlo simulation with 2000 replicates

Ethical requirements

All the experiments involving hens were conducted in compliance with regulations on animal experimentation (reference number of the Ethics committee: 036; project number: APAFIS#2339-2015101122029640 v4).

RESULTS

By bringing together the results of all type of predation tests, we provide an insight into the potential trophic interactions among selected henhouse-dwelling arthropods, showing that PRM could potentially be part of this partial food web. (Fig. 1).

Predation on PRM

In our experiments, ten indigenous arthropod taxa fed on PRM, with predation rates (Abbott-corrected prey mortality percent) ranging from 100 to 4%: M. penicilliger, Lam. nodosus, Cheyletus spp., M. muscaedomesticae, Lyc. campestris, linyphiid spider, Par. fimetorum, A. casalis, U. fimicola, and Dendrolaelaps spp. (in order of decreasing predation rates, Fig. 1). The last two did so only occasionally (8 and 4% predation rates, respectively). In the presence of Uropodina spp. or of Pro. parascolyti, PRM mortality did not exceed natural mortality in the corresponding predator-free controls.

Preference of putative predators for PRM or astigmatic mites

Of the seven arthropods tested as predators, four fed significantly more on PRM than on astigmatic mites in our conditions (Cheyletus spp.: ꭕ 2 = 27, P = 4.9 × 10 -4 ; A. casalis: ꭕ 2 = 13, P = 9.9 × 10 -4 ; M. muscaedomesticae: ꭕ 2 = 21, P = 4.9 × 10 -4 ; U. fimicola: ꭕ 2 = 8.5, P = 7.9 × 10 -3 ; Fig. 2). Neither Cheyletus spp. nor U. fimicola fed on astigmatic mites in the presence of PRM juveniles.

Effect of stage and physiological status of PRM preys on predation preference

Cheyletus spp. and A. casalis fed significantly more on juveniles than on adult females (Cheyletus spp.: ꭕ 2 = 42, P = 4.9 × 10 -4 ; A. casalis: ꭕ 2 = 12, P = 4.9 × 10 -4 ; Fig. 3). Predation by M. muscaedomesticae and the linyphiid spider did not differ significantly between the two stages.

The bug Lyc. campestris fed significantly more on freshly fed PRM juveniles than on unfed ones (ꭕ 2 = 8, P = 8.9 × 10 -3 , Fig. 4A). Cheyletus spp. fed significantly more on unfed adult PRM females than on freshly fed ones (ꭕ 2 = 17, P = 9.9 × 10 -4 ; Fig. 4B).

Predator-to-predator interactions

Predation occurred within all pairs of predators tested in our conditions (Fig. 5). Predation was significantly asymmetric between the pairs M. penicilliger-A. casalis (ꭕ 2 = 48, P = 4.9 x 10 -4 ), and Lam. nodosus-A. casalis (ꭕ 2 = 13.1, P = 9.9 x 10 -4 ). A. casalis did not kill any of the two other predators. It was more balanced and insignificant in the pair A. casalis-Cheyletus spp., but mortality was twice as high in Cheyletus spp. as in A.

casalis. Predation interactions in the other pairs appear to be accidental.

DISCUSSION

Potential food web involving PRM

Although the predation interactions we observed occurred in the absence of any physical barrier between predator and prey, our results suggest that PRM in barn layer farms could be embedded in a food web containing at least ten arthropod taxa able to feed on PRM. Eight of these taxa are recorded to feed on PRM for the first time here: Lyc. campestris, Lam. nodosus, a linyphiid spider, M. penicilliger, M. muscaedomesticae, Par. fimetorum, U. fimicola, and Dendrolaelaps spp. The last two did so only occasionally. We also confirmed predation on PRM by natural populations of Cheyletus spp. and A. casalis which are actually commercially available for release in layer buildings as biocontrol agent of PRM, [START_REF] Knapp | Use of predatory mites in commercial biocontrol: current status and future prospects[END_REF] although in our experimental design predation by A. casalis can be considered moderate (19% predation rate).

Predators' preference for different forms of PRM

Among the four tested arthropod taxa that showed substantial predation on PRM (predation rate > 10%), the largest-sized taxa (M. muscaedomesticae and the linyphiid spider) did not feed differentially on juveniles and adult PRM, whereas the smallest-sized (Cheyletus spp. and A. casalis) preferred juveniles, consistent with previous findings that larger predators use a wider range of prey sizes. [START_REF] Cohen | Body sizes of animal predators and animal prey in food webs[END_REF]58 Interestingly, PRM consistently form multi-layered clusters with juveniles remaining in the center and adult females staying at the top. [START_REF] Entrekin | Aggregation of the Chicken Mite, Dermanyssus gallinae (Acari: Dermanyssidae)[END_REF] This arrangement was considered a protective behavior to preserve juveniles from predators. Furthermore, predators with different prey-stage preferences do not affect the population dynamics of prey and predator in the same manner. [START_REF] Gurtin | On predator-prey interactions with predation dependent on age of prey[END_REF] Mathematical modelling of prey-predator population dynamics suggested that only predators with no prey-stage preference allow Lotka-Volterra periodic prey-predator oscillations to occur and be maintained. [START_REF] Gurtin | On predator-prey interactions with predation dependent on age of prey[END_REF] Predators that prefer juvenile stages induce a rapid increase in the prey-predator oscillations' amplitude, which ultimately results in the extinction of both predator and prey. [START_REF] Gurtin | On predator-prey interactions with predation dependent on age of prey[END_REF] Selective predation focused on young preys is thus more likely to lead to extinction than non-selective predation. [START_REF] Gurtin | On predator-prey interactions with predation dependent on age of prey[END_REF][START_REF] Hastings | Age-dependent predation is not a simple process. I. Continuous time models[END_REF] The predatory bug Lyc. campestris was the only tested taxon that fed more on freshly fed juveniles than on unfed ones. This is perhaps because of the slower movement owing to blood ingestion, as preference for slow moving preys over faster ones was reported in other predatory bugs. 41 An alternative explanation (though not exclusive) is that Lyc. campestris could be prone to feed indirectly on vertebrate blood, as does an African jumping spider. [START_REF] Jackson | A spider that feeds indirectly on vertebrate blood by choosing female mosquitoes as prey[END_REF] The preference of Cheyletus spp. for unfed females over freshly fed ones may be explained by the size and/or the motility of prey instead of deleterious effects of the blood meal as these predators did not show any preference for unfed/fed juveniles. In fact, individuals of predator taxa for which no preference for unfed/fed PRM (juveniles or adults) was detected frequently consumed both preys (fed and unfed PRM), showing no aversion to taking a potentially deleterious meal of vertebrate fresh blood. However, in these taxa, we may have missed preferences between unfed and fed PRM by hungry predators that may have recurrently chosen a particular status in the first act of predation.

Implications of interactions other than predator-PRM

The most abundant poultry mite taxon that can constitute an alternative prey to PRM predators is Astigmata. A significant preference for PRM over astigmatic mites was observed in A. casalis, Cheyletus spp., Lam. nodosus, and M. muscaedomesticae. The former two are known to feed and successfully develop on several species of astigmatic mites. [START_REF] Brady | Litter Mites and their Effects on Poultry[END_REF][START_REF] Barker | Bionomics of Androlaelaps casalis (Berlese) (Acarina: Laelapidae) a predator of mite pests of stored cereals[END_REF][START_REF] Silva Gl Da | Mite fauna (Acari) associated to commercial laying hens and bird nests in Vale do Taquari, Rio Grande do Sul, Brazil[END_REF][START_REF] Solomon | Experiments on predator-prey interactions of storage mites[END_REF][START_REF] Cebolla | Prey range of the predatory mite Cheyletus malaccensis (Acari: Cheyletidae) and its efficacy in the control of seven stored-product pests[END_REF] Cheyletus spp. is a sit-and-wait predator ("ambush" foraging mode) and it rarely actively stalks its prey. [START_REF] Wharton | Predatory behaviour of the mite cheyletus aversor[END_REF] As astigmatic mites are very slowly moving, [START_REF] Brady | Litter Mites and their Effects on Poultry[END_REF] these mites are less likely to encounter the cheyletid predator than the fast moving unfed PRM individual in the experimental arenas. This could partially explain the complete absence of predation on astigmatic mites by this predator in our experiments. As our results show that intraguild predation is possible between several of the PRM predators, the natural regulation of PRM populations might be affected by the predator diversity in farms. Geden et al. [START_REF] Geden | Predation by predators of the house fly in poultry manure: Effects of predator density, feeding history, interspecific interference, and field conditions[END_REF] reported the negative impact of intraguild predation between two naturally-occurring predatory mites, Parasitus sp. and M. muscaedomesticae, on the control of the house fly.

CONCLUSION AND PERSPECTIVES

Prior to the present study, information about native predators of PRM that have been recorded from poultry farms to date was limited to three species A. casalis , C. eruditus and C. malaccensis. 13,14,[START_REF] Toldi | Life cycle of the predatory mite Cheyletus malaccensis (Acari: Cheyletidae) fed on Poultry Red Mite Dermanyssus gallinae (Acari: Dermanyssidae)[END_REF] Our study revealed that six additional taxa native to poultry farms are substantial predators of PRM, all with predation rates on PRM higher than that of A. casalis. Among these, two (Par. fimetorum and M. penicilliger) are only found occasionally in poultry farms, and thus, do not seem to establish sustainably in this environment. The other four (the bug Lyc. campestris, the pseudoscorpionid Lam. nodosus, the mite M. muscaedomesticae, and the linyphiid spider) might be worth considering as additional candidates for augmentative biological control. They could also provide an unknown regulating ecosystem service that could be worth promoting via conservation biological control practices. Intricate interactions of varying degrees between arthropods in poultry ecosystems likely make actual food-webs different from the present reconstruction, as inferred from barrier-free tests. Ecosystem services provided by pest enemies are strongly dependent on spatial heterogeneity and layout coupled with the enemies' dispersal capabilities. [START_REF] Collard | Spatial organisation of habitats in agricultural plots affects per-capita predator effect on conservation biological control: An individual based modelling study[END_REF][START_REF] Rusch | Flow and stability of natural pest control services depend on complexity and crop rotation at the landscape scale[END_REF] Special attention should be paid to spatial mapping of prey-predator meeting points in henhouses and identification of their determining factors to anticipate impediments to biological control owing to limited penetration by predatory mites into pest-infested areas. 67

15

  Bartley K, Turnbull F, Wright HW, Huntley JF, Palarea-Albaladejo J, Nath M, et al., Field evaluation of poultry red mite (Dermanyssus gallinae) native and recombinant prototype vaccines, Veterinary parasitology 244:25-34 (2017). 16 Price DR, Küster T, Øines Ø, Oliver EM, Bartley K, Nunn F, et al., Evaluation of vaccine delivery systems for inducing long-lived antibody responses to Dermanyssus gallinae antigen in laying hens, Avian Pathology 48:S60-S74 (2019). 17 Ali W, George DR, Shiel RS, Sparagano OAE, and Guy JH, Laboratory screening of potential predators of the poultry red mite (Dermanyssus gallinae) and assessment of Hypoaspis miles performance under varying biotic and abiotic conditions, Vet Parasitol 187:341-344 (2012). 18 Lesna I, Sabelis MW, van Niekerk TGCM, and Komdeur J, Laboratory tests for controlling poultry red mites (Dermanyssus gallinae) with Hinton JL and Moon RD, Arthropod populations in high-rise, cagedlayer houses after three manure cleanout treatments, J Econ Entom 96:1352-1361 (2003). 27 Polis GA and Holt RD, Intraguild predation: the dynamics of complex trophic interactions, Trends Ecol Evol 7:151-154 (1992). 28 Rosenheim JA, Wilhoit LR, and Armer CA, Influence of intraguild predation among generalist insect predators on the suppression of an herbivore population, Oecologia 96:439-449 (1993). 29 Rosenheim JA, Kaya HK, Ehler LE, Marois JJ, and Jaffee BA, Intraguild predation among biological-control agents: theory and evidence, Biological control 5:303-335 (1995). 30 Schausberger P and Walzer A, Combined versus Single Species Release of Predaceous Mites: Predator-Predator Interactions and Pest Suppression, Biol Control 20:269-278 (2001). 31 Croft BA and MacRae IV, Biological control of apple mites by mixed populations of Metaseiulus occidentalis (Nesbitt) and Typhlodromus pyri Scheuten (Acari: Phytoseiidae), Environ Entomol 21:202-209 (1992). 32 Rosenheim JA and Harmon JP, The influence of intraguild predation on the suppression of a shared prey population: an empirical reassessment, Trophic and Guild in Biological Interactions Control, Springer, Dordrecht, pp. 1-20 (2006). 33 Ong TWY and Vandermeer JH, Coupling unstable agents in biological control, Nat Commun 6:1-9 (2015).

9. 4

 4 Figure 4. Mortality (%) in D. gallinae with two physiological statuses. A: unfed and freshly fed juveniles; B: unfed and freshly fed females. 24-h contact period. n = number of tested predators. Black: predation on freshly fed prey, white: predation on unfed prey, grey: predation on both preys. Bars are labeled with * when mortality rates differ significantly between preys within the corresponding modality. Significant difference means predator preference for the prey with the highest mortality rate (black + grey for

ACKNOWLEDGEMENTS

We thank all the poultry farmers whose farms were sampled, Dr Jean-Claude Streito for identifying the bugs, and Dr Mark Judson for identifying the pseudoscorpions. In addition, we thank Tristan Gambin, Jordan Dijoux, and Dylan Tallon for expert lab technical assistance.

The FEADER (Fonds Européen Agricole pour le Développement Rural) European funds, the French Rhone-Alpes-Auvergne Region [grant number RRHA 160116CR0820011], as well as the CNPO (Comité National pour la Promotion de l'OEuf, France) supported this work. Ghais Zriki was supported via a PhD fellowship from the CeMEB LabEx and the French Occitanie Region.

Sikorska D, Garnis J, Dąbrowski ZT, Sikorski P, Gozdowski D, and Hopkins RJ, Thus far but no further: predatory mites do not migrate effectively into strawberry plantations, Exp Appl Acarol 77:359-373 (2019). 29

TABLES

Table1

Figure 3.

Mortality (%) of D. gallinae females and juveniles in the presence of single predatory individuals over a 24-h contact period. n = number of tested predators. Black: mortality of unfed females, white: mortality of unfed juveniles, grey: mortality of both preys. Bars are labeled with * when mortality rates differ significantly between preys within the corresponding modality. Significant difference means predator preference for the prey with the highest mortality rate (black + grey for females, white + grey for juveniles). (P < 0.05, Chi-squared tests).


Abdelghani El Mousaoui

Pierre Argoul

Mohammed El Rhabi

Abdelilah Hakim

Marrakech Abdelghani El Mousaoui

Mohammed El Paris

Rhabi

Abdelilah Rouen

Marrakech Hakim

Kinetic BGK model for a crowd: Crowd characterized by a state of equilibrium

Keywords: Discrete kinetic theory, Crowd dynamics, BGK model, Semi-Lagrangian schemes MSC 2010: 35A01, 35A02, 97M70, 97N40

Introduction

Mathematical representations of crowd motion from the microscopic to macroscopic scale have been an active field of study for the last three decades. An overview of the most important models at microscopic, macroscopic, and mesoscopic scale is reviewed in [START_REF] Cristiani | Multiscale Modeling of Pedestrian Dynamics[END_REF]. Indeed, the most popular crowd simulation models are the individual models, namely the heuristic rule-based models [START_REF] Reynolds | Steering behaviors for autonomous characters[END_REF], mechanical models [START_REF] Helbing | A mathematical model for the behavior of pedestrians[END_REF][START_REF] Helbing | Social force model for pedestrian dynamics[END_REF][START_REF] Helbing | Self-organization phenomena in pedestrian dynamics[END_REF], and cellular automate [START_REF] Nishinari | Extended floor field ca model for evacuation dynamics[END_REF], continuous models are based on fluid dynamics [START_REF] Helbing | A fluid dynamic model for the movement of pedestrians[END_REF][START_REF] Henderson | On the fluid mechanics of human crowd motion[END_REF][START_REF] Rinaldo | Existence of nonclassical solutions in a pedestrian flow model[END_REF], and the kinetic (Gas-kinetic) models are intermediate models between the two discrete and continuous models [START_REF] Hoogendoorn | Gas-kinetic modeling and simulation of pedestrian flows[END_REF][START_REF] Bouchut | On zero pressure gas dynamics[END_REF][START_REF] Bouchut | Duality solutions for pressureless gases, monotone scalar conservation laws, and uniqueness[END_REF]. Handerson was the first to apply this type of "kinetic gas" model to empirical pedestrian crowd data [START_REF] Henderson | The statistics of crowd fluids[END_REF][START_REF] Henderson | Sexual differences in human crowd motion[END_REF]. In 2011 Bellomo et al. [START_REF] Bellomo | On the modeling of crowd dynamics: Looking at the beautiful shapes of swarms[END_REF][START_REF] Bellomo | Modeling complex living systems: a kinetic theory and stochastic game approach[END_REF][START_REF] Agnelli | A kinetic theory approach to the dynamics of crowd evacuation from bounded domains[END_REF][START_REF] Bellomo | Toward a mathematical theory of behavioral-social dynamics for pedestrian crowds[END_REF][START_REF] Bellomo | Behavioral crowds: Modeling and monte carlo simulations toward validation[END_REF][START_REF] Burini | A kinetic theory approach to the modeling of complex living systems[END_REF] have developed the kinetic approach for crowds in a recent approach called, the kinetic theory for active particles. This approach considers the crowd as a complex system. The microscopic state of each particle is characterized by a geometric variable, position x = (x, y) and a mechanical variable velocity v = (v x , v y ). In addition there is a microscopic state related to their socio-biological behavior, called activity, it noted u. The representation of the system is defined by a distribution function noted f (t, x, v, u), where f (t, x, v, u)dxdv represents the number of active pedestrians who at the moment t are in the elementary volume [x, x+dx]× [v, v+dv] and who have activity u. Pedestrian movements are governed by the partial derivative equation (PDE) of transport applied to f , Γ characterizing the different interactions between pedestrians and their environment (1.1) ∂f (t, x, v, u) ∂t + v.∇ x f (t, x, v, u) = Γ[f ](t, x).

In our previous work [START_REF] Elmoussaoui | Discrete kinetic theory for 2d modeling of a moving crowd: Application to the evacuation of a non-connected bounded domain[END_REF], we considered the model (1.1). The term Γ[f ](t, x) which models the interactions between pedestrians with various obstacles, is treated from a probabilistic point of view.

In this paper, we are interested in one of the simplest ways to model the term Γ[f ](t, x). It consists of describing a pedestrian tendency to equilibrium similar to the BGK operator which replaces the collision operator of the Boltzmann equation. Specifically, the case where the pedestrian system is characterized by an equilibrium configuration f e and a relaxation term τ [ρ]. In an emergency evacuation case, pedestrians try to achieve a desired velocity noted v d to reach a target. τ [ρ] is the relaxation term describes the adaptation of the density f to the equilibrium density f eq (v). Therefore the interactions term takes the following simple form:

(1.2) Γ[f ](t, x) = 1 τ [ρ] (f eq (v → v d ) -f ) .
This paper develops a special theory for pedestrian motion and these interactions namely deceleration avoidance. Consequently the BGK model (1.1), (1.2) proposed in this work does not use the assumptions of conservation of momentum and energy. Only the conservation of mass must be verified in our study. The equilibrium state function of pedestrians is developed and based on Henderson works [START_REF] Henderson | The statistics of crowd fluids[END_REF][START_REF] Henderson | Sexual differences in human crowd motion[END_REF]. A mathematical framework for a theoretical study of the proposed model is determined. The rest of this paper is organized as follows: section 2

provides the discrete velocity model derived from the continuous BGK equation. This model describes the motion of pedestrians reaching an equilibrium configuration in a domain Ω. Then an approximation of the Maxwellian discrete density representing this state of equilibrium is presented. In section 3, a result of existence and uniqueness for the discrete velocity model is demonstrated. Then we prove the convergence of this solution to a solution of the continuous BGK equation. Section 4 is devoted to numerical simulations to validate the proposed model, and to show its ability to describe the main features of the dynamics of pedestrians. (ped/m 2 ). This group of N pedestrians present in the domain Ω at the initial time t 0 . They are in a normal traffic situation, i.e. they have the ability to direct themselves towards all possible directions. After, all pedestrians have the tendency of a comfortable destination noted v d (state of equilibrium). The group of pedestrians wants to reach a target located at the point x c ( see Fig. 1). The state of the crowd is represented by the density f (t, x, v). They move with a velocity v ∈ D v .

Mathematical model

The average crowd quantities obtained by integrating f in the velocity space D v :

• density:

ρ(t, x) = Dv f (t, x, v)dv,
• average velocity:

ξ(t, x) = 1 ρ(t, x) Dv vf (t, x, v)dv,
• total energy:

E = 1 2 Dv ||v|| 2 f (t, x, v)dv.
The evolution of the particle density f (t, x, v) is described by the following equation:

(2.1) ∂f (t, x, v) ∂t + v.∇ x f (t, x, v) = 1 τ [ρ] (f eq (v) -f ) .
The coefficient τ can depend on the density ρ(t, x), this term expresses that the distribution f would not go instantly to the desired velocity distribution f eq , but would need some time called relaxation time τ . For reasons of simplicity, we assume that this relaxation term is a constant i.e. τ [ρ] = τ . According to this hypothesis, the model (2.1) takes the following form:

(2.2) ∂f (t, x, v) ∂t + v.∇ x f (t, x, v) = 1 τ (f eq (v) -f (t, x, v)) .
We consider the model (2.1) with following initial data:

(2.3) f (t = 0, x, v) = f 0 (x, v), x ∈ Ω, v ∈ D v .
The system defined by the two equations (2.2), (2.3) represents the Bhatnagar-Gross-Krook model. In this paper it describes the temporal evolution of the distribution of particles (pedestrians). This model is less expensive than the Boltzmann equation because it is sufficient to update the macroscopic fields at each time step. On the other hand it provides qualitatively correct solutions for macroscopic moments. These two aspects, namely the relatively low cost of calculation and the correct description of the hydrodynamic limit, explain the interest in the BGK model during the last decades.

It also shares important features with Boltzmann's original equation, such as the following conservation laws:

conservation of mass Dv f eq (v)dv = Dv f (t, x, v)dv. (2.4)
Thus, the BGK equation is a kinetic collision equation that takes into account only the overall effect of pedestrian interactions.

Remark 1. The BGK model for pedestrian motion contains corrections due to interactions such as avoidance and deceleration. Therefore this model does not obey the conservation of momentum and energy. In our case, the only law of conservation which must be respected by the model is the conservation of the mass. It is expressed by the following equation:

(2.5) Dv f eq (v)dv = Dv f (t, x, v)dv.

2.2.

Maxwellian approximation: equilibrium density f eq (v). In 1971 L.F.Henderson [START_REF] Henderson | The statistics of crowd fluids[END_REF], suggested that the motion of people in a crowd represents a system similar to a gas molecule collection. Specifically, he suggested that the classical Maxwell-Boltzmann theory of a molecular system could also describe the velocity distribution of individuals.

Henderson measured the speed distribution function for 3 crowd categories: a crowd of university students walking from the library to the university, an adult crowd of all ages using a pedestrian crossing on a street, and children in a playground air. Analyses and estimates are made under certain numbers of assumptions about the crowd, namely, the movement is defined in every moment t by position (x, y) and the velocity v = (v x , v y ). All the individuals in the crowd have the same mass. The two figures 2 (a), (b) show that Henderson's empirical results agree with the classical Maxwell-Boltzmann theory. The distribution of the v x component of the velocity is given by the following equation:

(2.6) f eq (v x ) = 1 N dN vx dv x = 1 √ 2πv m exp Å - 1 2 v 2 x v 2 m ã ,
where,

• v m is the square root of the average module value of speed.

• N the total number of pedestrians.

• N vx is the number of pedestrians with speed v x .

In a similar way for the equation of the v y component distribution of velocity, Henderson found the result for the distribution of v = (v x , v y ), is given by the following equation :

(2.7)

f eq (v) = 1 N dN v dv = 1 2πv 2 m exp Å - 1 2 v 2 v 2 m ã ,
where N v is the number of pedestrians with velocity v. According to the above and the last remark, we define in our case the density of equilibrium for pedestrians, with the following formulation:

(2.8)

f eq (v, ρ) = ρ(t, x) 2πv 2 m exp Å - 1 2 ||v -v d || 2 v 2 m ã ,
where, v m = » 1 2 Dv ||v|| 2 f (t, x, v)dv; v d is the desired direction. f eq models the equilibrium state that each pedestrian wishes to achieve. Indeed the tendency of all pedestrians to reach a comfortable destination defined by the direction v d . The density of equilibrium that we have defined satisfies the conservation of the mass defined by the equation (2.5), indeed:

Dv f eq (v)dv = Dv ρ(t, x) 2πv 2 m exp Å - 1 2 ||v -v d || 2 v 2 m ã dv = ρ(t, x) = Dv f (t, x, v)dv.

2.3.

Model of discrete velocities. The number of pedestrians N is generally insufficient to justify the hypothesis of continuity of the particle distribution function f (t, x, v) with respect to velocity. Thus for numerical simulations, a discrete velocity approximation of the BGK equation is introduced. We refer to the discrete velocity models for the Boltzmann equation developed by Rogier and Schneider [START_REF] Rogier | A direct method for solving the boltzmann equation[END_REF], Buet [START_REF] Buet | A discrete-velocity scheme for the boltzmann operator of rarefied gas dynamics[END_REF], Heintz and Panferov [START_REF] Heintz | A new consistent discrete-velocity model for the bolzmann equation[END_REF], and Mieussens [START_REF] Mieussens | Convergence of a discrete-velocity model for the boltzmann-bgk equation[END_REF]. The proposed approximation in this work has the same conservation properties as the continuous BGK model.

Let K be a set of multi-indices of Z 2 , defined by

K = {k = (k 1 , k 2 ), |k| ≤ B}, with B is a scalar. We define V ⊂ R 2 , a set of N v discrete velocities, V = {v k = k∆v /k ∈ K},
where ∆v is a scalar. The distribution f of the continuous velocities is then replaced by the following

N v vector, f K = (f k (t, x)) k∈K where each component f k (t, x) is an approximation of the function f (t, x, v k ).
Thanks to the previous discretization of the velocity, we define the local density by the following equation:

(2.9)

ρ (t, x) = k∈K f k (x, t) .
The discrete kinetic model associated with the equation BGK (2.2), is defined by the set of the following equations:

(2.10)

∂f k (t, x) ∂t + v k .∇ x f k (t, x) = 1 τ (f eq,k (v k , ρ) -f k (t, x)) , k ∈ K v k ∈ V,
where f eq,k is an approximation of the equilibrium density defined by (2.8) .

The main problem is to define this approximation of the discrete density f eq,k such that the property of conservation of mass is satisfied. We used the natural approximation used by Yang and Huang [START_REF] Yang | Rarefied flow computations using nonlinear model boltzmann equations[END_REF] and which has been developed

by Luc [START_REF] Mieussens | Discrete-velocity models and numerical schemes for the boltzmann-bgk equation in plane and axisymmetric geometries[END_REF]:

(2.11) f eq,k (v k , ρ) = f eq (v k ), k ∈ K, hence, (2.12) f eq (v k , ρ) = ρ(t, x) 2πv 2 m exp Å - 1 2 ||v k -v d || 2 v 2 m ã ,
where

v m = 1 ρ(t, x) k∈K ||v k || 2 f k (t, x).
We considered the mathematical model (2.10) with an initial data defined by: (2.13)

f k (t = 0, x) = f 0,k (x) k ∈ K. x ∈ Ω.

Boundary condition on ∂Ω

In our study, we are interested in the adaptation of an equilibrium situation of the pedestrians inside the domain Ω. We assume that their target is inside the domain. The disk diameter 2r D , occupied by pedestrians is always less than the distance between the target in x c and the edge of the domain ∂Ω, i.e.

2r D < d (x c , ∂Ω) .
According to this hypothesis, the theoretical study of our problem is reduced to all plane R 2 in places Ω.

Theoretical study of the proposed mathematical model

Some important mathematical results concerning the BGK equation have been obtained during the last decade. For example, Perthame has proved in [START_REF] Perthame | Global existence to the bgk model of boltzmann equation[END_REF] the existence and stability of a distribution solution throughout the space. This result has been extended to a bounded domain with various boundary conditions by Ringeisen [START_REF] Ringeisen | Contributions à l'étude Mathématique des Equations Cinétiques[END_REF]. More recently, Perthame and Pulvirenti have proved the existence and uniqueness of a "mild" solution with weighted estimates in L ∞ [START_REF] Perthame | Weighted L ∞ bounds and uniqueness for the boltzmann bgk model[END_REF]. We also mention the result of Issautier [START_REF] Issautier | Convergence of a weighted particle method for solving the boltzmann (bgk) equation[END_REF] which proved that the "mild" solution of Perthame and Pulvirenti is strong, if certain assumptions of regularity on the initial condition are made. However, it is important to note that in all these results, the authors assume the relaxation time is constant (i.e. τ = 1).

In our study we are interested in the existence and uniqueness of the BGK model with a different source or a different equilibrium density to that defined in the case of fluid dynamics. Namely a density that is suitable for pedestrian movement. In addition, we assume that τ = 1.

It is interesting to study the convergence of such an approximation to the continuous BGK equation. We refer to Mischler's proof for the convergence of a discrete velocity model for the Boltzmann equation [START_REF] Mischler | Convergence of discrete-velocity schemes for the boltzmann equation[END_REF]. There are essentially two distinct points to prove:

• The existence and uniqueness of a discrete velocity model solution.

• Convergence of the discrete kinetic equation towards the continuous equation.

Existence and uniqueness of the model solution.

To define a discrete velocity model of approximation 2.10, we consider the following notations:

We consider V n a grid of N n velocities defined by:

V n = {v n k = k∆v n /k ∈ K n }. where, K n = k = (k 1 , k 2 ) ∈ Z 2 , |k| ≤ B n .
∆v n , B n are two real suites, assumed such that:

∆v n -→ n→+∞ 0 ∆v n B n -→ n→+∞ +∞.
We also define velocity's cells

I n k par I n k = [v n k1 , v n k1 + 1 2 ∆v n [×[v n k2 , v n k2 + 1 2 ∆v n [
. The discrete velocity model approximation 4.1 is then given by the following systems:

(3.1)    ∂f n k (t, x) ∂t + v n k .∇ x f n k (t, x) = f n eq,k (v k , ρ) -f n k (t, x)) D [0, T ] × R 2 x k ∈ K, f n k (t = 0, x) = f n 0,k (x) k ∈ K,
where f n 0,k is an approximation of the initial density f 0 ,

D [0, T ] × R 2 x = ß W : D [0, T ] × R 2 x -→ R, W continuous, linear ™ , (3.2) D [0, T ] × R 2 x = ß f ∈ C ∞ [0, T ] × R 2 x : supp(f), compact ™ . (3.3)
Our goal is to show the existence and uniqueness of the model 3.1. We consider the characteristic curves associated with the problem 3.1, that are given by:

(3.4) γ k (t) = (γ k1 (t), γ k2 (t)) = x + tv n k1 , y + tv n k2 , k = (k 1 , k 2 ) ∈ K.
These curves are solutions of the following equations:

(3.5)    dγ k (t) dt = v n k k ∈ K, γ k (0) = (x, y) T k ∈ K. with v n k = v n k1 , v n k2
. Along these curves, the system solution 3.1 satisfies the following system of ordinary differential equations:

(3.6)    df n k (t, γ k (t)) dt = f n eq,k (v k , ρ) -f n k (t, γ k (t)) k ∈ K, f n k (t = 0, x) = f n 0,k (x) k ∈ K. We pose f n k (t, x) = f n k (t, γ k (t)) f or k ∈ K,
where, f n k is the value of f n along these characteristic curves. We introduce the "mild" form of the system (3.6) obtained by integration along the characteristic curves (3.4), for k ∈ K,

(3.7) f n k (t, x) = f n 0,k + t 0 Ä f n eq,k (v k , ρ) -f n k (s, x) ä ds .
For a given time, we define the following functional space:

L 1 R 2 x = f (t) = (f k ) k∈K , ||f (t)|| 1 = k∈K R 2 x ||f k (t, x)||dx < ∞ .
For a time T > 0, consider the following Banach space:

X T = C [0, T ], L 1 R 2 x , with the following norm: ||f || X T = sup t∈[0,T ] ||f (t)|| 1 .
Our theoretical study is based on the following theorem:

Theorem 1. (Local existence) Let f n 0 = Ä f n 0,k ä k∈K ∈ L ∞ R 2 x ∩ L 1 R 2
x with f n 0 ≥ 0, then there is a time T > 0 and a constant R such as ∀ t < T , the problem 3.1 admits a unique solution

f n = (f n k ) k∈K ∈ C [0, T ], L 1 R 2 x
and which satisfies the following estimates:

(3.8) sup n sup [0,T ] R 2 x k∈K f n k (t, x) dx ≤ Θ (T ) , and
(3.9) ρ(t, x) = k∈K f n k (t, x) ≤ R.
If furthermore k∈K ||f n 0,k || ∞ < 1 then the solution of the model has a physical meaning i.e. :

(3.10) ρ(t, x) ≤ 1.
Since this theorem is independent of n, for simplicity, the exponent n is omitted in this section.

Proof of the theorem

We introduce the following function:

ψ k (t, x) = f k (t, x) exp (λt) f or k ∈ K λ > 0.
Therefore, the system 3.6 equivalent to the following system:

(3.11)      d ψ k (t, x) dt = λ ψ k (t, x) + ψ eq,k (t, x) -ψ k (t, x) k ∈ K, ψ k (t = 0, x) = f 0,k (x) k ∈ K.
For all t ∈ [0, T ], we integrate the equation 3.11, we deduce the following "mild" formulation:

ψ k (t, x) = f 0,k (x) + t 0 Ä λ ψ k (s, x) + ψ eq,k (s, x) -ψ k (s, x) ä ds ∀ k ∈ K.

Consider the following operator

A = Ä Aψ k ä k∈K : Aψ k (t, x) = f 0,k (x) + t 0 Ä λ ψ k (s, x) + ψ eq,k (s, x) -ψ k (s, x) ä ds ∀ k ∈ K.
To show that the system (3.11) has a solution, it is enough to show that the operator A has a unique fixed point in the Banach space X T . Indeed, we introduce the set defined by:

B T,a0,λ,R = " ψ = Ä ψ ä k∈K ∈ X T : ψ k ≥ 0, || " ψ|| X T ≤ a 0 ||f 0 || 1 , k∈K ψ k (t, x) ≤ R exp(λt), t ∈ [0, T ], x ∈ R 2 x . Lemma 1. Let " ψ ∈ B T,a0,λ,R , 140 (1)
There is λ 0 , such that ∀λ ≥ λ 0 , we have :

Ä Aψ ä k∈K ≥ 0.
(

) If k∈K ψ k (t, x) ≤ R exp(λt) then, there are two constant R 0 , T , such that: ∀ R ≥ R 0 and t ∈ [0, T ], k∈K " Aψ k (t, x) ≤ R exp(λt). 2
(

) Let C 2 = 1 2πv 2 m exp -1 2 ||v k -v d || 2 v 2 m -1 such that (C 2 + λ) T = 1 3
, then there is a constant a 0 such that:

(3.12) || Aψ|| X T ≤ a 0 ||f 0 || 1 (4) Let " ψ 1 , " ψ 2 ∈ B T,a0,λ,R
, then there exist C 3 such that:

|| ' Aψ 1 (t) -' Aψ 2 (t)|| 1 ≤ (λT + C 3 T + T ) || " ψ 1 -" ψ 2 || X T
Proof of the lemma:

(1) Since f 0,k (x) ≥ 0 k ∈ K, then Ä Aψ ä k∈K ≥ 0 if, Ä λ ψ k (s, x) + ψ eq,k (s, x) -ψ k (s, x) ä ≥ 0 k ∈ K.
We have

λ ψ k (s, x) + ψ eq,k (s, x) -ψ k (s, x) = λ ψ k (s, x) + k∈K ψ k (s, x) 2πv 2 m exp Å - 1 2 ||v -v d || 2 v 2 m ã -ψ k (s, x) = Å λ + 1 2πv 2 m exp Å - 1 2 ||v k -v d || 2 v 2 m ã -1 ã ψ k (s, x) + l∈K l =k ψ l (s, x) 2πv 2 m exp Å - 1 2 ||v -v d || 2 v 2 m ã .
Since l∈K l =k "

ψ l (s,x) 2πv 2 m exp -1 2 ||v-v d || 2 v 2 m ≥ 0 then λ ψ k (s, x) + ψ eq,k (s, x) -ψ k (s, x) ≥ 0 if: Å λ + 1 2πv 2 m exp Å - 1 2
||v k -v d || 2 v 2 m ã -1 ã ≥ 0, hence, if: λ ≥ λ 0 = 1 - 1 2πv 2 m exp Å - 1 2
||v k -v d || 2 v 2 m ã .
(2) we have k∈K ψ k (t, x) ≤ R exp(λt), hence:

k∈K " Aψ k (t, x) = k∈K f 0,k (x) + t 0 k∈K λ ψ k (s, x) + k∈K ψ eq,k (s, x) - k∈K ψ k (s, x) ds = k∈K f 0,k (x) + t 0 k∈K ψ k (s, x) k∈K 1 2πv 2 m exp Å - 1 2 ||v k -v d || 2 v 2 m ã -1 + k∈K λ ψ k (s, x)ds. ≤ k∈K ||f 0,k || ∞ + C 1 R λ (exp(λt) -1) + R exp(λt) -R.
with

C 1 = k∈K 1 2πv 2 m exp -1 2 ||v k -v d || 2 v 2 m -1 .
145

Assuming that the choice of v k and v d , is made such that

C 1 > 0. Hence k∈K " Aψ k (t, x) ≤ R exp(λt) if, k∈K ||f 0,k || ∞ + C 1 R λ (exp(λt) -1) -R ≤ 0 (3.13) hence R ≥ R 1 = k∈K ||f 0,k || ∞ (3.14) and, t ≤ T = 1 λ ln 1 + λ C 1 R R - k∈K ||f 0,k || ∞ (3.15) (3) Since " ψ ∈ B T,a0,λ,R , we have || " ψ|| X T ≤ a 0 ||f 0 || 1 , moreover, k∈K R 2 x | " Aψ k (t, x) |dx = k∈K R 2 x |f 0,k (x) + t 0 Ä λ ψ k (s, x) + ψ eq,k (s, x) -ψ k (s, x) ä ds|dx (3.16) ≤ ||f 0 || 1 + λT || " ψ|| X T + R 2 x | t 0 k∈K ψ k (s, x) Å 1 2πv 2 m exp Å - 1 2 ||v k -v d || 2 v 2 m ã -1 ã ds|dx (3.17) ≤ ||f 0 || 1 + λT || " ψ|| X T + T C 2 || " ψ|| X T (3.18) ≤ (1 + λa 0 T + T a 0 C 2 ) ||f 0 || 1 , (3.19) with C 2 = 1 2πv 2 m exp -1 2 ||v k -v d || 2 v 2 m -1 . Whence, || Aψ|| X T ≤ a 0 ||f 0 || 1 , if (C 2 + λ) T = 1, i.e. λ = 1
T -C 2 and a 0 , satisfies the following equation: 1 + λa 0 T + T a 0 C 2 = a 0 . Indeed, the constant a 0 is given by: (3.20)

a 0 = 1 1 -(C 2 + λ) T
where , (C 2 + λ) T = 1.

(4)

|| ' Aψ 1 (t) -A " ψ 2 (t)|| 1 = k∈K R 2 x | ' Aψ 1 k (t, x) -' Aψ 2 k (t, x) |dx = k∈K R 2 x | t 0 λ Ä ψ 1 k (s, x) -ψ 2 k (s, x) ä ds + t 0 Ä ψ 1 eq,k (s, x) -ψ 2 eq,k (s, x) ä ds + t 0 Ä ψ 2 k (s, x) -ψ 1 k (s, x) ä ds|dx.
By definition of the equilibrium density, we have,

ψ 2 eq,k (s, x) -ψ 1 eq,k (s, x) = 1 2πv 2 m exp Å - 1 2 ||v k -v d || 2 v 2 m ã l∈K ψ 1 l (s, x) - l∈K ψ 2 l (s, x)
We pose

C 3 = k∈K 1 2πv 2 m exp -1 2 ||v k -v d || 2 v 2 m
Hence by a change of variable, we find:

|| ' Aψ 1 (t) -A " ψ 2 (t)|| 1 ≤ λT || " ψ 1 -" ψ 2 || X T + C 3 T || " ψ 1 -" ψ 2 || X T + T || " ψ 1 -" ψ 2 || X T ≤ (λT + C 3 T + T ) || " ψ 1 -" ψ 2 || X T .
Which ends the proof of the lemma.

According to 1, 2, 3 of the lemma, for any

λ 0 ≤ λ (λ = 1 T -C 2 ), R ≥ R 1 , t ≤
T and a 0 verifies (3.20), we have: if " ψ ∈ B T,a0,λ,R then Aψ ∈ B T,a0,λ,R . According to 4 of the lemma, we have

|| ' Aψ 1 -A " ψ 2 || X T ≤ (λT + C 3 T + T ) || " ψ 1 -" ψ 2 || X T . Moreover λ < 1 T -(C 3 + 1) := λ 1 then: (λT + C 3 T + T ) < 1.
Hence the operator A : B T,a0,λ,R → B T,a0,λ,R is a contraction. Banach's fixed point theorem refers to the local existence of the model solution.

From the foregoing, there exist λ, (λ

0 ≤ λ < λ 1 , λ = 1 T -C 2 )
, T , a 0 and R, such that the problem 3.1 has a unique positive solution

f n = (f n k ) k∈K ∈ C [0, T ], L 1 R 2 x
and which satisfies :

R 2 x k∈K f n k (t, x) dx ≤ a 0 ||f 0 || 1 exp(-λt) ( since || " ψ|| X T ≤ a 0 ||f 0 || 1 ).
Hence,

sup [0,T ] R 2 x k∈K f n k (t, x) dx ≤ a 0 ||f 0 || 1 .
Since a 0 ||f 0 || 1 does not depend on n and a 0 = a 0 (T ), then:

sup n sup [0,T ] R 2 x k∈K f n k (t, x) dx ≤ Θ (T ) .
where, Θ (T ) = a 0 (T )||f 0 || 1 , from where the estimation (3.8).

In addition the solution satisfies

ρ (t, x) ≤ R, t ∈ [0, T ] x ∈ Ω, R ≥ R 1 ,
where

R 1 = k∈K ||f 0,k || ∞ . Moreover if k∈K ||f 0,k || ∞ < 1, (R 1 ≤ 1), then we choose R = R 1 such that ρ (t, x) ≤ 1, ∀ t ∈ [0, T ], ∀ x ∈ Ω,
Hence the estimation (3.10). That ends the proof of the theorem 1.

f n = (f n k ) k∈K ∈ C [0, T max [, L 1 R 2
x , moreover satisfies the following estimate:

(3.21) sup n sup [0,T ] R 2 x k∈K f n k (t, x) dx ≤ Θ (T ) .
In order to prove the convergence of this solution, we define the following functions: Then the discrete model 3.1 can be linked to the BGK equation 2.1 by the following equation:

f n (t, x, v) = k∈K f n k (t, x)1 n k (v), (3.22) f n eq (t, x, v) = k∈K f n eq,k (t, x)1 n k (v), (3.23) C n (v) = k∈K v n k 1 n k (v) (3.24) f 0 (0, x, v) = k∈K f 0,n k (x)1 n k (v) (3.
(3.26)        ∂f n (t, x, v) ∂t + C n (v).∇ x f n (t, x, v) = f n eq (v) -f n (t, x, v) D ]0, T max [×R 2 x × R 2 v , f 0 (0, x, v) = k∈K f 0,n k (x)1 n k (v),
where,

D ]0, T max [×R 2 x × R 2 v = ß W : D ]0, T max [×R 2 x × R 2 v -→ R, W continuous, linear ™ , (3.27) D ]0, T max [×R 2 x × R 2 v = ß f ∈ C ∞ ]0, T max [×R 2 x × R 2 v : supp(f), compact ™ . (3.28)
We now denote our convergence result in the following theorem:

Theorem 2. let v n such that ∆v n -→ n→+∞ 0 ∆v n B n -→ n→+∞ +∞, and
C n : R 2 v → R 2 v v → C n (v)
such that:

(3.29) C n is locally, uniformly bounded in L ∞ loc R 2 v 2 , (3.30) C n (v) -→ n→+∞ v simply.
Then we can extract a sub-sequence noted (f n ) n which converges weakly in

L 1 ]0, T max [×R 2 x × R 2 v , ∀T max ≥ 0 to a solution of the equation (2.

1).

Proof

We refer to the works of Perthame [START_REF] Perthame | Global existence to the bgk model of boltzmann equation[END_REF] and [START_REF] Mieussens | Convergence of a discrete-velocity model for the boltzmann-bgk equation[END_REF], we divide the proof into 4 steps: step 1: weak convergence of f n According to (3.21), (3.22), it is clear that f n satisfies the following estimate:

(3.31) sup n sup [0,T ] R 2 x ×R 2 v f n (t, x, v) dxdv ≤ Θ 1 (T ) ,
which implies that the following (f n ) n is equi-integrable.

According to Dunford-Pettis theorem [START_REF] Dunford | Linear operators part I: general theory[END_REF], we can extract a subsequence noted (f n ) n which converges weakly in

L 1 ]0, T max [×R 2 x × R 2 v , towards f ∈ L 1 ]0, T max [×R 2 x × R 2 v , i.e.: (3.32) f n n→+∞ f in L 1 ]0, T max [×R 2 x × R 2 v . Hence f n -→ n→+∞ f in D ]0, T max [×R 2 x × R 2 v
, for all T max ≥ 0. We thus obtain the convergence of the transport term of (3.26)

towards ∂f + v.∇f in D ]0, T max [×R 2
x × R 2 v ., i.e.:

(3.33)

∂f n + C n (v).∇ x f n -→ n→+∞ ∂f + v.∇ x f in D ]0, T max [×R 2 x × R 2 v . Indeed, let ϕ ∈ D ]0, T max [×R 2 x × R 2 v , ]0,Tmax[×R 2 x ×R 2 v (∂ t f n + C n (v).∇ x f n ) ϕdtdxdv = - ]0,Tmax[×R 2 x ×R 2 v f n (∂ t ϕ + C n (v).∇ x ϕ) dtdxdv = - ]0,Tmax[×R 2 x ×R 2 v f n (∂ t ϕ) dtdxdv - ]0,Tmax[×R 2 x ×R 2 v f n (C n (v).∇ x ϕ) dtdxdv.
We have :

]0,Tmax[×R 2 x ×R 2 v f n (C n (v).∇ x ϕ) dtdxdv = ]0,Tmax[×R 2 x ×R 2 v f n (C n (v) -v) .∇ x ϕdtdxdv + ]0,Tmax[×R 2 x ×R 2 v f n v.∇ x ϕdtdxdv
According to (3.29), (3.30)

]0,Tmax[×R 2 x ×R 2 v f n (C n (v) -v) .∇ x ϕdtdxdv ≤ ||f n || L 1 (]0,Tmax[×R 2 x ×R 2 v ) ||C n (v) -v|| L ∞ (R 2 v ) ||∇ x ϕ|| L ∞ (]0,Tmax[×R 2 x ×R 2 v ) . From (3.30), ||C n (v) -v|| L ∞ (R 2 v ) -→ n→+∞ 0, consequently, (3.34) ]0,Tmax[×R 2 x ×R 2 v f n (C n (v) -v) .∇ x ϕdtdxdv -→ n→+∞ 0. Moreover, since v.∇ x ϕ ∈ L ∞ ]0, T max [×R 2 x × R 2 v , then (3.35) ]0,Tmax[×R 2 x ×R 2 v f n v.∇ x ϕdtdxdv -→ n→+∞ ]0,Tmax[×R 2 x ×R 2 v f v.∇ x ϕdt.
In addition,

∂ t ϕ ∈ L ∞ ]0, T max [×R 2 x × R 2 v , then 185 (3.36) ]0,Tmax[×R 2 x ×R 2 v f n (∂ t ϕ) dtdxdv -→ n→+∞ ]0,Tmax[×R 2 x ×R 2 v f (∂ t ϕ) dtdxdv.
From, (3.34), (3.35) and (3.36), we have

]0,Tmax[×R 2 x ×R 2 v f n (∂ t ϕ + C n (v).∇ x ϕ) dtdxdv -→ n→+∞ ]0,Tmax[×R 2 x ×R 2 v f (∂ t ϕ + v.∇ x ϕ) dtdxdv,
hence the result (3.33).

For the convergence of the non-linear part, we first have the convergence of ρ n , according to (3.32):

R 2 v f n (t, x, v) dv -→ n→+∞ R 2 v f (t, x, v) dv weakly in L 1 ]0, T max [×R 2 x ,
i.e.

(3.37)

ρ n (t, x) n→+∞ ρ(t, x) in L 1 ]0, T max [×R 2 x .
Step 2 : weak convergence of f n eq

We need the following lemma:

Lemma 2. Suppose that f n = (f n k ) k∈K satisfies the inequality (3.21), then ∀ T ≥ 0, there is C (T ) such that 190 (3.38) sup n sup [0,T ] R 2 x ×R 2 v f n eq (t, x, v) dxdv ≤ C (T ) ,

Proof of the lemma 2

We have

f n eq (t, x, v) = k∈K f n eq,k (t, x)1 n k (v), with, f n eq,k (t, x) = 1 2πv 2 m exp -1 2 ||v k -v d || 2 v 2 m , ρ (t, x) = k∈K f k (x, t) . R 2 x ×R 2 v f n eq (t, x, v) dxdv = R 2 x ×R 2 v l∈K f l (x, t) k∈K 1 2πv 2 m exp Å - 1 2 ||v k -v d || 2 v 2 m ã 1 n k (v)dxdv = R 2 x R 2 v l∈K f l (x, t) dx R 2 v k∈K 1 2πv 2 m exp Å - 1 2 ||v k -v d || 2 v 2 m ã 1 n k (v)dv = R 2 x l∈K f l (x, t) dx I n k k∈K 1 2πv 2 m exp Å - 1 2 ||v k -v d || 2 v 2 m ã dv
from where,

sup n sup [0,T ] R 2 x ×R 2 v f n eq (t, x, v) dxdv ≤ sup n sup [0,T ] R 2 x l∈K f l (x, t) dx.C (I n k ) according to (3.21) ≤ C (T ) ,
where C (T ) := Θ(T ).C (I n k ). According to this lemma (f n eq ) is weakly compact in

L 1 ]0, T max [×R 2 x × R 2 v , so we can extract a subsequence noted (f n eq ) such that f n eq n→+∞ g in L 1 ]0, T max [×R 2 x × R 2 v
, where from the step (1) and (2), the low limit of f n satisfies the equation

(3.39) ∂f + v.∇f = g -f in D ]0, T max [×R 2 x × R 2 v .
The next step is to show that g = f eq

Step 3: strong convergence of ρ n .

According to the compactness lemma on averages obtained by Mischler [START_REF] Mischler | Convergence of discrete-velocity schemes for the boltzmann equation[END_REF], and [?], on sets of bounded velocities f n is strongly compact therefore. Indeed, we have for any R C (extracting again subsequences):

(3.40)

|v|≤R C f n (t, x, v) dv -→ n→+∞ |v|≤R C f (t, x, v) dv, strongly in L 1 ]0, T max [×R 2 x .
From the uniform estimates (3.31), we thus obtain

(3.41) f n -→ n→+∞ f in L 1 ]0, T max [×K x × R 2 v . ∀ K x compact of R 2 x .
Hence from the above and the equation (3.37)

(3.42) ρ n (t, x) = R 2 v f n (t, x, v) dv -→ n→+∞ ρ(t, x) = R 2 v f (t, x, v) dv L 1 (]0, T max [×K x ) . ∀K x compact.
Step 4: Passing to the limit According to the step (3),

(3.43) ρ n (t, x) -→ n→+∞ ρ(t, x) L 1 (]0, T max [×K x ) ∀K x compact,
then, we can extract a sub-sequence again, as

(3.44) ρ n (t, x) -→ n→+∞ ρ(t, x) a.e [0, T max [×R 2 x .
On the other hand, by hypothesis we have:

(3.45) ||v n k -v d || -→ n→+∞ ||v -v d || ∀k ∈ K.
Hence from (3.44), (3.45),

(3.46) f n eq -→ n→+∞ f eq a.e. [0, T max [×R 2 x × R 2 v .
Combine this results with those of the step (2) we have:

(3.47) g = f eq .
As a result, the left side of the equation (3.26) converges to (f eq -f ) weakly in

L 1 ]0, T max [×R 2 x × R 2 v .
We conclude that f is solution of the equation BGK (2.1).

Results and numerical simulations

4.1. Numerical method. The general idea of the semi-Lagrangian method used is to fix a grid in the velocity space and to transform the kinetic equation into a set of linear hyperbolic equations with source terms. We refer to [START_REF] Dimarco | Towards an ultra efficient kinetic scheme. part i: Basics on the bgk equation[END_REF][START_REF] Dimarco | Self-alignment driven by jump processes: Macroscopic limit and numerical investigation[END_REF] for the detailed description of this numerical method. Here we recall only the basic principles. We summarize the semi-Lagrangian numerical method used in this work as follows:

(1) the discretization of the BGK model equation in the velocity space.

(2) A splitting procedure of time between transport and relaxation operators for each of the system evolution equations (4.1). ( 3) The exact resolution of the transport part which means without using a spatial mesh, the initial data of this step are given by the solution of the relaxation operator.

(4) The resolution of the relaxation part on the grid with initial data defined by the value of the distribution function at the center of the cells after the transport step.

We introduce a Cartesian grid V in R 2 in two-dimensional velocity space, and a set K multi-indices of Z 2 , such that:

K = {(-1, 0); (-1 , 1); (0, 1); (1, 1); (1
   ∂f k (t, x) ∂t + v k .∇ x f k (t, x) = 1 τ (f eq,k (v k ) -f k (t, x)) k ∈ K, f k (t = 0, x) = f 0,k (x) k ∈ K.
We describe the first step of the [t 0 → t 1 ] method, starting from t 0 = 0 and then we generalize at an arbitrary time step.

With splitting, the first step is reduced to the N v linear transport equation resolution of the form

(4.2) ∂f k (t, x) ∂t + v k .∇ x f k (t, x) = 0 k ∈ K.
In order to solve this part, we consider for each equation of the system (4.2), the initial data defined by:

(4.3) Û f k (x, t = 0) = f 0,k (x) x ∈ Ω, k ∈ K.
Thanks to this reconstruction, the exact system solution (4.2) at time t 1 = t 0 + ∆t = ∆t is given by:

(4.4) Û f * k (x) = Û f k (x -v k ∆t) x ∈ Ω, k ∈ K.
To complete a step of time, we must calculate the solution of the interaction part of the equation (4.1) on the points of the grid, (4.5)

∂ t f k = 1 τ (f eq,k -f k ) k ∈ K,
where the initial data is given by the resolution of the transport step to time

t 1 = t 0 + ∆t, Ä Û f * k (x) ä k∈K .
To solve (4.5), we define the value of the equilibrium distribution at the instant

t 1 , Ä f 1 eq,k ä k∈K , f 1 eq,k (v k ) = ρ(t, x) 2πv 2 r,m,s exp Ç - 1 2 ||v k -v d || 2 v 2 r,m,s å , k ∈ K, with ρ(t, x) = k∈K Û f * k (x).
Finally, the solution of relaxation equation (4.5) is given by: (4.6)

f 1 k = e (-∆t τ ) Û f * k + 1 -e (-∆t τ ) f 1 eq,k , k ∈ K.
Given a density value

Ä Û f n k (x)
ä k∈K at the moment t n for x ∈ Ω, k ∈ K, then the density at instant t n+1 can be calculated as follows:

• Û f * k (x) = Û f n k (x -v k ∆t). • f n+1 k = e (-∆t τ ) Û f * k + 1 -e (-∆t τ ) f n+1 eq,k .
The local density at the instant t n+1 defined by:

ρ(t n+1 , x) = k∈K f n+1 k (x).
By Referring to [START_REF] Dimarco | Towards an ultra efficient kinetic scheme. part i: Basics on the bgk equation[END_REF][START_REF] Stracquadanio | High order semi-Lagrangian methods for BGK-type models in the kinetic theory of rarefied gases[END_REF][START_REF] Michael Lentine | An unconditionally stable fully conservative semi-lagrangian method[END_REF], ?], this scheme is unconditionally stable , however for reasons of precision, the time step is chosen to satisfy the condition: ∆t ∆x < 1 because the maximum speed of the pedestrians is fixed at one. In conclusion, we summarize the procedures where the steps of this scheme in the following algorithm:

Require: f 0 i (x) i∈K : initial data. for m = 0 : N t -1 do
• Resolution of the transport part the N v equations of the system 4.2 with (f m i ) i∈K the initial data, we get

Ä Û f m i ä i∈K .
• Computation of equilibrium density Ä f m eq,i ä i∈K .

• Relaxation term resolution (4.5) with

Ä Û f m i ä i∈K initial data. end for • (f m i (x)) i∈K .
The output parameters of the algorithm are the numerical solution of the original problem defined by the two equations (4.1) .

4.2.

Convergence test for the semi-Lagrangian scheme. Our goal in this paper is to validate the proposed mathematical model for pedestrian motion. This choice is motivated by considering a small number of pedestrians (N ≤ 100 pedestrians).

In order to analyze the convergence of the spatial discretization of the semi-Lagrangian method. We solve the kinetic equation (4.1), with initial data that is a piecewise constant function in two-dimensional space. They are concentrated densities in a center disk x 0 with radius R d :

(4.7) ∀k ∈ K, f k (t = 0, x) = ® f 0 k (x) , ||x -x 0 || ≤ R d , 0
This initial data corresponds to a small number of pedestrians N = 100.

We performed the simulation with a temporal discretization step, ∆t = 2.10 -2 . In order to estimate the accuracy of the method, we used as a reference a solution f * computed with a space step ∆x = ∆y = 1 2 6 = 0.0156, then we estimate the solution with different steps ∆x = 1 2 n-1 , n = 1, ..., 5, with ∆x = ∆y. To evaluate the convergence, we calculated the difference between the norm L 1 of the reference solution and the density ρ(x) estimated at time t = 2.5 s with different steps ∆x and ∆y, (4.8)

Error(∆x) =||ρ(∆x, t = 2.5) -ρ * (∆x, t = 2.5)|| L 1 , with ρ * (., t = 2.5) = k∈K f k, * (., t = 2.5).

We illustrate the results obtained in Fig. 3. . The Fig. 3, shows that the norm of error L 1 decreases linearly with respect ∆x. In order to demonstrate the convergence order, we compute the so called experimental order of convergence (EOC), (4.9)

EOC := log 2 ||ρ(∆x, t = 2.5) -ρ * || L 1 ||ρ( ∆x 2 , t = 2.5) -ρ * || L 1 ,
for considered spacial steps ∆x = 1 2 n-1 , n = 1, ..., 5, the obtained results are given in the following table (see Table 265 1 )

Table 1. Experimental order of convergence (EOC) for the semi-Lagrangian scheme.

n the order of refinement ∆x = We observe that pedestrian density diffuses into space and pedestrians point to the following directions:

4.3.

Movement of a group of pedestrians towards a single desired direction v d . In this paragraph our aim is to show the adaptation of all pedestrians to a desired configuration namely the movement in a desired direction v d . We consider the BGK model (4.1) with the same initial condition (4.7) corresponding to 100 pedestrians. The 270 equilibrium density modeling the trend towards the desired configuration is defined as follows:

(4.10)

     f eq,k (v k ) = ρ(t, x) 2πv 2 r,m,s exp Ç - 1 2 ||v k -v d || 2 v 2 r,m,s å , k ∈ K v d = [1, 1] T , K = {(-1, 0); (1, -1); (1, 0); (1, 1); (0, 1)} v k ∈ V = {v k = k∆v /k ∈ K, ∆v = 1} .
We assume that all pedestrians have the same speed v = 1.00 m /s. We considered the same time step as in the first case ∆t = 2.10 -2 ∆x = ∆y = 5.10 -2 .

We then represent the evolution of pedestrian density for 6 different instants and for two values of relaxation time τ = 5.10 -3 (see Fig. 5) and τ = 5.10 -2 (see Fig. The figures show that the crowd is divided into two groups in order to reach their state of equilibrium. It consists the movement of the two groups towards the two desired directions v d,1 , v d,2 .

Conclusion and looking forward

In this paper we have developed the kinetic approach for the dynamics of a crowd, based on the BGK equation.

The existence and uniqueness of the proposed discrete velocity model solution has been demonstrated thanks to the Banach fixed point theorem. Thus, the convergence of this discrete model towards the continuous BGK model is proven. Numerical simulations using the semi-Lagrangian method are performed. The mathematical model proposed is capable of describing the tendency of a crowd towards a situation of equilibrium, namely the tendency towards a desired direction.

As already laid out in the introduction, it should be clear that the aim of this paper is only to adapt the BGK model in kinetic theory to the movement of a crowd. In other words, the derivation of the equilibrium function f eq describing a pedestrian tendency to achieve a desired direction. The most important perspective, is related to the derivation of a equilibrium function describing the panic conditions, evacuation problems or lane formation. The relaxation time τ is assumed constant this is a strong assumption considering the application. This coefficient τ may depend on f eq or on the density ρ.

2. 1 . 2 D

 12 Boltzmann equation: the BGK model for a crowd. We consider a system composed of N particles (pedestrians) randomly distributed in a two-dimensional bounded domain Ω ⊂ R 2 . At the moment t = t 0 , the pedestrians are distributed in a disk D 0 of radius r D and center M 0 (x 0 , y 0 ). The initial overall density is then ρ 0 = N πr

Figure 1 .

 1 Figure 1. Density of pedestrians characterized by (a) normal state: pedestrians have the ability to direct towards n possible directions. (b) the equilibrium state: all pedestrians look for a certain destination comfortable noted v d to achieve the target located at the point xc .

Figure 2 .

 2 Figure 2. The density function of the first component vx of the speed (a) for 693 students walking outside the library at the University of Sydney. Curve a represents the measured distribution and curve b represents the Maxwell-Boltzmann distribution, vx = 1.53m.s -1 ,vr,m,s = 0.201m.s -1 , (b) for 628 pedestrians on a pedestrian crossing in Sydney, vx = 1.44m.s -1 ,vr,m,s = 0.228m.s -1

Remark 2 . 2 x 3 . 2 .

 2232 In our proof we have demonstrated the local existence of problem 3.1, in [0,T]. Through an iteration process. We can successively solve the equation 3.1 with initial conditions in t 0 = 0 until T , t 1 = T until t 2 = 2T , ... By concatenation we build a maximum solution on [0, T max [, with T max = sup j t j . This solution belongs C [0, T max [, L 1 R Convergence of the discrete kinetic equation towards the continuous one. According to the previews section, we have shown the existence and uniqueness of the discrete model solution(3.1)

25 ) with 1 n

 251 k the indicator function on the velocity cells I n k .

  , 0); (1, -1); (0, -1); (-1, -1)} which means we discretize the square [-1, 1] × [-1, 1]. In what follows, all the simulations are done in a square space domain Ω = [0, 20] × [0, 20]. Thanks to the discrete velocity approximation above, the continuous distribution function f is replaced by 8 vector where each component is supposed to be an approximation of the distribution function f , i.e. f k (t, x) ≈ f (t, x, v k ), and the original kinetic equation (1.1) is replaced by a set of 8 evolution equations for f k of the following form:(4.1)

Figure 3 .

 3 Figure 3. Precision of the schema in space with the error given by (4.8).

1 ,

 1 -1), (-1, 0); (-1, 1); (1, -1); (1, 0); (1, 1); (0, 1); (0, -1) ™ , V = {v k = k∆v /k ∈ K, ∆v = 1} , we represented the solution at 6 different moments, with the same initial condition defined by (4.7), and space step ∆x = ∆y = 5.10 -2 and ∆t = 2.10 -2 , the results obtained are represented in Fig.4 a-f. (a) t = 0s (b) t = 2s (c) t = 2.44s (d) t = 3.16s (e) t = 4.2s (f) t = 6s.

Figure 4 .

 4 Figure 4. The evolution of the density ρ(x, y) during the instants, (a): t = 0s, (b): t = 2s, (c): t = 2.44s (d): t = 3.16s, (e): t = 4.2s, (f): t = 6s.

  ß (±1, 0); (1, ±1); (-1, ±1); (0, ±1) ™ .

  6.) (a) t = 0s (b) t = 1.04s (c) t = 1.76s (d) t = 2.52s (e) t = 3.16s (f) t = 6s.

Figure 5 .

 5 Figure 5. Evolution of the local density during the times : (a): t = 0s, (b): t = 1.04s, (c): t = 1.76s (d): t = 2.52.16s, (e): t = 3.16s, (f): t = 6s, with a relaxation time τ = 5.10 -3 .

Figure 6 . 4 . 4 .

 644 Figure 6. Evolution of the local density during the times : (a): t = 0s, (b): t = 1.04s, (c): t = 1.76s (d): t = 2.52.16s, (e): t = 3.16s, (f): t = 6s, with a relaxation time τ = 5.10 -2 .

Figure 7 .

 7 Figure 7. Evolution of local density during the times : (a): t = 0s, (b): t = 1.12s, (c): t = 2.08s (d): t = 3.2s, (e): t = 4.24s, (f): t = 6s, with a relaxation time τ = 5.10 -3 .

  (a) t = 0s (b) t = 1.12s (c) t = 2.08s (d) t = 3.2s (e) t = 4.24s (f) t = 6s.

Figure 8 .

 8 Figure 8. Evolution of local density during the times : (a): t = 0s, (b): t = 1.12s, (c): t = 2.08s (d): t = 3.2s, (e): t = 4.24s, (f): t = 6s, with a relaxation time τ = 5.10 -2 .

  According to the results in Table1, the experimental order of convergence EOC ≈ 1 showing that the method is first order accurate in space. Then, to show the possible directions for pedestrians

		1 2 n-1	Error(∆x)	EOC
	1	1	0.6138	--
	2	0.5000	0.2964	1.0502
	3	0.2500	0.1450	1.0315
	4	0.1250	0.0715	1.0200
	5	0.0625	0.0362	0.9820

Acknowledgments

The authors are grateful to the reviewers for their insightful remarks and corrections. Their feedback has greatly influenced the quality of this paper

France, e-mail: mohammed.elrhabi@unilasalle.fr; Abdelilah Hakim, LAMAI, FST Marrakech, Université Cadi Ayyad, Morocco, e-mail: a.hakim@uca.ma


Samir Jaber

Giuseppe Citerio

Acute respiratory failure and mechanical ventilation in the context of the COVID-19 pandemic: why a special issue in ICM?