Kévin Colin

Xavier Bombois

Laurent Bako

Federico Morelli

Data Informativity for the Identification of particular Parallel Hammerstein Systems

Keywords: System Identication, Data Informativity, Hammerstein Systems, Prediction Error Method, Consistency

  L'archive ouverte pluridisciplinaire

INTRODUCTION

In the prediction error (PE) identication framework, for the estimate to be consistent, the prediction error must be dierent for dierent values of the to-be-identied parameter vector. If the data used for the identication ensures this property, we say that the data are informative.

The data informativity has been extensively studied in the case of Linear Time-Invariant (LTI) systems. This has been done both for Single Input Single Output (SISO) systems [START_REF] Ljung | System identication: Theory for the user[END_REF][START_REF] Bazanella | Necessary and sucient conditions for uniqueness of the minimum in Prediction Error Identication[END_REF][START_REF] Gevers | Identiability and informative experiments in open and closed-loop identication[END_REF][START_REF] Saukoski | System and circuit design for a capacitive MEMS gyroscope[END_REF][START_REF] Gevers | Identication and the Information Matrix: How to Get Just Suciently Rich[END_REF] and for Multiple Inputs Multiple Outputs (MIMO) systems [START_REF] Bazanella | Closed-loop identication of MIMO systems: a new look at identiability and experiment design[END_REF][START_REF] Colin | Closed-loop Identication of MIMO Systems in the Prediction Error Framework: Data informativity Analysis[END_REF]b;[START_REF] Colin | Identication-Based Approach for Electrical Coupling Compensation in a MEMS Gyroscope[END_REF]d].

While the PE identication framework is generally used for LTI systems, it can also be used as an ecient tool to identify certain classes of non-linear systems. This is, e.g., the case for block-oriented systems with static nonlinearities (Hammerstein/Wiener systems) [START_REF] Giri | Block-oriented nonlinear system identication[END_REF].

Adapting the PE framework to this type of systems entails a number of challenges. As an example, the identication problem boils down to a complex non-convex optimization problem. Consequently, a good initialization of this optimization problem is crucial and the best linear approximation framework can be used for this purpose [START_REF] Schoukens | Parametric identication of parallel Hammerstein systems[END_REF][START_REF] Schoukens | Parametric identication of parallel Wiener-Hammerstein systems[END_REF][START_REF] Schoukens | Identication of block-oriented nonlinear systems starting from linear approximations: A survey[END_REF]. Another issue (the one we will consider in the present paper) is that we have to ensure that the data used for the identication are informative (i.e. yields a different prediction error for all values of the to-be-identied parameter vector). Up to our knowledge, this problem has never been studied in the literature.

In this paper, we will tackle data informativity for a particular type of block-oriented systems with two branches described by y(t) = G 0 (z)u(t) + P 0 (z)u n (t) + v(t) where u is the excitation, y the output, v the measurement noise and n ∈ N. This system is a parallel Hammerstein system with one monomial nonlinearity. For instance, the linear dynamics of the micro-electromechanical structure (MEMS) gyroscope with capacitive instrumentation can be modeled quite accurately with such representation with n = 2 [START_REF] Saukoski | System and circuit design for a capacitive MEMS gyroscope[END_REF][START_REF] Kempe | Inertial MEMS: principles and practice[END_REF].

For this type of Hammerstein systems, we derive results in order to verify whether a given set of data is informative or not. We do that by rewriting the system as an equivalent system with one output and two inputs, i.e., u 1 = u and u 2 = u n and we use informativity results for Multiple Input Single Output (MISO) systems to tackle the data informativity problem for the considered parallel Hammerstein system.

Notations. For all matrices A, A T denotes its transpose. The notation 0 n×m refers to the matrix of size n × m full of zeros. For quasistationary signals x [START_REF] Ljung | System identication: Theory for the user[END_REF]], we dene the operator Ē[x(t)] = lim

N →+∞ 1 N N t=1 E[x(t)]
where E is the expectation operator. Finally, for discretetime systems, z denotes the forward-shift operator.

PREDICTION ERROR FRAMEWORK

Consider the parallel Hammerstein system S with one input u and one output y described by

S : y(t) = G 0 (z)u(t) + P 0 (z)u n (t) + H 0 (z)e(t) (1)
where G 0 (z) and P 0 (z) are stable transfer functions, H 0 (z) is a stable, inversely stable and monic transfer function, e a white noise and n an integer that will be assumed to be known and such that n ≥ 2. We will suppose that u and e are independent (i.e. the identication experiment is performed in open loop). Therefore, u n and e are also independent.

As already mentioned in the introduction, we want to identify a model of S by using the PE identication framework. For this purpose, one could consider a parametrized nonlinear model structure

M = {(f (u(t), θ), H(z, η)) | θ = ( θT , η T ) T ∈ D θ } with f (u(t), θ) a nonlinear function in u(t) dened by f (u(t), θ) = G(z, θG )u(t) + P (z, θP )u n (t)
where G(z, θG ), P (z, θP ) and H(z, η) are parametrized rational transfer functions 1 , θ the parameter vector and θ = ( θT G , θT G ) T . The set D θ ⊂ R m restricts θ = ( θT G , θT P , η T ) T to those values for which G(z, θG ) and P (z, θP ) are stable and H(z, η) is stable and inversely stable.

Instead of choosing this model structure, one can notice that the system S in (1) is equivalent to the following MISO system with the input vector u = (u, u n ) T : S :

y(t) = (G 0 (z), P 0 (z))u(t) + H 0 (z)e(t) (2)
Therefore, it can be identied within the following MISO (linear) Box-Jenkins (BJ) model structure M described by

M = {(G(z, θG ), P (z, θP ), H(z, η))| θ = ( θT G , θT P , η T ) T ∈ D θ } (3)
In the sequel, we will denote by µ G (resp. µ P ) the dimension of θG (resp. θP ). Moreover, we will consider the so-called full-order assumption for M , i.e., ∃θ 0 = ( θT 0,G , θT 0,P , η T 0 ) T ∈ D θ such that (G(z, θ0,G ), P (z, θ0,P ), H(z, η 0 )) = (G 0 (z), P 0 (z), H 0 (z)). Finally, we will assume that the model structure M is globally identiable at θ 0 , by considering that there is no pole-zero cancellation at θ 0 [START_REF] Bazanella | Necessary and sucient conditions for uniqueness of the minimum in Prediction Error Identication[END_REF].

Assume that we have a set of N input-output data

Z N = {x(t) = (u T (t), y(t)) T | t = 1, • • • , N } collected on S.
From each (G(z, θG ), P (z, θP ), H(z, η)) ∈ M , we construct the one-step ahead predictor ŷ(t, θ) given by ŷ(t, θ) = W u (z, θ)u(t) + W y (z, θ)y(t) = W (z, θ)x(t) (4) where

W u (z, θ) = H -1 (z, η)(G(z, θG ), P (z, θP )) (5) W y (z, θ) = 1 -H -1 (z, η) (6) W (z, θ) = (W u (z, θ), W y (z, θ)). (7)
Based on the dataset Z N , we compute the optimal parameter vector denoted θN minimizing a least-square criterion on the prediction error (t, θ) = y(t) -ŷ(t, θ): θN = arg min

θ∈D θ V N (θ, Z N ) (8) V N θ, Z N = 1 N N t=1 2 (t, θ) (9)
We want θN to be a consistent estimate of the true parameter vector θ 0 , i.e., θN → θ 0 with a probability equal to 1 when N → ∞. For that, it is important that M is globally identiable at θ 0 and that the data are informative with respect to (w.r.t.) M . The denition of the latter adapted to our problem is given below: Denition 1. Consider the framework dened above with u T = (u, u n ) and with the data

x(t) = (u T (t), y(t)) T 1 H(z, η) is moreover assumed monic, i.e., H(z = ∞) = 1.
collected by applying a quasi-stationary input u to the true system S in (1)-( 2). Consider the model structure M dened in (3) yielding the predictor ŷ

(t, θ) = W (z, θ)x(t). Dene the set ∆ W = {∆W (z) = W (z, θ ) - W (z, θ ) | θ and θ in D θ }.
The data x(t) are said to be informative w.r.t. the model structure M when, for all ∆W (z) ∈ ∆ W , we have

Ē ||∆W (z)x(t)|| 2 = 0 =⇒ ∆W (z) ≡ 0 1×3 (
10) where ∆W (z) ≡ 0 1×3 means that ∆W (e jω ) = 0 1×3 at all or almost all ω.

DATA INFORMATIVITY FOR MISO SYSTEMS IN OPEN-LOOP WITH TWO INPUTS

As already mentioned in the previous section, the nonlinear SISO system S with the input u in (1) can be rewritten as the MISO linear system with the input vector u = (u, u n ) T in (2). Since we want a consistent estimate of (G 0 (z), P 0 (z), H 0 (z)), the objective of this paper is to develop results in order to verify if a given experiment with an excitation u will yield informative data x(t) = (u T (t), y(t)) T w.r.t. the MISO BJ model structure M . To derive these results, we will recall some data informativity results for MISO systems. These results are valid for arbitrary input vectors, i.e., u must not be necessarily equal to (u, u n ) T . For this purpose, let us introduce the following notations:

X(z, θ) = (G(z, θG ), P (z, θP )) and ∆X(z) = X(z, θ ) -X(z, θ ) where θ = ( θT G , θT P ) T . Based on D θ , we dene the set D θ = { θ | θ = ( θT , η T ) T ∈ D θ }. We also dene the set ∆ X = {∆X(z) = X(z, θ ) - X(z, θ ) | θ and θ ∈ D θ }.
Theorem 1. Consider data x(t) = (u T (t), y(t)) T collected on a MISO system (2) in open loop. Consider also a model structure M for this MISO system (see (3)) and the set ∆ X dened above. Then, the data x(t) are informative w.r.t. M if and only if, for all ∆X(z) ∈ ∆ X , we have

Ē ||∆X(z)u(t)|| 2 = 0 =⇒ ∆X(z) ≡ 0 1×2 (11)
Proof. See [START_REF] Colin | Data Informativity for the Open-Loop Identication of Multivariate System in the Prediction Error Framework[END_REF].

In [START_REF] Colin | Data Informativity for the Open-Loop Identication of Multivariate System in the Prediction Error Framework[END_REF]], we have developed ecient conditions to verify whether a given input vector yields informative data (i.e., satises the condition in Theorem 1). As we will see in the sequel, these conditions may be more dicult to verify in the case of an input vector u = (u, u n ) T . Therefore, we present the following lemmas which give simpler (but more conservative) conditions for data informativity.

Lemma 1. Consider Theorem 1. Assume that we collect the data x(t) = (u T (t), y(t)) T by applying a quasistationary input vector u to the system S in (2). Assume that the power spectrum matrix Φ u (ω) of the input u is strictly positive denite for almost all frequencies ω. Then, the data x(t) are informative w.r.t. a full-order BJ model structure M for S (see ( 3)). This result holds whatever the orders of the transfer functions G 0 , P 0 , and H 0 in (2) may be.

Proof. From Parseval theorem, the left-hand side of ( 11) is equivalent to 1 2π π -π ∆X(e jω )Φ u (ω)∆X * (e jω )dω = 0.

Since Φ u (ω) is positive denite for almost all frequencies ω, we have that Ē ||∆X(z)u(t)||2 = 0 always implies ∆X(z) ≡ 0 1×2 irrespectively of the complexity of ∆X ∈ ∆ X , i.e., a full-order BJ model structure M describing any system S. This concludes the proof.

Lemma 1 cannot be used to verify the data informativity for multisine input vectors (since Φ u is never strictly positive denite in this case). Lemma 2 presents a data informativity condition for this type of excitation vectors.

Lemma 2. Consider Theorem 1 and suppose that the data x(t) = (u T (t), y(t)) T are generated with a multisine input vector u = (u 1 , u 2 ) T . Denote by Φ u1 (resp. Φ u2 ) the power spectral density (PSD) of u 1 (resp. u 2 ). Suppose that Φ u1 (resp. Φ u2 ) is non-zero at s 1 (resp. s 2 ) frequencies in ]-π, π]. Finally, suppose that Φ u1 and Φ u2 are both nonzero at s 1,2 frequencies in this interval (s 1,2 can be equal to zero). Then, the data x(t) are informative w.r.t. M (see (3)) if s 1 ≥ µ G and s 2 -s 1,2 ≥ µ P where µ G and µ P are dened below (3).

Proof. Let us rst observe that the multisine u 2 can be decomposed as follows

u 2 = u ( u1) 2 + u (⊥u1) 2 where • u ( u1) 2
is the multisine whose PSD shares the same frequencies with the PSD of u 1 , i.e., the multisine part of u 2 that is totally correlated to u 1 . The PSD of u

( u1) 2 is non-zero in s 1,2 frequencies. • u (⊥u1) 2
is the multisine whose PSD does not share any frequency with the PSD of u 1 , i.e., the multisine part of u 2 that is not correlated to u 1 (and to u ( u1) 2 too by construction). The PSD of u

(⊥u1) 2 is non-zero in s 2 -s 1,2 frequencies.
Consequently, the left hand side of ( 11) is equivalent to the following equation system, for all ∆X(z) = (∆G(z), ∆P (z)

) ∈ ∆ X ,    Ē ||∆G(z)u 1 (t) + ∆P (z)u ( u1) 2 (t)|| 2 = 0 Ē ||∆P (z)u (⊥u1(t)) 2 (t)|| 2 = 0 (12)
We have to prove that (12) implies ∆X = (∆G, ∆P ) ≡ 0 1×2 . The PSD of u

(⊥u1(t)) 2
is non-zero in at least µ P dierent frequencies in ] -π, π]. Therefore, the second equation of (12) implies ∆P (z) ≡ 0 [START_REF] Ljung | System identication: Theory for the user[END_REF][START_REF] Gevers | Informative data: How to get just suciently rich?[END_REF]. By injecting the latter in the rst equation of ( 12), we obtain Ē ||∆G(z)u 1 (t)|| 2 = 0. Since the PSD of u 1 is non-zero in at least µ G frequencies in the set ]-π, π] and so Ē ||∆G(z)u 1 (t)|| 2 = 0 implies that ∆G(z) ≡ 0, which concludes the proof from Theorem 1.

As already mentioned, Lemmas 1 and 2 pertain to an arbitrary input vector u. They can therefore also be used in the case where the input vector u = (u 1 , u 2 ) T is of the form u 1 = u and u 2 = u n . This fact will be used in the sequel to derive data informativity results for the Hammerstein system S in (1).

WHITE GAUSSIAN NOISE EXCITATION

As shown in Appendix A, if the input signal u of ( 1) is chosen as a zero-mean Gaussian white noise, the PSD of u = (u, u n ) T is strictly positive denite at (almost) all frequencies (whatever the value of n ≥ 2). Consequently, we have the following result.

Theorem 2. Consider Theorem 1. Assume that we collect the data x(t) = (u T (t), y(t)) T with u = (u, u n ) T by applying a zero-mean white Gaussian noise u to the system S in (1). Then, the data x(t) are informative with respect to a full-order BJ model structure M describing any system S.

Proof. See Appendix A. This type of stochastic excitation is interesting since it can allow to identify any system of the type (1).

MULTISINE EXCITATION

In this section, u is a sum of m cosinusoids given by

u(t) = m l=1 A l cos(ω l t + φ l ) (13)
where A l > 0 and φ l are respectively the amplitude and the phase-shift of the cosinusoid at the non-zero frequency ω l belonging to the normalized frequency 2 interval ]0, π[. Since u is a multisine, u n is also a multisine, excited at more cosinusoids than u. In [START_REF] Colin | Data Informativity for the Open-Loop Identication of Multivariate System in the Prediction Error Framework[END_REF]], we gave a condition to verify if a given multisine input vector u will yield informative data w.r.t. M by verifying the rank of a matrix depending on the model structure complexity and on the amplitudes, phase-shifts and frequencies of the cosinusoids in u. To apply it to our problem, we have to compute the amplitudes, phase-shifts and frequencies of u n which can become computationally expensive when n and m increase. Fortunately, using Lemma 2, we can derive a sucient condition for data informativity that only requires the knowledge of the number of frequencies present in u n (and not its full expression).

