Xavier Erny
email: xavier.erny@univ-evry.fr

Eva Löcherbach
email: eva.locherbach@univ-paris1.fr

Dasha Loukianova
email: dasha.loukianova@univ-evry.fr

Mean eld limits for interacting Hawkes processes in a diusive regime

Keywords: MSC 2010 subject classications: 60K35, 60G55, 60J35 Multivariate nonlinear Hawkes processes, Mean eld interaction, Piecewise deterministic Markov processes

We consider a sequence of systems of Hawkes processes having mean eld interactions in a diusive regime. The stochastic intensity of each process is a solution of a stochastic dierential equation driven by N independent Poisson random measures. We show that, as the number of interacting components N tends to innity, this intensity converges in distribution in the Skorokhod space to a CIR-type diusion. Moreover, we prove the convergence in distribution of the Hawkes processes to the limit point process having the limit diusion as intensity. To prove the convergence results, we use analytical technics based on the convergence of the associated innitesimal generators and Markovian semigroups.

Introduction

Hawkes processes were originally introduced by [START_REF] Hawkes | Spectra of some self-exciting and mutually exciting point processes[END_REF] to model the appearance of earthquakes in Japan. Since then these processes have been successfully used in many elds to model various physical, biological or economical phenomena exhibiting self-excitation or -inhibition and interactions, such as seismology ( [START_REF] Helmstetter | Subcritical and supercritical regimes in epidemic models of earthquake aftershocks[END_REF], (Y. [START_REF] Kagan | Statistical Distributions of Earthquake Numbers: Consequence of Branching Process[END_REF], [START_REF] Ogata | Seismicity Analysis through Point-process Modeling: A Review[END_REF], [START_REF] Bacry | Second order statistics characterization of Hawkes processes and non-parametric estimation[END_REF]), nancial contagion [START_REF] Aït-Sahalia | Modeling nancial contagion using mutually exciting jump processes[END_REF]), high frequency nancial order books arrivals ( [START_REF] Lu | High dimensional Hawkes processes for limit order books Modelling, empirical analysis and numerical calibration[END_REF], [START_REF] Bauwens | Modelling nancial high frequency data using point processes[END_REF], [START_REF] Hewlett | Clustering of order arrivals, price impact and trade path optimisation[END_REF]), genome analysis [START_REF] Reynaud-Bouret | Adaptive estimation for Hawkes processes; application to genome analysis[END_REF]) and interactions in social networks ( [START_REF] Zhou | Learning triggering kernels for multi-dimensional Hawkes processes[END_REF]). In particular, multivariate Hawkes processes are extensively used in neuroscience to model temporal arrival of spikes in neural networks ( [START_REF] Grün | Analysis of Parallel Spike Trains[END_REF], [START_REF] Okatan | Analyzing Functional Connectivity Using a Network Likelihood Model of Ensemble Neural Spiking Activity[END_REF], [START_REF] Pillow | Spatio-temporal correlations and visual signalling in a complete neuronal population[END_REF], [START_REF] Reynaud-Bouret | Goodness-of-Fit Tests and Nonparametric Adaptive Estimation for Spike Train Analysis[END_REF]) since they provide good models to describe the typical temporal decorrelations present in spike trains of the neurons as well as the functional connectivity in neural nets.

In this paper, we consider a sequence of multivariate Hawkes processes (Z N ) N ∈N * of the form Z N = (Z N,1 t , . . . Z N,N t ) t≥0 . Each Z N is designed to describe the behaviour of some interacting system with N components, for example a neural network of N neurons. More precisely, Z N is a multivariate counting process where each Z N,i records the number of events related to the i-th component, as for example the number of spikes of the i-th neuron. These counting processes are interacting, that is, any event of type i is able to trigger or to inhibit future events of all other types j. The process (Z N,1 , . . . , Z N,N ) is informally dened via its stochastic intensity process λ N = (λ N,1 (t), . . . , λ N,N (t)) t≥0 through the relation P(Z N,i has a jump in ]t, t + dt]|F t ) = λ N,i (t)dt, 1 ≤ i ≤ N, where F t = σ Z N s : 0 ≤ s ≤ t . The stochastic intensity of a Hawkes process is given by

λ N,i (t) = f N i   N j=1 t -∞ h N ij (t -s)dZ N,j (s)   . (1)
Here, h N ij models the action or the inuence of events of type j on those of type i, and how this inuence decreases as time goes by. The function f N i is called the jump rate function of Z N,i . Since the founding works of [START_REF] Hawkes | Spectra of some self-exciting and mutually exciting point processes[END_REF] and [START_REF] Hawkes | A Cluster Process Representation of a Self-Exciting Process[END_REF], many probabilistic properties of Hawkes processes have been well-understood, such as ergodicity, stationarity and long time behaviour (see [START_REF] Brémaud | Stability of Nonlinear Hawkes Processes[END_REF], [START_REF] Daley | An Introduction to the Theory of Point Processes: Volume I: Elementary Theory and Methods[END_REF], [START_REF] Costa | Renewal in Hawkes processes with self-excitation and inhibition[END_REF], [START_REF] Raad | Renewal Time Points for Hawkes Processes[END_REF] and [START_REF] Graham | Regenerative properties of the linear Hawkes process with unbounded memory[END_REF]). A number of authors studied the statistical inference for Hawkes processes ( [START_REF] Ogata | The asymptotic behavior of maximum likelihood estimators for stationary point processes[END_REF] and [START_REF] Reynaud-Bouret | Adaptive estimation for Hawkes processes; application to genome analysis[END_REF]). Another eld of study, very active nowadays, concerns the behaviour of the Hawkes process when the number of components N goes to innity. During the last decade, large population limits of systems of interacting Hawkes processes have been studied in [START_REF] Fournier | On a toy model of interacting neurons[END_REF], (Delattre, Fournier and Homann, 2016) and [START_REF] Ditlevsen | Multi-class Oscillating Systems of Interacting Neurons[END_REF].

In (Delattre, Fournier and Homann, 2016), the authors consider a general class of Hawkes processes whose interactions are given by a graph. In the case where the interactions are of mean eld type and scaled in N -1 , namely h N ij = N -1 h and f N i = f in (1), they show that the Hawkes processes can be approximated by an i.i.d. family of inhomogeneous Poisson processes. They observe that for each xed integer k, the joint law of k components converges to a product law as N tends to innity, which is commonly referred to as the propagation of chaos. [START_REF] Ditlevsen | Multi-class Oscillating Systems of Interacting Neurons[END_REF] generalize this result to a multi-population frame and show how oscillations emerge in the large population limit. Note again that the interactions in both papers are scaled in N -1 , which leads to limit point processes with deterministic intensity.

The purpose of this paper is to study the large population limit (when N goes to innity) of the multivariate Hawkes processes (Z N,1 , . . . , Z N,N ) with mean eld interactions scaled in N -1/2 . Contrarily to the situation considered in (Delattre, Fournier and Homann, 2016) and [START_REF] Ditlevsen | Multi-class Oscillating Systems of Interacting Neurons[END_REF], this scaling leads to a non-chaotic limiting process with stochastic intensity. As we consider interactions scaled in N -1/2 , we have to center the terms of the sum in (1) to make the intensity process converge according to some kind of central limit theorem. To this end, we consider intensities with stochastic jump heights. Namely, in this model, the multivariate Hawkes processes

(Z N,i ) 1≤i≤N (N ∈ N * ) are of the form Z N,i t = ]0,t]×R+×R 1 {z≤λ N s } dπ i (s, z, u), 1 ≤ i ≤ N, (2)
where (π i ) i∈N * are i.i.d. Poisson random measures on R + × R + × R of intensity ds dz dµ(u) and µ is a centered probability measure on R having a nite second moment σ 2 . The stochastic intensity of Z N,i is given by λ

N,i t = λ N t = f X N t-,
where

X N t = 1 √ N N j=1 [0,t]×R+×R h(t -s)u1 {z≤f(X N s-)} dπ j (s, z, u).
Moreover we consider a function h of the form h(t) = e -αt so that the process (X N t ) t is a piecewise deterministic Markov process. In the framework of neurosciences, X N t represents the membrane potential of the neurons at time t. The random jump heights u, chosen according to the measure µ, model random synaptic weights and the jumps of Z N,j represent the spike times of neuron j. If neuron j spikes at time t, an additional random potential height u/ √ N is given to all other neurons in the system. As a consequence, the process X N has the following dynamic

dX N t = -αX N t dt + 1 √ N N j=1 R+×R u1 {z≤f(X N t-)} dπ j (t, z, u).
Its innitesimal generator is given by

A N g(x) = -αx g (x) + N f (x) R g x + u √ N -g(x) µ(du),
for suciently smooth functions g. As N goes to innity, the above expression converges to

Āg(x) = -αx g (x) + σ 2 2 f (x)g (x),
which is the generator of a CIR-type diusion given as solution of the SDE

d Xt = -α Xt dt + σ f ( Xt )dB t . (3)
It is classical to show in this framework that the convergence of generators implies the convergence of X N to X in distribution in the Skorokhod space. In this article we establish explicit bounds for the weak error for this convergence by means of a Trotter-Kato like formula. Moreover we establish for each i, the convergence in distribution in the Skorokhod space of the associated counting process Z N,i to the limit counting process Zi which has intensity (f ( Xt )) t . Conditionally on X, the Zi , i ≥ 1, are independent. This property can be viewed as a conditional propagation of chaosproperty, which has to be compared to (Delattre, Fournier and Homann, 2016) and [START_REF] Ditlevsen | Multi-class Oscillating Systems of Interacting Neurons[END_REF] where the intensity of the limit process is deterministic and its components are truly independent, and to [START_REF] Carmona | Mean eld games with common noise[END_REF], [START_REF] Dawson | Nonlinear Differential Equations and Applications NoDEA[END_REF] and [START_REF] Kurtz | Particle representations for a class of nonlinear SPDEs[END_REF] where all interacting components are subject to common noise. In our case, the common noise, that is, the Brownian motion B of (3), emerges in the limit as a consequence of the central limit theorem.

To obtain a precise control of the speed of convergence of X N to X we use analytical methods showing rst the convergence of the generators from which we deduce the convergence of the semigroups via the formula

Pt g(x) -P N t g(x) = t 0 P N t-s Ā -A N Ps g(x)ds. (4)
Here Pt g(x) = E x g( Xt ) and P N t g(x) = E x g(X N t ) denote the Markovian semigroups of X and X N . This formula is well-known in the classical semigroup theory setting where the generators are strong derivatives of semigroups in the Banach space of continuous bounded functions (see Lemma 1.6.2 of [START_REF] Ethier | Markov Processes. Characterization and Convergence[END_REF]). In our case, we have to consider extended generators (see [START_REF] Davis | Markov Models and Optimization[END_REF] or [START_REF] Meyn | Stability of Markovian Processes III: Foster-Lyapunov Criteria for Continuous-Time Processes[END_REF])), i.e. A N g(x) is the point-wise derivative of t → P N t g(x) in 0. The proof of formula (4) for our extended generators is given in the Appendix (Proposition 5.6).

It is well-known that under suitable assumptions on f, the solution of (3) admits a unique invariant measure λ whose density is explicitly known. Thus, a natural question is to consider the limit of the law L(X N t ) of X N t when t and N go simultaneously to innity. We prove that the limit of L(X N t ) is λ, for (N, t) → (∞, ∞), under suitable conditions on the joint convergence of (N, t). We also prove that there exists a parameter α * such that for all α > α * , this converges holds whenever (N, t) → (∞, ∞) jointly, without any further condition, and we provide a control of the error (Theorem 1.6).

The paper is organized as follows: in Section 1, we state the assumptions and formulate the main results. Section 2 is devoted to the proof of the convergence of the semigroup of X N to that of X (Theorem 1.4.(i)), and Section 3 to the study of the limit of the law of X N t as N, t → ∞ (Theorem 1.6). In Section 4, we prove the convergence of the systems of point processes (Z N,i ) 1≤i≤N to ( Zi ) i≥1 (Theorem 1.7). Finally in the Appendix, we collect some results about extended generators and we give the proof of (4) together with some other technical results that we use throughout the paper.

1. Notation, assumptions and main results

Notation

The following notation are used throughout the paper:

• If X is a random variable, we note L(X) its distribution.

• If g is a real-valued function which is n times dierentiable, we note ||g|| n,∞ = n k=0 ||g (k) || ∞ .

• If g : R → R is a real-valued measurable function and λ a measure on (R, B(R)) such that g is integrable with respect to λ, we write λ(g) for R gdλ. • We write C n b (R) for the set of the functions g which are n times continuously dierentiable such that ||g|| n,∞ < + ∞, and we write for short C b (R) instead of C 0 b (R). Finally, C n (R) denotes the set of n times continuously dierentiable functions that are not necessarily bounded nor have bounded derivates.

• If g is a real-valued function and I is an interval, we note ||g|| ∞,I = sup x∈I |g(x)|.

• We write C n c (R) for the set of functions that are n times continuously dierentiable and that have a compact support.

• We write D(R + , R) for the Skorokhod space of càdlàg functions from R + to R, endowed with the Skorokhod metric (see Chapter 3 Section 16 of [START_REF] Billingsley | Convergence of Probability Measures[END_REF]), and D(R + , R + ) for this space restricted to non-negative functions. • α is a positive constant, L, σ and m k (1 ≤ k ≤ 4) are xed parameters dened in Assumptions 1, 2 and 3 below. Finally, we note C any arbitrary constant, so the value of C can change from line to line in an equation. Moreover, if C depends on some non-xed parameter θ, we write C θ .

