Franco Fusco

Guillaume Allibert

Olivier Kermorgant
email: olivier.kermorgant@ls2n.fr

Philippe Martinet
email: philippe.martinet@inria.fr

Benchmarking Nonlinear Model Predictive Control with Input Parameterizations

Model Predictive Control (MPC) while being a very effective control technique can become computationally demanding when a large prediction horizon is selected. To make the problem more tractable, one technique that has been proposed in the literature makes use of control input parameterizations to decrease the numerical complexity of nonlinear MPC problems without necessarily affecting the performances significantly. In this paper, we review the use of parameterizations and propose a simple Sequential Quadratic Programming algorithm for nonlinear MPC. We benchmark the performances of the solver in simulation and compare them with state-of-the-art solvers. Results show that parameterizations allow to attain good performances with (significantly) lower computation times.

I. INTRODUCTION

Model Predictive Control (MPC) is probably one of the most attractive control design methodologies, as illustrated in the state of the art made by Mayne [START_REF] Mayne | Model predictive control: Recent developments and future promise[END_REF]. The ability of MPC to handle nonlinearities, enforce constraints and provide performance/cost trade-offs are just some of the reasons for its popularity. Historically, MPC was developed for systems with slow dynamics and required long computation times. Nonetheless, in the last two decades a lot of work has been done to deploy predictive control also on fast systems. Some of these pose additional restrictions and challenges, as the optimization problem must be solved on embedded hardware with limited computing power. One way to improve the efficiency of MPC schemes is to reduce the complexity of the optimization problem, and several methods have been proposed in order to speed up the solving process like [START_REF] Kouvaritakis | Who needs QP for linear mpc anyway?[END_REF], [START_REF]Efficient robust predictive control[END_REF], [START_REF] Ferreau | An online active set strategy to overcome the limitations of explicit MPC[END_REF], [START_REF] Bemporad | Model predictive control based on linear programming -the explicit solution[END_REF] to cite a few.

Control parameterizations can also be an effective solution [START_REF] Alamir | Stabilization of Nonlinear Systems Using Receding-Horizon Control Schemes: a Parametrized Approach for Fast Systems[END_REF], [START_REF] Van Donkelaar | Model predictive control with generalized input parametrization[END_REF], [START_REF] Muehlebach | Parametrized infinite-horizon model predictive control for linear time-invariant systems with input and state constraints[END_REF], [START_REF] Limaverde Filho | Trajectory tracking for a quadrotor system: A flatness-based nonlinear predictive control approach[END_REF]. The objective is to drastically reduce the number of control variables in the optimization problem while trying to keep the performance loss as small as possible. One of the major advantages of using a parameterization is the ability to decouple the prediction horizon from the control one. This is a very important improvement since for many systems the use of large prediction horizon is mandatory to ensure system stability in closed loop [START_REF] Mayne | Model predictive control: Recent developments and future promise[END_REF]. Even disregarding the issue of system instabilities, it is often necessary to select a large prediction horizon so as to obtain good closed-loop behavior. This is especially true when dealing with small sampling periods.

In this paper, we investigate the use of different input parameterizations in MPC to control non-linear systems. In particular, we review few existing input parameterizations coupled with a simple single-shooting scheme that can solve the MPC optimization using Sequential Quadratic Programming (SQP) quickly and while enforcing control input constraints such as bounds on the maximal command and its rate. Thanks to the parameterization, the sub-problems arising at each SQP iteration are low-dimensional and can thus be solved extremely quickly. Indeed, we show that using the parameterized approach it is possible to attain sub-millisecond performances despite a very dense timesampling and outperform existing state-of-the-art algorithms both on a laptop and on a Raspberry Pi device.

The remainder of this paper is organized as follows. In the next section, we firstly detail the adopted nonlinear MPC formulation and continue by recalling the studied control parameterizations. We propose a simple solver for the parameterized MPC optimization problem in section II-C and benchmark its performances in section III. In particular, we analyze how they change in relation to the chosen number of parameters and compare against state of the art solvers from an existing framework. Finally, the last Section reports our conclusions and proposes future perspectives.

II. NONLINEAR MODEL PREDICTIVE CONTROL

In this section, we recall the formulation of Direct Shooting methods for constrained Nonlinear MPC, which is at the basis of the employed solver. We then discuss how the introduction of a parameterization can help in reducing the complexity of the described problem and recall some existing parameterizations. We finally propose an algorithm to compute optimal parameters, which is a SQP strategy based on the Gauss-Newton Hessian approximation.

A. Direct Shooting Nonlinear MPC

The problem that Model Predictive Control tries to solve can be stated informally as follows: given an evolution model of the controlled system, find an optimal control trajectory which can steer the current state of a system towards a desired configuration, while meeting feasibility requirements such as actuation limits. A variety of mathematical formulations can be found in the literature [START_REF] Binder | Introduction to model based optimization of chemical processes on moving horizons[END_REF], [START_REF] Diehl | Fast direct multiple shooting algorithms for optimal robot control[END_REF] .

Direct methods work by transcribing the problem as a finite-dimensional nonlinear optimization. First of all, the controlled system is modeled as a discrete-time one, characterized by a transition function f that, given the state x k ∈ R n and control u k ∈ R m at the discrete time step k, returns the corresponding state that should be reached by the system at the next time step k + 1, i.e.,

x k+1 = f (x k , u k ) (1)
The optimal control problem can be written as:

min Np k=1 x k -x k 2 Q k + Np-1 k=0 u k -u k 2 R k (2a)
subject to:

x k+1 -f (x k , u k ) = 0 ∀k = 0, • • • , N p -1 (2b) c x (x k ) ≤ 0 ∀k = 1, • • • , N p (2c) c u (u k ) ≤ 0 ∀k = 0, • • • , N p -1 (2d)
wherein N p is known as prediction horizon, while x k and u k represent the desired state and control inputs at the discrete time step k. v A is used to denote the norm of vector v weighted by the positive semidefinite matrix A. x 0 is assumed to be a known constant, corresponding to the actual state of the system when performing the optimization.

Model constraints (2b) can be enforced in two ways. The first one is to use both u and x (respectively the sequence of controls and states) as decision variables and to explicitly deal with these constraints during the optimization. This leads to a sub-family of direct methods known as multipleshooting techniques [START_REF] Bock | A multiple shooting algorithm for direct solution of optimal control problems[END_REF], [START_REF] Albersmeyer | The lifted newton method and its application in optimization[END_REF], [START_REF] Giftthaler | A family of iterative gauss-newton shooting methods for nonlinear optimal control[END_REF]. The second option consists in using a single-shooting method, in which only the control input sequence u is used as free variables for the search. Predicted states can be evaluated iteratively starting from x 0 and the selected control inputs using [START_REF] Mayne | Model predictive control: Recent developments and future promise[END_REF]. In this paper, we have focused our attention on the later.

B. Complexity Reduction via Parameterization

In the optimization scheme presented above, performances highly depend on the prediction horizon. Ideally, one would wish to sample the system at a high frequency to reduce the inaccuracies introduced by discretization. At the same time, predicting over a longer period of time should lead to better results since it enables the algorithm to intelligently consider maneuvers that are globally optimal and steer the system to the desired state. It is thereby beneficial to increase the value of N p in both cases. However, the number of decision variables (and therefore the complexity of the nonlinear optimization) grows with the prediction horizon as well.

One way to simplify the optimization would be to allow only the first N c ≤ N p control samples to be freely allocated in the optimization [START_REF] Allibert | Predictive control for constrained image-based visual servoing[END_REF], [START_REF] Perez-Morales | Multi-sensor-based predictive control for autonomous parking in presence of pedestrians[END_REF], with the remaining N p -N c samples set to the same value of u Nc-1 (fig. 1(a)). More generally, the "free samples" from u could be spread uniformly along the prediction horizon rather than at the beginning, with remaining values obtained by "holding" the free samples constant (fig. 1(b)) or alternatively using linear interpolation (fig. 1(c)). These strategies are known in the literature as Move Blocking. For clarity, they will be referred to as "simple", "ZOH" and "LERP" in the sequel.

A more general option relies on the choice of a set of N b basis functions

φ i : R + → R (i = 1, • • • , N b ). The control
sequence can be generated from these functions by linear combination as:

u k = N b i=1 φ i (k∆t)η i (3)
where the combination coefficients η i are to be determined by the optimization and ∆t is the sampling time used in the prediction model. Some examples in this sense are damped Laguerre polynomials [START_REF] Muehlebach | Parametrized infinite-horizon model predictive control for linear time-invariant systems with input and state constraints[END_REF] and Haar wavelets [START_REF] Lee | Improving computational efficiency of model predictive control algorithm using wavelet transformation[END_REF]. In this paper, we propose to use a set of exponentially damped polynomials in the form:

φ i (t) = et iτ i e -t/τ (4)
All functions from this basis have a unique maximum in R + located at t = iτ , such that φ i (iτ ) = 1. It is thus reasonably easy to tune the shape constant τ , which simply controls the locations of the maxima. To the best of our knowledge, this is the first time such basis is employed in a control parameterization.