Theorem 3. Consider that the system S in (1) is excited with the multisine (13) where As already mentioned, to use the (sucient) data informativity condition of Theorem 3, we only need to know how many frequencies are present in u n and to compare them with the ones in u (see ( 13)).

ω l (l = 1, • • • , m) ∈ ]0,
For this purpose, one could compute the Fast Fourier Transform (FFT) of the sequence {u n (t) | t = 1, • • • , M } for a suciently large value of M . Another procedure to this end will be derived in the next subsection. In this procedure, we will formally suppose that φ l = 0 (l = 1, • • • , m) in ( 13). See Remark 1 for further details.

Method for the computation of s 2

Let us thus consider (13) with A l > 0, φ l = 0 and ω

l ∈]0, π[ (l = 1, • • • , m) and let us observe that u n (t) = u n-1 (t)u(t), u n-1 (t) = u n-2 (t)u(t), • • • , u 2 (t) = u(t)u(t).
Therefore, to get the cosinusoids of u n , we have rst to determine the ones of u 2 from u, then the ones of u 3 from u 2 , • • • , the ones of u n-1 from u n-2 and nally the ones of u n from u n-1 .

Let us now study how to obtain the frequencies of u p from the ones of u p-1 (p = 2, • • • , n). For this purpose, we know that u p-1 is a multisine and so u p (t) = u p-1 (t)u(t) can be written as the sum of the products of each cosinusoid in u by each cosinusoid in u p-1 . By using the fact that cos(α) cos(β) = 1/2(cos(α + β) + cos(α -β)) ∀(α, β) ∈ R 2 , all the products in this sum can be written as the sum of two cosinusoids.

By doing this, we see that the set of frequencies in u p can be determined by adding and substracting ω l (l = 1, • • • , m) to the frequencies of each cosinusoid present in u p-1 . Of course, in the obtained set, we have to remove the duplicates (e.g., the frequency -ω is equivalent to the frequency ω and the frequency ω + 2kπ with k ∈ Z is equivalent to ω ). We then obtain a set of frequencies in the interval [0, π]. Let us denote by m the number of frequencies in this set. Then, s 2 = 2m if the set does not contain neither the frequency 0 nor π while s 2 = 2m -1 if it contains either 0 or π and s 2 = 2m -2 if it contains both 0 and π.

The above procedure supposes that no terms in the sum of products can cancel out. This is the reason why we suppose that A l > 0 and φ l = 0. Consequently, all the terms in the summation will be characterized by a positive amplitude and a zero phase-shift and no cancellations can then occur.

Let us illustrate the above procedure on an example.

Example 1. Consider n = 3, m = 2 and u(t) = cos(ω 1 t) + cos(ω 2 t) with ω 1 = 0.2 and ω 2 = 0.3. To compute s 2 , we need the frequencies of the cosinusoids in u 3 . For this purpose, we need rst to determine the ones in u 2 from u.

For each frequency in u (i.e., 0.2 and 0.3), we add and subtract ω 1 = 0.2 and ω 2 = 0.3. With ω 1 we obtain the terms 0.4, 0, 0.5, 0.1 and with ω 2 we have 0.5, -0.1, 0.6, 0. By removing the duplicates, we obtain the following frequencies in u 2 : 0, 0.1, 0.4, 0.5 and 0.6.

For each frequency in u 2 (i.e., 0, 0.1, 0.4, 0.5 and 0.6), we add and subtract ω 1 = 0.2 and ω 2 = 0.3. With ω 1 , we obtain the terms 0.2, -0.2, 0.3, -0.1, 0.6, 0.2, 0.7, 0.3, 0.8 and 0.4 and with ω 2 , we have 0.3, -0.3, 0.4, -0.2, 0.7, 0.1, 0.8, 0.2, 0.9 and 0.3. By removing the duplicates, we obtain the 8 following frequencies in u 3 : 0.1, 0.2, 0.3, 0.4, 0.6, 0.7, 0.8 and 0.9, all belonging to [0, π]. We see that m = 8 and that s 2 = 16. Since the frequencies ±0.2 and ±0.3 are also in Φ u , we have that s 1,2 = 4 and thus that s 2 -s 1,2 = 12.

Remark 1. As already mentioned, the above procedure is not guaranteed to yield the right s 2 if the phase-shifts are not all equal to zero. To show this let us consider the following u(t) = 1/2 cos(0.4t) + cos(0.6t) + cos(t + π/2) and let us apply the above procedure for n = 2: we obtain the set of frequencies 0, 0.1, 0.2, 0.8, 0.9, 1, 1.1, and 1.2. However, the frequency 1 does not appear in u 2 due to the cancellations of the terms at this frequency. Note, however, that such a situation is rare in practice and we can therefore carefully use the procedure in this section when the phase-shifts are non-zero. However, a verication using the FFT can be always useful.

Result without computation of s 2

In some cases, we even do not have to compute s 2 to check the data informativity.

Lemma 3. Consider that the excitation of (1) (with n ≥ 2) is given by u(t) = m l=1 A l cos(ω l t) with A l > 0 and ω l ∈]0, π[ (l = 1, • • • , m) and such that the m frequencies nω l are • all dierent between them (modulo 2π).

• all dierent from π and 0 (modulo 2π).

• all dierent from each ω j (j = 1, • • • , m) (modulo 2π).

Then, the data x(t) = (u T (t), y(t)) are informative with respect to M if m ≥ max(µ G /2, µ P /2) where µ G and µ P are dened below (3).

Proof. Following the procedure in Section 5.1, we see that nω 1 , • • • , nω m are frequencies of the multisine u n . Since these m positive frequencies are all dierent (modulo 2π) and dierent from π (modulo 2π), this implies that the PSD of u n will be non-zero at at least 2m frequencies in ]π, π], i.e., s 2 ≥ 2m. Due to the fact that the m frequencies nω l are (modulo 2π) dierent from the frequencies ω j (j = 1, • • • , m) in u, we have also that s 2 -s 1,2 ≥ 2m. The result then follows from Theorem 3.

Remark 2. The idea of Lemma 3 is to give a lower bound for s 2 -s 1,2 which is 2m. However, this bound is conservative. In the example of Section 5.1 where n = 3 and m = 2, we have seen that s 2 -s 1,2 = 12 while the lower bound is equal to 2m = 4.

Synthesis of the results

Let us now summarize the dierent results of this section by giving a general approach to verify the informativity w.r.t. M for a given multisine u. First, we verify that m ≥ µ G /2. If it is not the case, the data x(t) are certainly not informative with respect to M .

If m ≥ µ G /2, we can check whether the (conservative) condition of Lemma 3 is satised. If it is not the case, we compute s 2 and s 1,2 (using the FFT approach or the procedure of Section 5.1) and we verify the (less conservative) condition of Theorem 3. If this condition is still not validated, we need to compute the full expression of u n (with all amplitudes and phase-shifts) and use the results of [START_REF] Colin | Data Informativity for the Open-Loop Identication of Multivariate System in the Prediction Error Framework[END_REF]] on this expression and the one of u (see ( 13)).

REAL LIFE EXAMPLE: MEMS GYROSCOPE

MEMS gyroscope description

The MEMS gyroscope is an inertial sensor used to measure angular rates by using Coriolis eect (see [START_REF] Kempe | Inertial MEMS: principles and practice[END_REF][START_REF] Saukoski | System and circuit design for a capacitive MEMS gyroscope[END_REF] for details of its working principle). By focusing on its main dynamics and from physics laws, the MEMS can be modeled by (1) with n = 2 and where3

• G 0 (z) illustrates a parasite electrical bond between the excitation and measurement capacitive circuits:

G 0 (z) = 10 -2 9.47z -1 + 6.69z -2 -16.21z -3 1 -0.685z -1 + 0.175z -2 -0.0415z -3 • P 0 (z) describes the mechanical motion of the MEMS which is a resonance with a high quality factor:

P 0 (z) = 10 -3 2.39z -1 -5.47z -2 1 -0.743z -1 + z -2
We will consider that e is Gaussian with a variance of 10 -3 and H 0 (z) = 1 for the sake of simplicity and to verify easily the consistency of the estimator. All data used in this example are simulated from the above true system in order to verify the consistency. The sampling frequency f s is 62500Hz. Finally, we choose the full-order model structure M with the same complexity as (G 0 (z), P 0 (z), H 0 (z)). Consequently, it is globally identiable at θ 0 and the number of parameters to be identied in G(z, θ) (resp. in P (z, θ)) is equal to µ G = 6 (resp. µ P = 4).

White Gaussian noise excitation

From Theorem 2, the data are informative with respect to any rational model structure M when u is a white Gaussian noise. We do 100 Monte-Carlo simulations to illustrate the consistency with 100 realizations of the noise e, with a white Gaussian excitation u of variance 1 and a data number N = 5000 for each identication. By computing the mean of the 100 computed parameter vectors, we obtain the following model:

     G(z, θG ) = 10 -2
9.48z -1 + 6.68z -2 -16.21z -3 1 -0.685z -1 + 0.175z -2 -0.0416z -3 P (z, θP ) = 10 -3 2.32z -1 -5.42z -2 1 -0.743z -1 + z -2 The closeness between the identied model and the true one suggests consistency of the estimator.

Multisine excitation

We consider the following multisine excitation for u u(t) = cos(ω 1 t) + cos(ω 2 t) + 0.5 cos(ω 3 t) where ω 1 , ω 2 and ω 3 are the normalized frequencies4 are given by ω 1 = 0.08, ω 2 = 1.19 and ω 3 = 1.92, corresponding to the true frequencies ω1 = 503 rad/s, ω2 = 74362 rad/s and ω3 = 120001 rad/s respectively. First, we have indeed m = 3 ≥ µ G /2. In this case, the sinusoid frequencies have been chosen such that we can verify the condition of Lemma 3. We have indeed that m = 3 ≥ max(µ G /2, µ P /2) and so the data will be informative w.r.t. M with this excitation. Let us verify it by doing 100 Monte-Carlo simulations to illustrate the consistency with 100 realizations of the noise e, with a data number N = 5000 for each identication. By computing the mean of the 100 computed parameter vectors, we obtain the following model:      G(z, θG ) = 10 -2 9.53z -1 + 6.65z -2 -16.21z -3 1 -0.691z -1 + 0.176z -2 -0.0420z -3 P (z, θP ) = 10 -3 2.40z -1 -5.50z -2 1 -0.743z -1 + z -2 Here again, the closeness between the identied model and the true one illustrates the consistency of the estimator.

CONCLUSION

In this paper, we have studied the data informativity with respect to a particular parallel Hammerstein system with an input monomial nonlinearity. We have considered two commonly used signals in Prediction Error Identication: white Gaussian noise and multisine. In the white Gaussian noise case, we can identify any model structure. For multisine excitation, we give some advice on the number of cosinusoids to use. A real life example has been considered to illustrate the results. For future works, we want to study the data informativity property for most complex Hammerstein/Wiener systems.

The normalized frequency ω l is obtained from the true frequency ωl by ω l = ωl /fs.

The given transfer functions G 0 (z) and P 0 (z) are the ones obtained with the approach in[Colin et al. 2019e].

Recall that the normalized frequencies ω l are obtained from the true ones ωl by ω l = ωl /fs.

Appendix A. PROOF OF THEOREM 4

First, since u is a white Gaussian noise with a non-zero variance, u n is also a white noise with a non-zero variance. We will need the expectation value of u n (t) when u(t) is Gaussian. It is given in the next lemma. Lemma 4. [START_REF] Papoulis | Random Variables and Stochastic Processes[END_REF]]. Consider a zero-mean Gaussian variable X with variance equal to σ 2 . Then,

where the operator !! is dened for odd integer p by

We are going to prove that the power spectrum Φ u (ω) of u is positive denite at almost all frequencies ω, which will conclude the proof from Lemma 1. For that, let us calculate Φ u (ω) by taking the Fourier transform of the correlation matrix R u (τ ) given by

. By using the fact that u and u n are white and that u is zero-mean, we have that, from Lemma 4,

For the calculation of Ē[u n (t)u n (t + τ )], we give here the details since it is not as simple as the previous ones.

First, when τ = 0, we have that

We deduce the expression for Ē[u n (t)u n (t + τ )]:

• when n is even

• when n is odd

By taking the Fourier transform of R u (τ ), the power spectrum matrix Φ u (ω) is given by

when n is odd.

Since the power spectrum matrix Φ u (ω) is positive semi-denite at all frequencies ω, let us prove that the determinant of Φ u (ω) is non-zero for almost all ω to prove that it is strictly positive denite at almost all frequencies ω. When

Let us prove that (2n -1)!! -((n -1)!!) 2 = 0 by proving that (2n -1)!!/((n -1)!!) 2 > 1 when n is non-zero and even. The proof that (2n -1)!! -(n!!) 2 > 1 when n is odd is based on the same principle. Let us rst observe that

The numerator and the denominator of the latter are the product of n/2 factors. Since the minimal factor of the numerator is strictly greater than the maximal factor of the denominator for n ≥ 2, then (2n -1)!! > ((n -1)!!) 2 . Therefore, det(Φ u (ω)) > 0 for almost all ω. With Lemma 1, the conclusion follows.


Brigitte Bidégaray-Fesquet

Von Neumann Stability Analysis of Finite Difference Schemes for Maxwell-Debye and Maxwell-Lorentz Equations

.

Before describing and analysing the schemes one by one, we give below the principle of the von Neumann analysis which allows us to study their stability.

Principles of the von Neumann Analysis

Schur and von Neumann polynomials

We define two families of polynomials: Schur and simple von Neumann polynomials.

Definition 1 A polynomial is a Schur polynomial if all its roots r satisfy |r| < 1.

Definition 2 A polynomial is a simple von Neumann polynomial if all its roots r belong to the unit disk (|r| ≤ 1) and all the roots of modulus 1 are simple roots.

Introduction

We address the stability study of finite difference schemes for Maxwell-Debye and Maxwell-Lorentz models. To this aim we selected the same schemes as those already studied by Petropoulos [START_REF] Petropoulos | Stability and phase error analysis of FD-TD in dispersive dielectrics[END_REF], who after having correctly defined characteristic polynomials associated to each scheme, merely computed its roots with a numerical algorithm. This implies having to specify values for the physical parameters which occur in the models as well for the time and space steps chosen for the discretization. The analysis has therefore to be carried out anew for each new material or discretization. We perform here a von Neumann analysis on the characteristic polynomials in their literal form, which yields once and for all stability conditions which are valid for all materials.

Maxwell-Debye and Maxwell-Lorentz Models

Le us consider Maxwell equations without magnetisation (Faraday) ∂ t B(t, x) = -curl E(t, x),

(Ampère) ∂ t D(t, x) = 1 µ 0 curl B(t, x), (1)
where x ∈ R N . This system is closed by the constitutive law of the material

D(t, x) = ε 0 ε ∞ E(t, x) + ε 0 t 0 E(t -τ, x)χ(τ )dτ, (2)
where ε ∞ is the relative permittivity at the infinite frequency and χ the linear susceptibility. If we discretize the integral equation [START_REF] Bidégaray-Fesquet | NAUtil: von Neumann Analysis Utilities[END_REF], we obtain what is called a recursive scheme (see e.g. [START_REF] Luebbers | A frequencydependent finite-difference time-domain formulation for dispersive materials[END_REF], [START_REF] Young | On the dispersion errors related to FD 2 TD type schemes[END_REF]). We can also differentiate Eq. ( 2) to obtain a time-differential equation for D which depends on the specific form of χ. For a Debye medium, this differential equation reads

t r ∂ t D + D = t r ε 0 ε ∞ ∂ t E + ε 0 ε s E, (3)
where t r is the relaxation time and ε s the static relative permittivity. We can derive an equivalent form dealing with the polarisation polarisation P(t, x) = D(t, x) -ε 0 ε ∞ E(t, x), namely

t r ∂ t P + P = ε 0 (ε s -ε ∞ )E. ( 4
)
1

For a Lorentz medium with one resonant frequency ω 1 , we have similarly

