Chao Lyu

Yann Capdeville
email: yann.capdeville@univ-nantes.fr

Gang Lu

Liang Zhao
email: zhaoliang@mail.iggcas.ac.cn

Removing the CFL stability criterion of the explicit time-domain very high degree spectral-element method with eigenvalue perturbation

Keywords: 2D, acoustic, finite element, wave propagation Manuscript Focus Area: Seismic Inversion, Seismic Modeling and Wave Propagation

Removing the CFL stability criterion of the explicit timedomain very high degree spectral-element method with eigenvalue perturbation

INTRODUCTION

The forward simulation of the acoustic and elastic wave equation is a basic and critical part for imaging at all scales from the exploration geophysics to the global Earth scale. In the past two decades, the explicit time-domain spectral-element method (SEM) [START_REF] Komatitsch | The spectral element method : an effective tool to silmulate the seismic response of 2d and 3d geological structures[END_REF][START_REF] Chaljub | Spectral element analysis in seismology[END_REF] has been widely used, especially in the academia. The adjoint methods by the two-way SEM modeling have made significant progress in recent regional and global seismic tomographies [START_REF] Tromp | Seismic tomography, adjoint methods, time reversal and banana-doughnut kernels[END_REF][START_REF] Fichtner | Full seismic waveform tomography for upper-mantle structure in the Australasian region using adjoint methods[END_REF][START_REF] Zhu | Structure of the european upper mantle revealed[END_REF][START_REF] French | Broad plumes rooted at the base of the earth's mantle beneath major hotspots[END_REF][START_REF] Fichtner | The collaborative seismic earth model: Generation 1[END_REF][START_REF] Capdeville | Elastic full waveform inversion based on the homogenization method: theoretical framework and 2-d numerical illustrations[END_REF][START_REF] Tromp | Seismic wavefield imaging of earth's interior across scales[END_REF][START_REF] Lyu | Intrinsic non-uniqueness of the acoustic full waveform inverse problem[END_REF]. However, the popular explicit time-domain SEM still suffers from two main constraints. On the one hand, each material interface has to be explicitly meshed to ensure the accuracy of the model definition and the precision of the simulation. On the other hand, only a very small time step ∆t can be used to fulfill the CFL stability condition, because the Gauss-Lobatto-Legendre (GLL) points cluster toward the edges of the elements especially for the very high degree SEM and the model with fine structures.

The main intent of the CFL stability condition is governed by the minimum ratio between the size of the grid spacing and the P-wave velocity as follows:

∆t ≤ C ∆x min V p , ( 1
)
where C is the courant constant, ∆x min is the minimum distance between any two GLL points, V p is the P-wave velocity. In this article, we use the ∆t CFL as the maximum stable time step of the secondorder Newmark scheme [START_REF] Newmark | A method of computation for structural dynamics[END_REF]. And please note that the selection of the available ∆t close to the ∆t CFL of an explicit time-domain SEM in the heterogeneous model has been actually discussed under the control of the CFL stability condition by the local element stability criteria [START_REF] Cottereau | Stability of an explicit high-order spectral element method for acoustics in heterogeneous media based on local element stability criteria[END_REF]).

The errors of numerical wave simulation come from the spatial-and time-dispersion errors, and these two kinds of dispersion errors have been proven not relevant [START_REF] Stork | Eliminating nearly all dispersion error from fd modeling and rtm with minimal cost increase[END_REF][START_REF] Koene | Eliminating time dispersion from seismic wave modeling[END_REF][START_REF] Lyu | Efficiency of the spectral element method with very high polynomial degree to solve the elastic wave equation[END_REF]. That means they accumulate independently during the numerical simulation.

Many studies have been focused on how to achieve both high efficiency and accuracy with coarser grids and larger time steps in the forward numerical simulation. The recent introduction of the non-1 2 3 4 5 6 7 8 9 The CFL's removal of SEM periodic homogenization has solved the meshing difficulty by replacing the complex models with their equivalent smooth versions [START_REF] Capdeville | 2-d non-periodic homogenization to upscale elastic media for p-sv waves[END_REF][START_REF] Cupillard | Non-periodic homogenization of 3-d elastic media for the seismic wave equation[END_REF]. The accuracy of the simulated waveforms can be ensured to some degree, which makes the meshing by SEM unrestricted by physical discontinuities anymore. This means that the SEM can flexibly reduce the spatial-dispersion error by tuning the number of elements and the internal polynomial degree N in each tensorial direction for each element. [START_REF] Lyu | Efficiency of the spectral element method with very high polynomial degree to solve the elastic wave equation[END_REF] validates the high efficiency of using the very high degree SEM (N = 12 ∼ 24) both for the required memory and the computation time, which makes SEM with very high degree attractive and competitive for solving the wave equation.

Except for the spatial-dispersion error, the selection of the time step is very critical for the simulated accuracy and computational efficiency as well. And please note that the time-dispersion error is only a function of frequency, time step, and propagation time, independent of the propagation path and the complex medium responses [START_REF] Stork | Eliminating nearly all dispersion error from fd modeling and rtm with minimal cost increase[END_REF][START_REF] Koene | Eliminating time dispersion from seismic wave modeling[END_REF]. For a model with sufficiently accurate spatial meshing, a converged waveform without spatial and time-dispersion errors can be obtained by using an extremely small ∆t but with low efficiency. While with the high computational efficiency, the time-dispersion error of using a relatively large time step ∆t ≤ ∆t CFL would also contaminate the numerical waveforms for the finite difference (FD) method [START_REF] Gao | Third-order symplectic integration method with inverse time dispersion transform for long-term simulation[END_REF] and as well as SEM especially with N = 4 [START_REF] Lyu | Efficiency of the spectral element method with very high polynomial degree to solve the elastic wave equation[END_REF] and even worse for long duration simulations owing to its accumulation over time. [START_REF] Stork | Eliminating nearly all dispersion error from fd modeling and rtm with minimal cost increase[END_REF] demonstrate that the time-dispersion error can be removed by time-varying filters and interpolation after the wave modelling. Based on this idea, [START_REF] Wang | Finite-difference time dispersion transforms for wave propagation[END_REF] propose that the forward time dispersion transform (FTDT) can predict the time-dispersion error, and the inverse time dispersion transform (ITDT) can eliminate the time-dispersion error from synthetic traces. [START_REF] Koene | Eliminating time dispersion from seismic wave modeling[END_REF] mentions that the obtained traces should not only be filtered by the ITDT method as proposed by [START_REF] Wang | Finite-difference time dispersion transforms for wave propagation[END_REF], but also the source time function needs to be filtered by the FTDT method before the simulation, which has been validated in the numerical methods including FD, pseudo-spectral (PS) method and SEM.

Their combination have been used to obtain better image by the reverse-time migration in real application. All of them introduce a great improvement in the simulated accuracy and computational efficiency with a time step ∆t ≤ ∆t CFL by the combination of the FTDT and ITDT methods.

Even so, their adopted time steps are still limited by the CFL stability condition. The spatial The CFL's removal of SEM filtering method [START_REF] Ecer | Digital filtering techniques for parallel computation of explicit schemes[END_REF][START_REF] Sarris | Extending the stability limit of the fdtd method with spatial filtering[END_REF][START_REF] Chang | A spatial filter-enabled high-resolution subgridding scheme for stable fdtd modeling of multiscale geometries[END_REF]Gao et al., 2019) and the eigenvalue-type method [START_REF] He | Explicit time-domain finite-element method stabilized for an arbitrarily large time step[END_REF][START_REF] Gaffar | An explicit and unconditionally stable fdtd method for electromagnetic analysis[END_REF][START_REF] Li | Model order reduction and stability enforcement of finite-difference time-domain The CFL's removal of SEM equations beyond the cfl limit[END_REF][START_REF] Li | Overcoming the fdtd stability limit via model order reduction and eigenvalue perturbation[END_REF][START_REF] Jin | Explicit and unconditionally stable fdtd method without eigenvalue solutions[END_REF][START_REF] Gao | Removing the stability limit of the explicit finitedifference scheme with eigenvalue perturbation[END_REF] are the two available methods to extend or remove the CFL stability limit of explicit scheme at present.

For the spatial filtering method, in fluid dynamics, [START_REF] Ecer | Digital filtering techniques for parallel computation of explicit schemes[END_REF] find that a time step ∆t > ∆t CFL would introduce unstable high-wavenumber wavefields. Therefore, the larger time steps of an explicit scheme could be used if one can filter out the unstable high-wavenumber wavefields. When numerically solving Maxwell's equations by a time step beyond the CFL condition, [START_REF] Sarris | Extending the stability limit of the fdtd method with spatial filtering[END_REF] adapt a low-pass filter to remove the unstable high-wavenumbers components. In the seismic field, Gao et al. (2019) extend the stability limit of explicit scheme with spatial filtering for solving the wave equations, and the FTDT and ITDT methods are then used to reduce the large time-dispersion error when using a time step larger than the ∆t CFL .

For the eigenvalue-type method, in the electromagnetic wave numerical simulation field, [START_REF] He | Explicit time-domain finite-element method stabilized for an arbitrarily large time step[END_REF] propose an explicit time-domain finite-element method (FEM) with unconditional stability, and this method eliminates the dependence of the time step size on the spatial interval. [START_REF] Gaffar | An explicit and unconditionally stable fdtd method for electromagnetic analysis[END_REF] address that the source of instability in an explicit time-domain FD using a time step ∆t > ∆t CFL is that the unstable eigenmodes whose eigenvalues are larger than what can be accurately simulated. [START_REF] Jin | Explicit and unconditionally stable fdtd method without eigenvalue solutions[END_REF] propose that these unstable modes can be discarded before performing an explicit time marching by using partial eigenvalues solution of the global system matrix. As an alternative, [START_REF] Li | Model order reduction and stability enforcement of finite-difference time-domain The CFL's removal of SEM equations beyond the cfl limit[END_REF] and [START_REF] Li | Overcoming the fdtd stability limit via model order reduction and eigenvalue perturbation[END_REF] adopt the so-called eigenvalue perturbation method to use a time step ∆t > ∆t CFL , which normalizes the unstable eigenvalues to one to extend the CFL stability limit. In the seismic field, [START_REF] Etgen | Stability of explicit depth extrapolation through laterally-varying media: 64th Annual International Meeting[END_REF] indicates if the singular values are not less than one, especially if they are larger than one, the numerical scheme will develop unstable behavior. [START_REF] Gao | Removing the stability limit of the explicit finitedifference scheme with eigenvalue perturbation[END_REF] also illustrates that for a given time step ∆t ≤ ∆t CFL , the eigenvalues of the updated matrix would be distributed along the unit circle; However, for a time step ∆t > ∆t CFL , some eigenvalues would be distributed outside of the unit circle. This is the reason of unstable phenomena of the acoustic wave numerical simulation. [START_REF] Gao | Removing the stability limit of the explicit finitedifference scheme with eigenvalue perturbation[END_REF] indicate that combining the FTDT method, the eigenvalue perturbation and the ITDT method can simulate highly accurate waveforms when applying a time step several times of the ∆t CFL in the FD framework of solving To our knowledge, the combination of FTDT, the eigenvalue perturbation and ITDT methods has received a lot of attention in the FD framework, but the related research in the SEM framework is rarely studied. No method has been proposed to overcome the CFL stability limit for the explicit time-domain SEM in the acoustic wave simulation from previously published literature.

