A King

W Ludwig

M Herbig

J.-Y Buffière

A A Khan

N Stevens

T J Marrow

Three dimensional in-situ observations of short fatigue crack growth in magnesium

Keywords: fatigue, x-ray diffraction (XRD), magnesium alloys, grain boundary structure

Three-dimensional (3D) short fatigue crack growth behaviour in a cast magnesium alloy Elektron 21 was studied using a combination of X-ray diffraction contrast tomography (DCT) and microtomography at the European Synchrotron Radiation Facility (ESRF). A 3D map of grain shapes and crystallographic orientations was produced in a miniature fatigue specimen using DCT. A focused ion beam (FIB) instrument was used to introduce small notches in selected grains. Synchrotron microtomography was used to study the evolution of fatigue cracks from these notches during interrupted in-situ fatigue cycling. Stage I crack growth occurred preferentially on the basal plane of the Magnesium hexagonal close packed crystal, with local crack growth rates between 4 and 40 nm/cycle. Retardations in the local crack growth rate were observed in certain grains, and at certain grain boundaries. The observed interactions between the crack and the polycrystalline microstructure can be explained using Schmid factors and a development of the tilt-twist model of Zhai and Wilkinson [1].

Introduction

Cracks are three-dimensional, and they interact with the microstructure of materials, which is also three-dimensional. During constant stress amplitude loading, crack growth rates generally increase exponentially with the crack length since the mechanical crack driving force increases with crack length, unless significant relaxation occurs in response to the effect of the crack (i.e. for displacement controlled tests), or in the case of complex residual stress gradients for example. Hence the majority of the lifetime of a crack in a structural component may be spent when the crack is quite short, and uncertainties in short crack propagation behaviour can have a strong influence on lifetime prediction. There is therefore a need to good understanding of the behaviour of short cracks, which may be short relative to the microstructure length scale (i.e. grain size, for example, in which case interactions between cracks and microstructure are important), or the component scale, in which case there are issues of elastic plastic behaviour and constraint to be addressed [START_REF] Suresh | Fatigue of Materials[END_REF][START_REF] Krupp | Fatigue crack propagation in metals and alloys: microstructural aspects and modelling concepts[END_REF][START_REF] Davidson | Small fatigue cracks[END_REF][START_REF] Schijve | Fatigue of structures and materials[END_REF]. Such cracks are in the regime within which the experimental study and modelling of crack behaviour are most difficult, and in which detailed understanding and data are least sound. In the last twenty years, significant experimental and theoretical efforts have been devoted to the characterisation and modelling of short fatigue crack propagation [e.g. [START_REF] Newman | The merging of fatigue and fracture mechanics concepts : a historical perspective[END_REF][START_REF] Mcdowell | Microstructure-sensitive computational modelling of fatigue crack formation[END_REF]. However, crack growth rate predictions for such defects remain far from completely satisfactory. The consequent uncertainty in the prediction of crack behaviour and lifetime can lead to conservatism in engineering design and may result in unnecessary and uneconomic repair of structures in which the early stages of cracking are identified in service.

Studies of short crack behaviour are generally restricted to the crack length at the surface, which provides little or no information about the development of the crack shape in the bulk. Some understanding of crack shapes can be obtained by methods such as serial metallographic sectioning, which can be used to reveal the crack shape at the end of a fatigue test [START_REF] Clement | Short crack behavior in nodular cast-iron[END_REF][9][10]. Similarly, interrupted testing and fracture can determine the typical development of crack shape [START_REF] Beretta | Propagation of short cracks in a high strength steel for railway axles[END_REF]. ""Beach-marking"", i.e. defining the position of the crack tip by varying the cyclic stress amplitude, has also been used to record the development of individual cracks [START_REF] Nadot | Fatigue life assessment of nodular cast iron containing casting defects[END_REF]. These observations have shown that physically short cracks (i.e. those which have a size greater than the microstructure scale) may have an approximately semi-circular shape. In some cases, particularly when microstructure has a strong influence on crack development, a simple relationship between the crack length as measured on the surface and length in depth cannot necessarily be assumed [START_REF] Marrow | X-ray tomography of short fatigue cracks in ductile iron[END_REF]. Indeed, very little is found by such techniques concerning the crack shape at the very early stage of growth when its size is comparable to that of the surrounding microstructure.

Computed X-ray tomography [START_REF] Stock | X-ray microtomography of materials[END_REF][START_REF] Maire | On the application of X-ray microtomography in the field of materials science[END_REF] is a relatively new technique that can be used to visualise the internal structure of materials. In its classical form, this non-destructive technique uses two-dimensional (2D) radiographic projections of a sample to reconstruct a three-dimensional (3D) map of the X-ray attenuation coefficient of the material. In the last ten years, significant progress has been made in terms of resolution with both the availability of new third generation synchrotron X-ray sources and new detectors [START_REF] Maire | On the application of X-ray microtomography in the field of materials science[END_REF]. A spatial resolution (i.e. pixel size) better than 1 µm can now be achieved in 3D, which is well suited to the investigation of fatigue cracks [START_REF] Ludwig | Study of the interaction of a short fatigue crack with grain boundaries in a cast Al alloy using X-ray microtomography[END_REF][START_REF] Marrow | High resolution X-ray tomography of short fatigue crack nucleation in austempered ductile cast iron[END_REF][START_REF] Buffière | Three dimensional visualisation of fatigue cracks in metals using high resolution synchrotron X-ray microtomography[END_REF][START_REF] Herbig | 3-D growth of a short fatigue crack within a polycrystalline microstructure studied using combined diffraction and phase-contrast X-ray tomography[END_REF]. In particular, the partially coherent radiation from a third generation synchrotron provides phase contrast enhancement under appropriate conditions, allowing features smaller than the pixel size to be detected [START_REF] Ludwig | Study of the interaction of a short fatigue crack with grain boundaries in a cast Al alloy using X-ray microtomography[END_REF].

X-ray diffraction contrast tomography (DCT) is an emerging technique for mapping grain shape and crystal orientation in polycrystalline materials, non-destructively and in 3D [START_REF] Ludwig | X-ray diffraction contrast tomography: a novel technique for three-dimensional grain mapping of polycrystals. I. Direct beam case[END_REF][START_REF] Johnson | X-ray diffraction contrast tomography: a novel technique for three-dimensional grain mapping of polycrystals. II. The combined case[END_REF][START_REF] Ludwig | Three dimensional grain mapping by x-ray diffraction contrast tomography and the use of Friedel pairs in diffraction data analysis[END_REF]. It is a natural complement to synchrotron microtomography for the investigation of fatigue crack growth, because it allows the grain microstructure to be characterised prior to fatigue loading.

In contrast, high-resolution post-test techniques such as serial sectioning that use focussed-ion-beams to reconstruct the microstructure have limited application to large volumes and also can investigate only the damaged and deformed microstructure after failure, i.e. the damage chronology is lost. By combining computed tomography with diffraction contrast tomography of the undamaged microstructure the interaction of the growing crack with the grains in the bulk of the material can be studied. This enables in-depth studies of early crack growth, including observations of the interaction of cracks with the grains and grain boundaries of polycrystals [START_REF] Herbig | 3-D growth of a short fatigue crack within a polycrystalline microstructure studied using combined diffraction and phase-contrast X-ray tomography[END_REF][START_REF] King | Observations of intergranular stress corrosion cracking in a grain-mapped polycrystal[END_REF]. In particular, Herbig et al studied this combination of techniques to study the propagation of a fatigue crack in a body centred cubic titanium sample containing around 1000 grains [START_REF] Herbig | 3-D growth of a short fatigue crack within a polycrystalline microstructure studied using combined diffraction and phase-contrast X-ray tomography[END_REF]. The local crack propagation rate was measured, and the crystallographic character of the crack plane was described in terms of combinations of expected slip planes. However, interpretation was made difficult by the fact that several slip systems, with high multiplicity, are active in the material studied [START_REF] Herbig | 3-D growth of a short fatigue crack within a polycrystalline microstructure studied using combined diffraction and phase-contrast X-ray tomography[END_REF]. To favour a more quantitative interpretation, in this work a magnesium alloy was used.

In hexagonal close packed magnesium the limited number of available slip systems restricts the number of crystallographic planes that are available for easy crack growth. The primary slip system is <10-10>(0002) (the close packed directions on the hexagonal basal planes). This limits easy plasticity, and crack growth rates are affected by the crystal orientation relative to crack plane [24]. Plasticity in magnesium alloys at room temperature is dominated by this slip system and {10-12} twinning [25]. Twinning can allow the reorientation of grains or parts of grains during loading. Twinning has been observed to interact with fatigue crack initiation and propagation in single crystals, or in rolled or extruded microstructures with strong texture which inhibit basal slip [24,26,27]. The cast material used in the present work has a much weaker texture, and is therefore contains a higher proportion of grains in which basal slip is available for deformation. Primatic slip <10-10>{10-12} has also been observed in fatigue specimens at higher stresses, but has a greater resistance to dislocation movement [27]. This paper reports observations, obtained by DCT and in-situ high-resolution X-ray tomography of miniature fatigue specimens, of the early stages of nucleation and growth of short fatigue cracks in magnesium alloy Elektron21, nucleated within selected grains at fine notches machined by a focused ion beam (FIB). The principal objective was to study the interaction of stage I crystallographic cracks with grain boundaries.

Experimental Details

The magnesium alloy Elektron 21 used in this work was provided by Magnesium Elektron Ltd. as cast plates, measuring 200 mm × 200 mm × 25 mm. They were in the T6 condition (8 hours at 520°C, hot water quenched, 16 hours at 200°C). The nominal chemical composition is given in Table 1. The tensile properties of the asreceived microstructure given in Table 2 are the average of three tensile specimens tested. (60 mm long, 10 mm wide and 2 mm thick with a short gauge section ~4 mm long and ~ 2 mm wide) [START_REF] Khan | In-situ Observation of Damage Mechanisms by Digital Image Correlation during Tension and Low Cycle Fatigue of Magnesium alloys[END_REF]. Miniature fatigue specimens were designed for the experiment, as shown in Figure 1.

The sample length is 14 mm, and the gauge diameter is 1 mm. Their size and microstructure were selected to suit the two imaging techniques, which have certain requirements regarding sample size and microstructure. It has been previously shown that fatigue cracks and stress corrosion cracks can be successfully imaged, with the assistance of phase contrast, using reconstructed volumes that have a voxel size of 0.7 µm (a voxel is the 3D analogue of a pixel) [START_REF] Ludwig | Study of the interaction of a short fatigue crack with grain boundaries in a cast Al alloy using X-ray microtomography[END_REF][START_REF] Marrow | High resolution X-ray tomography of short fatigue crack nucleation in austempered ductile cast iron[END_REF][START_REF] Buffière | Three dimensional visualisation of fatigue cracks in metals using high resolution synchrotron X-ray microtomography[END_REF][START_REF] Herbig | 3-D growth of a short fatigue crack within a polycrystalline microstructure studied using combined diffraction and phase-contrast X-ray tomography[END_REF][START_REF] King | Observations of intergranular stress corrosion cracking in a grain-mapped polycrystal[END_REF]. Using a 2048 2 pixel detector, this implies that the sample diameter should be less than 1.4 mm. The diameter should also be small enough to allow at least 10% transmission of the incoming X-ray beam through the thickest part of the sample. DCT is capable of mapping a few thousand grains in a single scan in suitable microstructures [START_REF] Johnson | X-ray diffraction contrast tomography: a novel technique for three-dimensional grain mapping of polycrystals. II. The combined case[END_REF][START_REF] Ludwig | Three dimensional grain mapping by x-ray diffraction contrast tomography and the use of Friedel pairs in diffraction data analysis[END_REF]. Thus a further requirement is that the grain size be large enough so as not to exceed this number of grains in the volume of interest. For this study a sample gauge diameter of 1 mm was used, and the T6 material was subsequently heat treated at 520°C to increase the average grain size to around 55 µm, to give around 2500 grains in the mapped volume. The studied sample is weakly textured, with a tendency for the {0002} poles to be aligned with the tensile axis of about two multiples of a random texture, based on the grain orientations determined by DCT. Several small nickel spheres (diameter 150 µm, just visible in figure 1) were fixed to the sample in the vicinity of the gauge volume. These served as fiducial reference markers to aid the alignment of different datasets, and when orienting the sample in the FIB to introduce the notches in the desired region.

The imaging experiments were performed at the European Synchrotron Radiation Facility (ESRF) in Grenoble on the high-resolution tomography beam line (ID19) on separate experimental periods within beam time allocation MA-423. For both techniques, images were acquired using a FReLoN camera coupled to a scintillator screen with magnifying optics [START_REF] Labiche | The fast readout low noise camera as a versatile x-ray detector for time resolved extended x-ray absorption fine structure and diffraction studies of dynamic problems in materials science, chemistry, and catalysis[END_REF].

For DCT, it is necessary to capture both diffraction spots and the absorption contrast image of the sample on the detector. A pixel size of 1.75 µm was used, giving a 3.5 mm field of view. A sample to detector distance of 5.1 mm gave sufficient solid angle coverage to acquire diffraction spots from the first five {hk.l} reflections. A beam of 20.5 keV with energy bandwidth E/E=10 -4 was used. A set of 7200 diffraction contrast images was recorded during a full 360° rotation of the sample.

The procedure for reconstruction of the grain shapes and their crystal orientations is described elsewhere [START_REF] Johnson | X-ray diffraction contrast tomography: a novel technique for three-dimensional grain mapping of polycrystals. II. The combined case[END_REF][START_REF] Ludwig | Three dimensional grain mapping by x-ray diffraction contrast tomography and the use of Friedel pairs in diffraction data analysis[END_REF]. The reconstructed volumes were re-scaled to match the 0.7 µm voxel size of the tomographic reconstructions for ease of comparing the two datasets.

For tomography, optics giving an effective pixel size of 0.7 µm were used. A multilayer monochromator was used to select a beam with an energy of 20.5 keV and a bandwidth of E/E=10 -2 for optimal flux. The distance between the sample and detector was set to 40 mm in order to obtain some phase contrast in the reconstructed 3D images, a mode of imaging which is known to enhance the visibility of sub-pixel sized defects [START_REF] Ludwig | Study of the interaction of a short fatigue crack with grain boundaries in a cast Al alloy using X-ray microtomography[END_REF]. A series of 1000 radiographs taken during a 180° rotation of the sample was used to reconstruct the 3D volumes using a filtered back projection algorithm. The specimen was first characterised by DCT to provide a 3D grain map. Five fine notches to initiate cracks were made at different locations at the centre of gauge section by FIB milling. For three of these, an attempt was made to place notches in specific grains, selected by inspection of the DCT data, where the notch was oriented parallel to the basal plane and was therefore inclined to the tensile axis. Part of the sample was removed by gentle manual polishing to expose the selected grains at the surface. Two further notches, normal to the loading axis, were placed in arbitrary grains. Each notch was approximately 90 µm long, with an estimated depth of 30 µm, with a mouth opening of ~10 µm.

The fatigue loading was performed using a cyclic loading machine designed for insitu tomographic measurements (i.e. compact size and minimal vibration levels). The load train is supported by an X-ray transparent Perspex tube to allow images to be acquired through a 180º rotation of the sample [START_REF] Buffière | Three dimensional visualisation of fatigue cracks in metals using high resolution synchrotron X-ray microtomography[END_REF]. The tension-tension fatigue loading used a sinusoidal waveform at a frequency of 5 Hz with a maximum stress of 130 MPa and R-ratio (minimum/maximum load) of 0.25. The cycling was interrupted at intervals for tomographic scanning at 500, 1000, 2000, 2500, 3000, 3500, 4000, 5000 cycles, then at intervals of 250 to 500 cycles up to 10250 cycles when the test was suspended. The scans were performed with the sample under maximum load to increase the crack opening, and thus its visibility.

Results and analysis

Although five notches were prepared, only one notch exhibited significant fatigue crack development. This notch was aligned perpendicular to the stress axis, and was inserted in an arbitrary grain. It is the results from this notch that are presented here.

Figure 2 shows sections through the reconstructed volume at the position of the notch, after 500 and after 10500 cycles. The crack can be observed to have begun to propagate from the notch in the first observation at 500 cycles. Both the crack and small, highly attenuating Nd-rich precipitates are resolved. Although it is possible to use automatic segmentation of these images to visualise the 3D shape of the crack, this is insufficient for the purpose of crack dimension measurements due to errors caused by the relatively low attenuation contrast between the crack and the magnesium matrix and also because of phase contrast artefacts. Consequently, the position of the fatigue crack front was traced manually by inspection of individual 2D slices from the 3D datasets to determine the coordinates of the crack front as it developed throughout the experiment. The error in determining the position of the crack front was estimated to be five voxels (i.e. 3.5 µm). An algorithm (described in the Appendix) was implemented in Matlab (Matlab 7.4, The Mathworks Inc, Natick, MA) to determine the crack growth velocity as a function of position, using these data.

These tomographic reconstructions were aligned with the grain map determined by DCT, as were the manually determined crack front data. Figure 3 shows a rendering of a region of interest around the crack in which the grains have been rendered partially transparent. This emphasises the 3D nature of the crack, which deviates strongly from the x-y plane.

2D and 3D representation of data

The 3D representations of the data contain a very large amount of information, and hence are complex to view, complicating the task of describing the behaviour of the crack as it interacts with the microstructure. It is useful to obtain a 2D visualisation of the crack. An example is shown in figure 4, which is a projected view of the crack looking along the tensile axis (z), which is normal to the x-y plane. The free surface of the sample and the FIB notch are on the left hand side of the view. The projected crack has an approximately semi-elliptical shape. The small "holes" that appear in the crack (labelled with arrows in (a)) are regions where the opening displacement is too small for the crack to be segmented using simple threshold methods. These may be removed by image processing methods. Each pixel in the projected image of the crack was defined as a crack element (with a size of 0.7 x 0.7 µm), which could have various parameters associated with it such as the grain orientation at its corresponding X-Y-Z position in the 3D dataset, the crack plane or the crack growth rate, for example. Although the crack exhibits some branching (see Fig 3), the branches do not tend to overlap vertically. In these cases, the element was assigned to the largest overlapping branch. Such areas where there may be more than one crack element occupy less than 20% of the total surface. This projected view of the crack may be coloured to display the parameters at each element. For example, Figure 4 is coloured according to the orientations of the grains the crack passes through. The colour is defined by the grain basal plane orientation relative to the sample. The colour scheme, which is also used in the 3D rendering in figure 3, is illustrated for a stereographic projection in figure 4(c). The angular rotation of the basal plane normal about the z-axis determines the colour, while the saturation of the colour, from white to near-black, depends on the angle between the basal plane normal and the z-axis. Grains in which the basal plane is perpendicular to the tensile axis are lightest.

To the left and right of the notch (figure 3), the crack grows along approximately flat planes that are inclined to the notch, whereas as the crack propagates away from the notch into the sample (towards the rear of figure 3), it splits into several "fingers" lying on approximately parallel planes whose normal vectors are inclined at approximately 30° to the tensile axis. Excluding the area of the notch, approximately 50% of the crack surface is inclined by more than 30° from the x-y plane. The variations in the crack plane can be quantified. To do this, at each element in the projected crack, the 3D positions of the voxels of the crack surface associated with a region of interest of 8x8x8 µm 3 around the element were selected. The local crack plane was obtained by fitting to these points. The local crack plane, displayed on the 2D projection of the crack with the same orientation-colour map used for figures 3 and 4, is shown in figure 5.

The 3D information about the crack plane and the basal plane orientations can be visually compared using these 2D images. For example, comparison of Figure 4 and 5 reveals a significant proportion of corresponding regions of similar colour, indicating that the crack has a tendency to propagate on the basal plane. This is shown with greater sensitivity in Figure 6, which presents the angular deviation between the crack plane and the local basal plane. The low values identify the regions of basal propagation. Higher values tend to be observed around the periphery of the crack, and where the crack intersects certain grain boundaries. Excluding the notch, approximately 54% of the crack surface is inclined by less than 30° from the basal plane.

In the single crystal work of Iida et al, twinning was observed to interact with crack propagation [24]. The crack morphology obtained in this work does not reveal evidence of twinning. However, the resolution is insufficient to detect fine twins, and it is possible that these may occur.

Crack growth rate

Data, extracted from the tomographic reconstructions, for the total surface length (2c) and the projected depth (a) for the observed crack are summarised in figure 7. The extent of propagation is insufficient for fit anything more complex than a linear relationship, and the average crack increment per cycle, N, (dc/dN and da/dN) are 8 and 7 nm/cycle respectively.

When visualised in 3D the crack growth behaviour is complex. Figure 4(a) shows the projections of the manually segmented crack fronts at increments of ~1000 cycles.

Because the surface of the crack is non-planar, the crack growth rate cannot be accurately visualised directly from these projections. However, they provide a qualitative indication of the crack propagation, showing that the crack front is retarded at the surface in the region of grains G and H. It is clear that surface observations of crack growth rates can be inadequate to describe its 3D development.

The algorithm described in the Appendix was used to calculate the local growth rate at each point along the observed crack fronts. The data were smoothed for visualisation by presenting the median value of the nearest ten crack growth rate data points for each position on the crack surface. The resulting data for local growth rate as a function of position are shown in Figure 8. This shows the crack growth rate varies considerably with position, varying between 4 and 40 nm/cycle. Several regions of fast crack growth are clearly observed on either side of the originating FIB notch. The average rate is around 9 nm/cycle, which is consistent with the values obtained from the projected dimensions of the crack. The crack growth rate shows no tendency to increase with distance from the notch, and indeed tends to be lower with increasing distance. Slow crack growth is observed particularly in grains I, J and L at the top and right of figure 4.

Discussion

The average fatigue crack growth rate extracted from the observations is around 7-9 nm/cycle. Comparative data on short fatigue crack propagation in magnesium alloys (i.e. crack length<1mm) are quite limited. Wang et al [30,31] The combined use of the CT and DCT datasets for the fatigue crack and microstructure helps to reveal its characteristic behaviour. There is a tendency for the crack to be oriented with the basal plane, and those regions where this is the case tend to exhibit the higher growth rates (see, for example, grains D and F). Lower rates appear to be associated with the non-crystallographically oriented regions that developed towards the end of the experiment, particularly in grains D, I, J and K. The crack therefore has the general appearance of a predominantly stage I crystallographic fatigue crack. Such cracks are known to be exhibit variable crack growth rates [eg 9] as is observed here.

Closer examination of the crack growth rates shows regions of local growth rate variation. Referring to the grain labels in figure 4, some regions of interesting behaviour are as follows: Fast, crystallographic crack growth occurs in grains A, D and F, compared with slower crystallographic growth in B. Fast, non-crystallographic growth is seen in grain H. Slower crack growth occurs around grain C, between grains E and L, and along the crack front in grains I, J and K. The unique data obtained by combined DCT and CT allows the crack tip environment at these locations to be examined to determine whether these may be understood in terms of the Schmid factor in the grains, the character of the grain boundaries, and also the orientation of the slip planes relative to the crack initiating notch.

Grain Schmid factor

The Schmid factor describes the fraction of the shear stress resolved onto each slip system to an applied tensile stress. Stage I short fatigue crack growth is driven by slip, and therefore the maximum Schmid factor may represent the driving force for crack growth [33]. In single crystals, short fatigue crack growth has been observed to occur on the slip planes with the highest Schmid factor [33]. The Schmid factor of the most highly stressed slip system in each grain through which the crack passed is shown in Figure 9. For this analysis, the applied stress has been assumed to be uniform and parallel to the loading axis, so each grain is effectively treated as an isolated single crystal. This neglects both the effect of the notch or crack on the local stress field, and also the effect of the elastic anisotropy of the grains. Both of these factors will complicate the local stress field, but are neglected as a first approximation. It should be noted that because of the symmetry of the hexagonal system there is only one plane available for easy slip, the basal plane. Therefore the Schmid factor can vary between 0 and 0.5 depending on the orientation of the grain.

In other crystal systems, for example cubic structures, the higher multiplicities of slip planes means that that some will always have non-zero Schmid factors. For example in face-centred cubic systems it can be shown that the minimum Schmid factor is ~0.27.

From Figure 9 it can be seen that the Schmid factor is high (greater than 0.25) in most of the grains (i.e. A-G, K). In particular, the crack grows quickly across grain F (Schmid factor 0.46), before dramatically slowing on reaching grains I and J (Schmid factors 0.05 and 0.25 respectively). In this region of reduced crack growth rate, the crack plane tends to be normal to the applied load, which is a characteristic of stage II crack growth.

Interestingly, in grain H, which lies between grain F and the sample surface, quite fast crack growth is observed despite its low Schmid factor (0.16). Examination of Figure 6 shows that the crack does not grow on the basal plane in this grain. It appears, from the successive positions of the crack front, that the faster crack growth in grain F causes grain H to develop a ligament at the sample surface. Such behaviour is illustrated in [34]. The ligament subsequently fails at a higher growth rate (presumably driven by the increased crack opening around it as the overall crack length increased). The crack path through the ligament appears to be defined by the earlier failure of grain F, as indicated in Figure 5.

The implication of this simplistic analysis is that the Schmid factor may provide a first order approximation to the tendency for stage I crack growth. Local arrangements of the crack geometry may, however, have additional effects, such as those observed for grain H.

Tilt/twist character of the grain boundaries

It has been argued for many years that stage I fatigue cracks should be retarded or arrested at grain boundaries where the slip plane deviates [35,36]. Recently, Zhai and Wilkinson used two angular measures (twist and tilt), to characterise the misorientation between slip planes at grain boundaries in a 2D study of crack growth [START_REF] Zhai | A crystallographic mechanism for fatigue crack propagation through grain boundaries[END_REF]. In their nomenclature, twist is the angle between the traces of the two slip planes on the boundary plane, and tilt is the angle between the traces of the slip planes on the surface of the sample. In their investigation, the boundary plane was assumed to be normal to the surface plane. The twist angle was identified as being the more important parameter to describe the resistance to crack propagation. Re-initiation at multiple sites along the intersection between the crystallographic crack plane and the grain boundary would tend to cause crack branching, with a consequent reduction in crack tip strains due to mechanical shielding from the intact ligaments of material between these parallel cracks.

In this work, the orientation of the grain boundary can be obtained in 3D by identifying those voxels lying on the boundary between two grains, and then fitting a plane through these. The twist angle can then be calculated readily. The trace of the slip plane on the surface of the sample has no meaning within the 3D volume. A more general definition of the tilt angle is required. This was done at each grain boundary by determining the necessary twist angle, as a rotation about the grain boundary normal, to bring the traces of the slip planes on the boundary plane into alignment. The tilt angle is then the remaining angle between the slip planes, about a rotation axis that lies in the plane of the grain boundary. These stages are illustrated in figure 10.

For basal plane slip in the hexagonal crystal system, the maximum tilt and twist values can vary between 0 and 90º. Values for those boundaries traversed by the crack are shown in Figure 11. It is can be seen that several regions of slow crack propagation (i.e. between grains E and L and around grain C) are associated with boundaries of either high twist or tilt values (Between E and L, tilt of 60°, twist of 88°, between C and D tilt of 86°, twist of 71°). High twist values are also observed between grain F and grains I and J (78° and 72° respectively) where the crack decelerates. These observations support the hypothesis that stage I crack propagation is retarded at boundaries with large misorientations of the slip plane.

Compatibility between the notch and the slip plane

The notch initiates the crack growth and its strain field influences the early stage I crack growth. It is interesting to assess the compatibility of the slip planes with the notch, particularly in the grains within which the notch lies. Here we define a quantity, the notch twist angle, to describe this. The notch is horizontal, with the shape shown in the figures. The notch tip opening angle is here assumed to be 0°. The ideal crack growth direction is assumed to be radial from the edges of the notch. Thus a plane may be defined that is normal to the growth direction at each point. The angle between the trace of the notch (in this case horizontal) and the trace of the slip plane in the grain is calculated on this plane. The analysis is equivalent to that for the tilt/twist calculation at grain boundaries, with the plane of the notch replacing one of the slip planes. The twist angle obtained will be referred to as the notch twist angle. The results of this calculation for the grains in which the crack propagated are shown in figure 12. The large variations of the notch twist angle within grains D and F are a consequence of the slip plane normal being close to the radial crack growth direction.

The effect of the notch twist angle might be responsible for the slow growth in grain B, despite its high Schmid factor of 0.45. In parts of grains D and F (marked with arrows in figure 12), the notch plane is close to the basal plane. Particularly in grain F, rapid early crack propagation is observed starting from the compatible region. However, throughout grain B the basal plane is twisted with respect to the notch, as shown by the relatively high value of the notch twist angle (>30°). The crack initiating from the notch within grain A therefore splits into several parallel fingers in order to propagate locally on the basal planes. There is a twist of around 15° between the plane of the notch and the basal plane in grain A, followed by a further 14° twist between grain A and grain B. The effect on the crack growth rate is thus similar to that of two successive twist boundaries. In contrast, for grains F and D the crack propagates more rapidly on a single plane for most of the grain, only deviating near the edges.

This description is useful only close to the notch. In grain B the crack growth rate tends to increase as the crack develops in size. Here, it is judged that the high Schmid factor allows fast crack growth as the incompatibility to the notch, and hence the effect of the notch twist angle, becomes less significant with distance.

Perspectives

This is a preliminary study, making use of the new technique of diffraction contrast tomography. It demonstrates that the crack growth rates are sensitive to the local environment of the grains and the geometry of the crack, in particular that stage I crack propagation is retarded at grain boundaries with large misorientations of the basal slip plane. It reveals the potential for studies to observe the behaviour of longer cracks, and their interaction with obstacles in the microstructure. In future studies, it would be interesting to observe how the crack interacts with low Schmid factor grains (e.g. grain I). If the crack arrests at such a grain boundary and then subsequently reinitiates, studies of this type might identify the mechanism by which the obstacle is overcome, whether by local twinning, prismatic slip, or another failure mode. This would extend the work of Iida et al to a polycrystalline context [24]. The behaviour at boundaries with high twist and low tilt angles and the effects of texture in wrought material are also of interest. A logical direction to pursue will be finite element modelling to estimate the true local crack tip strains and the correlation with growth rates. It also may be interesting to use 3D image correlation to study the strains around the crack tip, using the precipitates as markers.

Conclusion

The combined use of diffraction contrast tomography and synchrotron microtomography provides unique opportunities for the in-situ study of short fatigue cracks, since it gives high-resolution access to the grain structure and grain boundary characteristics, the crack plane and the local crack growth rate. This study shows that by using these microstructure parameters, a qualitative explanation for the observed 3D behaviour of a short stage I fatigue crack can be obtained.

Appendix: A Crack Advance Measurement Algorithm Abstract

By meshing two subsequent measured crack front positions it is possible to map the connections between points on the two lines, which allows the crack velocity to be estimated as a function of position. The meshing algorithm proceeds by producing triangles with sides as short as is practically achievable in order to connect the two lines in the likely direction of crack growth.

Method

Given a set of points defining two measured lines, at successive time steps t and t+1, a search may be made for the nearest point on line t+1 to each point on line t, as shown in figure A1.

A similar search is then performed in reverse to find the nearest point on line t to each point on line t+1, as shown in figure A2. Next, the common pairs of points in both searches that are mutual nearest neighbours are identified, as shown in Figure A2. These connections are used as seeds for the triangulation of the mesh.

Figure A3. Triangulation proceeding from pairs identified as close contacts.

The lines connecting these mutually nearest neighbours are identified and are used as an initial set to divide up the space along the two crack fronts. The regions between these lines are then triangulated by comparing the four possible lines that could be drawn to add a new triangle to the interior of each region, and adding the one which is the shortest, as shown in Figure A3. This is repeated for each region bounded by two close connections, or for regions at the starts or ends of the lines beginning with one close connection and advancing away from it.

The velocity at each point on the crack is then estimated using the longest line that begins at a particular point, as shown in figure A4, and the time between the two observations of the crack.

To increase the accuracy of the method, and to give connections that pass between the two lines closer to the normal angle (which is assumed to be most representative of the local growth rate), the initial point set is resampled to contain the original points plus a set of additional points distributed in between them in those regions where the average point spacing is greater than the average for the set. The facility exists to vary the criterion for adding points so as to resample segments greater than 1/n times the average, where n is user definable. In the analysis of data in this paper, n=1.

Figure 1 .

 1 Figure 1. Miniature fatigue sample mounted on beam line.

Figure 2 .

 2 Figure 2. Reconstructed sections after (a) 500 and (b) 10,500 fatigue cycles, showing the crack growing from the FIB notch, and the microstructure of the material. The white particles are Nd-rich precipitates.

Figure 3 .

 3 Figure 3. 3D rendering showing the crack (white) and the grains (coloured, semi-transparent). Colour denotes basal plane orientation according to the key in figure 4.

Figure 4 .

 4 Figure 4. (a) Plan view of crack from along the loading (z) axis. The crack is coloured according to the basal plane orientation of the grains intersected by the crack. Black lines mark the grain boundaries. Superimposed white lines are the manually segmented crack fronts (increments of ~1000 cycles). The black lines are the grain boundary intersections with the surface. (b) The labels A-L that are used in the text to describe the grains. (c) The colour map used to display orientation (stereographic projection).

Figure 5 .

 5 Figure 5. The local orientation of the crack plane. Colour represents the plane orientation relative to the sample, as in Figure 4.

Figure 6 .

 6 Figure 6. Deviation between the crack plane and the crystallographic basal plane.

Figure 7 .

 7 Figure 7. Development of crack surface length and depth with number of fatigue cycles extracted from tomography observations.

Figure 8 .

 8 Figure 8. Local crack growth rate. Note that while the maximum observed rate was 40 nm/cycle, the colour scale varies between 0 and 20 nm/cycle to better reveal local variations.

Figure 9 .

 9 Figure 9. Schmid factor for the most activated basal slip system in each grain.

Figure 10 .

 10 Figure 10. The calculation of twist and tilt angles in 3D. (a) Twist is the angle between the traces of the slip planes on the grain boundary. (b) Tilt is the remaining angle after applying the twst angle to one of the slip planes.

Figure 11 .

 11 Figure 11. (a) Tilt and (b) twist angles for boundaries traversed by the crack.

Figure A1 .

 A1 Figure A1. Initial search for nearest points on line t+1 to points on line t.

Figure A2 .

 A2 Figure A2. Initial search for nearest points on line t to points on line t+1.

Figure A4 .

 A4 Figure A4. Longest triangle sides identified as giving crack velocity between two time steps.

Figure 1 .

 1 Figure 1. Miniature fatigue sample mounted on beam line.

Figure 2 .

 2 Figure 2. Reconstructed sections after (a) 500 and (b) 10,500 fatigue cycles, showing the crack growing from the FIB notch, and the microstructure of the material. The white particles are Nd-rich precipitates.

Figure 3 .

 3 Figure 3. 3D rendering showing the crack (white) and the grains (coloured, semi-transparent). Colour denotes basal plane orientation according to the key in figure 4.

Figure 4 .

 4 Figure 4. (a) Plan view of crack from along the loading (z) axis. The crack is coloured according to the basal plane orientation of the grains intersected by the crack. Black lines mark the grain boundaries. Superimposed white lines are the manually segmented crack fronts (increments of ~1000 cycles). The black lines are the grain boundary intersections with the surface. (b) The labels A-L that are used in the text to describe the grains. (c) The colour map used to display orientation (stereographic projection).

Figure 5 .

 5 Figure 5. The local orientation of the crack plane. Colour represents the plane orientation relative to the sample, as in Figure 4.

Figure 6 .

 6 Figure 6. Deviation between the crack plane and the crystallographic basal plane.

Figure 7 .

 7 Figure 7. Development of crack surface length and depth with number of fatigue cycles extracted from tomography observations.

Figure 8 .

 8 Figure 8. Local crack growth rate. Note that while the maximum observed rate was 40 nm/cycle, the colour scale varies between 0 and 20 nm/cycle to better reveal local variations.

Figure 9 .

 9 Figure 9. Schmid factor for the most activated basal slip system in each grain.

Figure 10 .

 10 Figure 10. The calculation of twist and tilt angles in 3D. (a) Twist is the angle between the traces of the slip planes on the grain boundary. (b) Tilt is the remaining angle after applying the twst angle to one of the slip planes.

Figure 11 (

 11 Figure 11 (a) Tilt and (b) twist angles for boundaries traversed by the crack.

Figure 12 .

 12 Figure 12. The notch twist angle for all grains traversed by the crack, assuming radial growth from the edges of the notch.

Figure A1 .

 A1 Figure A1. Initial search for nearest points on line t+1 to points on line t.

Figure A2 .

 A2 Figure A2. Initial search for nearest points on line t to points on line t+1.

Figure A3 .

 A3 Figure A3. Triangulation proceeding from pairs identified as close contacts.

Figure A4 .

 A4 Figure A4. Longest triangle sides identified as giving crack velocity between two time steps.

Table 1 : Chemical Composition in weight percent.

 1

	Zn	Zr	Nd	Gd	Other rare earth Others Mg
	0.31 0.5	2.5	1.3	0.1	<0.01	Balance

Table 2 : Mechanical Properties

 2

	0.2% Proof Stress	Tensile Strength	Elongation
	(MPa)	(MPa)	(%)
	175	265	7

  report data for cast AM50 Mg-Al alloy in vacuum (yield stress 130 MPa, tensile strength 200 MPa). At constant stress amplitude (140 MPa) at room temperature the crack propagation could be described by a constant growth rate of approximately 300 nm/cycle for a crack length of less than 300 µm. At stress amplitude of 120 MPa, the crack growth rate averaged 70 nm/cycle below 150 µm, and 200 nm/cycle between 150 µm and 300 µm.Gall et al [32] present measurements made in a cast AM60B Mg-Al alloy (yield stress 150 MPa, tensile strength 285 MPa). Growth rates at a stress amplitude of 90 MPa were constant for crack lengths between 50 and 300 µm, varying between 7 to 12 nm/cycle in vacuum and 18-30 nm/cycle in water vapour. The observed average growth rates of around 8 nm/cycle in ambient air in this work are therefore consistent with such data. In both cases transgranular cracking was the predominant mode, having a crystallographic character associated with slip planes.

Table captions Table 3 : Chemical Composition in weight percent.Table 4 : Mechanical Properties

 captions34

Acknowledgements

This work was supported by the UK Research Council funded project "Crack Nucleation and Short Crack Behaviour by Non-Destructive In Situ Observation" (EP/C002946/1). AAK is grateful to the Higher Education Commission of Pakistan for supporting his studies at Manchester University. Much of the analysis was carried out while AK was employed by GKSS (now HZG), Geesthacht, Germany. The material was kindly supplied by Magnesium Elektron Ltd. The work was conducted while TJM was at the Materials Performance Centre, School of Materials, University of Manchester.


Katell Mevel

Béatrice Desgranges

Jean-Claude Baron

Brigitte Landeau

Vincent De

La Sayette

Fausto Viader

Francis Eustache

Gaël Chételat

Which SPM method should be used to extract hippocampal measures in early Alzheimer's disease?

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.


Magali Lécureux-Mercier

Conservation laws with a non-local flow Application to Pedestrian traffic

HAL is

Introduction

In the last decades, the crowds' dynamics has attracted a lot of scientific interest. A first reason of this interest is to understand how crowd disasters happened in some panic events, for example at the end of football play, during concerts, in the case of fire, or in place of pilgrimage (e.g. on Jamarat Bridge in Saudi Arabia, see [START_REF] Helbing | Dynamics of crowd disasters: An empirical study[END_REF]). Another reason lies in architecture of buildings such as subway stations or stadiums, where a lot of pedestrians are crossing. The goals here are consequently twofold: in one hand we want to understand the behavior of pedestrians in panic and adapt the regulation of traffic in order to avoid deaths; in the other hand, we want to modelize the interaction of several kinds of pedestrians with different objectives and in particular study how the geometry influences the general pattern.

In a macroscopic setting, a population is described by its density ρ which satisfies the conservation law

∂ t ρ + Div(ρ V (t, x, ρ)) = 0 , ρ(0, •) = ρ 0 , (1)
where V (t, x, ρ) is a vector field describing the velocity of the pedestrians depending on the time t ≥ 0, the space x ∈ R N and the density ρ. According to the choice of V , various behaviors can be observed. Several authors already studied pedestrian traffic in two dimensions space (N = 2). Some of these models are local in ρ, that is to say V depends on the local density ρ(t, x) [START_REF] Coscia | First-order macroscopic modelling of human crowd dynamics[END_REF][START_REF] Bellomo | On the modelling of traffic and crowds -a survey of models, speculations, and perspectives[END_REF][START_REF] Hughes | A continuum theory for the flow of pedestrians[END_REF][START_REF] Hughes | The flow of human crowds[END_REF][START_REF] Maury | A macroscopic crowd motion model of the gradient-flow type[END_REF][START_REF] Maury | Handling congestion in crowd motion modeling[END_REF] ; other models use not only the local density ρ(t, x) but the entire distribution of ρ, for example they depend on the convolution product ρ(t) * η [START_REF] Di Francesco | On the Hughes' model for pedestrian flow: the one-dimensional case[END_REF][START_REF] Piccoli | Time-evolving measures and macroscopic modeling of pedestrian flow[END_REF]. Here, in the line of preceding papers [START_REF] Colombo | Control of the continuity equation with a non local flow[END_REF][START_REF] Colombo | A class of non-local models for pedestrian traffic[END_REF][START_REF] Colombo | An analytical framework to describe the interactions between individuals and a continuum[END_REF][START_REF] Colombo | Nonlocal crowd dynamics models for several populations[END_REF], we present nonlocal macroscopic models for pedestrian traffic, we study these models and compare their properties. Our first aim is to modelize the behavior of pedestrians in different situations: crowd behaves indeed differently in panic or in a normal situation where courtesy rules do apply. We also introduce models in the case of a population interacting with an individual, and in the case of several populations with different objectives. For instance, we want to include in our study the case of two populations crossing in a corridor. 1 Technion, Israel Institute of Technology, Amado Building, 32000 Haifa, Israel 1 Second, we want to study the introduced models and prove existence and uniqueness of solutions under various sets of hypotheses. We will use two kinds of arguments: the first one comes from Kružkov theory [START_REF] Kružkov | First order quasilinear equations with several independent variables[END_REF][START_REF] Lécureux-Mercier | Improved stability estimates on general scalar balance laws[END_REF], the second one from the optimal transport theory. We want to prove existence and uniqueness of solutions for the various models presented below. Let us concentrate on the case of pedestrians in panic:

∂ t ρ + Div(ρ v(ρ * η) ν(x)) = 0. (2)
All the proofs of existence and uniqueness in this note are based on the following idea: let us fix the nonlocal term and, instead of (2), we study the Cauchy problem

∂ t ρ + Div(ρ v(r * η) ν(x)) = 0 , ρ(0) = ρ 0 , (3)
where r is a given function. Then, we introduce the application

Q : r → ρ X → X , (4)
where the space X has to be chosen so that (a) X is equipped with a distance d that makes X complete;

(b) the application Q is well-defined: the solution ρ ∈ X exists and is unique (for a fixed r);

(c) the application Q is a contraction.

Once we have fullfilled these conditions, we can prove existence of a solution using a fixed point argument. Note that, in the modelization of pedestrian traffic, the space dimension N has to be equal to two, but our results are in fact true for all N ∈ N. For instance, they can be adapted in dimension N = 3 to modelize the behavior of fishes or birds. Consequently, we keep here a general N , even if we essentially think to the case N = 2.

This note is organized as follows: in Section 2, we describe some nonlocal models and their properties. In Section 3 we study one of these models through Kružkov theory and in Section 4 we study the same model through optimal transport theory.

Pedestrian Traffic Modelization

One-Population model

Pedestrian in panic

The first model we present corresponds to pedestrians in panic and was studied in [START_REF] Colombo | Control of the continuity equation with a non local flow[END_REF][START_REF] Colombo | Nonlocal crowd dynamics models for several populations[END_REF], in collaboration with R. M. Colombo and M. Herty. A panic phenomenon appears under special circumstances in crowded events. In these cases, the people are no longer rational and try, no matter how, to reach their target. Let us denote ρ(t, x) the density of pedestrians at time t and position x ∈ R N . We consider the Cauchy problem:

∂ t ρ + Div ρ v(η * ρ(t, x)) ν(x) = 0 ; ρ(0, •) = ρ 0 . (5)
Here, v is a real function describing the speed of the pedetrians. This function does not depend on the local density ρ(t, x) but on the averaged density ρ(t) * η(x) = R N ρ(t, x -y) η(y) dy.

The vector field, ν(x) describes the direction that the pedestrian located in x will follow, independently from the distribution of the pedestrians' density. Note that we are working here on the all of R N and not on a subset of R N ; thus, we are not working on a restriction to a room, for example. However, we can still introduce the presence of walls and obstacles in the choice of the vector field ν. Let us denote Ω ⊂ R N the space where the pedestrians are authorized to walk, e.g. a room. If we choose ν(x) in a nice way (for example we can require that on the walls, i.e. for all x ∈ ∂Ω, ν(x) coincides with the entering normal to Ω), then we can conclude to the invariance of the room. More precisely, if the initial density has support on some closed set Ω ⊂ R N , then the solution will have support contained in Ω for all time. This remark allows us to avoid considering any boundaries and to have solutions on all R N . Using the Kružkov theory on classical scalar conservation laws, we are able to prove:

Theorem 2.1 (see [6]). Let ρ 0 ∈ (L 1 ∩ L ∞ ∩ BV)(R N , R + ). Assume v ∈ (C 2 ∩ W 2,∞ )(R, R), ν ∈ (C 2 ∩ W 2,1 )(R N , R N ), η ∈ (C 2 ∩ W 2,∞ )(R N , R).
Then there exists a unique weak entropy 5) with initial condition ρ 0 . Furthermore we have the estimate

solution ρ = S t ρ 0 ∈ C 0 (R + , L 1 (R N , R + )) to (
ρ(t) L ∞ ≤ ρ 0 L ∞ e Ct , (6)
where the constant C depends on v, ν and η.

For the definition of weak entropy solutions see Section 3; the proof is defered to Section 3.2.1. Note that in Theorem 2.1, the hypotheses are very strong. Let us denote P(R N ) the set of probability measures on R N and M + (R N ) the set of positive measures on R N . In collaboration with G. Crippa, using now some tools from optimal transport theory, we obtained the better result: Theorem 2.2 (see [START_REF] Crippa | Existence and uniqueness of measure solutions for a system of continuity equations with non-local flow[END_REF]). Let

ρ 0 ∈ M + (R N ). Assume v ∈ (L ∞ ∩ Lip)(R, R), ν ∈ (L ∞ ∩ Lip)(R N , R N ), η ∈ (L ∞ ∩ Lip)(R N , R + ).
Then there exists a unique weak measure solution ρ ∈ L ∞ (R + , M + (R N )) to [START_REF] Colombo | Control of the continuity equation with a non local flow[END_REF] with initial condition ρ 0 .

If furthermore ρ 0 ∈ L 1 (R N , R + ) then for all t ≥ 0, the solution ρ satisfies also

ρ(t) ∈ L 1 (R N , R + ).
For the definition of weak measure solution see Section 4; the proof is defered to Section 4.2.

Note that for the model ( 5), there is a priori no uniform L ∞ bound on the density. Indeed, heuristically, considering the case in which the density is maximal, equal to 1, on the trajectory of a pedestrian located in x. If the averaged density around x is strictly less than 1 (because, for example there is no one behind this pedestrian), then the speed v(ρ * η) will be strictly positive, which means the pedestrian in x will try to go forward, even though there is a queue in front of him. Consequently, we expect the density to become larger than one. This behavior is not really unexpected in the case of panic. In fact, in some events the density attained up to 10 persons per square meter, which is obviously too much and a cause of deaths (see [START_REF] Helbing | Dynamics of crowd disasters: An empirical study[END_REF]). Consequently, it is quite satisfactory to recover this behavior. One of our goal in this context is then to introduce a cost functional allowing to characterize the cases in which the density is too high, and to find extrema of this functional. Let us introduce

J T (ρ 0 ) = T 0 Ω f (S t ρ 0 ) dx ,
where Ω ⊂ R N is the room, ρ 0 is the initial condition and S t ρ 0 is the semi-group generated by Theorem 2.1. We choose the function f ∈ C 1 (R, R + ) so that it is equal to zero for any density ρ less than a fixed threshold ρ c and so that it is stictly increasing on [ρ c , +∞[. Consequently, the functional J T above allows to characterize the solutions with too high density and in particular it vanishes if the set {(t, x) ∈ [0, T ] × R N : S t ρ(x) ≥ ρ c } has measure zero. We are then interested in finding the minima of this cost functional.

Using the Kružkov theory we prove the following differentiability result:

Theorem 2.3 (see [START_REF] Colombo | Control of the continuity equation with a non local flow[END_REF][START_REF] Colombo | Nonlocal crowd dynamics models for several populations[END_REF]).

Let ρ 0 ∈ (W 2,∞ ∩ W 2,1 )(R N , R + ), r 0 ∈ (W 1,1 ∩ L ∞ )(R N , R) and denote ρ = S t ρ 0 . Assume v ∈ (C 4 ∩ W 2,∞ )(R, R), ν ∈ (C 3 ∩ W 2,1 )(R N , R N ), η ∈ (C 3 ∩ W 2,∞ )(R N , R + ).
Then there exists a unique weak entropy solution r = Σ ρ t r 0 to the Cauchy problem

∂ t r + Div(r v(ρ * η) ν(x)) = -Div(ρ v ′ (ρ * η) ν(x)) , r(0) = r 0 . (7)
Furthermore, the semi-group S t obtained in Theorem 2.1 is Gâteaux-differentiable, that is to say, for all

ρ 0 ∈ (W 2,1 ∩ W 2,∞ )(R N , R + ), r 0 ∈ (W 1,1 ∩ L ∞ )(R N , R), lim h→0 S t (ρ 0 + hr 0 ) -S t ρ 0 h -Σ ρ t r 0 L 1 = 0 .
The proof is defered to Section 3.2.2. It is not possible to obtain the same result by the use of optimal transport theory. Indeed, it seems already not possible to find a good definition of Gâteaux differentiability on the set of probability measures equipped with the Wasserstein distance of order 1.

Orderly crowd

In opposite to the previous model, for a model of orderly crowd it is required to have a uniform L ∞ bound on the density. In collaboration with R. M. Colombo and M. Garavello [START_REF] Colombo | A class of non-local models for pedestrian traffic[END_REF], we studied the equation

∂ t ρ + Div   ρv(ρ) ν(x) - ∇(ρ * η) 1 + ∇(ρ * η) 2    = 0 ; (8)
with

ρ 0 ∈ (L 1 ∩ L ∞ ∩ BV)(R N ; R).
In this model, the speed v depends on the local density ρ(t, x), which allows to prove some uniform bound in L ∞ . The prefered direction of the pedestrians is still ν(x), but they deviate from their optimal path trying to avoid entering regions with higher densities. Indeed, (ρ * η) is an average of the crowd density around x and -∇(ρ * η) is a vector going in the direction opposite to the area of maximal averaged density. Due to the nonlinearity of the flow with respect to the local ρ, it is no longer possible to use optimal transport theory. Hence, we use Kružkov theory, to prove existence and uniqueness of solutions. We obtain the theorem: Theorem 2.4 (see [START_REF] Colombo | A class of non-local models for pedestrian traffic[END_REF]).

Assume that v ∈ C 2 ([0, 1], R + ) satisfies v(1) = 0, that ν ∈ (C 2 ∩ W 2,1 ∩ W 1,∞ )(R N , R N ), and that η ∈ (C 3 ∩ W 3,1 ∩ W 2,∞ )(R N , R). Then, for any ρ 0 ∈ (L 1 ∩ L ∞ ∩ BV)(R N , [0, 1]), there exists a unique weak entropy solution ρ ∈ C 0 (R + , L 1 (R N , [0, 1])).
The proof of this theorem relies on Kružkov theory (see Section 3). Note that for this model the density is uniformly bounded in L ∞ , contrarily to the panic model, in which the L ∞ norm can grow exponentially in time.

Some difficulties now appear in proving that the pedestrians remain in the authorized area. Let us introduce the following invariance property:

(P) Let Ω ⊂ R N be region where the pedestrians are allowed to walk. The model ( 8) is invariant with respect to Ω if

Supp(ρ 0 ) ⊂ Ω ⇒ Supp(ρ(t)) ⊂ R N for all t ≥ 0 . (9)
To obtain that (P) is satisfied, we have to require the prefered direction to be strongly entering the room (see [4, Proposition 3.1 & Appendix A]). Some interesting phenomena show up through numerical computations. First, when considering a crowd walking along a corridor, we observe the formation of lanes (see Figure 1). Furthermore, this phenomenon seems very stable with respect to initial conditions and geometry. Indeed considering the room and initial distribution as in Figure 2, adding some various obstacles, we still have lanes, at least in large space (see Figures 3). A further remarkable property of the model ( 8) is that it captures the following wellknown, although sometimes counter intuitive phenomenon (Braess paradox): the evacuation time through an exit can be reduced by carefully positioning suitable "obstacles" that direct the outflow (see for instance [START_REF] Helbing | Pedestrian, crowd and evacuation dynamics[END_REF] and the references therein). Indeed, looking at Figure 4, we observe that the time of exist with obstacles is slightly shorter than the one without any obstacle.

Several populations

A natural wish now is to extend the previous models to the case of several population with different objectives.

Panic

With several populations, we have to consider several densities and several equations. For two populations, extending the idea of the equation ( 5) to the case of several populations, we obtain the system

∂ t ρ 1 + Div ρ 1 v(ρ 1 * η 1 + ρ 2 * η 2 ) ν 1 (x) = 0 , ∂ t ρ 2 + Div ρ 2 v(ρ 1 * η 1 + ρ 2 * η 2 ) ν 2 (x) = 0 . ( 10
)
Here we consider that the speed v depends on the average of the total density ρ 1 +ρ 2 . Furthermore, the difference of goals of the populations 1 and 2 is reflected in the choice of different prefered directions ν 1 and ν 2 .

We are able to prove theorems similar to the ones of section 2.1.1. In collaboration with R. M. Colombo [START_REF] Colombo | Nonlocal crowd dynamics models for several populations[END_REF], we obtained: Theorem 2.5 (see [START_REF] Colombo | Nonlocal crowd dynamics models for several populations[END_REF]).

Let ρ 0 = (ρ 0,1 , ρ 0,2 ) ∈ (L 1 ∩ L ∞ ∩ BV)(R N , R 2 + ). Assume that, for i ∈ {1, 2}, v i ∈ (C 2 ∩ W 2,∞ )(R, R), ν i ∈ (C 2 ∩ W 2,1 )(R N , R N ), η i ∈ (C 2 ∩ W 2,∞ )(R N , R).
Then there exists a unique weak entropy solution ρ = S t ρ 0 ∈ C 0 (R + , L 1 (R N , R 2 + )) to ( 5) with initial condition ρ 0 . Furthermore we have the estimate

ρ(t) L ∞ ≤ ρ 0 L ∞ e Ct , (11)
where the constant C depends on v, ν and η Similarly as for one population, it is possible with several populations to prove the Gâteauxdifferentiability thanks to Kružkov theory (see [START_REF] Colombo | Nonlocal crowd dynamics models for several populations[END_REF]Theorem 2.2]).

As in Theorem 2.5, the hypotheses of Theorem 2.5 are very strong. In collaboration with G. Crippa [START_REF] Crippa | Existence and uniqueness of measure solutions for a system of continuity equations with non-local flow[END_REF], using tools from optimal transport theory, we obtained the better result: Theorem 2.6 (see [START_REF] Crippa | Existence and uniqueness of measure solutions for a system of continuity equations with non-local flow[END_REF]).

Let ρ 0 ∈ P(R N ) 2 . Assume v i ∈ (L ∞ ∩ Lip)(R, R), ν i ∈ (L ∞ ∩ Lip)(R N , R N ), η i ∈ (L ∞ ∩ Lip)(R N , R + ).
Then there exists a unique weak measure solution ρ ∈ L ∞ (R + , P(R N ) 2 ) to [START_REF] Daamen | Experimental research of pedestrian walking behavior[END_REF] with initial condition ρ 0 .

If furthermore

ρ 0 ∈ L 1 (R N , R + ) then for all t ≥ 0, we have ρ(t) ∈ L 1 (R N , R + ).
The idea of the proof for this theorem is given in Section 4.

Interaction continuum / individuals. Note that in the framework of Theorem 2.6, we are dealing with measure solutions. This context allows us to describe a coupling between a group of density ρ 1 and an individual located in p(t). Indeed, let us assume that

ρ 1 ∈ L ∞ (R + , L 1 (R N , R + )) and ρ 2 = δ p(t)
is a Dirac measure. Then we have the following ODE/PDE coupling:

     ∂ t ρ 1 + div ρ 1 v ρ 1 * η 1 (x) + η 2 (x -p(t)) ν 1 (x) = 0 , ṗ(t) = v ρ 1 * η 1 p(t) + η 2 (0) ν 2 (p(t)) .

Orderly crowd

We now extend [START_REF] Coscia | First-order macroscopic modelling of human crowd dynamics[END_REF] to the case of several populations. We consider here not only that ρ 1 and ρ 2 have two different goals, but also that ρ 1 is repelled by ρ 2 , and that ρ 2 is repelled by ρ 1 . Let us denote ε i the parameter of self-interaction and ε o the parameter of interaction with the other population. We obtain:

               ∂ t ρ 1 + Div ρ 1 v 1 (ρ 1 ) ν 1 (x) -ε i ∇ρ 1 * η 1 1 + ∇ρ 1 * η 1 2 -ε o ∇ρ 2 * η 2 1 + ∇ρ 2 * η 2 2 = 0 , ∂ t ρ 2 + Div ρ 2 v 2 (ρ 2 ) ν 2 (x) -ε o ∇ρ 1 * η 1 1 + ∇ρ 1 * η 1 2 -ε i ∇ρ 2 * η 2 1 + ∇ρ 2 * η 2 2 = 0 . ( 12
)
Remark 2.7. Note that the interaction between the equations is only in the nonlocal term.

We obtain again existence and uniqueness of weak entropy solutions thanks to Kružkov theory (see [START_REF] Colombo | Nonlocal crowd dynamics models for several populations[END_REF]Theorem 3.2]). The obtained theorem is similar to Theorem 2.4 so we don't rewrite it.

A classical situation considered in the engineering literature is that of two groups of people moving in opposite directions and crossing each other. The numerical integration (see Figure 5), shows formation of lanes that are not superimposing as described in the engineering literature ( [START_REF] Helbing | Self-organizing pedestrian movement[END_REF][START_REF] Daamen | Experimental research of pedestrian walking behavior[END_REF] or [START_REF] Appert-Rolland | Two-way multi-lane traffic model for pedestrians in corridors[END_REF] ). Figure 5: Numerical integration of ( 12). Above, the population ρ 1 moving to the right and, below, ρ 2 moving to the left. Note the lanes that are formed. First, due to the self-interaction (ε i = 0.3) within each populations and then due to the crossing between the two populations (ε o = 0.7). The latter lanes of different populations do not superimpose. Picture from [START_REF] Colombo | Nonlocal crowd dynamics models for several populations[END_REF].

Another situation developping interesting features is the following: two populations are initially uniformly distributed in the same region. At time t =0, the first populations starts moving towards an exit (first line of Figure 6), on the right while the second moves only to let the first one pass (second line of Figure 6). Figure 6: Numerical integration of ( 12) when with ε i = 0, ε o = 0.3. Above, the population ρ 1 and, below, ρ 2 . Note first the formation of small clusters separating the two populations, then lanes and, finally, a sort of fingering. Picture from [START_REF] Colombo | Nonlocal crowd dynamics models for several populations[END_REF].

Interaction group / isolated agent

Let us finally introduce a model that would take into account the situation in which an isolated individual interact with a crowd. We think for example to the situation in which a predator is running after a group of preys. We can also think to the case of a leader willing to carry a group of followers to a given region. Let ρ ∈ R + be the density of the group and p ∈ R k be the position of an isolated agent (e.g. a leader or a predator). We describe the interaction by the coupling (see [START_REF] Colombo | An analytical framework to describe the interactions between individuals and a continuum[END_REF]):

     ∂ t ρ + Div ρ V t, x, ρ, p(t) = 0 , ṗ = ϕ t, p, Aρ(t) (p(t)) , (t, x) ∈ R + × R N , (13)
with initial conditions ρ(0, x) = ρ 0 (x) , p(0) = p 0 .

Using once again Kružkov theory and tools on the stability of ordinary differential equations, in collaboration with R. M. Colombo, we proved existence and uniqueness of solutions (see [7, Theorem 2.2]).

For example, we can consider:

• Followers / Leader (see Figure 7): here the vector p ∈ R 2 is the position of the leader and ρ is the density of the group of followers. The function v(ρ) describes essentially the speed of the followers and is, as usual, a decreasing function, vanishing in ρ = 1; the direction of a follower located in x is given by the vector p(t) -x, directed toward the leader. The velocity of the leader is increasing with respect to the averaged density ρ * η, computed in the position of the leader. Indeed, we expect the leader to wait for the followers to join him when the density of followers is small around him, and to accelerate when the density of followers becomes bigger; this means low speed when the averaged density is small and high speed when the averaged density is maximal. The direction of the leader is a chosen vector field ψ(t).

   ∂ t ρ + Div ρ v(ρ) (p(t) -x)e -p-x = 0 , ṗ = (1 + ρ * η(p(t)) ) ψ(t) .
Solution at time 0.000 13) in a predator/prey interaction. Picture from [START_REF] Colombo | An analytical framework to describe the interactions between individuals and a continuum[END_REF].

• Sheeps / Dogs: in this example, there are n dogs, located in p i (t) ∈ R 2 for i ∈ {1, . . . , n} and a group of sheeps of density ρ(t, x). The dogs want to constrain the sheeps to a given area by running around them and the sheeps are afraid by the dogs. As above, the speed of the sheeps is given by the decreasing function v(ρ) that vanishes in ρ = 1.

The direction of a sheep located in x is a sum of two terms. The first one is their prefered direction ν(x) ; the second one is given by the vector n i=1 (x -p i (t))e -p i -x representing the repulsive effects of the dogs on the sheeps, and its presence is is less felt if it is too far away. Each dog run around the flock tacking the direction perpendicular to the direction of maximal averaged density.

           ∂ t ρ + Div ρ v(ρ) ν(x) + n i=1 (x -p i )e -p i -x = 0 , ṗi = (ρ * ∇η) ⊥ (p i (t)) 1 + ρ * ∇η(p i (t)) 2 ,
for all i ∈ {1, . . . n} .

• Preys / Predator: here a predator, located in p(t) ∈ R 3 is pursuing a group of preys of density ρ(t, x). For example, we can think to an hawk pursuing a swarm of little birds or to a shark pursuing a shoal of little fishes. As above, the speed of the preys is given by the decreasing function v(ρ) that vanishes in ρ = 1. The direction of a prey located in x is given by the vector (x -p(t))e -p(t)-x so that the predator has a repulsive effect on the preys, and its presence is is less felt if it is too far away. The acceleration of the predator is directed toward the maximal averaged density of preys, as felt from its position, so it is directed as ∇(ρ * η)(p(t)).

       ∂ t ρ + Div ρ v(ρ) 1 + e -x-p(t) x -p(t) = 0 , d 2 p dt 2 = (ρ * ∇)η(p(t)) .
Note that Theorem 2.6 deals with measure solutions. Consequently, if ρ 1 ∈ L 1 and ρ 2 = δ p(t) , we recover the above coupling PDE/ODE. However, in Theorem 2.6, there are restricitions on the flow as we only deal with continuity equations. In particular, with a wellchosen nonlinearity with respect to ρ in the first equation of [START_REF] Helbing | Dynamics of crowd disasters: An empirical study[END_REF], we are able to provide for (13) an a priori uniform bound on the L ∞ norm, which is not possible for (10): in the three examples above, if v(1) = 0, then 0 ≤ ρ 0 ≤ 1 implies 0 ≤ ρ(t) ≤ 1 for all t ≥ 0.

Using Kružkov theory

General Theory on scalar conservation laws

In this section, we consider classical scalar conservation laws.

∂ t u + Divf (t, x, u) = F (t, x, u) , u(0, x) = u 0 (x) , ( 14
)
where f is the flow and F is the source. Note that, Div stands for the total divergence whereas div stands for the partial divergence. Thus Divf (t, x, u(t, x)) = divf (t, x, u(t, x))

+ ∂ u f (t, x, u(t, x)) • ∇u(t, x).
Let us recall the definition of weak entropy solution:

Definition 3.1. The function ρ ∈ L ∞ ([0, T ] × R N ,
∈ C 0 c (] - ∞, T [×R N , R + ), R * + R N (u -k) ∂ t ϕ + f (t, x, u) -f (t, x, k) ∇ x ϕ + F (t, x, u) -divf (t, x, k) ϕ ×sign(u -k) dx dt + R N ρ 0 (x) -k ϕ(0, x) dx ≥ 0 . (15)
We also recall the well-known Kružkov theorem [17, Theorem 4], giving existence and uniqueness of weak entropy solutions.

Theorem 3.2 (Kružkov). Let u 0 ∈ (L ∞ ∩ L 1 )(R N , R). For all A, T ≥ 0, we denote Ω A T = [0, T ] × R N × [-A, A] and Ω = R + × R N × R. Under the conditions f ∈ C 0 (Ω; R N ), F ∈ C 0 (Ω; R) and (K)            f ∈ C 0 (Ω; R N ) , F ∈ C 0 (Ω; R)
, and for all T, A > 0 :

f , F have continuous derivatives: ∂ u f , ∂ u ∇f , ∇ 2 f , ∂ u F , ∇F , for all T, A > 0 , ∂ u f ∈ L ∞ (Ω A T ) , F -divf ∈ L ∞ (Ω A T ) , ∂ u (F -divf ) ∈ L ∞ (Ω A T )
there exists a unique weak entropy solution u ∈ L ∞ ([0, T ]; L 1 (R N ; R)) of ( 14) that is rightcontinuous in time.

Let v 0 ∈ (L 1 ∩ L ∞ )(R N ; R).
Let u be the solution associated to the initial condition u 0 and v be the solution associated to the initial condition v 0 . Let M be such that

M ≥ sup( u L ∞ (R + ×R N ;R) , v L ∞ ([0,T ]×R N ;R) ). Then, for all t ∈ [0, T ], with γ = ∂ u F L ∞ (Ω M T ) , we have (u -v)(t) L 1 ≤ e γt u 0 -v 0 L 1 . (16)

Estimate on the dependence with respect to flow and source

Let us now consider a pair of different flows f, g

∈ C 2 ([0, T ] × R N × R; R N ) and different sources F, G ∈ C 1 ([0, T ] × R N × R; R), let us denote u, v the solutions of ∂ t u + Div f (t, x, u) = F (t, x, u) , u(0, •) = u 0 , (17)
∂ t v + Div g(t, x, v) = G(t, x, v) , v(0, •) = v 0 , (18)
with initial conditions u 0 , v 0 ∈ (L 1 ∩L ∞ ∩BV)(R N , R). We want here to estimate

(u -v)(t) L 1 with f -g, F -G, u 0 -v 0 .
To do that, we use the doubling variable method due to Kružkov [17]. This strategy was already employed by Lucier [START_REF] Lucier | A moving mesh numerical method for hyperbolic conservation laws[END_REF] and Bouchut & Perthame [START_REF] Bouchut | Kružkov's estimates for scalar conservation laws revisited[END_REF] to study the case in which the flows f, g depend only on u and the source are identically zero (F = G = 0).

A key ingredient in the previous results is that the solution of

∂ t u + Div(f (u)) = 0 with initial condition u 0 ∈ (L ∞ ∩ L 1 ∩ BV)(R N , R) satisfies TV(u(t)) ≤ TV(u 0 ).
Let us remind the definition of total variation. Definition 3.3. For u ∈ L 1 loc (R N ; R) we denote the total variation of u:

TV(u) = sup R N u divΨ ; Ψ ∈ C 1 c (R N ; R N ) , Ψ L ∞ ≤ 1 .
The space of function with bounded variation is then defined as

BV(R N ; R) = u ∈ L 1 loc ; TV(u) < ∞ . Remark 3.4. If u ∈ (C 1 ∩ W 1,1 )(R N , R) then TV(u) = ∇u L 1 .
When f and F are not depending on u, we have

u 0 ∈ (L ∞ ∩ BV)(R N , R) ⇒ ∀t ≥ 0 , u(t) ∈ (L ∞ ∩ BV)(R N , R) and, with γ = ∂ u F L ∞ (Ω M T ) , TV(u(t)) ≤ TV(u 0 )e γt .
Theorem 3.5 (see [START_REF] Lucier | A moving mesh numerical method for hyperbolic conservation laws[END_REF]). Let f, g : R → R N be globally Lipschitz, let u 0 , v 0 ∈ (L 1 ∩L ∞ )(R N ; R) be initial conditions of

∂ t u + Div f (u) = 0 , ∂ t v + Div g(v) = 0 .
Assume furthermore that v 0 ∈ BV(R N ; R). Then for all t ≥ 0,

(u -v)(t) L 1 ≤ u 0 -v 0 L 1 + C t TV(v 0 ) Lip (f -g) ,
where C > 0 is a given constant.

We want to generalize this theorem. Let us introduce the set of hypotheses:

(FS)      for all U, T > 0 , ∂ u F ∈ L ∞ (Ω U T ) T 0 R N (F -div f )(t, x, •) L ∞ ([-U,U ];R) dx dt < +∞ .
Denoting

Ω V T = [0, T ] × R N × [-V, V ],
we obtain: Theorem 3.6 (see [START_REF] Lécureux-Mercier | Improved stability estimates on general scalar balance laws[END_REF]). Assume that (f, F ), (g, G) satisfy (K), that (f -g, F -G) satisfies (FS). Let u 0 , v 0 ∈ (L 1 ∩ L ∞ ∩ BV)(R N ; R). Let u and v be the solutions of ( 17) and ( 18) associated to (f, F ) and (g, G) with initial conditions u 0 and v 0 . Assume furthermore that TV(u(t)) < ∞ for all t ≥ 0.

Let V = max( u L ∞ , v L ∞ ) and κ = ∂ u F L ∞ (Ω V T )
. Then for all t ∈ [0, T ]:

(u -v)(t) L 1 ≤ e κt u 0 -v 0 L 1 + e κt sup τ ∈[0,t] TV(u(τ )) t 0 ∂ u (f -g)(τ ) L ∞ (R N ×[-V ;V ]) dτ + t 0 e κ(t-τ ) R N ((F -G) -div (f -g))(τ, x, •) L ∞ ([-V,V ]) dxdτ .
The proof of this theorem is based on the Kružkov doubling variables method and is detailed in [START_REF] Lécureux-Mercier | Improved stability estimates on general scalar balance laws[END_REF].

In some particular cases, we recover known estimates:

• f (u), g(u), F = G = 0 : κ = 0 and (u -v)(t) L 1 ≤ u 0 -v 0 L 1 + t TV(u 0 ) ∂ u (f -g) L ∞ (Ω V T ) . • f (t, x), F (t, x) : κ = 0, u(t, x) = u 0 (x) + t 0 (F -divf )(τ, x) dτ , and
(u -v)(t) L 1 ≤ u 0 -v 0 L 1 + t 0 R N ((F -G) -div (f -g))(τ, x) dxdτ .

Total variation estimate

To complete Theorem 3.6, we need an estimate on the total variation, defined in Definition 3.3.

To do so we need a few more hypotheses. Indeed, there is no reason why the total variation should be bounded for all time. In fact, let us consider the equation ∂ t u+ ∂ x cos x = 0, u 0 = 0. The solution of this Cauchy problem is u(t, x) = t sin x whose total variation is 0 at time 0 and +∞ for any time t > 0.

Our goal now is to give a general estimate on the total variation. Let us introduce the following set of hypotheses (TV) :

           for all A, T > 0 ∇∂ u f ∈ L ∞ (Ω A T ; R N ×N ) , ∂ u F ∈ L ∞ (Ω A T ; R) , T 0 R N ∇(F -divf )(t, x, •) L ∞ ([-A,A];R N ) < ∞ .
We obtain

Theorem 3.7. Let us assume (f, F ) satisfies (K)-(TV). Denote W N = π/2 0 (cos θ) N dθ, M = u L ∞ ([0,T ]×R N )
, and

κ 0 = (2N + 1) ∇∂ u f L ∞ (Ω M T ) + ∂ u F L ∞ (Ω M T ) .
Then, the weak entropy solution u of ( 14) satisfies u(t) ∈ BV(R N ; R) for all t > 0, and

TV u(T ) ≤ TV(u 0 ) e κ 0 T + N W N T 0 e κ 0 (T -t) R N ∇(F -divf )(t, x, •) L ∞ ([-Ut,Ut];R) dx dt ,
where U t = sup y∈R N u(t, y) .

The proof of this theorem is based on a good representation formula for the total variation and on the doubling variable method and is detailed in [START_REF] Lécureux-Mercier | Improved stability estimates on general scalar balance laws[END_REF]. Remark 3.8. In some cases, we recover known estimates.

• When f depends only on u and F = 0, we have a result similar to the one that was already known: TV(u(t)) ≤ TV(u 0 ).

• When f and F do not depend on u, the equation reduces in fact to the ODE ∂ t u = (F -divf )(t, x), whose solution writes

u(t, x) = u 0 (x) + t 0 (F -divf )(τ, x) dτ .
Meanwhile, the bound above reduces to

TV(u(t)) ≤ TV(u 0 ) + N W N t 0 R N (F -divf )(τ, x) dτ
which is essentially what we expected, up to the coefficient N W N .

Proof. Let us admit the following proposition:

Proposition 3.9. Let µ ∈ C ∞ c (R + ; R + ) such that µ L 1 = 1 and µ ′ < 0 on R * + . Let u ∈ L 1 loc (R N , R). We define µ λ (x) = 1 λ N µ x λ .
If there exists C 0 > 0 such that ∀λ > 0,

I(λ) = 1 λ R N R N u(x + y) -u(x) µ λ (y)dxdy ≤ C 0 , then u ∈ BV(R N , R) and C 1 TV(u) = lim λ→0 I(λ) ≤ C 0 .
with

C 1 = R N |y 1 |µ( y )dy.
Hence, we introduce:

F(T, λ) = T 0 R N B(x 0 ,R+M (T 0 -t))
u(x + y) -u(x) µ λ (y)dx dy dt , our goal being to estimate this quantity using the Kružkov doubling variable method. Doubling variable method. Let us denote u = u(t, x) and v = u(s, y) for (t, x), (s, y) ∈

R * + × R N . For all k, l ∈ R and all test-function ϕ = ϕ(t, x, s, y) ∈ C 1 c (R * + × R N ) 2 ; R + ,
R * + R N (u -k) ∂ t ϕ + f (t, x, u) -f (t, x, k) ∇ x ϕ + F (t, x, u) -divf (t, x, k) ϕ ×sign(u -k) dx dt ≥ 0 (19)
and

R * + R N (v -l) ∂ s ϕ + f (s, y, v) -f (s, y, l) ∇ y ϕ + (F (s, y, v) -divf (s, y, l))ϕ ×sign(v -l) dy ds ≥ 0 . (20)
Let us choose k = v(s, y) in [START_REF] Lucier | A moving mesh numerical method for hyperbolic conservation laws[END_REF] and integrate in (s, y). Similarly, we choose l = u(t, x) in [START_REF] Maury | A macroscopic crowd motion model of the gradient-flow type[END_REF] and we integrate in (t, x). Furthermore, we choose ϕ(t, x, s, y) = Ψ(t -s, x -y)Φ(t, x) and we summate to get

R * + R N R * + R N sign(u -v) (u -v) Ψ ∂ t Φ + f (t, x, u) -f (t, x, v) • (∇Φ) Ψ + f (s, y, v) -f (s, y, u) -f (t, x, v) + f (t, x, u) • (∇Ψ) Φ + F (t, x, u) -F (s, y, v) + divf (s, y, u) -divf (t, x, v) ϕ dx dt dy ds ≥ 0 . (21)
We choose now µ and ν unit approximations and

Ψ(t -s, x -y) = 1 λ N µ x -y λ 1 η ν t -s η ,
so that Ψ converges to a Dirac in space and time when λ, η → 0; and we choose Φ such that Φ → 1 [0,T ]×B(x 0 ,R+M (T -t)) when θ, ε → 0. The last step consists in estimating the integral (21) when η, θ, ε → 0, keeping the λ fixed to make appear F.

Estimate on F. By the doubling variable method, we obtain

∂ T F(T, λ) ≤∂ T F(0, λ) + Cλ ∂ λ F(T, λ) + N λ F(T, λ) + C ′ F(T, λ) + λ T 0 A(t)dt , (22)
where

A(t) = M 1 R N ∇(F -div f )(t, x, •) L ∞ (u) dx , M 1 = R N y µ( y )dy, C ′ = N ∇∂ u f L ∞ + ∂ u F L ∞ . C = ∇∂ u f L ∞ .
We use the hypothesis

u 0 ∈ BV(R N , R) to get 1 λ ∂ T F(0, λ) ≤ M 1 TV(u 0
), and we integrate with respect to time, obtaining:

0 ≤ M 1 C TV(u 0 ) + ∂ λ F(T, λ) + α(T ) λ F(T, λ) + 1 C T 0 A(t)dt ,
where

α(T ) = N + C ′ C - 1 C T → T →0 -∞ .
Next, we choose T so that α < -1 and we integrate on [λ, +∞[. We obtain:

F(T, λ) ≤ λ -α -1 M 1 C TV(u 0 ) + λ C(-α -1) T 0 A(t)dt .
Besides, note that, thanks to the shape of µ, we can find K > 0 so that

∂ λ F(T, λ)+ N λ F = T 0 R N B(x 0 ,R+M (T 0 -t)) u(x + y) -u(x) µ ′ λ (y)dxdydt ≤ K 2λ F(T, 2λ) .
Hence, we finally obtain

1 λ ∂ T F(T, λ) ≤ 1 λ ∂ T F(0, λ) + CK + C ′ (-α -1)C M 1 TV(u 0 ) + T 0 A(t)dt + T 0 A(t)dt .
By Proposition 3.9 we have then u(t) ∈ BV(R N , R).We can now divide [START_REF] Piccoli | Time-evolving measures and macroscopic modeling of pedestrian flow[END_REF] by λ and make λ → 0. Applying a standard Gronwall-type argument, we obtain the desired inequality.

Scalar conservation law with a non-local flow

Let us consider equation ( 5) with initial condition

ρ 0 ∈ (L 1 ∩ L ∞ ∩ BV)(R N , R).
We want to prove Theorem 2.1 and Theorem 2.3. Note that this proof adapt smoothly to prove Theorem 2.4 and 2.5. We use here the scheme proposed in Section 1, choosing X = C 0 ([0, T ], L 1 (R N , R)) at point (a). Point (b) is then satisfied by Kružkov theorem (see Theorem 3.2) and point (c) is satisfied by the stability estimate given by Theorem 3.6 & 3.7. This scheme of proof can be also used for the other models: the proof of Theorem 2.4 for an orderly crowd is based on it, as well as Theorem 2.5 giving existence and uniqueness for the several populations models [START_REF] Daamen | Experimental research of pedestrian walking behavior[END_REF] and [START_REF] Helbing | Pedestrian, crowd and evacuation dynamics[END_REF]. For more details on these various theorems, refer to [START_REF] Colombo | Nonlocal crowd dynamics models for several populations[END_REF][START_REF] Colombo | A class of non-local models for pedestrian traffic[END_REF].

We use also this scheme of proof in order to prove existence and uniqueness of solution to equation [START_REF] Helbing | Dynamics of crowd disasters: An empirical study[END_REF], where we need furthermore stability results on the ODE (see [START_REF] Colombo | An analytical framework to describe the interactions between individuals and a continuum[END_REF]).

Proof of Theorem 2.1

Let us introduce the space

X = C 0 ([0, T ], L 1 (R N ; R + )) .
Let r ∈ X and ρ ∈ X be the solution of (3) given by Kružkov Theorem (see Theorem 3.2). Note that, u = 0 being a solution to (3), ρ 0 ≥ 0 implies ρ(t) ≥ 0 for all t ≥ 0 thanks to the maximal principle on scalar conservation laws [START_REF] Kružkov | First order quasilinear equations with several independent variables[END_REF]Theorem 3]. Remark 3.10. For the model ( 8) of orderly crowd, we can use the same king of argument to have a uniform bound on the L ∞ norm. Indeed, if u = 1 is a solution to the equation

∂ t ρ + Div(ρv(ρ) w(x)) = 0 (23)
then ρ 0 ≤ 1 implies ρ(t) ≤ 1 for all t ≥ 0. This fact allows us in this case to consider the space

X = C 0 ([0, T ], L 1 (R N , [0, 1])).
The same fact is also true for each density of the model ( 12) of orderly crowd with several populations. Indeed, the interaction between the equations is only in the non-local term; so, if we fix the nonlocal term, we have to deal with the same kind of scalar conservation law as [START_REF] Villani | Topics in optimal transportation[END_REF]. In particular, for each density ρ i the maximum principle gives us

0 ≤ ρ 0,i ≤ 1 ⇒ 0 ≤ ρ i (t) ≤ 1 for all t ≥ 0 .
So finally, if there are k populations, we have the following uniform bound on the total density density

0 ≤ k i=1 ρ i (t) ≤ k .
Let us introduce the application

Q : r ∈ X → ρ ∈ X .
For r 1 , r 2 , we get by Theorem 3.6 and Theorem 3.7,

Q(r 1 ) -Q(r 2 ) L ∞ ([0,T ],L 1 ) ≤ f (T ) r 1 -r 2 L ∞ ([0,T ],L 1 ) ,
where f is continuous, increasing, f (0) = 0 and f → T →∞ ∞. Consequently, we can choose T small enough so that f (T ) = 1/2 and apply the Banach fixed point theorem for a small time.

Iterating the process in time, we obtain global in time existence.

Gâteaux derivative of the semi-group

We prove now that the semi-group obtained in Theorem 2.1 (or 2.5) is Gâteaux-differentiable with respect to the initial condition. The following result is obtained by the use of Kružkov theory.

Note that the Gâteaux-differentiability is also possible for the more non-linear models of orderly crowd [START_REF] Coscia | First-order macroscopic modelling of human crowd dynamics[END_REF] and [START_REF] Helbing | Pedestrian, crowd and evacuation dynamics[END_REF]. Indeed an important tool below is that the equation ( 5) preserves the regularity of the initial condition, which is no longer true for (8). Definition 3.11. We say the application S :

L 1 (R N ; R) → L 1 (R N ; R) is Gâteaux differen- tiable in L 1 , in ρ 0 ∈ L 1 , in the direction r 0 ∈ L 1 if there exists a linear continuous application DS(ρ 0 ) : L 1 → L 1 such that S(ρ 0 + hr 0 ) -S(ρ 0 ) h -DS(ρ 0 )(r 0 ) L 1 → h→0 0 .
We expect the Gâteaux derivative to be the solution of the linearized problem [START_REF] Colombo | An analytical framework to describe the interactions between individuals and a continuum[END_REF]. First, we prove that the linearized problem admits an entropic solution (see [START_REF] Colombo | Control of the continuity equation with a non local flow[END_REF][START_REF] Colombo | Nonlocal crowd dynamics models for several populations[END_REF]).

Theorem 3.12.

Assume that v ∈ (C 4 ∩ W 2,∞ )(R, R), ν ∈ (C 3 ∩ W 2,1 ∩ W 2,∞ )(R N , R N ), η ∈ (C 3 ∩W 2,∞ )(R N , R + ). Let ρ ∈ C 0 ([0, T ex [; W 1,∞ ∩W 1,1 (R N , R)), r 0 ∈ (L 1 ∩L ∞ )(R N ; R).
Then the linearized problem [START_REF] Colombo | An analytical framework to describe the interactions between individuals and a continuum[END_REF] with initial condition r 0 admits a unique entropic solution

r ∈ C 0 ([0, T ex [; L 1 (R N ; R)) and we denote Σ ρ t r 0 = r(t, •). If furthermore r 0 ∈ W 1,1 (R N , R), then for all t ∈ [0, T ex [, r(t) ∈ W 1,1 (R N , R).
The proof is similar to the one of Theorem 2.1 so we omit it.

Theorem 3.13. Assume that v ∈ (C 4 ∩ W 2,∞ )(R, R), ν ∈ (C 3 ∩ W 2,1 ∩ W 2,∞ )(R N , R N ), η ∈ (C 3 ∩ W 2,∞ )(R N , R + ). Let ρ 0 ∈ (W 1,∞ ∩ W 2,1 )(R N , R + ), r 0 ∈ (W 1,1 ∩ L ∞ )(R N , R).
Then, for every t ≥ 0, the local semigroup of the pedestrian traffic problem is L 1 Gâteaux differentiable in the direction r 0 and

DS t (ρ 0 )(r 0 ) = Σ Stρ 0 t r 0 .
We prove this theorem following the same sketch of proof as above. Proof. Let ρ, ρ h the solutions of the problem (5) with initial condition ρ 0 , ρ 0 + hr 0 .

Let r be the solution of the linearized problem (7), r(0) = r 0 and let z h = ρ + hr which then satisifies

∂ t z h +Div z h (v(ρ * η) + hv ′ (ρ * η)(r * η)) ν(x) = h 2 Div(rv ′ (ρ * η)(r * η) ν(x)) , z h (0) = ρ 0 +hr 0 .
Then we use Theorem 3.6 & 3.7 to compare ρ h et z h . We get:

1 h ρ h -z h L ∞ ([0,T [,L 1 ) ≤F (T ) 1 h ρ h -ρ 2 L ∞ (L 1 ) + 1 h ρ h -z h L ∞ (L 1 ) + hC(β)T e C(β)T r L ∞ (W 1,1 ) r L ∞ (L 1 ) ,
with F increasing and F (0) = 0. Then, we choose T small enough so that F (T ) = 1/2. Besides, applying Theorem 3.6 to compare ρ and ρ h and using Gronwall lemma, we obtain that 1 h ρ h -ρ L ∞ (L 1 ) remains bounded when h → 0. Finally, we can take the limit h → 0 and obtain the desired result.

Extrema of a cost functional

We are now able to characterize the extrema of a given cost functional. Let J be a cost functional such that

J(ρ 0 ) = R N f (S t ρ 0 ) ψ(t, x)dx . Proposition 3.14. Let f ∈ C 1,1 (R; R + ) and ψ ∈ L ∞ (R + × R N ; R). Assume S : [0, T ] × (L 1 ∩ L ∞ )(R N ; R) → (L 1 ∩ L ∞ )(R N ; R) is L 1 Gâteaux differentiable. If ρ 0 ∈ (L 1 ∩ L ∞ )(R N ; R) is solution of the problem: Find min ρ 0 J(ρ 0 ) such that S t ρ 0 is solution of problem (Traffic) . then, for all r 0 ∈ (L 1 ∩ L ∞ )(R N ; R) R N f ′ (S t ρ 0 ) Σ ρ 0 t r 0 ψ(t, x) dx = 0 .

Using the Optimal transport theory

In this section, we want to prove Theorem 2.2, giving existence and uniqueness of weak measure solurtion to ( 5) and [START_REF] Daamen | Experimental research of pedestrian walking behavior[END_REF]. In a more general settings, we consider the system:

∂ t ρ i + div(ρ i V i (x, ρ 1 * η i,1 , . . . , ρ k * η i,k )) = 0 , (t, x) ∈R + × R N (24)
ρ i (0) = ρi , i ∈{1, . . . , k} .

in which the coupling between the equations is only present through the nonlocal term. Fixing the nonlocal term, we obtain a system of decoupled continuity equations:

     ∂ t ρ 1 + div(ρ 1 b 1 (t, x)) = 0 , . . . ∂ t ρ k + div(ρ k b k (t, x)) = 0 ,
where b 1 , . . . , b k are regular with respect to x.

Existence and uniqueness of weak measure solutions

Let us remind the following definitions of weak measure solution, but also of Lagrangian solution.

Definition 4.1. ρ ∈ L ∞ ([0, T ], M + (R N ) k ) is a measure solution of (24) if, ∀ϕ ∈ C ∞ c (] - ∞, T ] × R N , R) T 0 R N ∂ t ϕ + V i (x, ρ * η i ) • ∇ϕ dρ i t (x) dt + R N ϕ(0, x) dρ i (x) = 0 . Definition 4.2. ρ ∈ L ∞ ([0, T ], M + (R N ) k ) is a Lagrangian solution with initial condition ρ ∈ M + (R N ) k if there exists an ODE flow X i : [0, T ] × R N → R N , solution of    dX i dt (t, x) = V i (X i (t, x), ρ t * η i (X i (t, x))) , X i (0, x) = x ;
and such that ρ

i t = X i t ♯
ρi where X i t : R N → R N is the map defined as

X i t (x) = X i (t, x) for any (t, x) ∈ R + × R N . Definition 4.3. Let µ be measure on Ω and T : Ω → Ω ′ a measurable map. Then T ♯ µ is the push-forward of µ if for any ϕ ∈ C 0 c (Ω ′ ), Ω ′ ϕ(x) dT ♯ µ(x) = Ω ϕ T (y) dµ(y) .
If T ♯ dµ (x) = f (x) dx and dµ (y) = g(y) dy; and T is a C 1 -diffeomorphism. Then we have the change of variable formula g(x) = f (T (x)) det(∇T (x)) .

Instead of proving directly Theorem 2.2 (or Theorem 2.6), we prove first the existence and uniqueness of lagrangian solutions.

4.4. Let ρ ∈ M + (R N , R k ). Let us assume that V ∈ (L ∞ ∩ Lip)(R N × R k , R N ×k )
and that η ∈ (L ∞ ∩ Lip)(R N , R k×k ). Then there exists a unique Lagrangian solution to system (24) with initial condition ρ. Remark 4.5. Note that a Lagrangian solution is also a weak measure solution.

Besides, we recover the results of Theorem 2.1, except the continuity in time. Indeed, Lagrangian solutions conserve the regularity of the initial condition along the time: if ρ ∈ L 1 (R N , (R + ) k ) then ρ ∈ L ∞ (R + , L 1 (R N , (R + ) k )) and we have

ρ(t) L 1 = ρ L 1 .
If moreover ρ ∈ (L 1 ∩ L ∞ )(R N , (R + ) k ), then ρ(t) ∈ L ∞ for all t ≥ 0 and we have the estimate

ρ(t) L ∞ ≤ ρ L ∞ e Ct ,
where C depends on ρ M , V and η.

We conclude the proof of Theorem 2.2 stating the uniqueness of measure solutions, so that the Lagrangian solution obtained is also the unique measure solutions. Proposition 4.6. Under the same set of hypotheses as Theorem 4.4, the weak measure solutions for (24) are unique.

Proof of Theorem 4.4

Once again, the proof is based on the scheme described in Section 1. We consider below probability measures instead ; for any i, we could also consider positive measures of fixed total mass ρi We equip this space with the distance d(µ, ν) = sup t∈[0,T ] W 1 (µ t , ν t ). The Wasserstein distance W 1 is defined as follow: Definition 4.7. Let µ, ν ∈ P(R N ). Let us denote P x : R d ×R d → R d the projection on the first coordinate; that is, for any (u, v) ∈ R d × R d , P x (u, v) = u. In a similar way, P y : R d × R d → R d is the projection on the second coordinate; that is, for any (u, v) ∈ R d × R d , P y (u, v) = v. We denote Ξ (µ, ν) the set of plans, that is Ξ(µ, ν) = γ ∈ P(R N × R N ) : P x♯ γ = µ and P y ♯ γ = ν .

The Wasserstein distance of order one between µ and ν is Let ρ = (ρ 1 , . . . , ρ k ), σ = (σ 1 , . . . , σ k ) ∈ P(R N ) k . The Wasserstein distance of order one between ρ and σ, denoted W 1 (ρ, σ), as

W 1 (ρ, σ) = k i=1 W 1 (ρ i , σ i ) .
Note besides that, thanks to Proposition 4. where C = Lip x (V ) + Lip r (V ) r t M Lip(η), C ′ = Lip r (V )Lip(η) ρ M . We conclude taking γ 0 in an optimal way so that Q(0) = W 1 (ρ 0 , σ 0 ) and using the inequality W 1 (ρ t , σ t ) ≤ Q(t) .

The stability estimate allows us to apply Banach fixed point Theorem for T small enough.

Measure solutions are Lagrangian solutions

Proof of Proposition 4.6. Let ρ be a measure solution of (24). Let b = V (x, ρ * η) and denote σ the Lagrangian solution associated to ∂ t σ + div(σb) = 0 with σ(0) = ρ.

Then δ = ρ -σ is a measure solution of ∂ t δ + div(δb) = 0, with δ(0) = 0. That is to say, for any

ϕ ∈ C ∞ c (] -∞, T ] × R N , R), T 0 R N ∂ t ϕ + b i (t, x) • ∇ϕ dδ t dt = 0 .
Let ψ ∈ C 0 c (] -∞, T ] × R N , R). We can find ϕ ∈ C 1 c (] -∞, T ] × R N , R) so that ψ = ∂ t ϕ + b i (t, x) • ∇ϕ. Hence, for any ψ ∈ C 0 c (] -∞, T ] × R N , R), we have T 0 R N ψ dδ t dt = 0, which implies δ ≡ 0 a.e. so ρ = σ a.e.. Consequently, we have b i (t, x) = V i (x, σ * η i ), and ρ is a Lagrangian solution of (24).

Conclusion

In the Kružkov framework, we are able to prove existence and uniqueness of weak entropy solution for the equation ∂ t ρ + div(ρV (x, ρ, ρ * η)) = 0. Furthermore, we can prove uniform bound in L ∞ if V = v(ρ) W (x, ρ * η), with v(1 =)0. However, the required hypotheses are very strong: we need indeed V ∈ C 2 ∩ W 2,1 ∩ W 2,∞ .

In the optimal transport theory framework, we can treat only equations such that ∂ t ρ + div(ρV (x, ρ * η)) = 0 . For this equation, we have only L ∞ bound that are exponentially growing in time. The hypotheses are nevertheless weaker since we only ask V ∈ Lip ∩ L ∞ , but we are no longer able to prove the Gâteaux-differentiability of the semi-group.

Figure 1 :

 1 Figure 1: Solution to (8) at times t = 0, 2.529, 5.043, 7.557, 10.071, 15.014. First 3 lanes are formed, then the middle lane bifurcates forming the fourth lane. Picture from [4].

Figure 2 :

 2 Figure 2: Initial datum and room geometry. Picture from [4].

Figure 3 :

 3 Figure 3: Solution to (8) with different geometries, computed at time t = 2.521, 5.043 and 7.563. Picture from [4].

Figure 4 :

 4 Figure 4: Solution to (8) with ε = 0.2, at times t = 4.438, 6.253, 11.396. On the first line, no obstacle is present. On the second line, 4 columns direct the crowd flow. The evacuation time in the latter case is shorter than in the former one. Picture from[START_REF] Colombo | A class of non-local models for pedestrian traffic[END_REF].

Figure 7 :

 7 Figure 7: Numerical solution of (13) in a predator/prey interaction. Picture from[START_REF] Colombo | An analytical framework to describe the interactions between individuals and a continuum[END_REF].

  we have, by definition of weak entropy solution (seeDefinition 3.1)

Theorem

  Let us first introduce a spaceX = L ∞ ([0, T ], P(R N ) k ) .

W 1

 1 (µ, ν) = inf γ∈ Ξ (µ,ν) R N ×R N |x -y| dγ(x, y) .

  R) is called weak entropy solution to the Cauchy problem[START_REF] Helbing | Self-organizing pedestrian movement[END_REF] if it satisfies for all k ∈ R and any test-function ϕ

  8, we have Ct + tC ′ e Ct sup

	(r i t -s i t ) * η(z) =	R N	η(z -ζ)(dr i t (ζ) -ds i t (ζ)) ≤ Lip(η)W 1 (r i t , s i t ) .
	Integrating, we get
	Q(t) ≤ Q(0)e

τ W 1 (r t , s t ) .

Let us recall the following duality formula: Proposition 4.8 (cf. Villani [23, p. 207]). Let f, g be two probability measures. The Wasserstein distance of order one between f and g satisfies

Let us consider the equation

Let ρ i = X t♯ ρ be the Lagrangian solution of (25). Then the application

is well-defined.

(c) Stability estimate.

Proposition 4.9. Let ρ, σ ∈ P(R N ) and r, s

If ρ and σ are Lagrangian solutions of

We have the estimate:

where

Proof. Let X, Y be the ODE flows associated to ρ, σ. Let γ 0 ∈ Ξ(ρ, σ). Define

Then Q is a Lipschitz function and

By triangular inequality, we obtain Q ′ (t) ≤(Lip x (V ) + Lip r (V )Lip(r t * η))Q(t) + Lip r (V )

R N ×R N (r t -s t ) * η(Y t (y)) dγ 0 (x, y) .