∂ 2 t D + ν∂ t D + ω 2 1 D = ε 0 ε ∞ ∂ 2 t E + ε 0 ε ∞ ν∂ t E + ε 0 ε s ω 2 1 E, (5)
where ν is a damping coefficient and

∂ 2 t P + ν∂ t P + ω 2 1 P = ε 0 (ε s -ε ∞ )ω 2 1 E. (6)
Denoting by J the time derivative of P, Maxwell system (1) can be cast as

∂ t B(t, x) = -curl E(t, x), ε 0 ε ∞ ∂ t E(t, x) = 1 µ 0 curl B(t, x) -J(t, x). (7)

Yee Scheme

To discretize Maxwell equations in a passive medium (J = 0), we use Yee scheme [START_REF] Yee | Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media[END_REF], which consists in staggering space and time discretization grids for the different fields. We denote by c ∞ = 1/ √ ε 0 ε ∞ µ 0 the light speed at infinite frequency. If the space step δx is the same in all directions and δt is the time step, the CFL condition is c ∞ δt/δx ≤ 1 in space dimension N = 1 and c ∞ δt/δx ≤ 1/ √ 2 for N = 2 or 3. In dimension 1, we can for example only consider fields E ≡ E x et B ≡ B y which discrete equivalents are E n j ≃ E(nδt, jδx) (with similar notations for D ≡ D x ) and B

n+ 1 2 j+ 1 2 ≃ B((n + 1
2 )δt, (j + 1 2 )δx). Yee scheme for the initial Maxwell system [START_REF] Bidégaray-Fesquet | Stability of FD-TD schemes for Maxwell-Debye and Maxwell-Lorentz equations[END_REF] in variables E, B and D therefore reads

1 δt (B n+ 1 2 j+ 1 2 -B n-1 2 j+ 1 2 ) = - 1 δx (E n j+1 -E n j ), 1 δt (D n+1 j -D n j ) = - 1 µ 0 δx (B n+ 1 2 j+ 1 2 -B n+ 1 2 j- 1 2
).

In the same way, for Maxwell system [START_REF] Strikwerda | Finite Difference Schemes and Partial Differential Equations[END_REF] in variables E, B and J, we have the Yee discretization

1 δt (B n+ 1 2 j+ 1 2 -B n-1 2 j+ 1 2 ) = - 1 δx (E n j+1 -E n j ), ε 0 ε ∞ δt (E n+1 j -E n j ) = - 1 µ 0 δx (B
It may be difficult to localise roots of a polynomial with complicated coefficients. On the other hand, we can turn this difficult problem into the solving of many simpler small problems. To this aim, we construct a polynomial series with strictly decreasing degree. To a polynomial φ defined by

φ(z) = c 0 + c 1 z + • • • + c p z p ,
where c 0 , c 1 . . . , c p ∈ C and c p = 0, we associate its conjugate polynomial φ * which reads

φ * (z) = c * p + c * p-1 z + • • • + c * 0 z p .
Given a polynomial φ 0 , we can define a series of polynomials by recursion

φ m+1 (z) = φ * m (0)φ m (z) -φ m (0)φ * m (z) z .
This series is finite since it is clearly strictly degree decreasing: degφ m+1 < degφ m , if φ m ≡ 0. Besides, we have the following two theorems at our disposal. To localise roots of φ 0 in the unit disk or not, we only have to check conditions at each step m (non zero leading coefficient, |φ m (0)| < |φ * m (0)|, . . . ) until we obtain a negative answer or a polynomial of degree 1.

The proofs of the above results are based on Rouché theorem and are given in [START_REF] Strikwerda | Finite Difference Schemes and Partial Differential Equations[END_REF].

Stability Analysis

The models we consider are linear. They can therefore be analysed in the frequency domain. Hence we assume that the scheme deals with a variable U n j with space dependency in the form

U n j = U n exp(iξ • j),
where ξ et j ∈ R N , N = 1, 2, 3. Let G be the matrix such that U n+1 = GU n and we assume it does not depend on time, nor on δx and δt separately but only on the ratio δx/δt. Let φ 0 be the characteristic polynomial G, then we have the following sufficient stability condition.

Theorem 3 A sufficient stability condition is that φ 0 is a simple von Neumann polynomial.

This condition is not a necessary one. The stability is linked to the fact that U n = G n U 0 and corresponds to the boundedness of the iterates G n of the matrix G. The case of multiple unit modulus roots can give rise to iterates of G which are bounded (e.g. for the identity matrix) or not. For example 1 0 0 1 n = 1 0 0 1 is bounded, and 1 1 0

1 n = 1 n 0 1 is not bounded.
This case occurs for the schemes we are dealing with and have to be treated separately, without the help of the von Neumann analysis, which handles characteristic polynomials and not the matrices they stem from, which induces a loss of information.

Debye Type Media

For Debye type media, we study two schemes. The first one is due to Joseph et al. [START_REF] Joseph | Direct time integration of Maxwell's equations in linear dispersive media with absorption for scattering and propagation of femtosecond electromagnetic pulses[END_REF] and consists in coupling Maxwell equations in variables E, B and D with the Debye model linking E and D. The second is due to Young [START_REF] Young | Propagation in linear dispersive media: Finite difference time-domain methodologies[END_REF] and couples Maxwell equations in variables E, B and J with the Debye model linking E, P and J.

Joseph et al. Model

Model Setting

Maxwell system ( 8) is closed by a discretization of the Debye model (3), namely

ε 0 ε ∞ t r E n+1 j -E n j δt + ε 0 ε s E n+1 j + E n j 2 = t r D n+1 j -D n j δt + D n+1 j + D n j 2 . ( 10
)
System ( 8)-( 10) deals with the variable

U n j = (c ∞ B n-1 2 j+ 1 2 , E n j , D n j /ε 0 ε ∞ ) t = (B n-1 2 j+ 1 2 , E n j , D n j ) t
and reads

B n+ 1 2 j+ 1 2 -B n-1 2 j+ 1 2 = - c ∞ δt δx (E n j+1 -E n j ), D n+1 j -D n j = - c ∞ δt δx (B n+ 1 2 j+ 1 2 -B n+ 1 2 j- 1 2
),

E n+1 j -E n j + ε s ε ∞ δt 2t r (E n+1 j + E n j ) = D n+1 j -D n j + δt 2t r (D n+1 j + D n j ).
We see that this formulation contains dimensionless parameters:

λ = c ∞ δt/δx CFL constant, δ = δt/2t r normalised time step, ε ′ s = ε s /ε ∞ normalised static permittivity.
We write this system into the explicit form

B n+ 1 2 j+ 1 2 = B n-1 2 j+ 1 2 -λ(E n j+1 -E n j ), D n+1 j = D n j -λ(B n-1 2 j+ 1 2 -B n-1 2 j-1 2 ) + λ 2 (E n j+1 -2E n j + E n j-1 ), (1 + δε ′ s )E n+1 j = (1 -δε ′ s )E n j + (1 + δ)λ 2 (E n j+1 -2E n j + E n j-1 ) +(1 + δ)λ(B n-1 2 j+ 1 2 -B n-1 2 j-1 2 ) + 2δD n j ,
which yields the amplification matrix

G =    1 -λ(e iξ -1) 0
-(1+δ)λ(1-e -iξ ) 1+δε ′ s (1-δε ′ s )+(1+δ)λ 2 (e iξ -2+e -iξ ) 1+δε ′ s 2δ 1+δε ′ s -λ(1 -e -iξ ) λ 2 (e iξ -2 + e -iξ ) 1    .
We set σ = λ(e iξ -1) and q = |σ| 2 = -λ 2 (e iξ -2 + e -iξ ) = 4λ 2 sin 2 (ξ/2). With these notations G reads

G =    1 -σ 0 (1+δ)σ * 1+δε ′ s (1-δε ′ s )-(1+δ)q 1+δε ′ s 2δ 1+δε ′ s σ * -q 1    .

Computation of the Characteristic Polynomial

The characteristic polynomial of G is equal to

P (Z) = 1 1 + δε ′ s Z -1 σ 0 -(1 + δ)σ * (1 + δε ′ s )Z -(1 -δε ′ s ) + (1 + δ)q -2δ -σ * q Z -1 = 1 1 + δε ′ s (Z -1) (1 + δε ′ s )Z -(1 -δε ′ s ) + (1 + δ)q -2δ q Z -1 + q 1 + δ -δ 1 Z -1 = 1 1 + δε ′ s (Z -1) [1 + δε ′ s ]Z 2 -[2 -(1 + δ)q]Z + [1 -δε ′ s -(1 -δ)q] +q {[1 + δ]Z -[1 -δ]} .
The characteristic polynomial is proportional to

φ 0 (Z) = [1 + δε ′ s ]Z 3 -[3 + δε ′ s -(1 + δ)q]Z 2 + [3 -δε ′ s -(1 -δ)q]Z -[1 -δε ′ s ].

Von Neumann Analysis

From the polynomial φ 0 , we perform the recursive construction of the above-mentioned series of polynomials. We therefore define

φ * 0 (Z) = [1 + δε ′ s ] -[3 + δε ′ s -(1 + δ)q]Z + [3 -δε ′ s -(1 -δ)q]Z 2 -[1 -δε ′ s ]Z 3 .
The condition |φ 0 (0)| < |φ * 0 (0)| is valid. We define by recursion

φ 1 (Z) = 1 Z {φ * 0 (0)φ 0 (Z) -φ 0 (0)φ * 0 (Z)} = 2δ{2ε ′ s Z 2 -[4ε ′ s -(ε ′ s + 1)q]Z + [2ε ′ s -(ε ′ s -1)q]}, φ * 1 (Z) = 2δ{2ε ′ s -[4ε ′ s -(ε ′ s + 1)q]Z + [2ε ′ s -(ε ′ s -1)q]Z 2 }.
Since ε s ≥ ε ∞ , we have ε ′ s ≥ 1 and the quantity ε ′ s -1 is nonnegative. If q = 0 or ε ′ s = 1, we have exactly |φ 1 (0)| = |φ * 1 (0)|, and these specific cases have to be treated separately (see below). In the opposite case, condition

|φ 1 (0)| < |φ * 1 (0)| reverts to (ε ′ s -1)q < 4ε ′ s .
It is reasonable to assume we will not obtain a better result than with the raw Yee scheme (λ ≤ 1) and therefore q ∈ [0, 4]. In that case, and provided q = 0, we do have

|φ 1 (0)| < |φ * 1 (0)|. Moreover the degree of polynomial φ 1 is 2. Last φ 2 (Z) = 1 Z (φ * 1 (0)φ 1 (Z) -φ 1 (0)φ * 1 (Z)) = 4δ 2 (ε ′ s -1)q (4ε ′ s -(ε ′ s -1)q)Z -(4ε ′ s -(ε ′ s + 1)q) .
Always in the case when ε ′ s > 1 and q ∈]0, 4], the leading coefficient

4ε ′ s -(ε ′ s -1)q = ε ′ s (4 -q) + q = 0,
and is degree of φ 2 is 1. The root of φ 2 is

Z = 4ε ′ s -(ε ′ s + 1)q 4ε ′ s -(ε ′ s -1)q .
The modulus of this root is strictly lower than 1 if q = 4 and therefore φ 2 and hence φ 0 are Schur polynomials thanks to Theorem 1. If q = 4, the root of φ 2 is -1 and φ 2 and hence φ 0 are simple von Neumann polynomials thanks to Theorem 2. In both cases, we obtain the stability with the only assumption λ < 1, provided we treat the above-mentioned special cases.

Case q = 0

The case when q = 0 corresponds to the characteristic polynomial

φ 0 (Z) = (Z -1) 2 ([1 + δε ′ s ]Z -[1 -δε ′ s ])
which is not a simple von Neumann one. We shall therefore study the amplification matrix directly, which is then simply

G =     1 0 0 0 1 -δε ′ s 1 + δε ′ s 2δ 1 + δε ′ s 0 0 1     .
We clearly see that the eigenvectors corresponding to the eigenvalue 1 are in two stable eigensubspaces.

The other eigenvalue has a modulus strictly lower than 1. Iterates of this matrix are therefore bounded. This conclusion is valid for all ε ′ s ≥ 1.

3.1.5 Case ε ′ s = 1 The case when ε ′ s = 1 gives rise to a breaking of condition |φ 1 (0)| < |φ * 1 (0)|.
We shall therefore study directly the nature of φ 1 without carrying recursion over. We have