The interest in using a basis of functions is that the overall control sequence can "inherit" some desirable properties from its generating functions. As an example, if all φ i are smooth (such as in ( 4)) then the generated sequences naturally feature reduced variations between consecutive samples. In addition, changes in a single combination parameter usually affects a whole portion of the control sequence, rather than a specific one. From a practical point of view, this implies that even with a reduced number of variables it is possible to generate a rather wide variety of control profiles.

All the strategies described so far can be unified under the same concept of control parameterization, which is a linear mapping from a given vector of parameters η ∈ R Nη to the control sequence u,

u =   u 0 • • • u Np-1   =   Π 0 • • • Π Np-1   η = Πη (5)

C. SQP Solver for Parameterized MPC

In this section we detail a Sequential Quadratic Programming solver that can be employed for the solution of the MPC optimization problem (2) when using parameterizations. The proposed algorithm works by formulating a quadratic subproblem, whose solution allows to update the current guess of the parameters η. These two steps (solution of the subproblem and update of the parameters) is repeated until termination conditions are met. As mentioned before, we consider in this paper a single-shooting strategy in which (2b) is directly satisfied, while concerning other constraints we pose limits on the magnitude of the control actions and its variation. Concerning the definition of the quadratic subproblem, we firstly consider a current guess η of the parameters, from which the control sequence u is computed using [START_REF] Bemporad | Model predictive control based on linear programming -the explicit solution[END_REF]. Afterwards, predicted states x are evaluated recursively as function of the initial state x 0 and u using (1). States and controls in the objective are then linearized around η + δ η , with δ η representing variation of the current parameters.

Injecting the linearized states into (2a) leads to

min 1 2 δ T η Hδ η + g T δ η (6a)
i.e., a QP with Gauss-Newton Hessian approximation, with

H = Np k=1 ∂x k ∂η T Q k ∂x k ∂η + Np-1 k=0 Π T k R k Π k (6b) g T = Np k=1 (x k -x k ) T Q k ∂x k ∂η + Np-1 k=0 (u k -u k ) T R k Π k (6c)
The Jacobians of the states with respect to the parameters can be evaluated using the recursive relation

∂x k+1 ∂η = ∂f ∂x ∂x k ∂η + ∂f ∂u Π k (7)
Constraints on the magnitude of the control actions and its variation can be easily translated into constraints on the parameters variation δ η thanks to [START_REF] Bemporad | Model predictive control based on linear programming -the explicit solution[END_REF]. Nevertheless, it is worth noting that insights on the used parameterizations can lead to more efficient solutions. First of all, in simple, ZOH and LERP parameterizations the parameters represent a set of control samples, and it is easy to understand that other values cannot exceed (be less than) the largest (smallest) free samples. As a consequence, control bounds can be readily translated into parameter limits. Similarly, variation constraints need not to be enforced at each discrete sample k, but rather at the free samples only.

When function bases are used to parameterize the control inputs, in the general case there is no way to reduce the number of control constraints. However, if the chosen functions are smooth, it is reasonable to assume that values of consecutive samples will not be too dissimilar. In particular, if a control sample located at a given discrete instant j is within bounds, chances are that neighboring samples will either be within bounds or exceed control limits by a small margin and a similar argument holds for control variations. Therefore, only a subset of the original constraints can be considered at the cost of (small) constraints violations, granting a considerable simplification of the problem.

The sub-problem defined by (6a) and the constraints is a dense QP optimization that can be solved by state-of-the-art solvers, such as qpOASES [START_REF] Ferreau | qpoases: A parametric active-set algorithm for quadratic programming[END_REF], HPIPM [START_REF] Frison | HPIPM: a high-performance quadratic programming framework for model predictive control[END_REF] or OSQP [START_REF] Stellato | OSQP: an operator splitting solver for quadratic programs[END_REF]. Once the solution is computed, a line search algorithm can update the current guess of the parameters as η = η + sδ η with s ∈ [0, 1]. The algorithm can then continue by repeating the previous steps over and over until convergence.

It must be noted that in a classical formulation the number of decision variables and constraints are proportional to the prediction horizon N p . Instead, in the parameterized approach proposed here, they are all proportional to N η . As shown in the next section, the use of a reduced number of parameters allows for considerable speed-ups in the solution process, without degrading the control performance.

III. SIMULATIONS

In this section we benchmark an implementation of the proposed solver in a simulated environment. The parameterization and algorithms detailed in previous sections were implemented as a C++ library, using Eigen [START_REF] Guennebaud | Eigen v3[END_REF] for linear algebra and qpOASES [START_REF] Ferreau | qpoases: A parametric active-set algorithm for quadratic programming[END_REF] to solve the dense QP sub-problems. Different parameterizations with increasing number of parameters have been tested to determine the performances of the solver and a comparison with state of the art methods from the acados framework [START_REF] Verschueren | acados-a modular open-source framework for fast embedded optimal control[END_REF] is included.

The system to be controlled is an inverted pendulum with moving base (sometimes referred to as cart-pole system), which is often used to investigate new approaches in nonlinear control theory. The objective here is not to propose a new control strategy for this type of machine but simply to compare the efficiency of the proposed approach with the state of the art in MPC on a simple system. The state is given by the quadruplet (p, θ, v, ω) (position of the cart, angle of the pendulum, linear and angular velocities). The reader is referred to [START_REF] Verschueren | acados-a modular open-source framework for fast embedded optimal control[END_REF] for the equations adopted in our simulations.

Different sets of simulations are presented in the sequel: first, we focus on a task in which the pendulum has to be kept in balance while also tracking a time-varying trajectory with the base. In this first set of simulations, all parameterizations presented in previous sections are tested multiple times with different values of N η . Afterwards, the results obtained with two parameterizations with N η = 5 are compared to four algorithms from acados. Finally, to test the effectiveness of the parameterized approach in a more challenging scenario, a swing-up task is considered.

In all simulations, the time horizon used for the predictions was set to 1 s, and a very high control frequency was employed to test the algorithms under challenging conditions. In particular, all presented simulations were performed with both ∆t = 10 ms and ∆t = 2 ms, leading to a prediction horizon of N p = 100 and N p = 500 samples respectively. Tests were run on two different platforms: (1) a Laptop computer featuring an Intel Core i7-8750H @2.2 GHz CPU and 16GB of RAM; (2) a Raspberry Pi 400 featuring a Broadcom BCM2711 Cortex-A72 (ARM v8) @1.8 GHz processor and 4GB of RAM. In this way, it was possible to get some insights on the performances of the implemented solver both on a workstation and an embedded device.

A. Performance benchmarking

In this first set of tests, the objective was to gather insight about the performances of parameterizations for an increasing number of parameters. Each simulation lasts 30 s, with the pendulum already starting upward and tracking a time-varying position while keeping balance. Weighting matrices in [START_REF] Alamir | Stabilization of Nonlinear Systems Using Receding-Horizon Control Schemes: a Parametrized Approach for Fast Systems[END_REF] were defined as Q k = diag(200, 50, 7, 2)/N p and R k = 1/N p , and only control bounds were included, with u max = -u min = 0.1 N.

Concerning the tested parameterizations, the simple, ZOH and LERP ones were tested for N η = 2, • • • , 20, with free samples spread uniformly along the prediction horizon in the case of ZOH and LERP parameterizations. In addition, a further parameterization using a basis of exponentially damped polynomials was tested (later referred to as "Poly"), with N b = N η = 2, • • • , 10 parameters. In each of these cases, the parameter τ was set to 1/(N η -2), so that the maxima of the basis functions are regularly spread along the prediction horizon. In addition, control bounds were imposed on 2N η samples distributed along the prediction horizon.

For each parameterization, prediction horizon and number of parameters, the same simulation was repeated 10 times to be able to obtain statistical data on the amount of time required to solve the nonlinear optimization problems. Results are shown in figs. . The results show that the simple parameterization consistently requires less computation time than other strategies, but that it fails to converge to a good solution, with tracking errors that remain more or less the same independently from the number of parameters used. This can be justified by the fact that parameters are able to affect only a small portion of the prediction horizon, with most of the control inputs being handled by a single parameter. ZOH and LERP parameterizations have very similar performances in terms of optimization times, with ZOH usually taking slightly more than ZOH. However, the former features worse tracking performances for an equal number of parameters. Finally, the Poly parameterization takes the longest computation times, which can be explained by the fact that computing the control sequence from parameters requires more flops and that the number of constraints is larger than in other cases. Nonetheless, it generally reaches good tracking performances with less computation times than the ZOH parameterization.

It is interesting to notice that in multiple cases the average optimization times are less than a millisecond. In particular, the results show that almost all the parameterizations could have been employed in real-time experiments in the case of a laptop. The average optimization times on the Raspberry are larger, but still compatible with real-time requirements in the case of N p = 100.

B. Comparison against existing solvers: trajectory tracking

To better place the obtained results within the state of the art, the same simulations performed in the previous section were implemented also using the tools contained in acados. All tuning parameters of the problem, such as the weights in the objective function, were taken as in the previous section to be as fair as possible in the comparison. Termination conditions were also chosen to ensure that the results from the different strategies had the same level of optimality. Four solvers were tested: one based on qpOASES and three HPIPM variants featuring full, partial and no condensing. As before, the simulations were repeated 10 times per solver, with N p set to 100 or 500.