Assumptions

Let X N satisfy

     dX N t = -αX N t dt + 1 √ N N j=1 R+×R u1 {z≤f(X N t-)} dπ j (t, z, u), X N 0 ∼ ν N 0 , (5)
where ν N 0 is a probability measure on R. Under natural assumptions on f, the SDE (5) admits a unique non-exploding strong solution (see Proposition 5.8).

The aim of this paper is to provide explicit bounds for the convergence of X N in the Skorokhod space to the limit process ( Xt ) t∈R+ which is solution to the SDE

d Xt = -α Xt dt + σ f Xt dB t , X0 ∼ ν0 , (6)
where σ 2 is the variance of µ, (B t ) t∈R+ is a one-dimensional standard Brownian motion, and ν0 is a suitable probability measure on R.

To prove our results, we need to introduce the following assumptions.

Assumption 1.

√ f is a positive and Lipschitz continuous function, having Lipschitz constant L.

Under Assumption 1, it is classical that the SDE (6) admits a unique non-exploding strong solution (see remark IV.2.1, Theorems IV.2.3, IV.2.4 and IV.3.1 of [START_REF] Ikeda | Stochastic Dierential Equations and Diusion Processes[END_REF]).

Assumption 1 is used in many computations of the paper in one of the following forms:

• ∀x ∈ R, f (x) ≤ ( f (0) + L|x|) 2
, or, if we do not need the accurate dependency on the parameter,

• ∀x ∈ R, f (x) ≤ C(1 + x 2 ).
Assumption 2.

• R x 4 dν 0 (x) < ∞ and for every

N ∈ N * , R x 4 dν N 0 (x) < ∞.
• µ is a centered probability measure having a nite fourth moment, we note σ 2 its variance.

Assumption 2 allows us to control the moments up to order four of the processes (X N t ) t and ( Xt ) t (see Lemma 2.1) and to prove the convergence of the generators of the processes (X N t ) t (see Proposition 2.3). Assumption 3. We assume that f belongs to C 4 (R) and that for each

1 ≤ k ≤ 4, ( √ f ) (k) is bounded by some constant m k . Remark 1.1. By denition m 1 = L, since m 1 := ||( √ f ) || ∞ and L is the Lipschitz constant of √ f .
Assumption 3 guarantees that the stochastic ow associated to (6) has regularity properties with respect to the initial condition X0 = x. This will be the main tool to obtain uniform, in time, estimates of the limit semigroup, see Proposition 2.4.

Example 1.2. The functions f

(x) = 1 + x 2 , f (x) = √ 1 + x 2 and f (x) = (π/2 + arctan x) 2 satisfy Assumptions 1 and 3. Assumption 4. X N 0 converges in distribution to X0 .
Obviously, Assumption 4 is a necessary condition for the convergence in distribution of X N to X.

Main results

Our rst main result is the convergence of the process X N to X in distribution in the Skorokhod space, with an explicit rate of convergence for their semigroups. This rate of convergence will be expressed in terms of the following parameters

β := max 1 2 σ 2 L 2 -α, 2σ 2 L 2 -2α, 7 2 σ 2 L 2 -3α (7)
and, for any T > 0 and any xed ε > 0,

K T := (1 + 1/ε) T 0 (1 + s 2 )e βs 1 + e (σ 2 L 2 -2α+ε)(T -s) ds. (8) Remark 1.3. If α > 7/6 σ 2 L 2 , then β < 0, and one can choose ε > 0 such that σ 2 L 2 -2α + ε < 0, implying that sup T >0 K T < ∞.
Recall that Pt g(x) = E x g( Xt ) and P N t g(x) = E x g(X N t ) denote the Markovian semigroups of X and X N .

Theorem 1.4. If Assumptions 1 and 2 hold, then the following assertions are true.

(i) Under Assumption 3, for all T ≥ 0, for each

g ∈ C 3 b (R) and x ∈ R, sup 0≤t≤T P N t g(x) -Pt g(x) ≤ C(1 + x 2 )K T ||g|| 3,∞ 1 √ N .
In particular, if α > 7 6 σ 2 L 2 , then

sup t≥0 P N t g(x) -Pt g(x) ≤ C(1 + x 2 )||g|| 3,∞ 1 √ N .
(ii) If in addition Assumption 4 holds, then

(X N ) N converges in distribution to X in D(R + , R).
We refer to Proposition 2.4 for the form of β given in (7). Theorem 1.4 is proved in the end of Subsection 2.2. (ii) is a consequence of Theorem IX.4.21 of [START_REF] Jacod | Limit Theorems for Stochastic Processes[END_REF], using that X N is a semimartingale. Alternatively, it can be proved as a consequence of (i), using that X N is a Markov process.

Below we give some simulations of the trajectories of the process (X N t ) t≥0 in Figure 1.

Remark 1.5. Theorem 1.4.(ii) states the convergence of X N to X in the Skorokhod topology. Since X is almost surely continuous, this implies the, a priori stronger, convergence in distribution in the topology of the uniform convergence on compact sets. Indeed, according to Skorohod's representation theorem (see Theorem 6.7 of [START_REF] Billingsley | Convergence of Probability Measures[END_REF])), we can assume that X N converges almost surely to X in the Skorokhod space, and this classically entails the uniform convergence on every compact set (see the discussion at the bottom of page 124 in Section 12 of [START_REF] Billingsley | Convergence of Probability Measures[END_REF]).

Under our assumptions, P admits an invariant probability measure λ, and we can even control the speed of convergence of P N t g(x) to λ(g), as (N, t) goes to innity, for suitable conditions on the joint convergence of N and t.

Fig 1. Simulation of trajectories of (X N t ) 0≤t≤10 with X N 0 = 0, α = 1, µ = N (0, 1), f (x) = 1 + x 2 , N = 100 (left picture) and N = 500 (right picture).
Theorem 1.6. Under Assumptions 1 and 2, X is recurrent in the sense of Harris, having invariant probability measure λ(dx) = p(x)dx with density

p(x) = C 1 f (x) exp - 2α σ 2 x 0 y f (y) dy .
Besides, if Assumption 3 holds, then for all g ∈ C 3 b (R) and x ∈ R,

P N t g(x) -λ(g) ≤ C||g|| 3,∞ (1 + x 2 ) K t √ N + e -γt ,
where C and γ are positive constants independent of N and t, and where K t is dened in (8). In particular,

P N t (x, •) converges weakly to λ as (N, t) → (∞, ∞), provided K t = o( √ N ).
If we assume, in addition, that α > 7 6 σ 2 L 2 , then P N t (x, •) converges weakly to λ as (N, t) → (∞, ∞) without any condition on the joint convergence of (t, N ), and we have, for any g ∈ C 3 b (R) and x ∈ R,

P N t g(x) -λ(g) ≤ C||g|| 3,∞ (1 + x 2 ) 1 √ N + e -γt .
Theorem 1.6 is proved in the end of Section 3. Finally, using Theorem 1.4.(ii), we show the convergence of the point processes Z N,i dened in (2) to limit point processes Zi having stochastic intensity f ( Xt ) at time t. To dene the processes Zi (i ∈ N * ), we x a Brownian motion (B t ) t≥0 on some probability space dierent from the one where the processes X N (N ∈ N * ) and the Poisson random measures π i (i ∈ N * ) are dened. Then we x a family of i.i.d. Poisson random measures πi (i ∈ N * ) on the same space as (B t ) t≥0 , independent of (B t ) t≥0 . The limit point processes Zi are then dened by

Zi t = ]0,t]×R+×R 1 {z≤f( Xs)} dπ i (s, z, u). (9)
Theorem 1.7. Under Assumptions 1, 2 and 4, for every k ∈ N * , the sequence (Z N,1 , . . . , Z N,k ) N converges to ( Z1 , . . . , Zk ) in distribution in D(R + , R k ). Consequently, the sequence (Z N,j ) j≥1 converges to ( Zj ) j≥1 in distribution in D(R + , R) N * for the product topology.

Let us give a brief interpretation of the above result. Conditionally on X, for any k > 1, Z1 , . . . , Zk are independent. Therefore, the above result can be interpreted as a conditional propaga- tion of chaos property (compare to [START_REF] Carmona | Mean eld games with common noise[END_REF] dealing with the situation where all interacting components are subject to common noise). In our case, the common noise, that is, the Brownian motion B driving the dynamic of X, emerges in the limit as a consequence of the central limit theorem. Theorem 1.7 is proved in the end of Section 4.

Remark 1.8. In Theorem 1.7, we implicitly dene Z N,i := 0 for each i ≥ N + 1.

Proof of Theorem 1.4

The goal of this section is to prove Theorem 1.4. To prove the convergence of the semigroups of (X N ) N , we show in a rst time the convergence of their generators. We start with useful a priori bounds on the moments of X N and X.

Lemma 2.1. Under Assumptions 1 and 2, the following holds.

(i) For all ε > 0, t > 0 and

x ∈ R, E x (X N t ) 2 ≤ C(1 + 1/ε)(1 + x 2 )(1 + e (σ 2 L 2 -2α+ε)t ), for some C > 0 independent of N, t, x and ε. (ii) For all ε > 0, t > 0 and x ∈ R, E x ( Xt ) 2 ≤ C(1 + 1/ε)(1 + x 2 )(1 + e (σ 2 L 2 -2α+ε)t ), for some C > 0 independent of t, x and ε. (iii) For all N ∈ N * , T > 0, E (sup 0≤t≤T |X N t |) 2 < +∞ and E (sup 0≤t≤T | Xt |) 2 < +∞. (iv) For all T > 0, N ∈ N * , sup 0≤t≤T E x (X N t ) 4 ≤ C T (1 + x 4 ) and sup 0≤t≤T E x ( Xt ) 4 ≤ C T (1 + x 4 ). (v) For all 0 ≤ s, t ≤ T and x ∈ R, E x Xt -Xs 2 ≤ C T (1 + x 2 )|t -s| and E x X N t -X N s 2 ≤ C T (1 + x 2 )|t -s|.
We postpone the proof of Lemma 2.1 to the Appendix. The inequalities of points (i) and (ii) of the lemma hold for any xed ε > 0. This parameter ε appears for the following reason. We prove the above points using the Lyapunov function x → x 2 . When applying the generators to this function, there are terms of order x that appear and that we bound by x 2 ε + ε -1 to be able to compare it to x 2 .

Convergence of the generators

Throughout this paper, we consider extended generators similar to those used in [START_REF] Meyn | Stability of Markovian Processes III: Foster-Lyapunov Criteria for Continuous-Time Processes[END_REF] and in [START_REF] Davis | Markov Models and Optimization[END_REF], because the classical notion of generator does not suit to our framework (see the beginning of Section 5.1). As this denition slightly diers from one reference to another, we dene explicitly the extended generator in Denition 5.1 below and we prove the results on extended generators that we need in this paper. We note A N the extended generator of X N and Ā the one of X, and D (A N ) and D ( Ā) their extended domains. The goal of this section is to prove the convergence of A N g(x) to Āg(x) and to establish the rate of convergence for test functions g ∈ C 3 b (R). Before proving this convergence, we state a lemma which characterizes the generators for some test functions. This lemma is a straightforward consequence of Itô's formula and Lemma 2.1.(i).

Lemma 2.2. C 2 b (R) ⊆ D ( Ā), and for all g ∈ C 2 b (R) and x ∈ R, we have

Āg(x) = -αxg (x) + 1 2 σ 2 f (x)g (x). Moreover, C 1 b (R) ⊆ D (A N )
, and for all g ∈ C 1 b (R) and x ∈ R, we have

A N g(x) = -αxg (x) + N f (x) R g x + u √ N -g(x) dµ(u).
The following result is the rst step towards the proof of our main result.

Proposition 2.3. If Assumptions 1 and 2 hold, then for all

g ∈ C 3 b (R), Āg(x) -A N g(x) ≤ f (x) g ∞ 1 6 √ N R |u| 3 dµ(u).
Proof. For g ∈ C 3 b (R), if we note U a random variable having distribution µ, we have, since

E [U ] = 0, A N g(x) -Āg(x) ≤f (x) N E g x + U √ N -g(x) - 1 2 σ 2 g (x) =f (x)N E g x + U √ N -g(x) - U √ N g (x) - U 2 2N g (x) ≤f (x)N E g x + U √ N -g(x) - U √ N g (x) - U 2 2N g (x)
.

Using Taylor-Lagrange's inequality, we obtain the result.

Convergence of the semigroups

Once the convergence A N g(x) → Āg(x) is established, together with a control of the speed of convergence, our strategy is to rely on the following representation

Pt -P N t g(x) = t 0 P N t-s Ā -A N Ps g(x)ds, (10)
which is proven in Proposition 5.6 in the Appendix. Obviously, to be able to apply Proposition 2.3 to the above formula, we need to ensure the regularity of x → Ps g(x), together with a control of the associated norm || Ps g|| 3,∞ . This is done in the next proposition.

Proposition 2.4. If Assumptions 1, 2 and 3 hold, then for all t ≥ 0 and for all

g ∈ C 3 b (R), the function x → Pt g(x) belongs to C 3 b (R) and satises Pt g ∞ ≤ C||g|| 3,∞ (1 + t 2 )e βt , (11)
with

