M Flesch
email: m.flesch@vermon.com

M Pernot

J Provost

G Ferin

A Nguyen-Dinh

M Tanter

Thomas Deffieux

4D in-vivo Ultrafast Ultrasound Imaging using a Row-Column Addressed matrix and Coherently-Compounded Orthogonal Plane Waves

come   L'archive ouverte pluridisciplinaire

Introduction

The concept of ultrafast ultrasound imaging with plane waves has been introduced more than 15 years ago in order to map, at several thousand of images per second, the rapid transient events occurring in the human body [START_REF] Tanter | Ultrafast imaging in biomedical ultrasound[END_REF]. Recent progress of ultrafast imaging has allowed the development of novel modalities such as shear wave elastography [START_REF] Bercoff | ShearWaveTM elastography: A new real time imaging mode for assessing quantitatively soft tissue viscoelasticity[END_REF], ultrasensitive Doppler imaging and functional ultrasound imaging [START_REF] Mace | Functional ultrasound imaging of the brain: Theory and basic principles[END_REF], to name a few.

Recently, 4D ultrafast imaging was developed and experimentally applied to ultrafast Doppler and Shear-wave imaging in vivo in several organs [START_REF] Provost | 3D ultrafast ultrasound imaging in vivo[END_REF] [START_REF] Gennisson | 4-D ultrafast shear-wave imaging[END_REF] [START_REF] Provost | 3-D ultrafast doppler imaging applied to the noninvasive mapping of blood vessels in Vivo[END_REF], including the heart and the breast, using a 4D ultrafast ultrasound imaging research platform with 1024 active channels and a 32 x 32 matrix probe [START_REF] Ratsimandresy | A 3 MHz two dimensional array based on piezocomposite for medical imaging[END_REF] in contrast with the 128 channels used typically on regular commercial 2D scanner. Electronics with such a large channel count provide a high flexibility and as such are ideal for research applications [START_REF] Jensen | SARUS: A synthetic aperture real-time ultrasound system[END_REF]), but their complexities limit clinical and commercial implementations.

Large fully sampled 2D matrix arrays transducers are also difficult to manufacture because of the high density of elements to individually connect, which require complex and expensive manufacturing processes, especially at high frequency where the inter-element pitch is smaller. Moreover, the high electrical impedance and low sensitivity of the small elements compromise matrix array performances. Finally, the large coaxial cable count required to connect all individual elements to electronic channels strongly limits the ergonomics and the portability of the probe. On the scanner side, synchronizing electronic boards with more than 1024 independent channels remains prohibitive in terms of development and production costs. Scaling to an even higher number of elements such as 64 x 64 or 128 x 128 is currently out of reach outside of dedicated research platforms.

Consequently, the road for 4D ultrafast imaging toward the clinic may still rely on alternative methods to efficiently reduce the thousands individual probe elements into a more manageable number of system channels (regularly 128 channels), while still maintaining high spatial and temporal resolution.

Alternative methods have been extensively developed for volumetric ultrasound imaging to reduce the number of electronic channels, including mechanical scanning, 2D matrix with integrated analog frontend electronics such as multiplexers [START_REF] Hara | A New 80V 32x32ch Low Loss Multiplexer LSI for a 3D Ultrasound Imaging System[END_REF] or µ-beamformers [START_REF] Savord | Fully sampled matrix transducer for real time 3D ultrasonic imaging[END_REF]. However, these approaches provide relatively low volume rate and are not compatible with ultrafast imaging.

Sparse arrays have been proposed as a solution to maintain a reasonable amount of output channels. Sparse arrays are made of a limited number of elements that are generally randomly positioned over the array footprint [START_REF] Davidsen | Sparse Geometries for Two-Dimensional Array Transducers in Volumetric Imaging[END_REF]. However, this geometry with large spacing between elements is not straightforward for the transmission of plane waves or diverging waves.

In this paper, we study another configuration of under-sampled arrays, i.e., the RCA matrix, originally proposed in 2003 [START_REF] Morton | Theoretical assessment of a crossed electrode 2-D array for 3-D imaging[END_REF]. RCA matrices are composed of two sets of orthogonal vertical and horizontal overlapping transducers which effectively reduces the number of elements from N x N to N + N individual channels, as shown in figure 1, and ensure compatibility with commercial scanners at the cost of a reduced image quality. Compared to sparse arrays, RCA matrices have two advantages: a much larger active surface and a geometry that is inherently adapted to plane wave transmissions, which make them a good candidate for ultrafast plane wave imaging using coherent compounding.

Recently, the fabrication and assessment of RCA matrices using novel approaches based on various technologies such as lead zirconate titanate (PZT) [START_REF] Seo | A 256 x 256 2-D array transducer with row-column addressing for 3-D rectilinear imaging[END_REF], polyvinylidene fluoride layer (PVdF) [START_REF] Yen | A dual-layer transducer array for 3-D rectilinear imaging[END_REF], or capacitive micromachined ultrasonic transducers (cMUT) [START_REF] Logan | 2-D CMUT wafer bonded imaging arrays with a rowcolumn addressing scheme[END_REF]) have been the object of active research.

Moreover, RCA matrices have been leveraged in various imaging sequences such as conventional beam focusing [START_REF] Morton | Theoretical assessment of a crossed electrode 2-D array for 3-D imaging[END_REF], synthetic aperture techniques [START_REF] Rasmussen | 3-D imaging using row-column-addressed arrays with integrated apodization -Part i: Apodization design and line element beamforming[END_REF], vector Doppler imaging [START_REF] Lehrmann | 3-D Vector Velocity Estimation with Row-Column Addressed Arrays[END_REF], or plane wave emissions with multiplexed reception [START_REF] Chen | A column-row-parallel ultrasound imaging architecture for 3d plane-wave imaging and Tx 2nd-order harmonic distortion (HD2) reduction[END_REF]), but, to our knowledge, no in vivo study on the full volumetric Doppler imaging capability of the RCA scheme has been conducted or compared against a 2D matrix.

In this work, we introduce an OPW compounding strategy dedicated to RCA matrices and have evaluated its performance for 4D ultrafast Doppler imaging against a 2D matrix approach. We first evaluate OPW in terms of PSF in simulation and then proceed to demonstrate its feasibility in phantoms and in a healthy volunteer.

Materials & Methods

OPW compounding

Ultrafast ultrasound imaging using a 1D or 2D array relies both on the parallel receive of channels data and the decomposition of the transmit focusing on an efficient and versatile basis with coherent compounding. A typical emission basis is a set of tilted plane waves or diverging waves, which allows a good repartition of the propagated energy in the medium and a good decorrelation between emissions. The method benefits from dynamic focusing both in transmit and receive and thus converges to quality of synthetic aperture imaging. The key point of this approach is the flexibility between image quality and frame rate, which can be adjusted by varying the number of emissions compounded for each image.

In the case of a RCA matrix, the focusing in transmit or receive can only be performed using a cylindrical law due to the geometry of the array. To achieve spatial separation, it is thus necessary to use complementary transmit and receive focusing: typically a cylindrical transmit focusing with the rows and an orthogonal cylindrical receive focusing with the columns or vice versa. However, it is still possible to decompose the cylindrical transmit focusing on a more efficient basis such as tilted plane waves and to use parallel receive with conventional delay and sum beamforming for dynamic receive focusing. Regarding the transmit, the emissions basis has to be adapted to the constraints of the RCA matrix which obviously rules out spherically diverging or converging waves or plane waves with random directions. Herein, we propose to use horizontally and vertically tilted plane waves, which can easily be transmitted using the RCA array by applying transmit delays rows or columns.

Contrary to the synthetic aperture approach, since the transmit focusing is not perfect but an approximation of the cylindrical focusing, the symmetry between emission and reception is broken and transmitting with rows then receiving with columns or transmitting with columns then receiving with rows lead to different results. As such, in order to enforce isotropy of the PSF at least in the two main axis we propose to coherently combine both sequences in a scheme we denote HV+VH where HV stands for Horizontal emission with Vertical reception and VH for Vertical emission with Horizontal reception as shown in figure 2.

Such dual scheme yields to a number of Na emissions where Na is the total number of angles of plane waves for both directions. While a lower image quality is to be expected when compared against focused or synthetic aperture approaches, it has the benefit of allowing for high volume rates. Table 1 summarizes the number of emissions required by the aforementioned approaches. Where N is the number of elements along a dimension of the array and Na is the number of plane or diverging waves used.

In figure 2, the plane wave compounding strategies for 2D matrix and RCA matrix are illustrated with their respective PSF. Figure 2a shows the regular compounding strategy dedicated to ultrafast ultrasound imaging using 2D matrices that is based on our previous work [START_REF] Provost | 3D ultrafast ultrasound imaging in vivo[END_REF]. Figure 2b shows the proposed OPW strategy dedicated to ultrafast imaging using RCA matrices.

A c c e p t e d M a n u s c r i p t The HV PSF and VH PSF have an extended shape in the transmit focusing direction. This is mainly caused by the partial focusing in transmit relying on the compounding of a few plane waves. As noted earlier, to reduce this effect, we choose to compound the volume obtained for the two complementary sequences since their respective PSF are rotated at 90 degrees. The OPW compounding strategy relies on summing the HV PSF and the VH PSF side lobes that are not located in the same location in space and thus reducing the final side lobe levels on the combined (HV + VH) PSF.

Since the aim of this study is to perform 4D ultrafast Doppler imaging, we compared the performances of the RCA matrix approach against the performances of the 2D matrix approach for a fixed frame rate, i.e., and equal number of emissions figure 2.

As described above, RCA matrices can only transmit a plane wave tilted either along the X or the Y axis, unlike the 2D matrix. Hence, an identical number of tilted planes were distributed uniformly in the [θX, θY] space for the 2D matrix, whereas tilted planes were distributed only in the [θX, 0]∪[0, θY] space in the case of the RCA matrix. The same angular range was used in both cases.

A c c e p t e d M a n u s c r i p t

PSF simulation

The numerical study focuses on the two configurations stated above: a 32 x 32 2D matrix and a 32 + 32 RCA matrix. The RF (Radio Frequency) signals were generated using Field II [START_REF] Jensen | Calculation of Pressure Fields from Arbitrarily[END_REF]. A dedicated delay-and-sum beamforming was applied on the RF signals to compute the Bscan volume of a unique scatterer. The scatterer was centered in front of the transducer at a depth equal to 1.5 times the transducer aperture (i.e., 14 mm). A Hanning apodization with a F-number of 1.5 was used in reception. The excitation of the transducer consisted of 2 periods of a 5 MHz sinusoid with Hanning weighting. As described in figure 2a, 17 tilted plane waves were equally distributed within an angular range of -8 to +8 degrees for the 2D and the RCA matrix, corresponding to a volume rate of 2.9 KHz and transmitted in sequential order by increasing angle.

Ultrasound scanner and probe

The experimental setup was based on a customized, programmable, 1024-channel ultrasound system [START_REF] Provost | 3D ultrafast ultrasound imaging in vivo[END_REF]) used to drive a PZT piezocomposite based 32 x 32 2D matrix array transducer centered at 5 MHz with a 300-µm pitch (Vermon, Tours, France).The system was composed of four 256-transmit/128-receive-channels units initially designed for 2D imaging on an Aixplorer system (Supersonic Imagine, Aix-en-Provence, France), assembled into a 1024-channel device with synchronized boards. The total channel count was thus 1024 in transmit and 512 in receive. Because the receive channels were multiplexed to 1 of 2 transducer elements, each emission was repeated twice to synthetize a total of 1024 receive channels.

Emulated RCA matrix

Figure 3 shows the basic principle of the RCA emulation: in transmission, the probe is steered by flat delays along a line or a column whereas in reception, the RF signals are summed across the line or the column in post treatment. This approach allows for the direct comparison of the 32 + 32 RCA scheme to the optimal 32 x 32 2D matrix configuration, with all other parameters kept constant (transducer aperture, central frequency, bandwidth, field of view, etc.). The hypothesis that a line of RCA matrix share the same acoustic properties that a set of aligned elements on 2D matrix could be justified in term of active area. Indeed the kerf (distance between transducer elements) is significantly smaller than the pitch. The backpropagation distance, d Rx was calculated from the projection of the voxel 𝑣 onto each row element 𝑟 𝑖 (figure 4b) due to the long elements used in the RCA configuration [START_REF] Rasmussen | 3-D imaging using row-column-addressed arrays with integrated apodization -Part i: Apodization design and line element beamforming[END_REF]:

𝑑 𝑅𝑥 (𝑦, 𝑧, 𝑟 𝑖 ) = √𝑧 2 + (𝑦 -𝑟 𝑖 ) 2
The delays 𝜏 are then given by:

𝜏(𝑥, 𝑦, 𝑧, 𝜃, 𝑟 𝑖 ) = 1 𝑐 • 𝑑 𝑇𝑥 (𝑥, 𝑧, 𝜃) + 1 𝑐 • 𝑑 𝑅𝑥 (𝑦, 𝑧, 𝑟 𝑖 ),
where 𝑐 is the speed of sound.

A beamformed volume was computed by repeating the delay and sum algorithm for several voxels. A compounded beamformed volume was computed by coherent summation of the beamformed volumes obtained for each titled plane wave. The process was then applied to both the VH and the HV sequences (figure 2).

Full-view, ultrafast Doppler imaging volumes were obtained by first removing clutter from the beamformed baseband data using multidimensional spatiotemporal filtering [START_REF] Demené | Spatiotemporal Clutter Filtering of Ultrafast Ultrasound Data Highly Increases Doppler and fUltrasound Sensitivity[END_REF], applied in four dimensions (i.e., 3D space + time) [START_REF] Provost | 3-D ultrafast doppler imaging applied to the noninvasive mapping of blood vessels in Vivo[END_REF]. This technique consists in eliminating clutter signal by removing the low-frequency principal components of the temporal realizations of each voxel. Power Doppler volumes were then obtained by integrating the energy of the clutter-free baseband signal over time. Voxel-wise pulsed Doppler was obtained by calculating the short-time Fourier transform of the clutter-free baseband signals in the temporal direction using overlapping windows of 16 ms.

All beamforming and Doppler processing were performed on graphical processing units and were coded using CUDA C language within a Matlab (2016a, The MathWorks Inc., Natick, MA, USA) interface.

A typical acquisition including all data transfers and calculations such as volume beamforming and Doppler processing required a total of fewer than 5 minutes. Calculation times are dependent on imaging depth and sampling.

In vitro 4D ultrafast CNR and resolution

The experiment was conducted on a commercial phantom (ATS Inc., Small Parts Ultrasound Phantom Model 551, Bridgeport, CT. USA), consisting of anechoic cysts (cylinders with a radius equal to 2 mm) and echogenic wires (0.05 mm in diameter). CNR and resolution were quantified for both 2D and the RCA matrix. CNR was defined as the ratio between energy outside and inside the anechoic cyst. (1)

The denominator represents the anechoic inclusion signal and is computed as the mean of the voxels include within the cyst shape. The numerator represents the background signal. It is computed as the mean of the voxels outside the cyst [START_REF] Tiran | Multiplane wave imaging increases signal-to-noise ratio in ultrafast ultrasound imaging[END_REF]. The lateral resolution was defined as the lateral width of the main peak at -6 dB. A Hanning apodization with a F-number of 1.5 was fixed in reception.

The excitation of the transducer consisted of 2 periods of a 5-MHz sinusoid. The final volumes was obtained by compounding either 1, 4, 16, or 36 tilted plane waves equally distributed within an angular range of -16 to +16 degrees for the 2D and the RCA matrix and transmitted in sequential order by increasing angle figure 5. The pulse repetition frequency was equal to 16 kHz.

Figure 5 : Tilted plane waves represented by their angles of emission

A c c e p t e d M a n u s c r i p t

In vitro 4D ultrafast Doppler

The experiment was performed using a custom-made Doppler phantom with ATS blood mimicking fluid passing through a 4 mm diameter vein-mimicking tube and a hydraulic pumping systems ATS 700-D 115V. The excitation of the transducer consisted of 2 periods of a 5 MHz sinusoid. This Doppler phantom was imaged using 3 sequences with different number of angle transmit: 2, 4 and 5 as shown in table 2.

Each of these sequences of 525 frames were repeated 5 times. The pulse repetition frequency was fixed at 21 KHz, therefore the volume rate was inversely proportional to the number of angle. The beamformed volume was a parallelepiped of 10 x 10 x 16 mm. Power Doppler CNR was defined as the sum of all the voxels energy included inside the real tube shape divided by the energy outside as defined by the equation ( 1).

Angles 4 Angles 5 Angles θx [degrees

] -2 0 -2 2 0 0 -2 2 0 0 0 θy [degrees] 0 -2 0 0 -2 2 0 0 0 -2 2
Table 2 : Tilted plane wave angle configurations expressed in degree

In vivo 4D ultrafast Doppler

A specific sequence was developed for imaging the humeral artery of a healthy volunteer (table 3). The humeral artery was imaged using fives plane waves tilted at θX = [-4 4 0 0 0] and θY =[0 0 0 4 4] degrees (according to the convention fixed at figure 2a). The excitation of the transducer consisted of 2 periods of a 5 MHz sinusoid. A total of 1500 compounded volumes were acquired, at a volume rate of 1.

PSF simulation

The simulated PSF are presented along the different planes XY (figure 6b) and XZ (figure 6c). These PSF are computed for an equal number of tilted plane wave transmitted, and thus the same volume rate of 2.9 KHz. The tilted plane wave distributions are represented in the figure 6a. The PSF side lobe level associated to the 2D matrix is -37 dB while the HV + VH has a level of -16.5 dB. The lateral width at -6 dB in the 2D matrix case is 0.6 mm whereas the HV + VH has a 1 mm width.

As expected, The PSF of HV and VH are respectively elongated in the X and Y directions. In the elongated axe, the side lobe level is -11.5 dB and the lateral width is 2 mm whereas in the short axe (perpendicular to elongated axe), the side lobe level is -6.75 dB and the lateral width is 0.77 mm. The coherent summation of HV and VH allowed improving the PSF of HV and VH.

In figure 7, the PSF from the OPW sequence on RCA matrix has been computed for different set of tilted angle and thus different volume rate. The PSF quality improved rapidly with the number of plane waves up to 17. Beyond 17 angles, the PSF continues improving but at a slower rate (for this angular range configuration from -8 to +8 degrees) Figure 9a shows the image CNR measured experimentally as a function of the number of transmitted plane waves for both 2D matrix (blue curve) and emulated RCA matrix (red curve) at 13 mm depth.

The error bars are computed as the standard deviation over the pixel value inside the cyst. For a large number of transmits, the CNR was two times lower with the RCA than with the 2D matrix, which is consistent with the numerical simulation results.

Figure 9b, shows the lateral resolution as a function of depth for both 2D matrix (blue curve) and emulated RCA matrix (red curve) with fixed steering parameters: 16 plane waves transmitted and an angular range of -16 to +16 degrees. The lateral resolution decreased with depth and was 0.5mm with the 2D matrix and 0.6mm with the RCA. These results show that in this experimental context, the power Doppler CNR with the RCA is approximately the square root of the CNR of the 2D array. Nevertheless, by increasing the number of compounded angle, RCA can achieve power Doppler with enough CNR to recover the vein shape.

Regarding the Doppler spectrum, RCA gave similar frequency information but with less energy than the 2D matrix.

A c c e p t e d M a n u s c r i p t These in vivo results showed that ultrafast OPW compounding imaging using emulated RCA matrix offer sufficient performances to image an artery on power Doppler and gave accurate quantification on pulsed Doppler.

Discussion

In this paper we have investigated the performance of a RCA scheme with OPW compounding which allows the reconstruction of a large volume with only N + N channels at high volume rate. The study was deliberately focused on ultrafast Doppler imaging and comparisons with an equivalent 2D array.

Ultrafast imaging based on coherent plane wave compounding enables improving the image or volume quality by transmitting a limited number of plane waves that approximate a synthetic cylindrical focusing. In the RCA approach, however, when transmitting plane waves on one array, the symmetry between the emissions and receive focusing is broken leading to anisotropic PSF. OPW compounding was introduced in this study to compensate the anisotropy of the PSF. Jensen et al. recently published several studies describing a synthetic aperture approach for RCA matrices [START_REF] Rasmussen | 3-D imaging using row-column-addressed arrays with integrated apodization -Part i: Apodization design and line element beamforming[END_REF]. The synthetic aperture approach relied on a row-by-row element transmission with a reception of all the columns. To sweep through the whole aperture, the sequence is repeated for the N rows. Synthetic aperture allows the reconstruction of any kind of focusing from the acquired dataset including transmit + receive focusing and is thus very flexible. However, since the number of emissions is equal to the number of N elements, synthetic aperture imaging approach does not scale for high volume rate imaging, at least not for hundreds of rows but is relevant for real-time Bmode imaging.

The resolution and CNR performance of the RCA with OPW compounding were evaluated and compared to the 2D array performances. As expected, it consistently lags behind the full matrix array for a same aperture, framerate and number of emissions. As a rule of thumb, CNR of the RCA is around the square root of the 2D array (half in decibel).

Nevertheless, we have shown here the feasibility of using the RCA approach to image in vivo arterial blood flows in a human volunteer. Volumetric reconstruction of phantoms and a human artery vessel was performed and compared to the one obtained with the 2D matrix. An accurate Doppler spectrum was also successfully derived from these acquisitions but with less than 2D matrix sensitivity.

An important advantage of RCA arrays over 2D matrix, is the possibility to manufacture large arrays for example with more than 128 + 128 elements. A 2D matrix (128 x 128 elements) with an equivalent footprint would be too complex to manufacture and would require an electronics with very high channel count which is currently not available. A state of the art RCA array of 128 + 128 elements could therefore outperform in specific applications many of the current 2D matrix thanks to a larger aperture and/or a higher frequency.

Vector flow imaging could benefit strongly of ultrafast imaging performed with RCA arrays. Jensen et al presented initial in vitro results of vector flow imaging with a RCA prototype [START_REF] Lehrmann | 3-D Vector Velocity Estimation with Row-Column Addressed Arrays[END_REF].

Correia et al showed the possibility to obtain the complete in vivo volumetric reconstruction of vector flow using ultrafast plane waves imaging with a 2D matrix [START_REF] Correia | 4D ultrafast ultrasound flow imaging: in vivo quantification of arterial volumetric flow rate in a single heartbeat[END_REF]. This work could be extended to the use of RCA arrays and would allow a simple and less expensive way to perform 4D flow imaging.

Conclusion

The RCA scheme with OPW compounding is a promising technique for 4D ultrafast Doppler imaging with a low number of acquisition channels. The imaging quality was compared to an equivalent 2D matrix with same volume rate. CNR and resolution were found to be lower than using the 2D matrix with the same volume rate. However, the CNR and resolution could be improved by adding more angles at the cost of lower volume rate. We demonstrated the feasibility of using an RCA with OPW sequence in vivo on human artery. OPW is a promising approach to 4D clinical vascular application with a low number of channels.

A c c e p t e d M a n u s c r i p t

Figure 1 :

 1 Figure 1 : RCA matrix basic principle

Figure 2 :

 2 Figure 2 : Plane wave compounding for 2D matrix a) and OPW compounding for RCA matrix b)

Figure 3 :

 3 Figure 3 : A 2D matrix emulated in a RCA column for transmit a) and a 2D matrix emulated in a RCA line for receive b)

Figure 4 :

 4 Figure 4 : Forward a) and backward b) distances illustration for an HV sequence on a RCA matrix The forward distance d Tx calculation is based on the regular plane wave beamforming algorithm. It corresponds to the projection of the voxel 𝑣 on the tilted plane wave as shown in figure 4a (Montaldo et al. 2009). 𝑑 𝑇𝑥 (𝑥, 𝑧, 𝜃) = [𝑧 • 𝑐𝑜𝑠(𝜃) + 𝑥 • 𝑠𝑖𝑛 (𝜃)]

  7 KHz, corresponding to an acquisition time of 0.88 second. The beamformed volume was a parallelepiped of 10 x 10 x 24 mm cube. (ISPTA = 376 mW/cm² & MI = 0

Figure 6 :

 6 Figure 6 : Tilted plane waves represented by their angles of emission a) and theirs PSF cut view in two different planes: XY b), and XZ c), for an identical volume rate of 2.9 KHz

Figure 8 :

 8 Figure 8 : 2D matrix and RCA Bscan of echogenic wires a) and anechoic cysts b) A c c e p t e d M a n u s c r i p t

Figure 9 :

 9 Figure 9 : CNR a) and resolution b) for both 2D matrix (blue curve) and emulated RCA matrix (red curve)

Figure 10 :

 10 Figure 10 : Power Doppler cross views of the phantom flow with different number of compounded angles for the 2D matrix and the emulated RCA matrix

Figure 14 :

 14 Figure 14 : In vivo power Doppler on a volunteer human humeral artery using a 2D matrix a) and an emulated RCA matrix b)

Table 1 :

 1 Principal imaging schemes for 2D Matrices and RCA Matrices.

Table 4 :

 4 The main lobe lateral width and side lobe levels are synthesized in the table 4. Since HV and VH PSF are identical with a 90° rotation, only HV is represented in the table. PSF Main lobe lateral width and side lobe level

	2D matrix	HV + VH	on X	HV	on Y

Acknowledgment

Funding for this project was provided by the French National Research Agency (ANR -DS0413 2015 -ULTRAFAST 4D)


