Hongliang Lü

David Boilley
email: boilley@ganil.fr

Yasuhisa Abe

Caiwan Shen

Synthesis of superheavy elements: Uncertainty analysis to improve the predictive power of reaction models

Keywords: numbers: 25.70.Jj, 24.10.-i, 24.60.Dr, 06.20.-f

Background: Synthesis of super-heavy elements is performed by heavy-ion fusion-evaporation reactions. However, fusion is known to be hindered with respect to what can be observed with lighter ions. Thus some delicate ambiguities remain on the fusion mechanism that eventually lead to severe discrepancies in the calculated formation probabilities coming from different fusion models.

Purpose: In the present work, we propose a general framework based upon uncertainty analysis in the hope of constraining fusion models.

Method: To quantify uncertainty associated with the formation probability, we propose to propagate uncertainties in data and parameters using the Monte-Carlo method in combination with a cascade code called KEWPIE2, with the aim of determining the associated uncertainty, namely the 95% confidence interval. We also investigate the impact of different models or options, which cannot be modeled by continuous probability distributions, on the final results. An illustrative example is presented in detail and then a systematic study is carried out for a selected set of cold-fusion reactions.

Results: It is rigorously shown that, at the 95% confidence level, the total uncertainty of the empirical formation probability appears comparable to the discrepancy between calculated values.

Conclusions:

The results obtained from the present study provide a direct evidence for predictive limitations of the existing fusion-evaporation models. It is thus necessary to find other ways to assess such models for the purpose of establishing a more reliable reaction theory, which is expected to guide future experiments on the production of super-heavy elements.

I. INTRODUCTION

Experimentally the synthesis of super-heavy elements (SHE) is realized by means of the so-called fusionevaporation reaction which, from a theoretical viewpoint based upon the Bohr independence hypothesis [1], is divided into two separate phases:

• The fusion phase, during which a nucleus-nucleus collision may lead to the formation of an excited heavy nucleus.

• The de-excitation phase, where the newlyformed heavy nucleus has to be chilled by evaporation of light particles, including γ-ray emission against nuclear fission.

For the latter, it is often described using a cascade code.

In our case, we shall employ the KEWPIE2 code [2,3].

Regarding the fusion process, it is decomposed into two successive steps [4]:

• The capture step, corresponding to passing over the Coulomb (or Bass [5]) barrier, after which two colliding nuclei stick together.

• The formation step, from the contact configuration to the compound nucleus (CN). During this step, the system has to overcome an inner potential barrier, which results in the creation of a compound system, where the incident particle combines with the target nucleus and the total energy is shared among all the confined nucleons. Here, it should be mentioned that the formation step was introduced to explain the so-called fusion hindrance phenomenon that only occurs in heavy reaction systems (with a charge product of the projectile-target combination Z p Z t 1600 -1800). Regarding lighter systems, the inner barrier should no longer exist and the colliding nuclei automatically fuse after the capture step.

It should be noted that the main contribution to the fusion hindrance in heavy-ion collisions is due to the formation step, consisting of a diffusion process over an inner potential barrier. Over the past two decades, qualitative features of the formation dynamics have been well understood within the framework of the Langevin equation [6][7][8][9][10][11][12]. This formalism is not unique and some se-rious ambiguities remain on the inner-barrier height and the dissipation strength. Moreover, during an experiment, fusion-fission events that have reached the compound state and quasi-fission ones that re-separate after crossing the Coulomb barrier are usually very difficult to distinguish from each other. As a consequence, this leads to a lack of reliable data on the fusion cross-sections or formation probabilities that could be employed to assess the fusion models. Another delicate issue is that, at the present time, the two separate steps of the fusion process cannot be treated within a unified framework. This remains a severe problem that could eventually affect the final results.

Precise theoretical predictions of evaporation-residue (ER) cross-sections are crucially important for conducting experiments on the synthesis of SHE, because in most cases production cross-sections are extremely low (of the order of one picobarn). For a recent experimental review, see Ref. [13]. A small change of the cross-section could mean months of experiments. To this extent, the first crucial question one might ask is, how exactly can we predict evaporation-residue (ER) cross-sections with current models?

Actually, there have been a number of theoretical calculations on ER cross-sections, which appear to agree remarkably well with the measured ER cross-sections. This, of course, can be regarded as a significant theoretical achievement considering that the measured crosssections usually span at least six orders of magnitude. In the present study, we only focus on the so-called coldfusion reaction, where only one neutron can be evaporated from the compound nucleus. In this case, typical calculations on the ER cross-sections for a selected set of cold-fusion reactions are displayed in Fig. 1 taken from Ref. [14]. The experimental measurements are nicely fitted by almost all models within less than one order of magnitude. This is also the case for other models not included in this comparison such as the one based on a previous version of the code used in this study [15]. Nevertheless, if we take a closer look at the theoretical formation probabilities, as shown in Fig. 2, the calculated values can differ by two or three orders of magnitude, even though all the fusion-evaporation models seem able to reasonably fit the measured data. The large uncertainty in the calculated formation probability is not a surprise because of some serious ambiguities in the reaction mechanism of fusion dynamics, as mentioned before. However, the capture and de-excitation steps are both considered to be better known from a theoretical viewpoint and can be relatively well described on the basis of knowledge about lighter reaction systems. How come the better-known parts can accommodate large discrepancies between the calculated formation probabilities and thus reproduce the measured data? What if the experimental data are not available? To this extent, the predictive capability of the fusion model appears to be quite limited.