φ 1 (Z) = 4δ{Z 2 -[2 -q]Z + 1},
which determinant is q(q -4) and is therefore negative if q ∈]0, 4[. The roots of φ 1 are therefore complex conjugate, distinct and their modulus is 1 (their product is equal to 1). The polynomial φ 1 is therefore a simple von Neumann one, and φ 0 also.

Case q = 4

The last case we have to treat is ε ′ s = 1 and q = 4, where -1 is a double root of φ 1 . It is also a double root of φ 0 which reads

φ 0 (Z) = (Z + 1) 2 ([1 + δ]Z -[1 -δ]).
Hence we study directly the amplification matrix which reads simply

G =   1 -σ 0 σ * 1-δ 1+δ -q 2δ 1+δ σ * -q 1   .
There is no trivial splitting in two distinct eigensubspaces. We compute the eigenvectors associated to the eigenvalue -1. To this aim we solve

(G + Id)V =   2 -σ 0 σ * 2 1+δ -q 2δ 1+δ σ * -q 2   V = 0 ⇐⇒   2 -σ 0 0 1 -1 σ * -q 2   V = 0
and we only find one eigendirection, that of V = (σ, 2, 2) t . A minimal two-dimensional eigensubspace is therefore associated to the eigenvalue -1 and iterates G n are linearly increasing with n. Hence we conclude to instability when q = 4 and ε ′ s = 1.

Synthesis for the Debye-Joseph et al. Model

The scheme ( 8)- [START_REF] Young | On the dispersion errors related to FD 2 TD type schemes[END_REF] for the one-dimensional Maxwell-Debye equation is stable with the condition

δt ≤ δx/c ∞ if ε s > ε ∞ and δt < δx/c ∞ if ε s = ε ∞ .
We have already seen in this first example different types of arguments to conclude to stability: the generic case (q ∈]0, 4[ and ε s > 1) gives rise to a Schur polynomial via Theorem 1, the cases q ∈]0, 4[ and ε ′ s = 1 or q = 4 and ε ′ s > 1 to a simple von Neumann polynomial via Theorem 2 and last, the case q = 0 to a (not simple) von Neumann polynomial, but with a double eigenvalue that operates on two stable and distinct eigensubspaces. We have also encountered an instable case when ε ′ s = 1 and q = 4 which nevertheless corresponds to a (non simple) von Neumann polynomial.

Young Model

Model Setting

Maxwell system (9) is closed by two discretizations of Debye equation (4), namely

t r P n+ 1 2 j -P n-1 2 j δt = - P n+ 1 2 j + P n-1 2 j 2 + ε 0 (ε s -ε ∞ )E n j , (11)
and

t r J n+ 1 2 j = -P n+ 1 2 j + ε 0 (ε s -ε ∞ ) E n+1 j + E n j 2 . ( 12
)
Although we make use of J

n+ 1 2 j
in the description of the scheme, this is not a genuine variable and the system ( 9)-( 11)-( 12) deals with the variable

U n j = (c ∞ B n-1 2 j+ 1 2 , E n j , P n-1 2 j /ε 0 ε ∞ ) t = (B n-1 2 j+ 1 2 , E n j , P n-1 2 j
) t , and reads

B n+ 1 2 j+ 1 2 -B n-1 2 j+ 1 2 = -λ(E n j+1 -E n j ), E n+1 j -E n j = -λ(B n+ 1 2 j+ 1 2 -B n+ 1 2 j-1 2 ) + 2δP n+ 1 2 j -δα(E n+1 j + E n j ), P n+ 1 2 j -P n-1 2 j = -δ(P n+ 1 2 j + P n-1 2 j
) + 2δαE n j . In this system, apart from the notations λ, δ which we have already defined, we have introduced the dimensionless parameter α = ε ′ s -1 which is , as we already mentioned, a non negative parameter. We rewrite this system in the explicit form

B n+ 1 2 j+ 1 2 = B n-1 2 j+ 1 2 -λ(E n j+1 -E n j ), (1 + δα)E n+1 j = (1 -δα)E n j -λ(B n-1 2 j+ 1 2 -B n-1 2 j-1 2 ) + λ 2 (E n j+1 -2E n j + E n j-1 ) +2δ 1 -δ 1 + δ P n-1 2 j + 4δ 2 α 1 + δ E n j , (1 + δ)P n+ 1 2 j = (1 -δ)P n-1 2 j
+ 2δαE n j , from which stems the amplification matrix

G =    1 -σ 0 σ * 1+δα (1+δ)(1-δα)+4δ 2 α-(1+δ)q (1+δ)(1+δα) 1-δ 1+δ 2δ 1+δα 0 2δα 1+δ 1-δ 1+δ    .

Computation of the Characteristic Polynomial

The characteristic polynomial G is

P (Z) = Z -1 σ 0 -σ * 1+δα Z -(1+δ)(1-δα)+4δ 2 α-(1+δ)q (1+δ)(1+δα) -1-δ 1+δ 2δ 1+δα 0 -2δα 1+δ Z -1-δ 1+δ .
To reduce computations, we set Y = Z -1, which yields

(1 + δ) 2 (1 + δα)P (Z) = Y σ 0 -(1 + δ)σ * (1 + δ)(1 + δα)Y + 2δα(1 -δ) + (1 + δ)q -2δ(1 -δ) 0 -2δα (1 + δ)Y + 2δ = Y σ 0 -(1 + δ)σ * (1 + δ)(1 + δα)Y + (1 + δ)q (1 + δ)(1 -δ)Y 0 -2δα (1 + δ)Y + 2δ
.

We see that (1 + δ) is a factor in both sides and therefore

(1 + δ)(1 + δα)P (Z) = Y (1 + δα)Y + q (1 -δ)Y -2δα (1 + δ)Y + 2δ + q 1 -2δ(1 -δ) 0 (1 + δ)Y + 2δ = Y [(1 + δ)(1 + δα)]Y 2 + [2δ(1 + α) + (1 + δ)q]Y + [2δq] +q(1 + δ) {[1 + δ]Y + [2δ]} = [(1 + δ)(1 + δα)]Y 3 + [2δ(1 + α) + (1 + δ)q]Y 2 +[(1 + 3δ)q]Y + [2δq].
The characteristic polynomial is proportional to

φ 0 (Z) = [(1 + δα)(1 + δ)]Z 3 -[3 + δ + δα + 3δ 2 α -(1 + δ)q]Z 2 +[3 -δ -δα + 3δ 2 α -(1 -δ)q]Z -[(1 -δα)(1 -δ)].

Von Neumann Analysis

Condition |φ 0 (0)| < |φ * 0 (0)| is valid without any assumption. We define by recursion

φ 1 (Z) = 2δ{[2(1 + α)(1 + δ 2 α)]Z 2 -[4(1 + α)(1 + δ 2 α) -(2 + α + δ 2 α)q]Z +[2(1 + α)(1 + δ 2 α) -α(1 -δ 2 )q]}.
The case when δ 2 > 1 does not allow to fulfil the condition |φ 1 (0)| < |φ * 1 (0)|. We will assume therefore for the von Neumann analysis that δ < 1, which bounds the time step with respect to the time delay t r . This is reasonable from the point of view of modelling: we cannot approximate the delay equation with too large a time step. Such an assumption was however not necessary for the Joseph et al. scheme. The equality case |φ 1 (0)| = |φ * 1 (0)| is obtained when q = 0, α = 0 (i.e. ε ′ s = 1) or δ = 1. These cases shall be treated separately again. If α > 0, q > 0 and δ < 1, then

|φ 1 (0)| < |φ * 1 (0)| is equivalent to α(1 -δ 2 )q < 4(1 + α)(1 + δ 2 α), which is clearly true if q ∈]0, 4].
Besides, the degree of polynomial φ 1 is 2. In the general case (α > 0, q ∈]0, 4] and δ < 1), we then compute φ 2

φ 2 (Z) = 4δ 2 α(1 -δ 2 )q{[4(1 + α)(1 + δ 2 α) -α(1 -δ 2 )q]Z -[4(1 + α)(1 + δ 2 α) -(2 + α + δ 2 α)q]} .
We split the study according to the sign of φ 2 (0) (φ * 2 (0) is clearly always positive) and in both cases |φ 2 (0)| < |φ * 2 (0)|, for q ∈]0, 4]. The root of φ 2 therefore belongs to the interval ] -1, 1[ and φ 0 is a Schur polynomial. Hence we obtain the stability with the assumptions λ < 1 and δ < 1, provided we treat the above-mentioned specific cases.

Case q = 0

The case when q = 0 corresponds to the characteristic polynomial

φ 0 (Z) = (Z -1) 2 (Z - (1 -δ)(1 -δα) (1 + δ)(1 + δα) )
and is not a simple von Neumann one. The amplification matrix reads

G =    1 0 0 0 (1+δ)(1-δα)+4δ 2 α (1+δ)(1+δα) 1-δ 1+δ 2δ 1+δα 0 2δα 1+δ 1-δ 1+δ    .
We clearly see that the eigenvectors corresponding to the eigenvalue 1 are in two stable eigensubspaces. Iterates of this matrix are therefore bounded. This conclusion is once more valid in the limit cases δ = 1 and α = 0.

Case ε ′

s = 1 We notice that φ 0 is the same as that for the Joseph et al. model for ε ′ s = 1. Polynomial φ 0 is therefore a simple von Neumann polynomial for q = 4 (see above). The value of δ does not play any rôle here. This corresponds to different amplification matrices, operating on different sets of variables, the link between both formulations being not straightforward. Case q = 4 has therefore to be treated anew.

Case q = 4

If q = 4, only the case when α = 0 has not been treated by the general study and -1 is once more a double eigenvalue

(G + Id)V =   2 -σ 0 σ * 2 -q 2δ 1-δ 1+δ 0 0 2 1+δ   V = 0 ⇐⇒   2 -σ 0 σ * 2 -q 0 0 0 2 1+δ   V = 0,
and the only eigendirection is that of V = (σ, 2, 0) t , which gives rise to linearly increasing iterates G n and to instabilities. The fact that δ = 1 or not does not play any rôle in this argument.

Case δ = 1

There remains to study the case δ = 1 for q ∈]0, 4] and α > 0. Then Z = 0 is a trivial root of φ 0 , which simply reads

φ 0 (Z) = 2(1 + α)Z{Z 2 -[2 - q 1 + α ]Z + 1}.
The discriminant of the second order factor is ∆ = q (1+α) 2 (q -4(1 + α)) < 0. We therefore have two distinct complex conjugate eigenvalues of modulus 1. Polynomial φ 0 is a simple von Neumann polynomial.

Synthesis for the Debye-Young Model

The scheme ( 9)-( 11)-( 12) for the one-dimensional Maxwell-Debye equation is stable with the condition

δt ≤ min(δx/c ∞ , 2t r ) if ε s > ε ∞ and δt < δx/c ∞ if ε s = ε ∞ .
If ε s > ε ∞ , the stability condition is more restrictive for the Young scheme than for Joseph et al. scheme. The obtained bound is also related to the good approximation of Debye equation.

Anharmonic Lorentz Type Media

For Lorentz type media, we study three schemes. The first one is due to Joseph et al. [START_REF] Joseph | Direct time integration of Maxwell's equations in linear dispersive media with absorption for scattering and propagation of femtosecond electromagnetic pulses[END_REF] and consists in coupling Maxwell equations in the variables E, B and D with the Lorentz model linking E and D. The second and third are due to Kashiwa et al. [START_REF] Kashiwa | A treatment by the FD-TD method of the dispersive characteristics associated with orientation polarization[END_REF] and Young [START_REF] Young | Propagation in linear dispersive media: Finite difference time-domain methodologies[END_REF] respectively and both couple Maxwell equations in the variables E, B and J with the Lorentz model linking E, P and J. They differ in the choice of the time discretization of J.

We restrict here to the anharmonic case for which the damping ν is non-zero. The harmonic case (ν = 0) is treated with the same schemes but the analysis happens to be much more technical. To keep proofs readable in the general case we postpone the harmonic case to the next section.

Joseph et al. Model

Model Setting

Maxwell system ( 8) is closed by a discretization of the Lorentz equation ( 5), namely

ε 0 ε ∞ E n+1 j -2E n j + E n-1 j δt 2 + νε 0 ε ∞ E n+1 j -E n-1 j 2δt + ε 0 ε s ω 2 1 E n+1 j + E n-1 j 2 = D n+1 j -2D n j + D n-1 j δt 2 + ν D n+1 j -D n-1 j 2δt + ω 2 1 D n+1 j + D n-1 j 2 . ( 13
)
The explicit version of system ( 8)-(13) does not use explicitly the variable D n-1 j . Indeed we can use the explicit formula to compute D n+1 j -D n j and the implicit one to compute D n j -D n-1 j and therefore the system deals with the variable

U n j = (c ∞ B n-1 2 j+ 1 2 , E n j , E n-1 j , D n j /ε 0 ε ∞ ) t = (B n-1 2 j+ 1 2 , E n j , E n-1 j , D n j ) t
and reads

B n+ 1 2 j+ 1 2 = B n-1 2 j+ 1 2 -λ(E n j+1 -E n j ), D n+1 j = D n j -λ(B n-1 2 j+ 1 2 -B n-1 2 j-1 2 ) + λ 2 (E n j+1 -2E n j + E n j-1 ), (1 + δ + ωε ′ s )E n+1 j = 2E n j + (1 + δ + ω)λ 2 (E n j+1 -2E n j + E n j-1 ) -(1 -δ + ωε ′ s )E n-1 j -2δλ(B n-1 2 j+ 1 2 -B n-1 2 j-1 2
) + 2ωD n j .

In this system, apart from the already used notations λ and ε ′ s , we have denoted

δ = δtν/2 normalised time step, ω = ω 2 1 δt 2 /2
square of the normalised frequency.

The amplification matrix of the system is

G =      1 -σ 0 0 2δσ * 1+δ+ωε ′ s 2-q(1+δ+ω) 1+δ+ωε ′ s -1-δ+ωε ′ s 1+δ+ωε ′ s 2ω 1+δ+ωε ′ s 0 1 0 0 σ * -q 0 1      .

Computation of the Characteristic Polynomial

The characteristic polynomial of G is equal to

P (Z) = Z -1 σ 0 0 -2δσ * 1+δ+ωε ′ s Z -2-q(1+δ+ω) 1+δ+ωε ′ s 1-δ+ωε ′ s 1+δ+ωε ′ s -2ω 1+δ+ωε ′ s 0 -1 Z 0 -σ * q 0 Z -1 . Therefore (1 + δ + ωε ′ s )P (Z) = Z -1 σ 0 0 -2δσ * (1 + δ + ωε ′ s )Z -2 + (1 + δ + ω)q 1 -δ + ωε ′ s -2ω 0 -1 Z 0 -σ * q 0 Z -1 = Z -1 0 0 -2δσ * 1 -δ + ωε ′ s -2ω -σ * 0 Z -1 + Z Z -1 σ 0 -2δσ * (1 + δ + ωε ′ s )Z -2 + (1 + δ + ω)q -2ω -σ * q Z -1 = (Z -1) 2 [1 -δ + ωε ′ s ] +Z(Z -1){[(1 + δ + ωε ′ s )Z -2 + (1 + δ + ω)q][Z -1] + 2ωq} + qZ[2δ(Z -1) + 2ω].
The characteristic polynomial is proportional to

φ 0 (Z) = [1 + δ + ωε ′ s ]Z 4 -[4 + 2δ + 2ωε ′ s -(1 + δ + ω)q]Z 3 + [6 + 2ωε ′ s -2q]Z 2 -[4 -2δ + 2ωε ′ s -(1 -δ + ω)q]Z + [1 -δ + ωε ′ s ].

Von Neumann Analysis

We successively compute

φ 1 (Z) = 2δ{[2(1 + ωε ′ s )]Z 3 -[6 + 4ωε ′ s -(2 + ω(1 + ε ′ s ))q]Z 2 +[6 + 2ωε ′ s -2q]Z -[2 + ω(ε ′ s -1)q]}, φ 2 (Z) = 4δ 2 ω{[4ε ′ s (2 + ωε ′ s ) -4(ε ′ s -1)q -ω(ε ′ s -1) 2 q 2 ]Z 2 -[8ε ′ s (2 + ωε ′ s ) -4((ε ′ s -1) -ε ′ s (2 + ωε ′ s ))q + 2(ε ′ s -1)q 2 ]Z +[4ε ′ s (2 + ωε ′ s ) -4(ε ′ s -1)(2 + ωε ′ s )q + (2 + ω(1 + ε ′ s ))(ε ′ s -1)q 2 ]}, φ 3 (Z) = 64δ 4 ω 2 (ε ′ s -1)(1 + ωε ′ s )q(2 -q) × × {[4ε ′ s (2 + ωε ′ s ) -(ε ′ s -1)(6 + 2ωε ′ s ))q + (ε ′ s -1)(1 + ω)q 2 ]Z -[4ε ′ s (2 + ωε ′ s ) -2((ε ′ s -1) + ε ′ s (2 + ωε ′ s ))q + (ε ′ s -1)q 2 ]}. The root of φ 3 is Z = (2 -q)(2ε ′ s (2 + ωε ′ s ) -(ε ′ s -1)q) (2 -q)(2ε ′ s (2 + ωε ′ s ) -(ε ′ -1)q) + 2(2 + ωε ′ s )q + (ε ′ s -1)(1 + ω)q 2 ,
form on which we easily see that Z remains of modulus < 1 if q ∈]0, 2[, which corresponds to the condition for a multi-dimensional Yee scheme (λ ≤ 1/ √ 2). We notice from now on that we shall treat cases q = 0 and q = 2 apart because φ 3 ≡ 0. In the general case, we have to check the intermediate properties. First, the degree of polynomials φ 1 and φ 2 is 3 and 2 respectively. The degree of polynomial φ 3 is 1 provided ε ′ s > 1. We shall treat the case ε ′ s = 1 apart. In the general case (q = 0 and ε ′ s > 1), there remains to check the estimates between

φ 0 (0) = 1 -δ + ωε ′ s , φ * 0 (0) = 1 + δ + ωε ′ s , φ 1 (0) = 2δ -2 -ω(ε ′ s -1)q , φ * 1 (0) = 2δ 2(1 + ωε ′ s ) , φ 2 (0) = 4δ 2 ω 4ε ′ s (2 + ωε ′ s ) -4(ε ′ s -1)(2 + ωε ′ s )q + (2 + ω(1 + ε ′ s ))(ε ′ s -1)q 2 , φ * 2 (0) = 4δ 2 ω 4ε ′ s (2 + ωε ′ s ) -4(ε ′ s -1)q -ω(ε ′ s -1) 2 q 2 .
It is clear that for q ∈]0, 2[, we have

|φ 0 (0)| < |φ * 0 (0)| and |φ 1 (0)| < |φ * 1 (0)|. A simple calculation shows that φ 2 (0) = 4δ 2 ω (2 -q) 2 (2(1 + ω)(ε ′ s -1) + ω(ε ′ s -1) 2 ) + 8 + 4ω + 4ω(ε ′ s -1)q , φ * 2 (0) = 4δ 2 ω (2 -q)(4(ε ′ s -1) + ω(ε ′ s -1) 2 q) + 8 + 4ω + 8ω(ε ′ s -1) ,
form on which we readily see that both quantities are positive. Besides

φ * 2 (0) -φ 2 (0) = 8δ 2 ω(1 + ωε ′ s )(ε ′ s -1)q(2 -q) > 0.
We therefore checked all the assumptions. In the general case, φ 0 is a Schur polynomial.

Case q = 0

The case q = 0 gives anew rise to a separate study. We have

φ 0 (Z) = (Z -1) 2 (Z -1) 2 + δ(Z 2 -1) + ωε ′ s (Z 2 + 1) .
The corresponding amplification matrix is

      1 0 0 0 0 2 1 + δ + ωε ′ s - 1 -δ + ωε ′ s 1 + δ + ωε ′ s 2ω 1 + δ + ωε ′ s 0 1 0 0 0 0 0 1       .
Once more, 1 is a double root but in two distinct eigensubspaces. There remains to check that the other factor of the polynomial, namely

ψ 0 (Z) = [1 + δ + ωε ′ s ]Z 2 -2Z + [1 -δ + ωε ′ s ],
is a Schur (or a simple von Neumann) one. We do have |ψ 0 (0)| < |ψ * 0 (0)| and we compute

ψ 1 (Z) = 4δ{[1 + ωε ′ s ]Z -1}.
The modulus of both remaining eigenvalues is strictly less that 1 and iterates of the amplification matrix are bounded. This holds whatever the value of ε ′ s .

Case q = 2

In the specific case q = 2, φ 0 reads

φ 0 (Z) = [1 + δ + ωε ′ s ]Z 4 -[2 + 2ω(ε ′ s -1)]Z 3 + [2 + 2ωε ′ s ]Z 2 -[2 + 2ω(ε ′ s -1)]Z + [1 -δ + ωε ′ s ] = (Z 2 + 1){[1 + δ + ωε ′ s ]Z 2 -[2 + 2ω(ε ′ s -1)]Z + [1 -δ + ωε ′ s ]},
which has ±i as simple roots. This is therefore a good candidate to be a simple von Neumann polynomial. The remains to study the other factor of the polynomial

ψ 0 (Z) = [1 + δ + ωε ′ s ]Z 2 -[2 + 2ω(ε ′ s -1)]Z + [1 -δ + ωε ′ s ]
which has not ±i as roots. We notice that |ψ 0 (0)| < |ψ * 0 (0)| and compute

ψ 1 (Z) = 4δ{[1 + ωε ′ s ]Z -[1 + ω(ε ′ s -1)]},
which is a Schur polynomial for all ε ′ s . Polynomial ψ 0 is therefore a Schur polynomial and φ 0 is a simple von Neumann polynomial.

Case ε

′ s = 1 If ε ′ s = 1, the polynomial φ 3 is identically zero and φ 2 (Z) = 16δ 2 ω(2 + ω){Z 2 -[2 -q]Z + 1}, φ ′ 2 (Z) = 16δ 2 ω(2 + ω){2Z -[2 -q]}.
The root of φ ′ 2 does have a < 1 modulus if q ∈]0, 2]. The polynomial φ 0 is therefore a simple von Neumann polynomial.