β = max( 1 2 σ 2 L 2 -α, 2σ 2 L 2 -2α, 7 2 σ 2 L 2 -3α).
Moreover, for all T > 0,

sup 0≤t≤T || Pt g|| 3,∞ ≤ Q T ||g|| 3,∞ (12)
for some Q T > 0, and for all i ∈ {0, 1, 2} and x ∈ R, s → ( Ps g

) (i) (x) = ∂ i ∂x i ( Ps g(x)
) is continuous.

The proof of Proposition 2.4 requires some detailed calculus to obtain the explicit expression for β, so we postpone it to the Appendix.

Proof of Theorem 1.4. Step 1. The proof of point (i) is a straightforward consequence of Proposition 2.3, since, applying formula (10),

Pt g(x) -P N t g(x) = t 0 P N t-s Ā -A N Ps g(x)ds ≤ t 0 E x Ā Ps g X N t-s -A N Ps g X N t-s ds ≤C 1 √ N t 0 Ps g ∞ E x f X N t-s ds ≤C 1 √ N ||g|| 3,∞ t 0 (1 + s 2 )e βs 1 + E x X N t-s 2 ds ≤C 1 + 1 ε 1 √ N ||g|| 3,∞ (1 + x 2 ) t 0 (1 + s 2 )e βs 1 + e (σ 2 L 2 -2α+ε)(t-s) ds,
where we have used respectively Proposition 2.4 and Lemma 2.1.(i) to obtain the two last inequalities above, and ε is any positive constant.

Step 2. We now give the proof of point (ii) of the theorem. With the notation of Theorem IX.4.21 of [START_REF] Jacod | Limit Theorems for Stochastic Processes[END_REF], we have K N (x, dy) := N f (x)µ( √ N dy), b N (x) = -αx + K N (x, dy)y = -αx, and c N (x) = K N (x, dy)y 2 = σ 2 f (x). Then, an immediate adaptation of Theorem IX.4.21 of [START_REF] Jacod | Limit Theorems for Stochastic Processes[END_REF] to our frame implies the result.

Proof of Theorem 1.6

In this section, we prove Theorem 1.6. We begin by proving some properties of the invariant measure of Pt . In what follows we use the total variation distance between two probability measures ν 1 and ν 2 dened by

ν 1 -ν 2 T V = 1 2 sup g: g ∞≤1 |ν 1 (g) -ν 2 (g)|.
Proposition 3.1. If Assumptions 1 and 2 hold, then the invariant measure λ of ( Pt ) t exists and is unique. Its density is given, up to multiplication with a constant, by

p(x) = C 1 f (x) exp - 2α σ 2 x 0 y f (y)
dy .

In addition, if Assumption 3 holds, then for every 0 < q < 1/2, there exists some γ > 0 such that, for all t ≥ 0,

|| Pt (x, •) -λ|| T V ≤ C 1 + x 2 q e -γt .
Proof. In a rst time, let us prove the positive Harris recurrence of X implying the existence and uniqueness of λ. According to Example 3.10 of [START_REF] Khasminskii | Stochastic stability of dierential equations[END_REF] it is sucient to show that S(x) :=

x 0 s(y)dy goes to +∞ (resp. -∞) as x goes to +∞ (resp. -∞), where

s(x) := exp 2α σ 2 x 0 v f (v) dv .
For x > 0, and using that f is subquadratic,

s(x) ≥ exp C x 0 2v 1 + v 2 dv = exp C ln(1 + x 2 ) = (1 + x 2 ) C ≥ 1,
implying that S(x) goes to +∞ as x goes to +∞. With the same reasoning, we obtain that S(x) goes to -∞ as x goes to -∞. Finally, the associated invariant density is given, up to a constant, by

p(x) = C f (x)s(x)
.

For the second part of the proof, take V (x) = (1 + x 2 ) q , for some q < 1/2, then

V (x) = 2qx(1 + x 2 ) q-1 , V (x) = 2q(1 + x 2 ) q-2 [2x 2 (q -1) + (1 + x 2 )].
As q < 1 2 , V (x) < 0 for x 2 suciently large, say, for |x| ≥ K. In this case, for |x| ≥ K,

ĀV (x) ≤ -2αqx 2 (1 + x 2 ) q-1 ≤ -2αq x 2 1 + x 2 V (x) ≤ -2qα K 2 1 + K 2 V (x) = -cV (x).
So we obtain that, for suitable constants c and d, for any x ∈ R, ĀV (x) ≤ -cV (x) + d.

(13)

Obviously, for any xed T > 0, the sampled chain ( XkT ) k≥0 is Feller and λ-irreducible. The support of λ being R, Theorem 3.4 of [START_REF] Meyn | Stability of Markovian Processes III: Foster-Lyapunov Criteria for Continuous-Time Processes[END_REF] implies that every compact set is petite for the sampled chain. Then, as (13) implies the condition (CD3) of Theorem 6.1 of [START_REF] Meyn | Stability of Markovian Processes III: Foster-Lyapunov Criteria for Continuous-Time Processes[END_REF], we have the following bound: introducing for any probability measure µ the weighted norm

µ V := sup g:|g|≤1+V |µ(g)|, there exist C, γ > 0 such that Pt (x, •) -λ V ≤ C(1 + V (x))e -γt .

This implies the result, since || • ||

T V ≤ || • || V .
Now the proof of Theorem 1.6 is straightforward.

Proof of Theorem 1.6. The rst part of the theorem has been proved in Proposition 3.1. For the second part, for any g ∈ C 3 b (R),

P N t g(x) -λ(g) ≤ P N t g(x) -Pt g(x) + Pt g(x) -λ(g) ≤ K t √ N (1 + x 2 )||g|| 3,∞ + ||g|| ∞ || Pt (x, •) -λ|| T V ≤||g|| 3,∞ C K t √ N (1 + x 2 ) + e -γt (1 + x 2 ) q ,
where we have used Theorem 1.4 and Proposition 3.1. Since (1 + x 2 ) q ≤ 1 + x 2 , q being smaller than 1/2, this implies the result.

4. Proof of Theorem 1.7

We prove the result using Theorem IX.4.15 of [START_REF] Jacod | Limit Theorems for Stochastic Processes[END_REF].

Let k ∈ N * , let us note Y N := (X N , Z N,1 , . . . , Z N,k ) and Ȳ := ( X, Z1 , . . . , Zk ). Using the notation of Theorem IX.4.15 of [START_REF] Jacod | Limit Theorems for Stochastic Processes[END_REF] with the semimartingales Y N (N ∈ N * ) and Ȳ and denoting e j (0 ≤ j ≤ k) the j-th unit vector, we have:

• b N,0 (x) = b 0 (x) = -αx and b N,i (x) = b i (x) = 0 for 1 ≤ i ≤ k, • c N,0,0 (x) = c 0,0 (x) = σ 2 f (x 0 ) and c N,i,j (x) = c i,j (x) = 0 for (i, j) = (0, 0), • g * K N (x) = f (x 0 ) k j=1 R g( u √ N e 0 + e j )dµ(u) + (N -k) R g( u √ N e 0 )dµ(u), • g * K(x) = f (x 0 ) k j=1 g(e j ).
The only condition of Theorem IX.4.15 that is not straightforward is the convergence of g * K N to g * K for g ∈ C 1 (R k+1 ). The complete denition of C 1 (R k+1 ) is given in VII.2.7 of [START_REF] Jacod | Limit Theorems for Stochastic Processes[END_REF], but here, we just use the fact that C 1 (R k+1 ) is a subspace of C b (R k+1 ) containing functions which are zero around zero. This convergence follows from the fact that any g ∈ C 1 (R k+1 ) can be written as g(x) = h(x)1 {|x|>ε} where h ∈ C b (R k+1 ) and ε > 0. This allows to show that, for this kind of function g,

(N -k)f (x 0 ) R g u √ N e 0 dµ(u) ≤ (N -k)f (x 0 )||h|| ∞ R 1 {|u|>ε √ N } dµ(u) ≤ f (x 0 )C N -K N 2 ≤ Cf (x 0 )N -1 ,
where the second inequality follows from the fact that we assume that µ is a probability measure having a nite fourth moment. Theorem IX.4.15 of [START_REF] Jacod | Limit Theorems for Stochastic Processes[END_REF] implies that for all k ≥ 1, (Z N,1 , ..., Z N,k ) converges to ( Z1 , ..., Zk ) in distribution in D(R + , R k ).

This implies the weaker convergence in D(R + , R) k for any k ∈ N * . Then, the convergence in D(R + , R) N * is classical (see e.g. Theorem 3.29 of [START_REF] Kallenberg | Foundations of Modern Probability[END_REF]).

Appendix

Extended generators

There are dierent denitions of innitesimal generators in the literature. The aim of this subsection is to dene precisely the notion of generator we use in this paper. Moreover we establish some properties of these generators and prove formula (10). In the general theory of semigroups, one denes the generators on some Banach space. In the frame of semigroups related to Markov processes, one generally considers (C b (R), || • || ∞ ). In this context, the generator A of a semigroup (P t ) t is dened on the set of functions D

(A) = {g ∈ C b (R) : ∃h ∈ C b (R), || 1 t (P t g -g) -h|| ∞ -→ 0 as t → 0}.
Then one denotes the previous function h as Ag. In general, we can only guarantee that D(A) contains the functions that have a compact support, but to prove Proposition 5.6, we need to apply the generators of the processes (X N t ) t and ( Xt ) t to functions of the type Ps g, and we cannot guarantee that Ps g has compact support even if we assume g to be in C ∞ c (R). This is why we consider extended generators (see for instance [START_REF] Meyn | Stability of Markovian Processes III: Foster-Lyapunov Criteria for Continuous-Time Processes[END_REF] or [START_REF] Davis | Markov Models and Optimization[END_REF]). These extended generators are dened by the point-wise convergence on R instead of the uniform convergence. Moreover, they verify the fundamental martingale property, which allows us to dene the generator on C n b (R) for suitable n ∈ N * and to prove that some properties of the classical theory of semigroups still hold for this larger class of functions.

Let (X t ) t be a Markov process taking values in R. We set

D(P ) = {g : R → R, measurable, s.t. ∀x ∈ R, ∀t ≥ 0, E x |g(X t )| < ∞}.
For g ∈ D(P ), x ∈ R, t ≥ 0, we dene

P t g(x) = E x [g(X t )] .
Denition 5.1. We dene D (A) to be the set of g ∈ D(P ) for which there exists a measurable function Ag : R → R, such that Ag ∈ D(P ), t → P t Ag(x) is continuous in 0, and ∀x ∈ R, ∀t ≥ 0,

(i) E x g(X t ) -g(x) = E x t 0 Ag(X s )ds; (ii) E x t 0 |A(g(X s ))|ds < ∞.
Remark 5.2. Using Fubini's theorem and (ii) we can rewrite (i) in the following form:

P t g(x) -g(x) = t 0 P s Ag(x)ds. (14)
Then ( 14) implies immediately that if g ∈ D (A), then

lim t→0 1 t (P t g(x) -g(x)) = Ag(x). ( 15
)
Note also that it follows from the Markov property and the denition of Ag that the process g(X t )g(X 0 ) -t 0 Ag(X s )ds is a P x -martingale w.r.t. to the ltration generated by (X t ) t .

The following result is classical and stated without proof. It is a straightforward consequence of ( 14) and ( 15). Proposition 5.3. Suppose that A is the extended generator of the semigroup (P t ) t , g ∈ D (A), and the map s → P s Ag(x) is continuous on R + for some x ∈ R. Then

d dt P t g(x) = P t Ag(x).
Moreover, if for all t ≥ 0, P t g ∈ D (A), then d dt P t g(x) = AP t g(x) = P t Ag(x). In what follows, we give some sucient conditions to verify the continuity and the derivability of the map s → P s h(x). These conditions are not intended to be optimal, they are stated such that it is easy to check them both for X N and X.

Proposition 5.4. Let (X t ) t be a Markov process with semigroup (P t ) t and extended generator A.

1

. Let h ∈ D(P ), x ∈ R. Suppose that (i) the map t → X t is continuous in L 2 , i.e. lim |t-s|→0 E x |X s -X t | 2 = 0; (ii) for all T > 0, sup 0≤t≤T E x (|X t | 4 ) < +∞; (iii) there exists C > 0, such that ∀x, y ∈ R, |h(x) -h(y)| ≤ C(1 + x 2 + y 2 )|x -y|.
Then the map s → P s h(x) is continuous on R + . 2. Suppose moreover that (i), (ii) and (iii) are satised with (iii)' g ∈ D (A) such that for some C > 0, and for all x, y ∈ R, we have that |Ag(x)-Ag(y)| ≤

C(1 + x 2 + y 2 )|x -y|.
Then the map s → P s g(x) is dierentiable on R + , and d dt P t g(x) = P t Ag(x).

Proof. The proof of point 1. follows from the following chain of inequalities

|P t h(x) -P s h(x)| ≤ E x |h(X t ) -h(X s )| ≤ CE x (1 + X 2 t + X 2 s )|X t -X s | ≤ C[E x (1 + X 4 t + X 4 s )] 1/2 [E x |X t -X s | 2 ] 1/2 ≤ C sup u≤s∨t [E x X 4 u ] 1/2 X s -X t 2 -→ 0 |t-s|→0
.