Average time performances (with 95 % confidence intervals) at each iteration of the simulation are plotted in fig. 4 for all acados' algorithms and for LERP and Poly parameterizations 1 (with N η = 5). In addition, for the case N p = 100, we include the results obtained using a parameterization with N η = N p = 100. This corresponds to a MPC problem with full degrees of freedom and is therefore labeled as "full" in the results. In addition, a summary of numerical values is included in the top-half of table I. In particular, for each algorithm we report the average optimization time over all iterations and all simulations in the different conditions. In addition, we evaluated a numerical index that quantifies how faster an algorithm is with respect to a given "baseline". This index is evaluated as:

speed gain = 1 N iters Niters i=1 T bl opt [i] T solver opt [i] (8)
where N iters is the number of iterations per simulation and T solver opt

[i] and T bl opt [i] are the average optimization times required at the i-th iteration to solve the SQP problem using respectively the given solver and the baseline. The two parameterized approaches consistently take shorter computation times than acados' solvers (they are always at the bottom of the plots shown in fig. 4). The Poly and in particular the LERP parameterizations both converge to a good solution faster than every solver from acadosespecially on the Raspberry board, where they run more than 10 times faster on average. Furthermore, even though the values have been omitted for brevity, the average tracking errors obtained with these two parameterizations are slightly smaller than those of acados' algorithms, showing that shorter optimization times have not been obtained by sacrificing optimality. Also note that the full parameterization takes much longer than any algorithm in acados. Since its working principles are close to those of the qpOASES-based solver in acados, this suggests that there is room for improvement in the efficiency of our software implementation.

C. Comparison against existing solvers: swing-up

The last set of simulations considered here involves a swing-up task, with the pendulum starting downward and having to reach the upward position. Weighting matrices were changed to Q k = diag(10 3 , 10 3 , 10 -2 , 10 -2 )/N p and R k = 10 -2 /N p , and the control bounds relaxed to u max = -u min = 12 N. To make the problem more challenging, we also added control variation constraints.

Considering the proposed solver, we decided to find for each parameterization the minimum number of N η that was sufficient to complete the task within 3 s. It turns out that the simple and ZOH parameterizations require at least 16 and 33 parameters respectively. Instead, LERP and Poly parameterizations were successful with as few as 5 and 3 parameters each. Moreover, the full parameterization and the qpOASES-based solver in acados both took too long to converge and were thus not included in the analysis. Similarly to what done in the previous section, we show plots of the average optimization times during the simulations in fig. 5, while numerical results are included in the bottom-half of table I. Once again, the parameterized methods LERP and Poly are confirmed to be faster than existing solvers, by quite a large margin. In particular, it is worth remarking that on the Raspberry Pi board with N p = 100, only the parameterized approaches converged within the sampling period of ∆t = 10 ms consistently.

Fig. 1 .

 1 Fig. 1. Example of control sequences u with a reduced number of degrees of freedom. Even if Np = 10, only 4 variables are freely assigned, while the remaining ones are defined in function of the free samples.

  2 and 3 for experiments run on a laptop and a Raspberry Pi respectively. For each parameterization, two types of graphs are included: average run-times (per iteration) in microseconds vs the number of parameters (figs. 2(a), 2(c), 3(a) and 3(c)) and average run-times vs average position tracking errors (figs. 2(b), 2(d), 3(b) and 3(d))

Fig. 2 .

 2 Fig. 2. Performances (on a Laptop) of the proposed parameterizations with increasing number of parameters. In (a) and (c) the average optimization times are reported for increasing numbers of parameters, while (b) and (d) show the relation between optimization times and tracking errors (each curve corresponds to a parameterization, and each data-point to a specific number of parameters). The prediction horizon for (a) and (b) is Np = 100, while for (c) and (d) it is Np = 500.

Fig. 3 .

 3 Fig. 3. Performances of the proposed parameterizations on a Raspberry Pi board. The meaning of each plot is analogous to those depicted in fig. 2.

Fig. 4 .

 4 Fig. 4. Optimization times for the tracking task. HPIPM (f.c.), HPIPM (p.c), HPIPM: acados solvers based on HPIPM (with full, partial and no condensing); qpOASES: acados solver based on qpOASES; lerp-5, poly-5: proposed solver with Nη = 5. full: proposed solver with Nη = Np.

Fig. 5 .

 5 Fig. 5. Optimization times for the swing-up task.

Simple and ZOH were not considered here since with 5 parameters they did not converge to satisfactory solutions in terms of position error.

IV. CONCLUSIONS

In this work we reviewed the benefits of using input parameterizations in Model Predictive Control to achieve fast and reliable performances. They allow to reduce the number of optimization variables and of the constraints, making the problem solvable even on embedded hardware with very limited computational resources. We implemented a single-shooting SQP algorithm to solve the nonlinear MPC optimization and showed that it can consistently outperform state-of-the-art solvers.

The results obtained so far are very promising, and we believe that it is worth investigating parameterizations more in depth. We would like to focus our attention on different types of basis functions and extend the analysis to nonlinear parameterizations. Furthermore, several works in the field suggest that multiple-shooting methods can provide more stable and fast solutions and ad-hoc condensing strategies for ZOH parameterizations have already been proposed [START_REF] Chen | Efficient move blocking strategy for multiple shooting-based non-linear model predictive control[END_REF]. We are confident that other parameterized approaches would also benefit from the use of lifted techniques. Furthermore, we would like to optimize our software library, add new features and possibly incorporate it into existing frameworks.


N Itare

Kosai Raoof
email: kosai.raoof@univ-lemans.fr

Comparison of different processing for DOA estimation of an Unmanned Aerial Vehicle with few sensors

During the last years, many methods have been developed and studied in the field of source localization. A particular area of interest is the tracking of UAVs (Unmanned Aerial Vehicles) because of the numerous threats that can appear near sensitive sites. Delay and Sum Beamforming (DSB) is one of the methods that can be very useful to face this challenge. Indeed, this technique has a good robustness to noise which makes it an interesting tool. Recently, some processes have been studied to enhance the performance of DSB by taking into account the signature of UAVs. Signals obtained from a microphone antenna can be filtered according to the signature of an UAV before beamforming. Beamforming can also be performed from the measured signals, then the harmonic signature can be considered using the time-frequency representation of the focused signal. A pitch tracking algorithm can provide the fundamental frequency of the signals for consideration of the UAV's signature. Another interesting approach is the Steered Response Power (SRP) which can perform well in noisy environment. The use of generalized cross correlation with different spectral weightings provides a wide range of options. This study aims at comparing the performance of beamforming with time-frequency representation and SRP-PHAT on an experimental measurement with a UAV in flight.

Introduction

Airspace traffic has to be regulated because of the increase of users, particularly unmanned aerial vehicle (UAV) flyers. The amount of recreational flyers is nowadays more important and UAVs can be extremely useful in other activities like deliveries or aerial imagery. These flies can be dangerous near airplane traffic or airports. Privacy is also an issue near sensitive sites or more generally when flyers use UAV for surveillance. Several surveys have been conducted about the threats and the different methods that can be used to counter them [START_REF] Humphreys | Statement on the security threat posed by unmanned aerial systems and possible countermeasures[END_REF], [START_REF] Altawy | Security, Privacy, and Safety Aspects of Civilian Drones: A Survey[END_REF]. Among these methods, the use of acoustical methods enables to estimate the direction of arrival (DOA) of an UAV thanks to the sound emitted during the flight. Beamforming is often used because of its robustness to noise and its real time implementation. Considering the UAV's acoustic signature enables to enhance the performance of beamforming. This can be made using the time-frequency representation of the temporal beamformer's output [START_REF] Itare | Tracking of an Unmanned Aerial Vehicle with few sensors using time-frequency representation[END_REF]. Thus, energy is calculated with time-frequency bins associated with the harmonic signature of the drone. Other algorithms based on the time difference of arrival (TDOA) can estimate the DOA of the drone [START_REF] Brutti | Comparison Between Different Sound Source Localization Techniques Based on a Real Data Collection[END_REF]- [START_REF] Dibiase | Robust Localization in Reverberant rooms, from Microphone arrays: signal processing, Techniques and Applications[END_REF]. These algorithms exploit the generalized cross correlation (GCC) between pairs of microphone with specific weightings (SCOT [START_REF] Carter | The smoothed coherence transform[END_REF], PHAT- [START_REF] Donohue | Audio signal delay estimation using partial whitening[END_REF], …) . GCC with a phase transform also called GCC-PHAT is known to be robust to reverberation. The Steered Response Power with Phase Transformation (SRP-PHAT) makes use of GCC-PHAT to estimate the source DOA estimation. The aim of this study is to evaluate the performance of beamforming with timefrequency representation and SRP-PHAT on an experimental measure.

The work is presented as follows: Part 2 describes the two DOA estimation approaches used, Part 3 tackles the comparison between the two algorithms on the experimental case, and Part 4 concludes and gives some perspectives.

2.

DOA estimation algorithms

Steered Response Power -Phase Transformation