Anatole Jimenez
email: anatole.jimenez@espci.fr

Bruno Osmanski
email: bruno.osmanski@iconeus.fr

Denis Vivien
email: vivien@cyceron.fr

Mickael Tanter
email: mickael.tanter@espci.fr

Thomas Gaberel
email: thomas.gaberel@hotmail.fr

Thomas Deffieux
email: thomas.deffieux@espci.fr

Toward Whole-Brain Minimally-Invasive Vascular Imaging

Keywords: Blood flow imaging, Doppler imaging, contrastfree microvascular imaging I

INTRODUCTION

Traumatic brain injuries and aneurysm-related hemorrhages are significant causes of death and severe disabilities. Brain lesions significantly arise from secondary ischemic damage in comatose patients [1] [START_REF] Macdonald | Delayed neurological deterioration after subarachnoid haemorrhage[END_REF]. The maintenance of the cerebral blood supply of those patients is the cornerstone of neurocritical care, and a reliable and early ischemia detection method is critical for appropriate intervention. Unfortunately, existing clinical modalities do not fulfill this task [START_REF] Wintermark | Comparative overview of brain perfusion imaging techniques[END_REF]. Although ultrasensitive Doppler (UD) can provide good sensitivity to cerebral blood volume (CBV) [START_REF] Hingot | Early Ultrafast Ultrasound Imaging of Cerebral Perfusion correlates with Ischemic Stroke outcomes and responses to treatment in Mice[END_REF] [5], it remains challenging to perform through the skull. In this work, we investigate how a burr hole needed for intracranial pressure (ICP) monitoring [START_REF] Forsyth | Routine intracranial pressure monitoring in acute coma[END_REF] could be used to monitor whole brain vascularization using a preclinical prototype.

II. METHODS

A. The preclinical prototype

We developed a preclinical prototype that consists of a phased array probe mounted on a motorized rotary unit, the motor's rotation axis aligned on the probe's depth axis (Fig 1a). The probe is a 64-element linear phased array with a 10-mm diameter circular aperture (5 MHz central frequency, 0.145 mm pitch, 9.1 mm aperture (L)). This prototype is not implantable. The probe was driven by an ultrasound research system (Vantage 128, Verasonics, USA). Raw radio frequency (RF) data were saved and processed offline.

We performed 3D Doppler scans by acquiring UD slices at different rotation angles. An UD image was acquired with full aperture using N-diverging waves corresponding to different virtual sources at position X = (x, z) set horizontally [START_REF] Papadacci | High Contrast Ultrafast Imaging of the Human Heart[END_REF] at minimum ψmin and for φwide-wide:

𝑿 𝒌 = (𝐿 -𝑧 tan 𝜓 ) -½ + , sin , 𝑘 = 1, … , 𝑁 (1)
In phase Quadrature (IQ) images were reconstructed on a polar grid (elevation φ, radius R), (dφ, λ/2) = (0.5°, 0.15 mm) from the RF data. The blood flow signal was retrieved from a block of frames using spatiotemporal filtering based on singular value decomposition (SVD) [START_REF] Demené | Spatiotemporal Clutter Filtering of Ultrafast Ultrasound Data Highly Increases Doppler and fUltrasound Sensitivity[END_REF]. The filtered frames were squared and averaged into a final ultrasensitive power Doppler image.

Acoustic parameters were set to N = 12 virtual sources of φwide = 90° wide and ψmin = 10° giving, using (1), X = ({-4.0, -3.3, -2.5, -1.8, -1.1, -0.4, 0.4, 1.1, 1.8, 2.5, 3.3, 4.0}, 9.1) mm at a 6-kHz pulse repetition frequency for Rim = 100 mm deep.

The wavelength resolution was λ=0.3 mm. A 1-s long UD Funded by AAPG2020 BrainCam Agence Nationale de Recherche acquisition at a 500-Hz frame rate, removing the first 10% eigenvalues from SVD, was set for Doppler imaging.

We used 89 slices with a step dθ = 2° for a full scan. Slice overlapping in elevation was estimated based on the mean piezo elements' height ℎ = ∫ sin 𝑥 𝑑𝑥 = 6.4 mm with max height h = 10 mm, focal distance f = 55 mm, and elevation at focus elf derived from diffraction 𝑒𝑙 = =2.6 mm (Fig 1b).

For deep imaging, high power transmission is needed. Because of the small probe's aperture, the safety heating limit of +2°C for implantable device was surpassed. We thus implemented a pause of 11 s between each slice acquisition, limiting heating to +1.6°C. The total duration was 18 min, which remains acceptable for our application B. 3D Doppler scan processing 3D reconstruction was performed by interpolating data on a cartesian grid. UD slices were first corrected with a depth attenuation compensation. Then we interpolated the data, from a polar slice with rotation, a spherical grid (azimuth θ, elevation φ, radius φ), (dθ, dφ, λ/2) = (2°, 0.5°, 0.15 mm) to cartesian grid (lateral x, elevation y, depth z), (0.25, 0.25, 0.25) mm 3 , using Delaunay triangulation (scatteredInterpolant, MatLab).

Relative 3D Doppler scan in percentage was computed from two successive 3D Doppler scans (𝑟𝑒𝑙𝑎𝑡𝑖𝑣𝑒𝑠𝑐𝑎𝑛 / = × 100). Scans were first registered with intensitybased rigid registration (imregtform, MatLab) using large vessels from thresholding. When registration was not possible, ROI comparisons were carried out and the nonuniform background amplitudes stemming from different depth of the probe were compensated.

Signal, image processing and 2D slice rendering were performed with MatLab (MathWorks, MA, USA). 3D rendering was performed with Amira software (Thermo Fisher Scientific Inc., MA, USA).

III. IN VITRO FEASIBILITY

3D Doppler scan feasibility with our non-implantable preclinical prototype was first assessed in vitro on a phantom with tubes (MODEL ATS 524, CIRS, Norfolk, USA). The setup consisted of the prototype being placed above a water tank with the phantom immersed inside. The phantom's tubes were filled with blood mimicking fluid and connected to a pump to create flow at 4 cm/s. The phantom could be moved to image flow at different locations relative to the probe. 3D Doppler scan (18 min) was performed for each configuration (Fig 2a).

Then, the Contrast-to-Noise Ratio (CNR) of the Doppler signal in tubes was computed for all configurations. All CNR slices were averaged on a single slice, roughly fitting an 80°wide sector field of view (Fig 2b).

IV. IN VIVO PROTOCOL

In vivo imaging was carried out on six pigs with the collaboration of GIP Cyceron (Caen, France). The pig animal model was chosen because it offers a relatively large brain (5 x 6 x 7 cm3) and is readily available. The animal was anesthetized and sedated, and a conventional 15-mm diameter burr hole was drilled through the skull at Bregma +20 mm and midline -12 mm up to the dura mater leaving brain protective layers intact. 2D real-time Bmode was displayed for positioning. Perfusion CT with iodinated contrast injection was also performed for stroke assessment with a CT scanner (Model, GE Healthcare, Chicago, USA).

Global blood flow variations were induced. CBV was decreased through mechanical hyperventilation inducing hypocapnia [START_REF] Laffey | Hypocapnia[END_REF] and later increased through acetazolamide injection, a vasodilator [START_REF] Vagal | The acetazolamide challenge: techniques and applications in the evaluation of chronic cerebral ischemia[END_REF]. A stroke model was carried out to induce ischemia. Iron chloride application on the carotid makes blood clots form in the vessel by oxidation. Then compressing the carotid, blood clots migrate into the brain vasculature [START_REF] Bonnard | Ferric Chloride-induced Thrombosis Mouse Model on Carotid Artery and Mesentery Vessel[END_REF]. A neurosurgeon assessed the ischemic core with perfusion CT [START_REF] Wintermark | Cerebral perfusion CT: technique and clinical applications[END_REF] before and after iron chloride application.

V. RESULTS

A. Coverage

First, the coverage of the brain through a burr hole was assessed. Using an anatomical CT scan, we registered a pig brain atlas (The Pig Imaging Group, University of Illinois Urbana-Champaign, USA) on the 3D Doppler scan. We measured about 60% brain coverage and up to 80% with good probe positioning. The field of view of an 80°-wide cone was retrieved in vivo (Fig 3).

B. Global cerebral blood volume variations

Then, global CBV variations were measured. First, hyperventilation induces a CBV decrease, and later, acetazolamide induces an increase. Two 3D Doppler scans were performed at baseline and when blood variations were established, allowing us to estimate relative 3D Doppler scans in percentage. In this example, we observed, as expected, a global decrease of 2.5% following mechanical hyperventilation (Fig. 4a) and a global increase of 6% following acetazolamide injection (Fig. 4b). The relative Doppler found in the whole brain was in accordance with induced CBV variations in 4/6 pigs for decrease and 4/5 pigs for increase

C. Stroke

