F Monteghetti
email: florian.monteghetti@onera.fr

G Haine

D Matignon
email: denis.matignon@isae-supaero.fr

Stability of Linear Fractional Differential Equations with Delays: a coupled Parabolic-Hyperbolic PDEs formulation

Keywords: Systems with time-delays, stability of delay systems, fractional-order systems, Lyapunov methods, semigroup and operator theory, diffusive representation

Fractional differential equations with delays are ubiquitous in physical systems, a recent example being time-domain impedance boundary conditions in aeroacoustics. This work focuses on the derivation of delay-independent stability conditions by relying on infinitedimensional realisations of both the delay (transport equation, hyperbolic) and the fractional derivative (diffusive representation, parabolic). The stability of the coupled parabolic-hyperbolic PDE is studied using straightforward energy methods. The main result applies to the vectorvalued case. As a numerical illustration, an eigenvalue approach to the stability of fractional delay systems is presented.

INTRODUCTION

Time-delay and fractional operators are ubiquitous in physics. The former model lossless convection phenomena (think wave equation), while the latter model "viscous" losses (i.e. lossy mechanisms that exhibit a long-memory effect, think heat equation). As a result, they have been studied by many different communities; below is a crude overview of the literature, biased towards stability results.

The theory of time-delay systems has been enjoying a sustained development since the 1960s: let us only cite the classic monograph of [START_REF] Bellman | Differential-Difference Equations[END_REF]. Stability criteria roughly split into two: frequency-domain and timedomain ones. Frequency-domain approaches rely on studying the characteristic equation to locate (or merely count) unstable poles (Michiels and Niculescu, 2014, Chap. 1). In the time domain, a popular method consists in designing a Lyapunov-Krasovskii functional and approximating the corresponding sufficient stability condition as a linear matrix inequality (LMI) that is tractable numerically: see (Fridman, 2014, Chap. 3) and (Briat, 2014, Chap. 5) for an introduction; [START_REF] Seuret | Complete quadratic Lyapunov functionals for distributed delay systems[END_REF] and [START_REF] Baudouin | Stability analysis of a system coupled to a transport equation using integral inequalities[END_REF] for advanced LMIs of arbitrary accuracy. Both visions are unified by remarking that time-delay operators can be realised (in the sense of systems theory) as transport equations (of hyperbolic nature). This fact is commonly used by the time-delay community as a means of computing zeros as eigenvalues (Michiels and Niculescu, 2014, § 2.2), as well as by the PDE community to show well-posedness of the corresponding infinite-dimensional

The first author has been financially supported by the French ministry of defence (Direction Générale pour l'Armement) and ONERA (the French Aerospace Lab).

Cauchy problem (Engel and Nagel, 2000, § VI.6), (Curtain and Zwart, 1995, § 2.4).

The theory of fractional calculus is also well-established. The stability of fractional systems of commensurate order can be fully characterised algebraically: no poles in the sector | arg(s)| ≤ α π 2 , see [START_REF] Matignon | Stability properties for generalized fractional differential systems[END_REF]. Moreover, fractional operators turn out to be convolution operators with a so-called "diffusive" kernel, a fact that enables to recast them into the output of an infinite-dimensional ODE (of parabolic nature), related to the heat equation. See [START_REF] Matignon | An introduction to fractional calculus[END_REF] for an overview chapter on both aspects. Applications of this parabolic realisation include, but are not limited to: numerical simulation [START_REF] Lombard | Diffusive approximation of a time-fractional burger's equation in nonlinear acoustics[END_REF]), well-posedness and stability of fractional ODE [START_REF] Matignon | Asymptotic stability of linear conservative systems when coupled with diffusive systems[END_REF]) and PDE [START_REF] Matignon | Asymptotic stability of Webster-Lokshin equation[END_REF]).

The theory of fractional delay differential equations is more recent and limited. A necessary and sufficient stability condition in the frequency domain is established in [START_REF] Bonnet | Analysis of fractional delay systems of retarded and neutral type[END_REF] for the commensurate delay case. (A sufficient condition for the non-commensurate case is proven in [START_REF] Deng | Stability analysis of linear fractional differential system with multiple time delays[END_REF], provided that the delay and fractional operators are not composed.) A corresponding numerical method and its MATLAB implementation are available, see [START_REF] Fioravanti | A numerical method for stability windows and unstable root-locus calculation for linear fractional timedelay systems[END_REF] and [START_REF] Avanessoff | Yalta: a matlab toolbox for the H ∞ -stability analysis of classical and fractional systems with commensurate delays[END_REF]. More recently, a frequency-sweeping approach that counts unstable poles has been proposed for systems of retarded type, see [START_REF] Zhang | Complete stability of linear fractional order time delay systems: A unified frequency-sweeping approach[END_REF]. This work focuses on the derivation of delay-independent stability results for linear fractional differential equations with delays, making use of the language of systems theory and PDEs. The method of proof is based on realisations in the sense of (infinite-dimensional) systems theory: hyperbolic realisation for the time-delay and parabolic realisation for the fractional derivative. This article is organised as follows. Sec. 2 presents a motivation for this work. Sec. 3 works out a scalar toy model to introduce the method of proof, and contrasts it with a frequency-domain method. Sec. 4 proves the main theoretical result of this paper, namely Thm. 7. Sec. 5 proposes an eigenvalue approach to study the stability of fractional delay systems, and provides numerical illustrations. Sec. 6 covers possible extensions and open questions.

AN APPLICATION IN AEROACOUSTICS

The purpose of this section is to present an application of fractional delay systems in aeroacoustics. Additional material, including numerical simulations that rely on the hyperbolic-parabolic realisation, will be presented at the conference. Context. In order to reduce the noise emitted by jet engines, their inlets are lined with a sound absorbing material (known as a liner ). The most widespread material is the so-called "single degree of freedom" (SDoF) liner, depicted in Fig. 1. Based on the Helmholtz resonator, it consists of a perforated plate and a hard-backed cavity. Physically, both components have a different purpose: the (short and narrow) perforations, through visco-thermal losses, account for the majority of sound absorption; the (long and large) cavities drive the resonant frequencies, through a lossy back-and-forth propagation.

Boundary condition. Acoustic liners are typically modelled as single-input single-output systems: p = ẑ ûn , where p, ûn , and ẑ respectively denote the Laplace transforms of acoustic pressure, (inward) normal velocity, and impedance [START_REF] Allard | Propagation of Sound in Porous Media[END_REF]). For the SDoF liner depicted in Fig. 1, a broadly-applicable model reads [START_REF] Monteghetti | Design of broadband time-domain impedance boundary conditions using the oscillatory-diffusive representation of acoustical models[END_REF])