In this paper, following the validity of the combination of the FTDT and ITDT methods to reduce the time-dispersion error in the SEM framework in [START_REF] Koene | Eliminating time dispersion from seismic wave modeling[END_REF], and the successful integration of the FTDT and ITDT methods with eigenvalue perturbation in the FD framework in [START_REF] Gao | Removing the stability limit of the explicit finitedifference scheme with eigenvalue perturbation[END_REF], we combine these three methods to remove the CFL limit in the explicit timedomain SEM simulation of solving the acoustic wave equation. We first illustrate the theoretical analyses of the combined methods in the SEM framework with the detailed introductions of the eigenvalue perturbation, FTDT and ITDT methods. Then, the associated numerical experiments are performed in the 2D homogeneous and heterogeneous cases by the SEM with time steps ∆t > ∆t CFL .

METHODOLOGY

In this section, we first briefly introduce the numerical simulation of 2D acoustic wave equation by the SEM, then illustrate how the CFL stability limit is removed by the eigenvalue perturbation and finally clarify how the resultant time-dispersion error yielded by ∆t > ∆t CFL is eliminated by the FTDT and ITDT methods.

2D acoustic wave equation

The acoustic waves propagating is governed by the equation of motion where q is the velocity potential, κ(x) is the acoustic bulk modulus, ρ(x) is the (mass) density, u(x, t) is the displacement vector, and f (x, t) is a scalar source term. Generally, the acoustic medium is fully described by only two parameters, namely, the density ρ(x) and the sound speed

     1 κ q = ∇ • u + f u = 1 ρ ∇q, (2)
V (x) so that κ(x) = ρ(x)V 2 (x)
, and the density is kept constant during the numerical experiments below. And the boundary conditions are free normal stress all around the domain, and the initial velocity potential q and its first time derivative q are zero.

Removing the CFL limit of the explicit time-domain SEM

Like in a standard finite-element method, for generating spatial meshing discretization, the original domain is discretized into nonoverlapping quadrilateral elements. Then the ordinary differential equation in SEM could be built by assembling the element mass and stiffness matrices, which is based on the numerical integration of the weak form of the equation 2 [START_REF] Komatitsch | The spectral element method : an effective tool to silmulate the seismic response of 2d and 3d geological structures[END_REF]. The resultant matrix expression of the acoustic wave equation 2 by the SEM can be written as:

M ∂ 2 Q ∂t 2 + KQ = F, (3)
where M is the diagonal global mass matrix, K is the sparse global stiffness matrix, F is the source vector, Q and ∂ 2 Q ∂t 2 are the unknown velocity potential and its acceleration vector.

And then, for the temporal discretization, we only consider the second-order accurate temporal finite-difference operator, which is the most commonly used time-marching scheme in the SEM [START_REF] Komatitsch | The spectral element method : an effective tool to silmulate the seismic response of 2d and 3d geological structures[END_REF][START_REF] Chaljub | Spectral element analysis in seismology[END_REF]. The form of equation 3 is changed as follows:

Q n+1 -2Q n + Q n-1 = -∆t 2 M -1 K A Q n + ∆t 2 M -1 F n , ( 4
)
where n is related to the time moment n∆t. Taking a structural meshing as an example, the size

of the matrix A = -∆t 2 M -1 K is (N x × N z ) 2 , N
x and N z are the number of grid points of SEM along the x and z directions for 2D cases. The left side of equation 4 in the frequency domain can be expressed as

(e iω∆t -2 + e -iω∆t ) Q = [2 cos(ω∆t) -2] Q = -4 sin 2 (ω∆t/2) Q, ( 5
)
1 2 3 4 5 6 7 8 9
This paper presented here as accepted for publication in Geophysics prior to copyediting and composition. © 20 Society of Exploration Geophysicists 21

Q = F + [Q(t)],
and F + represents the forward Fourier transform. The coefficient range of Q should be in [-4, 0]. The matrix A is a semi-negative definite matrix, like the updated matrix in FD, only has nonpositive real eigenvalues [START_REF] Gaffar | An explicit and unconditionally stable fdtd method for electromagnetic analysis[END_REF][START_REF] Li | Overcoming the fdtd stability limit via model order reduction and eigenvalue perturbation[END_REF][START_REF] Gao | Removing the stability limit of the explicit finitedifference scheme with eigenvalue perturbation[END_REF].

It has been demonstrated that the over-sampling nature of the spatial grids in the time-domain FD method, which is directly related to the construction of the matrix A in the FD framework, provides the theoretical justification for extending the CFL limit [START_REF] Li | Model order reduction and stability enforcement of finite-difference time-domain The CFL's removal of SEM equations beyond the cfl limit[END_REF][START_REF] Gao | Removing the stability limit of the explicit finitedifference scheme with eigenvalue perturbation[END_REF]. And a time step ∆t > ∆t CFL makes some eigenvalues of the matrix A outside of [-4, 0], which leads to the unstable phenomena. Fortunately, [START_REF] Li | Model order reduction and stability enforcement of finite-difference time-domain The CFL's removal of SEM equations beyond the cfl limit[END_REF] proposes that an enforcement process may be performed by perturbing the unstable eigenvalues to ensure the stability, by retaining the stable eigenvalues without significantly affecting the accuracy. And this method is called the eigenvalue perturbation [START_REF] Li | Model order reduction and stability enforcement of finite-difference time-domain The CFL's removal of SEM equations beyond the cfl limit[END_REF]. All of these features have only direct relationships with the iterative matrix A no matter it comes from the FD, the SEM or the other numerical wave propagation methods, so the iterative matrix A constructed by the SEM should also be able to be perturbed. As a preprocessing step, the eigenvalue perturbation of the iterative matrix A consists of the following three steps:

(1) compute the eigenvalue decomposition of the iterative matrix A = WΛW -1 , where the matrix W includes all the eigenvectors,

Λ = diag(λ 1 , λ 2 , • • • , λ m ) is a diagonal matrix con-
taining all the eigenvalues of the matrix A.

(2) perturb the unstable eigenvalues λ i < -4 in Λ to λi = -4 to construct a new diagonal matrix Λ in order that the stability of equation 4 is remained.

(3) reconstruct a new updated iterative matrix  by the new perturbed eigenvalues.  = W ΛW -1 .

The first prepossessing step can be time-consuming especially for large models (Gao et al., 2019).

Nevertheless, some local eigenvalues approximation method can be further used to find the relevant eigenvalues to the poor CFL conditions [START_REF] Yan | Fast explicit and unconditionally stable fdtd method for electromagnetic analysis[END_REF][START_REF] Cottereau | Stability of an explicit high-order spectral element method for acoustics in heterogeneous media based on local element stability criteria[END_REF] for a reduced numerical cost. Moreover, this eigenvalue decomposition only needs to be performed only once throughout the whole workflow for a given model and a given spatial grid. After the eigenvalue perturbation, equation 4 can be rewritten as where  is the new perturbed updated matrix, and the size of the full  is the same as the sparse A matrix. The second-order time derivative in the wave equation is typically computed by using the FD approximation. [START_REF] Koene | Eliminating time dispersion from seismic wave modeling[END_REF] mentions that the FTDT and ITDT methods could also be used in SEM to remove the time-dispersion error for the second-order accurate Newmark scheme [START_REF] Newmark | A method of computation for structural dynamics[END_REF]. Here the second-order accurate Newmark scheme is used to update the Q n , Qn and Qn to the Q n+1 , Qn+1 and Qn+1 as follows:

Q n+1 -2Q n + Q n-1 = ÂQ n + ∆t 2 M -1 F n , (6)
           Q n+1 = Q n + ∆t Qn + 0.5∆t 2 Qn Qn+1 = Qn + 0.5∆t Qn + Qn+1 Qn+1 = 1 ∆t 2 ÂQ n+1 + M -1 F n+1 . (7)
Note that the temporal central difference format and the above Newmark format are equivalent.

From the Newmark scheme, it is known that

Q n = Q n-1 + ∆t Qn-1 + 0.5∆t 2 Qn-1 . (8)
Then we substract the first equation of equations 7 and equation 8, after replacing ( Qn -Qn-1 ) by the second equation of equations 7, their equivalence could be obtained as follows:

           Q n+1 -2Q n + Q n-1 = ∆t( Qn -Qn-1 ) + 0.5∆t 2 ( Qn -Qn-1 ) = 0.5∆t 2 ( Qn + Qn-1 ) + 0.5∆t 2 ( Qn -Qn-1 ) = 0.5∆t 2 Qn (9)
Thus, the corresponding FTDT and ITDT methods could be used in the above Newmark scheme to remove the time-dispersion error, which are derived from the temporal centre difference format.

It's also worth noting that, in the computation of the classical SEM, the global stiffness matrix K is not usually assembled, and the product vector KQ, as a whole part, is calculated element by element based on the tensor product [START_REF] Komatitsch | The spectral element method : an effective tool to silmulate the seismic response of 2d and 3d geological structures[END_REF]. However, here we need assemble the global stiffness matrix K to construct the updated matrix A, and then obtain the new updated matrix  through eigenvalue perturbation. Eliminating the time-dispersion error from a time step beyond the CFL stability limit