Given a microphone antenna with N microphones, the signal 𝑝 𝑛 (𝑡) received by the n th microphone can be described with Equation [START_REF] Humphreys | Statement on the security threat posed by unmanned aerial systems and possible countermeasures[END_REF]. In this equation, 𝑠(𝑡) is the signal emitted by the source, 𝛺 𝑠 is the direction of the source, ℎ 𝑛 (𝛺 𝑠 , 𝑡) is the impulse response between the source and the n th microphone, and 𝑏 𝑛 (𝑡) is the noise due to the environment or the microphone system.

𝑝 𝑛 (𝑡) = 𝑠(𝑡) * ℎ 𝑛 (𝛺 𝑠 , 𝑡) + 𝑏 𝑛 (𝑡).

The generalized cross correlation between microphones 𝑛 and 𝑚 is given by:

𝑅 𝑛𝑚 () = ∫ Φ 𝑛𝑚 (f) +∞ -∞ 𝐺 𝑛𝑚 (𝑓)𝑒 2𝜋𝑗𝑓𝜏 𝑑𝑓, (2)
with Φ 𝑛𝑚 (f) a weighting of the cross-spectrum 𝐺 𝑛𝑚 (𝑓) between microphones 𝑛 and 𝑚 . The weighting corresponding to the phase transformation is given by:

Φ 𝑛𝑚 (f) = 1 |𝐺 𝑛𝑚 (𝑓)|
. SRP-PHAT is the extension of GCC-PHAT with multiple microphones, it is computed with:

𝑃(Ω) = ∑ ∑ 𝑅 𝑛𝑚 (𝜏 𝑛𝑚 (Ω)) 𝑁 𝑚=1 𝑁 𝑛=1 , (3)
where 𝜏 𝑛𝑚 (Ω) is the TDOA between microphones 𝑛 and 𝑚. The TDOA can be calculated using a plane or spherical wave model given the direction Ω. The direction Ω which maximizes 𝑃(Ω) gives the estimation of the source DOA.

Two beamforming-based processing

Temporal Delay and Sum Beamforming

Delay and Sum Beamforming (DSB) enables to focus a signal in a direction Ω by calculating the TDOA between the antenna's microphones and a reference microphone. The delays can also be calculated with a plane or spherical wave model. Each signal is delayed with its own TDOA and then all the delayed signals are summed. The DOA estimate is given by the direction Ω which maximizes the energy of the focused signal.

Beamforming adapted to the acoustic signature

The time-frequency representation (TFR) of the focused signal in the direction Ω enables to consider the acoustic signature of the source. This TFR can be obtained with the Short Time Fourier Transform (STFT). An interesting property of UAVs' signatures is their harmonic characteristics. By using a fundamental frequency detection algorithm, it is possible to automatically select time-frequency bins from the TFR corresponding to chosen harmonics linked to the spectral properties of the signals produced by the UAV. The DOA estimation is therefore given by the direction Ω that maximizes the energy of these selected bins instead of considering all the spectral content. The algorithm used here for the fundamental frequency detection is the Spectral Harmonic Correlation (SHC) [START_REF] Zahorian | A spectral/temporal method for robust fundamental frequency tracking[END_REF]. The choice of this frequency is very important to obtain a good estimate of the DOA. A procedure has been implemented for selecting the fundamental frequency among those given by SHC, which is detailed in part 3.2.

Comparison on an experimental trajectory 3.1 Array geometry and measurement set-up

The measurements were realized with an array of 10 BSWA Technology MPA 416 1/4 in. microphones (20 Hz-20 kHz) in the arrangement shown in Figure 1. Three different intermicrophonic spacings were used to give the [220,5 3430] Hz array bandwidth:

‖𝑥 1 ‖ = ‖𝑥 4 ‖ = ‖𝑥 7 ‖ = 0.05 𝑚, ‖𝑥 2 ‖ = ‖𝑥 5 ‖ = ‖𝑥 8 ‖ = 0.2 𝑚, ‖𝑥 3 ‖ = ‖𝑥 6 ‖ = ‖𝑥 9 ‖ = 1.1 𝑚. (1)
A PXI-1036 chassis from National instruments is used with a laptop for the recording of the acoustic signals with a sampling frequency of 20 kHz. A source is located from its direction Ω = (𝜑, 𝜃), which is respectively the azimuth and the elevation, and with its distance 𝑟 to the reference microphone (placed at the origin). During the recording, the drone starts with a circular trajectory followed by a smaller one with both a height to the ground around 14 m. The drone used is a Phantom IV from DJI.

TFR parameters

According to the study of the acoustic signature of the drone used [START_REF] Blanchard | Caractérisation de drones en vue de leur localisation et de leur suivi à partir d'une antenne de microphones[END_REF], [START_REF] Blanchard | Acoustic localization and tracking of a multi-rotor unmanned aerial vehicle using an array with few microphones[END_REF], two types of harmonics are present. Weak harmonics are produced by the rotation of the rotors and strong harmonics are produced by the aerodynamic phenomenon. The rotor frequency 𝑓 0 is related to the blade passing frequency 𝑓 𝑏𝑝 with the number of blades 𝑁 𝑏 , 𝑓 𝑏𝑝 = 𝑁 𝑏 × 𝑓 0 . Since there are two blades, weak harmonics (of 𝑓 0 ) are odd and strong harmonics are even. The spectrogram of the reference microphone is presented in Figure 2. Because of the noise present during the recording, weak harmonics are not visible in the spectrogram and 𝑓 𝑏𝑝 is detected by SHC instead of 𝑓 0 (𝑓 𝑏𝑝 detections are visible with red points in Figure 2). In order to select the time-frequency bins , where 𝑖 = 1, … , 𝑁 ℎ , 𝑁 ℎ is the number of harmonics chosen, and 𝑄, the quality factor. The signal is cut into portions and for each portion, 𝑓 𝑏𝑝 is detected by taking the maximum of the SHC calculated in a given frequency range. To verify that this 𝑓 𝑏𝑝 gives a good DOA estimation, a validation procedure is added to the TFR. If the DOA estimate is higher or lower of 10° than the previous one in azimuth or in elevation, then another 𝑓 𝑏𝑝 is chosen and another DOA is estimated. This condition enables to select the appropriate frequency content by favouring the continuity of the drone trajectory. The frequency range for the SHC calculation is extended and other DOAs are estimated if no previous DOA estimate meets the conditions.

Results

The DOA estimation is computed with classical DSB, DSB with TFR and SRP-PHAT on signal portions of 3000 points. DSB is performed with 2048 points both for classical and TFR, using a plane wave model with a resolution of 4° in azimuth and 2° in elevation for the spherical search grid with a radius of 1 m. SHC is computed with 8192 points and 5 harmonics. For the energy calculation with the TFR, 5 harmonics are chosen and Q=8. SRP-PHAT is performed with 2048 points for the cross-spectrum calculation with the same search grid as for DSB. Figure 3 shows the evolution of azimuth and elevation during time in comparison with data provided by a GPS embedded in the drone. The drone has an incertitude around 3 meters which can explain the constant bias between the GPS data and the results. It is visible that selecting frequencies from the TFR enables to enhance the performance of beamforming both in azimuth and elevation. SRP-PHAT gives also better results than classical DSB but some fluctuations are still present. Table 1 presents the mean errors and standard deviations for the three approaches.

Beamforming with TFR and SRP-PHAT give results very close in azimuth but the TFR approach performs better in elevation. The second portion gives a good DOA estimate for SRP-PHAT (d) but not for classical DSB (b).

Figure 6 shows energy maps in the azimuth/elevation plane as well as SHC calculated between 100 Hz and 250 Hz for both signal portions using the TFR approach. This makes it possible to see which blade passing frequency is chosen for the computation of the energy contributed by its harmonics. The first portion [Fig. 6 (a)] gives a clearer map than classical DSB [fig. 5 (a)] with a good estimate of the DOA (a) using the max of SHC for 𝑓 𝑏𝑝 (d). The second portion gives a poor estimate of the DOA (b) using the max of SHC for 𝑓 𝑏𝑝 (e). Taking another 𝑓 𝑏𝑝 in the SHC (f) enables to change the content selected for the energy calculation and thus gives a good estimate of the DOA (c). In this case, the third maximum (156 Hz) provided by SHC (e) is selected given the cartography (c). The TFR is very interesting in its ability to choose different DOA estimates depending on the blade passing frequency chosen. A good example is between 28 s and 32 s in the trajectory where a low frequency content is present associated to a car acceleration (see Figure 2). During this time, classical DSB is unable to provide a good localization (Figure 3) while the TFR approach performs well choosing the appropriate frequency content in the energy calculation. SRP-PHAT is also more effective than classical DSB but with more fluctuations than the TFR approach.

Conclusions

This study compares three approaches for estimating the DOA of an UAV. The first uses beamforming, the second, beamforming with a time-frequency representation and the last concerns the steered response power with a phase transformation. The performance of these approaches is evaluated on an experimental trajectory where the drone is flying in circles. Results

show that estimations with the TFR and SRP-PHAT are better than classical DSB. However, TFR approach gives better estimations thanks to frequency content selection. The detection of the blade passing frequency thanks to the SHC is an important step in the process and a bad choice can result in a poor estimate of the DOA. A procedure to avoid this case is presented and enables to better follow the trajectory of the drone. In the presence of strong noise or perturbing sources, classical DSB is not able to estimate the DOAs but SRP-PHAT and the TFR approach still give correct DOAs. To enhance the performance of SRP-PHAT, it could be interesting to compute different weightings as SCOT, PHAT-, or others.

