Ronen Shekel

Ohad Lib

Rodrigo Gutiérrez-Cuevas

Sébastien M Popoff

Alexander Ling

Yaron Bromberg

Shaping Single Photons through Multimode Optical Fibers using Mechanical Perturbations

The capacity of information delivered by single photons is boosted by encoding high-dimensional quantum dits in their transverse shape. Transporting such high-dimensional quantum dits in optical networks may be accomplished using multimode optical fibers, which support the low-loss transmission of multiple spatial modes over existing infrastructure. However, when photons propagate through a multimode fiber their transverse shape gets scrambled because of mode mixing and modal interference. This is usually corrected using freespace spatial light modulators, inhibiting a robust all-fiber operation. In this work, we demonstrate an all-fiber approach for controlling the shape of single photons and the spatial correlations between entangled photon pairs, using carefully controlled mechanical perturbations of the fiber. We optimize these perturbations to localize the spatial distribution of a single photon or the spatial correlations of photon pairs in a single spot, enhancing the signal in the optimized spot by over an order of magnitude. Using the same approach we show a similar enhancement for coupling light from a multimode fiber into a single-mode fiber.

I. INTRODUCTION

Quantum technologies are revolutionizing the fields of communication 1 , sensing 2 , and computing 3 . Due to their low decoherence, photons are excellent candidates for transmitting quantum information between distant points and are often referred to as "flying qubits" 4,5 . Photons also enjoy the advantage of existing mature optical fiber infrastructure offering low-loss transmission. In particular, multimode fibers (MMFs) can greatly boost the information capacity for both classical 6,[START_REF] Ma | Scaling information pathways in optical fibers by topological confinement[END_REF] and quantum 8-10 channels, by encoding high-dimensional dits in the transverse shape of the photons. In the quantum regime, single photons have been delivered through MMFs for high-dimensional quantum key distribution (QKD) [START_REF] Ding | High-dimensional quantum key distribution based on multicore fiber using silicon photonic integrated circuits[END_REF] , and to increase the collection efficiency of freespace QKD [START_REF] Jin | Demonstration of analyzers for multimode photonic time-bin qubits[END_REF] . In addition, photon pairs have been sent through the same MMF for creating programmable and high-dimensional quantum circuits [START_REF] Defienne | Two-photon quantum walk in a multimode fiber[END_REF][START_REF] Leedumrongwatthanakun | Programmable linear quantum networks with a multimode fibre[END_REF][START_REF] Makowski | Large reconfigurable quantum circuits with spad arrays and multimode fibers[END_REF] .

However, when classical or quantum light propagates through a complex medium such as an MMF, the information it carries is scrambled due to the random interference between the different modes of the MMF. In particular, the spatial distribution of a single photon in a superposition of multiple fiber modes exhibits a speckle pattern at the output facet, exactly like classical coherent light. When pairs of spatially entangled photons propagate in MMFs, the spatial distribution of the intensity at the output of the fiber is homogeneous, yet the a) yaron.bromberg@mail.huji.ac.il two-photon spatial correlations exhibit a speckle pattern, coined two-photon speckle [START_REF] Beenakker | Two-photon speckle as a probe of multi-dimensional entanglement[END_REF][START_REF] Peeters | Observation of two-photon speckle patterns[END_REF] .

The speckle and two-photon speckle patterns can be refocused to a spot by tailoring the wavefront at the input of the fiber using a spatial light modulator (SLM) [START_REF] Carpenter | 110x110 optical mode transfer matrix inversion[END_REF][START_REF] Plöschner | Seeing through chaos in multimode fibres[END_REF][START_REF] Defienne | Nonclassical light manipulation in a multiple-scattering medium[END_REF][START_REF] Wolterink | Programmable two-photon quantum interference in 10 3 channels in opaque scattering media[END_REF][START_REF] Peng | Manipulation of the spontaneous parametric down-conversion process in space and frequency domains via wavefront shaping[END_REF][START_REF] Defienne | Adaptive quantum optics with spatially entangled photon pairs[END_REF][START_REF] Lib | Real-time shaping of entangled photons by classical control and feedback[END_REF][START_REF] Cao | Controlling light propagation in multimode fibers for imaging, spectroscopy and beyond[END_REF] . The wavefront can be found by measuring the transmission matrix of the fiber, from which the desired input is calculated [START_REF] Popoff | Measuring the transmission matrix in optics: an approach to the study and control of light propagation in disordered media[END_REF] , or by using a "blind" optimization scheme to maximize the signal at the focal spot [START_REF] Vellekoop | Focusing coherent light through opaque strongly scattering media[END_REF] . Both methods rely on manipulating light outside the fiber, requiring careful alignment of the shaped wavefront into the fiber, thus inhibiting a robust all-fiber operation.

Recently, we demonstrated an all-fiber approach for shaping classical light through MMFs. The transmission matrix of an MMF is highly sensitive to mechanical perturbations of the fiber 28 , as depicted in Fig. 1. Turning this nuisance into a feature, we have recently shown that by applying carefully controlled mechanical perturbations on the fiber, we could shape classical light in both the spatial and spectral domains [START_REF] Resisi | Wavefront shaping in multimode fibers by transmission matrix engineering[END_REF][START_REF] Finkelstein | Spectral shaping in a multimode fiber by all-fiber modulation[END_REF] . We refer to this system of mechanical perturbations along the fiber as the fiber piano.

In this work, we use the fiber piano to refocus the spatial distribution of a single photon and the spatial correlations of photon pairs that propagate through an MMF. In both cases, we achieve over an order of magnitude enhancement of the signal in the focal spot. We further demonstrate the ability to enhance the coupling of a single photon from an MMF into a single-mode fiber (SMF), which could potentially increase the collection efficiency in free-space quantum communications [START_REF] Jin | Demonstration of analyzers for multimode photonic time-bin qubits[END_REF][START_REF] Cao | Distribution of high-dimensional orbital angular momentum entanglement over a 1 km few-mode fiber[END_REF] . These results pave the way toward an all-fiber control of heralded FIG. 1. When a photon propagates through a multi-mode fiber, its spatial distribution is scrambled, producing a speckle pattern. This speckle is highly sensitive to mechanical perturbations of the fiber, so different mechanical fiber conformations produce different output patterns. By optimizing controlled mechanical perturbations of the fiber, we localize the spatial distribution at the output.

single photons and photon pairs through MMFs.

II. RESULTS

To demonstrate our all-fiber approach for shaping photons, we place 37 piezoelectric actuators along a graded index MMF with a diameter of 50µm. Each actuator creates a three-point contact that induces a local bend, as depicted in Fig 2 . We couple into the MMF photons generated by type-0 spontaneous parametric downconversion (SPDC) (see methods). To study the shaping of a single photon passing through the MMF, we utilize the SPDC light to generate heralded single photons: we use a tunable beamsplitter in front of the MMF to probabilistically route one of the photons to a detector (D1 in Fig. 2a), and its twin photon to the MMF. We refer to this configuration as the single photon configuration. The image plane of the output facet of the MMF is then scanned by a second detector (D2). Coincidence events between the two detectors, i.e., simultaneous detection of a photon by each of the detectors, reveal the spatial distribution of the heralded photon that passes through the MMF.

To study the correlations between pairs of photons that propagate together through the MMF, we tune the beamsplitter before the fiber to couple both photons to the MMF. We refer to this configuration as the two-photon configuration. We place two detectors (D2 and D3) at the image plane of the MMF output facet. One detector is fixed in place while the other scans the transverse plane, and the coincidence map reveals the two-photon spatial correlations. In both configurations, the collection of light to the detectors is done using fibers with a diameter of 50µm, which is slightly smaller than the size of a single speckle grain. Narrowband filters and a polarizing beam splitter are placed before the detectors FIG. 2. Experimental setup. A continuous-wave pump laser (λp = 403.8nm) is focused using a lens (L1) upon a 4mm type-0 PPKTP crystal, generating vertically polarized photon pairs via SPDC. The pump beam is discarded using a dichroic mirror (DM), and the photon pairs are imaged using lenses L2 and L3 to the input facet of a graded index MMF. Along the MMF 37 piezoelectric actuators are placed (b), inducing 3-point local bends on the fiber (c). The output facet of the MMF is then imaged using lenses L4 and L5 to the coincidence measurement plane, where there are two singlephoton detectors. The photons are split using a polarizing beam splitter (PBS), so one photon reaches a fixed detector (D3) and the other a scanning detector (D2). In the heralded single-photon configuration, we place before the fiber a beam splitter with a tunable reflectivity, implemented by a λ 2 waveplate and a PBS. The PBS probabilistically routes one photon to the heralding detector (D1) and one photon into the fiber. (M -mirror; t-BS -tunable beam splitter).

to ensure the detection of approximately a single polarization and spectral mode. For full details, see Methods and Supplementary information.

As depicted in Fig. 3, when either a single heralded photon (a) or photon pairs (f) propagate through the MMF, the coincidence maps exhibit a speckle pattern (b, g). The single counts, which are equivalent to the classical intensity, also exhibit a speckle pattern, but with a low contrast (c, h). The reason for the low contrast in the single counts is that the two-photon quantum state we generate is in a superposition of ≈ 15 spatial modes that sum incoherently [START_REF] Peeters | Observation of two-photon speckle patterns[END_REF] . For a detailed discussion on the contrast of the single counts see Supplementary information.

To refocus the spatial distribution of the heralded single photons and the spatial correlations of the photon pairs, we use a particle swarm optimization (PSO) algorithm to find a set of mechanical perturbations that maximizes the coincidence rate of the two detectors. We define the enhancement as the ratio between the peak coincidence rate after optimization to the disorder average of the coincidence counts rate at the same point without optimization. For the single photon experiment we obtain an enhancement factor of 16 ± 3, and for the twophoton experiment we get an enhancement of 12 ± 2, as depicted in Fig. 3d,i.

We note that optimization of the coincidence rate at the target position slightly enhances also the single counts at that position (3e, j), yet by a factor ten times smaller. The small enhancement in the single counts shows that most of the enhancement of the correlations

FIG. 3.
We optimize the coincidence counts of single heralded photons (a-e) and photon pairs (f-j) passing through an MMF. Before the optimization process, the coincidence maps exhibit a speckle pattern (b,g). After the optimization, the coincidence maps are strongly localized (d, i), with an enhancement factor of 16 ± 3 and 12 ± 2 for the one photon and twophoton configurations, respectively. As expected, the single counts (c, e, h, j) exhibit a low contrast and show only a slight enhancement compared with the coincidences. The collection of light to the detectors is done using fibers with a diameter of 50µm, which is slightly smaller than the size of a single speckle grain. All scale bars represent 200µm. (ph -photons).

between the photons stems from unscrambling of the twophoton speckle pattern. The reason for this is that the two-photon signal is coherent, allowing to set the phases of different modes to interfere constructively at a desired spot, and is therefore easier to enhance by wavefront shaping than incoherent signals, similar to classical wavefront shaping. A detailed discussion on the contributions of the single counts, the polarization and the loss to the enhancement is given in the Supplementary information.

Next, to demonstrate the flexibility of our approach to tailoring the shape of single photons and two-photon correlations, we perform the optimization to focus the coincidence maps on two spots. We optimize the coincidence counts between the fixed detector and each one of the two spots, simultaneously (see methods). Relative to the average coincidence counts at each spot, we achieve an enhancement of 4.1 ± 0.4 (4.7 ± 0.5) on the upper spot and 5.7 ± 0.5 (6.4 ± 0.8) on the lower spot in the single photon (two-photon) setting. The coincidence maps after the optimization process are shown in Fig. 4.

Finally, we consider the coupling of single photons from an MMF to an SMF. While sending single photons through MMFs is appealing when the information is encoded in the spatial degree of freedom, SMFs are more attractive in the case of temporal and polariza-

FIG. 4.
To demonstrate a higher degree of control using mechanical perturbations, we change the cost function to optimize the coincidence counts between the fixed detector and each of two spots. This is done both in the single photon configuration (a) and in the two-photon configuration (b). The collection of light to the detectors is done using fibers with a diameter of 100µm, which is on the order of a single speckle grain. The two scale bars represent 200µm. (ph -photons).