Finally, stroke detection was assessed following induced ischemia. For this pig, an ischemic stroke formed with perfusion dropping by 60% (14.0 to 4.9 a.u. CT's blood volume). The assessed ischemic core was found on both 3D Doppler scans performed before and after iron chloride application. In the core, fewer vessels were observed, and a decrease in Doppler signal of 25% was measured (0.161 to 0.105 normalized Doppler signal).

The formation of the ischemic core was explored with a higher frame rate. Right after iron chloride application, waiting for the ischemia to establish, a 2D Doppler movie of 60 frames in 10 min was acquired. In the assessed ischemic core, CBV was diminished by about 38%.

VI. DISCUSSION

The proof of concept of minimally-invasive whole-brain perfusion monitoring was assessed with a preclinical nonimplantable prototype. A coverage of 60% of the brain was obtained, which allowed global blood flow monitoring and stroke detection.

A scan duration of 18 min was chosen to stay below the 2 degrees safety limitation for implantable device. The resulting scan duration is quite long for real-time imaging but remains acceptable in the context of perfusion and ischemia monitoring. Brain coverage was sufficient on the pig brain with 62 ± 12 % up to 80 % coverage. However, still missing superficial, extreme frontal, and occipital regions depending on the probe's relative orientation with respect to the brain. Nonetheless, the pig brain is still half the depth of the human one, and deep imaging over 40 mm could not be evaluated because of the skull's interface with the lower base.

Multiple 3D Doppler scans were performed to assess longitudinal CBV variations, close to the clinical application of monitoring. However, in our protocol, CT scans performed between the 3D Doppler scans required manually removing and replacing the probe, introducing shifts in the probe orientation and position, which would not have happened with an implantable device. If post hoc volume registration could handle most of those cases, this posed, in some cases, significant challenges for correctly assessing small CBV variations due to the non-uniform amplitude of the Doppler images due to the probe geometry.

The stroke monitoring could only be reliably performed in one over six pigs due to a combination of iron chloride model issues, probe repositioning issues, or premature animal death. For this case, the 2D Doppler movie clearly showed the ischemic core with a decrease of 38% in the ischemic core. Registered 3D Doppler scans with intensity normalization found a decrease of 25%.

In future experiments, we plan to couple Doppler with Bmode, elastography, and explore ULM scan without having to move the probe between the longitudinal scans.

CONCLUSION

The proof of concept of minimally-invasive whole-brain brain perfusion monitoring was performed with a preclinical prototype on pigs. Scan duration and overall brain coverage allowed global perfusion monitoring and stroke detection. This device, which could be further miniaturized for implantation inside a burr-hole, could offer a promising brain perfusion monitoring tool for neurocritical care alongside ICP monitoring.

Fig 1 .

 1 Fig 1. The preclinical prototype: (a) schematic view of the linear phased array mounted on a motorized rotary unit with imaging field of view of the slice at azimuth θ = 0° and θ = θmotor, (b) consecutive slices in elevation overlapping.

Fig 2 .

 2 Fig 2. In vitro: (a) cartesian transformed 3D Doppler scan of four 8-mm diameter inclined tubes with flowing water with the 80°-wide and 100-mm deep conic field of view, (b) Contrast-to-Noise Ratio averaged in a slice from all tube's configurations.

Fig 3 .

 3 Fig 3. Pig brain coverage: (a) 3D Doppler scan of pig brain registered on pig brain atlas with the 80°-wide and 100-mm deep field of view with main vessels indicated ICV: internal cerebral vein, TsS: thalamostriate system, TS: transverse sinuses, ICA: internal carotid arteries and RM: rete mirabilis, (b) MR atlas.

Fig 4 ,

 4 Fig 4, Global cerebral blood volume (CBV) variations: relative 3D Doppler scan of an induced CBV (a) decrease via mechanical hyperventilation, (b) increase via acetazolamide injection

ACKNOWLEDGMENT

Thanks to the Agence Nationale de la Recherche, which funded this project (Project ANR 20 CE19-0026). The GIP Cyceron team managed the experimental part with pigs, and Iconeus team advised us on the ultrasound sequence.


Andrea Natale

GRADIENT FLOWS OF INTERACTING LAGUERRE CELLS AS DISCRETE POROUS MEDIA FLOWS

published or not. The documents may come

Gradient flows of interacting Laguerre cells as discrete porous media flows

Introduction

Voronoi and Laguerre tessellations are a popular approach to describe the neighborhood relations within particle systems, and therefore to model the particle interactions and dynamics. For instance, they have been used as a model for biological cells [START_REF] Honda | Geometrical models for cells in tissues[END_REF], to describe the regions of influence of different agents in territorial models in ecology [START_REF] Votel | Equilibrium configurations for a territorial model[END_REF], or also as discretization tools in continuum mechanics and fluid dynamics [START_REF] Thomas O Gallouët | Convergence of a Lagrangian discretization for barotropic fluids and porous media flow[END_REF][START_REF] Leclerc | Lagrangian discretization of crowd motion and linear diffusion[END_REF]. This article focuses on a specific class of models in which the particles evolution in space is governed by the gradient flow of an energy depending on a Laguerre decomposition of a given domain. Our primary interest is to investigate the high number of particles limit of these models, and show how to interpret the particle dynamics as a discrete version of porous media flow, reproducing its Lagrangian gradient flow structure [START_REF] Lawrence C Evans | Diffeomorphisms and nonlinear heat flows[END_REF]. Adopting this point of view, we will establish quantitative estimates for the convergence of the discrete models to their continous counterparts.

1.1. Problem description. Given N ∈ N, a tessellation L = {L 1 , . . . , L N } of a measurable set A ⊂ R d is a collection of a finite number of measurable subsets L i ⊂ A, called cells, such that int(L i ) ∩ int(L j ) = ∅ ∀ i = j , and ∪ i L i = A.

We call T N (A) the set of tessellations of A composed of N cells, and T s N (A) := {L ∈ T N (B) : B ⊆ A} . Given a compact domain Ω ⊂ R d with Lipschitz boundary, we study the dynamics of N interacting cells, represented by a tessellation in T s N (Ω), and whose location is parameterized by a vector of cell centers (or particles) X = (x 1 , . . . , x N ) ∈ (R d ) N .

Specifically, for a given subset of admissible tessellations L N (Ω) ⊆ T s N (Ω) and a fixed parameter ε > 0, we consider the following energy:

(1.1) F ε (X) := inf

L∈L N (Ω) i L i |x -x i | 2 2ε dx + i C i (|L i |) .
Date: April 6, 2023.

In practice, we will focus on the two cases where L N (Ω) is either T N (Ω) (the union of the cells is fixed) or T s N (Ω) (the union of the cells is optimal). Loosely speaking, the first term in (1.1) measures how close the tessellation is representative of the particle distributions. The second term is the energy of the tessellation which we suppose to depend only on the cell volumes, C i : [0, ∞) → R being a given function which might be different for each cell. For example, a common choice for C i used to model biological cells, is the (non)linear spring model

C i (a) = C(a) = K(a) |a -ā| 2 2 ∀ i ,
where K(a) > 0 is a scaling factor called bulk modulus and which might depend on a, and ā > 0 is the target volume which is assumed common to all cells. The dynamics of the cell centers on the time interval [0, T ] is governed by the gradient flow of F ε , with respect to a weighted l 2 metric on (R d ) N , defined as follows:

(1.2) Ẋ, Ẏ m 0 := i m 0 i ẋi , ẏi for all Ẋ, Ẏ ∈ (R d ) N and where m 0 := {m 0 1 , . . . , m 0 N } ∈ R N >0 . More precisely the evolution of the cell centers is given by a curve X : [0, T ] → (R d ) N satisfying

(1.3) Ẋ(t) = -∇ m 0 F ε (X(t))
for all t ∈ (0, T ), with a given initial condition X(0) = X 0 = (x 0 i ) ∈ Ω N , where ∇ m 0 denotes the gradient with respect to (1.2).

1.2. Relation with Laguerre tessellations and other discrete models. A Laguerre tessellation (also called power diagram or weighted Voronoi tessellation) is a tessellation of the domain L(X, w) = {L i (X, w)} i ∈ T N (Ω) parameterized by a set of particles X = (x 1 , . . . , x N ) ∈ (R d ) N and associated weights w = (w 1 , . . . , w N ) ∈ R N , and in which the cells L i (X, w) ⊂ Ω are defined as follows L i (X, w) := {x ∈ Ω : |x -x i | 2 -w i ≤ |x -x j | 2 -w j ∀j = i} .

We will also refer to x i as the cell center of the cell L i (X, w). Note that if w is a constant vector then we retrieve the standard Voronoi tessellation of Ω.

The Voronoi tessellation can be shown to be the unique minimizer of the energy (1.1) when C i = 0 and L N (Ω) = T N (Ω). In this case, the remaining term in (1.1) is sometimes referred to as centroidal Voronoi tessellation energy [START_REF] Liu | On centroidal Voronoi tessellation-energy smoothness and fast computation[END_REF], and the resulting model governed by (1.3) coincides with the Voronoi liquid described in [START_REF] Ruscher | The voronoi liquid[END_REF], which is a fluid dynamic interpretation of the Lloyd's algorithm [START_REF] Lloyd | Least squares quantization in PCM[END_REF][START_REF] Mérigot | Non-asymptotic convergence bounds for Wasserstein approximation using point clouds[END_REF].

In general, the optimal tessellations for problem (1.1) with L N (Ω) = T N (Ω) is always a Laguerre tessellation, and when L N (Ω) = T s N (Ω) is contained in one (in particular each cell of the optimal tessellation is the intersection of a Laguerre cell with a ball); see Section 2. Our discrete model is therefore related to cell evolution models based on Laguerre tessellations [START_REF] Gareth | Modeling growth in biological materials[END_REF]. In Voronoi cell models, for example, one imposes the tessellation to be Voronoi and the energy of the system only contains the second term in (1.1) (plus extra energy terms often related to the cell perimeter or nodes distance). Estabilishing continuous limit for such models is not trivial, partly because the energies considered are usually more complex than those we treat here. Available results are therefore limited to 1d [START_REF] John A Fozard | Continuum approximations of individual-based models for epithelial monolayers[END_REF] or formal calculations [START_REF] Alt | Nonlinear hyperbolic systems of generalized Navier-Stokes type for interactive motion in biology[END_REF]. In this light, the energy (1.1) leads to a modified dynamics which is however more amenable to theoretical analysis, at least for the case of energies only depending on the cell area.

1.3. Relation with Lagrangian discretizations of porous media flow. In the following, we focus on the case where

(1.4) C i (a) :=    U m 0 i a a if a > 0 +∞ otherwise
where U : [0, ∞) → R is a smooth strictly convex function with superlinear growth, with U (0) = 0. For this energy, the particle dynamics generated by (1.3) can be reinterpreted as a spatially discrete version of the Lagrangian formulation of the porous medium equation, describing the evolution of a density ρ : [0, T ]×Ω → [0, ∞) as the solution of the PDE:

(1.5)

   ∂ t ρ -div [ρ∇U (ρ)] = 0 on (0, T ) × Ω , ∇U (ρ) • n ∂Ω = 0 on (0, T ) × ∂Ω ,
where n ∂Ω denotes the unit normal to the boundary ∂Ω, with given initial conditions ρ(0, •) = ρ 0 .

To make this precise, let us fix a smooth strictly-positive reference density ν : Ω → (0, ∞), and consider the energy F : Diff(Ω) → R on the space of diffeomorphisms of Ω, defined by

(1.6) F(ϕ) = Ω U ν det(∇ϕ) • ϕ -1 . Given Φ ∈ Diff(Ω) such that (1.7) ρ 0 = ν det(∇Φ) • Φ -1 ,
one can check that, at least formally, the flow ϕ : [0, T ] → Diff(Ω) of the vector field -∇U (ρ) solves the gradient flow system

(1.8) ∂ t ϕ(t) = -∇ F F(ϕ(t)) ϕ(0) = Φ and ρ(t) = ν det ∇ϕ(t) • ϕ(t) -1 ,
where F := L 2 ν (Ω; R d ) and therefore ∇ F is the gradient computed with respect to the L 2 inner product weighted by ν (see [START_REF] Lawrence C Evans | Diffeomorphisms and nonlinear heat flows[END_REF] for details on this interpretation, or also Section 3).

Let us now fix a tessellation T N = {T i } i ∈ T N (Ω), and let F N be the space of piecewise constant functions on the triangulation with values in R d , i.e.

F N := {ϕ X ∈ F : ϕ X (x) = x i ∈ R d for a.e. x ∈ T i }. Any element ϕ X ∈ F N can be identified with a collection of particles X = (x i ) i ∈ (R d ) N
given by the collection of images of the cells T i . Then, a general strategy to construct a particle discretization of equation (1.8)

is to look for solutions ϕ X : [0, T ] → F N of (1.9) ∂ t ϕ X (t) = -∇ F N F(ϕ X (t))
where F : F N → R is an approrpiate discrete version of F, and with ϕ X (0) being an approximation of Φ in F N . Different choices of F lead to different methods. We refer to [START_REF] Carrillo | Lagrangian schemes for Wasserstein gradient flows[END_REF] for a review of possible strategies in this context (see in particular [START_REF] José | Numerical study of a particle method for gradient flows[END_REF], for an approach that is particularly close to the one we study in this article). Importantly, our discrete model (1.3) is equivalent to (1.9) for an appropriate variational regularization of the energy (see again Section 3 for details). Such a regularization is related to recent approaches based on semi-discrete optimal transport recalled in Section 1.4.

1.4. Relation with semi-discrete optimal transport. The energy in (1.1) admits a reformulation based on semi-discrete optimal transport [START_REF] Merigot | Optimal transport: discretization and algorithms[END_REF]. In fact, denoting by W 2 (µ 1 , µ 2 ) the L 2 -Wasserstein distance between two positive measures µ 1 , µ 2 ∈ M + (Ω), one can show that if x i = x j for all i = j,

(1.10) F ε (X) = inf a∈R N >0 ,η∈C W 2 2 i a i δ x i , η 2ε + i U m 0 i a i a i ,
where we used for C i the expression given in (1.4), and where C is a convex subset of M + (Ω): in particular, C = {dx} (where dx is the Lebesgue measure on Ω) if

L N (Ω) = T N (Ω), and C = {f dx : f : Ω → [0, 1]} if L N (Ω) = T s N (Ω) (see Appendix A for a proof).
Equation (1.10) highlights a link with a different approach which consists in regarding the energy of the system (1.6) as a function of the density, and considering its Moreau-Yosida regularization on the space of positive measures measures M + (Ω) with fixed total mass, with respect to the W 2 distance. This yields:

(1.11) Fε (X) := inf

ρ∈M ac + (Ω) W 2 2 i m 0 i δ x i , ρ 2ε + Ω U (ρ) ,
where M ac + (Ω) is the set of absolutely continuous positive measures on Ω. The orginal idea of these type of energy regularizations stems from the work of Brenier in [START_REF] Brenier | Derivation of the Euler Equations from a Caricature of Coulomb Interaction[END_REF] for discretizing the incompressible Euler equations, and its more recent reformulation using semi-discrete optimal transport in [START_REF] Mérigot | Minimal geodesics along volume-preserving maps, through semidiscrete optimal transport[END_REF][START_REF] Thomas | A Lagrangian scheme à la Brenier for the incompressible Euler equations[END_REF]. The energy (1.11) was then used in [START_REF] Thomas O Gallouët | Convergence of a Lagrangian discretization for barotropic fluids and porous media flow[END_REF][START_REF] Leclerc | Lagrangian discretization of crowd motion and linear diffusion[END_REF] to discretize the same nonlinear diffusion models considered in this article as well as the compressible (barotropic) Euler equations, and in [START_REF] Sarrazin | Lagrangian discretization of variational mean field games[END_REF] in the context of mean field games.

Note that the advantage of using (1.10) with respect to (1.11) is that the first implies a piece-wise constant density reconstruction (see equation (1.12) below) which is easier to handle numerically than the minimizers of (1.11) whose structure strongly depends on U . Moreover the variational definition (1.11) is also easier to generalize to more complex energies (see Remark 3.1). Finally, we remark that another discretization strategy similar to ours, which also leads to piece-wise constant densities on Laguerre cells, was proposed in [START_REF] Benamou | Discretization of functionals involving the Monge-Ampère operator[END_REF], but the aim of this latter work was not to derive quantitative convergence estimates as we do here. 1.5. Continuous limit. For a given initial condition X 0 ∈ Ω N such that x 0 i = x 0 j for all i = j, let t ∈ [0, T ] → X(t) solve system (1.3) with energy (1.4) and

m 0 i = T i ν = Φ(T i ) ρ 0
where ρ 0 satisfies (1.7) as before. Consider the discrete density

(1.12) μN (t, •) := i m 0 i |L i (t)| 1 L i (t)
where L i (t) is the unique (see Section 2.2) optimal tessellation for problem (1.1) associated with the positions X(t), and where 1 L i (t) is the indicator function of the set L i (t). Our main result states that μN converges to sufficiently smooth solutions ρ of (1.5) as long as the error in the initial conditions, measured by

(1.13) δ 2 N := i T i |Φ(x) -x 0 i | 2 ν(x)dx ,
and ε go to zero with appropriate rates.

Theorem 1.1. Let U : [0, ∞) → R be a smooth strictly convex function, with U (0) = 0, verifying the assumptions of Lemma 4.1 and suppose that there exist R, α > 1 and β > 0, such that 1 and verifies (1.7), and ∇U (ρ) is of class C 2,1 in space, uniformly in time. Moreover, let μN : [0, T ] × Ω → [0, ∞) be the discrete solution defined in (1.12) via system (1.3), with energy defined using (1.4). Then

U (r) -inf U ≥ βr α ∀ r ≥ R . Suppose that ρ : [0, T ] × Ω → [0, ∞) is a strong solution to (1.5), such that ρ 0 : Ω → [ρ min , ∞) with ρ min > 0 is of class C 1,
(1.14) max t∈[0,T ] Ω U (μ N (t, •)|ρ(t, •)) ≤ C δ 2 N ε + ε p-1 ,
where p = min{2, α}, U (r|s) := U (r) -U (s) -U (s)(r -s) for all r ≥ 0 and s > 0, and where the constant C > 0 only depends on

sup t∈[0,T ] ∇U (ρ(t, •) C 2,1 , ρ 0 C 1,1 , ρ min , diam(Ω), U , T , and d.
The proof is contained in Section 4. Just as in [START_REF] Thomas O Gallouët | Convergence of a Lagrangian discretization for barotropic fluids and porous media flow[END_REF], it relies on a Grönwall argument applied on an appropriately constructed modulated energy (as in the classical approach to obtain weak strong stability results on (1.5); see, e.g., Chapter 5 in [START_REF] Constantine M Dafermos | Grèce Mathématicien[END_REF]). The error estimate we prove is actually stronger than (1.14), and it is given explicitly in equation (4.16) (see also Section 4.7 for an extension in the presence of external potentials). In particular, the same bound holds for the error in the L 2 ν norm between the exact Lagrangian flow (according to the interpretation described in Section 1.3) and the discrete one, given by the collection of the particle trajectories. Finally, note that if for example we set (1.15) x 0 i = arg min

y∈R d T i |Φ(x) -y| 2 ν(x)dx ,
then δ N is just the L 2 ν projection error of Φ onto F N , the space of piece-wise constant vector fields on T N . Hence, denoting by

h N = max i diam(T i ) the largest cell diameter, we have δ N ≤ ∇Φ ∞ (ρ 0 [Ω]) 1/2 h N , where ρ 0 [Ω] is the integral of ρ 0 over Ω.

Analysis of the discrete model

In this section, we describe in more detail the discrete model (1.3), in particular we explain the link with Laguerre tessellations and provide an explicit formula for the gradient of the energy with respect to the particle positions. While much of the analysis follows the same lines as in [START_REF] Leclerc | Lagrangian discretization of crowd motion and linear diffusion[END_REF], or is based on standard arguments from semi-discrete optimal transport (see Appendix A or [START_REF] Merigot | Optimal transport: discretization and algorithms[END_REF], for example), we provide the proofs of the main statements without using directly optimal transport theory, for the sake of clarity and to make the discussion self-contained.

2.1. Internal energy. We suppose that the functions C i defining the energy of the tessellation are given by

C i (a) :=    U m 0 i a a if a > 0 +∞ otherwise
where U : [0, ∞) → R is a smooth strictly convex function with superlinear growth, with U (0) = 0. From this, it is easy to deduce that C i is strictly convex, decreasing and

C i (a) → m 0 i U (0) as a → +∞. As a consequence C * i : R → (-∞, ∞] is also strictly convex. More precisely, dom(C * i ) ⊆ (-∞, 0] and C * i is an increasing diffeomorphism between (-∞, 0) and (-∞, -m i U (0)) and (C * i ) (0) = +∞ (since C i is finite on (0, ∞)).
The pressure function associated with U is the strictly increasing function

P : [0, ∞) → [0, ∞), defined by (2.1) P (r) = rU (r) -U (r) if r > 0 0 if r = 0
and satisfying P (r) = rU (r) for all r > 0. This is related to

C i and C * i by (2.2) C i (a) = -P m 0 i a and (C * i ) (w) = m 0 i P -1 (-w) ,
for all a > 0 and w < 0.

Dual formulation.

From now on we suppose that X = (x 1 , . . . , x N ) ∈ (R d ) N is given and that x i = x j for all i = j. First we rewrite the energy (1.1) as follows:

(2.3) F ε (X) = inf L∈L N (Ω),a∈R N sup w∈R N i L i |x -x i | 2 2ε dx + C i (a i ) + w i 2ε (a i -|L i |) .
We obtain the dual problem by swapping the inf and the sup,

D ε (X) := sup w∈R N inf L∈L N (Ω),a∈R N i L i |x -x i | 2 2ε dx + C i (a i ) + w i 2ε (a i -|L i |) = sup w∈R N inf L∈L N (Ω) i L i |x -x i | 2 -w i 2ε dx -C * i - w i 2ε . Let φ(w; •) : Ω → R be the function defined by φ(w; x) = min i |x -x i | 2 -w i . If L N (Ω) = T s N (Ω), then L ∈ T N ( Ω)
for some Ω ⊆ Ω, and

i L i |x -x i | 2 -w i 2ε dx ≥ 1 2ε Ω φ(w; x) dx ≥ 1 2ε Ω min(0, φ(w; x)) dx = i L s i (X,w) |x -x i | 2 -w i 2ε dx , where (2.4) L s i (X, w) := L i (X, w) ∩ B(x i , (w + i ) 1/2 ) , √ w i x i Figure 1.
An example of optimal tessellation with cells constructed via equation (2.4).

where

B(x i , (w + i ) 1/2 ) is the closed ball of radius w 1/2 i
if w i ≥ 0, and ∅ otherwise. On the other hand if L N (Ω) = T N (Ω), by the same argument we obtain

i L i |x -x i | 2 -w i 2ε dx ≥ i L i (X,w) |x -x i | 2 -w i 2ε dx . Therefore (2.5) D ε (X) = sup w∈R N D ε (X; w) ,
where

D ε (X; w) := i L * i (X,w) |x -x i | 2 -w i 2ε dx -C * i - w i 2ε .
and

L * i (X, w) = L i (X, w) if L N (Ω) = T N (Ω) L s i (X, w) if L N (Ω) = T s N (Ω)
.

A graphical representation of an optimal tessellation in the case

L N (Ω) = T s N (Ω) is given in Figure 1. Proposition 2.1. The function D ε (X; •) : R N → (-∞, ∞] is concave and C 1 on its effective domain dom(D ε (X; •)) ⊆ [0, ∞) N . In particular for all w ∈ (0, ∞) N , ∂ w i D ε (X; w) = 1 2ε (C * i ) - w i 2ε - |L * i (X, w)| 2ε ∀ i.
Furthermore there exists a unique w * ∈ (0, ∞) N such that ∂ w i D ε (X; w * ) = 0 for all i, or equivalently

(2.6) P m 0 i |L * i (X, w)| = w i 2ε ∀ i ,
and therefore maximising D ε (X; •).

Proof. The concavity and C 1 regularity of D ε (X; •) can be proven using standard arguments from the theory of semi-discrete optimal transport [START_REF] Merigot | Optimal transport: discretization and algorithms[END_REF]. In particular, consider the function

Q(w) := i L * i (X,w) |x -x i | 2 -w i 2ε dx
and observe that for any w,

w ∈ R N Q( w) ≤ i L * i (X,w) |x -x i | 2 -wi 2ε dx = Q(w) + i w i -wi 2ε |L * i (X, w)| .
This shows that the super-differential of Q at w is not empty since the vector

(-|L * i (X, w)|/(2ε)) i ∈ ∂ + Q(w)
, and therefore Q is concave. Furthermore, since |L * i (X, w)| is a continuous function of w (see, e.g., Proposition 38 in [START_REF] Merigot | Optimal transport: discretization and algorithms[END_REF]), Q (and therefore D ε (X; •)) is necessarily C 1 . Existence of maximizers can be shown oberving that the function D ε (X; •) is coercive since for any j,

i L * i (X,w i ) |x -x i | 2 -w i 2ε dx ≤ Ω |x -x j | 2 -w j 2ε dx ≤ diam(Ω) 2 2ε |Ω| - |Ω| 2ε w j .
Since (C * i ) (0) = +∞ we get that necessarily w i > 0, so the optimality conditions hold. Uniqueness of maximizers is a consequence of the strict convexity of the functions C * i . From Proposition 2.1, we can deduce the equivalence with the primal problem (1.1), and the existence and uniquess of solutions for this latter. In fact, denoting by w * the unique solution of

∂ w i D ε (X; w * ) = 0 for all i, then D ε (X; w * ) = i L * i (X,w * ) |x -x i | 2 2ε dx - w * i 2ε (C * i ) - w * i 2ε -C * i - w i 2ε = i L * i (X,w * ) |x -x i | 2 2ε dx + C i (|L * i (X, w * )|) ≥ F ε (X).
Since w * is the unique maximizer of the dual problem (2.5), this also implies that the tessellation L * (X, w * ) := {L * i (X, w * )} i is the unique solution to problem (1.1). In the following we will need the optimality conditions in Proposition 2.1 in the following alternative form: Lemma 2.2. Let L = {L i } i solve problem (1.1), with L N (Ω) equal to either T N (Ω) or T s N (Ω). Then for any smooth vector field u :

Ω → R d with u • n ∂Ω = 0 on ∂Ω, i L i x -x i ε • u(x) dx = i L i P m 0 i |L i | - |x -x i | 2 2ε divu(x) dx
Proof. Consider again the function φ(w;

•) : Ω → R defined by φ(w; x) = min i |x - x i | 2 -w i . If L N (Ω) = T N (Ω) we have i L i x -x i ε • u(x) dx = Ω 1 2ε ∇ x φ(w; x) • u(x) dx = i L i w i 2ε - |x -x i | 2 2ε divu(x) dx
and we conclude using Proposition 2.1, and the relation between (C * ) and P in (2.2). If L N (Ω) = T s N (Ω) we just need to replace φ(w; •) with min{φ(w; •), 0}. 2.3. Discrete dynamical system. Let us introduce the set of particle configurations where at least two particles share the same location:

∆ N = {X = (x 1 , . . . , x N ) ∈ (R d ) N : ∃ i, j such that x i = x j , i = j} . Proposition 2.3. The function X → F ε (X) is C 1 on (R d ) N \ ∆ N . Moreover, for any X ∈ (R d ) N \ ∆ N , (∇ m 0 F ε (X)) i = |L i | m 0 i x i -b i ε , where b i := 1 |L i | L i xdx ,
and where L = {L i } i is the unique minimizer of problem (1.1).

The proof is a slight adaptation of the arguments used in [START_REF] Leclerc | Lagrangian discretization of crowd motion and linear diffusion[END_REF] to prove an analogous result, and is therefore postponed to Section A.3 in the appendix.

In view of Proposition 2.3, the particle dynamics is governed by the system of ODEs

(2.7) ẋi (t) = - |L i (t)| m 0 i x i (t) -b i (t) ε
where we denote by L i (t) the optimal tessellation associated with the particle configuration X(t) = (x 1 (t), . . . , x N (t)) ∈ (R d ) N . By Proposition 2.3 the right-hand side of (2.7) is a continuous function of X(t) on the open set (R d ) N \ ∆ N , and therefore the system always admits solutions if X(0) ∈ (R d ) N \ ∆ N . We now show that such solutions are always defined for all times (again, adapting similar arguments from [START_REF] Leclerc | Lagrangian discretization of crowd motion and linear diffusion[END_REF]). For this, we will need the following lemma:

Lemma 2.4. Let L(X, w) the Laguerre tessellation of Ω associated with the position vector X ∈ (R d ) N and the vector of weights w. If |L i (X, w)| > 0 for all i = 1, . . . , N , then for all i, j = 1, . . . , N ,

|w i -w j | ≤ 2diam(Ω)|x i -x j | Proof. Since |L i (X, w)| > 0, there exists x ∈ Ω such that |x -x i | 2 -w i ≤ |x -x j | 2 -w j ,
for all j. Rearranging terms we obtain,

w j -w i ≤ |x j | 2 -|x i | 2 -2x • (x j -x i ) = (x j + x i -2x) • (x j -x i ) ≤ 2diam(Ω)|x j -x i |
Swapping the role of i and j we get the result.

Suppose that X(t) solves (2.7) on some interval [0, t * ) with X(0) ∈ Ω N \ ∆ N , and consider the decreasing function R(s) := P (1/s) for all s > 0. Then, with the same notation as above,

P m 0 i |L i (t)| -P m 0 j |L j (t)| ≥ C N |L i (t)| m 0 i - |L j (t)| m 0 j
where C N := |R (|Ω|/ m0 )| > 0, with m0 := min i m 0 i , and by the optimality conditions in Proposition 2.1 and then Lemma 2.4,

|L i (t)| m 0 i - |L j (t)| m 0 j ≤ |w i (t) -w j (t)| 2εC N ≤ diam(Ω) εC N |x i (t) -x j (t)| .
From these bounds we can deduce a lower bound on the distance between particles. In particular, using the fact that

x i -x j , b i -b j ≥ 0, we obtain (2.8) d dt |x i -x j | 2 2 = x i -x j , - x i -b i ε |L i | m 0 i + x j -b j ε |L j | m 0 j ≥ - |L i | m 0 i |x i -x j | 2 ε -|x j -b j ||x i -x j | |L i | m 0 i - |L j | m 0 j 1 ε ≥ - |Ω| m0 |x i -x j | 2 ε - diam(Ω) 2 ε 2 C N |x i -x j | 2 ,
where we omitted the time dependency of x i , b i and L i to simplify the notation. By a Grönwall inequality, this shows the long time existence of discrete solutions:

Lemma 2.5. If X(0) ∈ (R d ) N \ ∆ N the solutions to (2.7) are defined for all times t > 0. Moreover, if X(0) ∈ conv(Ω) N \ ∆ N , then X(t) ∈ conv(Ω) N \ ∆ N for all t > 0, where conv(Ω) denotes the convex hull of Ω.

Proof. Denote by π : R d → conv(Ω) the Euclidean projection of x onto conv(Ω) and by

d 2 (x) = |x -π(x)| 2 the square distance of x from conv(Ω). Then since b i (t) ∈ Ω, d dt d 2 (x) 2 = ẋi , x i -π(x i ) ≤ - |Ω| m0 ε d 2 (x).
Using Grönwall's lemma on this inequality and on (2.8), we obtain the result.

Lagrangian formulation of porous media and link with the discrete model

In this section we describe formally the gradient flow structure of the porous medium equation (1.5) in Lagrangian variables [START_REF] Lawrence C Evans | Diffeomorphisms and nonlinear heat flows[END_REF]. We then use this to reinterpret (1.3) as a discretization of such a system preserving its gradient flow structure. Note that at the Eulerian level, the gradient flow interpretation of porous media corresponds to the Wasserstein gradient flow formulation originally put forward by Otto [START_REF] Otto | The geometry of dissipative evolution equations: The porous medium equation[END_REF]. Note also that the discussion in this section is essentially independent of the convergence proof in the next section, but it sheds light on the link between to the continuous and discrete systems and justifies the construction of the modulated energy in Section 4.2.

Let us denote by M + (Ω) the set of positive measures on Ω. Given a measurable map Ψ : Ω → R d , and a measure ρ ∈ M + (Ω) the pushforward of ρ by Ψ is the measure

Ψ # ρ ∈ M + (R d ), satisfying f dΨ # ρ = f (Ψ(x))dρ(x)
for all f ∈ C 0 (R d ), the space of continuous functions vanishing at infinity.

Consider now two reference measures σ, ν ∈ M + (Ω) with smooth and strictly positive densities with respect to the Lebesgue measure on Ω, denoted dx in the following, and let us define the energy F :

L 2 ν (Ω; R d ) → R by F(ϕ) =    Ω U dϕ # ν dσ dσ if ϕ # ν σ +∞ otherwise
.

The gradient of this energy with respect to the L 2 ν metric, at a given configuration ϕ ∈ Diff(Ω), can be defined as follows. Consider a smooth curve s ∈

(-ε, ε) → ϕ(s) ∈ Diff(Ω) ⊂ L 2 ν (Ω; R d ) such that ϕ(0) = ϕ and dϕ(s)/ds| s=0 = δϕ. Then, (3.1) ∇ L 2 ν F(ϕ), δϕ = d ds s=0 F(ϕ s ) = d ds s=0 U dϕ # ν dσ • ϕ(s)dν = ∇U dϕ # ν dσ • ϕ, δϕ ,
where •, • denotes the inner product on L 2 ν (Ω; R d ). Then, we can interpret any smooth curve of diffeomorphisms ϕ : [0, T ] → Diff(Ω) satisfying

   ∂ t ϕ(t) = -∇U d[ϕ(t) # ν] dσ • ϕ(t) ϕ(0) = Φ
as the gradient flow of F with respect to the L 2 ν metric, starting at Φ ∈ Diff(Ω). Denoting ρ(t) = ϕ(t) # ν and ρ 0 = Φ # ν, then ρ(t) statisfies the continuity equation with velocity field -∇U (dρ(t)/dσ) tangent to the boundary, i.e.

           ∂ t ρ(t) -div ρ(t)∇U dρ(t) dσ = 0 ∇U dρ(t) dσ • n ∂Ω = 0 ρ(0) = ρ 0 .
In order to link this formulation with the discrete model, for any diffeomorphism ϕ ∈ Diff(Ω), let us denote

λ := dϕ -1 # σ dσ ,
or equivalently ϕ # (λσ) = σ. Then,

ϕ # ν = ϕ # 1 λ dν dσ λσ = 1 λ dν dσ • ϕ -1 σ , which implies
F(ϕ) = U 1 λ dν dσ λdσ .
This suggests the definition of the following regularized energy

(3.2) F ε (ϕ) := inf λ W 2 2 (ϕ # (λσ), σ) 2ε + U 1 λ dν dσ λdσ ,
where W 2 (•, •) denotes the L 2 -Wasserstein distance on between positive measures with equal mass (see Appendix A for a precise definition). The main reason for using the regularization (3.2) is that that this is well-defined on the whole space L 2 ν (Ω; R d ). In turn, this allows us to make sense of the variational structure of the system without regularity assumptions on the flow, and reproduce this at the discrete level. In order to do this, let us fix a tessellation T = {T i } i of Ω, and for any given vector X = (x i ) i ∈ (R d ) N consider the piece-wise constant map ϕ X : Ω → R d such that ϕ X (x) = x i ∈ R d for a.e. x ∈ T i . Then,

ϕ X # λσ = i a i δ x i where a i = T i λdσ and (3.3) F ε (ϕ X ) := inf λ W 2 2 ( i a i δ x i , σ) 2ε + U 1 λ dν dσ λdσ .
Optimizing over λ we find that for any cell T i there exist a constant c i such that almost everywhere on T i , dν/dσ = c i λ. This implies that

m 0 i := T i dν = c i T i λdσ = c i a i =⇒ 1 λ dν dσ = i m 0 i a i 1 T i , which replaced into (3.3) gives F ε (ϕ X ) := inf a W 2 2 ( i a i δ x i , σ) 2ε + i U m 0 i a i a i .
In the case where σ = dx this coincides with F ε (X) for the case L N (Ω) = T N (Ω). The proof of the equivalence is contained in Appendix A.2.

In the case where ϕ : Ω → Ω = Ω the constraint ϕ # (λσ) = σ is not appropriate, as in this case we only have dϕ # λσ dσ = 1 Ω .

Hence, we define the regularized energy by

F ε (ϕ) := inf λ, dη dσ ≤1 W 2 2 (ϕ # (λσ), η) 2ε + U 1 λ dν dσ λdσ .
In the case where σ = dx, by similar computations as above, we find

F ε (ϕ X ) = F ε (X), with L N (Ω) = T s N (Ω). Finally, let us denote by F N := {ϕ ∈ L 2 ν (Ω; R d ) : ϕ(x) = x i ∈ R d for a.e.
x ∈ T i } the space of piecewise constant flows on the reference triangulation equipped with the L 2 ν metric. Then, in both the cases described above, the discrete dynamics in (1.3) coincides with the gradient flow

∂ t ϕ X = -∇ F N F ε (ϕ X ) .
Remark 3.1 (Generalizations to models with advected quantites). The type of energy regularization considered here can be easily generalized to models where multiple scalar functions and densities are advected by the flow, i.e. to the case where

F(ϕ) = U a 1 • ϕ -1 , . . . , a n • ϕ -1 , dϕ # ν 1 dσ , . . . , dϕ # ν m dσ dσ ,
where now U : R n+m → R, a i : Ω → R and ν j ∈ M + (Ω) are given scalar functions and positive measures, respectively. In fact, as before, this can be written as a single function of λ, since by a change of variables

F(ϕ) = U a 1 , . . . , a n , 1 λ dν 1 dσ , . . . , 1 λ dν m dσ λdσ .
Formally, writing the Hamiltonian equations corresponding to such energies,

∂ 2 tt ϕ(t) = -∇ F F(ϕ(t)) ,
one recovers, with appropriate choices of U , a large class of compressible fluid models including, e.g., the thermal shallow water equations or the full compressible Euler equations (see, e.g., [START_REF] Holm | The Euler-Poincaré equations and semidirect products with applications to continuum theories[END_REF][START_REF] Khesin | Geometric hydrodynamics and infinitedimensional Newton's equations[END_REF]). Then the same discretization strategy described in this section leads naturally to simple Lagrangian schemes for all of these models as well.

Convergence towards smooth solutions

In this section we prove Theorem 1.1, i.e. the convergence of discrete solutions associated to (1.3) towards smooth solutions of the equation

∂ t ρ -div ρ∇U (ρ) = 0.
The proof follows similar lines as the one used in [START_REF] Thomas O Gallouët | Convergence of a Lagrangian discretization for barotropic fluids and porous media flow[END_REF] to analyze the system associated to a different energy regularization, given by (1.11). It relies on the construction of an appropriate relative entropy, based on the Lagrangian point of view described in Section 3, and on two main technical lemmas.

Preliminary lemmas.

The first lemma we will need provides us with a way to control the relative pressure (4.1) P (r|s) := P (r) -P (s) -P (s)(r -s) by U (r|s) := U (r) -U (s) -U (s)(r -s). In particular, we will make the following assumption: there exists a constant A > 0 such that

(4.2) |P (r)| ≤ A U (r) ∀ r > 0 .
This assumption is trivially verified for the important case of power energies, i.e. when U (r) = r γ /(γ -1) with γ > 1, which corresponds to P (r) = r γ . It implies the following lemma, which is extracted from Lemma 3.3 in [START_REF] Giesselmann | Relative energy for the Korteweg theory and related Hamiltonian flows in gas dynamics[END_REF]. Proof. We have P (r|s) = (r -s) 2 1 0 (1 -θ)P ((1 -θ)s + θr) dθ and similarly for U (r|s). Hence, using equation (4.2),

|P (r|s)| ≤ (r -s) 2 1 0 (1 -θ)|P ((1 -θ)s + θr)| dθ ≤ A U (r|s) .
The second lemma is necessary to deal with the fact that the particles may exit the domain Ω, when this latter is not convex. For this reason, we will need to use an extension of the continuous density ρ on the whole space. It will be clear in the following that such an extension needs to verify the continuity equation with respect to an appropriate velocity field also defined on the whole space. Here we just report a simplified version of the statement of Lemma 4.1 in [START_REF] Thomas O Gallouët | Convergence of a Lagrangian discretization for barotropic fluids and porous media flow[END_REF], which was used precisely for this purpose.

Lemma 4.2. Let u : [0, T ] × Ω → R d be such that u • n ∂M = 0 on [0, T ] × ∂Ω, and ρ 0 : Ω → [ρ min , ∞) with ρ min > 0. If u is of class C 2,1
in space, uniformly in time, and ρ 0 is of class C 1,1 , then there exist ũ :

[0, T ] × R d → R d and ρ : [0, T ] × R d → R such that:
(1) ũ is an extension of u, i.e. ũ(t)| Ω = u(t) for all t ∈ [0, T ], and there exists a constant C > 0 only depending on d such that

(4.4) sup t∈[0,T ] ũ(t) C 2,1 ≤ C sup t∈[0,T ] u(t) C 2,1 ;
(2) the couple (ρ, ũ) solves the continuity equation:

∂ t ρ + div(ρũ) = 0 on [0, T ] × R d ,

4.2.

Main assumptions and relative entropy. Suppose that ρ : [0, T ] × Ω → (0, ∞) is a sufficiently smooth solution to equation (1.5). In particular, we suppose that u := -∇U (ρ) and ρ(0, •) = ρ 0 satisfy the assumptions of Lemma 4.2. We will denote by ρ : [0, T ] × R d → [0, ∞) and ũ : [0, T ] × R d → R d the extensions of ρ and u, respectively, outside the domain. By construction these satisfy the continuity equation on the whole space, but in general outside the domain Ω, ṽ := -∇U (ρ) = ũ .

For any N > 0, let

X N : t ∈ [0, T ] → (x i (t)) i ∈ (R d
) N be a solution to the discrete model (1.3), with given initial conditions x i (0) = x 0 i ∈ Ω. We suppose that for a given diffeomorphism Φ : Ω → Ω and a smooth reference density ν : Ω → (0, ∞), the initial density can be written as follows:

ρ 0 = ν det(∇Φ) • Φ -1 .
Then, given a fixed tessellation

T N = {T i } i ∈ T N (Ω), we denote (4.5) δ N := T i |Φ(x) -x 0 i | 2 ν(x)dx , m 0 i = Φ(T i ) ρ 0 = T i ν .
Let us introduce the time-dependent positive measures µ N : [0, T ] → M + (Ω) and μN : [0, T ] → M + (Ω) defined as follows

(4.6) µ N (t) := N i=1 m 0 i δ x i (t) , μN (t) := N i=1 m 0 i |L i (t)| 1 L i (t) dx .
For simplicity, in the following we will use μN (t) to denote also its density with respect to the Lebesgue measure dμ N (t)/dx. In order to define the relative entropy between the smooth and discrete solutions, we first introduce the flow of u, which is the curve of diffeomorphisms ϕ : [0, T ] × Ω → Ω verifying

∂ t ϕ(t, x) = u(t, ϕ(t, x)) , ϕ(0, x) = Φ(x) .
Then let ϕ N (t) := (ϕ(t, Φ -1 (x 0 i ))) i ∈ Ω N be the collection of the exact trajectories of the particles located at x 0 i at time t = 0. The relative entropy of the discrete solution with respect to the continuous one is defined as follows:

F ε (X N |ρ; t) := i L i (t) |x -x i (t)| 2 2ε dx + Ω U (μ N (t)|ρ(t)) + ϕ N (t) -X N (t) 2 m 0 2 .

4.3.

Time derivative of the relative entropy. We now compute the time derivative of the relative entropy and isolate the terms that need to be estimated. It will be useful to define the following quantity:

H(t) := R d U (ρ(t))d(µ N (t) -μN (t)) ,
where μN is extended by zero on R d . We will keep using this convention in what follows.

Let us start by rewriting the relative entropy as follows:

(4.7) F ε (X N |ρ; t) = F ε (X N (t)) + Ω P (ρ(t)) - R d U (ρ(t))dµ N (t) + H(t) + ϕ N (t) -X N (t) 2 m 0 2 .
We compute the time derivative of the terms of the right-hand side separately. For the first term, we write

d dt F ε (X N (t)) = i L i (t) x i (t) -x ε • ẋi (t)dx = -ẊN (t), ẊN (t) -ũ(t, X N (t)) m 0 + i L i (t) x i (t) -x ε • (ũ(t, x i (t)) -ũ(t, x))dx + i L i (t) x i (t) -x ε • u(t, x) dx ,
and by Lemma 2.2 we can write the last term on the right-hand side as follows:

i L i (t) x i (t) -x ε • u(t, x)dx = i L i (t) |x -x i (t)| 2 2ε dx - Ω P (μ N (t))divu(t) .
For the second term, using the continuity equation

-∂ t ρ = ∇ρ • u + ρdivu, we obtain d dt Ω P (ρ(t)) = - Ω P (ρ(t))ρ(t)divu(t) -∇P (ρ(t)) • u(t) = Ω P (ρ(t)) -P (ρ(t))ρ(t) divu(t) .
Finally, for the third term, using again the continuity equation

-∂ t ρ = ∇ρ• ũ+ ρdivũ, this time on R d , we have d dt R d U (ρ(t))dµ N (t) = i ∇U (ρ(t, x i )) • ẋi m 0 i + R d U (ρ(t))∂ t ρdµ N (t) = i ∇U (ρ(t, x i )) • ( ẋi -ũ(t, x i ))m 0 i - R d P (ρ(t))divũ(t)dµ N (t) = -ṽ(t, X N (t)), ẊN (t) -ũ(t, X N (t)) m 0 - R d P (ρ(t))divũ(t)dµ N (t) .
Reinserting these expressions into the time derivative of (4.7) and rearranging terms we obtain

(4.8) d dt F ε (X N |ρ; t) + d dt H(t) + ẊN (t) -ũ(t, X N (t)) 2 m 0 = 5 j=1 I j (t) ,
where the terms in the sum on the right-hand side are defined as follows:

I 1 (t) := i L i (t) x i (t) -x ε • (ũ(t, x i (t)) -ũ(t, x))dx , I 2 (t) := i L i (t) |x -x i (t)| 2 2ε divu(t, x)dx - Ω P (μ N (t)|ρ(t))divu(t) , I 3 (t) := R d P (ρ(t))divũ(t)d(µ N (t) -μN (t)) , I 4 (t) := ṽ(t, X N (t)) -ũ(t, X N (t)), ẊN (t) -ũ(t, X N (t)) m 0 , I 5 (t) := u(t, ϕ N (t)) -ẊN , ϕ N (t) -X N (t) m 0 .
4.4. Uniform estimates. In the following, for any given Lipchitz function f ∈ C 0,1 (Ω), we will denote by Lip(f ) its Lipschitz constant, and for any time-dependent function g ∈ C([0, T ]; C 0,1 (Ω)) we denote by Lip T g := sup t∈[0,T ] Lip(g(t, •)), and similarly for vector-valued functions.

We estimate separately the terms on the right-hand side of (4.8). We have

I 1 (t) ≤ Lip T (ũ) i L i (t)
|x i (t) -x| 2 ε dx and using Lemma 4.1,

I 2 (t) ≤ (Lip T (u) + A) i L i (t) |x -x i (t)| 2 2ε dx + Ω U (μ N (t)|ρ(t))divu(t) ,
To bound I 3 , let us introduce h := P (ρ)div(ũ). Then, for any λ > 0,

I 3 (t) = i h(t, x i (t))m 0 i - m 0 i |L i (t)| L i (t) h(t, x)dx ≤ i Lip T (h) m 0 i |L i (t)| L i (t) |x i (t) -x|dx ≤ Lip T (h) i L i (t) |x i (t) -x| q λ q qε dx + ε p-1 λ p p i m 0 i |L i (t)| p |L i (t)|
where p, q > 1 are conjugate exponents, i.e. 1/p + 1/q = 1. Recall that we supposed that there exist R, α > 1 and β > 0, such that

U (r) -inf U ≥ βr α ∀ r ≥ R .
Then, choosing p = min{2, α} we get q ≥ 2, and therefore

I 3 (t) ≤ Lip T (h) 2diam(Ω) q-2 qλ q i L i (t) |x i (t) -x| 2 2ε dx+ Lip T (h) ε p-1 λ p p |Ω|R p + β -1 Ω U (μ N (t)) -β -1 |Ω| inf U
where we used the fact that, since X N (0) ∈ Ω N \ ∆ N , by Lemma 2.5 x i (t) ∈ conv(Ω) (the convex hull of Ω) for all times t ≥ 0. Hence, (4.9)

I 3 (t) ≤ Lip T (h) C 1 λ q F ε (X|ρ; t) + C 2 λ p ε p-1
where (4.10)

C 1 := 2diam(Ω) q-2 q , C 2 := β -1 F ε (X(0)) + |Ω|(R p -β -1 inf U ) p .
In the following we will use (4.9) with λ = 1. However, using the same arguments to bound H(t) and choosing λ appropriately, we get

(4.11) |H(t)| ≤ 1 2 F ε (X|ρ; t) + C 2 λp ε p-1 where λ = 2 C 1 max{Lip T (h), 1} 1 q .
Finally we observe that

I 4 (t) + I 5 (t) = ṽ(t, X N (t)) -ṽ(t, ϕ N (t)), ẊN (t) -ũ(t, X N (t)) m 0 + ũ(t, ϕ N (t)) -ũ(t, X N (t)), ẊN (t) -ũ(t, X N (t)) m 0 + ϕ N (t) -X N (t), ẊN (t) -ũ(t, X N (t)) m 0 ≤ 3 2 Lip(ṽ) 2 + Lip(ũ) 2 + 1 X N (t) -ϕ N (t) 2 m 0 + ẊN (t) -ũ(t, X N (t)) 2 m 0 2 .
Remark 4.3. The dependence on ε p-1 of the final estimate (1.14) is essentially due to the bounds on I 3 and H given in (4.9) and (4.11), respectively. Both bounds can be directly related to the natural requirement that two possible Eulerian reconstructions given in (4.6) converge (weakly) to the same limit.

4.5. Grönwall argument. Reinserting the estimates above into (4.8), we obtain

(4.12) d dt F ε (X N |ρ; t) + d dt H(t) + ẊN (t) -ũ(t, X N (t)) 2 m 0 2 ≤ C 3 F ε (X N |ρ; t) + C 4 ε p-1
Let G(t) := F ε (X N |ρ; t) + H(t) and observe that equation (4.11) implies

-F ε (X N |ρ; t) ≤ 2H(t) + 2C 2 λp ε p-1
and adding 2F ε (X N |ρ; t) on both sides we obtain

F ε (X N |ρ; t) ≤ 2G(t) + 2C 2 λp ε p-1
Substituting this into (4.12), we obtain

d dt G(t) + ẊN (t) -ũ(t, X N (t)) 2 m 0 2 ≤ 2C 3 G(t) + C 5 ε p-1
Hence by Grönwall's inequality we get (4.13) F ε (X N |ρ; t)

+ t 0 ẊN (s) -ũ(t, X N (s)) 2 m 0 2 ds ≤ exp(2C 3 t)G(0) + C 5 2C 3 ε p-1 (exp(2C 3 t) -1) -H(t) ,
and using again the bound on H(t) in (4.11) we find that for some constant C 6 > 0,

F ε (X N |ρ; t) 2 + t 0 ẊN (s) -ũ(t, X N (s)) 2 m 0 2 ds ≤ exp(2C 3 t)[G(0)+C 6 ε p-1 ]+C 6 ε p-1 .
4.6. Estimates on the initial datum. In order to conclude we only need to estimate G(0) and the initial energy F ε (0), since C 6 is an affine function of this latter, due to (4.10). Recall the definition of δ N in (4.5). Using Jensen's inequality and the expression for ρ 0 in (1.7), we obtain (4.14)

F ε (X N (0)) ≤ i Φ(T i ) |x -x 0 i | 2 2ε dx + i U m 0 i |Φ(T i )| |Φ(T i )| ≤ i T i |Φ(x) -x 0 i | 2 2ε ν(x) (ρ 0 • Φ)(x) dx + Ω U (ρ 0 ) = C 7 δ 2 N 2ε + Ω U (ρ 0 ) ,
where C 7 = ρ -1 min . Moreover, (4.15)

R d U (ρ 0 )d(µ N (0) -ρ 0 ) = i Φ(T i ) (U (ρ 0 (x i )) -U (ρ 0 ))ρ 0 ≤ Lip(U (ρ 0 )) i T i |x 0 i -Φ(x)|ν(x)dx ≤ Lip(U (ρ 0 )) 2 2 ρ 0 [Ω]ε + δ 2 N 2ε .
where ρ 0 [Ω] is the integral of ρ 0 over Ω. Hence, combining (4.14) and (4.15) we obtain

G(0) = F ε (X N (0)) - Ω U (ρ 0 ) - R d U (ρ 0 )d(µ N (0) -ρ 0 ) ≤ (C 7 + 1) δ 2 N 2ε + Lip(U (ρ 0 )) 2 2 ρ 0 [Ω]ε .
Combining the estimates above we finally find

(4.16) F ε (X N |ρ; t) 2 + t 0 ẊN (s) -ũ(t, X N (s)) 2 m 0 2 ds ≤ exp(2C 3 t) (C 7 + 1) δ 2 N 2ε + Lip(U (ρ 0 )) 2 2 ρ 0 [Ω]ε + C 6 ε p-1 + C 6 ε p-1 ,
where C 6 is an affine function of δ 2 N /ε, which concludes the proof of Theorem 1.1.

4.7.

External potentials. We now consider a slight modification of the original system where the total energy is given by

Z ε (X) := F ε (X) + i V (x i )m 0 i ,
where V : R d → R is a Lipschitz function. The gradient flow of this energy, i.e. the trajectories satisfying Ẋ = -∇ m 0 Z ε (X) solve the following modified system ODEs (4.17)

ẋi (t) = - |L i (t)| m 0 i x i (t) -b i (t) ε -∇V (x i (t))
In this case the limit PDE is given by (4.18)

   ∂ t ρ -div [ρ(∇U (ρ) + ∇V )] = 0 on (0, T ) × Ω (∇U (ρ) + ∇V ) • n ∂Ω = 0 on (0, T ) × ∂Ω .
The proof above also apply to this case with some minor changes. First of all, we observe that the velocity field is now u = -∇U (ρ) -∇V . We assume that this is sufficiently smooth so that the extension Lemma 4.2 applies. Then, using the same modulated energy as above, the only different term in equation (4.8) is I 4 (t) which should be replaced by

Ĩ4 (t) := ṽ(t, X N (t)) -∇V (X N (t)) -ũ(t, X N (t)), ẊN (t) -ũ(t, X N (t)) m 0
where as before ṽ = -∇U (ρ). This can be controlled exactly as above, leading to the same convergence result as in Theorem 1.1, but with ∇U (ρ) replaced by ∇U (ρ) + ∇V .

Time discretization and numerical tests

5.1. Time discretization. In order to compute numerically the solution to discrete model (1.3) on a given time interval [t 0 , T ], we will consider the same explicit time discretization used in [START_REF] Thomas O Gallouët | Convergence of a Lagrangian discretization for barotropic fluids and porous media flow[END_REF] and originally proposed by Brenier in [START_REF] Brenier | Derivation of the Euler Equations from a Caricature of Coulomb Interaction[END_REF]. Given a time step τ = |T -t 0 |/N T > 0 with N T ∈ N, define the discrete solution (X n ) N T n=0 as follows: given X 0 , compute X n+1 = (x n+1 i ) N i=1 for n ≥ 0 by (5.1)

x n+1 i = b n i + exp - |L n i | m 0 i ε τ (x n i -b n i )
where L n i and b n i are the ith cell of the optimal tessellation at the nth step and its barycenter, respectively. This scheme can be obtained by following on each time interval [nτ, (n + 1)τ ] the gradient flow of the energy

F n ε (X) = i L n i |x -x i | 2 2ε dx + i U m 0 i |L n i | |L n i | ,
where L n i is fixed. As a consequence of the definition of the discrete energy (1.1), this latter is dissipated by the discrete process defined by (5.1):

F ε (X n+1 ) ≤ F n ε (X n+1 ) ≤ F n ε (X n ) = F ε (X n ) .

Numerical tests.

In this section we present some numerical tests to verify the convergence estimates of Section 4. All the experiments correspond to the case where L N (Ω) = T s N (Ω). The computation of the energy and optimal tessellation is perfermed using Newton's method applied on the system of optimality conditions for the vector of weights w ∈ R N given in (2.6), similarly to the case of semi-discrete optimal transport described in [START_REF] Kitagawa | Convergence of a Newton algorithm for semi-discrete optimal transport[END_REF]. Computationally, this is simpler than the case of the Moreau-Yosida regularization (1.11) considered in [START_REF] Thomas O Gallouët | Convergence of a Lagrangian discretization for barotropic fluids and porous media flow[END_REF][START_REF] Leclerc | Lagrangian discretization of crowd motion and linear diffusion[END_REF], as the optimality conditions in (2.6) do not require computing integrals of nonlinear functions over the cells. The code that generated the tests in this section is available at https: //github.com/andnatale/gradient_flows_of_interacting_cells, and is based on the open-source library sd-ot, which is available at https://github.com/sd-ot. 5.2.1. Barenblatt test case. We consider the case where U (r) = r γ /(γ -1) and P (r) = r γ with γ > 1, in which case the corresponding PDE (1.5) is the porous medium equation. For this energy, we have an exact solution on R d which is given by the Barenblatt profile:

(5.2) ρ(t, x) = 1 t α C 2 - k t 2β |x| 2 1 γ-1 + where α = d d(γ -1) + 2 , β = α d , k = β(γ -1) 2γ ,
The exact flow is given by

ϕ(t, x) = t t 0 β x .
Note that this case falls outside the hypotheses of our theorem, due to lack of a positive lower bound on the density. Note also that since the solution has a compact support the choice of the domain Ω, if sufficiently large, has no impact on the results. We solve the discrete system on the interval [t 0 , T ] with t 0 = 1/16, T = 1, and C = 1/3, and using ε = 10/N and τ = 10/N 2 . The initial conditions for the particle γ = 1.5 model are defined via equations (1.7), (1.13) and (1.15), where Φ is a radial map from a reference ball of given radius (on which we set ν = 1) to the support of ρ 0 , which can be computed explicitly from (5.2), and T N is a Voronoi tessellation of the reference ball. For all tests we will monitor the weighted l 2 error on the flow, at the final time T , defined as in the last section as ∆ϕ := X N (T ) -ϕ N (T ) m 0 , where ϕ N (t) = ϕ(t, X N (0)) and where ϕ is the exact flow associated with the vector field -∇U (ρ). The results in table 1 show generally a faster convergence than that predicted by Theorem 1.1 but confirm a dependence of the convergence rates on the growth rate of the internal energy function U . 5.2.2. Quadratic potential. We consider again the internal energy function U (r) = r 2 , but with an additional quadratic potential V (x) = |x-x| 2 /2, driving the particles towards x ∈ R d . As described in Section 4.7, the discrete model is now defined by the system of ODEs (4.17), and we can apply the same time discretization strategy described in Section 5.1, which leads to the scheme:

x n+1 i = c n i + exp (-λ n i τ ) (x n i -c n i )
, where

λ n i = |L n i | m 0 i ε + 1 and c n i = b n i + x -b n i λ n i .
In this case the density in the continuous model (4.18) converges exponentially towards the Barenblatt profile ρ ∞ (x) = max(( M 2π ) 1/2 -1 4 |x -x| 2 , 0) where M is the total mass. Here, we consider as initial condition a configuration where the particles are equally spaced within a cross of unit height and width, with barycenter at x = 0, and share the same mass, m 0 i = M/N for all i. In particular we set M = 0.12, N = 1.23 2 and3 show the particle distribution at different times and the energy evolution, respectively, and show the exponential decay of the density towards the equilibrium distribution. for all functions ψ ∈ C(Ω). Problem (A.1) always admits at least a solution γ, and we call this optimal transport plan from ρ to µ. Semi-discrete optimal transport refers to the case one of the two measures is discrete and the other is absolutely continuous. In this case, it can be that (A.1) admits a dual formulation which can be expressed in terms of Laguerre tessellations. Suppose that µ = i m i δ x i where X = (x i ) i ∈ (R d ) N and m i > 0, and that ρ ∈ M + (Ω) is absolutely continuous and satisfies ρ[Ω] = i m i , then (A.2) W 2 2 (ρ, µ) = max

• 10 4 , τ = 1/3 • 10 -2 , ε = 2/3 • 10 -2 . Figures
w∈R N i L i (X,w) |x -x i | 2 -w i dρ(x) + w i m i .
The maximum is always attained and the maximizer w ∈ R N is related to the optimal plan γ by Ω×Ω ψ(x, y)dγ(x, y) = i L i (X,w) ψ(x, x i )dρ(x) ∀ ψ ∈ C(Ω × Ω) .

A.2. Energy reformulation. Let us show the equivalence between (1.1) and (1.10). Suppose that X ∈ R d \ ∆ N . Then, by definition of the W 2 distance,

(A.3) F ε (X) ≥ inf a∈R N >0 ,η∈C W 2 2 i a i δ x i , η 2ε + i C i (a i ) =: E ε (X)
where F ε (X) is given by (1.1), and where C = {dx} if L N (Ω) = T N (Ω), and C = {f dx : f : Ω → [0, 1]} if L N (Ω) = T s N (Ω). Note that the minimization on the right hand side of (A.3) is implicitely taken under the constraint i a i = η[Ω], and we will keep using this convention in the following. Using the dual formulation (A.2) and exchanging inf and sup we find

E ε (X) ≥ sup w∈R N inf a∈R N >0 ,η∈C i L i (X,w) |x -x i | 2 -w i 2ε dη(x) + w i 2ε a i + C i (a i )
Optimizing over η and a, we find that the right-hand side is equal to D ε (X) = F ε (X) and thereofore F ε (X) = E ε (X).

A.3. Proof of Proposition 2.3. Let X, Y ∈ (R d ) N \ ∆ N and L ∈ L N (Ω) the optimal tessellation associated with X. Then

F ε (Y ) ≤ i L i |x -x i + x i -y i | 2 2ε + C(|L i |) ≤ F ε (X) + G ε (X), Y -X m 0 + |Ω| i |x i -y i | 2 2ε ,
where

(G ε (X)) i := |L i | m 0 i x i -b i ε , b i := 1 |L i | L i x dx .
This shows that G ε (X) ∈ ∂ + F ε (X), the Fréchet superdifferential of F ε at X. We now prove that G ε (X) is continuous, which is implies that F ε (X) is C 1 and that G ε (X) is its gradient at x with respect to the inner product •, • m 0 . For this, we will use the following expression for F ε (X) (shown in Section A.2):

(A.4) F ε (X) = min a∈R N ,η∈C W 2 2 (η, i a i δ x i ) 2ε + i C i (a i )
where C is a convex subset of M(Ω) defined as above. If X ∈ (R d ) N \ ∆ N , then problem (A.4) admits a unique solution which is linked to the solution L(X) = {L i (X)} i of problem (1.1) by

η(X) = i 1 L i (X) dx , a i (X) = |L i (X)| .
Since the function minimized in (A.4) is continuous with respect to a and η (with respect to the narrow topology) on the set {(a, η) : i a i = η[Ω]}, then the optimal η(X) and a i (X) are continuous functions of X on (R d ) N \ ∆ N . In particular, given

a sequence (X n ) n ⊂ (R d ) N \ ∆ N , such that X n = (x n i ) i → X ∈ (R d ) N \ ∆ N for n → ∞, we have η(X n ) η(X) , i a i (X n )δ x n i i a i (X)δ x i .
Denoting by γ n ∈ M(R d , R d ) the optimal transport plan from η(X n ) to i a i (X n )δ x n i , by the stability of optimal transport plans γ n γ, the optimal plan from η(X) to i a i (X)δ x . Now, for any ε > 0 and n sufficiently large we can assume that |x n i -x n j | ≥ 3ε for all i = j, and |x i -x n i | ≤ ε. Then let B(x i , δ) the closed ball of radius δ > 0 centered at x i , and consider a continuous function ψ : R d → R such that ψ(x) = 1 for x ∈ B(x i , ε), and ψ(x) = 0 for x ∈ R d \ B(x i , 2ε). Then, for n → +∞,

L i (X n ) xdx = R d ×R d xψ(y)dγ n (x, y) → R d ×R d xψ(y)dγ(x, y) = L i (X)
xdx , which shows that G ε (X) is continuous.

Lemma 4 . 1 .

 41 Let U and P be smooth functions on [0, ∞) verifying (2.1) and (4.2). Then (4.3) |P (r|s)| ≤ AU (r|s) ∀ r, s > 0 .

  and in particular the curve ρ : t ∈ [0, T ] → ρ(t)| Ω is the unique solution to the continuity equation on [0, T ] × Ω associated with u and initial conditions ρ(0) = ρ 0 ; ρ ≥ ρmin > 0, where ρmin only depend on ρ min , sup t∈[0,T ] u(t) C 2,1 , T and d; moreover, sup t∈[0,T ] ρ(t) C 1,1 only depends on ρ 0 C 1,1 , sup t∈[0,T ] u(t) C 2,1 , T , d and on ρ min .

Figure 2 .

 2 Figure 2. Scatter plot of the particle positions at different times (from left to right, t = 0, 0.05, 0.2, 8) for the quadratic potential test case. The color scale refers to the density, computed for each particle as m 0 i /|L i |.

Figure 3 .

 3 Figure 3. Energy evolution for the quadratic potential test case: (a) F ε (X); (b) i U (m 0 i /|L i |)|L i |; (c) internal energy of the equilibrium density U (ρ ∞ ).

Table 1 .

 1 Error and rate of convergence for the Barenblatt test case.

	1/	√	N	∆ϕ	rate	∆ϕ	γ = 2	rate	∆ϕ	γ = 4	rate
	1.00e-01 4.54e-02		8.23e-02		2.07e-01
	5.00e-02 3.81e-02 2.54e-01	5.32e-02 6.30e-01	1.21e-01 7.75e-01
	2.50e-02 3.03e-02 3.27e-01	3.40e-02 6.45e-01	6.84e-02 8.23e-01
	1.25e-02 2.28e-02 4.10e-01	2.11e-02 6.85e-01	3.88e-02 8.19e-01
	6.25e-03 1.64e-02 4.80e-01	1.28e-02 7.23e-01	2.12e-02 8.69e-01

Acknowledgements

This work was partly supported by the Labex CEMPI (ANR-11-LABX-0007-01).

Appendix A. Optimal transport tools and proof of Proposition 2.3

A.1. Optimal transport. Given two positive measures ρ, µ ∈ M + (Ω) with fixed total mass ρ[Ω] = µ[Ω], the L 2 -Wasserstein distance between ρ and µ, is defined via