Figure 1 :

 1 Figure 1: Microphone antenna used for the experimental measurements.

Figure 2 :

 2 Figure 2: Spectrogram of the reference microphone, red points are detections of 𝑓 𝑏𝑝 using SHC.

Figure 3 :

 3 Figure 3: Evolution of azimuth and elevation in time for classical DSB, DSB with TFR, and SRP-PHAT. The GPS trajectory is shown in black dotted lines.

Figure 4

 4 Figure 4 shows an example of an energy map in the azimuth/elevation plane of a portion of signals where classical DSB is very noisy and gives a poor estimate of the DOA in comparison with the TFR approach and SRP-PHAT. Because all the frequency content is selected with classical DSB, the signal to noise ratio (SNR) is very low. The TFR enables to enhance the SNR and therefore clarifies the energy map to find the DOA associated to the source. SRP-PHAT gives also a close estimation thanks to the GCC and weighting.

Figure 4 :

 4 Figure 4: Energy maps in the azimuth/elevation plane for the same signal portion with classical DSB, DSB with TFR, SRP-PHAT (from left to right). The maximum of each cartography gives the DOA estimate for the considered method. For comparison the three DOA estimates are shown on each map Two successive portions of signals have been chosen to demonstrate the importance of choosing the right 𝑓 𝑏𝑝 . Figure 5 presents the energy maps associated to both portions for classical DSB and SRP-PHAT. The first portion shows accurate DOA estimates for the three methods. Indeed, despite a map with a lot of energy everywhere, classical DSB gives a good estimate of the DOA (a). SRP-PHAT gives a localized energy maximum near the true DOA for both portions [(c) (d)].

Figure 5 :

 5 Figure 5: Energy maps in the azimuth/elevation plane for: the first chosen signal portion with classical DSB (a) and SRP-PHAT (c); the next signal portion with classical DSB (b) and SRP-PHAT (d).

Figure 6 :

 6 Figure 6: First signal portion chosen: energy map in the azimuth/elevation plane for the TFR approach (a) and SHC associated (d). Second signal portion: energy map in the azimuth/elevation plane for the TFR approach and for the first fbp estimate (b) with the SHC associated (e), for another (right) fbp estimate (c) with the SHC associated (f).

Table 1 :

 1 Comparison of mean errors and standard deviations for classical DSB, DSB with TFR, and SRP-PHAT


Guénolé Huitric

Michael Rodriguez-Fano

Lucas Gournay

Nicolas Godin

Marius Hervé

Gaël Privault

Julien Tranchant

Zohra Khaldi

Marco Cammarata

Eric Collet

Etienne Janod

Christophe Odin
email: christophe.odin@univ-rennes1.fr

Impact of the terahertz and optical pump penetration depths on generated strain waves temporal

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

INTRODUCTION

The study of Photo-Induced Phase Transitions (PIPT) is in constant development, owing to new experimental techniques and progress in the physical and theoretical interpretation 1 . Moreover, the recent possibility to generate almost single cycle terahertz pulses of high intensity offers appealing perspectives for the excitation and investigation of coherent lattice vibrations [2][3][4][5][6] .

It was soon recognized that, in addition to the microscopic excitation mechanisms, the macroscopic excitation profile is also strikingly important both for the generation of new phases, their propagation in space and time, and for the interpretation of collected data. The finite penetration depth and hence excitation profile of the pump will induce an inhomogeneous distribution of phase nuclei from the surface 7 , that may propagate in the depth of the sample (for instance, photoinduced metallic domains in an insulating phase 8 ). The initial spatial pattern is determined by the excitation profile and mechanism, and subsequent evolutions (nucleation, phase front propagation) are driven by corresponding spatial gradients. This spatial and dynamical mixture of sample properties modulates the sample dielectric function, whose effects are probed through time-resolved optical reflectivity or transmission measurements. A model of the full optical detection process was developed in the seminal work by Thomsen et al. 9 , pointing out the role and interplay between the pump and probe in the interpretation of the transient signal.

It is admitted that optical pump microscopic excitation processes are mainly mediated by electrons, while terahertz pulses can also interact directly with the lattice through infrared active phonon modes. Among possible transformations that are triggered by such ultrafast strong pulses are coherent phonons or strain-waves 10,11 . In addition to possible applications to probe material properties, such strain wave were shown to be able to drive a semi-conductor to metal phase transition 12 .

For these reasons, the comparison of ultrafast pump properties is of primary importance in the understanding of PIPT phenomena, in particular in metal-insulator phase transitions of strongly correlated metals.

Vanadium sesquioxide V 2 O 3 is the working horse in the study of Mott-Hubbard phase transitions 13,14 . At ambient pressure, it presents a first order symmetry-breaking Paramagnetic Metal (PM) to Antiferromagnetic Insulator (AFI) phase transition near 160 K, while a third phase (Paramagnetic Insulator, PI) can be obtained by changing pressure or chemical doping (Fig. 1a). At room temperature and pressure, V 2 O 3 is in a paramagnetic correlated metallic phase, with a corundum structure with two V 2 O 3 formula units in the rhombohedral R 3c unit cell (Fig. 1a).

The ultrafast response of such materials in pump-probe experiments are very important both for the understanding of the underlying physics, and for potential applications. For instance, optical pump experiments have shown the possibility of inducing metallic nucleation-growth from the low-temperature insulating phase 15 , while very high terahertz-pump fields (> 1 MV/cm) at 4 K achieved insulating to metal transition through electronic tunneling 16 .

In the following, we compare the optical relative transmission transients induced by a terahertz or an optical pump on a V 2 O 3 thin film, using the same optical probe, to get information on photoinduced static and dynamical strains within the acoustic timescale.

Materials and Methods

2

ACCEPTED MANUSCRIPT / CLEAN COPY

Sample preparation

A thin film of V 2 O 3 was deposited on a c-cut (0001) sapphire substrate by sputtering of a vanadium target in Ar discharge. Post-deposition annealing at 530 • C in a reducing atmosphere yields a 210±10 nm thick crystallized and single-phased V 2 O 3 layer.

The sample thickness of 210 nm was chosen on the one hand because it is higher than the penetration depth of the optical pump (80 nm), and lower than that of THz pump (1300 nm) used in pump-probe experiments presented in this work. On the other hand, a thickness of 210 nm results in an Optical Density OD ≈ 1 (probe penetration depth is 95 nm), which usually provides a good contrast in ultrafast transmission measurements.

Experimental set-up and simulations.

Schematics of the experimental setups are given in Fig. 1c where THz pump/optical probe and optical pump/optical probe techniques are used. The output of a femtosecond regenerative amplifier (Coherent) that delivers 100 fs pulses at 1 kHz repetition rate is split and used for both pump and probe beams. For the THz experiments, the main part (3 mJ) of the amplifier output is used to generate intense THz fields by optical rectification using the tilted-pulse-front technique 2,17 in LiNbO 3 . We obtain a fluence of 1 µJ/mm 2 and an electric fields reaching 250 kV/cm (measured by Electro-Optical Sampling, EOS). The remaining output of the amplifier seeds an Optical Parametric Amplifier, OPA (TOPAS Light Conversion) in order to obtain the 740 nm pulse used as probe. The photo-induced signal is recorded with a lockin detection. For the optical pump experiment, we use the output of an 3 ACCEPTED MANUSCRIPT / CLEAN COPY OPA as pump. The wavelength is set to 540 nm and its fluence fixed to 35 µJ/mm 2 . Due to experimental constraints, in this setup, a 700 nm probe wavelength is preferred. Both experiments are performed in transmission configuration at nearly normal incidence.

In this way, we measured the temporal relative variation of the optical transmission T r , ∆T r T r = [T r (on) -T r (o)] /T r (o), where "on" and "off" mean respectively with and without pump. Eventually, a timeindependent constant corresponding to the mean signal at negative times was subtracted. No other data treatments have been performed, nor any part of the signal removed.

Data reduction and all simulations were performed with home-made Python codes, that uses Numpy 18 The V 2 O 3 sample was characterized by terahertz Time Domain Spectroscopy. The time traces of the THz electric field, measured by Electro Optical Sampling with a 1 mm thick ZnTe crystal, are shown in Fig. 2a. From their Fourier transforms, we extracted the complex dielectric function of the film using Tinkham Formula [21][22][23] . The mean absolute value of the transmission is T ≈ 0.375, from which a static conductivity of σ o ≈ 1000 S/cm was calculated from Tinkham formula 21 . This estimation is close to those obtained with the same type of samples by THz-TDS 15,[24][25][26] . The frequency dependence of the real and imaginary parts of the dielectric function are shown in Fig. 2b, with the best fit using Drude's model for the conductivity Eq.1

ε ω = ε 1 + i ε 2 = ε ∞ + i σ ω ωε o where σ ω = σ o 1 + iωτ ( 1
)
where ω is the angular frequency (rad/s), ε o is the vacuum dielectric permittivity and ε ∞ the dielectric function at infinite frequency. We ob-

4