Synthesis for the Lorentz-Joseph et al. Model

The scheme (8)-(13) for the one-dimensional anharmonic Maxwell-Lorentz equations is stable with the condition δt ≤ δx/ √ 2c ∞ .

Kashiwa et al. Model

Model Setting

A modified version of Maxwell system ( 9) is closed by a discretization of Lorentz equation ( 6), namely

1 δt (B n+ 1 2 j+ 1 2 -B n-1 2 j+ 1 2 ) = - 1 δx (E n j+1 -E n j ), ε 0 ε ∞ δt (E n+1 j -E n j ) = - 1 µ 0 δx (B n+ 1 2 j+ 1 2 -B n+ 1 2 j-1 2 ) - 1 δt (P n+1 j -P n j ), 1 δt (P n+1 j -P n j ) = 1 2 (J n+1 j + J n j ), 1 δt (J n+1 j -J n j ) = - ν 2 (J n+1 j + J n j ) + ω 2 1 (ε s -ε ∞ )ε 0 2 (E n+1 j + E n j ) - ω 2 1 2 (P n+1 j + P n j ). ( 14
)
The system (14) deals with the variable

U n j = (c ∞ B n-1 2 j+ 1 2 , E n j , P n j /ε 0 ε ∞ , δtJ n j /ε 0 ε ∞ ) t = (B n-1 2 j+ 1 2 , E n j , P n j , J n j ) t
and reads

B n+ 1 2 j+ 1 2 = B n-1 2 j+ 1 2 -λ(E n j+1 -E n j ), [1 + δ + 1 2 ωε ′ s ]E n+1 j = [1 + δ - 1 2 ω(ε ′ s -2)]E n j + λ 2 (1 + δ + 1 2 ω)(E n j+1 -2E n j + E n j-1 ) -λ(1 + δ + 1 2 ω)(B n-1 2 j+ 1 2 -B n-1 2 j-1 2 ) + ωP n j -J n j , [1 + δ + 1 2 ωε ′ s ]P n+1 j = [1 + δ + 1 2 ω(ε ′ s -2)]P n j - 1 2 λω(ε ′ s -1)(B n-1 2 j+ 1 2 -B n-1 2 j-1 2 ) +ω(ε ′ s -1)E n j + 1 2 λ 2 ω(ε ′ s -1)(E n j+1 -2E n j + E n j-1 ) + J n j , [1 + δ + 1 2 ωε ′ s ]J n+1 j = [1 -δ - 1 2 ωε ′ s ]J n j -λω(ε ′ s -1)(B n-1 2 j+ 1 2 -B n-1 2 j-1 2 ) +2ω(ε ′ s -1)E n j + λ 2 ω(ε ′ s -1)(E n j+1 -2E n j + E n j-1 ) -2ωP n j ,
from which stems the amplification matrix

G =       1 -σ 0 0 σ * (D-1 2 ω(ε ′ s -1)) D (1-q)D-(2-q) 1 2 ω(ε ′ s -1) D ω D -1 D σ * 1 2 ω(ε ′ s -1) D (2-q) 1 2 ω(ε ′ s -1) D D-ω D 1 D σ * ω(ε ′ s -1) D (2-q)ω(ε ′ s -1) D -2ω D 2-D D      
where, along with earlier notations,

D = 1 + δ + 1 2 ωε ′ s .

Computation of the Characteristic Polynomial

The characteristic polynomial of G is

P (Z) = Z -1 σ 0 0 -σ * (D-1 2 ω(ε ′ s -1)) D Z - (1-q)D-(2-q) 1 2 ω(ε ′ s -1) D -ω D 1 D -σ * 1 2 ω(ε ′ s -1) D - (2-q) 1 2 ω(ε ′ s -1) D Z -D-ω D -1 D -σ * ω(ε ′ s -1) D -(2-q)ω(ε ′ s -1) D 2ω D Z -2-D D .
hence setting X = D(Z -1)

D 4 P (Z) = X Dσ 0 0 -σ * (D -1 2 ω(ε ′ s -1)) X + qD + (2 -q) 1 2 ω(ε ′ s -1) -ω 1 -σ * 1 2 ω(ε ′ s -1) -(2 -q) 1 2 ω(ε ′ s -1) X + ω -1 -σ * ω(ε ′ s -1) -(2 -q)ω(ε ′ s -1) 2ω X + 2(D -1) = X Dσ 0 0 -σ * D X + qD X 0 -σ * 1 2 ω(ε ′ s -1) -(2 -q) 1 2 ω(ε ′ s -1) X + ω -1 0 0 -2X X + 2D = X X + qD X 0 -(2 -q) 1 2 ω(ε ′ s -1) X + ω -1 0 -2X X + 2D + qD D X 0 1 2 ω(ε ′ s -1) X + ω -1 0 -2X X + 2D = X{(X + qD)(X + ω)(X + 2D) + (2 -q) 1 2 ω(ε ′ s -1)X(X + 2D) -2X(X + qD)} +qD{D(X + ω)(X + 2D) -2DX - 1 2 ω(ε ′ s -1)X(X + 2D)} = X 4 + [2D -2 + ωε ′ s + (D - 1 2 ω(ε ′ s -1))q]X 3 + D[2ωε ′ s + (3D -2 + 5 2 ω - 3 2 ωε ′ s )q]X 2 +D 2 [(2D -2 + 4ω -ωε ′ s )q]X + D 3 [2ωq] = X 4 + [2(δ + ωε ′ s ) + (1 + δ + 1 2 ω)q]X 3 + D[2ωε ′ s + (1 + 3δ + 5 2 ω)q]X 2 +D 2 [(2δ + 4ω)q]X + D 3 [2ωq].
The characteristic polynomial is proportional to

φ 0 (Z) = [1 + δ + 1 2 ωε ′ s ]Z 4 -[4 + 2δ -(1 + δ + 1 2 ω)q]Z 3 + [6 -ωε ′ s + (ω -2)q]Z 2 -[4 -2δ -(1 -δ + 1 2 ω)q]Z + [1 -δ + 1 2 ωε ′ s ].

Von Neumann Analysis

We successively compute

φ 1 (Z) = 2δ{[2 + ωε ′ s ]Z 3 -[6 + ωε ′ s -(2 + 1 2 ω(ε ′ s + 1))q]Z 2 + [6 -ωε ′ s -(2 -ω)q]Z -[2 -ωε ′ s + 1 2 ω(ε ′ s -1)q]}, φ 2 (Z) = 4δ 2 ω{[8ε ′ s -(ε ′ s -1)(2 -ωε ′ s )q - 1 4 ω(ε ′ s -1) 2 q 2 ]Z 2 -[16ε ′ s -8ε ′ s q + (ε ′ s -1)(1 - 1 2 ω)q 2 ]Z +[8ε ′ s -(ε ′ s -1)(6 + ωε ′ s )q -(ε ′ s -1) 2 (1 - 1 4 ω)q 2 ]}, φ 3 (Z) = 4δ 4 ω 2 (ε ′ s -1)q(4 -q)(2 + ωε ′ s ) × ×{[32ε ′ s -16(ε ′ s -1)q + (ε ′ s -1)(2 + ω)q 2 ]Z -[32ε ′ s -16ε ′ s q + (ε ′ s -1)(2 -ω)q 2 ]}.
We see that the specific cases q = 0 and ε ′ s = 1 which make φ 3 vanish shall be treated separately. The general case is treated by first checking that |φ 0 (0)| < |φ * 0 (0)|, which is obvious. We then notice that φ * 1 (0) > 0. The relation φ 1 (0) < φ * 1 (0) is equivalent to -ω(ε ′ s -1)q < 8, which always holds. As for relation -φ 1 (0) < φ * 1 (0), it can be cast as ω(ε ′ s -1)q < 4ωε ′ s , which holds true if q ≤ 4. We therefore have |φ 1 (0)| < |φ * 1 (0)| for q ∈]0, 4]. We carry on by studying the sign of

φ * 2 (0) = δ 2 ω[4ε ′ s -(ε ′ s -1)q][8 + (ε ′ s -1)ωq] > 0 and therefore we have to check that φ * 2 (0) + φ 2 (0) > 0 and φ * 2 (0) -φ 2 (0) > 0 φ * 2 (0) + φ 2 (0) = 2δ 2 ω[(ε ′ s -1)ωq 2 + 32 + 2(ε ′ s -1)(4 -q) 2 ] > 0, φ * 2 (0) -φ 2 (0) = 2δ 2 qω(ε ′ s -1)(2 + ε ′ s ω)(4 -q) > 0, if q ∈]0, 4[. Last let us study φ 3 φ * 3 (0) = 4δ 4 ω 2 (ε ′ s -1)q(4 -q)(2 + ε ′ s ω)[32 + (ε ′ s -1)(2(4 -q) + ωq 2 )] > 0, φ * 3 (0) + φ 3 (0) = 8δ 4 ω 2 (ε ′ s -1)q 2 (4 -q)(2 + ε ′ s ω)[8 + ω(ε ′ s -1)q] > 0, φ * 3 (0) -φ 3 (0) = 8δ 4 ω 2 (ε ′ s -1)q(4 -q) 2 (2 + ε ′ s ω)[4ε ′ s -(ε ′ s -1)q] > 0.
Hence we show that φ 3 and therefore φ 0 is a Schur polynomial if q ∈]0, 4[ and there remains to treat the specific cases.

Case q = 0

The case q = 0 has once more to be treated separately. We have

φ 0 (Z) = (Z -1) 2 (Z -1) 2 + δ(Z 2 -1) + 1 2 ωε ′ s (Z + 1) 2 .
The corresponding amplification matrix is

     1 0 0 0 0 D-ω(ε ′ s -1) D ω D -1 D 0 ω(ε ′ s -1) D D-ω D 1 D 0 2ω(ε ′ s -1) D -2ω D 2-D D     
.

Anew 1 is a double eigenvalue in two distinct eigensubspaces. The other factor of the characteristic polynomial, namely

ψ 0 (Z) = [1 + δ + 1 2 ωε ′ s ]Z 2 -[2 -ωε ′ s ]Z + [1 -δ + 1 2 ωε ′ s ],
should be a Schur (or a simple von Neumann) polynomial. We do have |ψ 0 (0)| < |ψ * 0 (0)| and we compute

ψ 1 (Z) = 4δ{[1 + 1 2 ωε ′ s ]Z -[1 - 1 2 ωε ′ s ]}.
Both remaining eigenvalues have a strictly lower to 1 modulus and iterates of the amplification matrix are bounded. This holds even if ε ′ s = 1.

Case q = 4

In the case when q = 4,

φ 0 (Z) = [1 + δ + 1 2 ωε ′ s ]Z 4 + [2δ + 2ω]Z 3 + [-2 -ωε ′ s + 4ω]Z 2 + [-2δ + 2ω]Z + [1 -δ + 1 2 ωε ′ s ] = (Z + 1) 2 {[1 + δ + 1 2 ωε ′ s ]Z 2 -2[1 -ω + 1 2 ωε ′ s ]Z + [1 -δ + 1 2 ωε ′ s ]}.
We have a double root Z = -1. We therefore have to study the amplification matrix which reads

G =       1 -σ 0 0 σ * (D-1 2 ω(ε ′ s -1)) D -3D+ω(ε ′ s -1) D ω D -1 D σ * 1 2 ω(ε ′ s -1) D -2 1 2 ω(ε ′ s -1) D D-ω D 1 D σ * ω(ε ′ s -1) D -2ω(ε ′ s -1) D -2ω D 2-D D       .
Only the vector (σ, 2, 0, 0) t is an eigenvector associated to the eigenvalue -1, and we have increasing iterates for G, whatever the study of the other factor of the characteristic polynomial. The value q = 4 gives rise to instabilities.

Case ε

′ s = 1 If ε ′ s = 1
, the polynomial φ 3 is identically zero and we shall study φ 2 for q ∈]0, 4[

φ 2 (Z) = 32δ 2 ω{Z 2 -[2 -q]Z + 1}.
This polynomial has two distinct complex conjugate roots with unit modulus. Polynomials φ 2 and therefore φ 0 are both simple von Neumann polynomials.

Synthesis for the Lorentz-Kashiwa Model

The scheme (14) for the one-dimensional anharmonic Maxwell-Lorentz equations is stable with the condition δt < δx/c ∞ .

Young Model

Model Setting

The Maxwell system ( 9) is closed by a discretization of ( 6), namely

1 δt (P n+1 j -P n j ) = J n+ 1 2 , 1 δt (J n+ 1 2 j -J n-1 2 j ) = - ν 2 (J n+ 1 2 j + J n-1 2 j ) + ω 2 1 (ε s -ε ∞ )ε 0 E n j -ω 2 1 P n j . (15)
The explicit version of system ( 9)-( 15) deals once more with the variable

U n j = (c ∞ B n-1 2 j+ 1 2 , E n j , P n j /ε 0 ε ∞ , δtJ n-1 2 j /ε 0 ε ∞ ) t = (B n-1 2 j+ 1 2 , E n j , P n j , J n-1 2 j ) t
and reads

B n+ 1 2 j+ 1 2 = B n-1 2 j+ 1 2 -λ(E n j+1 -E n j ), [1 + δ]E n+1 j = [1 + δ -2ω(ε ′ s -1)]E n j + λ 2 [1 + δ](E n j+1 -2E n j + E n j-1 ) -λ(B n-1 2 j+ 1 2 -B n-1 2 j-1 2 ) +2ωP n j -[1 -δ]J n-1 2 j , [1 + δ]P n+1 j = [1 + δ -2ω]P n j + 2ω(ε ′ s -1)E n j + [1 -δ]J n-1 2 j , [1 + δ]J n+ 1 2 j = [1 -δ]J n-1 2 j + 2ω(ε ′ s -1)E n j -2ωP n j ,
from which stems the amplification matrix

G =      1 -σ 0 0 σ * (1-q)(1+δ)-2ωα 1+δ 2ω 1+δ -1-δ 1+δ 0 2ωα 1+δ 1+δ-2ω 1+δ 1-δ 1+δ 0 2ωα 1+δ -2ω 1+δ 1-δ 1+δ      .

Computation of the Characteristic Polynomial

The characteristic polynomial of G is

P (Z) = Z -1 σ 0 0 -σ * Z -(1-q)(1+δ)-2ωα 1+δ -2ω 1+δ 1-δ 1+δ 0 -2ωα 1+δ Z -1+δ-2ω 1+δ -1-δ 1+δ 0 -2ωα 1+δ 2ω 1+δ Z -1-δ 1+δ .
Therefore setting X = (1 + δ)(Z -1),

(1 + δ) 4 P (Z) = X (1 + δ)σ 0 0 -σ * (1 + δ) X + (1 + δ)q + 2ωα -2ω 1 -δ 0 -2ωα X + 2ω -(1 -δ) 0 -2ωα 2ω X + 2δ = X (1 + δ)σ 0 0 -σ * (1 + δ) X + (1 + δ)q X 0 0 -2ωα X + 2ω -(1 -δ) 0 0 -X X + 1 + δ = X X + (1 + δ)q X 0 -2ωα X + 2ω -(1 -δ) 0 -X X + 1 + δ +(1 + δ) 2 q 1 0 0 -2ωα X + 2ω -(1 -δ) 0 -X X + 1 + δ = X{X 3 + [2δ + 2ωε ′ s + (1 + δ)q]X 2 + [2(1 + δ)ωε ′ s + 2(1 + δ)(δ + ω)q]X +[2(1 + δ) 2 ωq]} + (1 + δ) 2 q{X 2 + [2(δ + ω)]X + [2(1 + δ)ω]} = X 4 + X 3 [2δ + 2ωε ′ s + (1 + δ)q] + X 2 [(1 + δ)(2ωε ′ s + (1 + 3δ + 2ω)q)] +X[(1 + δ) 2 q(2δ + 4ω)] + [2(1 + δ) 3 ωq].
The characteristic polynomial is proportional to

φ 0 (Z) = [1 + δ]Z 4 -[4 + 2δ -2ωε ′ s -(1 + δ)q]Z 3 + 2[3 -2ωε ′ s + (ω -1)q]Z 2 -[4 -2δ -2ωε ′ s -(1 -δ)q]Z + [1 -δ].

Von Neumann Analysis

We successively compute

φ 1 (Z) = 4δ{Z 3 -[3 -ωε ′ s -q]Z 2 + [3 -2ωε ′ s + (ω -1)q]Z -[1 -ωε ′ s ]}, φ 2 (Z) = (4δ) 2 ω{[ε ′ s (2 -ωε ′ s )]Z 2 -[2ε ′ s (2 -ωε ′ s ) + (ε ′ s (ω -1) -1))q]Z +[ε ′ s (2 -ωε ′ s ) -(ε ′ s -1)q]}, φ 3 (Z) = (4δ) 4 qω 2 (ε ′ s -1){[2ε ′ s (2 -ωε ′ s ) -(ε ′ s -1)q]Z -[2ε ′ s (2 -ωε ′ s ) -(ε ′ s -1)q -(2 -ωε ′ s )q]}.
We see that we shall once more treat the cases q = 0 and ε ′ s = 1 separately since φ 3 is identically zero. Let us check the conditions in the general case. First, |φ 0 (0)| < |φ * 0 (0)| clearly holds as well as