The second assertion of the proof follows from point 1. and Proposition 5.3, observing that h := Ag satises point (iii).

Proof of (10)

In this section, we rst collect some useful results about the extended generators A N of X N and Ā of X. Then we give the proof of (10). We start with the following result.

Proposition 5.5. 1. For all g ∈ C 3 b (R), for all x, y ∈ R,

| Āg(x) -Āg(y)| ≤ C g 3,∞ (1 + x 2 + y 2 )|x -y| and | Āg(x)| ≤ C||g|| 2,∞ (1 + x 2 ).
In particular, for any g ∈ C 3 b , the map t → Pt g(x) is dierentiable on R + , and d dt Pt g(x) = Pt Āg(x) = Ā Pt g(x).

For all

g ∈ C 2 b (R), for all x, y ∈ R, |A N g(x) -A N g(y)| ≤ C g 2,∞ (1 + x 2 + y 2 )|x -y| and |A N g(x)| ≤ C||g|| 1,∞ (1 + x 2 ).
In particular, for any g ∈ C 2 b , the map t → P N t g(x) is dierentiable on R + , and

d dt P N t g(x) = P N t A N g(x).
Proof. The result follows from Proposition 5.4 together with Lemma 2.1 and Lemma 2.2. Finally, to show that Pt Āg(x) = Ā Pt g(x), we use Proposition 5.3 and Proposition 2.4.

We are now able to give the proof of the main result of this section. This result is a Trotter-Kato like formula that allows to obtain a control of the dierence between the semigroups P and P N , provided we dispose already of a control of the distance between their generators. It is an adaptation of Lemma 1.6.2 from [START_REF] Ethier | Markov Processes. Characterization and Convergence[END_REF] to the notion of extended generators.

Proposition 5.6. Grant Assumptions 1, 2 and 3. Let Ā and A N be the extended generators of respectively P and P N .

Then the following equality holds for each

g ∈ C 3 b (R), x ∈ R and t ∈ R + . Pt -P N t g(x) = t 0 P N t-s Ā -A N Ps g(x)ds. (16)
Proof. We x t ≥ 0,

N ∈ N * , g ∈ C 3 b (R), x ∈ R in the rest of the proof. Introduce for 0 ≤ s ≤ t the function u(s) = P N t-s Ps g(x). One can note that u = Φ • Ψ with Φ : R 2 → R; Φ(v 1 , v 2 ) = P N v1 Pv2 g(x)
and Ψ : R → R 2 ; Ψ(s) = (t -s, s). Let us show that Φ is dierentiable w.r.t. to both variables v 1 and v 2 . Indeed, for v 1 it is a consequence of the fact that h = Pv2 g ∈ C 3 b (R) by Proposition 2.4 and Proposition 5.5, from which we know that if

h ∈ C 2 b , then v 1 → P N v1 h(x) is dierentiable and ∂ dv 1 Φ(v 1 , v 2 ) = d dv 1 P N v1 h(x) = P N v1 A N h(x).
To show the dierentiability of Φ with respect to v 2 , let us write

Φ(v 1 , v 2 ) = E x Pv2 g(X N v1 ) . From Proposition 5.5, we know that since g ∈ C 3 b , v 2 → Pv2 g(X N v1 ) is a.s. dierentiable with derivative d dv 2 Pv2 g(X N v1 ) = Ā Pv2 g(X N v1 ) = Pv2 Āg(X N v1 ) = E X N v 1 ( Āg)( Xv2 ). Moreover, | Āg(x)| ≤ C||g|| 2,∞ (1 + x 2
) by Proposition 5.5. Now, using Lemma 2.1.(ii) we see that

sup v2≤T d dv 2 Pv2 g(X N v1 ) ≤ E X N v 1 sup v2≤T |( Āg)( Xv2 )| ≤ C T 1 + (X N v1 ) 2 .
By Lemma 2.1.(iii), we see that the last bound is integrable, hence by dominated convergence,

v 2 → Φ(v 1 , v 2 ) is dierentiable with derivative ∂ dv 2 Φ(v 1 , v 2 ) = P N v1 Ā Pv2 g(x) = P N v1 Pv2 Āg(x).
As a consequence, u is dierentiable on R + , and we have

u (s) = - ∂ ∂v 1 Φ(t -s, s) + ∂ ∂v 2 Φ(t -s, s) = -P N t-s A N Ps g(x) + P N t-s Ps Āg(x) =P N t-s Ā -A N Ps g(x).
Now we show that u is continuous. Indeed, if it is the case, then we will have

u(t) -u(0) = t 0 u (s)ds,
which is exactly the assertion.

In order to prove the continuity of u , we consider a sequence (s k ) k that converges to some s ∈ [0, t], and we write

P N t-s Ā -A N Ps g(x) -P N t-s k Ā -A N Ps k g(x) ≤ P N t-s -P N t-s k Ā -A N g s (x) (17)
+ P N t-s k Ā -A N Ps -Ps k g(x) , (18)
where g s = Ps g ∈ C 3 b (R). To show that the term (17) vanishes when k goes to innity, denote h s (x) = ( Ā -A N )g s (x). Using Proposition 5.5 and the fact that g s ∈ C 3 b (R), we have

|h s (x) -h s (y)| ≤ C(1 + x 2 + y 2 )|x -y|.
Proposition 5.4 applied to h s and to P N implies that u → P N u h s (x) is continuous. As a consequence the term (17) vanishes as k → ∞.

To nish the proof, we need to show that the term (18) vanishes. Denote g k = Ps -Ps k g. We have to show that

E x Ā -A N g k (X N t-s k ) → 0, when k → ∞.
In what follows we will in fact show that

E x Āg k (X N t-s k ) → 0 and E x A N g k (X N t-s k ) → 0, when k → ∞. (19)
To begin with, using Proposition 2.4, the functions g k belong to C 3 b (R), and for any i ∈ {0, 1, 2}, for all y ∈ R, g

k (y) vanishes as k goes to innity. Using again Proposition 2.4, we see that for all i ∈ {0, 1, 2, 3}, ||g

(i) k || ∞ is uniformly bounded in k. It follows that each sequence (g (i) k ) k , i ∈ {0, 1, 2},
is uniformly equicontinuous and thus converges to zero uniformly on each compact interval.

We next show that this implies that also the sequences (A N g k ) k and ( Āg k ) k converge to zero uniformly on each compact interval. For ( Āg k ) k , this is immediate, since Ā is a local operator having continuous coecients. For (A N g k ) k , it follows from the fact that A N g k (x) → 0 as k → ∞ for each xed x and the fact that by Lemma 5.7 given below, this sequence is uniformly (in k, for xed N ) equicontinuous on each compact.

We are now able to conclude. The sequence (X N t-s k ) k is almost surely bounded by sup 0≤r≤t |X N r | which is nite almost surely by Lemma 2.1.(iii). Hence, almost surely as k → ∞, Āg k (X N t-s k ) → 0 and A N g k (X N t-s k ) → 0. We now apply dominated convergence to prove (19). Using that by Proposition 5.5, for all

g ∈ C 3 b (R) and x ∈ R, | Āg(x)| ≤ C||g|| 2,∞ (1 + x 2 ),
we can bound the expression in the rst expectation by

C||g k || 2,∞ 1 + ( sup 0≤r≤t |X N r |) 2 ≤ 2C sup 0≤r≤t || Pr g|| 2,∞ 1 + ( sup 0≤r≤t |X N r |) 2 ,
whose expectation is nite thanks to Lemma 2.1(iii). The same arguments work for A N . This implies that (18) vanishes as k → ∞, and this concludes the proof.

We now prove the missing lemma Lemma 5.7. For all g ∈ C 2 b (R) and any M > 0,

sup x∈[-M,M ] | A N g (x)| ≤ C N g 2,∞ 1 + M 2 ,
for some constant C N > 0 that can depend on N, but not on M.

Proof. We have

A N g (x) = -αg (x) -αxg (x) + N f (x)E g x + U √ N -g(x) + N f (x)E g x + U √ N -g (x)
.

Since E g x + U √ N -g(x) ≤ g ∞ √ N E [|U |] , E |g x + U √ N -g (x)| ≤ g ∞ √ N E [|U |] ,
we obtain

sup x∈[-M,M ] | A N g (x)| ≤ |α| g 2,∞ (1 + M ) + √ N (|f (x)| ∨ |f (x)|) g 2,∞ E [|U |] . Assumption 3 implies that |f (x)| ≤ m 1 C √ 1 + x 2 for all x.
Together with the sub-quadraticity of f, this concludes the proof. 5.3. Existence and uniqueness of the process X N t t

Proposition 5.8. If Assumptions 1 and 2 hold, the equation ( 5) admits a unique non-exploding strong solution.

Proof. It is well known that if f is bounded, there is a unique strong solution of (5) (see Theorem IV.9.1 of [START_REF] Ikeda | Stochastic Dierential Equations and Diusion Processes[END_REF]). In the general case we reason in a similar way as in the proof of Proposition 2 in [START_REF] Fournier | On a toy model of interacting neurons[END_REF]. Consider the solution (X N,K t ) t∈R+ of the equation ( 5) where f is replaced by f K : x ∈ R → f (x) ∧ sup |y|≤K f (y) for some K ∈ N * . Introduce moreover the stopping time

τ N K = inf t ≥ 0 : X N,K t ≥ K . Since for all t ∈ 0, τ N K ∧ τ N K+1 , X N,K t = X N,K+1 t
, we know that τ N K (ω) ≤ τ N K+1 (ω) for all ω. Then we can dene τ N as the non-decreasing limit of τ N K . With a classical reasoning relying on Itô's formula and Grönwall's lemma, we can prove that

sup 0≤s≤t E X N,K s∧τ N K 2 ≤ C t 1 + x 2 , ( 20
)
where C t > 0 does not depend on K. As a consequence, we know that almost surely, τ N = +∞. So we can simply dene X N t as the limit of X N,K t , as K goes to innity. Now we show that X N satises equation ( 5). Consider some ω ∈ Ω and t > 0, and choose K such that τ N K (ω) > t. Then we know that for all s ∈ [0

, t], X N s (ω) = X N,K s (ω) and f (X N s-(ω)) = f K (X N,K s-(ω)).
Moreover, as X N,K (ω) satises equation ( 5) with f replaced by f K , we know that X N (ω) veries equation ( 5) on [0, t]. This holds for all t > 0. As a consequence, we know that X N satises equation ( 5). This proves the existence of a strong solution. The uniqueness is a consequence of the uniqueness of strong solutions of (5), if we replace f by f K in (5), and of the fact that any strong solution (Y N t ) t equals necessarily (X N,K t ) t on [0, τ N K ].

5.4. Proof of Lemma 2.1

Proof. We begin with the proof of (i). Let Φ(x) = x 2 and A N be the extended generator of (X N t ) t≥0 . One can note that, applying Fatou's lemma to the inequality (20), one obtains for all t ≥ 0, sup 0≤s≤t E (X N s ) 2 is nite. As a consequence Φ ∈ D (A N ) (in the sense of Denition 5.1). And, recalling that µ is centered and that σ 2 := R u 2 dµ(u), we have for all x ∈ R,

A N Φ(x) =-αxΦ (x) + N f (x) R Φ(x + u √ N ) -Φ(x) dµ(u) =-2αx 2 + N f (x) R 2x u √ N + u 2 N dµ(u) = -2αΦ(x) + σ 2 f (x) ≤ -2αΦ(x) + σ 2 L|x| + f (0) 2 ≤(σ 2 L 2 -2α)Φ(x) + 2σ 2 L|x| f (0) + σ 2 f (0).
Let ε > 0 be xed, and

η ε = 2σ 2 L f (0)/ε. Using that, for every x ∈ R, |x| ≤ x 2 /η ε + η ε , we have A N Φ(x) ≤ c ε Φ(x) + d ε , (21)
with

c ε = σ 2 L 2 -2α + ε and d ε = O(1/ε). Let us assume that c ε = 0, possibly by reducing ε > 0. Considering Y N t := e -cεt Φ(X N t ), by Itô's formula, dY N t = -c ε e -cεt Φ(X N t )dt + e -cεt dΦ(X N t ) = -c ε e -cεt Φ(X N t )dt + e -cεt A N Φ(X N t )dt + e -cεt dM t ,
where, denoting by πj (dt, dx, du) := π j (dt, dx, du) -dtdxdµ(u) the compensated measure of π j (1 ≤ j ≤ N ), (M t ) t≥0 is the P x -local martingale dened as

M t = N j=1 [0,t]×R+×R Φ(X N s-+ u √ N ) -Φ(x) 1 {z≤f(X N s-)} dπ j (s, z, u).
One can note that, since sup 0≤s≤t E (X N s ) 2 is nite for any t ≥ 0, (M t ) t≥0 is a locally square integrable P x -local martingale, and as a consequence, it is a P x -martingale.

Using (21), we obtain

dY N t ≤ d ε e -cεt dt + e -cεt dM t , implying E x Y N t ≤ E x Y N 0 + d ε c ε e -cεt + 1 . One deduces E x X N t 2 ≤ x 2 e (σ 2 L 2 -2α+ε)t + C ε e (σ 2 L 2 -2α+ε)t + 1 , (22)
for some constant C > 0 independent of t, ε, N.

The proof of (ii) is analogous and therefore omitted.

Consequently, the stopping times τ N K tend to innity as K goes to innity, and Fatou's lemma allows to conclude.