ẑ(s) = perforated plate a 0 + a α s ĥα (s) + a 1 s + coth(b 0 + b α s ĥα (s) + b 1 s) hard-backed cavity , (1)
where ĥα := s → s -α , with α = 1 /2. The coefficient a 0 (b 0 ) models frequency-independent losses in the perforation (cavity), while a α (b α ) models frequency-dependent losses.

A study of (1), based on complex analysis, leads to the following causal fractional delay differential equation in the time domain kernel, i.e. that admits the following representation in the time domain

p(t) = a 0 u n (t)+a 1 un (t)+a α d α C u n (t)+D 2 (u n )(t-τ ), (
D 2 (u n )(t) = n∈Z r n e snt + ˆ∞ 0 e -ξt dµ(ξ) u n (t),
where the residues r n , poles s n and diffusive weight µ are linked to b 0 , b α and b 1 . More generally, (2) can model a wide range of material, not limited to (1), provided that the fractional derivative d α C is replaced by a suitable diffusive operator, D 1 , expressed similarly to D 2 . PDE. To model sound absorption in a duct, ( 2) is used as a (time-domain impedance) boundary condition of a PDE on (p, u) such as the linearised homentropic Euler equations (u 0 and γ > 1 denote the base flow velocity profile and specific heat ratio, respectively)

∂ t p + ∇ • u + u 0 • ∇p + γ p∇ • u 0 = 0 ∂ t u+∇p+[u 0 •∇]u+[u•∇]u 0 +p[u 0 •∇]u 0 = 0.
(3)

A spatial discretisation of ( 2)-( 3) then leads to 4) where X ∈ R n denotes the DoF of the spatial discretisation and C • X the DoF that belong to the impedance boundary (a subset of X). In summary, an aeroacoustical problem (sound absorption in a duct) leads to the fractional delay differential equation (4).