|φ 1 (0)| < |φ * 1 (0)| under the condition ω < 2/ε ′ s . If ε ′ s > 1 and q = 0, the condition |φ 2 (0)| < |φ * 2 (0)| is equivalent to (ε ′ s -1)q < 2ε ′ s (2 -ωε ′ s ). If the worst case is q = 2, we must have (ε ′ s -1)2 < 2ε ′ s (2-ωε ′ s ), which is equivalent to ω < (ε ′ s +1)/ε ′ s 2 .
If the worst case is q = 4, we must have (ε

′ s -1)4 < 2ε ′ s (2 -ωε ′ s ), which is equivalent to ω < 2/ε ′ s 2 .
We wait until the study of φ 3 to choose between q ≤ 2 and q ≤ 4. The root of φ 3 is

Z = 2ε ′ s (2 -ωε ′ s ) -(ε ′ s -1)q -(2 -ωε ′ s )q 2ε ′ s (2 -ωε ′ s ) -(ε ′ s -1)q .
The denominator is positive under the same assumption found to ensure

|φ 2 (0)| < |φ * 2 (0)|. If we want |Z| < 1, the condition is hence (2 -ωε ′ s )q < 4ε ′ s (2 -ωε ′ s ) -2(ε ′ s -1)q. If the worst case is q = 2, we must have 2(2 -ωε ′ s ) < 4ε ′ s (2 -ωε ′ s ) -4(ε ′ s -1)
, which is equivalent to ω < 2/(2ε ′ s -1). If the worst case is q = 4, we must have 4(2

-ωε ′ s ) < 4ε ′ s (2 -ωε ′ s ) -8(ε ′ s -1)
, which is equivalent to 4ωε ′ s (ε ′ s -1) < 0, which is false. We therefore choose to take q ≤ 2 and the successive conditions found are

ω < 2 ε ′ s , ω < ε ′ s + 1 ε ′ s 2 , ω < 2 2ε ′ s -1 .
The more restrictive condition that we have encountered is ω < 2/(2ε ′ s -1), this is therefore our final stability condition in addition to q < 2, for which φ 0 is a Schur polynomial.

Which are the limiting case we have to study? The three conditions are equivalent if and only if ε ′ s = 1. In this case, if ω has its limit value ω = 2 and q its limit value q = 2, we have φ 2 (Z) ≡ 0. If ε ′ s = 1, q = 2 and ω = 2/(2ε ′ s -1), then the modulus of the root of φ 3 is 1, and we conclude that φ 3 and therefore φ 0 are simple von Neumann polynomials.

Case q = 0

If q = 0, the characteristic polynomial has the double eigenvalue 1

φ 0 (Z) = (Z -1) 2 ((Z -1) 2 + δ(Z 2 -1) + 2ωε ′ s Z). The corresponding amplification matrix is G =     1 0 0 0 0 1+δ-2ωα 1+δ 2ω 1+δ -1-δ 1+δ 0 2ωα 1+δ 1+δ-2ω 1+δ 1-δ 1+δ 0 2ωα 1+δ -2ω 1+δ 1-δ 1+δ     .
The double eigenvalue operates on two distinct eigensubspaces. We have to study the other factor of the polynomial

ψ 0 (Z) = [1 + δ]Z 2 -[2(1 -ωε ′ s )]Z + [1 -δ]. We clearly have |ψ 0 (0)| < |ψ * 0 (0)|. Moreover we compute ψ 1 (Z) = 4δ{Z -[1 -ωε ′ s ]}.
We recover the condition ω < 2/ε ′ s , under which we have a Schur polynomial. For ω = 2/ε ′ s , we have a simple von Neumann polynomial (-1 is a root), which allows to conclude.

Case ε ′

s = 1 In the case when ε ′ s = 1, only the lasts steps have to be considered. The only problems are the condition |φ 2 (0)| < |φ * 2 (0)| and a vanishing φ 3 . In this specific case,

φ 2 (Z) = (4δ) 2 ω(2 -ω){Z 2 -[2 -q]Z + 1}.
If ω < 2/ε ′ s , i.e. ω < 2, this polynomial is identically zero, and

φ ′ 2 (Z) = (4δ) 2 ω(2 -ω){2Z -[2 -q]}. Polynomial φ ′
2 is a Schur one if q ∈]0, 2] and hence φ 0 is a simple von Neumann polynomial. If ω = 2, polynomial φ 2 is identically zero and we compute

φ 1 (Z) = 4δ(Z + 1){Z 2 -[2 -q]Z + 1}.
For q ∈]0, 2] the roots of Z 2 -[2 -q]Z + 1 are complex conjugate, distinct, and their modulus is 1.

Synthesis for the Lorentz-Young model

The scheme ( 9)-(15) for the one-dimensional anharmonic Maxwell-Lorentz is stable with the condition

δt ≤ min δx √ 2c ∞ , 2 ω 1 2ε ′ s -1 .

Harmonic Lorentz Type Media

Harmonic Lorentz type media are treated thanks to the three above mentioned schemes. The computation of the amplification matrices and the characteristic polynomials remains unchanged. The harmonicity ν = 0 is expressed by the parameter δ = 0. This vanishing value makes φ 1 identically zero for the three schemes and the above analysis breaks down. We resume to the analysis of the three schemes.

Joseph et al. Model

General Case

Since φ 1 is identically zero, we want to apply Theorem 2 and study the derivative polynomial of φ 0 , which we denote by

ψ 0 (Z) = [4 + 4ωε ′ s ]Z 3 -[12 + 6ωε ′ s -3(1 + ω)q]Z 2 + [12 + 4ωε ′ s -4q]Z -[4 + 2ωε ′ s -(1 + ω)q].
We notice that ψ * 0 (0) > 0 and ψ 0 (0) = -(4 -q) -ω(2ε ′ s -q) < 0 for q ≤ 2. We therefore have to check that -ψ 0 (0) < ψ * 0 (0), which is equivalent to -(ω + 1)q < 2ωε ′ s and always holds. We therefore have

|ψ 0 (0)| < |ψ * 0 (0)|. Let us now compute ψ 1 (Z) = [4ωε ′ s (4 + 3ωε ′ s ) + 4(1 + ω)(2 + ωε ′ s )q -(1 + ω) 2 q 2 ]Z 2 +[-16ωε ′ s (2 + ωε ′ s ) + 8(ω 2 ε ′ s -2)q + 4(1 + ω)q 2 ]Z +[4ωε ′ s (4 + ωε ′ s ) + 4(2 -ωε ′ s + 6ω + 3ω 2 ε ′ s )q -3(1 + ω) 2 q 2 ].
Anew we have

ψ * 1 (0) = 4ωε ′ s (4 + 3ωε ′ s ) + (1 + ω)q[4(2 + ωε ′ s ) -(1 + ω)q] > 0.
Instead of studying the sign of ψ 1 (0), we will check that -ψ 1 (0) < ψ * 1 (0) and ψ 1 (0

) < ψ * 1 (0). The relation -ψ 1 (0) < ψ * 1 (0) is equivalent to 4ωε ′ s (2 + ωε ′ s ) + (4 -q)(1 + ω) 2 q + 4ω 2 (ε ′ s -1)q > 0,
which clearly holds. There remains ψ 1 (0) < ψ * 1 (0) which reverts to

4ω 2 ε ′ s 2 + 4ω(ε ′ s -2 -ωε ′ s ) + (1 + ω) 2 q 2 > 0.
Cast like this it is not easy to conclude, but we can write it has a polynomial of the variable ω

ω 2 (2ε ′ s -q) 2 + 2ωq(2(ε ′ s -2) + q) + q 2 > 0.
The reduced discriminant of this polynomial is

∆ ′ = -8(ε ′ s -1)q 2 (2 -q) < 0, if 0 < q < 2 and ε ′ s = 1.
The polynomial (in ω) is therefore always positive, which we were seeking. Hence we have |ψ 1 (0)| < |ψ * 1 (0)| and we can profitably carry on with the computation of ψ 2 (Z) which is the product

ψ 2 (Z) = 8[4ω 2 ε ′ s 2 + (4ωε ′ s -4ω 2 ε ′ s -8ω)q + (1 + ω) 2 q 2 ] × ×{[4ω 2 ε ′ s 2 + 8ωε ′ s + (4 + 4ωε ′ s + 8ω)q -(1 + ω) 2 q 2 ]Z -[4ω 2 ε ′ s 2 + 8ωε ′ s + (4 -2ωε ′ s )q -(1 + ω)q 2 ]} We notice that 8[4ω 2 ε ′ s 2 + (4ωε ′ s -4ω 2 ε ′ s -8ω)q + (1 + ω) 2 q 2 ] = 8[(4ε ′ s -q)ω -q] 2 + 16ωq[4(ε ′ s -1) + 2ε ′ s ] > 0
and we simplify by this factor denoting

ψ2 (Z) = [4ω 2 ε ′ s 2 + 8ωε ′ s + (4 + 4ω 2 ε ′ s + 8ω)q -(1 + ω) 2 q 2 ]Z -[4ω 2 ε ′ s 2 + 8ωε ′ s + (4 -2ωε ′ s )q -(1 + ω)q 2 ]. We see that ψ * 2 (0) = 4ω 2 [ε ′ s 2 + (ε ′ s -1)q] + (1 + ω) 2 q(4 -q) ≥ 0 and to prove |ψ 2 (0)| < |ψ * 2 (0)|, we only have to check that ψ * 2 (0) + ψ2 (0) > 0 and ψ * 2 (0) -ψ2 (0) > 0. We "notice" that ψ * 2 (0) + ψ2 (0) = ωq[(4ε ′ s -q)ω + 2(8 -q)] > 0, ψ * 2 (0) -ψ2 (0) = [2ωε ′ s + (1 + ω)q][(4ε ′ s -q)ω + 2(4 -q)]
> 0, which ends the proof in the general case for δ = 0.

Case q = 0

In the case when q = 0, along with the fact that 1 is a double root "which does not cause any trouble", φ 0 has the same roots as those of [1 + ωε ′ s ]Z 2 -2Z + [1 + ωε ′ s ] which are complex conjugate, distinct, and their modulus is 1, if ε ′ s ≥ 1.

Case q = 2

The same holds for q = 2 and this time, along with the roots ±i, we have the roots of the polynomial

[1+ωε ′ s ]Z 2 -2[1+ω(ε ′ s -1)]Z +[1+ωε ′ s ]
which has two complex conjugate, distinct roots with modulus 1, if ε ′ s ≥ 1.

Case ε

′ s = 1 Finally if ε ′ s = 1 (and q ∈]0, 2[),
we shall return to the polynomial φ 0 which can be cast as

φ 0 (Z) = [Z 2 -(2 -q)Z + 1][(1 + ω)Z 2 -2Z + (1 + ω)].
Each of the second degree polynomials has two distinct complex conjugate roots. We therefore have a simple von Neumann polynomial except in the particular case when the two polynomials are proportional and have the same roots, which are then double roots. This is reached if (2 -q) = 2/(1 + ω), namely q = 2ω/(1 + ω) or equivalently ω = q/(2 -q). In this case, von Neumann analysis is not useful anymore and we have to revert to the amplification matrix

G =     1 -σ 0 0 0 -2q + 2 -1 q 0 1 0 0 σ * -q 0 1     .
The two double eigenvalues of this matrix are (2 -q ± i q(4 -q))/2 and their each only have one corresponding eigenvector σ, q ∓ i q(4 -q) 2 , q(3 -q) ∓ i(2 -q) q(4 -q) 2 , -q ± i q(4 -q) 2 t .