As the theoretical analyses illustrated above, now we could adopt a time step ∆t > ∆t CFL to solve the equation 7 after applying the eigenvalue perturbation. However, it will result in an extremely large time-dispersion error due to the O(∆t 2 ) error accumulation for the second-order Newmark time scheme. Fortunately, we could eliminate the time-dispersion error by the combination of the FTDT and ITDT methods [START_REF] Wang | Finite-difference time dispersion transforms for wave propagation[END_REF][START_REF] Koene | Eliminating time dispersion from seismic wave modeling[END_REF][START_REF] Gao | Removing the stability limit of the explicit finitedifference scheme with eigenvalue perturbation[END_REF]. A simplified introduction of the FTDT and ITDT methods is given below.

After the second-order time discretization at the left side of the equation 4, the relation between the original ω 0 and time-dispersed ω frequencies is 2nd-order :

     -ω 2 0 q ≈ (e iω∆t -2+e -iω∆t ) ∆t 2 q ω 0 ≈ 2 ∆t sin( ω∆t 2 ), (10)
where q = F + [q(t)]. We thus could map the correct and dispersed phase rotations of a recording by adding dispersion in a forward Fourier transform to a trace q(t):

q(t) = F -1 ∞ -∞ q(t) e -i2 sin( ω∆t 2 ) t ∆t dt . (11)
We refer to this operation as the FTDT method. On the contrary, we could map a temporally dispersed signal back to the correct phase in a similar way as the FTDT method:

q(t) = F -1 ∞ -∞ q(t) e -i2 sin -1 ( ω 0 ∆t 2 ) t ∆t dt . (12)
And we refer to this operation as the ITDT method. For the sake of completeness, here we summarize the workflow about how to eliminate the time-dispersion error from a time step ∆t > ∆t CFL by the combination of the FTDT method, the eigenvalue perturbation, and the ITDT method in the SEM framework as follows:

(1) Apply the FTDT method to the source time function. This preprocessing of the source ensures that we impose a wavelet with an identical amount of time-dispersion error as that accumulating in the wavefields [START_REF] Koene | Eliminating time dispersion from seismic wave modeling[END_REF].

1 2 3 4 5 6 7 8 9
This paper presented here as accepted for publication in Geophysics prior to copyediting and composition. © 20 Society of Exploration Geophysicists 21 (2) Perform the acoustic wave numerical simulation with the new updated matrix  and a time step ∆t > ∆t CFL by the SEM, and record the desired traces.

(3) Apply the ITDT method to the recorded traces. This step can recover the original waveforms without the time-dispersion error.

NUMERICAL EXPERIMENTS

To further numerically validate the combined method, numerical experiments in the 2D homogeneous in the 2D homogeneous, moderate and strong heterogeneous models are followed. [START_REF] Seriani | Spectral element method for acoustic wave simulation in heterogeneous media[END_REF][START_REF] Capdeville | 2-d non-periodic homogenization to upscale elastic media for p-sv waves[END_REF]. The corresponding maximum time step available for the second-order temporal discretization is ∆t CFL ≈ 0.076 ms (Courant constant ≈ 0.6).

2D homogeneous model case