tion encoding. This is due to their lower dispersion and polarization mixing, better interface with superconducting detectors, and better components and infrastructure. However, in applications such as free-space quantum communications [START_REF] Bedington | Progress in satellite quantum key distribution[END_REF][START_REF] Liao | Satellite-to-ground quantum key distribution[END_REF] , coupling into SMFs limits the link efficiency, due to the need to map the distorted incoming wavefront to the single mode of the fiber. To demonstrate funneling of single photons through the MMF and coupling into an SMF, we couple heralded single photons into the MMF, and then optimize the coupling efficiency into a single polarization mode of an SMF using mechanical perturbations. We show that after optimization the collection efficiency into the SMF is enhanced by a factor of 11.3 ± 0.6. We show the optimization process for the coupling into the SMF in the Supplementary information.

III. DISCUSSION

We demonstrate the shaping of heralded single photons and photon pairs through MMFs via controlled mechanical perturbations, enhancing the signal at a focal spot by over an order of magnitude. In comparison to the standard approach of shaping the photons using SLMs, our approach is all-fiber, which can potentially offer a low-cost and robust solution. Most importantly, all-fiber shaping eliminates the need to achieve precise coupling of a large number of modes encoding the information of the quantum state into an MMF, which is particularly challenging in the high-dimensional case.

Building on our demonstration of all-fiber shaping, we believe that utilizing recent demonstrations of the generation and detection of spatially entangled photons directly inside MMFs [START_REF] Sulimany | All-fiber source and sorter for multimode correlated photons[END_REF][START_REF] Garay-Palmett | Fiber-based photon-pair generation: tutorial[END_REF] could enable all-fiber implementations of high-dimensional quantum technologies. For example, control of single photon propagation through MMFs could enable all-fiber high-dimensional QKD and entanglement distribution [START_REF] Valencia | Unscrambling entanglement through a complex medium[END_REF][START_REF] Goel | Inversedesign of high-dimensional quantum optical circuits in a complex medium[END_REF] , while control of photon pairs propagating through the same MMF could enable allfiber linear quantum networks [START_REF] Leedumrongwatthanakun | Programmable linear quantum networks with a multimode fibre[END_REF] .

Another advantage of the fiber piano is that it allows tailoring of the entire transmission matrix of an MMF that supports N guided modes, by increasing the number of controlled actuators in the experiment to approximately N 2 . This is in contrast to the degree of control that can be achieved with an SLM in front of the fiber, which can only control 2N degrees of freedom, namely the complex amplitudes of the excited guided modes. Therefore, we believe that scaling up our proof-of-concept demonstration will not only yield larger enhancements [START_REF] Resisi | Wavefront shaping in multimode fibers by transmission matrix engineering[END_REF] , but also provide an all-fiber approach for programming arbitrary transformations on multiple input modes. This is analogous to the advantage of using multiplane light converters (MPLCs) with respect to a single SLM in both classical [START_REF] Morizur | Programmable unitary spatial mode manipulation[END_REF][START_REF] Labroille | Efficient and mode selective spatial mode multi-plexer based on multi-plane light conversion[END_REF][START_REF] Fontaine | Laguerre-gaussian mode sorter[END_REF][START_REF] Būtaitė | How to build the "optical inverse" of a multimode fibre[END_REF] and quantum [START_REF] Brandt | High-dimensional quantum gates using full-field spatial modes of photons[END_REF][START_REF] Hiekkamäki | High-dimensional two-photon interference effects in spatial modes[END_REF][START_REF] Lib | Processing entangled photons in high dimensions with a programmable light converter[END_REF] experiments. We thus expect that while in this work our objective was to restore the correlations between specific modes, increasing the number of controlled actuators will allow for future applications in high-dimensional entanglement certification through MMFs [START_REF] Valencia | Unscrambling entanglement through a complex medium[END_REF][START_REF] Defienne | Two-photon quantum walk in a multimode fiber[END_REF][START_REF] Goel | Inversedesign of high-dimensional quantum optical circuits in a complex medium[END_REF] and the realization of programmable linear quantum networks [START_REF] Leedumrongwatthanakun | Programmable linear quantum networks with a multimode fibre[END_REF] .

Finally, similarly to adaptive optics with deformable mirrors which also relies on mechanical perturbations, we expect the optimization rate of our approach to be significantly improved by using a beacon laser for optimization [START_REF] Defienne | Two-photon quantum walk in a multimode fiber[END_REF][START_REF] Defienne | Adaptive quantum optics with spatially entangled photon pairs[END_REF] , fast electronics, and faster actuators, as discussed in the Supplementary information. We thus believe that our demonstration could be relevant to improving the link efficiency in free-space quantum communications through turbulence [START_REF] Zheng | Free-space to few-mode-fiber coupling under atmospheric turbulence[END_REF][START_REF] Ghalaii | Quantum communications in a moderate-to-strong turbulent space[END_REF] . This could be achieved by first coupling the distorted single photons into an MMF and then using mechanical perturbations to counter-effect the scattering in both the link and the MMF, yielding efficient coupling of the photons into an SMF.

IV. METHODS

The experimental setup is presented in Fig 2 . An L = 4mm long type-0 PPKTP crystal is pumped by a continuous-wave laser (140 mW, λ p = 403.8 nm). The wavefront of the pump beam is focused on the crystal using a lens with a focal length of L 1 = 200mm. The pump profile in the crystal plane is approximately Gaussian with a waist of w 0 ≈ 110µm. Pairs of vertically polarized photons centered at λ ≈ 807.6nm are generated by the SPDC process. Since the divergence of the pump beam (λ p /w 0 ) is much smaller than the divergence of the emitted pairs ( λ/L), the pairs are generated in an entangled superposition of multiple transverse modes [START_REF] Walborn | Spatial correlations in parametric down-conversion[END_REF] . The pump beam is separated from the entangled photons using a dichroic mirror and sent to a beam dump. The entangled photons are then imaged using two lenses with focal lengths L 2 = 200mm and L 3 = 60mm to the facet of an ≈ 1.8m long multimode bare fiber (Corning ClearCurve OM2, graded index, numerical aperture of 0.2 and diameter of 50µm). Between these lenses, the photons propagate through a λ/2 waveplate and a polarizing beam splitter (PBS), used for switching between the single photon and the two-photon configurations.

Along the fiber we place 37 piezoelectric actuators (CTS NAC2225-A01), with a 3cm separation between adjacent actuators, creating three-point contacts that induce local bends along the fiber. The range of curvatures induced by the bends was chosen such that they would induce a significant change in the speckle pattern while causing a relatively low loss of a few percent per actuator.

The far end of the fiber is imaged using a 20X objective lens (Olympus RMS20X) and another lens with focal length L 5 = 300mm, through another PBS, to the coincidence measurement plane. The coincidence measurements are performed after passing through 3nm FWHM bandpass filters (Semrock LL01-808-12.5), using 50µm multimode fibers (Thorlabs FG050LGA, step index, numerical aperture of 0.22) coupled to single photon detectors (Excelitas SPCM-AQRH). The coincidence counts are measured using a time tagger (Swabian, Time tagger 20), with a coincidence window of 1ns. All coincidence counts shown in this paper are after subtracting accidental counts.

In the single photon configuration, the heralding photon is collected between L2 and L3, using a 50µm multimode fiber (Thorlabs FG050LGA). When we optimize for two spots we use two fibers out of a fiber bundle with three 100µm diameter cores, with 200µm separation between core centers (FiberTech Optica FTO-CTFOLA). In the two-photon configuration when optimizing for two spots, we used also for the second fiber a 100µm multimode fiber. To enhance coincidence rates at the two spots, c 1 and c 2 , such that c 1 ≈ c 2 , we maximize the cost function

√ c 1 + √ c 2 -α • |c 1 -c 2 |
, where α ≈ 0.04 was adjusted to balance the tradeoff between a high signal and an equal signal at both spots.

When coupling heralded single photons from an MMF into an SMF, we connect the fibers using a standard connector, and pass through an inline fiber polarizer (Thorlabs ILP780PM-FC) before arriving at the detectors.

The Schmidt number K quantifies the entanglement of a quantum state. In the case of a maximally entangled state in the spatial domain, it is related to the number of different incoherent spatial modes that a single photon occupies when tracing out its twin photon. After propagating through the MMF, each of these modes produces a different speckle pattern. When looking at the single counts, these modes are summed incoherently, reducing the image contrast to zero for high Schmidt numbers. For a finite number of K modes, we can expect to see a speckle pattern with a reduced contrast C which is given by C

= 1 √ K 1 .
To get a good coupling efficiency of the SPDC light into the multimode fiber, the Schmidt number must be lower than the ≈ 300 modes per polarization that the fiber supports. This guides us in choosing the parameters of the pump beam on the crystal, and here we quantify the Schmidt number of the state created by the SPDC process.

To estimate the Schmidt number we measure the single counts at a given spot at the output of the fiber, for ≈ 1900 different random actuator configurations. We calculate the standard deviation of the measured single counts, which is equivalent to the contrast C of the speckle pattern at the measured spot. From the contrast we compute the Schmidt number K = 1 C 2 . However, besides the spatial modes, there are other factors that can reduce contrast, such as the photons' bandwidth, the spatial resolution in the acquisition, and polarization effects. Assuming these factors affect the contrasts of the single counts and coincidence in the same manner, the ratio between the number of estimated modes in the single counts and in the coincidence counts can help us get a better estimate of the number of spatial modes. This is due to the fact that for maximally entangled states the contrast of the two-photon speckle pattern is 1 2 . a) yaron.bromberg@mail.huji.ac.il

We calculate ≈ 1.2 modes in the coincidence counts, and arrive at an estimate of K ≈ 15 SPDC spatial modes.

B. Polarization and spectral modes

Using a PBS before the detectors ensures we choose a single polarization mode, where one photon is horizontally polarized, and the other is vertically polarized. This has the same effect as inserting a polarizer and then a beam splitter before the detectors, as it selects a single polarization mode.

A measure for the number of spectral modes in a given fiber is the spectral correlation width ∆λ, defined as the wavelength shift required for creating two uncorrelated speckle patterns at the output of the fiber 3 . To measure ∆λ, we send light from a superluminescent diode (spectral bandwidth of ≈ 14nm) through the MMF, and collect at the output light from a single speckle grain into a spectrometer. The variation in intensity for different wavelengths originates from the spectral correlation width, which we quantify by the full-width-half-max of the autocorrelation of the measured spectrum. We estimate ∆λ ≈ 2.8nm, which is comparable to the width of the bandpass filters we use for the SPDC light. We therefore conclude that the spectral correlation width of the fiber has a minor effect on the contrast of the two-photon speckle that we measure and on the enhancements that we obtain.

S2. THE OPTIMIZATION PROCESS

A. Contributions to the enhancement

When optimizing the correlations, we showed an enhancement of over an order of magnitude. There are several mechanisms by which mechanical perturbation can contribute to the enhancement, and here we quantify them. Two global mechanisms that might contribute to the enhancement are a global polarization rotation to the selected polarization mode and a global reduction of the loss induced by the actuators. These effects enhance the total signal at the output of the fiber but do not change the shape of the pattern. We quantify these effects using a camera before and after the optimization process and find that these factors change by less than 12%.

Another mechanism that can enhance the signal at the focused spot is amplifying the single counts. The single counts are indeed enhanced by a factor of 1.7 (1.5) in the single photon (two-photon) configuration. This slight enhancement is expected due to the incoherent sum of different SPDC modes, one of which we enhance.

The main contribution to the enhancement comes from unscrambling the two-photon speckle by creating constructive interference at the focused spot. To quantify this, we normalize the enhancement by the total coincidence counts before and after the optimization. The normalized enhancement we obtain is 12 ± 2 (7 ± 1) for the single photon (two-photon) configuration. The enhancement due to all of the other mechanisms mentioned above can be quantified by the ratio between the total coincidence counts after and before optimization, which is 1.4 (1.8) for the single photon (two-photon) configuration.