We nally prove (v). Indeed, by Itô's isometry and Jensen's inequality, for all 0 ≤ s ≤ t ≤ T, using the sub-quadraticity of f and (i),

E x (X N t -X N s ) 2 =E x      -α t s X N r dr + 1 √ N N j=1 ]s,t]×R+×R u1 {z≤f(X N r-)} dπ j (r, z, u)   2    ≤2α 2 (t -s) t s E x (X N r ) 2 dr + 2σ 2 t s E x f (X N r ) dr ≤2α 2 C t (1 + x 2 )(t -s) 2 + 2σ 2 C t (1 + x 2 )(t -s) ≤C T (t -s)(1 + x 2 ).
This proves that X N satises hypothesis (v). A similar computation holds true for X.

Proof of Proposition 2.4

Proof. To begin with, we use Theorem 1.4.1 of [START_REF] Kunita | Lectures on Stochastic Flows and Applications for the Indian Institute Of Science Bangalore[END_REF] to prove that the ow associated to the SDE (6) admits a modication which is C 3 with respect to the initial condition x (see also Theorem 4.6.5 of [START_REF] Kunita | Stochastic ows and stochastic dierential equations[END_REF]). Indeed the local characteristics of the ow are given by b(x, t) = -αx and a(x, y, t) = σ 2 f (x)f (y), and, under Assumptions 1 and 3, they satisfy the conditions of Theorem 1.4.1 of [START_REF] Kunita | Lectures on Stochastic Flows and Applications for the Indian Institute Of Science Bangalore[END_REF]:

• ∃C, ∀x, y, t, |b(x, t)| ≤ C(1 + |x|) and |a(x, y, t)| ≤ C(1 + |x|)(1 + |y|). • ∃C, ∀x, y, t, |b(x, t) -b(y, t)| ≤ C|x -y| and |a(x, x, t) + a(y, y, t) -2a(x, y, t)| ≤ C|x -y| 2 . • ∀1 ≤ k ≤ 4, 1 ≤ l ≤ 4 -k, ∂ k ∂x k b(x, t
) and ∂ k+l ∂x k ∂y l a(x, y, t) are bounded.

In the following, we consider the process ( X(x) t ) t , solution of the SDE (6) and satisfying X(x) 0 = x. Then we can consider a modication of the ow X(x) t which is C 3 with the respect to the initial condition x = X(x) 0 . It is then sucient to control the moment of the derivatives of X(x) t with respect to x, since with those controls we will have

Pt g(x) =E g X(x) t , Pt g (x) = E ∂ X(x) t ∂x g X(x) t , Pt g (x) =E   ∂ 2 X(x) t ∂x 2 g X(x) t + ∂ X(x) t ∂x 2 g X(x) t   , Pt g (x) =E   ∂ 3 X(x) t ∂x 3 g X(x) t + 3 ∂ 2 X(x) t ∂x 2 • ∂ X(x) t ∂x g X(x) t + ∂ X(x) t ∂x 3 g X(x) t   . ( 25
)
We start with the representation

X(x) t = xe -αt + σ t 0 e -α(t-s) f X(x) s dB s . This implies ∂ X(x) t ∂x = e -αt + σ t 0 e -α(t-s) ∂ X(x) s ∂x f X(x) s dB s . (26)
Writing U t = e αt ∂ X(x) t ∂x and

M t = t 0 σ f X(x) s dB s , (27)
we obtain U t = 1 + t 0 U s dM s , whence

U t = exp M t - 1 2 < M > t . (28)
Notice that this implies U t > 0 almost surely, whence

where m 1 is the bound of ( √ f ) introduced in Assumption 3. In particular we have

We introduce V t = ∂ 2 X(x) t ∂x 2 e αt and deduce from this that

V t =σ t 0 V s f X(x) s + e -αs U 2 s f (2)
X(x) s dB s , which can be rewritten as tn ∂x g ( X(x) tn ), n ≥ 1 is uniformly integrable. As a consequence, the second formula in (25) implies that ( Ptn g) (x) → ( Pt g) (x) as n → ∞, whence the desired continuity. The argument is similar for the second derivative, using (31) and (33). That concludes the proof.

dV t = V t dM t + Y t dB t , V 0 = 0, Y t = σe -αt U 2 t f (

  as in(27). Applying Itô's formula to Z t := V t /U t (recall that U t > 0), we obtain such that, by the precise form of Y t and since Z 0 = 0, Using Jensen's inequality, (29) and Burkholder-Davis-Gundy inequality, for all t ≥ 0, + t3 1 + t + e (6σ 2 m 2 1 -4α)t ≤ C(t + t 4 )e (6σ 2 m 2 1 -4α)t . (32) ≤ C(t 1/2 + t 2 )e 3σ 2 m 2 1 -2αt e 3σ 2 m 2 1 t ≤ C(t 1/2 + t 2 )e 6σ 2 m 2 1 -2αt , Introducing W t = e αt ∂ 3 X(x)Once again we can rewrite this asdW t = W t dM t + Y t dB t , W 0 = 0,whereY t = σ 3e -αt U t V t fwhich proves the rst assertion of the proposition. The proof of the second assertion, equation (12), follows similarly. Finally to prove the third assertion, we rst study the regularity of the rst derivative. Notice that t → surely continuous by equation (26). Now take any sequence t n → t. By (30), the family of random variables

	Z t = σ		t	e -αs U s		f	(2)		X(x) s		dB s -σ 2 (2) X(x) t	t	e -αs U s + e -2αt U 3 t	f	(2) f	X(x) s (3) X(x) t	f	,	X(x)
	0 whence, introducing Z t = Wt Ut ,								0
	E Z 4 t ≤ C E		t	e -αs U s Z t =	f 0	t	(2) Y s U s	X(x) s dB s -	dB s t 0 Y s U s	4 d < M, B > s .
	0 As previously, we obtain,
	+E 2 ≤C(1 + t) t 0 e -αs U s ≤ C E E (Z t ) t 0 E t ≤C(1 + t) e -2αs E V 2 f t e -2αs U 2 X(x) s s 2 Y s ds U s s + e -4αs E U 4 f (2) f (2) X(x) s s 0 +E t e -αs U s f X(x) s f (2) 0 ≤C(1 + t) t (s 1/2 + s 2 )e (6σ 2 m 2 1 -4α)s + e (6σ 2 m 2 X(x) s 2 ds ds X(x) s 1 -4α)s ds 4 ds 2 ds 4
								0	0
					≤ C t + t 3 ≤C(1 + t 3 )	0	t 0 e -4αs E U 4 s ds ≤ C t + t 3 t e (6σ 2 m 2 1 -4α)s dss	0	t	e (6σ 2 m 2 1 -4α)s ds
	≤ C t We deduce that
	E V 2 t 1/2 whence ≤ E Z 4 t 1/2 E U 4 t E  ≤C(1 + t 2 ) e  ∂ 2 X(x) 1 2 σ 2 m 2 1 t + e ( 7 2 2 σ 2 m 2 1 -2α)t ,  t  ≤ C(t 1/2 + t 2 )e (6σ 2 m 2 1 -4α)t . implying ∂x 2 Finally, dierentiating (31), we get E ∂ 3 X(x) t ∂ 3 x ≤ C(1 + t 2 ) e ( 1 2 σ 2 m 2 1 -α)t + e ( 7 2 σ 2 m 2 1 -3α)t .	1 2 σ 2 m 2 1 t	(33) (36)
	∂ 3 X(x) Finally, using Cauchy-Schwarz inequality, and inserting (30), (33) and (36) in (25), s ∂x 3 f X(x) s + 3 ∂ 2 X(x) s ∂x 2 ∂ X(x) s f (2) X(x) s ∂x + ∂ X(x) s ∂x 3 f X(x) (3) Pt g ∞ ≤ C||g|| 3,∞ (1 + t 2 ) e ( 1 2 σ 2 m 2 1 -α)t + e 2(σ 2 m 2 1 -α)t + e ( 7 2 σ 2 m 2 1 -3α)t ,	(34)
											∂	X(x)
	W t = σ	t	W s	f	X(x) s ∂x 3 , we obtain + 3e -αs U s V s t	f	(2)	X(x) s	2) + e -2αs U 3 s X(x) ∂ X(x)	f	(3)	X(x)
	0
										dZ t =	Y t U t	dB t -	Y t U t	d < M, B > t ,

t , with M t s ds.

t ∂x 3 = σ t 0 e -α(t-s) ∂ 3 X(x) s   dB s . s dB s . ≤C(1 + t 3 )(1 + t + e (6σ 2 m 2 1 -4α)t ) ≤ C(1 + t 4 ) 1 + e (6σ 2 m 2 1 -4α)t .

(35)

As a consequence,

E [|W t |] ≤E (Z t ) 2 1/2 E U 2 t 1/2 ≤ C(1 + t 2 ) 1 + e (3σ 2 m 2 1 -2α)t e t ∂x

is almost

Now we prove (iii). From

u1 {z≤f(X N s-)} dπ j (s, z, u),

Applying Burkholder-Davis-Gundy inequality to the last term above in (23), we can bound its expectation by

Now, bounding the expectation of ( 23) by ( 24), and using point (i) of the lemma we conclude the proof of (iii).

The assertion (iv) can be proved in classical way, applying Itô's formula and Grönwall's lemma. Let us explain how to prove this property for the process X N . The proof for X is similar. According to Itô's formula, for every t ≥ 0,

Let us recall that u is centered and has a nite fourth moment, and that f is subquadratic. Introducing the stopping times τ N K := inf{t > 0 :

) 4 , it follows from the above that for all t ≥ 0,

where C is a constant independent of t, N and K. This implies that for all t ≥ 0, sup


F Silva
email: silva@lma.cnrs-mrs.fr

T Hélie

V Wetzel

Port-Hamiltonian Representation of Dynamical Systems. Application to Self-Sustained Oscillations in the Vocal Apparatus

Phonation is a natural example of nonlinear dynamical system with self-sustained oscillations, here resulting from the controlled nonlinear coupling between the deformable vocal folds and the airflow expired from the lungs through the glottis and the vocal tract. As a proof of concept, we propose a minimal model of the full vocal apparatus using the port-Hamiltonian representation that emphasises on the structure of a system (on the separation between the behaviour of the subsystems and their interconnection) and on the power exchanges. Numerical results on bifurcations are qualitatively discussed in relation to voice pathologies.

Motivations The physics of voice production has motivated a wide variety of models, from full-featured numerical ones (mainly based on FEM for vocal folds and FVM for the airflow) to reduced order models focusing on the essential phenomenon underlying the phonation. A large body of work in the latter category relies on the description of the glottal aerodynamics from the late 1950 that is based on experiments in rigid static larynx-like ducts. However phonation intrinsically implies the vibration of the vocal folds periodically closing the glottis. There is a significant paradox in those simplified models: the vocal folds move as a consequence of the power exchanged with the glottal flow, but the description of the latter assumes that it does not receive or provide any power to the folds. This power imbalance in the models introduces a bias in the analysis of the coupling occurring in the larynx and of the instability leading to voice production.

Port-Hamiltonian systems (PHS)

The port-Hamiltonian theory combines the views of the (geometric) Hamiltonian mechanics and of port-based modelling approach, emphasising the separation between the behaviour of components and their interconnection. The lingua franca of this theory is energy : port-Hamiltonian systems are open passive systems that can store, dissipate and exchange power with their neighbourhood. Energy-storage is defined by the Hamiltonian H(x) as a function of state variables x. This dependency, formulated as the Hamiltonian gradient ∇ x H, expresses the effort of the energy-storing component. Conversely, the evolution of this component is described by the flux variable ẋ. Effort and flux are power dual variables, i.e., they jointly define the energy flow dH/dt = ẋ • ∇H(x).

Dissipating components are described by dissipation variable w and their constitutive laws z(w), such that the dissipated power is P diss = w • z(w) ≥ 0. Finally, external interactions are classically described in terms of inputs u (also called efforts) and output y (fluxes) that are dual with respect to the external power: P ext = y • u (≥ 0 when the system yields power).

The geometric structure (Dirac structure, see [START_REF] Van Der Schaft | Port-Hamiltonian Systems Theory: an Introductory Overview[END_REF]) accounting for the interconnection of the subsystems is then described by a matrix S relating efforts to fluxes:

  ẋ w y   = S(x, w)   ∇H z(w) u   . (1)
When the interconnection is conservative, the power balance is ensured provided that the matrix S is skew-symmetric (e • S • e = 0 for any effort e).

Minimal port-Hamiltonian model of the vocal apparatus

The vocal apparatus is modelled as the interconnection of vocal folds, glottal flow, vocal tract and a subglottal pressure supply. Every component is designed minimally with an emphasis on the dual pairing on the interconnection ports. For instance, each vocal fold is considered to be a single-d.o.f. oscillator with an elastic cover, and is submitted to pressure on the upstream and downstream faces, and interacts with the flow at the glottal face. Reciprocally, the glottal flow is modelled as the simplest kinematics of potential incompressible flow of inviscid air preserving the continuity of the normal velocity on the surface of the vocal folds (see details in Ref. [START_REF] Helie | Self-oscillations of a Vocal Apparatus: A Port-Hamiltonian Formulation[END_REF]). Downstream the glottis, the flow separates from the folds into a jet that spreads and dissipates its kinetic energy into heat. Finally, the vocal tract is represented as seen by the larynx, i.e. by means of its input impedance describing the acoustic feedback on the larynx.