For each eigenvalue the associated minimal eigensubspace is therefore two-dimensional, which corresponds to an unstable case. We can say that the scheme is stable for q ∈ [0, 2ω/(1 + ω)[. If we rewrite this in physical variables, we have

4c 2 ∞ δt 2 δx 2 sin 2 ( ξ 2 ) < 2ω 2 1 δt 2 2 + ω 2 1 δt 2 . If δx 2 < 4c 2 ∞ /ω 2 1
, this is not a bound on the time step. If δx if large enough, the bound on the time step is

δt 2 < δx 2 2c 2 ∞ - 2 ω 2 1 .

Synthesis for the Harmonic Lorentz-Joseph et al. Model

The scheme ( 8)-( 13) for the one-dimensional harmonic Maxwell-Lorentz equations is stable with the condition

δt ≤ δx √ 2c ∞ if ε s > ε ∞ and δt < δx 2 2c 2 ∞ - 2 ω 2 1 if ε s = ε ∞ , this last condition being meaningful only if δx > 2c ∞ /ω 1 .
It is therefore advisable not to use the Lorentz-Joseph et al. scheme in the harmonic case for ε s = ε ∞ .

Kashiwa Model

General Case

Anew the polynomial φ 1 is identically zero and we study the derivative of polynomial φ 0 , which we denote

ψ 0 (Z) = [4 + 2ωε ′ s ]Z 3 -[12 -3(1 + 1 2 ω)q]Z 2 + [12 -2ωε ′ s -(4 -2ω)q]Z -[4 -(1 + 1 2 ω)q].
The condition |ψ 0 (0)| < |ψ * 0 (0)| is equivalent to (8 -q) + 1 2 ω(4ε ′ s -q) > 0 which holds for q ≤ 4. Then we compute

ψ 1 (Z) = [4(4 + ωε ′ s )ωε ′ s + 4(2 + ω)q -(1 + 1 2 ω) 2 q 2 ]Z 2 -[32ωε ′ s + (16 -8ω -8ωε ′ s -4ω 2 ε ′ s )q + (ω 2 -4)q 2 ]Z +[4(4 -ωε ′ s )ωε ′ s + 4(2 -2ωε ′ s + 5ω)q -3(1 + 1 2 ω) 2 q 2 ].
We check that

ψ * 1 (0) = 1 4 [4ωε ′ s + (2 + ω)q][16 + 4ωε ′ s -(2 + ω)q] > 0, ψ * 1 (0) + ψ 1 (0) = q(4 -q)(ω + 2) 2 + 4q(ε ′ s -1)ω 2 + 8ε ′ s (4 -q) + 8q > 0, ψ * 1 (0) -ψ 1 (0) = 1 2 {[(4ε ′ s -q)ω -2q] 2 + 32qω(ε ′ s -1)} > 0,
under the only condition that q ∈ [0, 4], which ensures that |ψ 1 (0)| < |ψ * 1 (0)|. Last we compute ψ 2 (Z) which can be cast as

ψ 2 (Z) = 1 2 {[(-8 + 4ε ′ s + q)ω + 2q] 2 + 16α(4 -q)} ψ2 (Z) with ψ2 (Z) = [(2 + ω) 2 q(4 -q) + 4((ε ′ s -1)ω 2 q + 2(ε ′ s -1)ω(4 -q) + 8)]Z -[8 -(2 + ω)q][4ωε ′ s + (2 -ω)q].
We check that

ψ * 2 (0) ≥ 0, ψ * 2 (0) + ψ2 (0) = qω[(6ε ′ s -q)ω + 2(8 -q)] > 0, ψ * 2 (0) -ψ2 (0) = 2(4 -q)[4ε ′ s ω + (2 + ω)q] > 0.
These two last inequalities are strict only if q ∈]0, 4[, and we then have

|ψ 2 (0)| < |ψ * 2 (0)|.
The case ε ′ s = 1 is not specific in this general study.

Case q = 0

To treat the specific case when q = 0, we have to revert to the study of φ 0 which is here

φ 0 (Z) = [1 + 1 2 ωε ′ s ]Z 4 -4Z 3 + [6 -ωε ′ s ]Z 2 -4Z + [1 + 1 2 ωε ′ s ] = (Z -1) 2 {[1 + 1 2 ωε ′ s ]Z 2 -2[1 - 1 2 ωε ′ s ]Z + [1 + 1 2 ωε ′ s ]}.
The double eigenvalue Z = 1 is the same as in the anharmonic case and does not cause any trouble either (minimal eigensubspaces are still one-dimensional). The other factor of the polynomial has clearly two distinct complex conjugate roots of modulus 1. This configuration corresponds to a stability case for the scheme.

Case q = 4

The analysis performed in the anharmonic case remains valid here. The eigenvalue Z = -1 is double and the associated minimal eigensubspace is two-dimensional. Iterates of the amplification matrix are therefore linearly increasing and the case is unstable.

Synthesis for the Harmonic Lorentz-Kashiwa Model

The scheme (14) for one-dimensional harmonic Maxwell-Lorentz equations is stable with the condition δt < δx/c ∞ .

Young Model

general Case

Once more, the polynomial φ 1 is identically zero and we study the derivative of polynomial φ 0 , which we denote

ψ 0 (Z) = 4Z 3 + [-12 + 6ωε ′ s + 3q]Z 2 + [12 -8ωε ′ s + 4(ω -1)q]Z + [-4 + 2ωε ′ s + q].
The condition |ψ 0 (0)| < |ψ * 0 (0)| is equivalent to (4 -2ωε ′ s ) + (4 -q) > 0, which we assume (ω < 2/ε ′ s and q ∈ [0, 2]). We carry on computing

ψ 1 (Z) = [(2ωε ′ s + q)(8 -(2ωε ′ s + q))]Z 2 + [(2ωε ′ s + q)(4(2ωε ′ s + q) -16 -4ωq) + 16ωq]Z +[(2ωε ′ s + q)(8 -3(2ωε ′ s + q)) + 16ωq].
We see immediately in this formulation that ψ * 1 (0) > 0. Besides

ψ * 1 (0) + ψ 1 (0) = 4{-4ω 2 ε ′ s 2 + 8ωε ′ s + 4(ω(ε ′ s + 1) + 1)q -q 2 }, ψ * 1 (0) -ψ 1 (0) = 2{[2ω -q] 2 + 4ω(ε ′ s -1)[q + ω(ε ′ s + 1)}] > 0, Let us note f (ω) = -4ω 2 ε ′ s 2 + 8ωε ′ s + 4(ω(ε ′ s + 1) + 1)q -q 2 ,
we want to prove that this quantity is positive for ω ∈]0, 2/(2ε ′ s -1)]. We derive to obtain

f ′ (ω) = -4ε ′ s q -8ωε ′ s 2 + 8ε ′ s + 4q which vanishes at ω = (2ε ′ s -(ε ′ s -1)q)/2ε ′ s 2 .
This corresponds to values of ω between 1/ε ′ s 2 (value for q = 2) and 1/ε ′ s (value for q = 0), which always belong to the interval ]0, 2/(2ε ′ s -1)]. At this point we have a maximum of the function f (ω). To conclude, we only have to evaluate the limit for f (ω) as ω → 0 and the value of f (ω) at ω = 2/(2ε ′ s -1). If both values are positive, f (ω) will be positive on the whole interval.

lim ω→0 f (ω) = 4q -q 2 > 0, f (2/(2ε ′ s -1)) = (2ε ′ s -1) 2 (2 -q) + 2(2ε ′ s -1) + 4(ε ′ s -1)(2ε ′ s + 1) > 0, If ε ′ s = 1, q = 2
and ω = 2, we have ψ * 1 (0) = -ψ 1 (0). We will treat this case apart. We finally compute ψ 2 (Z) which can be cast as

ψ 2 (Z) = 8[(2ωε ′ s -q) 2 + 8ω(ε ′ s -1)q] ψ2 (Z) with ψ2 (Z) = [-4ω 2 ε ′ s 2 + 8ωε ′ s + 4(ω(ε ′ s + 1) + 1)q -q 2 ]Z -[(2ωε ′ s + (1 -ω)q)(-4 + 2ωε ′ s + q)].
We notice that ψ 2 is identically zero if ε ′ s = 1 and q = 2ωε ′ s = 2ω, which has to be treated separately. In the opposite case, we check that

ψ * 2 (0) ≥ f (ω) > 0, ψ * 2 (0) + ψ2 (0) = qω[8 -2ωε ′ s -q] > 0, ψ * 2 (0) -ψ2 (0) = [2ωε ′ s + q][8 -4ωε ′ s + (ω -2)q].
The quantity 8 -4ωε ′ s + (ω -2)q is minimum if q = 2 (since ω < 2) and is then equal to 4 -2ω(2ε ′ s -1). As in the anharmonic case if ω < 2/(2ε ′ s -1) this quantity is positive, and if ω = 2/(2ε ′ s -1) this quantity is zero. Yet we want to show that ψ 0 is a Schur or a simple von Neumann polynomial, we therefore have to revert to the study of φ 0 . Once more the cases ε ′ s = 1, q = 2 and ω = 2 have to be treated specifically.

Case q = 0

In the case when q = 0, we revert to the direct study of φ 0

φ 0 (Z) = Z 4 -[4 -2ωε ′ s ]Z 3 + 2[3 -2ωε ′ s ]Z 2 -[4 -2ωε ′ s ]Z + 1 = (Z -1) 2 [Z 2 -2(1 -ωε ′ s )Z + 1],
which has Z = 1 as a double root, which is no more a problem as in the anharmonic case. Both other roots are complex conjugate, distinct and have a unit modulus if ωε ′ s < 2. If ε ′ s = 1 and ω = 2, Z = -1 is also a double root. Then we have

G + Id =     2 -σ 0 0 σ * 2 1 -1 2 0 0 1 1 1 0 0 -2 1     ,
which has only (σ, 2, 0, 0) t as eigenvalue (associated to 0). This is a cause of instability for the scheme.

Case q = 2

The case q = 2 is treated by the general case except when ω = 2/(2ε ′ s -1). In this case 2ωε ′ s = 2 + ω, hence

φ 0 (Z) = Z 4 + ωZ 3 + 2[ω -1]Z 2 + ωZ + 1 = (Z + 1) 2 (Z 2 -(2 -ω)Z + 1).
We therefore have to study the stable subspaces through the amplification matrix for the eigenvalue -1. We have

G + Id =     2 -σ 0 0 σ * ω -2 2ω -1 0 2 -ω 2 -2ω 1 0 2 -ω -2ω 2    
which has a unique eigenvector (associated to 0), namely (σ, 2, -1, -2) t . This is an unstable case.

Case ε ′

s = 1 Only the case when ε ′ s = 1, q = 2ω remains to study, in which case ψ 2 vanishes. We compute φ 0 which is equal to

ψ 0 (Z) = (Z 2 -2[1 -ω]Z + 1) 2 .
The two complex conjugate roots 1 -ω ± i ω(2 -ω) are both double with modulus 1. We therefore have to study the associated stable subspaces. The only associated eigenvectors are (σ, ω ∓ i ω(2 -ω)) t respectively and the associated minimal eigensubspaces are two-dimensional, which leads to instability. If ε ′ s = 1, we should assume q < 2ω, which in physical variables reads δx > 2c ∞ /ω 1 which is not a stability condition. It should therefore be avoided to use the Lorentz-Young scheme in the harmonic case for ε s = ε ∞ when q can reach the value 2ω, i.e. if ω ≥ 1. Another way to see this condition is to give ω < 1 as a stability condition if ε ′ s = 1.

Synthesis for the Harmonic Lorentz-Young et al. Model

The scheme ( 9)-( 15) for the one-dimensional harmonic Maxwell-Lorentz equations is stable with the condition

δt < min δx √ 2c ∞ , 2 ω 1 2ε ′ s -1 if ε s > ε ∞ and δt < min δx √ 2c ∞ , √ 2 ω 1 if ε s = ε ∞ ,

Basic Polynomials in Dimension 1

The previous computations lead us to define basic polynomials associated to each one-dimensional scheme. We will see that these polynomials will prove useful in higher dimensions.

Debye (Joseph et al.)

P DJ (Z) = [1 + δε ′ s ]Z 3 -[3 + δε ′ s -(1 + δ)q]Z 2 + [3 -δε ′ s -(1 -δ)q]Z -[1 -δε ′ s ] = [1 + δε ′ s ]Y 3 + [2δε ′ s + (1 + δ)q]Y 2 + [(1 + 3δ)q]Y + [2δq].
Debye (Young)

P DY (Z) = [(1 + δα)(1 + δ)]Z 3 -[3 + δ + δα + 3δ 2 α -(1 + δ)q]Z 2 +[3 -δ -δα + 3δ 2 α -(1 -δ)q]Z -[(1 -δα)(1 -δ)] = [(1 + δ)(1 + δα)]Y 3 + [2δ(1 + α) + (1 + δ)q]Y 2 + [(1 + 3δ)q]Y + [2δq].
Lorentz (Joseph et al.)

P LJ (Z) = [1 + δ + ωε ′ s ]Z 4 -[4 + 2δ + 2ωε ′ s -(1 + δ + ω)q]Z 3 + [6 + 2ωε ′ s -2q]Z 2 -[4 -2δ + 2ωε ′ s -(1 -δ + ω)q]Z + [1 -δ + ωε ′ s ] = [1 + δ + ωε ′ s ]Y 4 + [2δ + 2ωε ′ s + (1 + δ + ω)q]Y 3 + [2ωε ′ s + (1 + 3δ + 3ω)q]Y 2 +[2(δ + 2ω)q]Y + [2ωq].
Lorentz (Kashiwa et al.)

P LK (Z) = [1 + δ + 1 2 ωε ′ s ]Z 4 -[4 + 2δ -(1 + δ + 1 2 ω)q]Z 3 + [6 -ωε ′ s + (ω -2)q]Z 2 -[4 -2δ -(1 -δ + 1 2 ω)q]Z + [1 -δ + 1 2 ωε ′ s ] = [1 + δ + 1 2 ωε ′ s ]Y 4 + [2(δ + ωε ′ s ) + (1 + δ + 1 2 ω)q]Y 3 +[2ωε ′ s + (1 + 3δ + 5 2 ω)q]Y 2 + [2(δ + 2ω)q]Y + [2ωq].
Lorentz (Young)

P LY (Z) = [1 + δ]Z 4 -[4 + 2δ -2ωε ′ s -(1 + δ)q]Z 3 + 2[3 -2ωε ′ s + (ω -1)q]Z 2 -[4 -2δ -2ωε ′ s -(1 -δ)q]Z + [1 -δ] = [1 + δ]Y 4 + [2δ + 2ωε ′ s + (1 + δ)q]Y 3 + [2ωε ′ s + (1 + 3δ + 2ω)q]Y 2 +[2(δ + 2ω)q]Y + [2ωq].
7 The Two-Dimensional Space Case

Maxwell Equations

In dimension 2, the field may be decoupled into two polarisations which lead to different schemes and also a different number of variables. We use here similar notations as those introduced in the one-dimensional case, namely λ x = c ∞ δt/δx, λ y = c ∞ δt/δy, σ x = λ x (e iξx -1), σ y = λ x (e iξy -1), q x = |σ x | 2 and q y = |σ y | 2 .

Polarisation

(B x , B y , E z )
The polarisation (B x , B y , E z ) is also called the transverse electric polarisation T E z .

B n+ 1 2 x,j,k+ 1 2 -B n-1 2 x,j,k+ 1 2 = -λ y E n z,j,k+1 -E n z,j,k , B n+ 1 2 y,j+ 1 2 ,k -B n-1 2 y,j+ 1 2 ,k = λ x E n z,j+1,k -E n z,j,k , D n+1 z,j,k -D n z,j,k = λ x B n+ 1 2 y,j+ 1 2 ,k -B n+ 1 2 y,j-1 2 ,k -λ y B n+ 1 2 x,j,k+ 1 2 -B n+ 1 2 x,j,k-1 2 . 7.1.2 Polarisation (B z , E x , E y )
The polarisation (B z , E x , E y ) is also called the transverse magnetic polarisation T M z .

B

n+ 1 2 z,j+ 1 2 ,k+ 1 2 -B n-1 2 z,j+ 1 2 ,k+ 1 2 = -λ x E n y,j+1,k+ 1 2 -E n y,j,k+ 1 2 + λ y E n x,j+ 1 2 ,k+1 -E n x,j+ 1 2 ,k , D n+1 x,j+ 1 2 ,k -D n x,j+ 1 2 ,k = λ y B n+ 1 2 z,j+ 1 2 ,k+ 1 2 -B n+ 1 2 z,j+ 1 2 ,k-1 2 , D n+1 y,j,k+ 1 2 -D n y,j,k+ 1 2 = -λ x B n+ 1 2 z,j+ 1 2 ,k+ 1 2 -B n+ 1 2 z,j-
G =      1 0 -σ y 0 0 1 σ x 0 (1+δ)σ * y 1+δε ′ s -(1+δ)σ * x 1+δε ′ s (1-δε ′ s )-(1+δ)(qx+qy) 1+δε ′ s 2δ 1+δε ′ s σ * y -σ * x -(q x + q y ) 1      associated to the variable (B n-1 2 x,j,k+ 1 2 , B n-1 2 y,j+ 1 2 ,k , E n z,j,k , D n z,j,k ) t .
The characteristic polynomial is proportional to the characteristic polynomial in dimension 1 for the same scheme with 1 as an extra root φ 0 (Z) = Y P DJ (Z).