B. Optimizing the single counts does not work

As mentioned above, propagation through an MMF of each SPDC spatial mode results in a different speckle pattern, and measuring the single counts will display an incoherent sum of these patterns. Conversely, measuring coincidence events effectively picks out a single mode and measures its spatial distribution. For this reason, trying to focus the single counts is equivalent to focusing all different modes simultaneously and would result in very low enhancement when looking at the coincidence counts for a specific mode.

To demonstrate this, we run the optimization using the single counts for feedback. As depicted in Fig. S1, we see that even though the single counts are slightly enhanced, the coincidence counts are not enhanced at all.

C. Optimization speed

The optimization speed is determined by the signal acquisition time and by the speed of the electronics and mechanical actuators. Since we use the coincidence counts as our signal, our acquisition time for each iteration is 5 seconds for the single photon configuration and 8 seconds for the two-photon configuration. In addition, our current electronics and mechanical actuators allow changing configurations every ≈ 0.3 seconds. This results in a slow optimization, which typically takes a few tens of minutes. However, using a beacon laser could reduce the acquisition time by many orders of magnitude 4,5 . In addition, faster electronics and mechanical actuators, such as those used in deformable mirrors, can result in rates approaching the KHz regime 6 , nearly three orders of magnitude faster than our current system. Furthermore, since our approach to modulating the transmission matrix of the fiber is general and not limited to mechanical perturbations, it could be directly transferred to other types of actuators, e.g., in-fiber acousto-optical or electro-optical modulators, which could reach even higher rates.

D. Optimization process for coupling into an SMF

We show the optimization process of coupling heralded single photons from an MMF into an SMF. As depicted in Fig. S2, within a few hundred iterations (corresponding to ≈ 20 minutes), the coupling into the SMF is enhanced by over an order of magnitude.

FIG. S2. We optimize the coupling of heralded single photons from an MMF into an SMF. The green dashed line and envelope correspond to the mean and standard deviation of the counts for the first 100 configurations, which are random. Within a few hundred iterations we achieve an enhancement of over an order of magnitude.

FIG. S1 .

 S1 FIG. S1. We run the optimization process using the single counts as feedback. While the single counts are enhanced (a), this is not sufficient for enhancing the correlations (b). Both scale bars represent 200µm. (ph -photons).

J. W. Goodman, Speckle phenomena in optics: theory and applications (Roberts and Company Publishers

, 2007). 2 C. Beenakker, J. Venderbos, and M. Van Exter, "Two-photon speckle as a probe of multi-dimensional entanglement," Physical review letters 102, 193601 (2009).

FUNDING INFORMATION

This research was supported by the Zuckerman STEM Leadership Program and the ISF-NRF Singapore joint research program (Grant No. 3538/20).

DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

DISCLOSURES

The authors declare no conflicts of interest.


Paola Goatin
email: paola.goatin@inria.fr

Macroscopic traffic flow modelling: from kinematic waves to autonomous vehicles

Keywords: Macroscopic traffic models, hyperbolic systems of conservation laws, finite volume schemes, PDE-ODE systems. AMS subject classification: 35L65, 90B20, 82B21

This survey offers a unified up-to-date presentation of macroscopic models of traffic flow, pointing out their main characteristics and possible drawbacks. The presentation is completed by several pictures illustrating the models' features. Some open problems and future research directions are also given to inspire the reader.

Introduction

Macroscopic traffic flow models consist of partial differential equations derived from fluid dynamics to describe the spatio-temporal evolution of locally aggregated quantities such as the traffic density ρ and the mean velocity v. Starting from the fundamental assumption that the number of cars on a road segment is conserved, one easily recovers a mass conservation law in one space dimension, in the form [START_REF] Lighthill | On kinematic waves. II. A theory of traffic flow on long crowded roads[END_REF] ∂ t ρ + ∂ x q = 0, x ∈ R, t > 0, where ρ = ρ(t, x) ≥ 0 is number of vehicles per space unit and q = q(t, x) ≥ 0 is the flow, expressed as number of vehicles per time unit. The fundamental hydrodynamic relation q = ρv allows to rewrite (1) as

(2) , where Vmax = 100 km/h and Rmax = 150 veh/km, corresponding to the initial datum ρ 0 (x) = Rmaxχ([0.1, 0.6]) and describing queue dissolution at a traffic light located at x = 0.6 and turning to green at t = 0. Left: (t, x) plot. Right: density profile at t = 0.003 h. We can observe the upstream shock and the downstream rarefaction, which start interacting after t = 0.005 h.

∂ t ρ + ∂ x (ρv) = 0, x ∈ R, t > 0, ρ v e 0 R max V max ρ q e 0 ρ cr R max ρ v e ρ cr 0 R max V max ρ q e 0 ρ cr R max
LWR model (3) allows to capture some basic features of traffic on road networks, such as congestion formation and propagation and queue dissolution, see Figure 2 for an example of solution to (3) consisting of an upstream shock interacting with a downstream rarefaction wave. Nevertheless, the assumed functional dependency of the average velocity on the traffic density, which holds at equilibrium, can be far from real behaviour in some situations. Indeed, it assumes that drivers can adjust instantaneously their velocities according to the densities they are experiencing, thus implying infinite acceleration of the vehicles. To overcome these limitations, second order models consist of an additional dynamical acceleration equation for v, thus resulting in 2 × 2 systems of balance laws, which capture traffic instabilities, capacity drop phenomena and scattered flow-density data, see e.g. [START_REF] Haut | A second order model of road junctions in fluid models of traffic networks[END_REF][START_REF] Seibold | Constructing set-valued fundamental diagrams from jamiton solutions in second order traffic models[END_REF]. The second order model proposed by Payne [START_REF] Payne | Models of freeway traffic and control[END_REF] and Whitham [START_REF] Whitham | Linear and nonlinear waves[END_REF] suffered of some drawbacks, such as the presence of negative velocities under certain circumtances [START_REF] Daganzo | Requiem for second-order fluid approximations of traffic flow[END_REF]. An alternative model was then proposed by Aw, Rascle [START_REF] Aw | Resurrection of "second order" models of traffic flow[END_REF] and Zhang [START_REF] Zhang | A non-equilibrium traffic model devoid of gas-like behavior[END_REF], whose non-conservative formulation reads as

(4)    ∂ t ρ + ∂ x (ρv) = 0, ∂ t v + (v -ρp (ρ))∂ x v = v e (ρ) -v τ , x ∈ R, t > 0,
for some function p satisfying p(ρ) > 0, p (ρ) > 0 and 2p (ρ) + ρp (ρ) > 0 for ρ > 0. This pseudo-pressure accounts for drivers' anticipation of downstream density changes. The source term expresses the willing of drivers to adapt their speed to the equilibrium one, τ ≥ 0 being the relaxation time parameter. We observe that the conservative form of (4) writes ( 5)

   ∂ t ρ + ∂ x (ρv) = 0, ∂ t (ρ(v + p(ρ))) + ∂ x (ρv(v + p(ρ))) = ρ v e (ρ) -v τ , x ∈ R, t > 0,
see Figure 3 for the corresponding speed-density and flux-density fundamental diagrams. Figure 4 shows a solution corresponding to the homogeneous version of ( 5), where we have set the right-hand side to zero (letting τ → +∞). Here, we model the presence of more aggressive drivers willing to run faster, thus causing a congestion upstream the group of slower cars preceding them and resulting in a more complex dynamics than the one in Figure 2, which cannot be captured by the LWR model [START_REF] Haut | A second order model of road junctions in fluid models of traffic networks[END_REF].

ρ v 0 R max V max V min ρ ρv 0 R max
Several extensions of the above mentioned models are possible. First of all, model parameters can be space and/or time dependent, to describe heterogeneous road conditions and variable speed limits. Besides, multi-lane dynamics can be handled by systems of balance laws with suitable lane-changing terms [START_REF] Klar | A hierarchy of models for multilane vehicular traffic. I. Modeling[END_REF][START_REF] Colombo | Well posedness for multilane traffic models[END_REF][START_REF] Holden | Models for dense multilane vehicular traffic[END_REF]. Moreover, macroscopic models can be extended to describe traffic flow on road networks by introducing suitable coupling conditions at junction nodes, see e.g. [START_REF] Garavello | Models for vehicular traffic on networks[END_REF] for an overview. Also, non-local dependencies in the flux function have been proposed recently, aiming at modeling the reaction of drivers to downstream traffic conditions [START_REF] Blandin | Well-posedness of a conservation law with non-local flux arising in traffic flow modeling[END_REF][START_REF] Friedrich | A Godunov type scheme for a class of LWR traffic flow models with non-local flux[END_REF]. These extensions will not be addressed in the present review and we invite the interested reader to refer to the corresponding literature.

Macroscopic traffic flow models

The principal macroscopic traffic models mentioned above can be cast in the following class of 2 × 2 hyperbolic systems of conservation laws, which is referred to as Generic Second Order traffic flow Models ). This situation describes queue dissolution at a traffic light located at x = 0.6 and turning to green at t = 0, with a group of more aggressive divers initially located behind slower ones, which causes a congestion formation behind the downstream group, where density almost reaches value 200. Left: (t, x) plot. Right: density profile at t = 0.001 h. We can observe the upstream shock and the downstream rarefaction, separated by an higher density region delimited by an upstream shock and a downstream contact discontinuity separating the two groups of drivers.

(GSOM) [START_REF] Lebacque | Generic second order traffic flow modelling[END_REF]:

(6) ∂ t ρ + ∂ x (ρv) = 0, ∂ t (ρw) + ∂ x (ρwv) = 0, x ∈ R, t > 0,
where the average speed of vehicles is a function of the density ρ = ρ(t, x) and a Lagrangian vehicle property w = w(t, x), namely v = V(ρ, w) for some speed function V satisfying the following hypotheses [START_REF] Fan | Comparative model accuracy of a data-fitted generalized Aw-Rascle-Zhang model[END_REF]:

V(ρ, w) ≥ 0, V(0, w) = w, (7a)
2V ρ (ρ, w) + ρV ρρ (ρ, w) < 0 for w > 0, (7b)