The full vocal apparatus is obtained by conservative interconnection of the previous subsystems. For the sake of conciseness, the resulting Hamiltonian and matrix S are not reported here but can be efficiently computed using the Python package PyPHS [START_REF] Falaize | PyPHS: Passive modeling and simulation in python[END_REF] designed for the symbolic manipulation of PHS systems. Numerical results Time-domain simulations are performed using the numerical scheme designed after the principles of the port-Hamiltonian theory [START_REF] Lopes | Approche passive pour la modélisation, la simulation et l'étude d'un banc de test robotisé pour les instruments de type cuivre[END_REF], i.e. preserving power balance, but also being more consistent than standard integrators (e.g., RK45, . . . ). Time-domain simulations and bifurcation analysys evidence the ability of the model to produce oscillating regimes above some pressure threshold, and quasi-periodic vibrations to possible intermittency when vocal folds are detuned.

P sub Q sub P sub r Q sub r P sub l Q sub l P - tot Q - F p l v l F p r v r P + Q + P sup l Q sup l P sup r Q sup r P ac Q ac

Figure 1 :

 1 Figure 1: Components of the vocal apparatus. The interconnection takes place via pairs of effort (P ) and flux (Q) variables.


Aude Lagier
email: aude.lagier@chuliege.be

Thierry Legou
email: thierry.legou@lpl-aix.fr

Fabrice Silva
email: silva@lma.cnrs-mrs.fr

Thomas Hélie
email: thomas.helie@ircam.fr

INVESTIGATING PHONATION THROUGH EXCISED HUMAN LARYNGES: RECENT DEVELOPMENTS AND ONGOING WORKS ON AN ANIMATED TESTBENCH

Many aspects of voice production can be investigated by means of in vivo observations and experiments on synthetic vocal folds replica. However, either mammals or human excised larynges provide an interesting trade-off between realistic physiological conditions, reproducibility and observability, paving the way to the development in the last years of experimental setups involving excised human larynges. Some of them focus on the control of the laryngeal articulation, in order to enable the simulation of coordinated muscular activations.

We here present the last developments of the excised human larynx test bench hosted at the Medecine Faculty in Marseille, with the motorized animation of some of the intrinsic laryngeal muscles. As an evidence of the possibilities of such experimental devices, we present a quantitative analysis of the contact pressure between vocal folds in presence of a Reinke's edema.

INTRODUCTION

From a physical point of view, phonation is seen as the sound production induced by the self-sustained oscillation of the vocal folds in interaction with the airflow in the glottis (see Ref. [START_REF] Titze | The Myoelastic Aerodynamic Theory of Phonation[END_REF] for a detailed description). This is enabled by the combination of convenient geometry (that speeds up the airflow in the glottal constriction) and elasticity of the tissues (that counteracts the resultant of the air pressure, with possible delay).

In vivo studies are made challenging by the difficult evaluation of the geometric, biomechanical and aerodynamical characteristics of the human larynx from videoendoscopic imaging, electroglottographic or even electromygraphic measurements through the vocal tract of the subject or his skin. This is even more intricate considering the uncontrolled muscles contraction that hinder repeatability. This paved the road to experiments on vocal folds replica: pioneering setup designed by Muller [START_REF] Muller | Elements of physiology[END_REF] was either simplified to investigate the peculiarities of the flow in idealized glottal channel (either static [START_REF] Berg | On the air resistance and the Bernoulli effect of the human larynx[END_REF] or motorized [START_REF] Triep | High-speed piv measurements of the flow downstream of a dynamic mechanical model of the human vocal folds[END_REF]) or rather augmented in order to improve the decoupled control of several geometrical and mechanical parameters (see, e.g., Ref. [START_REF] Ruty | An in vitro setup to test the relevance and the accuracy of low-order vocal folds models[END_REF]) involved in the self-sustained oscillations. They are usually made of an elastomer envelop containing water [START_REF] Ruty | An in vitro setup to test the relevance and the accuracy of low-order vocal folds models[END_REF] or a bioreactor [START_REF] Titze | Design and validation of a bioreactor for engineering vocal fold tissues under combined tensile and vibrational stresses[END_REF], or of a tissue-like hydrogel [START_REF] Yousefi-Mashouf | Multiaxial mechanical properties of hydrogel-based materials upon finite strains: towards the design of tailored vocal-fold composite replicas[END_REF].

Nevertheless, by now, none of these synthetic designs manages to reach to combine observability to realistic anatomical or physiological conditions. Either mammals or human excised larynges provide an interesting tradeoff between those requirements, and the last decades witnessed a large variety of experimental setups involving excised larynges (see Garcia and Herbst [8] and Luo et al [START_REF] Luo | Development of excised larynx[END_REF] for recent overviews). However, even these setups may be affected by oversimplified controls: as an example, the automated setup designed in Erlangen [START_REF] Birk | Automated setup for ex vivo larynx experiments[END_REF] suffers from an over-constrained mobility of the arytenoids cartilages as the chosen electromechanical actuators inhibits the essential rocking motion of the arytenoids on the cricoid facets.

The aim of this communication is thus to present the last developments of the excised human larynx test bench hosted at the Medecine Faculty in Marseille, with the motorized animation of some of the intrinsic laryngeal muscles (Sec. 2) designed to be physiologically realistic. As an evidence of the possibilities of such experimental devices, we present three studies: first the validation of the ex vivo laryngeal articulation (Sec. 3.1), then the control of the fundamental frequency f o using the tensioning of the crico-thyroid muscle (Sec. 3.2), and finally a quantitative analysis of the contact pressure between vocal folds in presence of a Reinke's edema (Sec. 3.3).

DESCRIPTION OF THE ANIMATED TESTBENCH 2.1 Larynx preparation

Human larynges are harvested from fresh cadavers, within 24 to 48 hours after death, in accordance with the French ethical laws. Larynges (including the upper three tracheal rings) are then rapidly frozen in 0.9% saline solution. Before experiment, the larynx is thawed in saline solution and cleaned with removal of the extraneous tissues and supraglottic structures (epiglottis and vestibular folds), within literature standards [START_REF] Alipour | Phonatory characteristics of the excised human larynx in comparison to other species[END_REF]. Ref. [START_REF] Hanna | Investigations of the acoustics of the vocal tract and vocal folds in vivo, ex vivo and in vitro[END_REF] evidences that dehydratation of the tissues during experiments seems to have more impact on phonation characteristics (such as the f o range) and on tissue damaging than successive freezing/thawing cycles.

The larynx is then secured on the experimental bench, the posterior edges of the thyroid cartilage lying horizontally upon adaptable steel workpieces (with hose ties on the superior horns). The gap (adjustable in order to deal with the interindividual variability in larynx dimensions) between the latter enables the cricoid cartilage and the trachea to move freely during the crico-thyroid tilt.

Actuation

The intrinsic laryngeal musculature being disabled on excised organs (except for preliminary ex vivo perfused larynx under neuromuscular electrical stimulation [START_REF] Mendelsohn | Preliminary study of the open quotient in an ex vivo perfused human larynx[END_REF]), it is necessary to set up substitutes enabling the reproduction of the effects of the muscle activation. It was the aim of the ANIMAGLOTTE project (funded by the Mission for Interdisciplinarity CNRS, 2018, coordinated by Th. Legou) not to impose the motions of the cartilages (as done, e.g., in Ref. [START_REF] Birk | Automated setup for ex vivo larynx experiments[END_REF] using needles pricked to the cartilages) but to mimic the forces resulting from muscular activation (see Ref. [START_REF] Titze | The Myoelastic Aerodynamic Theory of Phonation[END_REF] for a review of the muscles and their orientation). For most of the controls, this means that we focus on the approximation of the origin and the insertion of the essential intrinsic muscles of the larynx. The whole actuation setup is shown on Fig. 1 and is controlled using an Arduino board.

Crico-thyroid (CT) pars recta

This muscle links the anterior parts of the cricoid and the thyroid cartilages, is innervated by the external branch of the superior laryngeal nerve (contrarily to all other muscles, which are innervated by the recurrent laryngeal nerve) and is generally reported as the primary control for the fundamental frequency. Considering the thyroid as a reference frame, the contraction of the CT muscle induces the tilt of the cricoid that supports the arytenoids, which lengthens, thins and tenses the vocal folds. This is simulated on excised larynges by a pulling system: an inextensible nylon thread is stitched on the anterior arch of the cricoid, runs along the pars recta and exits through an orifice done into the lamina angle of the thyroid (see Fig. 2), and is attached on the mobile rod of a linear actuator (Firgelli L12-100, range 100 mm). The CT angle related to the tilt and the lengthening of the vocal folds using this pulling system have been reported in Ref. [START_REF] Lagier | Control of the glottal configuration in ex vivo human models: quantitative anatomy for clinical and experimental practices[END_REF], highlighting an intrinsic upper bound for this control corresponding to the situation where the anterior parts of the cricoid and the thyroid cartilages are in contact.

Thyroarytenoid muscle (TA)

These paired muscles run from the angle of the thyroid to the anterolateral surface of the arytenoids, they would shorten the vocal folds. However, being antagonist to the CT, the contraction of the muscular fibers essentially contributes to the thickening of the membranous vocal fold and thus to its medialization. As described in Ref. [START_REF] Lagier | Control of the glottal configuration in ex vivo human models: quantitative anatomy for clinical and experimental practices[END_REF], this effect is reproduced using a type 1 thyroplasty performed by means of Montgomery® implants inserted through bilateral openings done in the thyroid wings (see Fig. 3). The handle of each implant is attached to a linear motion actuator.

Figure 3. TA actuation: schematics of the medialization system (left, top view), and realization on an excised larynx (right, from Ref. [START_REF] Lagier | Control of the glottal configuration in ex vivo human models: quantitative anatomy for clinical and experimental practices[END_REF]).

Lateral cricoarytenoid muscles (LCA)

These paired muscles run obliquely from the anterolateral surfaces of the cricoid to the lateral tips of the arytenoids. Their activation contributes to the rotation of the arytenoids on their relative cricoid facets, the processes of the vocalis moving inwards which classifies the LCA as adductors muscles.

Following Refs. [START_REF] Hanna | Investigations of the acoustics of the vocal tract and vocal folds in vivo, ex vivo and in vitro[END_REF][START_REF] Berry | Bifurcations in excised larynx experiments[END_REF], sutures are done on the muscular processes of the arytenoids, the threads then course anteriorly, medially and inferiorly, passing below the thyroid wings to preserve the mean orientation towards the muscular origins on the anterolateral surfaces of the cricoid (see Fig. 4). Each thread is then attached to the wheel of a stepper motor 17HM15-0904S.

Posterior cricoarytenoid muscles (PCA)

These paired muscles attach to the posterior face of the cricoid and the muscular process of each arytenoid. Their activation contribute to rotate the latter (thus abducting the vocal folds) and to pull them posteriorly, acting then as direct antagonists to the LCA muscles. As for LCA actuation, threads are sewed to reproduce physiological force directions and are attached to wheels of servo (same as for LCA actuation, see Fig. 5).

Aerodynamics and Sensors

For experiments implying aerodynamics (needed for phonation), airflow may be supplied from a turbine (Werei Rietschle, Fürstenfeldbruck, Germany) through an electrovalve (Burkert 2875, now replacing the previous quarter-turn manual valve, see Ref. [START_REF] Legou | Test bench for human excised larynx studies[END_REF]), a flowmeter (TSI4043, full scale 200 L min -1 ), a heater/humidifier (Drager, Aquador, Dräger) and an intubation tube (Rüsch Elit, inner diameter 6.5 mm, with balloon cuff fitted in the tracheal rings), respectively. Finally, the subglottal aerodynamical pressure is measured by inserting a venous catheter (23-gauge) through the crico-thyroid membrane. The catheter is then connected to the piezo-resistive differential pressure sensor provided by the EVA workstation (SQLab, Aix-en-Provence, France).

Concerning acoustics, a microphone (AKG C520) may be placed 30 cm away for the larynx to record the radiated acoustic pressure. Additionally, a probe microphone (GRAS 40SC) can be inserted within the trachea along the balloon cuff in order to measure the acoustic component of the subglottal pressure. Vibrations of the vocal folds may also be monitored by electroglottographic electrodes placed on the anterior aspect of the thyroid laminae and connected to the Laryngograph microProcessor (Wallington, United Kingdom), and/or by a high-speed camera (Fastcam SA3, Photron, external trigger, with a micro-lens Nikon AF MicroNikkor 105 mm, f/2.8D, and a light source STORZ XENON NOVA 175) and recording up to 5000fps.

Finally, a piezo-resistive transducer (Model 105, Precision Measurement Company, Ann Arbor; full scale 207 kPa, 2.67 mm diameter and 0.508 mm thickness) can be positioned between the vocal folds to provide access to the contact pressure during the closed glottis episode.

SOME EXPERIMENTAL STUDIES

Ex vivo laryngeal articulation

We first present experiments testing the effectiveness of the actuation described in Sec. 2.2. As a first step, the ranges of laryngeal motions (described in terms of vocal folds angle, antero-posterior diameter, inter-arytenoid and mid-vocal folds distance) have been estimated from clinical in vivo fibroscopic videos (N = 28 healthy subjects performing breathing, M1 and M2 phonation; 24 females and 4 males) in order to provide normalized reference data.