ACCEPTED MANUSCRIPT / CLEAN COPY tained an electronic scattering rate Γ = 1/τ ≈ 14.3 ± 0.2 THz, and a static conductivity of σ o ≈ 1000 S/cm. These values are consistent with those obtained in Refs. 26,27 . These parameters will be used in the following to calculate the sample refractive index at terahertz frequencies.

2 Results and Analysis.

2.1 Modulation of optical transmission by terahertz or optical pump.

In the following we report the comparison of temporal relative transmission variation using either a THz pump or an optical pump. The THz pump EOS temporal trace, shown in Fig. 3c, has a maximum electric field ≈ 250 kV/cm and typical full duration ≈ 4 ps. The energy power spectrum has its maximum near 0.5 THz (Fig. 3d) with a Full Width at Half Maximum FWHM ≈ 0.5 THz. The optical pump at 540 nm has a duration of 100 fs. For experimental reasons, we could not use exactly the same probe wavelengths, 740 and 700 nm respectively for THz and optical pumps. However, it can be shown that the sensitivity function (see below Fig. 4) that determines how the probe interacts with the permittivity modulation induced by the pump are almost equal over the full film thickness, since the ratio of the complex index to wavelength ñ/λ are very close (difference of less than 6% for the real part, and 1% for the imaginary part).

The temporal evolution of the raw relative transmission variation

∆T r (t) T r
is presented on Fig. 3 for both pumps, with probes at almost the same wavelength.

Two features are immediately apparent. First, both signals increase within about 30 ps, reaching a maximum before decaying to an almost constant plateau value for delay times larger than 90 ps. In fact, for larger times, small oscillations modulate the plateau. For instance a bump appears near 150 ps in the terahertz pump profile. It is less clear for the optical pump profile. Second, the magnitude of the relative transmission variation with the optical pump is larger than the value with the terahertz pump, as expected from the largest available energy and fluence of the optical pulse.

Most significant is the difference in the profiles within the 35-90 ps time range. While the decay is rather uniform for the terahertz profile, a dip near 44 ps, and second maximum at 56 ps are clearly seen for the profile with the optical pump, making this profile strikingly different from the terahertz profile. Moreover, an additional strong sharp peak signal appears at very short time delays.

By noting that an acoustic timescale is the travel time of a wave from the top surface to the interface of the film with the substrate, that is in our case τ a = d/v ≈ 29 ps, and an electronic time-scale is typically a few 0.1 ps, we can split our analysis into two parts.

First, we analyse in more details the optical short-time trace, that was studied with a shorter time-step of 0.01 ps, as presented in Fig. 3e. This photo-induced signal can be divided in three contributions. A first ultrafast response is located just after the pump arrival and lasts for less than 0.5 ps (See Fig. 3e). This feature is related to the photon-electron interaction. At larger time delays, oscillations on a slowly growing background are perceptible. The extraction of the oscillating component is done by removing the two other contributions. The peak is fitted by the convolution of a Gaussian of full width at half height 120 fs, with an exponential decay of 80 fs time constant, while the increasing background is reproduced thanks to a first order polynomial. The Fourier Transform of the remaining oscillating component unveils a 7.2 THz mode (see Fig. 3g), in accordance with the frequency of the A 1 g optical phonon mode previously measured through Raman 28,29 and pump-probe experiments [29][30][31] . Concerning the THz pump experiment, neither the peak nor any oscillation were detected. The lack of oscillating component is due to the long THz excitation (∼1 ps) compared to the period of the phonon mode (∼0.14 ps), which washes out the oscillations due to the dephasing.

Second, as discussed before, for times larger than 1 ps, the features of the data fall in the acoustic timescale, the positions of peaks and dips corresponding more or less to multiples of the acoustic travel time n × τ a with n an integer. More precisely, we will show that the measured relative transmission variations as a function of time comes from the modulation of the complex dielectric function by photo-induced elastic strain-waves. For metals, the strain waves are usually generated by a thermoelastic mechanism [9][10][11]32 . In this process, part of the pump pulse energy is absorbed by the electrons, which thermalize in a few hundreds of femtoseconds. During and after the pump pulse, the electronic and lattice degrees of freedom exchange energy through electron-phonon coupling and equilibrate to a common temperature, typically within less than a few picoseconds, a process that is usually described with a two-temperature model 33 . This almost instantaneous heating of the lattice, at the elastic timescale, induces an inhomogeneous thermal stress related to the pump energy deposition profile. The related strain modulates the dielectric function as a function of depth and time, giving rise to a time-dependent optical transmission variation.

In the following, we use the simplest one-dimensional model to interpret our data with the minimum number of fitting parameters.

Thermoelastic Model

The relative transmission variation ∆T r (t) T r depends both on the pump excitation profile that determines the strain pulse shape η(x,t), and on the probe pulse transmission through the perturbed sample. The pioneering work of Thomsen et al. 9 with a semi-infinite sample was generalized to a multilayer system by Matsuda et al. 34 .

In the following, we consider a layer of thickness d on top of a semiinfinite substrate (Fig. 4). Note however that a more complete model, with a finite-width substrate, gives identical results.

Assuming that the permittivity, or refractive index, is modulated by the strain, and using the approximation ∆ ñ(x,t) = ∂ ñ ∂ η η(x,t), the relative change of transmission can be written as:

∆T r (t) T r = Re ∂ ñ ∂ η d 0 η(x,t) f (x, λ pr , ñ, ñs , d)dx (2)
where Re stands for the real part of a complex number, and f (x, λ pr , ñ, ñs , d) is a complex-valued "sensitivity function" 9 . This sensitivity function describes the interaction of the probe with the permittivity modulation, and accounts to both absorption and interferences effects. The parameters are the probe wavelength λ pr , the film and substrate complex refractive indexes ñ and ñs , and the film thickness d. All these values are known either from bibliography or direct measurement (see Table1).

In this equation, we made the assumption that the derivative of the film photoelastic complex constant ∂ ñ ∂ η is not a function of depth. It is an unknown that can only be estimated by fitting the data.

The strain profile η(x,t) = ∂ u(x,t) ∂ x , were u(x,t) is the displacement, can be obtained by solving the two coupled elastic-wave equations:

∂ 2 u ∂t 2 = v 2 ∂ 2 u ∂ x 2 -v 2 ∂ η th ∂ x for 0 ≤ x ≤ d (3) ∂ 2 u s ∂t 2 = v 2 s ∂ 2 u s ∂ x 2 for x>d (4)
The subscript s indicates the substrate. The last term in the right hand side of the first equation is the thermal strain, defined as:

η th (x,t) = η o exp(-x/ξ pump )Θ(t), η o = β T ∆T o 1 + ν p 1 -ν p (5)
Θ(t) is the Heaviside function, ξ pump the pump penetration depth, β T the linear expansion coefficient, ∆T o the initial lattice maximum temperature rise at the surface, and ν p the film Poisson ratio. The validity domain of such a system of equations is discussed for instance in Refs. 9,10 . The initial temperature jump is related to pump fluence F p , intensity reflection coefficient R p and sample heat capacity per unit volume C, by the formula

∆T o = (1-R p )F p
Cξ pump . It is worth mentioning however that the details of the energy exchange mechanisms that produce the initial stress on the lattice are not important in this model, as long as the initial photoinduced lattice perturbation is very fast compared to the acoustic timescale, leading to an approximate stepwise temporal rise. Thus, four dimensionless parameters determine the solutions. First, the strain wave reflection coefficient R defined as

R = Z s -Z Z s + Z (6)
where Z s and Z are respectively the substrate and thin film acoustic impedances (Z = ρv with ρ the density and v the sound velocity). Second the ratio K = d ξ pump of the film thickness d and the pump penetration depth ξ pump . if K 1, the excitation profile is almost constant within the film, whereas the exponential decay is strongly marked when K 1. Third, V = v v s the ratio of sound velocities in the layer and the substrate. And finally the strain at the top-layer surface η o .

Using Laplace transforms methods, the solutions of Eq.4, u(x,t) and η(x,t) in the film and the substrate, were calculated analytically, with the initial conditions of no displacement and strain, as well as zero velocity at t=0, boundary conditions of continuous displacement and stress at the film-substrate interface, and free surface at the air-layer interface (zero-stress condition).

Thus, by combining the solutions of Eq.4, and the formula Eq.2, it is possible to simulate and fit the data, as explained below.

We show in Fig. 4 the normalized excitation profiles for the two pumps and sensitivity functions at the two probe wavelengths, using values from Table1. As expected, the terahertz excitation is almost constant over the whole film thickness (ξ T Hz =1300 nm), in contrast to the optical excitation that has a penetration depth of ξ 540nm =80 nm. As already said above, the probe sensitivity function f (x, λ pr , ñ, ñs , d) are very close to each other, and the difference has a negligible effect on the transmission. Thus, we can conclude that if a change of transmission is observed, it is a consequence of the different pump excitations. Interestingly, we also note that the probe sensitivity function is almost constant, despite the ξ probe =95 nm penetration depth of the 700 nm or to constructive interferences within the film.

To get further insight into the interpretation of the data ∆T r (t) T r , we first discuss a rough approximation of Eq.2. Assuming a constant complexvalued sensitivity function f (x, ...) = f (...), Eq.2 can be easily simplified to:

∆T r (t) T r ≈ Re ∂ ñ ∂ η f (...) d 0 η(x,t)dx = Re ∂ ñ ∂ η f (...) d 0 ∂ u(x,t) ∂ x dx Thus ∆T r (t) T r ≈ Re ∂ ñ ∂ η f (...) [u(d,t) -u(0,t)] = F ∆u(t) (7)
where ∆u(t) = u(d,t)u(0,t) and F = Re ∂ ñ ∂ η f (...) . Within this approximation, the relative transmission variation is simply proportional to the variation of the film thickness ∆u(t).

The behaviour of ∆u(t) for both terahertz and optical pumps, as well as snapshots of the strain profile η(x,t) at typical times, is illustrated in Fig. 5, using values from Table1, and three values of the strain reflection coefficient R = 0 (matched layer), R = ±0.2. The behaviour of the strain spatial profile is interesting (Fig. 5c andd), and can be understood using three contributions. A time-independent contribution corresponding to the initial temperature profile, and a propagating and counter-propagating strain waves launched at the free surface and film-substrate surface by the initial stress impulse. The resulting strain is the sum of these three contributions, and evolves depending on the reflections at the free surface or film-substrate interface. At each reflection on the film-substrate interface, the amplitude of the strain-wave is multiplied by R, resulting in a magnitude after n round trips of duration 2τ a = 2d/v reduced by a factor R n . For times 2nτ a long enough such that R n 1, only the time-independent strain (static strain) survives. Thus, we expect the time-independent strain to determine the magnitude of the plateau. Of course, other phenomena are susceptible to modulate the shape of the transmission signal. For instance, the static strain being related to a temperature profile, and thus to a temperature gradient, heat diffusion, loss and transfer, will smooth and washed away the temperature gradient, and strain, till thermodynamical equilibrium is reached. Such phenomena are negligible on the studied timescale, and are not taken into account in this model.

As remarked by Schick et al. 41 who treated the case of a matched layer R = 0, the full strain presents a positive and negative contributions in the limit d ξ pump . In the opposite limit d ξ pump , the strain has almost the same sign over all the film thickness. It was explained by the fact that the static strain always compensate the strain wave. This conclusion remains valid in our model, at least for the small considered reflection coefficient magnitude |R| < 0.2.

Clearly, a comparison of the simulations shown in Fig. 5a,b with the raw data of Fig. 3 shows that the studied V 2 O 3 film is not matched to the sapphire substrate. If the layer were matched, ∆u(t) would reach Markers are placed at times at which corresponding strain wave snapshots are plotted as a function of position. Panel c and d: strain wave proles η(x,t) in the thin lm (left of the vertical black line) or in the substrate (right of the vertical black line) at a given time, as a function of position. Horizontal dotted gray lines represent zero strain η = 0 line for each time. For the simulation, the sound velocity in the substrate was set to 1.5v. a plateau after the strain wave launch at the film-substrate interface has travelled back and forth in the film, that is 2τ a = 2d/v ≈ 58 ps (black dotted lines of Fig. 5a,b). But experimental data show a variation over more than 90 ps. Considering now the unmatched case, we can determine the sign of the reflection coefficient. For R < 0, a dip below 11 ACCEPTED MANUSCRIPT / CLEAN COPY the plateau value appears near time 2τ a , which is inconsistent with Fig. 3. Thus we conclude that R is positive, R > 0, meaning that the sapphire acoustic impedance is larger than the one of the V 2 O 3 film, Z s > Z.

Analysis

To go further in the analysis of the data, we fitted data obtained from the terahertz and optical pumps with this simplest possible model using the full Eq.2, and the analytical solution of Eq.4.

First, the terahertz penetration depth deserves particular attention, because of the frequency distribution of energy |E ω | 2 of the terahertz excitation pump pulse, shown in Fig. 6. From the optical conductivity measured by TDS, and Drude's formula Eq.1, we calculated the complex refractive index (real and imaginary parts in Fig. 6a), and penetration depth 1/α ω (Fig. 6b). As shown by the figures, these parameters strongly depend on frequency, rising the question of choosing an effective penetration depth in the formula. We model the intensity absorption profile by a weighted average of the absorption profile at each frequency

I(x) ∼ dωP ω α ω e -α ω x ∼ 1 ξ T Hz e -x/ξ T Hz (8)
with probability P ω of absorbing energy within ω and ω + dω defined as

P ω = |E ω | 2 dω|E ω | 2
From the numerical integration of Eq.8, we obtained the effective penetration depth for the terahertz pump ξ T Hz ≈ 1300 nm. In fact, the mean absorption profile is almost indistinguishable to the one at single frequency 0.5 THz over all film thickness.

12 ACCEPTED MANUSCRIPT / CLEAN COPY Second, among the fitting parameters, the photoelastic constants ∂ n ∂ η and ∂ κ ∂ η , where n and κ are respectively the real and imaginary parts of the complex refractive index ñ, are usually difficult to estimate. We tried to fit these parameters to our data, without success. But as already mentioned before, the fact that the sensitivity-function is almost depth independent combines the photoelastic constants with the scaling factor, making difficult or impossible to distinguish both parameters that are strongly correlated. We also verified that a change of the photoelastic constants values could be almost compensated by a change of the scaling factor. Thus, we choose to separate the magnitude from the phase as

∂ ñ ∂ η = ∂ ñ ∂ η e -iφ n .
The photoelastic magnitude ∂ ñ ∂ η is incorporated to the overall scaling factor, while the phase φ n was calculated from the relationship ∂ n ∂ η = -1.63 ∂ κ ∂ η that was used in Refs. 24,25 . Therefore, we ended with only three fitting parameters : the acoustic reflection coefficient R, and two scaling factors to adjust the amplitudes to the terahertz or optical-pump results. All the other parameters were obtained from the bibliography or measured in the laboratory, and fixed to the corresponding values (Table1).

The best fit is presented in Fig. 7, and shows an overall very good agreement with the data for a three-parameter fit. We would like to stress again that the model fits the overall relative optical transmission variation data without any other corrections, including the rising, decreasing and pseudo-plateau parts. In particular, the result of the exponential excitation profile of the optical pump is well reproduced, as well as the small oscillation coming from the reflection of the strain pulse at the film-substrate interface. Interestingly, for the optical pump, it shows that the round shape of the rising part when time<25 ps is also given by the exponential excitation profile, in addition to other possible mechanisms such as non-instantaneous heating, heat diffusion or loss, that are not taken into account in the model. The value of fitted acoustic reflection coefficient R ≈ 0.18 ± 0.02, is larger than the expected value R ≈ 0.11 calculated from densities and 13 ACCEPTED MANUSCRIPT / CLEAN COPY sound velocities extracted from the bibliography (Table1). From R, we obtain the ratio of substrate to film acoustic impedance Z s /Z ≈ 1.44 . It gives an effective acoustic impedance of V 2 O 3 of Z ≈ 31.0 10 5 g/cm 2 /s instead of Z ≈ 35.9 10 5 g/cm 2 /s using density and sound velocity along (110) from bibliography, since our sample has a preferential orientation along (110). This small discrepancy of -14 % may have different origins.

First, a density change due to internal stress/strain created during the synthesis process might be invoked. However, such effects are usually smaller than 1% in V 2 O 3 films (see eg Ref. [START_REF] Kalcheim | Advanced Functional Materials[END_REF] ) and cannot account for the observed 14% differences in R values.

Second, values from the bibliography come usually from single crystal samples, that cannot be directly compared to thin films. As a matter of fact, the thin film sample has a heterogeneous structure, with typical grain sizes estimated to be around 25-30 nm, as revealed by the SEM image of Fig. 1b. This dense mixture of nanocrystals introduces orientational disorder, and changes the effective density because of the presence of holes or porosity. To different orientations correspond different elastic constants and sound velocities. As a result, the measurement over the beamsize would give an effective acoustic impedance by spatial averaging. For instance, in disordered inhomogeneous systems [START_REF] Fouque | Wave propagation and time reversal in randomly layered media[END_REF] , an effective density and elastic stiffness are obtained by averaging the density, and the inverse of the elastic stiffness. It results in an effective sound velocity that can be smaller than the components velocities. In our case, the presence of holes or less dense components decreases the effective density, and thus the acoustic impedance.

Discussion and Conclusion

In this article, the role of the pump profile on the temporal evolution of ultrafast optical time-resolved transmission response of a thin film of vanadium sesquioxide in the metallic phase was studied for two types of pumps: a terahertz pump and an optical pump. Using raw relative transmission data and only three fitting parameters, the temporal profile of the relative optical transmission for both pumps was analysed within the thermoelastic framework over a 250 ps time range, corresponding to the acoustic timescale of this V 2 O 3 film.

From the material point of view, we also measured the acoustic impedance mismatch between V 2 O 3 film and sapphire substrate (R ≈ 0.2). Such a mismatch was unexpected, and gives another important control parameter in addition to usual ones (X-rays diffraction, four-points resistivity measurements, SEM...), to improve the synthesis methodology and better control the quality of the films in the design of new electronic or opto-electronic devices.