V w (ρ, w) > 0, (7c) ∀w > 0 ∃ R(w) > 0 : V(R(w), w) = 0. (7d)
As in [START_REF] Fan | Comparative model accuracy of a data-fitted generalized Aw-Rascle-Zhang model[END_REF][START_REF] Chiarello | Micro-macro limit of a nonlocal generalized Aw-Rascle type model[END_REF], we observe that (7b) implies that Q(ρ, w) := ρV(ρ, w) is strictly concave and V ρ (ρ, w) < 0 for w > 0, if V is a C 2 function in ρ. For later use, we denote by ρ cr (w) ∈ ]0, R(w)[ the point of maximum of Q(ρ, w) and by Q max (w) = Q(ρ cr (w), w) the corresponding maximum. We also remark that in (7d) we can have R(w) = R max for all w > 0.

System (6) is strictly hyperbolic for ρ > 0, with eigenvalues

(8) λ 1 (ρ, w) = V(ρ, w) + ρV ρ (ρ, w), λ 2 (ρ, w) = V(ρ, w),
and corresponding eigenvectors

(9) r 1 (ρ, w) = -1 0 , r 2 (ρ, w) = V w (ρ, w) -V ρ (ρ, w) ,
with the first characteristic field being genuinely non-linear and the second linearly degenerate [19, Definition 5.2]. The associated Riemann invariants are

z 1 (ρ, w) = V(ρ, w), z 2 (ρ, w) = w.
Since shock and rarefaction curves of each family coincide, the system belongs to the Temple class [START_REF] Temple | Systems of conservation laws with invariant submanifolds[END_REF]. Notice that, setting V(ρ, w) = w -p(ρ) for a suitable "pressure" function p, system (6) corresponds to the Aw-Rascle-Zhang (ARZ) model [START_REF] Payne | Models of freeway traffic and control[END_REF]. We also remark that, taking w = w constant, we recover the classical Lighthill-Whitham-Richards (LWR) model ( 3) with v e (ρ) = V(ρ, w). We consider for [START_REF] Whitham | Linear and nonlinear waves[END_REF] an invariant domain of the form ( 10)

Ω := (ρ, w) ∈ R 2 : ρ ∈ [0, R(w max )], w ∈ [w min , w max ]
for some 0 < w min ≤ w max < +∞. Since V ρ (ρ, w) < 0 and V(0, w) = w, the range of v = V(ρ, w) is given by v ∈ [0, w] for any w ∈ [w min , w max ]. In particular, we have v = w in the vacuum. See Figure 5 for a representation of Ω in different coordinates. We recall that, for Temple class systems, general well-posedness results [START_REF] Baiti | The semigroup generated by a Temple class system with large data[END_REF][START_REF] Bianchini | Stability of L ∞ solutions for hyperbolic systems with coinciding shocks and rarefactions[END_REF][START_REF] Bressan | Stability of L ∞ solutions of Temple class systems[END_REF] (for initial data with possibly unbounded total variation) hold under the assumption of strict hyperbolicity. For system (6), we observe that

λ 1 (0, w) = V(0, w) = w = λ 2 (0, w),
showing that strict hyperbolicity is lost at vacuum. While existence of solutions can still be proved for data which are BV in the Riemann invariants [START_REF] Andreianov | A second-order model for vehicular traffics with local point constraints on the flow[END_REF][START_REF] Godvik | Existence of solutions for the Aw-Rascle traffic flow model with vacuum[END_REF], uniqueness does not hold if vacuum is involved. Let us consider for example a Riemann-like initial datum of the form

(11) (ρ 0 , w 0 )(x) = (ρ L , w L ) if x < 0, (0, w R ) if x > 0, with ρ L > 0 and w L ≥ w R > 0 such that V(ρ L , w L ) < V(0, ρ R ).
One solution consists of a first family rarefaction wave followed by a contact discontinuity

(ρ 1 , w 1 )(t, x) =            (ρ L , w L ) if x < λ 1 (ρ L , w L )t, λ -1 1 (x/t) if λ 1 (ρ L , w L )t ≤ x < λ 1 (ρ M , w L )t, (ρ M , w L ) if λ 1 (ρ M , w L )t ≤ x < w R t, (0, w R ) if x ≥ w R t, where (ρ M , w L ) is an intermediate state satisfying V(ρ M , w L ) = w R .
A second one can be constructed as a first family rarefaction wave connecting the left state directly to the vacuum, followed by a vacuum wave (not visible, since it has no physical meaning), as proposed in [START_REF] Andreianov | A second-order model for vehicular traffics with local point constraints on the flow[END_REF]:

(ρ 2 , w 2 )(t, x) =      (ρ L , w L ) if x < λ 1 (ρ L , w L )t, λ -1 1 (x/t) if λ 1 (ρ L , w L )t ≤ x < λ 1 (0, w L )t, (0, w L ) if x ≥ λ 1 (0, w L )t.
The two solutions are illustrated in Figure 6. Notice that they are both weak entropy solutions in the sense of Lax [START_REF] Lax | Hyperbolic systems of conservation laws[END_REF].

v w U R U L U M 0 w max w min w max w min ρ ρv U R = 0 U L U M R(w max ) R(w min ) Figure 6
. Two entropy admissible solutions to the Riemann problem ( 6)- [START_REF] Colombo | Well posedness for multilane traffic models[END_REF] in two different coordinate systems: (ρ 1 , w 1 ) is depicted in dotter red, (ρ 2 , w 2 ) in dashed blue.

The intrinsic instability at vacuum of second order models, which poses issues for their analytical and numerical treatment, has been addressed introducing the so called phase-transition models [START_REF] Colombo | Hyperbolic phase transitions in traffic flow[END_REF][START_REF] Goatin | The Aw-Rascle vehicular traffic flow model with phase transitions[END_REF]. The idea is to couple second order models describing the congested dynamics with the first order LWR equation ( 3) accounting for free flow conditions, by introducing a transition dynamics between free and congested flows. In particular, free flow and congested phases can be distinguished in model ( 6) by taking speed laws of the form v = min{V max , V(ρ, w)}, for some 0 < V max < w min as in [START_REF] Blandin | A general phase transition model for vehicular traffic[END_REF][START_REF] Colombo | A 2-phase traffic model based on a speed bound[END_REF], or

v = v e (ρ) if 0 ≤ ρ ≤ ρ f , V(ρ, w) if ρ f < ρ ≤ R max , for some ρ f ∈ ]0, R max [ such that V(ρ f , w) = v e (ρ f
) for all w ∈ [w min , w max ] to ensure continuity, as proposed by [START_REF] Fan | A collapsed generalized Aw-Rascle-Zhang model and its model accuracy[END_REF]. Such models turn out to be well posed for general BV data [START_REF] Colombo | Global well posedness of traffic flow models with phase transitions[END_REF].

Numerical approximations

The most widely used numerical scheme for traffic flow macroscopic simulations is the finite volume Godunov scheme [START_REF] Godunov | A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics[END_REF] in its Cell Transmission Model (CTM) version [START_REF] Daganzo | The cell transmission model: A dynamic representation of highway traffic consistent with the hydrodynamic theory[END_REF], where the fluxes across interfaces are given by the minimum of the sending capacity of the upstream cell and the receiving capacity at the downstream one.

Given a space step ∆x and a time step ∆t satisfying the CFL condition max

(ρ,w)∈Ω {|λ 1 (ρ, w)|, |λ 2 (ρ, w)|}∆t ≤ ∆x, let x j+1/2 =
j∆x, j ∈ Z, be the cells interfaces, and t n = n∆t, n ∈ N, the time mesh.

Denoting by U = (ρ, y) T = (ρ, ρw) T the vector of the conservative variables (where we set y = ρw), we construct a finite volume approximate solution of (6) of the form

U ∆x = (ρ ∆x , y ∆x ) T with ρ ∆x (t, x) = ρ n j and y ∆x (t, x) = y n j for (t, x) ∈ C n j = [t n , t n+1 [ ×[x j-1/2 , x j+1/2 [.
To this end, we approximate the initial data with piece-wise constant functions

ρ 0 j = 1 ∆x x j+1/2 x j-1/2 ρ 0 (x) dx, y 0 j = 1 ∆x x j+1/2
x j-1/2 y 0 (x) dx, ∀j ∈ Z, and we iterate in time according to the conservation formulas

ρ n+1 j = ρ n j - ∆t ∆x F ρ,n j+1/2 -F ρ,n j-1/2 , (12a) y n+1 j = y n j - ∆t ∆x F y,n j+1/2 -F y,n j-1/2 , (12b)
where F ρ,n j+1/2 and F y,n j+1/2 are the the flow respectively of ρ and y at x = x j+1/2 in the time interval [t n , t n+1 [. Let us remark that, since the variable w is advected with ρv, we get

(13) F y,n j+1/2 = w n j F ρ,n j+1/2 .
in (12b). To compute F ρ,n j+1/2 , we define the corresponding demand and supply functions [START_REF] Lebacque | The Godunov scheme and what it means for first order traffic flow models[END_REF][START_REF] Lebacque | Second order traffic flow modeling: supply-demand analysis of the inhomogeneous Riemann problem and of boundary conditions[END_REF] as

D j = D(ρ j , w j ) = Q(ρ j , w j ) if ρ j ≤ ρ cr (w j ), Q max (w j ) if ρ j > ρ cr (w j ), S j+1 = S(ρ j+1 , w j+1 ; w j ) = Q max (w j ) if ρ j+1/2 ≤ ρ cr (w j ), Q(ρ j+1/2 , w j ) if ρ j+1/2 > ρ cr (w j ),
where ρ j+1/2 is the density of the intermediate state in the solution of the Riemann problem corresponding to U j and U j+1 :

ρ j+1/2 such that V(ρ j+1/2 , w j ) = V(ρ j+1 , w j+1 ) if V(ρ j+1 , w j+1 ) ≤ w j ; 0 if V(ρ j+1 , w j+1 ) > w j .
We can thus set ( 14)

F ρ,n j+1/2 = min D(ρ n j , w n j ), S(ρ n j+1 , w n j+1 ; w n j )
in [START_REF] Holden | Models for dense multilane vehicular traffic[END_REF], [START_REF] Garavello | Models for vehicular traffic on networks[END_REF]. See [START_REF] Fan | A collapsed generalized Aw-Rascle-Zhang model and its model accuracy[END_REF] for a physical interpretation of demand and supply functions, and the role of the intermediate state.

Remark 2.1. Taking w 0 j = w constant, the scheme ( 12)-( 14) boils down to the standard CTM approximation of the LWR model (3):

ρ n+1 j = ρ n j - ∆t ∆x F n j+1/2 -F n j-1/2 , with F n j+1/2 = min D(ρ n j ), S(ρ n j+1
) and

D(ρ j ) = Q(ρ j , w) if ρ j ≤ ρ cr ( w), Q max ( w) if ρ j > ρ cr ( w), S(ρ j+1 ) = Q max ( w) if ρ j+1 ≤ ρ cr ( w), Q(ρ j+1 , w) if ρ j+1 > ρ cr ( w).

Macroscopic models for heterogeneous traffic flows

The macroscopic models presented above are based on the assumption that vehicles' characteristics are homogeneous, and traffic flow is mainly impacted by drivers' behaviour (modeled by the Lagrangian marker w). In real situations, the vehicles present on a road can differ a lot by technical characteristics, like size and speed. Think for example of cars and trucks: beside their size, trucks generally run at lower speed with respect to cars (because of the imposed speed limits, but also due to their lower acceleration). Also, the foreseen deployment of autonomous vehicles on public roads is expected to impact the overall traffic conditions due to their different dynamics. For these reasons, real traffic flows can be better described accounting for possible vehicle heterogeneities.

Multi-class models

When the number of vehicles in each class is of the same order, it is still meaningful to represent each class dynamics at a macroscopic level. Denoting by ρ c , c = 1, . . . , N , the density of the c-th class, the mass conservation law holding for each class leads to the N × N system of equations ( 15)

∂ t ρ c + ∂ x (ρ c v c ) = 0, c = 1, . . . , N,
coupled by the speed functions v c , which model each class behaviour depending on the different class densities:

(16) v c = v c (ρ 1 , . . . , ρ N ).
Different coupling laws [START_REF] Lebacque | Generic second order traffic flow modelling[END_REF] are proposed in the literature, depending on the targeted application: interaction dynamics can vary greatly if we consider mixed flows of cars and trucks obeying lane discipline [START_REF] Wong | A multi-class traffic flow model-an extension of LWR model with heterogeneous drivers[END_REF][START_REF] Zhang | Kinematic wave traffic flow model for mixed traffic[END_REF][START_REF] Benzoni-Gavage | An n-populations model for traffic flow[END_REF][START_REF] Chanut | Macroscopic model and its numerical solution for two-flow mixed traffic with different speeds and lengths[END_REF], or a disordered flow of cars and two-wheelers without any lane rule [START_REF] Nair | A porous flow approach to modeling heterogeneous traffic in disordered systems[END_REF][START_REF] Fan | A heterogeneous multiclass traffic flow model with creeping[END_REF][START_REF] Gashaw | Modeling and analysis of mixed flow of cars and powered two wheelers[END_REF]. Indeed, in the latter case, we may observe creeping behaviour of cycles, stemming form the fact that jam densities (for which v c = 0) are different for cars and cycles, the latter still being able to advance when cars are stopped. Taking v c (ρ 1 , . . . , ρ N ) = V c max Ψ(r) as in [START_REF] Benzoni-Gavage | An n-populations model for traffic flow[END_REF], where r = N c=1 ρ c /R c max is the total (normalized) traffic density, V 1 max > . . . > V N max are the classes' maximal speeds and Ψ : [0, 1] → [0, 1] is a decreasing function satisfying Ψ(0) = 1 and Ψ(1) = 0, the multi-class model ( 15) is defined on the simplex

S := (ρ 1 , . . . , ρ N ) ∈ R N : N i=1 ρ i ≤ 1, ρ i ≥ 0 for i = 1, . . . , N . System ∂ t ρ c + ∂ x (V c max ρ c Ψ(r)) = 0, c = 1, . . . , N,
is hyperbolic, but lacks of strict hyperbolicity at umbilic points on the boundaries of S (with ρ i = 0 for some i = 1, . . . , N ), hindering its analytical study. Nevertheless, numerical simulations show that the model can display observed behaviour, like overtaking of slower vehicles by faster ones, see Figure 7. In ) This situation describes queue dissolution at a traffic light located at x = 0.6 and turning to green at t = 0, with a group of more aggressive divers initially located behind slower ones, which causes a congestion formation behind the downstream group, where density reaches higher levels than the initial ones. Left: (t, x) plot of the total density (ρ 1 + ρ 2 ). Right: total density profile at t = 0.001 h. We can observe the upstream shock and the downstream rarefaction, separated by an higher density region. The global behaviour is similar to the one obtained using the ARZ model in Figure 4.

this case, instead of using the multi-class equivalent of Godunov flux, as proposed by [START_REF] Fan | A heterogeneous multiclass traffic flow model with creeping[END_REF] (and which displays oscillations due to unbounded total variation), we propose the following extension of the upwind scheme first introduced by [START_REF] Hilliges | A phenomenological model for dynamic traffic flow in networks[END_REF][START_REF] Bürger | A family of numerical schemes for kinematic flows with discontinuous flux[END_REF]

ρ c,n+1 j = ρ c,n 1,j - ∆t ∆x F c,n j+1/2 -F c,n j-1/2 , c = 1, . . . , N, with F c,n j+1/2 = V c max ρ c,n 1,j Ψ(r n j+1 ), c = 1, . . . , N.
Remark 3.1. The invariance of S for N > 2 and general existence results for N ≥ 2 are currently open problems.

Multi-scale models

If a class consists of a small number of vehicles, compared to the bulk traffic, it can be modeled microscopically by a system of Ordinary Differential Equations (ODEs) describing each vehicle's trajectory and coupled with the partial differential equations (PDEs) giving the macroscopic traffic evolution. Besides other approaches [START_REF] Lattanzio | Moving bottlenecks in car traffic flow: a PDE-ODE coupled model[END_REF][START_REF] Borsche | Mixed systems: ODEs -balance laws[END_REF], the capacity reduction due to the presence of a slow vehicle can be accounted by a moving flux inequality constraints, as proposed by [START_REF] Lebacque | Introducing buses into first-order macroscopic traffic flow models[END_REF][START_REF] Delle Monache | Scalar conservation laws with moving constraints arising in traffic flow modeling: an existence result[END_REF][START_REF] Garavello | A multiscale model for traffic regulation via autonomous vehicles[END_REF]. More precisely, equations ( 6) are coupled with the following laws

ξ(t) = min{V (t), v(t, ξ(t)+)}, (17a) lim x→ξ(t)± ρ(t, x)(v(t, x) -ξ(t)) ≤ F α ( ξ(t)), (17b) where V (t) ∈ W 1,1 loc (R + ; R + )
is the desired speed of the slow vehicle, and x = ξ(t) its position at time t > 0. The upper bound in (17b) accounts for the road capacity reduction due to the presence of the vehicle, which acts as a moving bottleneck: α ∈ ]0, 1[ is the reduction rate, and F α is the corresponding maximal flow attainable at the vehicle's position:

F α ( ξ(t)) := ρ 2 α V(ρ α , w α ),
where w α ∈ ]w min , w max [ satisfies R(w α ) = αR(w max ) and ρ α ∈ ]0, αR(w max )[ is the point of maximum of the relative flow function ψ(ρ) = ρ V(ρ, w α ) -ξ(t) , i.e.

V(ρ α , w α ) + ρ α V ρ (ρ α , w α ) = ξ(t).
Classical solutions of (6) remain admissible provided that the relative flux at x = ξ(t) does not exceed the upper bound F α ( ξ(t)). Otherwise, the enforcement of constraint (17b) causes the formation of a undercompressive shock, that is a non-classical jump discontinuity which satisfies the Rankine-Hugoniot conditions, thus guaranteeing conservation, but violates Lax entropy conditions. More precisely, let

U L = (ρ L , w L ), U R = (ρ R , w R )
∈ Ω and consider the Cauchy problem for ( 6)-( 17) corresponding to the initial data

(ρ, w)(0, x) = (ρ L , w L ) if x < 0, (ρ R , w R ) if x > 0, (18a) ξ(0) = 0, (18b)
and to the constant bottleneck speed ξ(t) = V for all t > 0. Let I(w L ) be the set

I(w L ) = {ρ ∈ [0, R(w max )] : ρV(ρ, w L ) = F α ( V ) + ρ V }.
Since the map ρ → ρ V(ρ, w L )) -V is strictly concave due to (7b), the set I(w L ) contains at most two elements. If I(w L ) = ∅, let ρ = ρ(w L ) and ρ = ρ(w L ) be the points defined by These are respectively the points with maximal and minimal density of the Lax curve of the first family passing through (ρ L , w L ) which satisfy the condition (17b) on the flux, see Figure 8. Let now R GSOM be the standard Riemann solver for [START_REF] Whitham | Linear and nonlinear waves[END_REF], see e.g. [START_REF] Fan | Data-fitted generic second order macroscopic traffic flow models[END_REF] for a detailed description, and let

ρ(w L ) = max I(w L ), ρ(w L ) = min I(w L ). 0 F α ( V ) ρ α ρ(w L ) ρ(w L ) ρ ρv R(w L ) R(w α )
ρ(U L , U R )(•) and w(U L , U R )(•)
be respectively the ρ and w components of R GSOM (U L , U R )(•), and let

f 1 (R GSOM (U L , U R )(•)) := ρ(U L , U R )(•) V (ρ(U L , U R )(•), w(U L , U R )(•))
be the first component of the flux function of the GSOM system.

An admissible Riemann solver R α GSOM for the constrained GSOM system ( 6)-( 17) can be constructed as follows.

1. If f 1 (R GSOM (U L , U R )( V )) > F α ( V ) + V ρ(U L , U R )( V ), then R α GSOM (U L , U R )(x/t) = R GSOM ((ρ L , w L ), (ρ, w L ))(x/t) if x < y(t), R GSOM ((ρ, w L ), (ρ R , w R ))(x/t) if x > y(t),
and y(t) = V t, see Figure 9.

0 ρ ρ ρ ρv R(w L ) U R U L U M U L U M U R (ρ, w L )(ρ, w L ) x t 0 Figure 9.
The constrained Riemann problem ( 6)-( 17), case 1. In this case ξ(t) = ξ(0) + V t, while the solution (ρ, w) is composed by two classical shocks of the first family separated by an undercompressive shock between (ρ, w L ) and (ρ, w L ), all followed by a contact discontinuity. The unconstrained solution would have consisted of a classical shock from U L to U M = (ρ M , w L ) with V(ρ M , w L ) = V(ρ R , w R ), followed by a contact discontinuity from U M to U R .

If

f 1 (R GSOM (U L , U R )( V )) ≤ F α ( V ) + V ρ(U L , U R )( V ) then R α GSOM (U L , U R )(x/t) = R GSOM (U L , U R )(x/t) and y(t) = min V , V ρ(U L , U R )( V ), w(U L , U R )( V ) t.
Remark that the solution R α GSOM is conservative for both the traffic density and the momentum. An alternative solution conserving only the number of vehicles is proposed in [START_REF] Villa | Moving bottlenecks for the Aw-Rascle-Zhang traffic flow model[END_REF]. Moreover, in case 1., the solution given by R α GSOM does not satisfy the Lax entropy condition at the jump discontinuity between the left state (ρ, w L ) and the right state (ρ, w L ), because ρ < ρ and thus λ 1 (ρ, w L ) > λ 1 (ρ, w L ). Therefore, (ρ, w L ) and (ρ, w L ) are connected by a non-classical shock.

In the case w min = w max , the constrained GSOM system ( 6)-( 17) boils down to the constrained LWR model described in [START_REF] Delle Monache | Scalar conservation laws with moving constraints arising in traffic flow modeling: an existence result[END_REF], see Figure 11. In this case, existence and stability of solutions satisfying a suitable entropy condition has been proved in [START_REF] Garavello | A multiscale model for traffic regulation via autonomous vehicles[END_REF][START_REF] Liard | Well-posedness for scalar conservation laws with moving flux constraints[END_REF]. Also, the above construction can be extended to several, possibly interacting, moving bottlenecks, see [START_REF] Goatin | Interacting moving bottlenecks in traffic flow[END_REF].

Non-classical shocks arising at the bottleneck positions cannot be captured by classical (entropic) finite volume schemes, and they require a specific treatment, see e.g. [START_REF] Daganzo | On the numerical treatment of moving bottlenecks[END_REF][START_REF] Delle Monache | A numerical scheme for moving bottlenecks in traffic flow[END_REF][START_REF] Čičić | Front tracking transition system model with controlled moving bottlenecks and probabilistic traffic breakdowns[END_REF]. The conservative reconstruction strategy proposed in [START_REF] Villa | Moving bottlenecks for the Aw-Rascle-Zhang traffic flow model[END_REF][START_REF] Chalons | A conservative scheme for non-classical solutions to a strongly coupled PDE-ODE problem[END_REF] consists of using the classical Godunov scheme ( 12)-( 14) away from bottleneck locations. Let C n m the cell where the moving bottleneck is located at time t n , i.e. ξ n ξ(t n ) ∈ [x m-1/2 , x m+1/2 [, and we set U n j = (ρ n j , w n j ). If

f 1 (R GSOM (U n m-1 , U n m+1 )( V )) > F α ( V ) + V ρ(U n m-1 , U n m+1 )( V ),
a non-classical shock is expected to arise in cell C n m , and the left and right numerical fluxes F ρ,n m±1/2 and F y,n m±1/2 are to be recomputed as follows. We replace the (constant) cell values ρ n m and y n m by the functions

ρ n rec = ρ(w n m-1 )1 [x m-1/2 ,x ρ m [ + ρ(w n m-1 )1 [x ρ m ,x m+1/2 [ , y n rec = ŷ(w n m-1 )1 [x m-1/2 ,x y m [ + y(w n m-1 )1 [x y m ,x m+1/2 [ , where xρ m = x m-1/2 + h d ρ,n m and xy m = x m-1/2 + h d y,n m with (19) d ρ,n m = ρ n m -ρ(w n m-1 ) ρ(w n m-1 ) -ρ(w n m-1 ) ∈ [0, 1] and d y,n m = y n m -y(w n m-1 ) ŷ(w n m-1 ) -y(w n m-1 ) ∈ [0, 1],
to guarantee conservation in both the ρ and y components, see Figure 10. Above, we have set ŷ(w n m-1 ) = ρ(w n m-1 ) w n m-1 and y(w n m-1 ) = ρ(w n m-1 ) w n m-1 . Since the non-classical shock travels with speed V > 0, we denote ∆t ρ m+1/2 and ∆t y m+1/2 respectively the time needed by the ρ and the y component of the discontinuity to reach the interface x m+1/2 :

x m-1/2 x m+1/2 ρ x xρ m ρ n m-1 ρ(w n m-1 ) ρ n m ρ(w n m-1 ) ρ n m+1 x m-1/2 x m+1/2
∆t ρ m+1/2 = ∆x 1 -d ρ,n m V , ∆t y m+1/2 = ∆x 1 -d y,n m V .
Therefore, we replace the classical Godunov flux F ρ,n m+1/2 by

F ρ,n m+1/2 = 1 ∆t min(∆t ρ m+1/2 , ∆t)ρ(w n m-1 )V(ρ(w n m-1 ), w n m-1 )
+ max(∆t -∆t ρ m+1/2 , 0)ρ(w n m-1 )V(ρ(w n m-1 ), w n m-1 ) , and F y,n m+1/2 by F y,n m+1/2 = w n m-1 F ρ,n m+1/2 . Also, we replace F ρ,n m-1/2 by

F ρ,n m-1/2 = min D(ρ n m-1 , w n m-1 ), S(ρ(w n m-1 ), w n m-1 ; w n m-1 )
and F y,n m-1/2 by F y,n m-1/2 = w n m-1 F ρ,n m-1/2 . Figure 11 displays de solution computed with the above scheme in the case w min = w max (and hence w constant), where we have applied the reconstruction technique also to the classical shock. Observe that both jump discontinuities are sharply captured.

Conclusions and perspectives

Macroscopic models can be used for estimation and prediction, but also as a basis to design control actions intended to optimize traffic flow on road networks, by improving throughput, reducing travel times and fuel consumption. Traditional optimization measures include variable message signs (such as traffic adaptive speed limits and dynamic routing), traffic-light control at intersections and ramp metering, or selective measures for specific vehicle classes (such as overtaking or circulation bans for trucks). We refer the reader to [START_REF] Treiber | Traffic Flow Dynamics[END_REF][START_REF] Ferrara | Freeway Traffic Modelling and Control[END_REF] for an overview of these approaches.

More recently, the expected deployment of Connected and Automated Vehicles (CAVs) has captured the attention of researchers as it offers the possibility to use CAVs as endogenous controls acting as moving bottlenecks, thus reducing traffic speed upstream. The benefits induced by CAV control have already been demonstrated experimentally [START_REF] Stern | Dissipation of stop-and-go waves via control of autonomous vehicles: Field experiments[END_REF][START_REF] Stern | Quantifying air quality benefits resulting from few autonomous vehicles stabilizing traffic[END_REF] and make the object of model-based theoretical studies [START_REF] Čičić | Traffic regulation via individually controlled automated vehicles: a cell transmission model approach[END_REF][START_REF] Piacentini | Traffic control via moving bottleneck of coordinated vehicles[END_REF][START_REF] Čičić | Stop-and-go wave dissipation using accumulated controlled moving bottlenecks in multi-class CTM framework[END_REF][START_REF] Piacentini | Traffic control via platoons of intelligent vehicles for saving fuel consumption in freeway systems[END_REF] as well as machine learning investigations [START_REF] Vinitsky | Lagrangian control through Deep-RL: Applications to bottleneck decongestion[END_REF][START_REF] Kreidieh | Dissipating stop-and-go waves in closed and open networks via deep reinforcement learning[END_REF]. In particular, the multi-scale models presented in Section 3.2 can be used to quantify the gain that can be obtained by controlling the desired speed of small fleets of CAVs in mixed traffic flows. As an example, we illustrate in Figure 12 the result of total fuel consumption minimization algorithms using a centralized control of a small fleet of CAVs [START_REF] Daini | Centralized traffic control via small fleets of connected and automated vehicles[END_REF], supporting the perspective of regulating traffic flow by means of a small number of controlled vehicles. ) veh/h respectively. The cost function to be minimized is the Total Fuel Consumption (TFC) as described in [START_REF] Piacentini | Traffic control via moving bottleneck of coordinated vehicles[END_REF]. We can observe on the left the non controlled situation, in the middle the solution corresponding to the optimal constant speed values on 1 hour time interval, on the right the solution given a Model Predictive Control algorithm corresponding to 15 min optimization updated every 5 min. The corresponding gain in TFC is 5.6% and 3.2% respectively, but we can also observe a consistent reduction of the congestion (in red).

Figure 1 .

 1 Figure 1. Speed-density (left) and flow-density (right) fundamental diagrams corresponding to two common choices of the speed function. Top: ve(ρ) = Vmax 1 -ρ Rmax (Greenshields' fundamental diagram). Bottom: ve(ρ) = min Vmax, w Rmax ρ -1 (triangular fundamental diagram).

Figure 2 .

 2 Figure 2. Solution of (3) with ve(ρ) = Vmax 1 -ρ Rmax , where Vmax = 100 km/h and Rmax = 150 veh/km, corresponding

Figure 3 .

 3 Figure 3. Speed-density (left) and flow-density (right) fundamental diagrams for the ARZ model (5) corresponding to the choice p(ρ) = w -αρ, with w ∈ [V min , Vmax] and α > 0. We observe that the model allows for zero flux (and zero average speed) for densities belonging to the interval 1/α • [V min , Vmax], and not for a sigle (maximal) density value. This may model different types of more or less aggressive drivers, jamming at higher or lower densities when traffic is stopped.

Figure 4 .

 4 Figure 4. Solution of (5) with τ = +∞ and p(ρ) = ρ, corresponding to the initial datum ρ 0 (x) = 150χ([0.1, 0.6]) and v 0 (x) = 50χ([0.1, 0.4]) + 10χ(]0.4, 0.6]). This situation describes queue dissolution at a traffic light located at x = 0.6 and turning to green at t = 0, with a group of more aggressive divers initially located behind slower ones, which causes a congestion formation behind the downstream group, where density almost reaches value 200. Left: (t, x) plot. Right: density profile at t = 0.001 h. We can observe the upstream shock and the downstream rarefaction, separated by an higher density region delimited by an upstream shock and a downstream contact discontinuity separating the two groups of drivers.

Figure 5 .

 5 Figure 5. Invariant domain Ω (10) in the (v, w)-coordinates (left) and the (ρ, ρv)-coordinates (rigtht). Vacuum states are highlighted in red.

Figure 7 .

 7 Figure 7. Solution of (15) with V 1 max = 140 km/h, V 2 max = 100 km/h and R 1 max = R 2 max = 200 veh/km, corresponding to the initial data ρ 0 1 (x) = 0.75Rmaxχ([0.1, 0.4]) and ρ 0 2 (x) = 0.75Rmaxχ(]0.4, 0.6]) This situation describes queue dissolution at a traffic light located at x = 0.6 and turning to green at t = 0, with a group of more aggressive divers initially located behind slower ones, which causes a congestion formation behind the downstream group, where density reaches higher levels than the initial ones. Left: (t, x) plot of the total density (ρ 1 + ρ 2 ). Right: total density profile at t = 0.001 h. We can observe the upstream shock and the downstream rarefaction, separated by an higher density region. The global behaviour is similar to the one obtained using the ARZ model in Figure4.

Figure 8 .

 8 Figure 8. The definition of ρ(w L ) and ρ(w L ).

Figure 10 .

 10 Figure 10. Non-classical shock numerical reconstruction.

Figure 11 .

 11 Figure 11. Solution of (6)-(17), with V(ρ, w) = Vmax 1 -ρ Rmax , where Vmax = 100 km/h and Rmax = 150 veh/km, corresponding to the initial datum ρ 0 (x) = Rmaxχ([0.1, 0.6]) and describing queue dissolution at a traffic light located at x = 0.6 and turning to green at t = 0. The presence of a slower vehicle at ξ(0) = 0.6, running at V (t) = 30 km/h and reducing by half the capacity of the road (α = 0.5), hinders the traffic flow upstream of it (compare with Figure 2). Left: (t, x) density plot and slow vehicle trajectory (in blue). Right: density profile at t = 0.003 h. We can observe the nonclassical shock at the slow vehicle position.

Figure 12 .

 12 Figure 12. Traffic flow optimization on a 50 km road stretch by centralized control of a fleet of 10 CAVs initially equally distributed in the interval [2, 20]. Here we have used model (6)-(17), with V(ρ, w) = Vmax 1 -ρ Rmax , where Vmax = 140 km/h and Rmax = 400 veh/km, and inflow and outflow boundary conditions are given by f in (t) = 1400χ([0, 0.5]) veh/h and f out (t) = 700χ([0, 1]) veh/h respectively. The cost function to be minimized is the Total Fuel Consumption (TFC) as described in[START_REF] Piacentini | Traffic control via moving bottleneck of coordinated vehicles[END_REF]. We can observe on the left the non controlled situation, in the middle the solution corresponding to the optimal constant speed values on 1 hour time interval, on the right the solution given a Model Predictive Control algorithm corresponding to 15 min optimization updated every 5 min. The corresponding gain in TFC is 5.6% and 3.2% respectively, but we can also observe a consistent reduction of the congestion (in red).


Paola Goatin
email: paola.goatin@inria.fr

Chiara Daini
email: chiara.daini@inria.fr

Maria Laura

Delle Monache
email: mldellemonache@berkeley.edu

Antonella Ferrara
email: antonella.ferrara@unipv.it

INTERACTING MOVING BOTTLENECKS IN TRAFFIC FLOW

Keywords: 2020 Mathematics Subject Classification. 90B20, 35L65 Conservation laws, PDE-ODE systems, macroscopic traffic flow models, moving bottlenecks, autonomous vehicles

We present a general multi-scale approach for modeling the interaction of controlled autonomous vehicles (AVs) with the surrounding traffic flow. The model consists of a scalar conservation law for the bulk traffic, coupled with ordinary differential equations describing the possibly interacting AV trajectories. The coupling is realized through flux constraints at the moving bottleneck positions, inducing the formation of non-classical jump discontinuities in the traffic density. In turn, AVs are forced to adapt their speed to the downstream traffic average velocity in congested situations. We analyze the model solutions in a Riemann-type setting, and propose an adapted finite volume scheme to compute approximate solutions for general initial data.

The work paves the way to the study of general optimal control strategies for AV velocities, aiming at improving the overall traffic flow by reducing congestion phenomena and the associated externalities.

1. Introduction. Moving bottleneck models have been introduced in the engineering literature starting from the end of the last century to describe the interaction of slow moving vehicles, such as buses or trucks, with bulk traffic, see [START_REF] Claudel | Lax-Hopf Based Incorporation of Internal Boundary Conditions Into Hamilton-Jacobi Equation. Part II: Computational Methods[END_REF][START_REF] Leclercq | Moving bottlenecks in Lighthill-Whitham-Richards model: A unified theory[END_REF][START_REF] Lebacque | Introducing buses into first-order macroscopic traffic flow models[END_REF]. Further approaches have then be proposed in the applied mathematics community [START_REF] Borsche | Mixed systems: ODEs -balance laws[END_REF][START_REF] Delle Monache | Scalar conservation laws with moving constraints arising in traffic flow modeling: an existence result[END_REF][START_REF] Lattanzio | Moving bottlenecks in car traffic flow: a PDE-ODE coupled model[END_REF][START_REF] Sylla | Influence of a slow moving vehicle on traffic: well-posedness and approximation for a mildly nonlocal model[END_REF][START_REF] Villa | Moving bottlenecks for the Aw-Rascle-Zhang traffic flow model[END_REF], always following the underlying idea of using multi-scale models consisting of Partial Differential Equations (PDEs) to model the overall traffic flow evolution, coupled with Ordinary Differential Equations (ODEs) accounting for the moving bottlenecks trajectories. Despite the intrinsic difficulties arising as a result of the flux constraint induced by the presence of the slow vehicles, corresponding existence and stability results have been provided, see also [START_REF] Delle Monache | Stability estimates for scalar conservation laws with moving flux constraints[END_REF][START_REF] Garavello | A multiscale model for traffic regulation via autonomous vehicles[END_REF][START_REF] Liard | Well-posedness for scalar conservation laws with moving flux constraints[END_REF][START_REF] Liard | On entropic solutions to conservation laws coupled with moving bottlenecks[END_REF], together with numerical schemes to compute approximate solutions capturing the emerging non-classical phenomena [START_REF] Bretti | Two algorithms for a fully coupled and consistently macroscopic PDE-ODE system modeling a moving bottleneck on a road[END_REF][START_REF] Chalons | A conservative scheme for non-classical solutions to a strongly coupled PDE-ODE problem[END_REF][START_REF] Daganzo | On the numerical treatment of moving bottlenecks[END_REF][START_REF] Delle Monache | A front tracking method for a strongly coupled PDE-ODE system with moving density constraints in traffic flow[END_REF][START_REF] Delle Monache | A numerical scheme for moving bottlenecks in traffic flow[END_REF].

More recently, moving bottleneck models have been applied to traffic flow regulation by means of Autonomous Vehicles (AVs), see e.g. [START_REF] Cičić | Coordinating vehicle platoons for highway bottleneck decongestion and throughput improvement[END_REF][START_REF] Cičić | Traffic regulation via individually controlled automated vehicles: a cell transmission model approach[END_REF][START_REF] Cičić | Energy-optimal platoon catch-up in moving bottleneck framework[END_REF][START_REF] Cičić | Stop-and-go wave dissipation using accumulated controlled moving bottlenecks in multi-class CTM framework[END_REF][START_REF] Cičić | Front tracking transition system model with controlled moving bottlenecks and probabilistic traffic breakdowns[END_REF][START_REF] Piacentini | Traffic control via moving bottleneck of coordinated vehicles[END_REF]. Indeed, AVs can act as endogenous actuators of moving variable speed limits, improving throughput by hindering upstream traffic. This possibility has already been the object of several theoretical studies (see e.g. [START_REF] Davis | Effect of adaptive cruise control systems on traffic flow[END_REF][START_REF] Guériau | How to assess the benefits of connected vehicles? a simulation framework for the design of cooperative traffic management strategies[END_REF][START_REF] Knorr | Reducing traffic jams via VANETs[END_REF][START_REF] Piacentini | VACS equipped vehicles for congestion dissipation in multi-class ctm framework[END_REF][START_REF] Talebpour | Influence of connected and autonomous vehicles on traffic flow stability and throughput[END_REF][START_REF] Wang | Cooperative car-following control: Distributed algorithm and impact on moving jam features[END_REF] and [START_REF] Ferrara | Freeway traffic modelling and control[END_REF]Chapter 11]), but also artificial [START_REF] Jang | Simulation to scaled city: Zero-shot policy transfer for traffic control via autonomous vehicles[END_REF] and real world [START_REF] Stern | Dissipation of stop-and-go waves via control of autonomous vehicles: Field experiments[END_REF] experiments. In particular, it has been demonstrated that even low AV penetration rates can have a significant impact on reducing congestion phenomena and the associated energy consumption. In this context, macroscopic approaches are needed to address the curse of dimensionality in control design for microscopic models [START_REF] Delle Monache | Feedback Control Algorithms for the Dissipation of Traffic Waves with Autonomous Vehicles[END_REF].

Even if some of the above works account for the presence of multiple AVs on the road, very few of them allow for AV interactions, and only in very specific modeling frameworks [START_REF] Gasser | Vehicular traffic flow dynamics on a bus route[END_REF][START_REF] Simoni | A fast simulation algorithm for multiple moving bottlenecks and applications in urban freight traffic management[END_REF]. In this paper, we aim at generalizing the control framework designed in [START_REF] Garavello | A multiscale model for traffic regulation via autonomous vehicles[END_REF][START_REF] Piacentini | Traffic control via moving bottleneck of coordinated vehicles[END_REF] for the Lighthill-Whitham-Richards model [START_REF] Lighthill | On kinematic waves. II. A theory of traffic flow on long crowded roads[END_REF][START_REF] Richards | Shock waves on the highway[END_REF] to multiple, possibly interacting, moving bottlenecks, as if they were distributed on several lanes and they were allowed both to merge (if on the same lane) and to overtake (if on different lanes). This can be achieved by a careful study of the corresponding Riemann-type problems, and a natural generalization of the numerical reconstruction technique presented in [START_REF] Chalons | A conservative scheme for non-classical solutions to a strongly coupled PDE-ODE problem[END_REF]. The proposed model can be used to develop a general simulation environment accounting for several AVs interacting with the surrounding traffic flow and among themselves, which in turns allows to design optimization-based control strategies having fleets of AVs as actuators.

The paper is organized as follows: Section 2 details the model framework, providing a mathematical description of the fully coupled PDE-ODE model and the interactions among AVs at the level of Riemann problems. Section 3 describes the numerical scheme and Section 4 shows the results of the numerical experiments.

2. PDE-ODE model with interacting moving bottlenecks. Following [START_REF] Delle Monache | Scalar conservation laws with moving constraints arising in traffic flow modeling: an existence result[END_REF][START_REF] Garavello | A multiscale model for traffic regulation via autonomous vehicles[END_REF], we consider the following strongly coupled PDE-ODEs system the maximal vehicle density attainable on the considered road section). Above, ρ 0 ∈ [0, R] and y 0 ∈ R, = 1, . . . , N , are respectively the initial traffic density and AV positions, while the function F α in (1c), α ∈ ]0, 1[, represents the road capacity reduction rate due to the presence of the AV at x = y (t), = 1, . . . , N , acting as a moving bottleneck which imposes a unilateral flux constraint at its position.

∂ t ρ (t, x) + ∂ x f (ρ (t, x)) = 0, , (1a)
ẏ (t) = min{u (t), v(ρ(t, y (t)+))}, (1b) f (ρ (t, y (t))) -ẏ (t)ρ (t, y (t)) ≤ F α ( ẏ (t)) := max ρ∈[0,R] (αf (ρ/α) -ρ ẏ (t)) , (1c) ρ(0, x) = ρ 0 (x), (1d)
y (0) = y 0 , (1e)
For simplicity, we assume that this parameter is the same for all AVs, and it is a function of the number of lanes: α = (M -1)/M , where M ∈ N denotes the number of lanes. (The extension to a more general framework would follow the same lines detailed in this work.) To determine the function F α , we consider the rescaled flux function f α : [0, αR] -→ R + ρ -→ ρv(ρ/α) = αf (ρ/α), which is strictly concave function and such that f α (0) = f α (αR) = 0. For every u ∈ [0, V ], we define the point ρu as the unique solution to the equation f α (ρ) = u. Moreover, for every u ∈ [0, V ], we define the function

ϕ u : [0, R] -→ R + ρ -→ f α (ρ u ) + u (ρ -ρu ) .
Hence, if ẏ(t) = u, the function F α in (1c) is defined by

F α : [0, V ] -→ R + u -→ ϕ u (0) = f α (ρ u ) -uρ u .
If ẏ(t) = v(ρ(t, y(t)+)), the inequality (1c) is trivially satisfied since the left-hand side is zero. Finally, the points 0 [START_REF] Delle Monache | Scalar conservation laws with moving constraints arising in traffic flow modeling: an existence result[END_REF][START_REF] Garavello | A multiscale model for traffic regulation via autonomous vehicles[END_REF] and Figure 1. We remark that ρV = ρV = ρV = ρ * V = 0. Solutions to (1) are defined as follows (see [START_REF] Garavello | A multiscale model for traffic regulation via autonomous vehicles[END_REF]Definition 3.1]):

≤ ρu ≤ ρu ≤ ρu ≤ ρ * u ≤ R are uniquely defined by ρu = min I u , ρu = max I u , I u = {ρ ∈ [0, R] : f (ρ) = ϕ u (ρ)} , v(ρ * u ) = u, see
Definition 2.1. The N + 1-tuple (ρ, y 1 , . . . , y N ) provides a solution to (1) if the following conditions hold.

1.

ρ ∈ C 0 R + ; L 1 loc (R; [0, R]) and TV (ρ(t)) < +∞ for all t ∈ R + ; 2. y ∈ W 1,1 loc (R + ; R) for = 1, . . . , N ; 3. For every κ ∈ R and for all ϕ ∈ C 1 c (R 2 ; R + ) it holds R + R (|ρ -κ|∂ t ϕ + sgn(ρ -κ)(f (ρ) -f (κ))∂ x ϕ) dx dt + R |ρ 0 -κ|ϕ(0, x) dx + 2 N =1 R + (f (κ) -ẏ (t)κ -min{f (κ) -ẏ (t)κ, F α ( ẏ (t))}) ϕ(t, y(t)) dt ≥ 0 ; 4. For a.e. t > 0, f (ρ (t, y (t)±))-ẏ (t)ρ (t, y (t)±) ≤ F α ( ẏ (t)) for = 1, . . . , N ; 5. For a.e. t > 0, ẏ (t) = min {u (t), v (ρ (t, y (t)+))} for = 1, . . . , N .
As long as AV trajectories do not intersect, one can apply the theory and the numerical schemes developed in [START_REF] Chalons | A conservative scheme for non-classical solutions to a strongly coupled PDE-ODE problem[END_REF][START_REF] Delle Monache | A front tracking method for a strongly coupled PDE-ODE system with moving density constraints in traffic flow[END_REF][START_REF] Delle Monache | Scalar conservation laws with moving constraints arising in traffic flow modeling: an existence result[END_REF][START_REF] Delle Monache | A numerical scheme for moving bottlenecks in traffic flow[END_REF][START_REF] Garavello | A multiscale model for traffic regulation via autonomous vehicles[END_REF]. In particular, we recall that the Constrained Riemann Solver is defined as follows.

Definition 2.2. The Constrained Riemann Solver R α for

∂ t ρ (t, x) + ∂ x f (ρ (t, x)) = 0, (2a) ẏ(t) = min{u, v(ρ(t, y(t)+))}, (2b)
f (ρ (t, y(t))) -ẏ(t)ρ (t, y(t)) ≤ F α ( ẏ(t)) , (2c)
ρ(0, x) = ρ L if x < 0, ρ R if x > 0, (2d)
y(0) = 0, (2e)
is the map

R α : [0, R] 2 × [0, V ] → L 1 loc (R, [0, R]) defined as follows. 1. If f (R(ρ L , ρ R )(u)) > F α (u) + u R(ρ L , ρ R )(u), then R α (ρ L , ρ R ; u)(x/t) = R(ρ L , ρu )(x/t) if x < u t, R(ρ u , ρ R )(x/t) if x ≥ u t, and
y(t) = u t. 2. If f (R(ρ L , ρ R )(u)) ≤ F α (u) + u R(ρ L , ρ R )(u), then R α (ρ L , ρ R ; u) = R(ρ L , ρ R ) and y(t) = min{u, v(ρ R )} t.
Above, R denotes the standard (i.e. without the constraint (2c)) Riemann solver for (2a)-(2d), i.e. the (right continuous) map (t, x) → R(ρ L , ρ R )(x/t) given by the standard weak entropy solution to (2a)-(2d).

In this paper, we are interested in describing what happens when two vehicles AV i and AV j , with i, j ∈ {1, . . . , N }, i = j, interact. Let us assume y i (t) < y j (t) with constant speeds u i (t) = u i > u j (t) = u j for t ∈ ] t-, ε t[, so that y i ( t) = y j ( t) for some t > 0, see Figure 2. Depending on the surrounding traffic density conditions, only three situations may occur in a Riemann-like setting (piece-wise constant density with at most one jump discontinuity):

• The constraint (1c) is enforced for = i (we will say that AV i is active), so that ρ(t, x) = ρui for x < y i (t) and ρ(t, x) = ρui for x > y i (t). Since, at x = y j (t), it holds f (ρ ui ) -u j ρui ≤ F α (u j ) (see Figure 2), the constraint is not enforced for = j (AV j is then inactive). • If AV j is active, we have ρ(t, x) = ρuj for x < y j (t) and ρ(t, x) = ρuj for

x > y j (t). At x = y i (t), it therefore holds f ρuj -

u i ρuj ≤ F α (u i ) , thus AV i is inactive. 0 Fα(uj) Fα(ui) ρuj ρui ρuj ρui ρ f f α R ϕ ui ϕ uj Figure 2.
Notation for interacting bottlenecks.

• If both AV i and AV j are inactive, we may have either a constant density or a classical shock coinciding with one of the AV trajectories. In any case, the density must satisfy ρ(t, x) ∈ [0, ρui ] ∪ [ρ uj , R]. We remark that AV i and AV j cannot be both active at the same time, unless other waves are present between them.

To detail the problem evolution after the interaction (at t > t), we distinguish if AVs are moving in different lanes or if they are located on the same lane.

Same lane interactions.

Let us assume that the two vehicles are in the same lane, so for t > t the upstream vehicle must adapt to the downstream vehicle speed: we will have y i (t) = y j (t) with u i (t) = u j (t) for t ≥ t. Three situations may occur, see Figure 3: • If, before the interaction, AV i was active and AV j inactive, at t > t the solution will consists of a classical shock joining ρui to ρuj , followed by a non-classical shock at x = y j (t) = y i (t) between ρuj and ρuj and a rarefaction wave from ρuj to ρui , see Figure 3a. In particular, after the interaction, both AVs will be active. • If, before the interaction, AV i was inactive and AV j active, at t > t the solution will consists only of a non-classical shock joining ρuj to ρuj , see Figure 3b. Again, after the interaction, both AVs will be active. • If, before the interaction, AV i and AV j were both inactive, at t > t the solution will consists at most of the classical shock that was already present before the interaction, see Figure 3c. After the interaction, both AVs will still be inactive. In any case, after the interaction, the two vehicles behave as AV j . Different lane interactions. If the two interacting vehicles travel in different lanes, for t > t they will follow their own trajectory and AV i will overtake AV j : we will have y i (t) ≥ y j (t) with u i (t) ≥ u j (t) for t ≥ t. Three situations may occur, see Figure 4: • If, before the interaction, AV i was active and AV j inactive, at t > t the solution will consists of a classical shock joining ρui to ρuj , followed by a non-classical shock at x = y j (t) from ρuj to ρuj , a classical shock from ρuj to ρui and a non-classical shock at x = y i (t) from ρui to ρui , see Figure 4a. In particular, after the interaction, both AVs will be active. • If, before the interaction, AV i was inactive and AV j active, at t > t the solution will consists of a non-classical shock at x = y j (t) from ρuj to ρuj , followed by a classical shock from ρuj to ρui , a non-classical shock at x = y i (t) from ρui to ρui and a classical shock from ρui to ρuj , see Figure 4b. In particular, after the interaction, both AVs will be active. • If, before the interaction, AV i and AV j were both inactive, at t > t the solution will consists at most of the classical shock that was already present before the interaction, see Figure 4c. After the interaction, both AVs will still be inactive. 3. Numerical scheme. We extend the reconstruction strategy introduced in [START_REF] Chalons | A conservative scheme for non-classical solutions to a strongly coupled PDE-ODE problem[END_REF] to include moving bottleneck interactions. The scheme is composed of two parts:

Numerical approximation of (1a)-(1c). To approximate the conservation equation (1a), we use a conservative finite volume scheme for the constrained hyperbolic PDE using a flux reconstruction technique at the constraint locations, which allows to capture the non-classical shocks sharply. Let ∆x and ∆t be the fixed space and time steps satisfying the standard Courant-Friedrichs-Lewy (CFL) condition [START_REF] Courant | On the partial difference equations of mathematical physics[END_REF]: max

ρ∈[0,R] |f (ρ)|∆t < ∆x,
and set x j-1/2 = j∆x, x j = (j + 1/2)∆x for j ∈ Z, and t n = n∆t for n ∈ N.

The initial data ρ 0 is approximated by the piece-wise constant function obtained by averaging it on the discretization cells C j = [x j-1/2 , x j+1/2 [, namely

ρ 0 j = 1 ∆x x j+1/2 x j-1/2 ρ 0 (x) dx, j ∈ Z.
Away from AV positions, equation (1a) is approximated by the standard Godunov scheme [START_REF] Godunov | A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics[END_REF], whose numerical fluxes at cell interfaces

F n j+ 1 2 = F (ρ n j , ρ n j+1
) can in this case be derived using the supply-demand formula [START_REF] Lebacque | The Godunov scheme and what it means for first order traffic flow models[END_REF]

F (ρ n j , ρ n j+1 ) = min{D(ρ n j ), S(ρ n j+1 )}, (3)
where

D(ρ) = f (min{ρ, ρ cr }), S(ρ) = f (max{ρ, ρ cr }), ρ cr = argmin ρ∈[0,R] f (ρ)
being the point of maximum of the flux function f . Let now the approximate -th AV position be y n ∈ C m for some m ∈ Z, = 1, . . . , N . If

f (R(ρ n m -1 , ρ n m +1 )(u n )) > F α (u n ) + u n R(ρ n m -1 , ρ n m +1 )(u n ), we assume a moving bottleneck at xm = x m -1/2 + d n m ∆x with d n m = ρu n -ρ n m ρu n -ρu n . If 0 ≤ d n m ≤ 1, then xm ∈ C m and we set ∆t n m = 1 -d n m u n ∆x, F n m -1 2 = F (ρ n m -1 , ρu n ), ∆tF n m + 1 2 = min (∆t n m , ∆t)f (ρ u n ) + max (∆t -∆t n m , 0)f (ρ u n ).
We can then update the density by means of the conservative formula

ρ n+1 j = ρ n j - ∆t ∆x F n j+ 1 2 -F n j-1 2 , j ∈ Z. (4)
Above, we set u n = u (t n ). If two or more AVs are located in the same cell, the moving bottlenecks are treated sequentially one after the other, starting from those which are not active (i.e. satisfy (1c)), then processing those which are active (violating (1c)).

Remark 1. For better resolution, we also apply the above reconstruction technique to classical shocks, as described in [3, Section 3.1].

Numerical approximation of (1b). To track the AV trajectories, at each time step, we update the positions y n using an explicit Euler scheme

y n+1 = y n + v(ρ n )∆t n , = 1, . . . , N.
If two AVs are in the same lane and in the same cell, and the upstream AV 1 moves faster than the downstream AV 2 , i.e. u n 1 > u n 2 , then we set

y n+1 1 = y n+1 2 = y n 2 + v(ρ n )∆t n and u n 1 = u n 2 .
(A more accurate computation of the interaction point could be implemented to improve simulation accuracy.) Remark 2. The above procedure seems not coherent with the observation that two vehicles travelling at close speeds u i ≈ u j would act as a single bottleneck with α ij = α i + α j -1. In this perspective, another option could be to reconstruct a single front moving at speed u ij = (u i + u j )/2 with α = α ij whenever two vehicles are in the same cell on two different lanes. In this case, the procedure should be the following:

If

f (R(ρ n m-1 , ρ n m+1 )(u ij )) > F αij (u ij ) + u ij R(ρ n m-1 , ρ n m+1 )(u ij ),
d n m = ρuij -ρ m ρuij -ρuij .
To reconstruct the numerical flux at the cell interface x m+ 1 2 , we compute

∆t m+ 1 2 = 1 -d n m u ij ∆x and we set ∆tF n m+ 1 2 = min{∆t m+ 1 2 , ∆t}f (ρ uij ) + max{∆t -∆t m+ 1 2 , 0}f (ρ uij ).
4. Numerical tests. In this section we illustrate the performances of the numerical scheme described in the previous Section 3 in capturing AV interactions correctly.

We first analyze the Riemann-like cases described in Section 2, to conclude with a more general test case. We consider a three lane stretch of road of length 50 km, parametrized by the interval [0, 50] with absorbing boundary conditions. We consider the quadratic flux function proposed by Greenshields [START_REF] Greenshields | A study of traffic capacity[END_REF] that assumes a linear decreasing dependence of the speed on the traffic density

f (ρ) = V ρ 1 - ρ R , (5)
where V = 140 km/h denotes the maximal speed and R = 400 veh/km the maximal (bump-to-bump) density on the road. The capacity reduction ratio is set to α = 0.6. In all simulations, we set ∆x = 0.2 and ∆t = 0.9∆x/V . For the Riemann-like cases, we consider the following initial data:

AV 1 : y 1 (0) = 7.5, u 1 = 50, AV 2 : y 2 (0) = 15, u 2 = 20, (a) ρ 0 (x) = ρu1 ≈ 210 if x < 7.5, ρu1 ≈ 47 if x > 7.5, (b) ρ 0 (x) = ρu2 ≈ 280 if x < 15, ρu2 ≈ 63 if x > 15, (c) ρ 0 (x) ≡ 20.
The results are depicted in Figures 5 and6. Except small spurious oscillations emerging from AV interactions in some cases, the solutions, and in particular nonclassical shocks, are accurately reproduced.

To give a more general example, we consider the following initial data: ρ 0 (x) ≡ 200 and

AV 1 : y 1 (0) = 2.5, u 1 = 120, AV 2 : y 2 (0) = 7.5, u 2 = 30, AV 3 : y 3 (0) = 10, u 3 = 55, AV 4 : y 4 (0) = 20, u 4 = 20, (6)
Moreover, AV 1 and AV 3 are on the same lane. Figure 7 displays the solution, accounting for the different interaction types (queuing and overtaking). We observe that AV 1 and AV 3 are initially inactive while AV 2 and AV 4 are active. Moreover, AV 1 and AV 3 overtake the preceding vehicles and finally merge. Also, AV 1 is always inactive until it merges with AV 3 , which becomes active after the interaction with AV 4 , while AV 4 is always active. On the contrary, AV 2 becomes inactive when it reaches an high traffic density region.

Conclusion.

We have presented a strongly coupled PDE-ODE model describing the interaction of a small number of autonomous vehicles with bulk traffic on a multi-lane highway. The model allows for AV queuing and overtaking, expanding the range of traffic control applications of previous works [START_REF] Cičić | Coordinating vehicle platoons for highway bottleneck decongestion and throughput improvement[END_REF][START_REF] Cičić | Traffic regulation via individually controlled automated vehicles: a cell transmission model approach[END_REF][START_REF] Cičić | Energy-optimal platoon catch-up in moving bottleneck framework[END_REF][START_REF] Cičić | Stop-and-go wave dissipation using accumulated controlled moving bottlenecks in multi-class CTM framework[END_REF][START_REF] Cičić | Front tracking transition system model with controlled moving bottlenecks and probabilistic traffic breakdowns[END_REF][START_REF] Piacentini | Traffic control via moving bottleneck of coordinated vehicles[END_REF]. A first study on traffic management opportunities offered by this framework is presented in [START_REF] Daini | Centralized Traffic Control via Small Fleets of Connected and Automated Vehicles[END_REF], showing that low penetration rates are sufficient to reach nearly optimal improvements of a selected performance index.

for = 1 ,FαFigure 1 .

 11 Figure 1. The definition of ρu , ρu , ρu and ρ * u .

Figure 3 .

 3 Figure 3. Possible configurations of AV interactions on the same lane.

Figure 4 .

 4 Figure 4. Possible configurations of AV interactions on different lanes.

( a )

 a Same lane interaction of an upstream active AV with an inactive one, causing the formation of a classical shock followed by a non-classical shock at the AVs' common position, see Fig.3a. (b) Same lane interaction of an inactive AV with an active preceding one, causing no change in the surrounding traffic density displaying only a non-classical shock at the AVs' common position, see Fig. 3b. (c) Same lane interaction of inactive AVs, resulting in no change in the surrounding traffic density, see Fig.3c.

Figure 5 .

 5 Figure 5. Numerical reconstruction of possible configurations of AV interactions on same lanes. Left: (t, x) representation. Right: density profile at given time instant.

( a )

 a Different lane interaction of an upstream active AV with an inactive one, causing the formation of two classical shocks alternated with two non-classical shocks at the AVs' positions, see Fig. 4a. (b) Different lane interaction of an inactive AV with an active preceding one, causing the formation of a classical shock separating two non-classical shocks at the AVs' positions, see Fig. 4b. (c) Different lane interaction of inactive AVs, resulting in no change in the surrounding traffic density, see Fig.4c.

Figure 6 .

 6 Figure 6. Numerical reconstruction of possible configurations of AV interactions on different lanes. Left: (t, x) representation. Right: density profile at given time instant.

Figure 7 .

 7 Figure 7. Numerical reconstruction of the solution of (1) corresponding to four interacting AVs with initial data (6). Top: (t, x) representation. Bottom: density profiles at given time instants. The first and third AV are on the same lane and merge after overtaking the others.

The second author was partially supported by ERASMUS+/KA1 "NORTH SOUTH TRAINEESHIP" academic year 2020/21.