Then the designed actuation has been tested on two excised larynx (2 female). As a first step, individual actuations (one muscle simulated at a time) were performed and recorded (see Fig. 6). It is visible that the medialization of the membraneous vocal folds by the Montgomery implants, the adduction of the vocal processes (with a net hourglass shape), the abduction by the simulation of CAP are effective. The CT actuation is less noticeable on the picture but has already been validated in Ref. [START_REF] Lagier | Control of the glottal configuration in ex vivo human models: quantitative anatomy for clinical and experimental practices[END_REF]. As a last comment, it is worth noticing that individual actuations were performed here with no effort applied on the other muscles. Those are expected to stabilize the tissues, i.e., to reduce the magnitude of the displacements obtained in this individual actuation experiment and to reduce the hysteresis visible for some of the muscles (left PCA, right TA and right LCA).

Coactivations of the muscles have also been tested, in a scenario mimicking sequentially breathing, prephonatory posturing, M3 phonation, and successive transitions to M2 and M1 phonations. All these data are still to be compared to the in vivo reference data, to the magnitude of the vocal process displacement reported in Ref. [START_REF] Hunter | Review of range of arytenoid cartilage motion[END_REF] and to the correspondence between laryngeal EMG data and the image-based estimation of vocal folds angle discussed in Ref. [START_REF] Poletto | Correspondence between laryngeal vocal fold movement and muscle activity during speech and nonspeech gestures[END_REF].

Control of the fundamental frequency

Another experiment focused on the control of the fundamental frequency. Air supply line being connected to the excised larynx preparated as in Ref. [START_REF] Lagier | Control of the glottal configuration in ex vivo human models: quantitative anatomy for clinical and experimental practices[END_REF] (i.e., the LCA and PCA actuations described in Sec. 2.2 being replaced by the double adduction technique). Traction and release cycles were performed on the CT actuation leading to changes in the fundamental frequency of the radiated acoustic pressure as displayed on Fig. 7. On this recording, the control is coarse-grained (±1mm steps) and not all commands seem to induce a step in f o . This may be due to a failure in the motor control (the linear actuator does not respond to the electrical input, which seems a common issue on such devices) or either a saturation effect related to the cricoid and the thyroid being in contact on their anterior edges. Another interesting situation is depicted in Fig. 8 where a bifurcation from the M1 laryngeal mechanism to the M2 one occurs about t = 28 s. The opposite situation when releasing the tension at t = 43 s.

Figure 8. Same quantities as in Fig. 7 for a different larynx and different controls on the subglottal pressure and on the CT actuation.

Contact pressure on a larynx with Reinke's edema

As a last example, we refer to Ref. [START_REF] Silva | Contact pressure between the vocal folds in reinke's edema: Experimental observations on an excised human larynx[END_REF] by some of the current authors. It is the first comparison of the contact pressure between vocal folds for two excised human larynges, one of them being pathological with a bilateral loose swelling of the vocal folds diagnosed as bilateral grade 1 Reinke's edema. Both larynges were tested on the testbench (with controls and sensors depicted in Fig. 9) with the airflow controlled in a kind of Mesa di Voce singing exercise (progressively increased, maintained, then decreased), scanning the subglottal pressure range for human phonation.

While most of the usual descriptors (onset and offset subglottal pressures; sound level; fundamental frequency; RMS values for EGG; and Harmonics-to-noise ratio) are comparable and consistent with classical knwoledge on patients suffering from Reinke's edema, the main observation is the extremely high level of contact pressure during stable phonation (up to 296 kPa) in the pathological larynx, compared to the healthy larynx (no more than 34 kPa). This suggests that the phonotrauma may contribute to maintain the lesion after its appearance.

CONCLUSION

The excised human larynx test bench designed and hosted at the Medecine Faculty in Marseille has shown to be a valuable tool to investigate phonation and its pathologies. To the best of our knowledge, this is the first implementation of this full set of actuators and sensors on excised larynx testbeds that offers, in addition to a laryngeal flow control, a physiologically realistic motorized animation of the intrinsic laryngeal musculature with less constrained motions of the cartilages and the soft tissues than other setups prescribing idealized translations and rotations.

Further work is still needed to provide a calibration procedure (in order to automatically determine the bounds and reference state of each actuator) and to map high-level controls to low-level actuation commands.

Figure 1 .

 1 Figure 1. Full experimental setup, showing the positioning of the actuators with respect to the excised larynx.

Figure 2 .

 2 Figure 2. CT actuation: schematics of the pulling system and of the resulting rotation of the cricoid cartilage around the CT joint (left), and realization on an excised larynx (right, from Ref. [15]).

Figure 4 .

 4 Figure 4. LCA actuation: schematics of the pulling system and of the resulting rotation of the arytenoids cartilages around the CA facets (left, top view), and realization on an excised larynx (right, threads as black lines -dashed when within tissues).

Figure 5 .

 5 Figure 5. PCA actuation: schematics of the pulling system (left, rear view), and realization on an excised larynx (right).

Figure 6 .

 6 Figure 6. Effect on the glottis of each individual actuation. Position of the glottis at rest (blue: before activation, green: after activation) and during the activation (background image and red contour).

Figure 7 .

 7 Figure 7. CT actuation commands (u.a., top) and spectrogram of the radiated acoustic pressure (bottom). Effect of the tension of the CT thread on the phonation fundamental frequency.

Figure 9 .

 9 Figure 9. Setup for the estimation of the contact pressure on healthy and pathological larynges.

ACKNOWLEDGEMENTS

Part of this work has been supported by the Mission for Interdisciplinarity CNRS (Call Instrumentation aux limites, 2018). We thank Laurent Sabatier (LMA), Henri Boutin and Djellal Chalabi (STMS) for helping with the imaging of the larynx preparation and its actuation. Furthermore, we would like to thank Pr. Pierre Champsaur for hosting the experiments, and Pr. Antoine Giovanni for supervising the initial studies on the testbench.


Shaoqing Zhou
email: sq.zhou@seu.edu.cn

Member, IEEE Kezhi Wang
email: kezhi.wang@northumbria.ac.uk

Fellow, IEEE Marco Di Renzo
email: marco.direnzo@centralesupelec.fr

Fellow, IEEE Mohamed-Slim Alouini
email: slim.alouini@kaust.edu.sa

W Xu

M.-S Alouini

Spectral and Energy Efficiency of IRS-Assisted MISO Communication with Hardware Impairments

Keywords: Intelligent reflecting surface, hardware impairments, downlink spectral efficiency, energy efficiency. I

In this letter, we analyze the spectral and energy efficiency of an intelligent reflecting surface (IRS)-assisted multiple-input single-output (MISO) downlink system with hardware impairments. An extended error vector magnitude (EEVM) model is utilized to characterize the impact of radio-frequency (RF) impairments at the access point (AP) and phase noise is considered for the imperfect IRS. We show that the spectral efficiency is limited due to the hardware impairments even when the numbers of AP antennas and IRS elements grow infinitely large, which is in contrast with the conventional case with ideal hardware. Moreover, the performance degradation at high SNR is shown to be mainly affected by the AP hardware impairments rather than the phase noise of IRS. We further obtain the optimal transmit power in closed form for energy efficiency maximization. Simulation results are provided to verify these results.

which independently induce phase adjustments on impinging signals to conduct reflecting beamforming.

Significantly different from existing technologies, IRS reconfigures wireless communication environment between transmitter and receiver via programmable and highly controllable intelligent reflection. Moreover, it avoids active radio-frequency (RF) chains and operates passively for short range coverage enhancement so that it can be densely deployed in a flexible way with affordable hardware cost and energy consumption.

Traditional communication theories may no longer be applied because the IRS-assisted wireless system consists of both active and passive components, instead of solely active entities [START_REF] Renzo | Smart radio environments empowered by reconfigurable AI meta-surfaces: An idea whose time has come[END_REF]. Researches on channel estimation, IRS beamforming design and system performance analysis are on the way. Two efficient uplink channel estimation schemes were proposed in [START_REF] Zheng | Intelligent reflecting surface assisted multi-user OFDMA: Channel estimation and training design[END_REF] for IRS-assisted multi-user systems with various channel setups. In [START_REF] Pan | Intelligent reflecting surface aided MIMO broadcasting for simultaneous wireless information and power transfer[END_REF], transmit precoding and passive IRS phase shifts were jointly optimized for simultaneous wireless information and power transfer systems. Ergodic spectral efficiency of an IRS-assisted massive multiple-input multiple-output system was characterized in [START_REF] Han | Large intelligent surface-assisted wireless communication exploiting statistical CSI[END_REF] under Rician fading channel. In [START_REF] Gao | Reconfigurable intelligent surface for MISO systems with proportional rate constraints[END_REF], spectral efficiency of an IRS-aided multi-user system was studied under proportional rate constraints and an iteratively optimizing solution was proposed with closed-form expressions. Secrecy rate was maximized in [START_REF] Shen | Secrecy rate maximization for intelligent reflecting surface assisted multi-antenna communications[END_REF] for an IRS-assisted multi-antenna system by alternately optimizing transmitting covariance and IRS phase shifts. IRS was also shown to be effective in enhancing the performance of cell-edge users [START_REF] Pan | Multicell MIMO communications relying on intelligent reflecting surface[END_REF].

In practice, precise phase control is infeasible at IRS due to hardware limitations and imperfect channel estimation. Corresponding researches are still in their infancy. Discrete phase shifts were considered for IRS-assisted multi-user communication in [START_REF] Di | Hybrid beamforming for reconfigurable intelligent surface based multi-user communications: Achievable rates with limited discrete phase shifts[END_REF], where a hybrid beamforming optimization algorithm was proposed for sum rate maximization. In addition to non-ideal IRS, the impacts of RF impairments at transmitter on the performance of an IRS system have not been clear. To capture the aggregate impacts of various RF impairments, a generalized model named extended error vector magnitude (EEVM) was proposed in [START_REF] Schenk | RF Imperfections in High-Rate Wireless Systems: Impact and Digital Compensation[END_REF] for cellular transmitters.

In this letter, we focus on an IRS-assisted multiple-input single-output (MISO) system with hardware impairments at both access point (AP) and IRS. Theoretical expression of spectral efficiency is derived for this non-ideal case. We discover that the performance is limited even with increasing numbers of elements at both the AP and IRS. The impact of phase noise at IRS diminishes at high SNR. Meanwhile, we obtain a closed-form solution to the optimal power design for maximizing energy efficiency. The optimal power increases with RF impairments.

II. SYSTEM MODEL

A. Signal Model

We consider a MISO downlink system where an IRS consisting of N reflecting elements is deployed to assist the communication from an M-antenna AP to a single-antenna user. The IRS is triggered by an attached smart controller connected to the AP. Denote the reflection matrix of IRS by Θ = diag{ζ 1 e jθ 1 , ζ 2 e jθ 2 , . . . , ζ N e jθ N }, where ζ n ∈ [0, 1] and θ n ∈ [0, 2π) for n = 1, 2, . . . , N are respectively the amplitude reflection coefficient and the phase shift introduced by the nth reflecting element. In practice, each reflecting element is usually designed to maximize the signal reflection. Without loss of generality, we set ζ n = 1 for all n [START_REF] Zheng | Intelligent reflecting surface-enhanced OFDM: Channel estimation and reflection optimization[END_REF]. The direct link between AP and user is blocked by obstacles, such as buildings or human body, which is common in the communication at high-frequency bands, like millimeter wave.

Thus it would be better to deploy IRS at positions where line-of-sight (LoS) communication is ensured for both AP-to-IRS and IRS-to-user links.

Considering the flat-fading model, the channel from the AP to IRS and that from the IRS to user are respectively denoted by H 1 and h H 2 . Both channels are assumed to be LoS, which are represented by

H 1 = αa N (φ a r , φ e r )a H M (φ a t , φ e t ), h H 2 = βa H N (ϕ a t , ϕ e t ), (1)
where α and β are the corresponding strength of path AP-to-IRS and IRS-to-user, φ a r (φ e r ) is the azimuth (elevation) angle of arrival (AoA) at IRS, φ a t (φ e t ) and ϕ a t (ϕ e t ) are the azimuth (elevation) angles of departure (AoD) at AP and IRS, respectively, and a X (ϑ a , ϑ e ) is the array response vector. We consider uniform square planar array (USPA) with √ X × √ X antennas. The array response vector can be written as

a X (ϑ a , ϑ e ) = [1, . . . , e j2π d λ (x sin ϑ a sin ϑ e +y cos ϑ e ) , . . . , e j2π d λ (( √ X-1) sin ϑ a sin ϑ e +( √ X-1) cos ϑ e ) ] T , (2)
where d and λ are the antenna spacing and signal wavelength, and 0 ≤ x, y < √ X are the antenna indices in the planar. Assume that the AP knows the channel state information (CSI) of both H 1 and h H 2 . Channel estimation methods for communication with IRS can be found in [11][12].

With the errors caused by imperfect RF chains at AP, we adopt the EEVM in [START_REF] Schenk | RF Imperfections in High-Rate Wireless Systems: Impact and Digital Compensation[END_REF] to model the transmit signal, which can be written as

x = χws + n RF , (3)
where s is the signal satisfying E[|s| 2 ] = P with P being the transmit power budget, w is the nor- Furthermore, there always exist some phase errors at IRS in implementation. The received signal with phase errors can be modeled as