To establish a reliable theory for guiding future experiments on the synthesis of SHE, it is necessary to find FIG. 1. [Color online] Theoretical calculations of ER crosssections for the one-neutron evaporation channel of coldfusion reactions. The calculated results can be found in Refs. [10,[16][17][18][19]. The figure is taken from Ref. [14].

FIG. 2. [Color online

] Calculated formation probabilities for the selected set of cold-fusion reactions. The theoretical results can be found in Refs. [10,[16][17][18][19][20]. The figure is adapted from Ref. [14].

new methods to constrain existing fusion models with a special focus on the formation step. Thus, we perform an uncertainty analysis of what are considered as the betterknown parts of the reaction, namely the capture and deexcitation steps, to constrain the lesser-known part. In this work we will include both experimental uncertainties and uncertainties due to the remaining ambiguities in modeling. For the latter, this includes input parameters and models or approximations. The main aim of the present paper is to employ some state-of-the-art statistical methods to quantify different uncertainty sources and look more closely at their impact on the formation probability.

Over the past few decades, the importance of uncertainty analysis in modeling has been increasingly recognized [21], especially with the enormous development of computer simulation. It is mainly concerned with the characterization and quantification of uncertainty in numerical applications by running a great number of trials to investigate the impact of minor differences in the input on final outcomes. This kind of study is essential when one has very poor information about the input. For instance, both the free parameters and the various theoretical models used for capture and evaporation phases are usually not unambiguously determined [2,3].

Before entering into details about uncertainty analysis, let us first briefly present the theoretical models used in this study.

II. THEORETICAL FRAMEWORK

It is commonly known that, based upon the Bohr independence hypothesis [1], the ER cross-section for a fusion-evaporation reaction can be explicitly expressed as

σ ER (E cm ) = Jmax J C ≥0 σ fus (E cm , J C )P sur (E * C , J C ), (1)
recalling that the relationship between incident energy in the center-of-mass frame E cm and total excitation energy of the CN E * C is simply given by E * C = E cm + Q with Q being the Q-value evaluated from ground-state masses. It should be noted that Eq. ( 1) takes into account all partial-wave contributions to the total crosssection, ending up with an upper limit, say J max , which can be approximately determined once the partial-wave cross-section vanishes.

A. A simplified model

In the case of SHE, the maximum spin J max is essentially related to the survival probability P sur (E * C , J C ) which does not vanish only in the vicinity of J C 0. Hence, one has P sur (E * C , J C ) P sur (E * C ), where P sur (E * C ) can be estimated by means of the Weisskopf-Ewing model, and the total ER cross-section can thus be approximated as

σ ER (E cm ) σ fus (E cm )P sur (E * C ). (2)
Moreover, for heavy fusing systems leading to the formation of SHE, the fusion probability needs to be rewritten as a product of the capture and formation probabilities. Accordingly, the above approximate formula can be further transformed into

σ ER (E cm ) σ cap (E cm )P form (E cm )P sur (E * C ), (3)
where P form (E cm ) stands for the formation probability due to the fusion hindrance. We should be aware that this simplified version of Eq. ( 1) has been widely employed to estimate ER cross-sections in many recent studies [10,14,19,20] due to its simplicity. As mentioned in the introduction, we only focus on the one-neutron evaporation channel. In this case, the ER cross-section is simply given by

σ 1n ER σ cap P form P 1n sur . (4)
It should be recalled that there is a huge discrepancy between different theoretical results on P form . On the one hand, it is difficult to reliably measure this quantity as the identification of quasi-fission events remains a delicate issue. On the other hand, as previously mentioned, serious ambiguities in the fusion mechanism still persist. Thus, the formation probability, which is very specific to heavy systems, is considered to be the most poorly known quantity in Eq. ( 4).

In the following we shall attempt to constrain the formation probability by deducing it from experimental data,

P form σ 1n ER σ cap P 1n sur . (5)
In this study, we will estimate the uncertainty of P form and compare it to the theoretical values calculated with various models as displayed in Fig. 2. Note that in the r.h.s. of the above equation, the numerator directly comes from experiments and the denominator from models. For the latter, we shall use the KEWPIE2 code specially dedicated to the synthesis of super-heavy elements [3].

B. Capture cross section

The capture cross-section corresponds to the Coulomb barrier crossing or to the sum of quasi-fission, fusionfission and ER cross-sections. There are several models that can be tested against some measurements. According to a review by Loveland [22] these models are able to describe the magnitudes of the capture cross sections within 50%. However, the ratio of calculated to measured cross sections for various systems and models spans from 0.55 to 2. 24.

In this work we did not perform an extended analysis of the capture cross-section but picked up two typical models that are widely used and implemented as default options in the KEWPIE2 code, namely the empirical barrier-distribution (EBD) method [10,23] and the Wentzel-Kramers-Brillouin (WKB) approximation with a proximity potential [24,25]. In the former model, the barrier distribution is supposed to be a Gaussian with a mean value and width adjusted to reproduce a systematics of 48 reactions, whereas in the latter proximity potential parameters were adjusted to reproduce experimental data without any coupling to other channels.

C. Survival probability

The survival probability of an excited heavy nucleus results from the competition between neutron evaporation and fission that are governed by the neutron-evaporation width Γ n and fission-decay width Γ f , respectively.

On the one hand, the neutron-evaporation width is estimated using the Weisskopf-Ewing theory [26][27][28]:

Γ n (E * C ) = (2s n + 1)µ n π 2 2 E * C -Sn 0 σ n inv ( n )ρ B (E * B ) ρ C (E * C ) n d n , (6)
where C denotes a compound nucleus at excitation energy E * C and B the daughter nucleus at excitation energy

E * B = E * C -S n -n
, where S n is the binding energy of neutrons in the compound nucleus and n the kinetic energy of evaporated neutrons in the center-of-mass frame. µ n is the reduced mass of the reaction system. For neutrons, the cross section for the time-reversed reaction σ n inv is given by

σ n inv ( n ) = g 0 (1 + g 1 n )π R2 , (7)
where g 0 = 0.76 + 1.93A [29]. On the other hand, the fission-decay width is usually calculated within the Bohr-Wheeler (BW) transitionstate method [30]:

-1/3 B , g 0 g 1 = 1.66A -2/3 B -0.05 and R = 1.7A 1/3 B fm
Γ BW f (E * C ) = 1 2πρ gs C (E * C ) E * C -B f 0 ρ sd C (E * sd )d f , (8)
where the excitation energy at the saddle point E * sd is equal to E * C -B ff with f being the kinetic energy of the collective motion. Here, B f represents the fission barrier that is known to have a great influence on the stability of super-heavy elements. A recent review [31] shows that this quantity is difficult to estimate and large discrepancies that can reach up to several MeV still remain in the predicted values.

In this work, we first use a common approximation: B f B LDM -∆E sh with B LDM and ∆E sh being respectively the liquid-drop fission barrier and the shellcorrection energy in the ground state. We also considered a table based upon a so-called microscopic-macroscopic theory [32]. Note that the shell-correction energy is the same in both cases.

From a dynamical point of view, the fission-decay width evaluated by Eq. ( 8) can actually be refined by introducing the Kramers factor [33],

K = 1 + β 2ω sd 2 - β 2ω sd , (9)
and the Strutinsky factor [34],

S = ω gs T gs . (10)
In the above factors, β stands for the reduced friction coefficient that takes into account the effect of viscosity on the fission process. Its value is not well known. ω gs and ω sd denote the potential curvature in the ground state and at the saddle point respectively. They are both fixed at 1.0 MeV. Finally, the total fission-decay width is given by

Γ f = K • S • Γ BW f . (11)
In the following study, the product K • S is referred to as the Kramers-Strutinsky factor. Some models take it into account; some others do not. The state density ρ enters both decay widths in Eqs. ( 6) and ( 8). In our model, the intrinsic state-density formula for a nucleus comprising two kinds of particles, namely protons and neutrons, is explicitly given by [35]

ρ int (E * ) = √ π 12 exp(β 0 E * + a/β 0 ) β 0 E * 3 g 2 0 4g n g p 1/2 × 1 -exp(-a/β 0 ) 1 - 1 2 E * β 0 exp(-a/β 0 ) 1/2 , (12)
where g n and g p are respectively the neutron and proton single-particle state densities at the Fermi energy. Here, g 0 = g n + g p and typically, one has the following approximation: g n g p . In this study, we shall consider various formulas for the level density parameter. In addition to the Tōke-Świątecki model [36] chosen as default value, we also selected Reisdorf and Pomorska et al.'s models [24,37]. Ignatyuk's prescription was also taken into consideration for the energy dependence of the level-density parameter a in the ground state [38]:

a gs (E * ) = a gs 1 + 1 -e -E * /E d ∆E sh E * , (13)
where E d represents the shell-damping factor that characterizes how fast the shell effect is disappearing with excitation energy. Its value is not well established. Moreover, the intrinsic state density should be enhanced by a certain factor as a function of excitation energy E * , namely

ρ(E * ) = ρ int (E * )κ coll (E * ). ( 14
)
For more details on the collective enhancement factor κ coll used in KEWPIE2, the reader is referred to Ref. [3]. This factor is not always implemented in models used to study the survival probability of super-heavy elements.

D. Average values due to the loss of beam energy

Here, we only concentrate on the survival probability with respect to the 1n-channel P 1n sur , which can be written in the following form: In practice, the measured production cross-sections are usually obtained by averaging over a certain range of excitation energy, due to the resolution of the detection setup and the loss of beam energy in the target. When the calculation is confronted with experimental data, this dispersion effect would be crucially important in the case of SHE formed via cold-fusion reactions, because the survival probability is rather peaked, as illustrated in Fig. 3. Therefore, to be close to experimental observables, it is common to take an average over a fixed energy width δE * , as proposed in Ref. [10]. Thus, the average survival probability can be defined as

P 1n sur = P 1n sur (E * , β, E d ). (15
P 1n sur = 1 δE * E * C +δE * /2 E * C -δE * /2 P 1n sur (E * , β, E d )dE * , (16)
which is also dependent on the energy loss δE * in addition to the model parameters. To carry out the integral, we employ the Gauss-Legendre quadrature method, with a number of abscissas chosen to be 12 to get a numerical precision of less than 1.0%.

III. UNCERTAINTY ANALYSIS

In this section, we concretely formulate the problem we are dealing with and specify remaining ambiguities in the input and options of the model.

A. Modeling the physical system

The main purpose of a numerical simulation is to establish an abstract mathematical model for describing some basic features of the physical system under study.

In a general way, a physical system can be abstracted as a set of multivariate real functions, namely

M : R N → R M x → y = M(x), (17)
together with some constraints or restrictions so that the subsequent derivations would make sense. Here, the input parameters of the model can be represented by a N -dimensional vector x ∈ D ⊂ R N , where D denotes the domain of a model function. The vector of M output quantities (or response variables) is thus a Mdimensional vector function. In rare cases, the model might be a simple closed-form function. But more generally, it could correspond to a black-box function, such as a computer program, which requires N input values and then yields some numerical results. It should be mentioned that the model itself can also be changed so that sometimes one has different model functions for the same issue.

In our case, KEWPIE2 is employed to model the better-known parts, namely the capture cross-section and survival probability. The main input parameters for the decay part are nothing else but the reduced friction parameter β and the shell-damping energy E d . It should be noted that the fission barrier has been considered to be model-dependent, instead of a free parameter. In addition to the input parameters, we have also incorporated several sub-models into the code [3].

The main objective here is to extract the empirical formation probability which is simply a scalar-valued model function. Based upon Eq. ( 5), the model function y = M(x) can be constructed. Accordingly, the empirical formation probability is given as follows:

P form = σ 1n exp σ cap P 1n sur = P form (δE * , β, E d , σ 1n exp ), (18)
where the measured data are assumed to be normally distributed. It should be noted that, sometimes, uncertainty intervals associated with the experimental data can be asymmetric. In this case, the method presented here is still valid, provided that the corresponding distribution function is available.

B. Identifying and quantifying uncertainty sources

The common uncertainty sources can be roughly summarized as follows [39]:

• Parameter uncertainty, which comes from the input parameters of the model whose values cannot be exactly inferred, neither theoretically nor experimentally. In our case, for example, the reduced friction coefficient and the shell-damping energy cannot be determined with certainty [38,[40][41][42]. Hence, one needs to construct some input probability distributions for them. • Model uncertainty, which results from the lack of knowledge about the underlying true physics. Thus, it is dependent upon how accurately a model describes the true system in a realistic situation, knowing that models are only rough approximations of reality. In the present study, for example, to estimate the capture cross-section, two approximate methods are employed. Here, to assess model uncertainty, we simply investigate the effects of changing models on the final outcomes.

• Experimental uncertainty, which is directly related to the precision and accuracy of measurements. In simple cases, experimental data can be described by a Gaussian distribution.

• Numerical uncertainty, which arises because of errors corresponding to the implementation of computational methods. We should be aware that most models are too complicated to be solved analytically, so that doing a certain approximation is often necessary. This contribution should be negligible whenever possible.

It should be noted that the first two uncertainty sources are essentially related to theoretical modeling and the distinction between them is sometimes blurred, since a model can also comprise a set of parameters. The next step consists in identifying different uncertainty sources that might come up in various contexts with the aim of quantifying them using probabilistic methods.

Regarding parameters, a sensitivity analysis introduced in our previous studies [2,3] allowed us to extract the most relevant ones, namely the friction coefficient β and the shell-damping energy, E d . Their typical ranges of values can be determined by doing an exhaustive literature search. Since the choice of ranges is somewhat subjective, we would like to introduce two kinds of cases, namely the pessimistic and optimistic hypotheses. In the former one, the ranges of values are summarized as follows:

• The reduced friction coefficient β ∼ 1.0 -9.0 zs -1 .

• The shell-damping energy E d ∼ 13.0 -25.0 MeV.

As its name indicates, such ranges are able to cover most of values that can be found in the literature [38,[40][41][42][43][44][45][46], but the associated uncertainties appear quite large. To justify this choice, Fig. 4 illustrates typical values of β according to previous studies. Here, the points or bars indicate the values extracted from experimental data. It can be seen that the chosen range of β can accommodate most of the extracted values. As regards the shell-damping energy E d , typical values can be found in some recent papers [42,47], which demonstrate the chosen range of values. Accordingly, the default values of the above two parameters are taken to be 5.0 zs -1 and 19.0 MeV, namely their mean values. By contrast, in the optimistic case, the above ranges are expected to be better constrained so that they could become much narrower. In view of this, some smaller ranges of values are worth taking into account, say reduced by 50%, namely

• The reduced friction coefficient β ∼ 3.0 -7.0 zs -1 .

• The shell-damping energy E d ∼ 16.0 -22.0 MeV. This optimistic hypothesis is also inspected in the following study.

Regarding the energy loss in the target δE * , it is assumed to lie in the range 2 -4 MeV [48]. Unlike the excitation energy, uncertainty in the beam energy can be safely neglected because the typical relative uncertainty is only a few percents at most.

To specify the probability density function (PDF) associated with each parameter, one needs to take into account all available information. Here, it should be noted that, in addition to the normalization condition, only the lower and upper limits are available (not that exactly). The maximum entropy principle [49,50] would give a complicated probability distribution as detailed in Ref. [START_REF] Bipm | Joint Committee for Guides in Metrology[END_REF]. However, for the sake of simplicity, we stick here to the uniform probability distribution denoted as U (l min , l max ) in the following study.

In addition to the parameters, it should be recalled that using different models can also have an impact on the final results. This influence needs to be taken into account during analysis.

The experimental data employed in this study are mainly taken from Refs. [START_REF] Münzenberg | [END_REF][53][54] and are assumed to be normally distributed, namely N (µ, σ 2 ) with the mean value µ and the variance σ 2 . Therefore, we are only interested in the ones with symmetric uncertainty bars, even though there are more measured data in the literature. Accordingly, the selected reaction systems and the corresponding maximum production cross-sections are summarized in Tab. I.

As a whole, it is clearly demonstrated that the maximal production cross-section drops considerably with increasing atomic number from Z = 104 to Z = 108, whereas the experimental uncertainty gradually rises due to increasing difficulty in performing such measurements.

C. Uncertainty propagation

To propagate input uncertainties through the model so as to quantify the output uncertainty, we use a Monte-Carlo simulation as recommended in Ref. [START_REF] Bipm | Joint Committee for Guides in Metrology[END_REF], for calculating the mean value µ Y and the standard deviation u(Y ) associated with the model response. More concretely, supposing that one has a set of samples for the input random vector, namely {x 1 , x 2 , . . . , x S }, the random model function is then evaluated for each of them as follows:

y r = M(x r ), r = 1, ..., S. (19)
Thus, the usual statistical estimators of these quantities are simply given by

µ Y = 1 S S r M(x r ), (20)
u 2 (Y ) = 1 S -1 S r [M(x r ) -µ Y ] 2 . ( 21
)
As regards the sample size, it can be approximately determined by checking the stability and convergence of the calculated result. At a 95% confidence level, the value of S is determined to be 40 000, which would give a relative numerical uncertainty of less than 1%. In the following study, the confidence level is kept fixed at 95% as usual.

To obtain the corresponding confidence interval, first, one has to estimate the cumulative distribution function (CDF) related to the output quantity, namely

F Y (y) = y -∞ f Y (t)dt. ( 22
)
The estimated CDF, denoted by FY (y), can be obtained as follows [START_REF] Bipm | Joint Committee for Guides in Metrology[END_REF]:

• Sorting the values {y r } of the output quantity provided by the Monte-Carlo simulation into increasing order. The sorted values are denoted by {y (r) }.

• Assigning uniform cumulative probabilities p r = (r -1/2)/S to the ordered values.

Finally, once the estimated CDF has been constructed, it is possible to determine the endpoints that define the required confidence interval. In the case of symmetric output distributions, the 95% confidence interval is determined by the 0.025-and 0.975-quantiles. If the output distribution is asymmetric the shortest confidence interval should be adopted. It generally does not match with the 0.025-and 0.975-quantiles anymore and can be obtained numerically from the estimated CDF. For details on this approximation and the practical algorithm, the reader is referred to Ref. [START_REF] Bipm | Joint Committee for Guides in Metrology[END_REF].

D. Summary

In the following, we calculate the empirical formation probability, P form (δE * , β, E d , σ 1n exp ) defined in Eq. ( 18), and perform an uncertainty analysis that is essentially based upon the recommendations given in Ref. [START_REF] Bipm | Joint Committee for Guides in Metrology[END_REF] which serves as the international standard for expression of uncertainty.

Regarding the PDFs for input parameters, they are summarized as follows:

• σ 1n exp is assumed to be normally distributed. • δE * ∼ U (2, 4) in MeV. Note that over each value of δE * , the survival probability is averaged (cf. Eq. ( 16)).

• β ∼ U (1.0, 9.0) in zs -1 .

• E d ∼ U (13.0, 25.0) in MeV.

As previously discussed, the above ranges of values of the last two parameters are considered under the pessimistic hypothesis, which simply means that they seem somewhat large but can cover most of the parameter values that can be found in the literature. Moreover, it would be instructive to take a look into the effects of using some narrower ranges of values on the final results. This is the optimistic case, where the values of β and E d are supposed to lie within the following ranges:

• β ∼ U (3.0, 7.0) in zs -1 .

• E d ∼ U (16.0, 22.0) in MeV.

In addition to the input parameters, the effects of changing model functions on the calculated results should also be examined. Here, we mainly focus on some theoretical corrections, namely the collective enhancement and the Kramers-Strutinsky factors, as well as the fissionbarrier models, the level-density parameters and two simple methods for estimating capture cross-sections, which are considered to be the main physical ingredients included in the model.

It should be recalled that the default parameter values and models were concretely presented in a recent paper on the description of KEWPIE2 [3]. To briefly summarize, in the present study, the empirical barrierdistribution method [10,55] for estimating the capture cross-section, the Tōke-Świątecki model [36] for the leveldensity parameter, the Thomas-Fermi model [56] for calculating fission barriers, and the collective enhancement and Kramers-Strutinsky factors [33,34,57,58] have been considered by default.

The Monte-Carlo approach is employed to propagate input distributions through the model, as the uncertainty associated with each of the input parameters is considerably large. In the present work, the GSL scientific library [START_REF] Gough | GNU scientific library reference manual[END_REF] was employed to generate uniform and normal random numbers for input distributions.

IV. RESULTS AND DISCUSSION

In this section, we start by showing an illustrative example for the cold-fusion reaction leading to the formation of element Hs, namely 208 Pb( 58 Fe,1n) 265 Hs, at its optimum energy. Then, a systematic study is carried out and discussed. A. Influence of input distributions Fig. 5 illustrates the estimated distribution related to the empirical formation probability obtained by means of the KEWPIE2 code together with available production cross-sections (cf. Eq. ( 18)). It should be mentioned that all the input parameters were considered under the pessimistic hypothesis. In this case, the shape of the distribution appears to be highly skewed to the left, that is, it has a long tail on the right hand side of the distribution and thus results in a quite asymmetric confidence interval. The mean value is hence located on the right-hand side of the peak value. Fig. 5 also displays the corresponding CDF used to determine both lower and upper bounds of the 95% confidence interval.

Separate uncertainty contributions coming from each of the input distributions are displayed in Fig. 6. Each point stands for the mean value of the empirical formation probability, together with its 95% confidence interval. First, the contribution of the experimental uncertainty is represented by the first point on the left. The upper bound of the confidence interval is higher than the lower one by a factor of about 2. The following three points are connected with the energy loss in the target, the reduced friction parameter β and the shell-damping energy E d , respectively. Overall, it can be clearly noticed that the theoretical contributions are dominant compared to the experimental ones. This is simply due to the fact that β and E d are closely related to the fission process that is known to be the dominating decay channel for SHE. It should also be mentioned that the relative uncertainty of β is about 46%, which appears to be much larger than that of E d estimated to be about 18%. Put differently, this means that uncertainty related to the empirical formation probability would be more sensitive to the shell-damping energy. Finally, when including all input distributions, the total confidence interval spans slightly more than one order of magnitude.

When reducing the intervals of the two dominant parameter values by 50%, that is, under the optimistic hypothesis, Fig. 6 clearly demonstrates that the total uncertainty is also decreased by almost the same factor. This result also indicates that further constraints on these two critical model parameters would be able to reduce the total uncertainty connected with the empirical formation probability. We then look into the impact of changing models on the empirical formation probability. Fig. 7 illustrates how correction factors, namely the Kramers-Strutinsky and collective enhancement ones, affect the mean value of the empirical formation probability. Removing the Kramers-Strutinsky factor, it is readily seen that the mean value slightly increases by a factor of about 2, whereas it becomes a bit lower after removing the collective enhancement factor. It should also be noted that the uncertainty amplitudes do not remain the same due to the fact that β is jointly removed along with the Kramers-Strutinsky factor so that the output distribution becomes narrower.

B. Influence of models

Fig. 8 shows the effects of changing level-density parameter formulas on the mean value of the empirical formation probability and the associated confidence interval. As already mentioned, in addition to the default model of Tōke and Świątecki [36], we also selected the Reisdorf [24] and Pomorska et al. models [37]. Here, one can notice that the mean value of the empirical probability decreases at most by a factor of around 2, which appears to be negligible compared to the total uncertainty interval. As a whole, the lengths of their confidence intervals seem to remain the same. and confidence interval of the empirical formation probability. Here, the empirical barrier-distribution (EBD) method [10,23] and the WKB approximation [24,25], were considered. As can be seen in Fig. 9, the results based upon these two simple capture models seem to differ by a factor of about 7, that is, less than one order of magnitude. However, special attention should always be drawn to the capture step because calculated crosssections could differ by one order of magnitude, especially in the case of deformed colliding nuclei due to nuclear structure effects. Fig. 10 demonstrates how the mean value and the confidence interval associated with the empirical formation probability evolve with various fission-barrier models. Here, three different models were considered. Two of them are based on a commonly used approximation with a macroscopic part calculated with the Thomas-Fermi (TF) model [56] or the Lublin-Strasbourg Drop (LSD) [START_REF] Ivanyuk | [END_REF] to which a shell-correction table [61] is added. The third model is Möller et al's microscopic-macroscopic model [32]. It is thus observed that the mean values can differ by almost one order of magnitude that is comparable to the total uncertainty interval. This is mainly due to the fact that nuclear fission is known to be the dominant decay channel for heavy nuclei. A small change in fission barrier can give rise to a significant variation in the calculated survival probability. It is worth noticing that if we had taken into account all fission barriers available in the scientific literature, the impact would have been larger.

C. An extended comparison

So far, we have illustrated an uncertainty analysis for a special cold-fusion reaction. In this subsection, an extended study is presented for the whole set of cold-fusion reactions listed in Tab. I. Two different cases related to the empirical formation probabilities and the associated 95% confidence intervals for the selected cold-fusion reactions leading to the synthesis of SHE ranging from Z = 104 to Z = 108 are shown in Figs. 11 and12. The former displays calculated results with two simple capture models and the latter shows those based upon two extreme fission-barrier models. It should be kept in mind that the uncertainty interval is always related to all input distributions within the pessimistic hypothesis.

First, one can immediately notice that the discrepancy due to the capture step can even reach two orders of magnitude for the lightest reaction system, but it gradually declines as the atomic number of the compound system goes up. This is due to the fact that, as the system becomes heavier, the incident energy for the one-neutron channel is approaching to the Coulomb barrier so that uncertainty associated with the capture model because of sub-barrier coupled-channels effects should be gradually reduced. To this extent, the uncertainty remaining in the estimated capture cross-section could also have a significant impact on the final results, but fortunately it can be directly measured and thus be constrained. Here, it should be kept in mind that we only considered two simple capture models in our calculations, whereas many other more sophisticated and accurate ones are currently available. Accordingly, one would expect that uncertainty associated with the capture model can be reduced either theoretically or experimentally.

Regarding the impact of fission barriers, as displayed in Fig. 12, the discrepancy practically seems to remain the same. Overall, the mean values based upon two extreme fission-barrier models with the same shell correction energy can differ by about one order of magnitude, which is comparable to the lengths of their corresponding uncertainty intervals. Unlike in the case of capture models, since the fission-barrier height cannot be directly measured, it is unlikely that uncertainty would be reasonably reduced. To this extent, other methods should to be invented so as to constrain fission-barrier models. For instance, we recently started to consider Bayesian inference together with its possible application [63,64].

Finally, it is also interesting to take a look into the extreme cases, that is, with two capture models and two extreme fission-barrier models. This is displayed in Fig. 13. As for some lighter systems, the maximal discrepancy can reach more than three orders of magnitude, whereas it gradually diminishes to around two orders of magnitude as the atomic number goes up. Considering the total uncertainty interval, the discrepancy between the calculated formation probabilities can be completely accommodated and hence, it appears unlikely to discriminate the various formation models. Put differently, as mentioned at the beginning of this paper, all fusion-evaporation models would be capable of reproducing the measured data, even though there are still serious ambiguities in the reaction mechanism for heavy-ion fusions.

V. CONCLUSION AND PERSPECTIVES

Although various models converge to experimental data when estimating excitation function of ER crosssections of super-heavy nuclei, a closer look at intermediate steps shows large discrepancies between them. It is especially the case for the formation step that is the lesser known part of the reaction. Consequently, when experimental data are not available, predictions diverge. The results obtained from the uncertainty analysis presented in this paper are helpful for comprehending the contradiction between theoretical formation probabilities and ER cross-sections. It was demonstrated that, although the formation probability is not yet quantitatively known, by somehow tuning the better-known factors, namely the capture cross-section and survival probability, one can still be able to obtain a reasonable fit to experimental data.

In this sense, the predictive power of fusionevaporation models remains quite limited in the case of super-heavy element synthesis. To assess them, it would be necessary to refine both the capture and fission-barrier models. In the former case, as previously mentioned, this can be done with the help of more accurate data as the capture cross-section can be directly measured, or using some more sophisticated capture models. In the latter case, however, there are still severe discrepancies between various types of calculations [31] and moreover, the fission-barrier height cannot be directly measured. Its extraction from experimental data is model dependent [63,64].

How to constrain both formation models and fission barriers simultaneously? There is a crucial need for experiments dedicated to a better understanding of the fusion hindrance. The synthesis of the same nucleus by hindered and non-hindered reactions would be helpful as the decay part would be the same. From a theoretical point of view, there is a real need to decrease the number of free parameters. Thus, fission barriers and inner fusion barriers responsible for the fusion hindrance should be calculated within the same model. This generally requires microscopic or microscopic-macroscopic approaches [31].

VI. ACKNOWLEDGMENT

The authors acknowledge JSPS and the French-Japanese International Associated Laboratory for Nuclear Structure Problems (LIA-FJNSP) for partial support of this work. Part of this work is supported by NSF

FIG. 4 .

 4 FIG. 4. Typical values of the reduced friction parameter that can be found in the literature. The figure is adapted from Ref. [41].

FIG. 5 .

 5 FIG. 5. [Color online] Estimated distribution of the empirical formation probability for the reaction 208 Pb( 58 Fe,1n) 265 Hs. The solid line represents the mean value and the dotted lines refer to the lower and upper bounds of the 95% confidence interval, respectively. The cumulative distribution function (CDF) of the empirical formation probability used to define the confidence interval is indicated.

FIG. 6 .

 6 FIG. 6. [Color online] Impact of input distributions on uncertainty associated with the empirical formation probability. The pessimistic case is considered in the upper panel and the optimistic input distributions for β and E d in the lower one. See text for details.

FIG. 7 .

 7 FIG. 7. [Color online]Impact of correction factors on the empirical formation probability. Note that the uncertainty interval is due to all parameters under the pessimistic hypothesis.

  FIG. 8. [Color online]Impact of level-density parameters on the empirical formation probability. Uncertainty interval is due to all parameters under the pessimistic hypothesis.

  FIG. 9. [Color online] Impact of capture models on the empirical formation probability. Uncertainty interval is due to all parameters under the pessimistic hypothesis.

Fig. 9

 9 FIG. 10. [Color online] Impact of fission barriers on the empirical formation probability. Uncertainty interval is due to all parameters under the pessimistic hypothesis.

  )

			Influence of energy loss in the target
	0.010					Without average
							δE* = 2 MeV
	0.008						δE* = 3 MeV δE* = 4 MeV
	0.006 sur P -1n
	0.004
	0.002
	0.000	6	8	10	12	14	16	18	20
					E* [MeV]
	FIG. 3. [Color online] Average of the survival probability
	for the 1n-channel over different energy intervals to take into
	account the loss of beam energy in the target. The solid
	red curve indicates the calculated survival probability without
	taking an average.

TABLE I .

 I Experimental data for the 1n-channel with the optimum energy and the corresponding maximal ER crosssection. The data are taken from Refs.[START_REF] Münzenberg | [END_REF][53][54].

	Reaction	ZC AC Ecm [MeV] max σ 1n exp [pb]
	208 Pb( 50 Ti,1n) 257 Rf 104 258	185.02	10 419 +1 284 -1 284
	209 Bi( 50 Ti,1n) 258 Db 105 259	187.50	2 200 +240 -240
	208 Pb( 54 Cr,1n) 261 Sg 106 262	202.01	2 520 +253 -253
	209 Bi( 54 Cr,1n) 262 Bh 107 263	205.81	163 +34 -34
	208 Pb( 58 Fe,1n) 265 Hs 108 266	219.27	69 +12 -12

of China under grant number 11275068 and 11547312. The authors greatly appreciate the warm hospitality and support provided by GANIL and RCNP at Osaka Univ., which enabled them to work together. H. L. would like to thank the staff at Huzhou Univ. for their hospital-ity. This work is part of H. L.'s Ph.D. thesis which was supported by the Centre National de la Recherche Scientifique (CNRS) and the Conseil Régional de Basse-Normandie.


Florent Altché
email: florent.altche@mines-paristech.fr

Xiangjun Qian
email: xiangjun.qian@mines-paristech.fr

Arnaud De La

Arnaud De

La Fortelle

Least Restrictive and Minimally Deviating Supervisor for Safe Semi-Autonomous Driving at an Intersection: An MIQP Approach

published or not. The documents may come

I. INTRODUCTION

While self-driving vehicles are not likely to reach consumer market before a few more years, most of the technology required for semi-autonomy already exists. For instance, adaptive cruise control [START_REF] Vahidi | Research advances in intelligent collision avoidance and adaptive cruise control[END_REF] or advanced emergency braking [START_REF] Coelingh | Collision Warning with Full Auto Brake and Pedestrian Detection -a practical example of Automatic Emergency Braking[END_REF] are becoming standard equipment in new cars to assist, and sometimes supplant, human drivers in order to improve safety and driving comfort. Such systems are efficient to prevent rear-end collisions, but they are not suited to more complex situations such as intersections, which represented more than 20% of road fatalities in the United States in 2014.

Until recently, most studies considered fully autonomous driving at intersections, assuming that vehicles were perfectly controlled [START_REF] Dresner | Multiagent traffic management: a reservation-based intersection control mechanism[END_REF], [START_REF] De | Analysis of reservation algorithms for cooperative planning at intersections[END_REF], and able to exchange their future intentions and to quickly and cooperatively react to any unexpected event. Some methods do exist to allow autonomous cars to handle unexpected maneuvers from human-driven vehicles [START_REF] Qian | Prioritybased coordination of autonomous and legacy vehicles at intersection[END_REF], but the safety of the human-driven vehicle itself is rarely considered.

In this article, we consider a road intersection with semiautonomous vehicles equipped with communication capability, driving along fixed paths through the intersection (see, Fig. 1. Example of paths in a general intersection with three-lane roads. The outer (dotted) square is the boundary of the supervision area, and the inner square defines the collision area. e.g., Fig. 1). At a given time, each vehicle is assumed to have a preferred desired control input for the next few fractions of a second: it can simply consist of tracking a desired speed or a more complicated human input. As vehicles are not fully autonomous, this target control may, or may not lead the system of vehicles into a so-called "inevitable collision state" [START_REF] Fraichard | Inevitable Collision States. A Step Towards Safer Robots?[END_REF], or in a deadlock configuration where no participant is able to cross the intersection. Therefore, an outside input is needed to ensure safe and non-blocking coordination of the vehicles across the intersection. Note that the fixed-paths assumption is classic in intersection problems, as vehicles generally drive inside a preassigned lane.

Colombo et al. [START_REF] Colombo | Least Restrictive Supervisors for Intersection Collision Avoidance: A Scheduling Approach[END_REF] introduced the idea of a supervisory instance (called supervisor) tasked with preventing the system of vehicles from entering undesirable states by overriding the controls of one or several vehicles only when necessary. The question of determining whether overriding is needed or not, called "verification problem", is NP-hard [START_REF] Reveliotis | On the Complexity of Maximally Permissive Deadlock Avoidance in Multi-Vehicle Traffic Systems[END_REF]. Under several simplifying assumptions, this verification problem is shown in [START_REF] Colombo | Least Restrictive Supervisors for Intersection Collision Avoidance: A Scheduling Approach[END_REF] to be equivalent to solving a scheduling problem, by assigning time slots during which only one vehicle is allowed inside the intersection. This allows to design a socalled least restrictive supervisor, which verifies the safety of the desired inputs and overrides them if necessary. However, the proposed supervisor only works with simple intersection geometries with one conflict point, i.e. only one vehicle can be inside the intersection at any given time. Moreover, no bound is given on the controls used for overriding, which can widely deviate from the desired ones. We call a supervisor minimally deviating if it minimizes the difference between the desired and overridden controls.

Several variations have been proposed based on the equiv-alence demonstrated in [START_REF] Colombo | Least Restrictive Supervisors for Intersection Collision Avoidance: A Scheduling Approach[END_REF]. Reference [START_REF] Bruni | Robust multi-agent collision avoidance through scheduling[END_REF] designs a supervisor which is robust to bounded uncertainties by adding safety margins. Reference [START_REF] Ahn | Semi-autonomous Intersection Collision Avoidance through Job-shop Scheduling[END_REF] leverages job-shop scheduling to develop a supervisor that considers several possible conflict points inside the intersection; however, vehicles are modeled as single integrators which is not realistic in a realworld setting. Campos et al. [START_REF] De Campos | Optimal and least restrictive supervisory control: Safety verification methods for humandriven vehicles at traffic intersections[END_REF] proposed a Pareto-optimal supervisor leading to a minimally deviating formulation by recursively finding the most constrained vehicle, reserving its optimal crossing time, and scheduling the crossing of the other vehicles using the previously scheduled ones as constraints. This method allows to minimize the deviation between the overridden and desired controls, but may be computationally intensive.

Inspired by previous work on using operations research for autonomous intersection management [START_REF] Peng | Coordinating Multiple Robots with Kinodynamic Constraints Along Specified Paths[END_REF] and studies demonstrating the interest of using pairwise priorities between vehicles [START_REF] Gregoire | Priority-based coordination of mobile robots[END_REF], we propose a novel design of a least restrictive and minimally deviating supervisor. Our main contribution is a new formulation of the supervisory problem as a finite horizon Model Predictive Control (MPC) problem, in which the safety constraints are modeled as mixedinteger linear inequalities on vehicles control and a quadratic objective function is used to penalize the deviation from the vehicles target control inputs. This formulation, extending our previous work on time-optimal coordination of fully autonomous vehicles [START_REF] Altché | Time-optimal Coordination of Mobile Robots along Specified Paths[END_REF], results in a mixed-integer quadratic program (MIQP) that can be solved using widely available software. Compared to previous approaches [START_REF] Colombo | Least Restrictive Supervisors for Intersection Collision Avoidance: A Scheduling Approach[END_REF], [START_REF] Bruni | Robust multi-agent collision avoidance through scheduling[END_REF], [START_REF] Ahn | Semi-autonomous Intersection Collision Avoidance through Job-shop Scheduling[END_REF], the advantage of our formulation is twofold: 1) intersections with complex geometries (crossing, merging and diverging paths) can be handled intuitively, 2) the computational burden can be significantly reduced as we only consider a finite time horizon, and use highly efficient MIQP solvers.

The rest of this paper is articulated as follows. In Section II, we describe the modeling of the system of vehicles and, inspired by [START_REF] Colombo | Least Restrictive Supervisors for Intersection Collision Avoidance: A Scheduling Approach[END_REF], [START_REF] De Campos | Optimal and least restrictive supervisory control: Safety verification methods for humandriven vehicles at traffic intersections[END_REF], formally define a "Supervision problem" and its optimal, minimally deviating version. In Section III, we present an infinite horizon MPC formulation which solves the discretized optimal supervision problem. In Section IV, we deduce an equivalent, finite time horizon MPC problem which guarantees infinite horizon safety and deadlock avoidance. These results are demonstrated in Section V using computer simulation; finally, Section VI concludes the study.

II. PROBLEM STATEMENT

A. System modeling

We consider the problem of safely coordinating a set N of N semi-autonomous vehicles at a generic road intersection, which could be comprised of multiple entry and exit lanes. Contrary to full autonomy, we assume that semi-autonomous vehicles are able to safely navigate outside intersections, for instance with adaptive cruise control techniques, but need a form of coordination to safely cross intersections. Therefore, we consider a finite supervision area comprising the intersection and a portion of the roads leading to, and exiting from the intersection, and only consider possible collisions inside this area. In the supervision area, vehicles are assumed to follow predetermined paths with good lateral accuracy. Fig. 1 shows an example intersection between two roads, each one comprising three entry and exit lanes, and the corresponding possible paths inside the intersection.

We model the longitudinal dynamics of each vehicle as a double integrator: for a vehicle i ∈ N following a path γ i , the curvilinear position s i , longitudinal speed v i and control input (or acceleration) u i are linked as

ṡi vi = v i u i . (1)
We choose the origin of s i so that s i = 0 when the front bumper of vehicle i enters the supervision area, and s i = s out i > 0 when its rear bumper exits the supervision area. Vehicles are assumed to only move forward along their paths, so that for all i ∈ N , v i ≥ 0. Furthermore, we consider bounds on the speed and acceleration of each vehicle, noted

v i ∈ V i = [0, v i ] and u i ∈ U i = [u i , u i ]
where u i < 0 < u i . We let x i = (s i , v i ) be the state of vehicle i, and we use the bold characters x = (x i ) i∈N and u = (u i ) i∈N to respectively denote the vectors of states and control inputs for the system of all vehicles.

We consider complex intersection geometries with merging, diverging, crossing and non-conflicting paths, in which both rear-end and side collisions may occur. For two vehicles i, j with a risk of rear-end collision, we note S ij = s ij , s ij the interval of positions s i where i can have rear-end collisions with j. For vehicles on the same path, S ij = [0, s out i ]; for vehicles on paths diverging at the curvilinear position

s diverge ij , S ij = [0, s diverge ij ]; for vehicles on merging paths, S ij is in the form [s merge ij , s out ij ].
If there is no risk of rear-end collision for i and j, we let S ij = ∅.

We call collision area the center of the intersection where side collisions can occur, as presented in Fig. 1. Note that side collisions cannot occur for vehicles on diverging paths, but can occur for merging paths. The portion of the path of a vehicle i comprised in the collision area is called the collision segment for vehicle i, and we note S ⊥ i = s ⊥ i , s ⊥ i the set of curvilinear coordinates corresponding to this segment. For two vehicles i and j, we let C ij be the set of positions (s i , s j ) ∈ S ⊥ i × S ⊥ j where a side collision between i and j can happen, and we say that two vehicles i and j are conflicting if C ij = ∅. To ensure a safe crossing of the intersection, we require that two conflicting vehicles cannot be inside their collision segments simultaneously. However, to reduce unnecessary delays and increase traffic efficiency, non-conflicting vehicles are allowed to be inside their collision segments at the same time.

At a given time τ , we say that the corresponding system state x(τ ) is safe if there exists a control t → u(t) ∈ i∈N U i defined over [τ, +∞[ corresponding to a safe infinite-horizon system trajectory t → x(t) over [τ, +∞[, i.e. for which no rear-end or side collision happens at any future time. On the contrary, if no such control exists, the state x(τ ) is called an inevitable collision state [START_REF] Fraichard | Inevitable Collision States. A Step Towards Safer Robots?[END_REF]. Formally, an infinite horizon trajectory is safe if, ∀t ≥ τ, ∀i, j ∈ N :

C ij = ∅ ⇒ (s i (t), s j (t)) / ∈ S ⊥ i × S ⊥ j (2) (s i (t), s j (t)) ∈ S ij × S ji ⇒ (s i (t) ≥ s j (t) + d ij ∨ s j (t) ≥ s i (t) + d ji ) (3)
where d ij and d ji account for a safe following distance and a potential offset of longitudinal positions between vehicles. Note that condition (3) corresponds to either i following j or j following i, accounting for the case of vehicles on merging paths, for which a relative order is not set in advance.

B. Discretized dynamics

At a given time, each vehicle has a certain desired control u i,des , for instance taking into account passengers comfort, energy consumption or time constraints. The vector of these desired controls defines a desired system control u des . Depending on the dynamic constraints and the current system state, this desired control may lead the system into an unsafe situation, i.e. a state where a collision is inevitable.

Following the ideas introduced in [START_REF] Colombo | Least Restrictive Supervisors for Intersection Collision Avoidance: A Scheduling Approach[END_REF], we define a Supervisor as an instance working in discrete time, which assigns a control to each vehicle to ensure that the system remains in a safe state for all future times. We let ∆t > 0 be the duration of a time step for the supervisor, and we assume that individual vehicles update their control in discrete time with the same time step duration, so that the dynamic equations for a vehicle i can be rewritten as:

s k+1 i = v k i ∆t + 1 2 u k i ∆t 2 (4)
v k+1 i = v k i + u k i ∆t (5)
where s k i , v k i and u k i are respectively the position, speed and control of the vehicle at the beginning of the time interval [k∆t, (k +1)∆t[. We note

x k i = (s k i , v k i )
the state of a single vehicle, and we let x k = (x k i ) i∈N and u k = (u k i ) i∈N . As in the continuous case, we formally introduce the safety of the discretized system: at a time instant κ, we say that the current system state x κ is safe if there exists a piecewiseconstant control (u k ) k≥κ that ensures the existence of a safe infinite-horizon system trajectory (x k ) k≥κ . Note that directly discretizing conditions (2) and ( 3) is not sufficient to ensure the infinite horizon safety of the discretized system. Therefore, we adopt the following sufficient conditions for system safety, adapted from ( 2) and ( 3), and we require that for all k ≥ κ and all i, j ∈ N ,

C ij = ∅ ∧ s k i ∈ S ⊥ i ⇒ s k+1 j / ∈ S ⊥ j (6) (s k i , s k j ) ∈ S ij × S ji ⇒ s k+1 i ≥ s k+1 j + d ij ∨ s k+1 j ≥ s k+1 i + d ji (7) (s k i , s k j ) ∈ S ij × S ji ∧ s k+1 i ≥ s k+1 j + d ij ⇒ s k+1 i -s k+1 j + ∆t 2 v k+1 i -v k+1 j ≥ d ij (8)
The additional condition [START_REF] Reveliotis | On the Complexity of Maximally Permissive Deadlock Avoidance in Multi-Vehicle Traffic Systems[END_REF] ensures the absence of longitudinal collisions between two time steps.

C. Supervision problem

At the beginning of a time step κ and for an initially safe system state x κ , we note u κ des the desired (constant) system control over [κ∆t, (κ + 1)∆t[. The Supervisor is tasked with returning a safe (constant) control u κ saf e for the set of vehicles so that the next system state x κ+1 is safe. In what follows, we note U κ saf e ⊂ i U i the set of safe constant system controls over [κ∆t, (κ + 1)∆t[.

As introduced in [START_REF] Colombo | Least Restrictive Supervisors for Intersection Collision Avoidance: A Scheduling Approach[END_REF] we define the (least restrictive) Supervision Problem (SP ) as follows:

Problem 1: SP . For a safe system state x κ and a desired system control u κ des , return u κ saf e = SP (x κ , u κ des ) such that u κ saf e ∈ U κ saf e and so that u κ saf e = u κ des if u κ des ∈ U κ saf e . In the supervised driving paradigm, the control actually implemented by the vehicles over [κ∆t, (κ + 1)∆t[ is u κ saf e . The supervision problem SP is least restrictive in the sense that an "overriding" only occurs when the desired system input would steer the system into an inevitable collision state. The fact that the state reached using u κ saf e is safe ensures the recursive feasibility of SP .

We now define the optimal Supervision Problem (SP * ), which is a generalization of [START_REF] De Campos | Optimal and least restrictive supervisory control: Safety verification methods for humandriven vehicles at traffic intersections[END_REF]. We say that a supervisor solving SP * is a minimally deviating supervisor for a given norm || • || over the set of bounded system controls Problem 2: SP * . For a system state x κ and a desired system control u κ des , return u * κ saf e = SP * (x κ , u κ des ) such that u * κ saf e = arg min u∈U κ saf e ||u κu κ des ||. Note that the optimal supervision problem is an instance of the least-restrictive supervision problem, in the sense that any solution to SP * is a solution to SP .However, not all solutions of SP are solutions of SP * . Previous literature [START_REF] Colombo | Least Restrictive Supervisors for Intersection Collision Avoidance: A Scheduling Approach[END_REF], [START_REF] Bruni | Robust multi-agent collision avoidance through scheduling[END_REF], [START_REF] Ahn | Semi-autonomous Intersection Collision Avoidance through Job-shop Scheduling[END_REF] mainly considers SP , by converting it to an equivalent scheduling problem. Campos et al. [START_REF] De Campos | Optimal and least restrictive supervisory control: Safety verification methods for humandriven vehicles at traffic intersections[END_REF] proposed to exhaustively enumerate all solutions of SP to find the solution to SP * , which may be computationally costly. In this paper, we propose an MPC-based approach that can solve SP * efficiently, even with a large number of vehicles and complex intersection geometries.

III. INFINITE HORIZON MPC FORMULATION FOR SP *

In this section, we present an infinite horizon MPC formulation to exactly solve the optimal supervision problem SP * .

For two conflicting vehicles i and j, we introduce the binary decision variable π ij which is equal to 1 if i goes before j, and 0 otherwise; we require that π ij + π ji = 1. Moreover, we define a set of indicator variables ε such that ε ⊥,in

i (k) = 1 if s k i ≥ s ⊥ i and ε ⊥,out i (k) = 1 if s k i ≥ s ⊥ i . Similarly, we let ε ,in ij (•, k) = 1 if s k • ≥ s ij (•) and ε ,out ij (•, k) = 1 if s k • ≥ s ij (•) for • = i or j.
Reference [START_REF] Altché | Time-optimal Coordination of Mobile Robots along Specified Paths[END_REF] shows how constraints (2) and (3), as well as the constraints on the indicator variables ε can be enforced using linear inequalities involving continuous and integer variables using a "big-M" formulation. To ensure safety and deadlock avoidance over an infinite time horizon we require, for every pair of conflicting vehicles i and j and for all time steps k:

π ij = 1 ∧ ε ⊥,out i (k) = 0 ⇒ ε ⊥,in j (k + 1) = 0 (s1) π ij = 1 ∧ ε ,in ji (k) = 1 ∧ ε ,out ij (k) = 0 ⇒ s k+1 i -s k+1 j ≥ d ij (s2) π ij = 1 ∧ ε ,in ji (k) = 1 ∧ ε ,out ij (k) = 0 ⇒ s k+1 i -s k+1 j + ∆t 2 v k+1 i -v k+1 j ≥ d ij (s3)
Note that (s1), (s2) and (s3) respectively correspond to ( 6), ( 7) and [START_REF] Reveliotis | On the Complexity of Maximally Permissive Deadlock Avoidance in Multi-Vehicle Traffic Systems[END_REF]. Conditions (s2) and (s3) are used for vehicles on merging paths. For non-conflicting vehicles having following constraints (i.e. on parallel or diverging paths), similar conditions are required but the value of π ij is implicitly chosen based on the initial positions of the vehicles.

Any trajectory verifying the above conditions (for all vehicles for which they are applicable) at all time steps is safe over an infinite time horizon. To ensure trajectories are dynamically feasible for the vehicles, we also enforce constraints ( 4) and ( 5), as well as bound constraints

v k i ∈ [0, v i ] and u k i ∈ [u i , u i ]
for all vehicles and at all time steps. As a result, any control verifying all the above conditions is a member of U κ saf e . To choose the optimal control (in the sense of SP * ) in U κ saf e , we introduce the following objective function:

O κ = i ∈ N w κ i u κ i,des -u κ i 2 (9)
where the constant w κ i ≥ 0 is a weighting coefficient for vehicle i at step κ which can be used, for instance, to distinguish high priority vehicles. From the previous considerations and noting X = (X k ) k≥κ where X k is the tuple of the variables defined above at step k, we deduce the following theorem about our Infinite Horizon Model Predictive Control (IH-MPC) formulation:

Theorem 1: Let u κ be a solution of the optimization problem:

u κ = Π u arg min X O κ (X) (IH-MPC)
subj. to ∀k ≥ κ, (4), ( 5), ( s1)

-(s3)
where Π u is a projection operator consisting of only keeping the components of the solution corresponding to u κ . Then u κ is a solution to the optimal supervision problem SP * at step κ for a (semi-)norm of R N . Proof: We note that any solution to the above problem ensures safety for all vehicles at all future times; therefore, the corresponding control u κ is in U κ saf e and is thus a solution to SP . Moreover, by construction the solution minimizes O κ , thus it is a solution to SP * for the seminorm associated with the bilinear form Π u O κ . Note that, if all w κ i = 1, this corresponds to the euclidean norm of R N .

One property of this receding horizon formulation is the guarantee that all vehicles can eventually exit the supervision area, and that no deadlock situation can happen, as stated in the following theorem:

Theorem 2: For all κ, there exists a feasible solution to IH-MPC in which all vehicles in N exit the supervision area in finite time.

Proof: The proof builds on previous results from [START_REF] Gregoire | Priority-based coordination of mobile robots[END_REF], defining the notion of "priority graph" for a set of vehicles. For a valuation of the priority variables π ij , the corresponding priority graph is built by choosing each vehicle as a vertex, and adding the (directed) edge i → j if and only if π ij = 1. It is shown in [START_REF] Gregoire | Priority-based coordination of mobile robots[END_REF] that no deadlock can arise from an acyclic priority graph; this corresponds to the intuitive notion that deadlocks only occur in the presence of circular wait lists.

Assume that some vehicles are in a deadlock situation: by definition, there exists a time step k ≥ κ when all deadlocked vehicles are stopped; moreover, the above results imply that the vehicles have been assigned cyclic priorities. The problem constraints ensure that, for any pair of vehicles (i, j) involved in the circular wait list with assigned priority π ij = 1, we necessarily have s k i ≤ s ⊥ i and s k j ≤ s ⊥ j . Therefore, assigning priority π ji = 1 is also feasible, and thus it is possible to reverse the edge i → j in the priority graph without violating any constraint.

Results on minimum feedback arc sets [START_REF] Younger | Minimum feedback arc sets for a directed graph[END_REF] show that it is possible to remove cycles in a directed graph by reversing a set of its edges, all of these edges participating in at least one cycle. This property ensures that, if all vehicles involved in circular wait lists are stopped at step k, there exists a feasible solution with acyclic priorities which proves the theorem.

IV. AN EQUIVALENT FINITE HORIZON FORMULATION

In the previous section, we introduced an optimal and least restrictive supervisor based on an infinite horizon MPC scheme. However, this formulation is not suitable for practical implementations because of the infinite number of variables. We now present a finite horizon formulation which can be implemented in widely available solvers and, given the proper choice of the time horizon, ensures the safety of the system over an infinite time horizon.

A. Finite Horizon formulation

We propose an approach based on a restriction of the infinite horizon formulation introduced in the previous section. At a time step κ, we consider the restriction of IH-MPC to the first following K steps (corresponding to a horizon T = K∆t), i.e. O κ is minimized under constraints (4), ( 5) and (s1) to (s3) for all steps k such that κ ≤ k ≤ κ + K, and the last considered variables are those at step k + K. We note FH-MPC this finite horizon Model Predictive Control problem.

B. Equivalence of FH-MPC and IH-MPC

By construction, a solution to FH-MPC is guaranteed to avoid collisions up to time (κ + K)∆t. However, due to the dynamic constraints on the vehicles, there is no a priori guarantee that the resulting trajectory can be extended to a safe infinite horizon trajectory, i.e. x κ+1 might not be safe. In what follows, we prove that provided the horizon T = K∆t is chosen long enough, the solution u κ of FH-MPC ensures that x κ+1 is safe.

First, note that the constraints in FH-MPC are a restriction of those of IH-MPC. Therefore, the restriction of any solution to IH-MPC to the interval [κ∆t, (κ+K)∆t] is a feasible (but potentially not optimal) solution of FH-MPC, and so: Proposition 1: If there exists a solution to IH-MPC, then the restriction of the solution to the first K steps is a feasible solution to FH-MPC.

We will now prove a reciprocal implication to Proposition 1, which is stated in Theorem 3. In what follows, we let v max > 0, u a > 0 and u b < 0 be respectively a global upper bound for v i , lower bound for u i and upper bound for

u i such that for all i, v i ≤ v max , u i ≥ u a and u i ≤ u b .
The key idea of the proof lies in the choice of a planning horizon long enough to allow any vehicle to fully stop. The structure of the proof is as follows: Lemma 1 gives a lower bound on the time horizon to allow a single vehicle to stop using discrete dynamics, although with a potential risk of rear-end collisions from following vehicles. In Proposition 2, we show that this time horizon is sufficient to ensure safety from lateral (but not necessarily from rear-end) collisions. Finally, in Proposition 3, we give a slightly higher bound on the time horizon ensuring that all vehicles in a line can all safely stop without rear-end collisions, and therefore that this horizon also ensures absence of rear-end collisions. There exists a discrete control (ũ k i ) κ≤k≤κ+K such that for all κ ≤ k ≤ κ + K, u k i ∈ [u i , u i ] and ũκ i = u κ i , and for which the corresponding dynamically feasible trajectory t → xi (t) = (s i (t), ṽi (t)) verifies si (t 1 + T + ∆t) ≤ s i (t 1 + T ) and ṽi = 0 over [t 1 + T, t 1 + T + ∆t].

Proof: Let (u k i ) be the control corresponding to trajectory s i , and let us define a control (w k i ) as:

w κ i = u κ i , w k i = min(u b , u k i ) for κ < k < κ + K, and
w κ+K i = u b . We construct (ũ k i ) iteratively as ũκ i = u k i and, for k ≥ κ + 1, ũk i = w k i if ṽk i + w k i ∆t ≥ 0 - v k i ∆t otherwise
, where ṽk i is the speed of vehicle i at time k∆t under control (ũ k i ).

As ṽκ+1 i = v κ+1 i ≤ v max ≤ (K -1)∆t|u b | from the hypothesis, there exists a minimal value of k 0 ≥ κ such that ṽk0 i ≤ |u b |∆t; moreover, the condition on K ensures that ṽκ+1 i -(K -2)|u b |∆ t ≤ |u b |∆ t , and so k 0 ≤ (κ + 1) + (K -2) = κ + K -1.
From the definition of (ũ k i ), we know that for all

k 0 + 1 ≤ k ≤ κ + K, ṽk i = 0. Since ũk i ≤ u k i for κ ≤ k ≤ k 0 -1, we also know that sk0 i ≤ s k0 i and ṽk0 i ≤ v k0 i .
Finally, ũk0 i is the minimal admissible control starting from xk0 i ; therefore, sk0+1 i ≤ s k0+1 i = s κ+K i which proves the above lemma.

We can now formulate the following proposition, in the particular case where no pair of vehicles has following constraints:

Proposition 2: Let ∆t > 0 and T = K∆t ≥ vmax |u b | + ∆t. Suppose that there exists a solution of FH-MPC at step κ leading the system of vehicles in N in a state x κ+1 at step κ + 1, and that vehicles have no following constraints. There exists a solution of FH-MPC at step κ + 1 starting from state x κ+1 with horizon T .

Proof: We will show that, using control u κ over [κ∆t, (κ + 1)∆t[, the safe crossing constraints (s1) of FH-MPC remain feasible at step κ + K + 1 without violating the kinodynamic constraints (4) and [START_REF] Qian | Prioritybased coordination of autonomous and legacy vehicles at intersection[END_REF]. First, if s κ+1 i ≥ s ⊥ i for all i ∈ N , the crossing constraints are trivially feasible for κ ≤ k ≤ κ + K + 1. Otherwise, let i be a vehicle such that [START_REF] Dresner | Multiagent traffic management: a reservation-based intersection control mechanism[END_REF] shows that i can fully stop before reaching s ⊥ i and so the crossing constraints involving i remain feasible.

s κ+1 i ≤ s ⊥ i . If s κ+K i ≤ s ⊥ i , proposition
If s ⊥ i ≤ s κ+K i ≥ s ⊥ i , condition (s1)
ensures that each vehicle j conflicting with i either verifies s κ+K j ≤ s ⊥ j (in which case j can fully stop), or s κ+K j ≥ s ⊥ j (in which case constraint (s1) is always satisfied). Finally, if s κ+K i ≥ s ⊥ i , crossing constraints involving i are inactive and thus satisfied. In all cases, the safe crossing constraints remain feasible.

When following constraints are involved, the condition on the time horizon becomes more restrictive. In the following proposition and noting • the ceiling function, we prove a bound ensuring that a line of vehicles can safely stop before the leader reaches its final computed position at the end of the time horizon, without risk of rear-end collisions:

Proposition 3: Let ∆t > 0 and suppose that p vehicles (noted 1, . . . , p from rear to front) are following one another at step κ. We consider a horizon

T = K∆t ≥ vmax |u b | + (p -1) 1 + umax |u b |
∆t + ∆t at time t 1 = κ∆t, and we assume that every vehicle i ∈ {1, . . . , p} has a safe discrete control (u k i ) κ≤k<κ+K such that, for all κ ≤ k < κ + K, u k i ∈ [u i , u i ]. We let t → x i (t) be the safe trajectory over [t 1 , t 1 + T ] for vehicle i under control (u k i ). For all i ∈ {1, . . . , p}, there exists a safe discrete control

(û k i ) κ≤k≤κ+K such that for all κ ≤ k ≤ κ+K, u k i ∈ [u i , u i ], ûκ i = u κ
i and for which the corresponding dynamically feasible and safe trajectory t

→ xi (t) = (ŝ i (t), vi (t)) verifies ŝi (t 1 + T + ∆t) ≤ s i (t 1 + T ) and vi = 0 over [t 1 + T, t 1 + T + ∆t].
Proof: We will prove by induction that, for i ∈ {1, . . . , p}, there exists a dynamically feasible control (û k i )

with ûκ i = u κ i and ûk i ≤ u k i for k ≥ κ such that the corresponding vehicle speed (v k i ) verifies vκ+Ki i ≤ |u b |∆t with K i ∆t = vmax |u b | + (i -1) 1 + umax |u b |
∆t. First, for the rearmost vehicle i = 1, the proof of lemma 1 provides the result with ûk i = ũk i . We now let i ≥ 2 and assume that every vehicle j ∈ {1, . . . , i -1} follows its corresponding control (û k j ). We note ŝk j and vk j the position and speed of vehicle j at step k under this control. Since ûk j ≤ u k j for these vehicles, we deduce from the monotony of the system that the original control solution for vehicle i (u k i ) κ≤k<κ+K prevents rearend collisions if vehicle i -1 applies (u k i-1 ). Therefore, any dynamically feasible extension of (u k i ) is safe over [κ∆t, (κ + K + 1)∆t[. As a result, the set U saf e i (u κ i , [κ, κ + K]) of all admissible controls (û k i ) κ≤k≤κ+K for vehicle i such that ûκ i = u κ i and u k i ≤ ûk i ≤ u k i for κ ≤ k < κ + K is not empty. We note (û k i ) a minimum element of this set (and so ûk i ≤ u k i ); we will prove that vκ+Ki i ≤ |u b |∆t. If for all k ≥ κ, vk i ≤ ṽk i-1 , we conclude that vehicle i stops before vehicle i -1 which proves the result from the induction hypothesis. Otherwise, we let k i 0 ≥ κ be the minimum k such that vk i ≥ ṽk i-1 , and we know that

vk i 0 i ≤ vk i 0 -1 i-1 + u max ∆t.
For all k ≥ k i 0 , we know from the monotony of the system that the control min(û k i-1 , u k i , u b ) prevents rear-end collisions; we deduce that, for k

≥ k i 0 , vk i ≤ vk-1 i-1 + u max ∆t. Therefore, vκ+Ki-1+1 i ≤
vKi-1 i-1 + u max ∆t and we deduce from the induction hypothesis that vκ+Ki-1+1 i ≤ u max ∆t + |u b |∆t. Therefore, we obtain the recursion relation

K i = K i-1 + 1 + umax |u b |
which yields the announced result.

Finally, we conclude that vehicle i can fully stop (without rear-end collisions) at step κ + K i + 1; therefore the set of p vehicles can safely stop before the beginning of step

κ + K if K ≥ K p = vmax |u b |∆t + (p -1) umax |u b | + 1.
Since the recursion ensures that for all i and k, ûk i ≤ u k i , we deduce that s κ+K i ≤ ŝκ+K i which proves the proposition. Finally, we deduce the following theorem: Theorem 3: Let ∆t > 0 and suppose that at most p vehicles are following one another, and consider a horizon

T = K∆t ≥ vmax |u b | +(p-1) 1 + umax |u b |
∆t+∆t. At every time step κ, formulation IH-MPC is equivalent to FH-MPC with horizon T , i.e. the restriction of an optimal solution from one formulation can be extended to an optimal solution for the other in the sense of O κ .

Proof: Consider an optimal solution u ∞ to IH-MPC and an optimal solution u T to FH-MPC. From proposition 1, we know that the restriction of u ∞ is a solution to FH-MPC, and therefore O κ (u ∞ ) ≥ O κ (u T ) since u T is optimal for FH-MPC. Moreover, from propositions 2 and 3, we know that u T can be extended into an infinite horizon solution u T,∞ by recursively solving the problem over successive time steps. Therefore, O κ (u T,∞ ) ≥ O κ (u ∞ ); since O κ only depends on the first time step of the solution, we deduce

O κ (u T,∞ ) = O κ (u T ) and so O κ (u ∞ ) = O κ (u T ).
Note that, from the previous equivalence, the liveness guarantee demonstrated in theorem 2 for the infinite horizon formulation also applies to the finite horizon one.

V. SIMULATION RESULTS

The use of the FH-MPC formulation to design an optimal supervisor has been validated by computer simulation on the example of the intersection of Fig. 1. The simulation is based on the commercial MIQP solver Gurobi [START_REF] Gurobi Optimization | Gurobi optimizer reference manual[END_REF] and the free traffic modeling tool SUMO [START_REF] Krajzewicz | Recent development and applications of SUMO -Simulation of Urban MObility[END_REF], and uses SUMO's path generation algorithm to compute collision segments

u i , u i ] = [-4, +4] m s -2 .
The duration of a time step is chosen as ∆t = 0.25 s, and the horizon T = 4 s, ensuring that the hypotheses of theorems 2 and 3 are verified. All weighting terms w κ i are chosen as w κ i = 1.

A. Microscopic simulation

To illustrate the safety and liveness of the proposed supervisor, we first provide a simple example of 6 vehicles in the intersection presented in Fig. 1. For clarity reasons, vehicles either come from the west or the south. Table I presents the parameters used for the first simulation (in meters and meters per second); note that vehicles 1 and 2 traveling from west to east and vehicles 4 and 5 traveling from south to north do so on the middle lane of the road, so that they do not conflict with vehicle 3, which turns right. For i ∈ {1, 2, 4, 5, 6}, s ⊥ i = 89 m and s ⊥ i = 111 m. Vehicles request a control corresponding to tracking their initial speed, with u κ i,des =1 ∆t (v 0 i -v κ i ). Fig. 2 shows the curvilinear positions of the vehicles under the above supervisor; the gray area corresponds to the [89, 111] m region. Note that vehicle 3, which does not conflict with any other vehicle, crosses at its target speed. The accompanying video submission 1 shows this example, as well as a larger instance of 25 vehicles to illustrate a more complex situation.

B. Computation time

The results from theorem 3 allows to use a relatively small time horizon in the simulator, which enables a fast iterative resolution of the FH-MPC formulation, even with a high number of vehicles. Figure 3 shows the distribution of the computation time needed to solve the FH-MPC formulation for a given number of vehicles, over 10 different instances, for a time step duration of 0.25 s. For up to 16 vehicles, this computation time is lower than the duration of a time step in 95% of the situations, which could allow near real-time applications.

VI. CONCLUSION

In this article, we have proposed a novel approach to design a least restrictive and minimally deviating supervisor for the safe navigation of semi-autonomous vehicles through an intersection. The supervisor is tasked with overriding the desired controls of individual vehicles, only when necessary. Moreover, the overridden control of all vehicles should be as close as possible, for a given metric, to the controls initially requested. We proposed an infinite horizon model predictive control scheme to solve this optimal supervision problem. We then derived a finite horizon MPC formulation and, for a long enough horizon, we theoretically proved that the finite horizon MPC is ensured to maintain safety and avoid deadlocks for all future times. Finally, we implemented the supervisor using the commercial mixed-integer quadratic programming solver Gurobi jointly with the open-source traffic simulation tool SUMO to demonstrate these properties on a small set of 6 vehicles as well as on a larger instance comprising 25 vehicles.

Contrary to previous work on semi-autonomous supervision, our approach is based on the simultaneous resolution of the discrete problem of finding an optimal vehicle ordering, and the continuous problem of optimizing trajectories within this ordering. This method, using the concept of pairwise priorities [START_REF] Gregoire | Priority-based coordination of mobile robots[END_REF] between vehicles, allows to consider more complex intersection geometries while taking into account second-order kinodynamics constraints, with potentially several vehicles behind one another or having merging or diverging paths. Since the supervisor only needs to consider a relatively short time horizon, it is able to run in real-time with at least a dozen vehicles.

This new approach to semi-autonomous intersection supervision opens up several possibilities for future research. Careful study of the hypotheses of the demonstrated theorems and the solutions of the optimization problem can help determine the best instant at which priorities should be assigned to vehicles. Using additional assumptions on vehicles behavior, some of the hypotheses can likely be relaxed to further improve the throughput of the intersection; different objective functions can be used to find solutions balancing respect of individual vehicles wishes with traffic efficiency for the intersection. Moreover, our formulation may be distributed among vehicles, which could provide a fast and efficient method for safe and fully autonomous intersection crossing.

Lemma 1 :

 1 Let ∆t > 0 be a time step duration and consider a horizon T = K∆t ≥ vmax |u b | +∆t at time t 1 = κ∆t. Consider a vehicle i for which there exists a discrete control (u k i ) κ≤k<κ+K such that, for all κ ≤ k < κ+K, u k i ∈ [u i , u i ] corresponding to a dynamically feasible trajectory s i (t) over [κ∆t, (κ + K)∆t].

Fig. 2 .

 2 Fig. 2. Curvilinear positions of the supervised vehicles.

Fig. 3 .

 3 Fig. 3. Distribution of computation times for a given number of vehicles, for a time step duration of 0.25 s.

Also available at https://youtu.be/GF5KKr2V7Sg