M Ẋ + K X = F 1 d α C (C • X) + F 2 D 2 (C • X)(• -τ ), (

A TOY MODEL, AS WORKED-OUT EXAMPLE

The purpose of this section is to derive sufficient delayindependent stability condition for the following scalar model (0 < α < 1) ẋ(t) = ax(t) + bx(t -τ ) -gd α C x(t) for t > τ (5) x(t) := x 0 (t) for 0 ≤ t ≤ τ .

(6) To derive general results, we consider complex parameters (a, b, g) ∈ C 3 . The interest of this case is twofold: firstly, it acts as a stepping stone to the more general result given in Thm. 7; secondly, it enables to contrast the frequencydomain and time-domain (energy) methods. Remark 1. Due to the delay τ in the evolution equation (5), the given initial datum is a function on the time interval [0, τ ]. The choice of [0, τ ], which differs from the traditional convention for time-delay systems, is justified by the fact that fractional derivatives (be they Caputo or Riemann-Liouville) are naturally defined and related to one another within a causal setting, which implies zero values before t = 0.

Analysis in the Laplace domain

Let us write g = |g| exp(i θ g ), and denote α := 1 -α. Theorem 2. Under the following algebraic condition (a) < -|b| ≤ 0 , (7) and θ g ∈ J α := [-α π 2 , α π 2 ], system (5)-( 6) is delayindependent stable. Remark 3. Let us first inspect the two limiting cases for J α : when α = 1, J 1 = {0} only, and the term (1 + |g|) appears in front of ẋ in the left-hand side, giving rise to a delay-independent stable delay equation without any fractional derivative; whereas when α = 0,

J 0 = [-π 2 , π 2 ],
and the term (a-g) appears in front of x in the right-hand side: since (-g) = -|g| cos(θ g ) ≤ 0, delay-independent stability holds.

Proof. For c ∈ R, let us define the Laplace transform L c (f ) := ´∞ c f (t) e -st dt, and denote f := L 0 (f ). Firstly, take the Laplace transform of (5), taking care of the fact that it is defined for t > τ : L τ (( 5)). Secondly, simplify L τ ( ẋ), L τ (x(•-τ )) and L τ (d α C x) to recover terms that uses only the Laplace transform L 0 . Eventually, the Laplace transform of ( 5)-( 6) reads x(s) = (s -a) ĥ(s) x 0 (s)

(a) +g ĥ(s) L 0 [d α C x 0 1 [0,τ ] ](s) (b) + gx 0 (0) ĥ(s) s α (c)
+gx 0 (τ ) ĥ(s)s α e -sτ (d)

+ x 0 (τ ) ĥ(s)e -sτ (e)
, where ĥ(s) := (s -a -b exp(-τ s) + g s α ) -1 and 1 is the characteristic function. The function ĥ is defined in the open right half-plane C + β for some real abscissa β; since the fractional power of s requires a cut to be performed from the branching point s = 0, necessarily β ≥ 0. Let us check its analyticity. Let s = x + i y = r exp(iθ) with {x ≥ 0 and y ∈ R} or {r ≥ 0 and |θ| ≤ π 2 }. We have

(s -a -b e -τ s + g s α ) = x -(a) -|b| exp(-τ x) cos(τ y -θ b ) + (g s α ) ≥ x -(a) -|b| + |g| |s| α cos(θ g + α θ) ≥ x -(a) -|b| ≥ -(a) -|b| ,
where we used the assumption θ g ∈ J α := [-α π 2 , α π 2 ], which yields cos(θ g + α θ) ≥ 0 for any |θ| ≤ π 2 . If additionally (a) < -|b| then ĥ is analytic in the closed right half-plane C + 0 . Let us now investigate the asymptotic behaviour of x for t → ∞. The limit of the terms (c), (d) and (e) can be studied using the final value theorem, which yields a null value. The terms (a) and (b) require a more involved treatment. In the time domain, (b) reads h

(d α C x 0 1 [0,τ ] ), where h ∈ A(0) (see § B.2) and d α C x 0 1 [0,τ ] is compactly supported in [0, τ ].
Using the decomposition (B.2) and the Lebesgue dominated convergence theorem enables to conclude that (h (

d α C x 0 1 [0,τ ] ))(t) → 0.
The term (a) yields p x 0 , where p ∈ A(0) and x 0 is compactly supported in [0, τ ]. Therefore the exact same arguments enable to conclude that (p x 0 )(t) → 0. 2

Energy analysis with parabolic-hyperbolic realisation

In this section, we use the parabolic-hyperbolic realisation to establish the following result. Proposition 4. Let g > 0. Under the algebraic condition (a) < -|b| ≤ 0, the system (5)-( 6) with x 0 (0) = 0 is delay-independent stable. Remark 5. Here, the parabolic-hyperbolic realisation leads to a less general result than Thm. 2. A better result might be achieved by using less stringent estimates; however, this proof is presented only as a step towards the vector-valued case (see the conditions on G in Thm. 7).

Proof. We begin by using (B.1) to recast (5) using only the Riemann-Liouville fractional derivative ẋ(t) = a x(t) + b x(t -τ ) -g D α RL x(t) + g w(t), (8) where w(t) = x 0 (0) h 1-α (t) = 0. For the sake of compactness, let us denote x τ (•) := x(• -τ ). The natural energy functional is E x := 1 2 |x| 2 , and its decay rate along the trajectories is given by Ėx = 2 (a)

E x + (x (b x τ -g D α RL x + g w)), (9)
whose sign is indefinite. However, by using suitable realisations, energy decay can be proven. The proof is split into three steps: hyperbolic realisation of the delay x τ in § 3.2.1; parabolic realisation of the fractional derivative D α RL in § 3.2.2; study of the global energy decay in § 3.2.3.

Transport PDE for the time-delay operator

Let us realise the delay term through a transport PDE on the 1D spatial domain z ∈ (0, ) with := c τ .

∂ t ψ(t, z) = -c ∂ z ψ(t, z) (0 < z < ) (10) ψ(t, z = 0) := x(t) (11) x(t -τ ) = ψ(t, z = ) . (12)
The initial data is ψ(t = τ, z) := x 0 (τ -z /c). It is natural to examine the energy

E ψ (t) := 1 2 ˆ 0 |ψ(t, z)| 2 dz . (13)
A careful computation (all quantities are complex-valued) leads to the following energy balance

d dt E ψ (t) = -c ˆ 0 (∂ z ψ(t, z) ψ(t, z)) dz (14) = - c 2 [ |ψ(t, z)| 2 ] 0 (15) = c 2 (|x(t)| 2 -|x(t -τ )| 2 ) . (16)

Diffusive realisation for the fractional derivative

For the Riemann-Liouville fractional derivative, the main idea is to use the diffusive representation, which enables to get an energy balance and help control the cross-product term (x D α RL x) that does not have a definite sign. A realisation of the Riemann-Liouville fractional derivative is provided in § A.2, see (A.3-A.4). The energy of the additional variables is naturally defined as

E ϕ (t) := 1 2 ˆ∞ 0 ξ | ϕ(ξ, t)| 2 µ α (ξ) dξ . ( 17
)
This realisation enjoys the following energy balance (which expresses the dissipativity of

D α RL ) d dt E ϕ = (x D α RL x) - ˆ∞ 0 |x -ξ ϕ(ξ, •)| 2 µ α (ξ) dξ . (18)

Lyapunov stability of the coupled system

The energy E of the coupled system can be built from the three energies encountered above: original E x , hyperbolic E ψ , and parabolic E ϕ .

E(t) := E x (t) + k E ψ (t) + g E ϕ (t) . (19)
The coefficient k > 0 is a degree of freedom to be tuned later. The global energy balance reads:

Ė = Ėx + k Ėψ + g Ė ϕ .
(20)

The cross term g (x D α RL x) from the original energy balance (9) cancels out with that from the parabolic energy balance (18). The further use of the hyperbolic energy balance ( 16) leads to

Ė = -X H Σ k X -g ˆ∞ 0 |x -ξ ϕ(ξ, •)| 2 µ α (ξ) dξ, ( 21
)
where X := (x, x τ ) and

Σ k :=    (a) + k c 2 b 2 b 2 -k c 2    . ( 22
)
To conclude the proof, it remains to choose k such that Σ k is positive definite. Let us first recall the following elementary lemma.

Lemma 6. A 2 × 2 hermitian matrix u v v w is positive definite if and only if u > 0 and |v| 2 < u w.

Proof. It suffices to note that a 2 × 2 hermitian matrix Σ is positive definite if and only if det(Σ) > 0 and tr(Σ) > 0.

(Product and sum of the two eigenvalues.) 2

The hermitian matrix Σ k is positive definite if and only if 0 < k < -2 (a)/c and |b| 2 < kc (-2 (a) -kc). For the optimal value k * = -(a)/c, the least stringent condition on the parameters is obtained: (a) < -|b| ≤ 0 . 2

A GENERAL RESULT ON STABILITY OF VECTOR-VALUED LINEAR FRACTIONAL DIFFERENTIAL EQUATIONS WITH DELAY

The general vector-valued model under study is ẋ(t) = Ax(t) + Bx(t -τ ) -Gd α C x(t) for t > τ (23) x(t) := x 0 (t) for 0 ≤ t ≤ τ.

(24) All the matrices involved above are n × n matrices. Theorem 7. Let G be a diagonalisable matrix with eigenvalues (g

1 , • • • , g n ) ≥ 0. Under the algebraic condition max a∈σ(A) (a) < - max b∈σ(B H B) |b| ≤ 0, ( 25
)
the system ( 23)-( 24) with x 0 (0) = 0 is delay-independent stable. (Here B H and σ denote the transconjugate of B and point spectrum, respectively.) Remark 8. If B is diagonalizable (on C), then an elementary computation shows that max b∈σ(B H B) |b| = max b∈σ(B) |b| 2 , which is the square of the spectral radius of B. The algebraic condition (25) then reads max

a∈σ(A) (a) < -max b∈σ(B) |b| ≤ 0.
When G = 0, several type of sufficient conditions are known, see for instance [START_REF] Kharitonov | Robust stability analysis of time delay systems: a survey[END_REF]; [START_REF] Engel | One-parameter semigroups for linear evolution equations[END_REF]; [START_REF] Fridman | Introduction to Time-Delay Systems[END_REF] and references therein. Such a condition typically reads: for any (s) ≥ 0 and |z| ≥ 1, det(s I -A -B z -1 ) = 0. The challenge in verifying such criterion lies in the need to localize zeros. As an alternative, we herein focus on a so-called energy approach that includes the case G = 0.

Proof. Without loss of generality, we can assume that G is diagonal (otherwise, consider a change of basis on ( 23)-( 24)). The proof is similar in spirit to that of Prop. 4. The delay x τ is realised as a vector-valued transport equation on ψ (see ( 10)), while the fractional derivative D α RL x is realised through its diffusive representation. The global energy is defined as (contrast with ( 19))

E(t) := i∈ 1,n E xi (t) + k E ψi (t) + g i E ϕi (t).
Along the trajectories, we get:

Ė = ẋ • x + kc 2 [ x 2 -x(• -τ ) 2 ] + i∈ 1,n g i Ė ϕi . ( 26
)
As G is diagonal, there is no coupling between the diffusive variables ϕ i , and (26) leads to an energy balance akin to ( 21): the coupling between the delay and fractional operators is straightforward. Since g i ≥ 0, it is sufficient for delay-independent stability to prove that

-Σ k :=    A + A H 2 + kc 2 I 1 2 B 1 2 B H - kc 2 I    < 0,
where I denotes the identity matrix on C n , c is the convection velocity and k > 0 is a parameter to be tuned. Let us denote A S = (A+A H ) /2 the symmetric part of A. We have for all x and y in C n

-Σ k x y • x y =A S x•x+ kc 2 x 2 + 1 2 By•x+ 1 2 B H x•y- kc 2 y 2 . Now, for any ε > 0, from 1 2 By • x + 1 2 B H x • y = (By • x) ≤ ε 2 By 2 + 1 2ε x 2 ,
and

By 2 = B H By • y ≤ max b∈σ(B H B) |b| y 2 , we can choose k * = k ε = ε max b∈σ(B H B) |b| /c > 0 to get -Σ k x y • x y ≤ max a∈σ(A) (a) + k ε c 2 + 1 2ε x 2 .
Taking the least stringent value of ε > 0, we derive max

a∈σ(A) (a) + max b∈σ(B H B)
|b| < 0 as a delay-independent stability sufficient condition. 2 Remark 9. If the three matrices A, B and G are assumed simultaneously diagonalizable, the proof of Thm. 7 can be performed directly in the Laplace domain, as for the toy model above. However, the Laplace transform does not appear fruitful to deal with the general case.

Remark 10. The proof breaks down if the delay and fractional operator are composed, i.e. d α x(t -τ ) in ( 23). Nonetheless, this case can be tackled numerically, see § 5.

AN EIGENVALUE APPROACH TO STABILITY OF FRACTIONAL DELAY SYSTEMS

The method of proof, based on parabolic-hyperbolic realisations, naturally suggests a numerical method to discretise (notice the composition for

τ α = 0) ẋ(t) = A x(t) + B x(t -τ ) -G d α C x(t -τ α ) (27) into the Cauchy problem on C n Ẋh (t) = A h X h (t), (28)
where X h := (x, ψ h , ϕ h ) is the extended state. A h is a finite-dimensional approximation of the semigroup generator A (whose existence is herein assumed). To study the stability of (27) using σ(A h ), we assume spectral convergence, without lack of approximation or spectral pollution. A numerical method enables to both check the stability conditions provided above and investigate more intricate systems where it is the only way to conclude. An efficient numerical method will be presented at the conference; some early results are presented below.

Spectral structure. Fig. 2 plots σ(A), with A and B chosen to verify the condition (25) of Thm. 7. The spectrum has two components: an essential one on (-∞, 0), associated with the diffusive representation of the fractional derivative; a discrete one, typical of time-delay operators. Setting g > 0 (g < 0) has a (de)stabilizing effect.

Delay-dependent stability. The scalar case is convenient to perform a parametric study with a complex-valued g. Thm. 2 states that for (a,b) that verify the condition (7), delay-independent stability is achieved for g within the sector J α . Nonetheless, the τ -dependent stability region is larger than this sector. This is illustrated in Fig. 3 that shows, for a given modulus |g| and delay τ , the angle θ g above which ( 5) is unstable, denoted θ max g .

Composition.

Although the presented proofs break down for τ α = 0, this case can be studied numerically. while the reverse is true for a larger value g = |a| (in contrast with the phenomenon highlighted in Fig. 2).

EXTENSIONS AND OPEN QUESTIONS

Multiple delay case

A straightforward extension of Thm. 7 is the case of multiple, non-commensurate delays. Each delay x(t -τ j ) is realised through its own transport equation: in (10), a diagonal velocity matrix C then comes into play.

Semigroup formulation

The proofs of stability given in Prop. 4 and Thm. 7 would be complete if one could invoke LaSalle's invariance principle (Cazenave and Haraux, 1998, Thm. 9.2.3). However, its rigorous application to infinite-dimensional systems requires to check a priori the precompactness of the trajectories in the energy space. This is not a mere formality: the diffusive representation does induce a lack of precompactness, as pointed out in [START_REF] Matignon | Asymptotic stability of linear conservative systems when coupled with diffusive systems[END_REF] for an ODE with fractional derivative, and in [START_REF] Matignon | Asymptotic stability of Webster-Lokshin equation[END_REF] for a PDE with fractional derivative. In both these references, the authors turn to the Arendt-Batty stability theorem to prove the asymptotic stability of the coupled system, which requires a refined analysis of the spectrum of the generator of the underlying semigroup. Similar techniques applied to the operator A that generates the coupled system (23) can be expected to yield well-posedness and, under the given algebraic conditions, asymptotic stability.

Time-domain method with the Caputo derivative

If x 0 (0) = 0, then d α C is different from D α RL , see (B.1). As a result, w = 0 in (9) and the given proofs for Prop. 4 and Thm. 7 break down (i.e. energy decay is not achieved using the diffusive representation of D α RL (A.3)-(A.4)). To circumvent this, one may consider directly using the diffusive representation of d α C given by (B.3)-(B.4). However, since ϕ(•, 0) = 0, the energy ( 17) is then infinite: the energy method still breaks down. A suitable energy space is proposed in § B.3; nonetheless, this norm does not yield an energy balance analogous to (18). These difficulties stem from the fact that the Caputo derivative is not a dissipative operator (see (Lozano et al., 2007, Def. 4.7)) if x 0 (0) = 0, hence the breakdown of the energy method. To the best of the authors' knowledge, a time-domain method of proof that do not need the assumption x 0 (0) = 0 is unavailable.

Eigenvalue approach to stability

As pointed out in § 5, the numerical method can be used to investigate the stability of more intricate cases. Moreover, the spectral properties of the numerical method have only been assumed so far. This could be further investigated.

  Fig. 1. Example of SDoF liner: perforated plate and hardbacked cavity.

Fig. 2 .

 2 Fig. 2. Influence of g on the spectrum σ(A) for (23) with (α, τ ) = ( 1 2 , 10), A = 1 2 -3 1 1 -3 , B = 1 4 1 1 1 1 , and G = g I. ( ): g = 0. ( ): g = +2. ( ): g = -2. ( ): essential spectrum (g = 0 only).

Appendix A. FRACTIONAL DERIVATIVES AND DIFFUSIVE REPRESENTATIONS

A.1 Fractional integral

Let u ∈ L 2 (0, T ), for any β ∈ (0, 1), define the causal convolution kernel h β (t) := Γ(β) -1 t β-1 for t > 0, then the fractional integral of order β ∈ (0, 1) of u is defined by I β u := h β u.

Since h β (t) = ´∞ 0 µ β (ξ) e -ξ t dξ, with specific weight µ β (ξ) = π -1 sin(β π) ξ -β , the fractional integral can be reformulated by the following input-output representation

with e ξ (t) := e -ξ t , and [e ξ u](t) = ´t 0 e -ξ (t-τ ) u(τ ) dτ .

The following infinite-dimensional dynamical system can be seen as a state-space realisation of the fractional integral of order β ∂ t ϕ(ξ, t) = -ξ ϕ(ξ, t) + u(t), ϕ(ξ, 0) = 0 , (A.1)

The functional state space for the diffusive variables ϕ(•, t) is H β = L 2 (R + , µ β dξ), see e.g. [START_REF] Matignon | Asymptotic stability of Webster-Lokshin equation[END_REF]. It defines a well-posed linear system, see [START_REF] Matignon | Standard diffusive systems are well-posed linear systems[END_REF].

A.2 Riemann-Liouville fractional derivative

For short, the Riemann-Liouville fractional derivative of order α ∈ (0, 1) of u ∈ H 1 (0, T ) is defined by y = D α u = D[I 1-α u], and a careful computation shows that the following input-output representation holds

The following infinite-dimensional dynamical system can be seen as a state-space realisation of the fractional derivative of order α

The functional state space for the diffusive variables [START_REF] Matignon | Asymptotic stability of Webster-Lokshin equation[END_REF]. Compared with, for instance, the Grünwald-Letnikov formula, the computational interest of (A.3)-(A.4) is that it is a time-local representation of the hereditary operator D α .

Appendix B. CAPUTO FRACTIONAL DERIVATIVE AND INITIAL CONDITIONS

B.1 Caputo and Riemann-Liouville fractional derivatives

A careful computation shows the following relation between the two fractional derivatives

showing more regularity for d α C than D α RL , as pointed out in [START_REF] Matignon | An introduction to fractional calculus[END_REF]. Since the infinite-dimensional dynamical system (A.3)-(A.4) provides a realisation for D α RL , it thus proves necessary to adapt to d α C . Hopefully, since only the initial condition on u does play a role, a slight modification of the initial state ϕ(., 0) in (A.3) will do. Indeed, following e.g. (Lombard and Matignon, 2016, Prop. 2.1), we find that ϕ 0 (ξ) := u(0 + ) /ξ is the appropriate initial datum. However, since ϕ 0 / ∈ H α , there is a theoretical need to extend the framework, see § B.3.

B.2 Making use of the Callier-Desoer class

The proof of Thm. 2 relies heavily on Thm. A.7.49 of [START_REF] Curtain | An Introduction to Infinite-Dimensional Linear Systems Theory[END_REF]. Let us denote ĥ(s) := (s -a -b exp(-τ s) + g s α ) -1 , we want to prove that ĥ ∈ A(0), the Callier-Desoer class of transfer functions, when (a) < -|b| ≤ 0 and

Then, it is enough to prove that inf (s)≥0 | f (s)| > 0 to ensure that 1/ f ∈ A(0) also. Finally, the product of transfer functions is stable in A(0).

Indeed, since inf (s)≥0 |s -a| = -(a) > 0, we rather compute (s -a) f (s) = s -a -b exp(-τ s) + g s α , and examine its real part with s = x + i y = r exp(iθ), when {x ≥ 0 and y ∈ R}, or {r ≥ 0 and |θ| ≤ π 2 }. Recall, from the proof of Thm. 2, that ((s -a) f (s)) ≥ -(a) -|b| > 0 , where have used

As a consequence, h ∈ A(0), i.e. can be decomposed into

Then, the Lebesgue dominated convergence theorem can be applied, and for any w ∈ L 1 loc (R + ), if lim t→∞ w(t) = 0, then (h w)(t) → 0, as t → ∞. We particularise the computation with w(t) ∝ t -α . 2

B.3 Diffusive representation of fractional derivatives of Caputo type with non-null initial condition

The following infinite-dimensional dynamical system can be seen as a state-space realisation of the Caputo fractional derivative of order α ∂ t φ(ξ, t) = -ξ φ(ξ, t) + u(t), φ(ξ, 0) = u(0 + ) /ξ , (B.3)

The main difficulty is to give a functional framework such that these equations make sense. Indeed, as already said, φ(ξ, 0) = u(0 + ) /ξ / ∈ H α . Let us now define

Then, one can show that for all inputs u ∈ H 1 (0, T ), φ(ξ, 0) = u(0 + ) /ξ ∈ Ȟα and 1(ξ)u(t) ∈ H 1 loc ((0, ∞), Ȟα,-1 ). Finally from [START_REF] Engel | One-parameter semigroups for linear evolution equations[END_REF] and (Tucsnak and Weiss, 2009, Thm. 4.1.6), there exists a unique solution to (B.3) satisfying φ ∈ C([0, T ); Ȟα ) ∩ C 1 ([0, T ); Ȟα,-1 ). From straightforward computation and the following Duhamel's formula φ(ξ, t) = e -ξt u(0 + ) ξ + ˆt 0 e -ξ(t-τ ) u(τ ) dτ , one can easily check, by integration by parts, that [e ξ u](t) = u(t) -ξ φ(ξ, t) , which implies for (B.4) that y(t) = ˆ∞ 0 ˆt 0 e -ξ(t-τ ) u(τ ) dτ µ 1-α (ξ) dξ ,

The last equality comes from Fubini theorem, and the result means that y = I 1-α [ u] := d α C u indeed. Then, since u ∈ H 1 (0, T ), we can easily see that for all α ∈ (0, 1), y ∈ L 2 (0, T ); thus, the output makes sense.


Florian Monteghetti
email: florian.monteghetti@onera.fr

Denis Matignon

Estelle Piot

Lucas Pascal

Asymptotic stability of the linearised Euler equations with long-memory impedance boundary condition

Keywords: impedance boundary condition, diffusive representation, stability, discontinuous Galerkin

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

This work focuses on the (dimensionless) homentropic linearised Euler equations (LEEs)

∂ t p + ∇ • u + u 0 • ∇p + γ p∇ • u 0 = 0 ∂ t u+∇p+[u 0 •∇]u+[u•∇]u 0 +p[u 0 •∇]u 0 = 0, (1)
defined on (0, ∞) × Ω, where Ω ⊂ R n is a bounded Lipschitz open subset, p (u) is the acoustical pressure (velocity), u 0 ∈ C ∞ (Ω) n is the (given) base flow, and γ > 1 is the specific heat ratio. On the boundary Γ := ∂Ω (with outward normal n), a socalled acoustical impedance boundary condition is prescribed :

p(t, x) = [z t u • n(•, x)] (t) (x ∈ Γ := ∂Ω), (2)
where the impedance (z ∈ D + (R) ∩ S (R), causal convolution kernel) models a mono-dimensional medium as a continuous linear time-invariant system.

A recent analysis of acoustical models in the time domain [START_REF] Monteghetti | Design of broadband time-domain impedance boundary conditions using the oscillatory-diffusive representation of acoustical models[END_REF] has shown that a wide range of sound absorbing materials and ground layers, assumed locally-reacting, can be modelled by kernels such as (" " is the weak derivative, a 0 ,a 1 ≥ 0):

z = a 0 δ + a 1 δ + D 2 + D 3 (• -τ ), (3)
where τ ≥ 0 and

D i ∈ L 1 loc (0, ∞) is a causal oscillatory-diffusive kernel (I i ⊂ Z countable, poles [s n,i ] < 0, r n,i > 0, µ i positive Borel measure): D i (t) = n∈I i r n,i e s n,i t oscillatory + ˆ∞ 0 e -ξt dµ i (ξ) diffusive , (4)
which models resonances and visco-thermal losses (e.g. fractional kernel D 2 ∝ t -1 /2 ). A key feature of such positive real kernels is that they can be realised (in the sense of systems theory) by a diagonal, dissipative, infinite-dimensional dynamical system. Note that, if τ > 0 in (3), then (2) is a delayed boundary condition, which models wave reflections.

The two sources of instability in (1,2) are the base flow u 0 and the impedance z.

Acoustical case

The acoustical assumption (u 0 = 0) removes hydrodynamic instabilities, but leaves room for purely acoustical ones triggered by the impedance boundary condition [START_REF] Hesthaven | Nodal discontinuous Galerkin methods[END_REF][START_REF] Khamis | Acoustic boundary conditions at an impedance lining in inviscid shear flow[END_REF]. Below, the delayed (τ = 0) and undelayed (τ > 0) cases are successively investigated by recasting the PDE (1,2) into a Cauchy problem on a Hilbert space H:

Ẋ(t) = A X(t), X(0) = X 0 ∈ D(A). (5)
To express A, a time-domain realisation of z in a state-space Θ is needed. The given asymptotic stability results (see Thms. ). Eq. ( 5) is then defined on

H = L 2 (Ω) × (L 2 (Ω)) n × L 2 (Γ; Θ). Example 2. Let z = a 0 δ + D 2 (not D 2 )
, and define the weighted spaces

Φ α(ξ) = L 2 (0, ∞; α(ξ) dµ 2 ). The diagonal, dissipative, infinite-dimensional real- isation of D 2 in Θ = Φ 1 leads to H = L 2 (Ω) × (L 2 (Ω)) n × L 2 (Γ; Φ 1 )
, and (5) then reads:

A X = A   p u ϕ   =   -∇ • u -∇p A ϕ + B u • n   V = H 1 × (H 1 (div) ∩ (H 1 /2 ) n ) × L 2 (Γ; Φ 1+ξ ) D(A)= X ∈ V [A ϕ+B u•n] ∈ L 2 (Γ; Φ 1 ) p |Γ = a 0 u•n+Cϕ (in L 2 (Γ)) ,
where, formally,

(A ϕ)(x, ξ) = -ξϕ(x, ξ) (state operator), (B u • n)(x, ξ) = 1(ξ) u • n(x) (control),
and

(C ϕ)(x) = ´∞ 0 ϕ(x, ξ) dµ 2 (ξ) (observation). Theorem 3. Assume that τ = 0 in (3). If [a 0 ] > 0, a 1 ≥ 0, [s n,i ] < 0, r n,i > 0 and µ i is a positive
Borel measure, then z admits dissipative realisation, and ( 5) has a unique strong solution X, such that X(t) ≤ X 0 for t ≥ 0 and X(t

) t → ∞ 0.
Proof (Sketch). We follow [START_REF] Matignon | Asymptotic stability of Webster-Lokshin equation[END_REF]. The dissipativity of the realisation of z implies that of A. Well-posedness follows from the m-dissipativity of A. With the Fredholm alternative, we show that ρ(A) ⊃ iR * (we use that H s (Ω) ⊂ L 2 (Ω), s > 0, is a compact embedding). Since 0 / ∈ σ p (A), asymptotic stability then follows from the Arendt-Batty theorem.

Remark 4. With an infinite-dimensional realisation of z, the embedding D(A) ⊂ H may not be compact, hence the need to finely inspect ρ(A), as the pre-compactness condition of LaSalle's invariance principle is not straightforward to verify.

Delayed impedance (τ > 0)

The delayed case (τ > 0) can also be recast into (5) using a hyperbolic realisation of the delay through a transport equation, which leads to an additional extension: H = H×L 2 (Γ; L 2 (0, τ ; Θ)). Asymptotic stability then becomes delay-dependent, which is typical of time-delayed linear systems (see [START_REF] Monteghetti | Stability of Linear Fractional Differential Equations with Delays: a coupled Parabolic-Hyperbolic PDEs formulation[END_REF] and references therein). The energy method of Thm. 3 leads to a sufficient stability condition for the pure delay case (i.e. D 3 = a τ δ, not a diffusive kernel).

= a 0 δ + a 1 δ + D 2 + a τ δ(• -τ ).
Proof (Sketch). Similar to Thm. 3. The energy norm on the hyperbolic variables, • L 2 (Γ;L 2 (0,τ ;C)) (here, Θ = C), is tuned so that A is dissipative. [START_REF] Monteghetti | Stability of Linear Fractional Differential Equations with Delays: a coupled Parabolic-Hyperbolic PDEs formulation[END_REF] 2 Aeroacoustical case The aeroacoustical assumption is u 0 = 0 in (1). In the case of a subsonic base flow (|u 0 | < 1), and under stringent assumptions on u and u 0 (which must be, in particular, potential), the energy functional of Cantrell and Hart [1, Eq. ( 64)] can be used to construct a contraction C 0 -semigroup. Without these assumptions, however, there is no energy balance, and the dissipativity of A is lost: well-posedness can only be achieved in a space like "e -µt L 2 (Ω)", for some µ > ω 0 (A) > 0, where ω 0 (A) is the growth rate of A. (This constitutes a difficulty of the LEEs, compared to e.g. the Galbrun equation, see [START_REF] Brazier | Derivation of an exact energy balance for Galbrun equation in linear acoustics[END_REF].) Current research focuses on the identification of instabilities with (2,3), see e.g. [START_REF] Khamis | Acoustic boundary conditions at an impedance lining in inviscid shear flow[END_REF].

Numerical method

Insights into the stability of ( 5) can be gained by a numerical approximation of the temporal growth rate ω 0 (A). A nodal discontinuous Galerkin method [START_REF] Hesthaven | Nodal discontinuous Galerkin methods[END_REF] is used to formulate Ẋh = A h X h +B h X(•-τ ), with X h = (p h , u h , ϕ h ). The time-domain impedance boundary condition (2,3) is enforced through a centred numerical flux that couples the acoustical unknowns (p h , u h ) with the memory variables ϕ h . If τ > 0, finite-dimensional criteria, which rely on e.g. linear matrix inequalities (LMIs) or spectral conditions, are used to assess stability.

  3 and 5), crucially rely on the dissipativity of this realisation. 1.1 Undelayed impedance (τ = 0) Impedances z of increasing complexity can be considered, with Θ either finite or infinite-dimensional: proportional (z = a 0 δ), for which no realisation is required; derivative (z = a 1 δ ), for which Θ = C. For the sake of brevity and clarity, only two simplified examples (compared with (3)) are given below before the statement of the general result.Example 1. Let ẑ(s) be a real rational function, bounded for [s] ≥ 0. If [ẑ(s)] ≥ 0 (passivity), then it can be realised by a dissipative ODE in Θ = R N , with a suitable energy norm (positive real lemma, see[4, § 3.1]

Theorem 5 .

 5 Let a τ ∈ C and a 1 > 0. If [a 0 ] > |a τ |, then the result of theorem 3 extends to the case z

Acknowledgment This research is supported jointly by ONERA and DGA.


High-order time-domain simulation of acoustic impedance models using diffusive representation Florian MONTEGHETTI 1* , Denis MATIGNON [START_REF] Staffans | Well-posedness and stabilizability of a viscoelastic equation in energy space[END_REF] , Estelle Piot 1 and Lucas PASCAL 1 [START_REF] Allard | Propagation of sound in porous media: modelling sound absorbing materials[END_REF] Dpt. for Aerodynamics and Energetics Modeling, ONERA -The French Aerospace Lab, 31055 Toulouse (FRANCE) [START_REF] Staffans | Well-posedness and stabilizability of a viscoelastic equation in energy space[END_REF] Dpt. of Applied Mathematics, ISAE-SUPAERO, 31055 Toulouse (FRANCE) * florian.monteghetti@onera.fr

Context: impedance boundary condition

Euler equation for isentropic acoustical perturbations p,u ∀x ∈ Ω,

∂ t p + ∇ • u + ∇ • pu = 0 ∂ t u + ∂ t pu + ∇ • u u + ∇ • pu u = -∇p (1)
Non-linear terms neglected in the spatial discretisation

n Ω Ω z Γ z Time-domain impedance boundary condition (TDIBC) ∀x ∈ Γ z , p (t, x) = [z(•, x) u • n(•, x)] (t) (2)
Causal and tempered distribution z ∈ D + (R) ∩ S (R) Integro-differential equation between p and u • n Intuition |ẑ| 1 (high impedance) ⇒ u • n = 0 ∂ n p = 0 (rigid wall) Application linear modelling of a locally-reactive propagation medium Ω z Acoustical impedance models

Ω z is an acoustic liner Γ z perforated plate (p) cavity (c) Impedance models 1 ẑ (s) = ẑp (s) + ẑc (s) (3) ẑp (s) = a 0 + a 1 /2 √ s + a 1 s (4) ẑc (s) = coth b 0 + b 1 /2 √ s + b 1 s (5)

Analysis of acoustical impedance models: oscillatory-diffusive representation

Recast of models (4,5) leads to the definition of ĥp and ĥc ẑp (s) = a 0 + a 1 /2 ĥp (s) s + a 1 s z p u n = a 0 u n + a 1 /2 h p un + a 1 un ẑc (s) = 1 + e -∆t s ĥc (s)

z c u n = u n + h c u n (• -∆t) (6)
Analysis in the complex plane reveals both poles s n and cut

C 0 (s) (s) cut C arg( ĥp ) 0 (s) (s) cut C arg( ĥc ) b 0 = 4, b 1 = 4, b 1/2 = 5b 1 0 (s) (s) cut C arg( ĥc ) poles s n b 0 = 0.1, b 1 = 4, b 1/2 = 5b 1 -π - π 2 0 π 2 π
Oscillatory-diffusive representation [START_REF] Staffans | Well-posedness and stabilizability of a viscoelastic equation in energy space[END_REF][START_REF] Hélie | Diffusive representations for the analysis and simulation of flared acoustic pipes with visco-thermal losses[END_REF][START_REF] Monteghetti | Design of broadband time-domain impedance boundary conditions using the oscillatory-diffusive representation of acoustical models[END_REF] in the time domain

h (t) = n∈Z r n e s n t + ˆC e -ξt dµ (ξ) (7)
Convolution is then recast into an observer of an infinite-dimensional state-space representation

h u n = n∈Z r n ψ n + ˆC ϕ ξ dµ (ξ) (8)
where variables ϕ ξ and ψ n follow a first order dynamic

φξ = -ξ ϕ ξ + u n ψn = s n ψ n + u n (9)
! Time-local realisation of a hereditary operator ! Analysis applicable to a wide class of acoustical models 4

Numerical impedance boundary condition

Numerical TDIBC stems from the discretisation of ( 8

N ϕ = 2 N ϕ = 6 ξ min ξ max f (kHz) [ ĥp (s)s
] 0 1 2 3 4 5 6 7 8 9 10

N ψ = 1 N ψ = 4 N ϕ = 2 f (kHz) [ẑ 2 (s)]

Spatial discretisation

! Discontinuous Galerkin method (high spatial order)

Numerical flux for the TDIBC 5 : centred, with fictive state

q z n q z Γ z q Ω Ω z f = A i n i q + q z (q) 2 q = p u Additional variables (ψ n , ϕ k ) n,k at each node of Γ z Spatial discretisation leads to the following global formulation Ẋ (t) = A X (t) + C X (t -∆t) (11) on extended state X = p i u i ψ n ϕ k ! TDIBC integrated at the semi-discrete level If C = 0, ( 11

Validation

Asymptotic stability study of (11) on a monodimensional case ! TDIBC does not impact stability limit as long as ξ max ≤ f DG max

Well-posedness

Admissibility conditions for z: reality, passivity, causality and stability [START_REF] Zemanian | Distribution theory and transform analysis[END_REF] Semigroup theory on extended state 8 p (•, t) u (•, t) ψ n ϕ ξ

  ) DoF: N ψ poles s n and weights rn / N ϕ poles ξ k and weight μk Optimisation of rn and μk in the frequency domain 3 (10) ! Broadband approximation through elementary linear least squares ! No optimisation against experimental data h (jω) h (jω) =

  ) is a delay differential equation (DDE)Time discretisation! High-order time integration with Runge-Kutta method (a ij , b j ) If C = 0, use of a continuous extension 6 (a ij , b i (θ))

Acknowledgments & Bibliography

This research has initially been supported by ISAE-SUPAERO and is now supported jointly by ONERA (The French Aerospace Lab) and DGA