We evidenced that the terahertz absorption profile is almost homogeneous over the film thickness (large penetration depth compared to film thickness). Conversely, the optical pump presents a marked exponential decay within the sample whose shape governs the transmission temporal profile at all times. These differences in excitation profiles are directly observed on the shape of the transmission transients. We showed that the superposition of strain waves to the static strain plays 14 ACCEPTED MANUSCRIPT / CLEAN COPY a central role in the interpretation of the data, the former having a contribution that decreases as time increases, the latter determining the plateau value within the acoustic timescale, and for time much smaller than heat transfer timescales. The crucial steps in the data interpretation was to use full raw relative optical transmission data with no time-varying "background" removal. Note also that the transmission sensitivity function being almost flat, it can catch the effect of static strain and dynamical contribution with a small distortion, in contrast to reflectivity sensitivity functions.

Among the few studies of the transmission of V 2 O 3 in the PM phase near room temperature, we only found one dealing with a terahertz pump and an optical probe, Giorgianni et al. 16 However, the transmission data were acquired on a very short time <6 ps, at 175 K, and with a very high field (> 1 MV/cm), making comparison with our data impossible.

Liu et al. 24 studied a 400 nm thick V 2 O 3 film at 220 K, with an optical pump and terahertz probe within 250 ps, and interpreted the data with a constant sensitivity function. The results are shown in the supplementary material of their article, where two acoustic echoes were evidenced. Similar results, but with a reflectivity setup ∆R/R and optical pumps, were also presented by Abreu et al. 25 or Mansart et al. 30 . Note that all these authors only presented the oscillatory part of the signal after removal of some time varying background. We would like to point out that the extraction of the "oscillating" contributions may obscures the presence or absence of the contribution of a static strain. Our analysis demonstrated the importance of this static strain in the interpretation of the plateau observed for both excitation processes.

Recently, a comparison of coherent acoustic phonon generation by terahertz and near infrared ultrashort pump pulses (1.55 eV) on Chromium metallic thin film (≈ 14 nm) with the same probe wavelength (400 nm) was carried out 32 . It is demonstrated that terahertz pulse pumping in metals generates coherent phonons by a thermoelastic process, through an ultrafast Joule effect. For the optical pump, this process arises by an intraband hot electron relaxation by heat transfer to the lattice degrees of freedom. Although the mechanisms differ, the oscillatory parts of the transient optical transmissions overlap almost perfectly, within the experimental time range <30 ps, suggesting similar response of the thin metallic film to both excitation processes. In contrast, our observations with a much thicker thin film of V 2 O 3 , and on much longer times, were clearly able to distinguish between the response to the terahertz and optical pump excitations.

For controlling materials at the macroscopic scale with strong ultrafast pulses, the strain-stress profile is a key feature, and we have shown that terahertz or optical pumps generate different profiles. With this respect, terahertz excitation allows for inducing a more spatially homogeneous excitation over the film depth than optical pump. Thus, we expect thermalization effects that depend on temperature gradients ∂ T ∂ x to be minimized. Moreover, in agreement with the matched layer model 41 , our unmatched layer model also predicts that the strain does not change sign when the excitation is almost uniform over the full film thickness, at least in the range |R| < 0.2 of our simulations, whereas tensile and compressive strain develop inhomogeneously in the film thickness when the intensity decreases exponentially. This change of sign of the strain, or not, within the sample may be of importance in the study of PIPT where the strain or stress-waves are at the origin of the driving mechanisms that trigger the phase transformation. The comprehension of such phenomena may help in tailoring strain profiles to study ultrafast phase transformations of volume changing materials such as V 2 O 3 .

Figure 1

 1 Figure 1 (a) Phase diagram of V 2 O 3 (PM: Paramagnetic Metal; PI: Paramagnetic Insulator ; AFI: Anti ferromagnetic Insulator) and structure in the PM phase. (b) Scanning Electron Microscope (SEM) image. The lm thickness was measured as 210±10 nm. (c) Experimental setups for the terahertz pump or optical pump/optical probe.

Figure 2

 2 Figure 2 Terahertz Time Domain Spectroscopy of V 2 O 3 at room temperature. (a) Electric eld waveform as a function of time : bare substrate (reference scan) and substrate with thin lm sample (sample scan). (b) Experimental real ε 1 and imaginary ε 2 parts of the dielectric function as a function of frequency. The continuous lines are the best t using a Drude model.

5 ACCEPTEDFigure 3

 53 Figure 3 (a) and (b): Relative optical transmission as a function of time ∆T r (t) T r . Raw data are presented. (a) THz pump .(b) Optical pump (540 nm). For t≤1.5 ps the timestep is 0.01 ps, while for t≥2 ps, the time-step is 0.5 ps. (c-d) Terahertz pump electric eld waveform as a function of time, and corresponding squared modulus of the Fourier transform |E ω | 2 . (e-g) Short-time analysis of the optical pump transmission acquisition shown in (b): (e) Raw data (points) and best t (continuous line); (f-g) Evidence of coherent phonon excitation: (f) Resulting transmission after subtraction of the peak and background as a function of time, and (g) the corresponding modulus of the Fourier transform.

7 ACCEPTEDFigure 4

 74 Figure 4 Left: THz pump and right: Optical pump. First line: Excitation proles and modulus of the sensitivity function in the thin lm. Parameters correspond to the real parameters of the experiment. All functions were normalized to one at x=0 for clarity. Second line: Sketch of the thin-lm and substrate interfaces and strain waves.

Figure 5

 5 Figure 5 Temporal evolution of ∆u(t) (a,b) and snapshots of the strain prole η(x,t) (c,d) within the lm and the substrate as a function of position at some selected times. Left column (a,c) and right columns (b,d) correspond respectively to THz and optical pump. Panels a and b represent ∆u(t) with matched layer R=0 (black dotted line) or unmatched R = ±0.2 (R = -0.2 orange dashed line and R = +0.2 green continuous line).Markers are placed at times at which corresponding strain wave snapshots are plotted as a function of position. Panel c and d: strain wave proles η(x,t) in the thin lm (left of the vertical black line) or in the substrate (right of the vertical black line) at a given time, as a function of position. Horizontal dotted gray lines represent zero strain η = 0 line for each time. For the simulation, the sound velocity in the substrate was set to 1.5v.

  Figure (a) Real and imaginary parts of the thin lm complex index at THz frequencies, as given by the TDS experiment. (b) Corresponding THz penetration depth as a function of frequency. The dashed line represents the terahertz pump pulse spectral density |E ω | 2 with arbitrary unit.

Figure 7

 7 Figure 7 Relative optical transmission as a function of time ∆T r (t) T r , and best t results with only three parameters. Experimental raw data are presented. (a) THz pump. (b) Optical pump. Only points every 0.5 ps are shown.

  and Scipy 19 libraries. Figures were drawn with

Table 1

 1 Parameters used in the simulations. i 2 =-1

	Parameter	Value	Reference
	Al 2 O 3 (sapphire)
	ρ	3.98 g/cm 3
	v 001 (sound velocity)	11.23x10 5 cm/s	Winey et al. 35
	Z s (acoustic impedance) 44.7x10 5 g/cm 2 /s
	n (visible)	1.76	Tropf et al. 36
	n (THz)	3.07	Grischkowsky et al. 37
	V 2 O 3 (PM)
	ρ	4.87 g/cm 3
	v 100 (sound velocity)	7.37x10 5 cm/s	Yelon and Keem 38
	Z (acoustic impedance)	35.9x10 5 g/cm 2 /s
	ñ @540nm	1.73 + 0.54 i	Qazilbash et al. 39
	ñ @700nm	1.74 + 0.58 i	Qazilbash et al. 39
	ñ @740nm	1.72 + 0.61 i	Qazilbash et al. 39
	ξ T Hz	1300 nm	From THz-TDS
	ν p (Poisson coefficient)	0.33
	C p (heat capacity)	3.25 J/K/cm 3	Keer et al. 40

8

ACCEPTED MANUSCRIPT / CLEAN COPY

ACCEPTED MANUSCRIPT / CLEAN COPY

ACCEPTED MANUSCRIPT / CLEAN COPY

Acknowledgements

The authors gratefully acknowledge Agence Nationale de la Recherche for financial support undergrant, ANR-19-CE30-0004 ELECTROPHONE, ANR-19-CE29-0018 MULTICROSS, the University Rennes 1 and CNRS Human Resources Direction, for funding.

We acknowledge Maciej Lorenc and Roman Bertoni for experimental support, and fruitful discussions. We would like to thank Pr. H. Cailleau, and Céline Mariette, for useful discussions.

Notes and references

ACCEPTED MANUSCRIPT / CLEAN COPY

Author Contributions

GH and CO supervised the project. GH performed pump-probe and THz-TDS experiments. A technical support has been provided by NG, LG, GP, MH, CO prior and during experiments. GH and CO performed data analysis. CO developed the analytical model and carried out the simulations. MR, ZK, JT, EJ synthetized and characterized the sample. MC initiated the THz project and designed the THz setup. LG built the THz source during his phD thesis. Financial support has been provided by EC and MC. CO and GH wrote the manuscript, with critical reading from EC and EJ. All authors contributed to discussions.

Conflicts of interest

There are no conflicts to declare.