y = h H 2 ΘH 1 x + u = h H 2 ΘH 1 χws + h H 2 ΘH 1 n RF + u, (4)
where Θ = diag{e j θ1 , e j θ2 , . . . , e j θN } with θn = θ n + θn being the practical phase shift of the nth reflecting element, θn is the phase noise due to the fact, e.g., only discrete phase shifts are possible at IRS, and u is the additive noise with zero mean and variance σ 2 u . Assume that θn is uniformly distributed as U[-δ θ , δ θ] with δ θ ∈ [0, π). Since the distortion noise is independent of channel noise, the received SNR is given by

SNR = P |h H 2 ΘH 1 χw| 2 (h H 2 ΘH 1 )C n RF (h H 2 ΘH 1 ) H + σ 2 u . (5)
Then we have the downlink spectral efficiency as

R = log 2 (1 + SNR). (6)

B. Power Consumption Model

Before we discuss the power consumption, it needs to be emphasized that the IRS does not consume any transmit power due to its nature of passive reflection. The total power consumption is modeled as [START_REF] Huang | Reconfigurable intelligent surfaces for energy efficiency in wireless communication[END_REF]

P T = µP + P C , (7)
where µ = ν -1 with ν being the efficiency of transmit power amplifier considering the RF impairments and P C is the total static hardware power dissipated in all circuit blocks. The establishment of (7) models well under two assumptions: 1) the transmit amplifier operates within its linear region; and 2) the power consumption P C does not rely on the rate of the communication link. Both assumptions are valid in typical wireless systems.

III. SPECTRAL AND ENERGY EFFICIENCY ANALYSIS

In this section, we quantitatively analyze the downlink spectral and energy efficiency and discover the impact of hardware impairments at both the AP and IRS. The ideal spectral and energy efficiency are retrieved as a special case of our analysis and it is presented for comparison.

A. Spectral Efficiency Analysis

Before analyzing the performance, we need to determine the transmit beamforming of AP and the reflecting beamforming of IRS. Since the hardware impairments are unknown and in order to facilitate the design in practice, the two parameters, w and Θ, are optimized by treating the hardware as ideal.

Maximum ratio transmission (MRT) is adopted for transmit beamforming, i.e.,

w = (h H 2 ΘH 1 ) H / h H 2 ΘH 1 . (8)
We identify the optimal reflecting beamforming of IRS by maximizing the received signal power as

Θ opt = arg max Θ |h H 2 ΘH 1 w| 2 (a) = arg max Θ h H 2 ΘH 1 2 = arg max Θ |a H N (ϕ a t , ϕ e t )Θa N (φ a r , φ e r )| 2 a H M (φ a t , φ e t ) 2 (b) = arg max Θ | 0≤x,y< √ N , n= √ N x+y+1 e j2π d λ (xp+yq)+jθn | 2 , (9)
where (a) is obtained by substituting w in ( 8), (b) makes use of a mapping from the two-dimensional index (x, y) to the index n, p = sin φ a r sin φ e rsin ϕ a t sin ϕ e t , and q = cos φ e rcos ϕ e t . Observing [START_REF] Di | Hybrid beamforming for reconfigurable intelligent surface based multi-user communications: Achievable rates with limited discrete phase shifts[END_REF], it is easy to get the optimal phase shift of the nth reflecting element as --→ log 2 1 +

θ opt n = -2π d λ (xp + yq), (10)
P MN 2 η 2 |αβ| 2 sinc 2 (δ ψ )sinc 2 (δ θ) MN 2 |αβ| 2 sinc 2 (δ θ)σ 2 + σ 2 u . ( 11
)
Special Case 1: For the ideal system without any impairments, we let η = 1, σ = 0 and δ ψ = δ θ = 0 in [START_REF] Zheng | Intelligent reflecting surface-enhanced OFDM: Channel estimation and reflection optimization[END_REF]. The downlink spectral efficiency reduces to

R ideal = log 2 1 + P σ 2 u MN 2 |αβ| 2 . ( 12
)
Special Case 2: For high SNR, [START_REF] Zheng | Intelligent reflecting surface-enhanced OFDM: Channel estimation and reflection optimization[END_REF] in Theorem 1 can be further simplified as

R → log 2 P + 2 log 2 η + 2 log 2 sinc(δ ψ ) -log 2 σ 2 . ( 13
)
Remark 1: It is concluded from ( 11) that the non-ideal spectral efficiency increases with η while decreases with parameters δ ψ , σ 2 and δ θ. The impact of the phase rotation at AP in terms of δ ψ is in general more significant than that of the phase noise at IRS in terms of δ θ.

Remark 2: The spectral efficiency in [START_REF] Zheng | Intelligent reflecting surface-enhanced OFDM: Channel estimation and reflection optimization[END_REF] increases with the transmit power approximately in a logarithmic manner similar to the ideal case in ( 12) but with a different scale. Contrary to the ideal case, the performance is ultimately upper bounded for increasing M and N, which is

R(M, N) ≤ R = log 2 (1 + η 2 P σ 2 sinc 2 (δ ψ ))
for all large M and N. Remark 3: An interesting observation from ( 13) is that the spectral efficiency at high SNR is merely limited by the RF impairments at AP rather than the phase noise at IRS, which can be explained from the perspective that the IRS reflecting beamforming simultaneously affects both the desired signal and the distortion noise under the considerations of hardware impairments at AP and LoS channel. It encourages us to use cheap IRS with low-resolution phase shifts without much consideration of performance degradation for large IRS.

B. Energy Efficiency Analysis

The energy efficiency is defined as the ratio of the spectral efficiency to the power consumption, i.e.,

EE

BR/P T where B is the channel bandwidth. We are interested in the performance at high SNR, which can be rewritten as

EE = B(log 2 P + 2 log 2 η + 2 log 2 sinc(δ ψ ) -log 2 σ 2 ) µP + P C . (14)
In the following Theorem 2, we give a closed-form expression of the optimal transmit power maximizing the EE in [START_REF] Xu | Secure communication for spatially sparse millimeter-wave massive MIMO channels via hybrid precoding[END_REF]. Theorem 2: The optimal transmit power to maximize the energy efficiency is the unique solution as

P ⋆ = µP -1 C W (µ -1 e C AP -1 P C ), (15)
where W (x) is the Lambert's W-function and C AP = 2 ln η + 2 ln sinc(δ ψ )ln σ 2 .

Note that for an ideal system without hardware impairments, the optimal transmit power can be similarly derived as

P ⋆ ideal = µP -1 C W (µ -1 e C-1 P C ), (16)
where C = ln(MN 2 |αβ| 2 )ln σ 2 u . We emphasize that the IRS has continuous phase in the ideal case, which increases the static hardware power consumption of IRS.

Remark 4: For P ⋆ in (15) and EE(P ⋆ ) in ( 14), the optimal transmit power increases with more severe RF impairments and the corresponding optimal energy efficiency decreases.

IV. SIMULATION RESULTS

In this section, simulation results are presented to validate the results in Section III. Assume that M = 16, N = 64, η = 0.9, δ ψ = π 18 , σ 2 = 0.1, δ θ = π 8 , α = 0.1, β = 0.5 and µ = 1.1. We plot the downlink spectral efficiency in Theorem 1, special cases and by simulations in Fig. 1. Both non-ideal case in [START_REF] Zheng | Intelligent reflecting surface-enhanced OFDM: Channel estimation and reflection optimization[END_REF] and ideal case in [START_REF] Ning | Channel estimation and transmission for intelligent reflecting surface assisted THz communications[END_REF] increase with the transmit power but by respective scales.

The simplified expression in [START_REF] Huang | Reconfigurable intelligent surfaces for energy efficiency in wireless communication[END_REF] appears to be fairly tight at high SNR.

Transmit Power P (dB) We further assume SNR = P σ 2 = 10 dB. Fig. 1 shows the spectral efficiency versus the number of IRS reflecting elements. As the number goes larger, the hardware impairments lead to limited growth of spectral efficiency, which is consistent with Remark 2, while the ideal case continues increasing logarithmically with the squared number of elements.

In Fig. 2, we give the energy efficiency with various degrees of RF impairments. The optimal transmit power derived in Theorem 2 matches the highest point of the curve well. When the RF impairments become worse, higher optimal transmit power is required while the corresponding energy efficiency decreases. Note that the ideal case may obtain poorer performance than the non-ideal case because of larger static hardware power consumption of continuous-phase IRS.

V. CONCLUSION

In this letter, we demonstrate the downlink spectral and energy efficiency of an IRS-assisted system with hardware impairments. The non-ideal spectral efficiency is upper bounded for large numbers of AP antennas and IRS elements. Specially, the impact of imperfect IRS diminishes at high SNR. The optimal transmit power for maximizing the energy efficiency increases as the RF impairments become more severe.

APPENDIX A PROOF OF THEOREM 1

Applying ( 8) and [START_REF] Schenk | RF Imperfections in High-Rate Wireless Systems: Impact and Digital Compensation[END_REF], we can rewrite the downlink spectral efficiency in [START_REF] Gao | Reconfigurable intelligent surface for MISO systems with proportional rate constraints[END_REF] as

R = log 2 1 + P |h H 2 ΘH 1 χ(h H 2 ΘH 1 ) H | 2 / h H 2 ΘH 1 2 h H 2 ΘH 1 2 σ 2 + σ 2 u R (c) = log 2       1 + P |αβ| 2 N n=1 e j θn 2 |a H M (φ a t , φ e t )χa M (φ a t , φ e t )| 2 M |αβ| 2 N n=1 e j θn 2 a M (φ t ) 2 σ 2 + σ 2 u       = log 2    1 + P η 2 |αβ| 2 M m=1 e jψ(m) 2 N n=1 e j θn 2 /M M|αβ| 2 N n=1 e j θn 2 σ 2 + σ 2 u    , (17)
where --→ sinc 2 (δ θ). The remainder proves that (21) has a unique solution. Define g(P ) µP (ln P + C AP -1). It follows d dP g(P ) = µ(ln P + C AP ) > 0. Thus g(P ) is monotonically increasing with respect to P , which implies that equation ( 21) has at most one solution, which is exactly (15).

  malized beamforming vector at AP, χ = diag{χ(1), χ(2), . . . , χ(M)} with χ(m) = η(m)e jψ(m) for m = 1, 2, . . . , M representing the RF attenuation and phase rotation of the mth RF chain with |η(m)| ≤ 1, and n RF = [n RF[START_REF] Dang | What should 6G be?[END_REF], n RF (2), . . . , n RF (M)] T represents the additive distortion noise with covariance matrix C n RF . The mapping of χ and n RF to particular type(s) of RF impairments, e.g., phase noise, I/Q imbalance and nonlinearity, could be found in[START_REF] Schenk | RF Imperfections in High-Rate Wireless Systems: Impact and Digital Compensation[END_REF] Ch. 7]. For notational simplicity, assume that ψ(m) is uniformly distributed as U[-δ ψ(m) , δ ψ(m) ] with δ ψ(m) ∈ [0, π), n RF (m) ∼ CN (0, σ 2 (m)), and the impairments of all RF chains fall into the same level, i.e., η(m) = η, δ ψ(m) = δ ψ , σ(m) = σ and C n RF = σ 2 I M .

  where x = ⌊(n -1)/ √ N⌋ and y = (n -1) mod √ N , and ⌊•⌋ represents rounding down the value and mod means taking the remainder. Now considering the design of Θ opt in (10) and w(Θ opt ) in (8), we characterize the impacts of both RF impairments at AP and phase noise at IRS on the downlink spectral efficiency in the following Theorem 1. Theorem 1: The downlink spectral efficiency for the massive IRS-assisted MISO with large M and N approaches R a.s.

= π/ 3 RFig. 1 .

 31 Fig. 1. Downlink spectral efficiency versus SNR and N .

9 , 1 η 1 η 1 ηFig. 2 .

 91112 Fig. 2. Energy efficiency versus P .

  |E[e jψ(m) ]|2 (e) = |E[cosψ(m)]| 2 (f ) = sinc 2 (δ ψ ),(18)where (d) applies the Strong Law of Large Numbers and the Continuous Mapping Theorem[START_REF] Xu | Secure communication for spatially sparse millimeter-wave massive MIMO channels via hybrid precoding[END_REF] which indicates that the convergence preserves for continuous matrix functions, (e) uses the symmetry of the odd function sin ψ(m) for ψ(m) ∈ [-δ ψ , δ ψ ], (f ) is obtained by substituting the probability density function of variable ψ(m), i.e., f X (x) = 1 2δ ψ for x ∈ [-δ ψ , δ ψ ], and sinc(x) = sin x x . Similarly, for large N → ∞,

  ) and (19) into (17) completes the proof.APPENDIX B PROOF OF THEOREM 2By calculating the partial derivative of EE in[START_REF] Xu | Secure communication for spatially sparse millimeter-wave massive MIMO channels via hybrid precoding[END_REF], we have∂ ∂P EE = B P -1 (µP + P C )µ(ln P + C AP ) (ln 2)(µP + P C ) 2 . (20)Letting the partial derivative be zero, we haveµP (ln P + C AP -1) = P C ,(21) t=lnP ⇐==⇒ µe t+C AP -1 (t + C AP -1) = e C AP -1 P C ,(g)⇒ t = W (µ -1 e C AP -1 P C ) -C AP + 1,(22)where (g) uses the fact the Lambert's W-function is the inverse function of f (W ) = W e W . Now rearranging (22) yields (15).