In polynomial P DJ (Z), the variable q means q = q x + q y in the two-dimensional case. The polynomial only depends on this sum and not on the separate values of q x and q y . The general case treated in dimension one concludes to a Schur polynomial, we therefore have a von Neumann polynomial here. We also see easily on the amplification matrix in the case q = 0 (if q x only or q y only vanish, we do not have a specific case), that the eigensubspaces associated to the eigenvalue 1 are indeed stable. As for the particular case ε ′ s = 1, which gave rise to two complex conjugate eigenvalues, different from 1, we may add this new eigenvalue with the same conclusion, namely stability with the condition q = q x + q y < 4 i.e. √ 2c ∞ δt < δx if δx = δy.

G =        1 σ y 0 -σ x 0 - (1+δ)σ * y 1+δε ′ s (1-δε ′ s )-(1+δ)qy 1+δε ′ s 2δ 1+δε ′ s (1+δ)σxσ * y 1+δε ′ s 0 -σ * y -q y 1 σ x σ * y 0 (1+δ)σ * x 1+δε ′ s (1+δ)σ * x σy 1+δε ′ s 0 (1-δε ′ s )-(1+δ)qx 1+δε ′ s 2δ 1+δε ′ s σ * x σ * x σ y 0 -q x 1        associated to the variable (B n-1 2 z,j+ 1 2 ,k+ 1 2 , E n x,j+ 1 2 ,k , D n x,j+ 1 2 ,k , E n y,j,k+ 1 2 , D n y,j,k+ 1 2
) t . Anew we have a proportional polynomial to that of the one-dimensional case with two extra roots

φ 0 (Z) = Y [(1 + δε ′ s )Y + 2δε ′ s ]P DJ (Z),
which are equal to 1 and (1 -δε ′ s )/(1 + δε ′ s ) respectively and which are not roots in dimension 1 in the general case.

Only the root 1 could induce stability problems and we have seen that only the case q = 0 would make it a multiple root. The matrix G is then block-diagonal, with two rank-two blocks identical to that of the one-dimensional case. The stability is once more ensured with the condition q = q x +q y < 4.

The Debye-Young Scheme

G =      1 0 -σ y 0 0 1 σ x 0 σ * y 1+δα -σ * x 1+δα (1+δ)(1-δα)+4δ 2 α-(1+δ)q (1+δ)(1+δα) 1-δ 1+δ 2δ 1+δα 0 0 2δα 1+δ 1-δ 1+δ      associated to the variable (B n-1 2 x,j,k+ 1 2 , B n-1 2 y,j+ 1
2 ,k , E n z,j,k , P n z,j,k ) t . The computation of the characteristic polynomial leads to the same polynomial as in one dimension but with 1 as an extra eigenvalue φ 0 (Z) = Y P DY (Z).

Anew we revert to the arguments of the one-dimensional case, the double eigenvalue 1 of the q = 0 case not being a problem.

Polarisation

(B z , E x , E y )
Coupling polarisation (B z , E x , E y ) with the Debye-Young scheme, we obtain the amplification matrix

G =         1 σ y 0 -σ x 0 - σ * y 1+δα (1+δ)(1-δα)+4δ 2 α-(1+δ)qy (1+δ)(1+δα) 1-δ 1+δ 2δ 1+δα σxσ * y 1+δα 0 0 2δα 1+δ 1-δ 1+δ 0 0 σ * x 1+δα σ * x σy 1+δα 0 (1+δ)(1-δα)+4δ 2 α-(1+δ)qx (1+δ)(1+δα) 1-δ 1+δ 2δ 1+δα 0 0 0 2δα 1+δ 1-δ 1+δ         associated to the variable (B n-1 2 z,j+ 1 2 ,k+ 1 2 , E n x,j+ 1 2 ,k , P n x,j+ 1 2 ,k , E n y,j,k+ 1 2 , P n y,j,k+ 1 2
) t . The computation of the characteristic polynomial leads to the same polynomial as in dimension 1 but with two extra eigenvalues

φ 0 (Z) = Y [(1 + δ)(1 + δα)Y + 2δ(1 + α)]P DY (Z)
which are 1 and (1-δ)(1-δα)/(1+δ)(1+δα). Anew we revert to the argument in the one-dimensional case, the triple eigenvalue 1 of case q = 0 not being a problem.

G =        1 0 -σ y 0 0 0 1 σ x 0 0 2δσ * y 1+δ+ωε ′ s -2δσ * x 1+δ+ωε ′ s 2-(1+δ+ω)q (1+δ+ωε ′ s ) -1-δ+ωε ′ s 1+δ+ωε ′ s 2ω 1+δ+ωε ′ s 0 0 1 0 0 σ * y -σ * x -q 0 1        associated to the variable (B n-1 2 x,j,k+ 1 2 , B n-1 2 y,j+ 1 2 ,k , E n z,j,k , E n-1 z,j,k , D n z,j,k ) t .
The computation of the characteristic polynomial leads to the same polynomial as in dimension 1 but with the extra eigenvalue 1

φ 0 (Z) = Y P LJ (Z).
The triple eigenvalue in the q = 0 case is not a problem.

G =             1 σ y 0 0 -σ x 0 0 - 2δσ * y 1+δ+ωε ′ s 2-(1+δ+ω)qy 1+δ+ωε ′ s -1-δ+ωε ′ s 1+δ+ωε ′ s 2ω 1+δ+ωε ′ s σxσ * y (1+δ+ω) 1+δ+ωε ′ s 0 0 0 1 0 0 0 0 0 -σ * y -q y 0 1 σ x σ * y 0 0 2δσ * x 1+δ+ωε ′ s σ * x σy(1+δ+ω) 1+δ+ωε ′ s 0 0 2-(1+δ+ω)qx 1+δ+ωε ′ s -1-δ+ωε ′ s 1+δ+ωε ′ s 2ω 1+δ+ωε ′ s 0 0 0 0 1 0 0 σ * x σ * x σ y 0 0 -q x 0 1             associated to the variable (B n-1 2 z,j+ 1 2 ,k+ 1 2 , E n x,j+ 1 2 ,k , E n-1 x,j+ 1 2 ,k , D n x,j+ 1 2 ,k , E n y,j,k+ 1 2 , E n-1 y,j,k+ 1 2 , D n y,j,k+ 1 2
) t . The characteristic polynomial is once more proportional to the one-dimensional polynomial

φ 0 (Z) = Y [(1 + δ + ωε ′ s )Y 2 + 2(δ + ωε ′ s )Y + 2ωε ′ s ]P LJ (Z) = Y [(1 + δ + ωε ′ s )Z 2 -2Z + (1 -δ + ωε ′ s )]P LJ (Z).
In the anharmonic case, and by the von Neumann technique, we check easily that

ψ 0 (Z) = [1 + δ + ωε ′ s ]Z 2 -2Z + [1 -δ + ωε ′ s ]
is a Schur polynomial. Besides the double root 1 if q = 0 is still no problem.

In the harmonic case, ψ 0 (Z) has two distinct complex conjugate roots with modulus 1. This is not a problem in itself, except if ε ′ s = 1 in which case ψ 0 (Z) is also a factor in P LJ (Z)

φ 0 (Z) = (Z -1)[(1 + ωε ′ s )Z 2 -2Z + (1 + ωε ′ s )]P LJ (Z) = (Z -1)[(1 + ωε ′ s )Z 2 -2Z + (1 + ωε ′ s )] 2 [Z 2 -(2 -q)Z + 1].
This is already the case when we have detected double eigenvalues in dimension 1, giving rise to instabilities for q = 2ω/(1 + ω). It is better to avoid this scheme in the case when ε ′ s .

7.5 The Lorentz-Kashiwa et al. Scheme

G =         1 0 -σ y 0 0 0 1 σ x 0 0 σ * y (D-1 2 ωα) D - σ * x (D-1 2 ωα) D (1-q)D-(2-q) 1 2 ωα D ω D -1 D σ * y 1 2 ωα D - σ * x 1 2 ωα D (2-q) 1 2 ωα D D-ω D 1 D σ * y ωα D -σ * x ωα D (2-q)ωα D -2ω D 2-D D         associated to the variable (B n-1 2 x,j,k+ 1 2 , B n-1 2 y,j+ 1 2 ,k , E n z,j,k , P n z,j,k , J n z,j,k ) t .
The calculation of the characteristic polynomial leads to the same polynomial as in dimension 1 with the extra root 1

φ 0 (Z) = Y P LK (Z).
The triple eigenvalue 1 of case q = 0 is not a problem. In the anharmonic case, and by the von Neumann technique, we check easily that

G =              1 -σ y 0 0 σ x 0 0 σ * y (D-1 2 ωα) D (1-qy)D-(2-qy) 1 2 ωα D ω D -1 D σxσ * y (D-1 2 ωα) D 0 0 σ * y 1 2 ωα D (1-qy) 1 2 ωα D D-ω D 1 D σxσ *
ψ 0 (Z) = [1 + δ + 1 2 ωε ′ s ]Z 2 -[2 -ωε ′ s ]Z + [1 -δ + 1 2 ωε ′ s ]
is a Schur polynomial. Besides, the root 1, which is a double one if q = 0, does not lead to any problem.

In the harmonic case, we have the extra roots 1 and two complex conjugate roots of modulus 1, which are not roots of P LK (Z). The stability is hence given under the same conditions as in the one-dimensional case.

7.6 The Lorentz-Young Scheme The triple eigenvalue 1 of the case q = 0 is not a problem. In the anharmonic case, and by von Neumann technique, we check easily that

ψ 0 (Z) = [1 + δ]Z 2 -[2 -2ωε ′ s ]Z + [1 -δ]
is a Schur polynomial. Besides, the root 1 which is a double one if q = 0 is not a problem.

In the harmonic case, we have the extra roots 1 and two complex conjugate roots of modulus 1, which are not roots of P LY (Z). The stability is therefore ensured under the same conditions as in the one-dimensional case.

Conclusion

We have studied the stability of numerical schemes for Maxwell-Debye and Maxwell-Lorentz equations in space dimension 1 and 2. In dimension 2, the characteristic polynomials of each scheme and in both polarisation happen to be proportional to the characteristic polynomials for the same scheme in space dimension 1. In all the cases, the extension to dimension 2 goes with an extra root 1 compared to the one-dimensional case. This is the only extra root in the T E z polarisation. For the T M z polarisation, there is one other extra root for the Debye equation and two other extra roots for the Lorentz equation, all these roots being on the unit circle. For the Yee scheme applied to the raw Maxwell equations, the stability condition is q ≤ 4 in dimensions 1, 2 et 3, recalling that q = q x + q y in dimension 2 (q = max(q x + q y , q x + q z , q y + q z ) in dimension 3). The results are gathered in two tables according to ε s = ε ∞ or not. For each model, we have at least one scheme for which the stability condition is the same as for the raw Maxwell equations (q < 4). In Young models, the extra conditions correspond to a fine enough discretization of Debye and Lorentz equations respectively,... because stability is not the only issue. Applications to classical materials show in general that the condition due to the Maxwell equations is the more restrictive one and not conditions due to the constitutive law of the material.

Model

Computations in dimension 3 are too tedious to be carried out by hand. They have been automated (see [START_REF] Bidégaray-Fesquet | NAUtil: von Neumann Analysis Utilities[END_REF]).

7. 3 . 1

 31 Polarisation (B x , B y , E z )Coupling polarisation (B x , B y , E z ) with Debye-Young scheme, we obtain the amplification matrix

7. 4 . 2

 42 Polarisation (B z , E x , E y ) Coupling polarisation (B z , E x , E y ) with the Lorentz-Joseph et al. scheme, we obtain the amplification matrix

7. 5 . 1

 51 Polarisation (B x , B y , E z ) Coupling polarisation (B x , B y , E z ) with the Lorentz-Kashiwa et al. scheme, we obtain the amplification matrix

7. 5 . 2

 52 Polarisation (B z , E x , E y ) Coupling polarisation (B z , E x , E y ) with the Lorentz-Kashiwa et al. scheme, we obtain the amplification matrix

2 )

 2 t . The computation of the characteristic polynomial leads to a polynomial proportional to the one-dimensional oneφ 0 (Z) = Y [(1 + δ + 1 2 ωε ′ s )Y 2 + (2(δ + ωε ′ s ))Y + (2ωε ′ s )]P LK (Z) = Y [(1 + δ + 1 2 ωε ′ s )Z 2 -(2 -ωε ′ s )Z + (1 -δ + 1 2 ωε ′ s )]P LK (Z).

7. 6 . 1 1 2 1 2z

 6111 Polarisation (B x , B y , E z )Coupling polarisation (B x , B y , E z ) with the Lorentz-Young scheme, we obtain the amplification matrix ,k , E n z,j,k , P n z,j,k , J n-,j,k ) t . The computation of the characteristic polynomial leads to the same polynomial as in dimension 1 but with 1 as an extra eigenvalue φ 0 (Z) = Y P LY (Z).

7. 6 . 2 2 )

 622 Polarisation (B z , E x , E y )Coupling polarisation (B z , E x , E y ) with the Lorentz-Young scheme, we obtain the amplification matrix t . The computation of the characteristic polynomial leads to a polynomial proportional to that of dimension 1φ 0 (Z) = Y [(1 + δ)Y 2 + (2(δ + ωε ′ s ))Y + (2ωε ′ s )]P LY (Z) = Y [(1 + δ)Z 2 -(2 -ωε ′ s )Z + (1 -δ)]P LY (Z).

  The Debye-Joseph et al. Scheme7.2.1 Polarisation (B x , B y , E z ) Coupling polarisation (B x , B y , E z ) with the Debye-Joseph et al. scheme, we obtain the amplification matrix

	1 2 ,k+ 1 2	.
	7.2

  7.2.2 Polarisation (B z , E x , E y ) Coupling polarisation (B z , E x , E y ) with the Debye-Joseph et al. scheme yields the amplification matrix

  7.4 The Lorentz-Joseph et al. Scheme 7.4.1 Polarisation (B x , B y , E z ) Coupling polarisation (B x , B y , E z ) with the Lorentz-Joseph et al. scheme, we obtain the amplification matrix

Table 1 :

 1 Stability of schemes for ε s > ε ∞ .

		Scheme		dimension 1	dimension 2 (δx = δy)
	Debye	Joseph et al.	q ≤ 4	δt ≤ δx c∞		δt ≤ δx √ 2c∞
	Debye	Young	q ≤ 4, δ ≤ 1	δt ≤ min( δx c∞ , 2t r )	δt ≤ min( δx √ 2c∞ , 2t r )
	Lorentz Joseph et al.	q ≤ 2	δt ≤ δx √ 2c∞		δt ≤ δx 2c∞
	Lorentz Kashiwa et al.	q < 4	δt < δx c∞		δt < δx √ 2c∞
	Lorentz	Young	q ≤ 2, ω ≤ 2 2ε ′ s -1	δt ≤ min( δx √ 2c∞ ,	ω 1	2 √ 2ε ′ s -1	) δt ≤ min( δx 2c∞ ,	ω 1	2 √ 2ε ′ s -1	)
	Harm.	Joseph et al.	q ≤ 2	δt ≤ δx √ 2c∞		δt ≤ δx 2c∞
	Harm. Kashiwa et al.	q < 4	δt < δx c∞		δt < δx √ 2c∞
	Harm.	Young	q < 2, ω ≤ 2 2ε ′ s -1 or q ≤ 2, ω < 2 2ε ′ s -1	δt < min( δx √ 2c∞ ,	ω 1	2 √ 2ε ′ s -1	) δt < min( δx 2c∞ ,	ω 1	2 √ 2ε ′ s -1	)

Model

Scheme dimension 1 dimension 2 (δx = δy) Harm.

Table 2: Stability of schemes for ε s = ε ∞ .