(
The spatial interval ∆x min is too small to be adapted in real numerical simulations. After using the combined method, we verify five larger time steps: ∆t = 1, 3, 5, 7 and 9 ms, which are all tens of times larger than ∆t CFL . A dense enough meshing containing 40 × 40 elements with 11 GLL points (G ≈ 10.7) and a small enough time step (∆t = 0.1 ms) are used to obtain the reference waveforms and snapshots without spatial and time-dispersion errors. And they will be used to benchmark the numerical accuracy of the waveforms and snapshots calculated by the five larger time steps. error, which is in agreement with the FD case in [START_REF] Gao | Removing the stability limit of the explicit finitedifference scheme with eigenvalue perturbation[END_REF]. Even though the time-dispersion error using ∆t = 7 ms is not very obvious after ITDT method, some oscillations suddenly arise near the source location in Figure 1b4 and in Figure2b for the ∆t = 7 ms case at about 2.95 s. And these phenomena become uncontrollable and unusable when the ∆t is chosen towards or above the Nyquist sampling time step ∆t Nyq = 1/(2f max ). For the Ricker wavelet we used, the maximum frequency is about 70 Hz, and the time step should be smaller than 7.1 ms depending on the Nyquist limit ∆t < ∆t Nyq . The waveform from ∆t = 9 ms in Figure 2b is not accurate any more especially for the high frequency components, due to the insufficient sampling of source frequency spectrum.

Thus, the time step cannot be infinitely large and still constrained by the Nyquist sampling theorem.

2D heterogenous cosine model case

Furthermore, we validate the combined method in a 2D smooth cosine heterogeneous model, which is a typical representation of the homogenized model. Its oscillating property with space needs to be accurately integrated, and a larger spatial sampling G than that in homogeneous models [START_REF] Capdeville | 2-d non-periodic homogenization to upscale elastic media for p-sv waves[END_REF] needs to be adapted. Our cosine model has the size of 2000 m×2000 m, and its velocity range V is from 3600 to 4400 m/s, which is a moderate variation of the velocity. The source time function is the same as the homogeneous case. The only receiver is located at (0 m, 2000 m). The associated spatial meshing includes 8 × 8 elements with 15 GLL points as shown in Figure 3a, which leads to almost no spatial-dispersion error owing to the good enough G ≈ 3.2 [START_REF] Seriani | Spectral element method for acoustic wave simulation in heterogeneous media[END_REF][START_REF] Capdeville | 2-d non-periodic homogenization to upscale elastic media for p-sv waves[END_REF]. The associated maximum time step for the second-order After the whole procedure of the combined method, the residuals between the numerical and reference waveforms in Figure 3b are unapparent even for ∆t = 7 ms (only with a little waveform oscillation), which implies that the obtained numerical waveforms have almost no time-dispersion error. The waveform from ∆t = 9 ms in Figure 3b is unusable due to the insufficient sampling of source frequency spectrum. This numerical experiment validates that through the combined method, in the SEM framework, we can adopt a time step, far beyond the CFL stability limit and below the Nyquist sampling theorem, to obtain the resultant waveforms without time-dispersion error in the general homogenized models with moderate velocity variation as well.

2D modified Marmousi model case

In the last numerical experiment, we validate the combined method in a 2D modified Marmousi model with a very long time simulation (about 1000 periods). Its spatially varying property also needs to be accurately integrated, and a larger spatial sampling G than that in homogeneous models needs to be adapted. Our model has the size of 2000 m × 2000 m, and its velocity range V is from 1500 to 5443 m/s, which is a strong variation of the velocity. The source, located at the center of the model, is a 15 Hz Ricker wavelet. The only receiver is located at (0 m, 0 m). The associated spatial meshing includes 10 × 5 elements with 21 GLL points as shown in Figure 4a, which leads to almost no spatial-dispersion error owing to the good precise G ≈ 2.9 for the very high order SEM [START_REF] Seriani | Spectral element method for acoustic wave simulation in heterogeneous media[END_REF][START_REF] Capdeville | 2-d non-periodic homogenization to upscale elastic media for p-sv waves[END_REF]. The associated maximum time step for the second-order Newmark temporal discretization is ∆t CFL ≈ 0.195 ms (Courant constant ≈ 0.61).

Similarly, we verify five larger time steps: ∆t = 3, 5, 7, 9 and 11 ms. They are slightly larger than the above two cases, because the used maximum frequency 52.5 Hz is smaller than the 70 Hz in the above cases, we could have a larger Nyquist sampling time step about 9.5 ms. And the reference waveform is calculated by the same meshing and a small enough time step ∆t = 0.05 ms almostly without time-dispersion error. After the calculation of using the combined method, the residuals between the numerical and reference waveforms in Figure 4b are unapparent even for ∆t = 9 ms (only with a little waveform oscillation), which implies that the obtained numerical waveforms have

1 2 3 4 5 6 7 8 9
This paper presented here as accepted for publication in Geophysics prior to copyediting and composition. © 20 Society of Exploration Geophysicists 21 In summary, for the homogeneous, moderate and strong heterogeneous models, we use time steps beyond the CFL limit even towards the Nyquist limit, and the resultant waveforms almostly have no time-dispersion error by the combined methods. For the time step above the Nyquist limit, the waveforms from ∆t = 9 ms in Figures 2b and3b and the waveform from ∆t = 11 ms in Figure 4b are not accurate any more especially for the high-frequency components, due to the insufficient sampling of source frequency spectrum.

DISCUSSION

As one of the methods to eliminate or extend the CFL stability limit of explicit scheme at present, the spatial filtering method is a more basic approach to understand the exact meaning of the CFL stability limit and how to extend this limit. However, the eigenvalue perturbation method performs better than the spatial filtering method for its feasibility for strongly inhomogeneous media.

It is well known that the GLL points cluster towards the edge of the element for very high degree SEM, which seriously restricts its application especially for the model with fine structures.

And the small available time steps of using the very high degree SEM for the homogenized model still constraint its real application. Based on the combined method, the time step is not limited by the minimum spatial GLL interval ∆x min anymore, which is a huge advantage because the very high-degree SEM with a low G can be used to perform simulation without the spatial-and timedispersion errors. As the homogenous numerical experiment shows, the very high-degree SEM can adopt tens of times larger than ∆t CFL to perform simulation owing to its extremely small space.

However, for FD method due to its general uniform spatial grid, the available time steps are only several times of the ∆t CFL . It seems that the resultant waveform calculated by the time steps near the Nyquist sampling in SEM is much stabler than the corresponding waveform in the FD, which can be observed between the Figure 3b in Gao et al. (2019) and the Figure 2b in this paper.

It's also worth noting that the eigenvalue perturbation of the iterative matrix A is a pre-processing procedure, and only needs to be calculated once throughout the whole workflow for a given model and a given spatial grid. This will greatly reduce the calculation time for some technical methods, such as the full waveform inversion and migration who require numerical simulations of multiple sources for the same model, which shows a potential application prospect. Although the combined method with the SEM is very promising, until now it is still very difficult to be applied in large models. The new updated iterative matrix  is full, so a lot of memory and large computational sources are needed. For a 2D model meshed by 20 × 20 spectral elements with 11 GLL points in each orthogonal direction, the memory cost of matrix  is about 6.0 GB. And the calculation of eigenvalue of matrix A increases as exponential as the number of points [START_REF] Gao | Removing the stability limit of the explicit finitedifference scheme with eigenvalue perturbation[END_REF], which greatly limits the application of this method in practice. And in order to build the updated iterative matrix Â, we need the global stiffness matrix K. However, for the classic SEM, the global stiffness matrix K is not usually built, and the product vector KQ is calculated element by element based on the tensor product as a whole part. For example, in the calculation of our Marmousi model case, the size of the matrix A is 20301 × 20301, which takes 4619 s for the eigenvalue related calculation.

While in the following SEM long time simulation (about 1000 periods) with the full stiffness matrix Â, it takes 3058 s. And the calculation of the reference waveform with fine time step ∆t = 0.05 s takes 1436 s, which is based on the traditional spectral-element method in the element by element way also shown in Table 1.

Next, we compare the computational complexity analysis of the acoustic FWI by SEM numerical simulation and eigenvalue calculation. For the acoustic FWI in the d dimensions with ns sources, nr receivers, and nt time steps, regarding the number of elements ne d with the element degree N , the corresponding computational complexity is nt × (ns + nr) × ne d × N d+1 . Note that N d+1 is usually where the main optimization effort is put into, following matrix-matrix multiplication strategy [START_REF] Deville | High-order methods for incompressible fluid flow[END_REF]. The conservative computational complexity of the eigenvalue calcualation is about O(n 3 ), where n is the number of degrees of freedom. In previous same model, the associated computational complexity is O(((ne

× N + 1) d ) 3 ) ≈ O((ne × N ) 3d
), for example, by the QR algorithm [START_REF] Parlett | The qr algorithm[END_REF]. And the ratio between SEM and eigenvalue calculation is about O( nt×(ns+nr) ne 2d ×N 2d-1 ). They are both in the usual numerical linear algebra statistics where the basic operations (+, -, *, /) are performed approximately in current machine arithmetic and cost O(1) each. Unless the computational complexity of solving eigenvalues increases significantly, obtaining all the eigenvalues and applying the eigenvalue perturbation for 3D large-scale models will remain Page 14 of 27 GEOPHYSICS 1 2 3 4 5 6 7 8 9

This paper presented here as accepted for publication in Geophysics prior to copyediting and composition.

© 20 Society of Exploration Geophysicists 21

Downloaded 07/29/21 to 193.52.109.12. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/page/policies/terms DOI:10.1190/geo2020-0623.1 The CFL's removal of SEM impossible. So to avoid globally calculating the stiffness matrix K and the eigenvalues of large matrices for real applications, some local eigenvalues approximation method can be further used to eliminate the relevant eigenvalues that destroy the CFL conditions to make use of the larger time step ∆t > ∆t CFL [START_REF] Yan | Fast explicit and unconditionally stable fdtd method for electromagnetic analysis[END_REF][START_REF] Cottereau | Stability of an explicit high-order spectral element method for acoustics in heterogeneous media based on local element stability criteria[END_REF]. In the global seismology, the eigenvalue calculation can be parallelized by the angular orders (wave number) for the spherically symmetric earth models as used in [START_REF] Capdeville | Coupling the spectral element method with a modal solution for elastic wave propagation in global earth models[END_REF] to calculate the limited normal modes [START_REF] Takeuchi | [END_REF]. For the more general cases, the MPI and OpenMP could be used for the parallelization computation of finding some eigenvalues, such as the FEAST package [START_REF] Polizzi | Density-matrix-based algorithm for solving eigenvalue problems[END_REF]: http://www.ecs.umass.edu/~polizzi/feast/. After obtaining the new iterative matrix by eigenvalue perturbation, one could store it for all subsequent forward simulations in the model. And another very important content about how to further implement the combined method in the numerical modeling with absorbing boundary conditions is also an issue to be addressed in the future.

CONCLUSION

In the explicit time-domain SEM framework for solving the 2D acoustic wave equation, based on the verifications of theoretical analyses and numerical experiments, we succeed in combining the FTDT method, the eigenvalue perturbation, and ITDT method to use time steps beyond the CFL limit even towards the Nyquist limit, and the resultant waveforms exhibit no time-dispersion error.

It lays an initial foundation that we probably could draw lessons from all the efficient eigenvalue calculation methods to perform the numerical wave simulation by the SEM with a very large time step ∆t > ∆t CFL in the complex models. The red star presents the location of the source and the blue triangle presents the receiver. The velocity V (x) = V 0 f (x), where V 0 is the constant velocity used in the homogeneous case.

And the heterogeneous mechanical properties have been chosen with the form f (

x) = 1 + a cos( 2π T h (k a •x)+cos( 2π λ h (k b •x)) , where k a = cos(t 1 ), sin(t 1 ) , k b = cos(t 2 ), sin(t 2 ) .
The amplitude a is adapted 0.5 and the two crossed directions are t This paper presented here as accepted for publication in Geophysics prior to copyediting and composition. © 20 Society of Exploration Geophysicists 21 This paper presented here as accepted for publication in Geophysics prior to copyediting and composition. © 20 Society of Exploration Geophysicists 21 This paper presented here as accepted for publication in Geophysics prior to copyediting and composition. © 20 Society of Exploration Geophysicists 21 This paper presented here as accepted for publication in Geophysics prior to copyediting and composition.

© 20 Society of Exploration Geophysicists 21

Firstly

  , we carry out numerical experiments in a 2D homogeneous square model with the size 2000 m × 2000 m by the SEM. The wave velocity is constant with V = 4000 m/s. The only receiver is chosen at (0 m, 2000 m). And the source, located at the center of the model, is a 20 Hz Ricker wavelet. In terms of the memory of the updated matrix A, we adopt only one element with 121 GLL points in x and z orthogonal directions for the spatial meshing, which leads to almost no spatial-dispersion error based on the empirical value of points per minimum wavelength G ≈ 3.2

Figure 1

 1 Figure 1 displays four snapshots at 1470 ms simulated by SEM with four different time steps, and Figure 2 shows the numerical waveforms with five time steps and the associated reference waveform. Figures 1a and 2a show that, as the time step ∆t increases the time-dispersion error becomes more and more obvious and serious in the snapshots and waveforms. After applying the

Figure Captions ••

 Captions Figure Captions

  1 = π 4 , t 2 = 3π 4 . The model spatial wavelength λ h = 50m. (b) Waveforms calculated by five different time steps ∆t, from 1 ms to 9 ms. The corresponding black waveforms shown with black solid lines are calculated by applying the combination of the FTDT method, the eigenvalue perturbation and the ITDT method, and the red dashed line is the reference waveform calculated by a fine time step. • Figure 4: Configuration of the 2D smooth modified Marmousi model and the waveforms calculated by five different time steps with a relatively long time simulation. (a) The modified Marmousi model. The gray solid lines display the spatial meshing by the 10 × 5 spectral elements. The red star presents the location of the source and the blue triangle presents the receiver. The velocity range is from 1500 to 5500 m/s. (b) Waveforms calculated by five

Figure1:

  Figure1: Snapshots obtained by four different time steps at 1470 ms. (a) (a1-a4) are obtained by time steps 1, 3, 5, and 7 ms with the FTDT method and eigenvalue perturbation respectively. (b) (b1-b4) are the corresponding snapshots in (a) after applying the ITDT method.249x132mm (300 x 300 DPI)

  Downloaded 07/29/21 to 193.52.109.12. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/page/policies/terms DOI:10.1190/geo2020-0623.1

Figure2:

  Figure2: Waveforms calculated by five different time steps. (a) are the waveforms obtained by time steps from 1 ms to 9 ms after using the FTDT method and eigenvalue perturbation. The five waveforms in black solid lines from bottom to top are calculated by time steps from 1ms to 9ms respectively, and each red dashed line displays the reference waveform calculated by a very small Δt. (b) are the corresponding black waveforms in (a) after applying the ITDT method, shown with black solid lines, and the red dashed line is the reference waveform. 203x127mm (300 x 300 DPI)

  Figure3: Configuration of the smooth 2D cosine heterogeneous model and the waveforms calculated by five different time steps with a relatively long time simulation. (a) The 2D cosine model. The gray solid lines display the spatial meshing by the 8 × 8 spectral elements. The red star presents the location of the source and the blue triangle presents the receiver. The velocity V(x) = V 0 f(x), where V 0 is the constant velocity used in the homogeneous case. And the heterogeneous mechanical properties have been chosen with the form f(x) = 1 + a(cos( 2π/T h (ka • x) +cos( 2π/λ h (kb • x))) , where ka = (cos(t 1 ), sin(t 1 )) , kb = (cos(t 2 ), sin(t 2 )) . The amplitude a is adapted 0.5 and the two crossed directions are t1 = π/4, t2 = 3/4π. The model spatial wavelength λ h = 50m. (b) Waveforms calculated by five different time steps ∆t, from 1 ms to 9 ms. The corresponding black waveforms shown with black solid lines are calculated by applying the combination of the FTDT method, the eigenvalue perturbation and the ITDT method, and the red dashed line is the reference waveform calculated by a fine time step. 317x159mm (300 x 300 DPI)

Figure4:

  Figure4: Configuration of the 2D smooth modified Marmousi model and the waveforms calculated by five different time steps with a relatively long time simulation. (a) The modified Marmousi model. The gray solid lines display the spatial meshing by the 10 × 5 spectral elements. The red star presents the location of the source and the blue triangle presents the receiver. The velocity range is from 1500 to 5500 m/s. (b) Waveforms calculated by five different time steps Δt, from 3 ms to 11 ms. The corresponding black waveforms shown with black solid lines are calculated by applying the combination of the FTDT method, the eigenvalue perturbation and the ITDT method, and the red dashed line is the reference waveform computed by a fine time step. 264x159mm (300 x 300 DPI)

  Downloaded 07/29/21 to 193.52.109.12. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/page/policies/terms

DOI:10.1190/geo2020-0623.1 The CFL's removal of SEM

  Downloaded 07/29/21 to 193.52.109.12. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/page/policies/terms DOI:10.1190/geo2020-0623.1 The CFL's removal of SEM waveform is calculated by the same meshing and a fine time step ∆t = 0.05 ms almost without the time-dispersion error.

	1
	2
	3
	4
	5
	6
	7
	8
	9

Newmark temporal discretization is ∆t CFL ≈ 0.6 ms (courant constant ≈ 0.608). Similar to the homogeneous case, we verify five larger time steps: ∆t = 1, 3, 5, 7 and 9 ms. And the reference This paper presented here as accepted for publication in Geophysics prior to copyediting and composition.

© 20 Society of Exploration Geophysicists 21

  Downloaded 07/29/21 to 193.52.109.12. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/page/policies/terms

DOI:10.1190/geo2020-0623.1 The CFL's removal of SEM almostly no time-dispersion error.

  Downloaded 07/29/21 to 193.52.109.12. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/page/policies/terms

	Page 13 of 27	GEOPHYSICS
	1
	2
	3
	4
	5
	6
	7
	8
	9

This paper presented here as accepted for publication in Geophysics prior to copyediting and composition.

© 20 Society of Exploration Geophysicists 21 DOI:10.1190/geo2020-0623.1 The CFL's removal of SEM

  This paper presented here as accepted for publication in Geophysics prior to copyediting and composition.© 20 Society of Exploration Geophysicists 21 Downloaded 07/29/21 to 193.52.109.12. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/page/policies/terms This paper presented here as accepted for publication in Geophysics prior to copyediting and composition. © 20 Society of Exploration Geophysicists 21 Downloaded 07/29/21 to 193.52.109.12. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/page/policies/terms This paper presented here as accepted for publication in Geophysics prior to copyediting and composition. © 20 Society of Exploration Geophysicists 21 Downloaded 07/29/21 to 193.52.109.12. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/page/policies/terms This paper presented here as accepted for publication in Geophysics prior to copyediting and composition. © 20 Society of Exploration Geophysicists 21 Downloaded 07/29/21 to 193.52.109.12. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/page/policies/terms This paper presented here as accepted for publication in Geophysics prior to copyediting and composition. © 20 Society of Exploration Geophysicists 21 Downloaded 07/29/21 to 193.52.109.12. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/page/policies/terms

	1 1 1 1 1	Lyu & al	20	The CFL's removal of SEM
	2 2 2 2 2
	3 3 3 3 3
	4 4 4 4 4
	5 5 5 5 5
	6 6 6 6 6
	7 7 7 7 7
	8 8 8 8 8
	9 9 9 9 9
	DOI:10.1190/geo2020-0623.1 DOI:10.1190/geo2020-0623.1 DOI:10.1190/geo2020-0623.1 DOI:10.1190/geo2020-0623.1 DOI:10.1190/geo2020-0623.1

Table Caption Table 1 :

 Caption1 Comparison of calculated time for the Marmousi model case. The computations are finished in the 2020 MacBook Pro16 with 2.4 GHz core and 64G memory in Matlab 2020b. This paper presented here as accepted for publication in Geophysics prior to copyediting and composition. © 20 Society of Exploration Geophysicists 21 Downloaded 07/29/21 to 193.52.109.12. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/page/policies/terms DOI:10.1190/geo2020-0623.1 The CFL's removal of SEM

	1
	2
	3
	4
	5
	6
	7
	8
	9

Table 1 :

 1 Downloaded 07/29/21 to 193.52.109.12. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/page/policies/terms DOI:10.1190/geo2020-0623.1 The CFL's removal of SEM different time steps ∆t, from 3 ms to 11 ms. The corresponding black waveforms shown with black solid lines are calculated by applying the combination of the FTDT method, the eigenvalue perturbation and the ITDT method, and the red dashed line is the reference waveform computed by a fine time step. Downloaded 07/29/21 to 193.52.109.12. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/page/policies/terms Comparison of calculated time for the Marmousi model case. The computations are finished in the 2020 MacBook Pro16 with 2.4 GHz core and 64G memory in Matlab 2020b. Downloaded 07/29/21 to 193.52.109.12. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/page/policies/terms

	1
	2
	3
	4
	5
	6
	7
	8
	9

This paper presented here as accepted for publication in Geophysics prior to copyediting and composition.

© 20 Society of Exploration Geophysicists 21

This paper presented here as accepted for publication in Geophysics prior to copyediting and composition. © 20 Society of Exploration Geophysicists 21 DOI:10.1190/geo2020-0623.1

ACKNOWLEDGEMENTS

This study is supported by the National Key R&D Program of China (Grant No. 2017YFC0601206), NSFC (grant no. 41888101, 42004045, 41625016, and 41774112), the China Scholarship Council (File No. 201804910289) and the French National Research Agency (ANR-16-CE31-0022-01, [2016][2017][2018][2019][2020][2021]. We thank the editor-in-chief Jeffrey Shragge, the assistant editors Joakim Blanch and

The CFL's removal of SEM Stig Hestholm, and five anonymous reviewers, for their comments which greatly helped to improve the manuscript. Chao Lyu also sincerely thank Yingjie Gao and Yan Jin for plentiful discussions on the CFL's removal.


Benoît Bonnet
email: benoit.bonnet@imj-prg.fr

Hélène Frankowska
email: helene.frankowska@imj-prg.fr

On the Properties of the Value Function Associated to a Mean-Field Optimal Control Problem of Bolza Type *

In this paper, we obtain several structural results for the value function associated to a mean-field optimal control problem of Bolza type in the space of measures. After establishing the sensitivity relations bridging between the costates of the maximum principle and metric superdifferentials of the value function, we investigate semiconcavity properties of this latter with respect to both variables. We then characterise optimal trajectories using set-valued feedback mappings defined in terms of suitable directional derivatives of the value function.

Introduction

During the past decade, variational problems formulated on mean-field approximations of collective dynamics have gained a tremendous amount of steam. Indeed, the mathematical field of multi-agent systems analysis has been for several years the stage on which modern theories such as mean-field control [START_REF] Bensoussan | Mean Field Games and Mean Field Type Control Theory[END_REF][START_REF] Carmona | Probabilistic Theory of Mean-Field Games with Applications[END_REF] and mean-field games [START_REF] Huang | Large Population Stochastic Dynamic Games : Closed-Loop McKean-Vlasov Systems and the Nash Certainty Equivalence Principle[END_REF][START_REF] Lasry | Mean Field Games[END_REF] are being developed.

The main rationale behind studying the former of these two fields of research can be heuristically summarised as follows. From a practical standpoint, it is usually easier to model collective dynamics by means of large families of microscopic entities -called the agents -, which evolution is prescribed e.g. by ODEs [START_REF] Choi | The Derivation of Swarming Models: Mean-Field Limit and Wasserstein Distances[END_REF] or time-dependent graphs [START_REF] Mesbahi | Graph Theoretic Multi-Agent Systems[END_REF]. However, such formulations often lead to numerically intractable problems, and rarely allow to capture global properties of the system in a meaningful way. For these reasons, an ever-expanding literature has been focusing on the development of mathematical tools for the investigation of adequate macroscopic approximations of such systems, see e.g. [START_REF] Bonnet | Intrinsic Lipschitz Regularity of Mean-Field Optimal Controls[END_REF][START_REF] Cardaliaguet | Notes on Mean-Field Games[END_REF][START_REF] Cardaliaguet | The Master Equation and the Convergence Problem in Mean Field Games[END_REF][START_REF] Piccoli | Control to Flocking of the Kinetic Cucker-Smale model[END_REF] and references therein. In this context, an intensive research effort has been directed towards the establishment of a theory of mean-field optimal control [START_REF] Bonnet | A Pontryagin Maximum Principle in Wasserstein Spaces for Constrained Optimal Control Problems[END_REF][START_REF] Bonnet | Necessary Optimality Conditions for Optimal Control Problems in Wasserstein Spaces[END_REF][START_REF] Bonnet | The Pontryagin Maximum Principle in the Wasserstein Space[END_REF][START_REF] Cavagnari | Generalized Dynamic Programming Principle and Sparse Mean-Field Control Problems[END_REF][START_REF] Cavagnari | Lagrangian, Eulerian and Kantorovich Formulations of Multi-Agent Optimal Control Problems: Equivalence and Gamma-Convergence[END_REF][START_REF] Fornasier | Mean-Field Optimal Control as Gamma-Limit of Finite Agent Controls[END_REF]. This denomination usually refers to optimal control problems formulated in the so-called Wasserstein spaces (see Definition 2 below) of optimal transport [START_REF] Ambrosio | Gradient Flows in Metric Spaces and in the Space of Probability Measures[END_REF], which constitute a convenient and natural framework for the generalisation of the classical geometric results of control theory to this setting.

In this paper, we transpose to this class of problems a series of result that were originally derived in [START_REF] Cannarsa | Some Characterizations of Optimal Trajectories in Control Theory[END_REF], dealing with fine properties of the value function associated to a general family of smooth optimal control problems. With this aim, we first investigate the so-called sensitivity relations -originally established in [START_REF] Flemming | Deterministic and Stochastic Optimal Control[END_REF] -which provide a link between the costate variables stemming from the Pontryagin maximum principle [START_REF] Pontryagin | The Mathematical Theory of Optimal Processes[END_REF] ("PMP" for short) and the superdifferentials of the value function. Among their many possible applications, the sensitivity relations are very useful for expressing sufficient conditions to ensure the optimality of Pontryagin's extremals. We proceed by providing a list of conditions under which the value function is semiconcave [START_REF] Cannarsa | Semiconcave Functions, Hamilton-Jacobi Equations and Optimal Control[END_REF] along suitably chosen interpolating curves. This notion of regularity is indeed very rich, and appears in a wide number of problems originating both in control theory and in the calculus of variations, see for instance [START_REF] Dolcetta | Approximate Solutions of the Bellman Equation of Deterministic Control[END_REF][START_REF] Rifford | Semiconcave Control-Lyapunov Functions and Stabilizing Feedbacks[END_REF]. Moreover in the context of Wasserstein spaces, it has been well established in the monograph [START_REF] Ambrosio | Gradient Flows in Metric Spaces and in the Space of Probability Measures[END_REF] that semiconcavity properties along interpolating curves are in fact the natural path leading to meaningful results on metric superdifferentials in the space of measures. By leveraging some of these results, it is then possible to characterise optimal feedbacks for mean-field optimal control problems, and to provide useful regularity properties for the latter.

The structure of the article is the following. In Section 2, we recollect notions of measure theory and optimal transport, and proceed by recalling several results which are more specific to meanfield optimal control problems in Section 3. In Section 4, we study the sensitivity relations and semiconcavity properties of the value function for a general Bolza problem in Wasserstein spaces. We then make use of these results to investigate optimal feedbacks and their regularity in Section 5.

Preliminaries

In this section, we recall preliminary notions of measure theory and optimal transport, for which we refer to [START_REF] Ambrosio | Functions of Bounded Variations and Free Discontinuity Problems[END_REF] and [START_REF] Ambrosio | Gradient Flows in Metric Spaces and in the Space of Probability Measures[END_REF] respectively.

Optimal transport and calculus in Wasserstein spaces

In the sequel, P(R d ) denotes the space of Borel probability measures over R d , endowed with the standard narrow topology (see [START_REF] Ambrosio | Gradient Flows in Metric Spaces and in the Space of Probability Measures[END_REF]Chapter 5]). Given p ∈ [1, +∞), let P p (R d ) be the subset of probability measure which p-momentum M p (µ) := R d |x| p dµ(x)

1/p is finite, and P c (R d ) be that of measures with compact support. For µ ∈ P(R d ), we denote by L p (R d , R d ; µ) the Banach space of maps from R d into itself that are p-summable with respect to µ. When p = 2, the latter is a Hilbert space endowed with the scalar product

ζ, ξ L 2 (µ) := R d ζ(x), ξ(x) dµ(x)
defined here for any ζ, ξ ∈ L 2 (R d , R d ; µ). We will also denote by L 1 the standard Lebesgue measure over R, and consider the associated spaces L p (I, R + ) of p-summable maps from an interval I ⊂ R into R + . Given a Borel map f : R d → R d , we define the pushforward f # µ ∈ P(R d ) of µ through f as the measure f # µ(B) := µ(f -1 (B)), for every Borel set B ⊂ R d .

Definition 1 (Transport plans). The set of transport plans between two measures µ, ν

∈ P(R d ) - denoted by Γ(µ, ν) -is the collection of elements γ ∈ P(R d × R d ) such that π 1 # γ = µ and π 2 # γ = ν, where π 1 , π 2 : R d × R d → R d represent
the projection operator onto the first and second factor respectively.

In the sequel, we will also need the following simplified version of [START_REF] Ambrosio | Gradient Flows in Metric Spaces and in the Space of Probability Measures[END_REF]Theorem 5.3.1].

Theorem 1 (Barycentric projection). Let p ∈ [1, +∞), µ, ν ∈ P p (R d ) and γ ∈ Γ(µ, ν). Then, there exists a unique map γ ∈ L p (R d , R d ; µ) -called the barycentric projection of γ on π 1 # γ = µ -, such that R 2d y, F (x) dγ(x, y) = γ, F L 2 (µ) , for every F ∈ L ∞ (R d , R d ; µ).

Definition 2 (Wasserstein distances). The quantity

W p (µ, ν) := min γ∈Γ(µ,ν) R 2d
|x -y| p dγ(x, y)

1/p
, defines a distance between any two measures µ, ν ∈ P p (R d ), and we denote by Γ o (µ, ν) the set of transport plans for which the minimum is attained.

The Wasserstein space of order p ∈ [1, +∞) is defined as the metric space (P p (R d ), W p ). In the sequel, we will often consider the (non complete) metric space (P c (R d ), W p ) of compactly supported measures equipped with the W p -metric. When p = 2, it has been known since [START_REF] Ambrosio | Gradient Flows in Metric Spaces and in the Space of Probability Measures[END_REF][START_REF] Otto | The Geometry of Dissipative Equations : The Porous Medium Equation[END_REF] that (P 2 (R d ), W 2 ) could be endowed -at least formally -with a differentiable manifold-like structure. We present below a revisited definition of the corresponding calculus, tailored for (P c (R d ), W 2 ), which is borrowed from [START_REF] Bonnet | Necessary Optimality Conditions for Optimal Control Problems in Wasserstein Spaces[END_REF]. Given R > 0, we will use the notation B R (µ) := ∪ x∈supp(µ) B(x, R).

Definition 3 (Local differentials). A functional φ : P 2 (R d ) → R ∪ {±∞} is locally differentiable at µ ∈ P c (R d ) if and only if there exists ∇φ(µ) ∈ L 2 (R d , R d ; µ) such that for all R > 0 and ν ∈ P(B R (µ)), it holds φ(ν) = φ(µ) + R 2d ∇φ(µ), y -x dγ(x, y) + o R (W 2 (µ, ν)), for every γ ∈ Γ o (µ, ν). We also say that φ(•) is locally continuously differentiable if the application (µ, x) ∈ P c (R d ) × R d → ∇φ(µ)(x) ∈ R d is continuous.
In [START_REF] Bonnet | Necessary Optimality Conditions for Optimal Control Problems in Wasserstein Spaces[END_REF], we proved the following chain rule for locally differentiable functionals.

Proposition 1 (Chain rule). Consider elements µ

∈ P c (R d ), F ∈ L ∞ (R d , R d ; µ) and let φ : P 2 (R d ) → R ∪ {±∞} be locally differentiable at µ. Then, it holds φ((Id + ǫF ) # µ) = φ(µ) + ǫ ∇φ(µ), F L 2 (µ) + o(ǫ),
for every ǫ ∈ R.

Non-local continuity equations

Consider the non-local continuity equation in a Wasserstein space

∂ t µ(t) + div x (v(t, µ(t))µ(t)) = 0, µ(0) = µ 0 , (1)
where the non-local velocity-field v :

[0, T ] × P c (R d ) × R d → R d satisfies the following assumptions.

Hypotheses (CE).

Suppose that for any compact set K ⊂ R d , the following holds.

(i) The mapping t → v(t, µ, x) ∈ R d is L 1 -measurable with respect to t ∈ [0, T ] for all (µ, x) ∈ P c (R d ) × R d , and there exists m(•) ∈ L 1 ([0, T ], R + ) such that |v(t, µ, x)| ≤ m(t) 1 + |x| + M 1 (µ) . (ii) There exists a map ℓ K (•) ∈ L 1 ([0, T ], R + ) such that |v(t, ν, y) -v(t, µ, x)| ≤ ℓ K (t) W 1 (µ, ν) + |x -y| for L 1 -almost every t ∈ [0, T ],
any µ, ν ∈ P(K) and all x, y ∈ K.

Under these standard Cauchy-Lipschitz assumptions, we have the following well-posedness result for (1) (see e.g. [START_REF] Bonnet | Differential Inclusions in Wasserstein Spaces: The Cauchy-Lipschitz Framework[END_REF]).

Theorem 2 (Well-posedness of ( 1)). Let r > 0, µ 0 ∈ P(B(0, r)) and suppose that hypotheses (CE) hold. Then, there exists a unique distributional solution µ(•) of (1), which also satisfies

supp(µ(t)) ⊂ B(0, R r ) and W 1 (µ(t), µ(τ )) ≤ t τ m r (s)ds, for all 0 ≤ τ ≤ t ≤ T , with R r > 0 and m r (•) ∈ L 1 ([0, T ], R + )
depending on r > 0 and m(•). Moreover, the curve µ(•) can be represented explicitly as

µ(t) = Φ (τ,t) [µ(τ )](•) # µ(τ ), for every τ, t ∈ [0, T ], where (Φ (τ,t) [µ](•)) τ,t∈[0,T ] denotes the non-local flow starting from µ ∈ P c (R d ) at time τ ∈ [0, T ], defined as the unique solution of    ∂ t Φ (τ,t) [µ](x) = v t, Φ (τ,t) [µ](•) # µ, Φ (τ,t) [µ](x) , Φ (τ,τ ) [µ](x) = x, ( 2
)
for all x ∈ R d .

Mean-field optimal control

In this section, we focus more in details on the following Bolza mean-field optimal control problem (P )

           min u(•)∈U T 0 L(µ(t), u(t))dt + ϕ(µ(T )) , s.t. ∂ t µ(t) + div x (v(t, µ(t), u(t))µ(t)) = 0, µ(0) = µ 0 .
Here

T > 0, v : [0, T ] × P c (R d ) × U × R d → R d is a controlled velocity-field, the set of admissible controls is U := L 2 ([0, T ], U ) where U ⊂ R m is a closed set, and µ 0 ∈ P c (R d ) is a fixed initial datum.
In the present paper, we make the following assumptions on (P ).

Hypotheses (OCP).

Suppose that for any compact set K ⊂ R d , the following holds.

(i) The map u ∈ U → v(t, µ, u, x) ∈ R d is continuous, and there exist constants m, ℓ K > 0 such that    |v(t, µ, u, x)| ≤ m 1 + |x| + M 1 (µ) (1 + |u| 2 ), |v(τ, ν, u, y) -v(t, µ, u, x)| ≤ ℓ K |t -τ | + W 1 (µ, ν) + |x -y| 1 + |u| 2 ),
for all t, τ ∈ [0, T ], any µ, ν ∈ P(K), every u ∈ U and each x, y ∈ K.

(ii) The map u ∈ U → L(µ, u) ∈ R + is lower semicontinuous, and there exist constants

l 1 , C 1 , L 1 K > 0 such that L(µ, u) ≥ C 1 |u| 2 and L(µ, u) ≤ l 1 (1 + M 1 (µ))(1 + |u| 2 ), |L(ν, u) -L(µ, u)| ≤ L 1 K W 1 (µ, ν)(1 + |u| 2 ),
for any µ, ν ∈ P(K) and every u ∈ U .

(iii) There exist C 2 ∈ R and L 2 K > 0 such that

ϕ(µ) ≥ C 2 and |ϕ(ν) -ϕ(µ)| ≤ L 2 K W 1 (µ, ν),
for any µ, ν ∈ P(K).

(iv) The maps (x, µ) → v(t, µ, u, x), L(u, µ), ϕ(µ) are locally continuously differentiable in the sense of Definition 3.

While the assumptions listed above are not minimal for the developments of this section (see e.g. more general conditions in [7, Section 4.2]), they will be necessary to establish several results in Section 4 and Section 5.

Theorem 3 (Existence). Suppose that (OCP) hold and that for any non-empty compact

K ⊂ R d and every (t, µ) ∈ [0, T ] × P c (R d ), the restricted sets (v(t, µ, u) |K , L(µ, u) + r) s.t. u ∈ U, r ≥ 0 ⊂ C 0 (K, R d ) × R,
are closed and convex. Then, there exists an optimal trajectory-control pair (µ * (•), u * (•)) for (P ).

Proof. Under hypothesis (OCP)-(ii), there exists for any minimising sequence (µ

* N (•), u * N (•)) a con- stant C U > 0 such that sup N ≥1 T 0 |u * N (t)| 2 dt ≤ C U .
The result then follows up to minor modifications in the proofs of [6, Theorem 6 and Theorem 7].

In the sequel, we will assume that (P ) admits such an optimal pair (µ * (•), u * (•)). By a simple adaptation of the arguments of [START_REF] Bonnet | Necessary Optimality Conditions for Optimal Control Problems in Wasserstein Spaces[END_REF], it can be shown that the latter satisfies the following Pontryagintype optimality conditions. Theorem 4 (Pontryagin Maximum Principle). Let (µ * (•), u * (•)) be an optimal pair for (P ). Then, there exists a unique state-costate curve ν * : [0, T ] → P c (R 2d ) which satisfies the following.

(i) It solves the Hamiltonian continuity equation

       ∂ t ν * (t) + div x,r ∇ ν H(t, ν * (t), u * (t))ν * (t) = 0, π 1 # ν * (t) = µ * (t)
for all times t ∈ [0, T ],

ν * (T ) = (Id, -∇ϕ(µ * (T ))) # µ * (T ),
where the Hamiltonian H is defined by

H(t, ν, u) := R 2d r, v(t, µ, u, x) dν(x, r) -L(µ, u), ( 3
)
with µ := π 1 # ν and ∇ ν H(t, ν * (t), u * (t)) being its gradient given in the sense of Definition 3.

(ii) The maximisation condition

H(t, ν * (t), u * (t)) = max u∈U H(t, ν * (t), u), ( 4
)
holds for L 1 -almost every t ∈ [0, T ].

Properties of the value function

In this section, we use some of the tools introduced above to study fine properties of the value function

V (τ, µ) :=            min u(•)∈U T τ L(µ(t), u(t))dt + ϕ(µ(T )) , s.t. ∂ t µ(t) + div x (v(t, µ(t), u(t))µ(t)) = 0, µ(τ ) = µ, associated to (P ), defined for (τ, µ) ∈ [0, T ] × P c (R d ).
A first key observation is that under our working assumptions, the mapping V (•, •) is absolutely continuous in its first argument and Lipschitz in the second one.

Proposition 2 (Regularity). Suppose that hypotheses (OCP) hold, and let r > 0. Then, there exist

M r (•) ∈ L 1 ([0, T ], R + ) and L r > 0 such that |V (τ 2 , µ 2 ) -V (τ 1 , µ 1 )| ≤ τ 2 τ 1 M r (t)dt + L r W 1 (µ 1 , µ 2 ),
for every τ 1 , τ 2 ∈ [0, T ] and all µ 1 , µ 2 ∈ P(B(0, r)).

Proof. This result follows from the Lipschitz dependence of solutions of non-local continuity equations on their initial data (see [START_REF] Bonnet | Differential Inclusions in Wasserstein Spaces: The Cauchy-Lipschitz Framework[END_REF]) together with the sublinearity and Lipschitz regularity of the cost functionals.

In what follows, we will first investigate the so-called sensitivity relations involving the costates of the PMP of Theorem 4, and then semiconcavity properties of V (•, •). Given u(•) ∈ U , we will use the notation (Φ u (τ,t) [•](•)) τ,t∈[0,T ] for the non-local flows generated by (t, µ, x) → v(t, µ, u(t), x) ∈ R d defined as in (2).

Sensitivity analysis and Pontryagin costates

In our context, the sensitivity relations link the state-costate curves ν * (•) of the maximum principle and the Hadamard superdifferential of the value function.

Definition 4 (Hadamard superdifferentials). Let

(τ, µ) ∈ [0, T ] × P c (R d ). Then, we say that a pair (δ, ξ) ∈ R × L 2 (R d , R d ; µ) belongs to the Hadamard superdifferential ð + V (τ, µ) of V (•, •) at (τ, µ) provided that for all (h, F ) ∈ R × L ∞ (R d , R d ; µ) lim sup ǫ→0 + τ +ǫh∈[0,T ] V (τ + ǫh, (Id + ǫF ) # µ) -V (τ, µ) ǫ ≤ δh + ξ, F L 2 (µ) .
On vector spaces, Hadamard superdifferentials are usually defined in terms of contingent directional derivatives. Yet, we showed in Proposition 2 above that under (OCP), the map V (t, •) is locally Lipschitz over P c (R d ) with local Lipschitz constants independent of t ∈ [0, T ]. This allows us to use the simpler upper Dini directional derivatives to define Hadamard superdifferentials.

Theorem 5 (Hadamard-type sensitivity). Suppose that hypotheses (OCP) hold. Let (µ * (•), u * (•)) be a solution of (P ) and ν * (•) be the corresponding state-costate curve given by Theorem 4. Then

H(t, ν * (t), u * (t)), -ν * (t) ∈ ð + V (t, µ * (t)), for L 1 -almost every t ∈ [0, T ], where the map ν * (t) ∈ L 2 (R d , R d ; µ * (t)) denotes the barycentric projection of ν * (t) onto its first marginal π 1 # ν * (t) = µ * (t).
Proof. The proof of this result is fairly long and technical, and will thus be published in greater details elsewhere. For the convenience of the reader, we sketch the main steps of the underlying arguments.

Given (τ, ǫ) ∈ [0, T ] × R + and (h, F ) ∈ R × L ∞ (R d , R d ; µ * (τ )), the perturbed curve μǫ (•) defined by μǫ (t) := Φ u * (τ +ǫh,t) [(Id + ǫF ) # µ * (τ )] • (Id + ǫF ) # µ * (τ ),
for all times t ∈ [τ, T ] satisfies

V (τ + ǫh, (Id + ǫF ) # µ * (τ )) ≤ T τ +ǫh L(μ ǫ (t), u * (t))dt + ϕ(μ ǫ (T )),
by definition of V (•, •). Following technical linearisation arguments in the spirit of [START_REF] Bonnet | A Pontryagin Maximum Principle in Wasserstein Spaces for Constrained Optimal Control Problems[END_REF][START_REF] Bonnet | The Pontryagin Maximum Principle in the Wasserstein Space[END_REF], one can show

Φ u * (τ +ǫh,t) [(Id + ǫF ) # µ * (τ )](x + ǫF (x)) = Φ * (τ,t) (x) + ǫw h,F t, Φ * (t,τ ) (x) + o τ,t,x (ǫ),
with sup (τ,t,x) |o τ,t,x (ǫ)| = o(ǫ) and where we denoted Φ * (τ,t) (x) := Φ u * (τ,t) [µ * (τ )](x). Therein, t ∈ [τ, T ] → w h,F (t, y) solves a suitable linearised Cauchy problem.

Combining the chain rule of Proposition 1 with the optimality of (µ * (•), u * (•)), the following inequality

1 ǫ V (τ + ǫh, (Id + ǫF ) # µ * (τ )) -V (τ, µ * (τ )) ≤ T τ ∇ µ L(µ * (t), u * (t)), w h,F (t, Φ * (t,τ ) (•)) L 2 (µ * (τ )) dt + ∇ϕ(µ * (T )), w h,F (T, Φ * (T,τ ) (•)) L 2 (µ * (T )) -hL(µ * (τ ), u * (τ )) + o τ (1),
holds for every ǫ > 0, provided that τ ∈ [0, T ] belongs to a suitable subset of full L 1 -measure. Upon performing lengthy computations involving the dynamics of the state-costate curve ν * (•) introduced in Theorem 4, this inequality can be further reformulated as

1 ǫ V (τ + ǫh, (Id + ǫF ) # µ * (τ )) -V (τ, µ * (τ )) ≤ H(τ, ν * (τ ), u * (τ )) h - R 2d r, F (x) dν * (τ )(x, r) + o τ (1) ≤ H(τ, ν * (τ ), u * (τ )) h -ν * (τ ), F L 2 (µ * (τ )) + o τ (1), (5)
for any ǫ > 0 small enough, and where we used the barycentric projection of ν * (τ ) defined as in Theorem 1. We end the proof by taking the limsup as ǫ → 0 + in ( 5) for L 1 -almost every t ∈ [0, T ].

The main application of the sensitivity relations is to provide sufficient optimality conditions for the PMP.

Theorem 6 (Sufficient optimality conditions). Let (µ * (•), u * (•)) be an admissible pair for (P ) and ν * (•) be an associated state-costate curve satisfying the PMP of Theorem 4 and the sensitivity relation of Theorem 5. Then, the pair (µ * (•), u * (•)) is optimal for (P ).

Proof. Observe that a pair (µ * (•), u * (•)) is optimal for (P ) if and only if

V (τ 2 , µ * (τ 2 )) = V (τ 1 , µ * (τ 1 )) + τ 2 τ 1 L(µ * (t), u * (t))dt, for all 0 ≤ τ 1 ≤ τ 2 ≤ T . By Proposition 2, the map t ∈ [0, T ] → V (t, µ * (t)
) is absolutely continuous, and by definition of the value function it necessary holds

d dt V (t, µ * (t)) ≥ L(µ * (t), u * (t)), (6)
for L 1 -almost every t ∈ [0, T ]. Following [7, Section 3.2], it can be shown that

Φ * (t,t+ǫ) (x) = x + ǫv(t, µ * (t), u * (t), x) + o t,x (ǫ),
for L 1 -almost every t ∈ [0, T ] and for any ǫ ∈ R, with T 0 sup x |o t,x (ǫ)|dt = o(ǫ). This combined with (2) and Proposition 2 then yields

V (t + ǫ, µ * (t + ǫ)) = V (t + ǫ, (Id + ǫv(t, µ * (t), u * (t))) # µ * (t)) + o t (ǫ) (7)
with T 0 |o t (ǫ)|dt = o(ǫ). By the sensitivity relation of Theorem 5, it further holds lim sup

ǫ→0 + 1 ǫ V (t + ǫ, (Id + ǫv(t, µ * (t), u * (t))) # µ * (t)) -V (t, µ * (t)) ≤ -ν * (t), v(t, µ * (t), u * (t)) L 2 (µ * (t)) +H(t, ν * (t), u * (t)) = L(µ * (t), u * (t)), (8)
for L 1 -almost every t ∈ [0, T ], where we used the definition (3) of the Hamiltonian. Upon combining [START_REF] Bonnet | Necessary Optimality Conditions for Optimal Control Problems in Wasserstein Spaces[END_REF] and ( 8), we can thus recover

d dt V (t, µ * (t)) ≤ L(µ * (t), u * (t)),
for L 1 -almost every t ∈ [0, T ], which together with (6) concludes the proof of our claim up to an integration over any time interval of the form [τ 1 , τ 2 ].

Regularity of the value function

Our aim now is to study semiconcavity estimates for V (•, •), taking inspiration from [2, Chapter 9]. We start our developments by recalling in which sense this notion is understood in Wasserstein spaces.

Definition 5 (Interpolation semiconcavity).

A map φ :

P 2 (R d ) → R ∪ {±∞} is locally semiconcave if it is finite on P c (R d )
, and for every r > 0 there exists C r > 0 such that for any λ ∈ [0, 1],

(

1 -λ)φ(µ 1 ) + λφ(µ 2 )-φ µ 1→2 λ ≤ C r λ(1 -λ)W 2 2,µ (µ 1 , µ 2 ),
for any µ 1 , µ 2 ∈ P(B(0, r)) and each µ ∈ Γ(µ, ν), where µ 1→2 λ := (1λ)π 1 + λπ 2 # µ, and

W 2,µ (µ, ν) := R 2d |x -y| 2 dµ(x, y) 1/2
.

Below, we make the following additional assumptions on the data of (P ).

Hypotheses (R). Suppose that the following holds.

(i) For every r > 0, there exists a constant C v r > 0 such that for all λ ∈ [0, 1], it holds

(1 -λ)v(t, µ 1 , u, x) + λv(t, µ 2 , u, y) -v(t, µ 1→2 λ , u, (1 -λ)x + λy)| ≤ C v r λ(1 -λ) W 2 2,µ (µ 1 , µ 2 ) + |x -y| 2 ,
for all (t, u) ∈ [0, T ] × U , any µ 1 , µ 2 ∈ P(B(0, r)) and each x, y ∈ B(0, r).

(ii) The maps µ ∈ P c (R d ) → L(µ, u), ϕ(µ) are locally semiconcave, uniformly with respect to u ∈ U .

Theorem 7 (Semiconcavity). Suppose that hypotheses (OCP) and (R) hold. Then for every r > 0, there exists C r > 0 such that for any λ ∈ [0, 1], it holds

(1 -λ)V (τ 1 , µ 1 ) + λV (τ 2 , µ 2 ) -V (1 -λ)τ 1 + λτ 2 , µ 1→2 λ ≤ C r λ(1 -λ) |τ 1 -τ 2 | 2 + W 2 2,µ (µ 1 , µ 2 ) ,
for every τ 1 , τ 2 ∈ [0, T ], any µ 1 , µ 2 ∈ P(B(0, r)) and each µ ∈ Γ(µ 1 , µ 2 ).

Proof. The proof of this result is again fairly long and technical, and will be exposed in greater details in a subsequent publication. The core idea is to use the refined superposition results from To recover the semiconcavity with respect to both variables, one must perform a suitable rescaling of the time variable in the spirit of [START_REF] Cannarsa | Some Characterizations of Optimal Trajectories in Control Theory[END_REF]Section 5], combined with the semiconcavity estimates already established with respect to the measure variable.

Optimal feedbacks

By leveraging several results of Section 4, it is possible to derive a general characterisation of optimal feedbacks for mean-field optimal control problems. Throughout this section, we assume that hypotheses (OCP) hold. We start by recalling the definition of an upper-semicontinuous set-valued map ("u.s.c." for short). Definition 6 (Regular set-valued maps). A multifunction G : S ⇒ X from a metric space (S , d S ) into a Banach space (X, • X ) is upper-semicontinuous if for any s ∈ S and every neighbourhood Ω of G(s), there exists η > 0 such that G(s ′ ) ⊂ Ω for every s ′ ∈ B(s, η). In particular, if G(•) is u.s.c. and single-valued, it is then continuous in the usual sense.

In order to characterise optimal feedbacks for (P ), we need to introduce two notions of inferior derivatives for functionals defined over the space of measures.

Definition 7 (Two notions of lower-derivatives). Let φ

: P 2 (R d ) → R ∪ {±∞} and µ ∈ P c (R d ), φ(µ) = {±∞}. Given F ∈ C 0 (R d , R d ), we define the lower Dini derivative of φ(•) at µ in the direction F as ∂ - D φ(µ)(F ) := lim inf ǫ→0 + φ((Id + ǫF ) # µ) -φ(µ) ǫ .
Given R > 0, we also define the regularised R-lower derivative of φ(•) at µ in the direction F as

∂ - o,R φ(µ)(F ) := lim inf ν→µ, ǫ→0 + ν∈P(B R (µ)) φ((Id + ǫF ) # ν) -φ(ν) ǫ .
Proposition 3 (The case of semiconcave functions). Suppose that φ : P 2 (R d ) → R ∪ {±∞} is locally semiconcave and let µ ∈ P c (R d ). Then, it holds that

∂ - D φ(µ)(F ) = ∂ - o,R φ(µ)(F ), for every R > 0 and all F ∈ C 0 (R d , R d ).
Proof. Adapt the arguments of [10, Section 3].

In the following theorem, we combine these notions of non-smooth analysis to provide a general characterisation of optimal feedbacks for (P ).

Theorem 8 (Optimal feedbacks

). An admissible pair (µ * (•), u * (•)) is optimal for (P ) if and only if

v(t, µ * (t), u * (t)) ∈ G(t, µ * (t)), (9)
for L 1 -almost every t ∈ [0, T ], where we defined the sets In conclusion, the results of Theorem 8 inform us that every optimal curve µ * (•) to a problem of the form (P ) can be written as the solution of a continuity inclusion defined in the sense of [START_REF] Bonnet | Differential Inclusions in Wasserstein Spaces: The Cauchy-Lipschitz Framework[END_REF].

G(t, µ) := v ∈ v(t, µ, U ) s.t. ∂ - D V (t, µ)(1, v) ≤ 0 , ( 10
)
for all (t, µ) ∈ [0, T ] × P c (R d ). If in addition hypotheses (R) hold, the restricted set-valued maps G |K : [0, T ] × P(K) ⇒ C 0 (K, R d )

Corollary 1 (Continuity inclusion).

A curve of measures µ * (•) is optimal for (P ) if and only if ∂ t µ * (t) ∈ -div x (G(t, µ * (t)µ * (t)), [START_REF] Cannarsa | Semiconcave Functions, Hamilton-Jacobi Equations and Optimal Control[END_REF] with G(•, •) being defined as in [START_REF] Cannarsa | Some Characterizations of Optimal Trajectories in Control Theory[END_REF], namely if and only if there exists an L 1 -measurable selection t ∈ [0, T ] → v * (t) ∈ G(t, µ * (t)) such that the pair (µ * (•), v * (•)) satisfies ∂ t µ * (t) + div x (v * (t)µ * (t)) = 0 in the sense of distributions.

Even though the theory developed in [START_REF] Bonnet | Differential Inclusions in Wasserstein Spaces: The Cauchy-Lipschitz Framework[END_REF] only accounts at the present moment for set-valued velocities that are Lipschitz with respect to the measure variable, the existence of solutions to (11) is guaranteed whenever there exists at least one optimal solution to (P ). In this context if we further assume that hypotheses (R) hold, then for every compact set K ⊂ R d the restricted set-valued maps G |K (•, •) are upper-semicontinuous. The importance of this regularity result comes from the fact that set-valued maps with closed graphs are more convenient when Euler constructions of solutions are involved. Also as previously mentioned, the applications G |K (•, •) are also continuous as soon as they are single-valued, e.g. when the value function is differentiable.

[ 6 ,

 6 Lemma 1] to prove that a norm-semiconcavity estimate similar to (R)-(i) holds for the non-local flows(Φ u (τ,t) [•](•)) τ,t∈[0,T ] .It is then possible to recover the local semiconcavity of V (τ, •) for any fixed τ ∈ [0, T ] by combining the latter estimates with hypotheses (OCP)-(ii),(iii) and (R)-(ii).

  (t))) for any R > 0. By definition of the latter together with the regularity (OCP)-(i) of the controlled velocity field v :[0, T ] × P c (R d ) × U × R d → R d , it can then be shown that for every K ⊂ R d , the graph of the restricted set-valued map G |K (•, •) is closed. Moreover under hypotheses (OCP), it can be noted that G |K (•, •) takes its values in a compact subset of C 0 (K, R d ). By classical results in set-valued analysis (see e.g.[START_REF] Aubin | Set-Valued Analysis[END_REF] Section 1.4]), every such multifunction with closed graph is upper-semicontinuous, which ends the proof of our claim.

  are upper-semicontinuous for every compact set K ⊂ R d .Proof. The proof of the first part of this result -namely that (µ * (•), u * (•)) is optimal if and only if (9) holds L 1 -almost everywhere in [0, T ] -is quite long, and can be obtained by performing computations similar to those in the proof of Theorem 5 above.The regularity of the set-valued map G(•, •) under hypotheses (R) can be established via the following arguments. By Theorem 7, the value function V (•, •) is locally semiconcave. Thus by Proposition 3, the lower-derivative ∂ -

D V (t, µ * (t))(1, v(t, µ * (t), u * (t))) coincides with the regularised Rlower derivative ∂ - o,R V (t, µ * (t))(1, v(t, µ * (t), u *

This material is based upon work supported by the Air Force Office of Scientific Research under award number FA9550-18-1-0254.