Bernard Collet

Michel Destrade

Love waves in a piezoelectric layered structure

Consider a layer made of a m3m insulator crystal, with faces cut parallel to a symmetry plane. Then bond it onto a semi-infinite mm2 piezoelectric substrate. For a X-or Y -cut of the substrate, a Love wave can propagate in the resulting structure and the corresponding dispersion equation is derived analytically. It turns out that a fully explicit treatment can also be conducted in the case of a Y -cut rotated about Z. In the case of a germanium layer over a potassium niobate substrate, the wave exists at any wavelength for X-and Y -cuts but this ceases to be the case for rotated cuts, with the appearance of forbidden ranges. By playing on the cut angle, the Love wave can be made to travel faster than, or slower than, or at the same speed as, the shear bulk wave of the layer. A by-product of the analysis is the derivation of the explicit secular equation for the Bleustein-Gulyaev wave in the substrate alone, which corresponds to an asymptotic behavior of the Love wave.

Introduction

Layered structures, especially film/coating substrate systems, play an important role in micro-electro-mechanical systems (MEMS) and in microelectronics packages. In order to achieve high performance, many surface acoustic wave (SAW) devices/sensors are made of layered architectures such as, for instance, a dielectric, or a piezoelectric, or a non-piezoelectric semiconductor (finite-thickness) layer deposited onto a (semi-infinite) piezoelectric substrate. For certain configurations, it is possible to have a one-component wave travel in the structure, in the direction of the interface: this guided (shear-horizontal) Love wave leaves the upper face of the layer free of mechanical tractions, its amplitude varies sinusoidally through the thickness of the layer and then decays rapidly with depth in the substrate, and it is such that all fields are continuous at the layer/substrate interface.

Love waves in piezoelectric layered acoustic devices are most suitable for high frequency filters because of their high phase velocity, and they also show great promise in bio-sensors applications with liquid environments because of their high sensitivity. Consequently, they have received much attention over the years. For example, Lardat et al. [START_REF] Lardat | Theory and performance of acoustical dispersive surface wave delay lines[END_REF] found under which conditions a piezoelectrically stiffened Love wave exist and presented experimental and analytical results on surface wave delay lines. Kessenikh et al. [START_REF] Kessenikh | Love surface waves in piezoelectrics[END_REF] investigated surface Love waves in piezoelectric substrates of classes 6, 4, 6mm, 4mm, 622, and 422 with an isotropic dielectric layer. Hanhua and Xingjiao [START_REF] Hanhua | Shear-horizontal surface waves in a layered structure of piezoelectric ceramics[END_REF] studied Love waves for a structure made of a 6mm piezoelectric layer and a 6mm piezoelectric substrate, with a common symmetry axis in the plane of the interface, while Darinskii and Weihnacht [START_REF] Darinskii | Supersonic Love waves in strong piezoelectrics of symmetry mm2[END_REF] had a similar structure, made of 2mm piezoelectric layer and substrate with common symmetry axes, one of which is aligned with the propagation direction and another is aligned with the normal to the interface. Jakoby and Vellekoop [START_REF] Jakoby | Properties of Love waves: applications in sensors[END_REF] reviewed the properties of Love waves and associated numerical methods for a piezoelectric/piezoelectric layered structure, when the substrate is made of ST-cut quartz; so did Ogilivy [START_REF] Ogilvy | The mass-loading sensitivity of acoustic Love waves biosensors in air[END_REF], with special emphasis on the mass-sensitivity loading of biosensors. The reader can find additional pointers to the literature on Love waves for piezoelectric sensors in those articles and in the references therein, as well as in the reviews by Farnell and Adler [START_REF] Farnell | Elastic wave propagation in thin layers[END_REF] or by Gulyaev [START_REF] Yu | Review of shear surface acoustic waves in solids[END_REF].

The present work is concerned with the propagation of Love waves in a composite structure, made of a m3m cubic insulator (non-piezoelectric semiconductor) layer of finite thickness h, bonded onto a Y -cut rotated about Z, mm2 orthorhombic, piezoelectric substrate, see Fig. 1. With this singular structure, it is possible to exploit the potentiality offered by an efficient direct analytic method. The steps to be followed are: (i) Determine the general solution for the layer, satisfying the boundary conditions on the upper face; (ii) Derive some fundamental equations in the substrate using the Stroh formulation of the equations of motion there; (iii) Match the solutions at the layer/substrate interface; and (iv) Deduce the dispersion equation for piezoelectric Love waves in explicit form. The dispersion relation admits a denumerable number of solutions (branches) and the solutions for the mechanical displacements, shear stresses, electric potentials, and electric inductions are deduced explicitly in the case of a metalized, mechanically free, upper surface, brought to the zero electric potential (short-circuit). In passing, the explicit secular equations for the Bleustein-Gulyaev wave speed and for the limiting wave speed (substrate only, no overlayer) are found as a cubic and as a sextic in the squared wave speed respectively, for rotated cuts. The special cases of a X-cut or Y -cut are treated separately. The effects of the angle of the cut on the phase velocity of the first modes are illustrated numerically for a specific layered structure, namely a germanium layer over a potassium niobate substrate and the appearance of a forbidden band of frequency is uncovered for a rotated cut, in sharp contrast with the non-rotated cuts where the waves exist for all frequencies.

The layer

First consider the upper layer, which is made of a cubic, m3m non-piezoelectric (semi-conductor) crystal, with mass density ρ. For two-dimensional motions (independent of x 3 ), the anti-plane equation of motion decouples from its in-plane counterpart and reads (0

≤ x 2 ≤ h) ĉ44 (û 3,11 + û3,22 ) = ρû 3,tt , (1)
where û3 is the anti-plane mechanical displacement, ĉ44 the transverse stiffness, and the comma denotes partial differentiation. For a solution in the form of an inhomogeneous wave traveling with speed v and wave number k in the x 1 direction, such as

û3 = Û3 (kx 2 )e ik(x 1 -vt) , (2)
(where Û3 is a function of kx 2 alone), the equation of motion (1) reduces to

Û′′ 3 + v 2 v2 -1 Û3 = 0, where v := ĉ44 /ρ. (3)
The general solution to this second-order differential equation is either (i)

Û3 (kx 2 ) = Û3 (0) cos v 2 v2 -1kx 2 +A sin v 2 v2 -1kx 2 , when v > v, (4)
or (ii)

Û3 (kx 2 ) = Û3 (0) cosh 1 - v 2 v2 kx 2 +A sinh 1 - v 2 v2 kx 2 , when v < v, ( 5
)
where A is a constant to be determined from the boundary condition at x 2 = h. This latter condition is that the upper face of the layer be free of mechanical tractions, so that σ23 = ĉ44 û3,2 is zero there. Then it follows from (4) and ( 5) that: A = tan v 2 /v 2 -1kh in Case (i), and that: A = -tanh 1 -v 2 /v 2 kh in Case (ii). Consequently the mechanical field is now entirely known in the layer. The procedure to find the electrical field is similar, and even simpler because the layer is not piezoelectric. Taking the electric potential φ in the form φ = φ(kx 2 )e ik(x 1 -vt) ,

(where φ is a function of kx 2 alone), Poisson's equation: ∆ φ = 0 reduces to φ′′ -φ = 0, [START_REF] Farnell | Elastic wave propagation in thin layers[END_REF] with general solution

φ(kx 2 ) = φ(0) (cosh kx 2 + B sinh kx 2 ) , (8)
where B is a constant to be determined from the boundary condition on the upper face of the layer x 2 = h. A thin metallic film covers that face, and it is grounded to potential zero: φ(kh) = 0. It then follows from ( 8) that B = -coth kh, and consequently, that the electrical field is now entirely known in the layer.

For the problem at hand, only the values of the fields at the interface x 2 = 0 between the layer and the substrate are needed. Because the mechanical displacement û3 and the electrical potential φ are in the forms (2) and [START_REF] Ogilvy | The mass-loading sensitivity of acoustic Love waves biosensors in air[END_REF], respectively, the mechanical traction σ23 and the electrical displacement D2 must also be of a similar form, due to the constitutive relations: σ23 = ĉ44 û3,2 and D2 = -ǫ 11 φ ,2 , where ǫ11 is the dielectric constant of the layer. Hence,

σ23 = ik t23 (kx 2 )e ik(x 1 -vt) , D2 = ik d2 (kx 2 )e ik(x 1 -vt) , (9)
(say), where t23 and d2 are functions of kx 2 alone. In particular, the conclusions drawn from the calculations conducted above are that at the layer/substrate interface,

t23 (0) = -iĉ Û3 (0), d2 (0) = -iǫ φ(0), (10)
where

ĉ := ĉ44 v 2 v2 -1 tan v 2 v2 -1kh, when v > v, ĉ := -ĉ 44 1 - v 2 v2 tanh 1 - v 2 v2 kh, when v < v, ǫ := ǫ11 coth kh. ( 11
)
3 The substrate

The substrate occupies the half-space x 2 ≤ 0 and is made of the superstrong piezoelectric crystal, potassium niobate KbNO 3 , with orthorhombic mm2 symmetry. Cut the crystal along a plane containing the Z axis and making an angle θ with the XY plane. Let x 1 x 2 x 3 be the coordinate system obtained after the rotation

  m n 0 -n m 0 0 0 1   , where m = cos θ, n = sin θ. (12)
Yet again, for a two-dimensional motion (independent of x 3 ) the anti-plane strain and stress decouple from their in-plane counterparts. Hence, with u 1 = u 2 = 0, u 3 = u 3 (x 1 , x 2 , t) for the mechanical displacement, and φ = φ(x 1 , x 2 , t) for the electric potential, the constitutive relations yield σ 11 = σ 22 = σ 33 = σ 12 = 0 for the stress components and D 3 = 0 for the electrical displacement, and they reduce to

σ 23 = c 44 u 3,2 + c 45 u 3,1 + e 14 φ ,1 + e 24 φ ,2 , σ 31 = c 45 u 3,2 + c 55 u 3,1 + e 15 φ ,1 + e 14 φ ,2 , D 1 = e 14 u 3,2 + e 15 u 3,1 -ǫ 11 φ ,1 -ǫ 12 φ ,2 , D 2 = e 24 u 3,2 + e 14 u 3,1 -ǫ 12 φ ,1 -ǫ 22 φ ,2 , (13)
where the c ij , e ij , ǫ ij are related to the corresponding quantities cij , ẽij , ǫij in the crystallographic coordinate system XY Z through tensor transformations [START_REF] Auld | Acoustic fields and waves in solids[END_REF] as:

c 44 = m 2 c44 + n 2 c55 , c 55 = n 2 c44 + m 2 c55 , c 45 = mn(c 44 -c55 ), e 24 = m 2 ẽ24 + n 2 ẽ15 , e 15 = n 2 ẽ24 + m 2 ẽ15 , e 14 = mn(ẽ 24 -ẽ15 ), ǫ 11 = m 2 ǫ11 + n 2 ǫ22 , ǫ 22 = n 2 ǫ11 + m 2 ǫ22 , ǫ 12 = mn(ǫ 22 -ǫ11 ). ( 14
)
Now consider the x 2 -cut, x 1 -propagation of a Shear-Horizontal interface acoustic wave that is, a motion with speed v and wave number k where the displacement field u 3 and the electric potential φ are of the form,

{u 3 , φ}(x 1 , x 2 , t) = {U 3 (kx 2 ), ϕ(kx 2 )}e ik(x 1 -vt) , (15)
(say), with

U 3 (-∞) = 0, ϕ(-∞) = 0. ( 16
)
It follows from the constitutive equations ( 13) that the traction σ 32 and the electric induction D 2 are of a similar form,

{σ 32 , D 2 }(x 1 , x 2 , t) = ik{t 32 (kx 2 ), d 2 (kx 2 )}e ik(x 1 -vt) , (17)
(say), where

t 32 (-∞) = 0, d 2 (-∞) = 0. ( 18
)
The non-trivial part of the equations of piezoacoustics, σ ij,j = ρu i,tt , D i,i = 0 (where ρ is the mass density of the crystal), can be written as a second-order differential system [START_REF] Mozhaev | Sectors of nonexistence of surface acoustic waves in potassium niobate[END_REF]:

T U ′′ 3 ϕ ′′ + 2iR U ′ 3 ϕ ′ -[Q -(ρv 2 )J] U 3 ϕ = 0 0 , (19)
where the prime denotes differentiation with respect to kx

or as a first-order differential system in the form

ξ ′ = iNξ, where ξ(kx 2 ) := [U 3 , ϕ, t 31 , d 2 ] t , (21)
and N has the Stroh block structure

N = N 1 N 2 N 3 + (ρv 2 )J N t 1 , (22)
with [START_REF] Ting | Anisotropic elasticity: theory and applications[END_REF]]

N 1 = -T -1 R, N 2 = T -1 , N 3 = RT -1 R -Q. ( 23
)
Here the components of N are easily computed from ( 20) and ( 23), by hand or using a Computer Algebra System. Seeking a solution to [START_REF] Zgonik | Material constants of KNbO 3 relevant for electro-and acousto-optics[END_REF] in the form:

[U 3 , ϕ] t = [U 3 (0), ϕ(0)] t e ikqx 2
, where q is constant, yields the propagation condition

det [q 2 T + 2qR + Q -(ρv 2 )J] = 0, (24)
a quartic in q. Explicitly,

c D 44 q 4 + 2 ǫ 12 ǫ 22 c D 16 q 3 + (c D 45 -ρv 2 )q 2 + 2 ǫ 12 ǫ 22 (c D 26 -ρv 2 )q + ǫ 11 ǫ 22 (c D 55 -ρv 2 ) = 0, (25) where c
Out of the four possible roots, only two have a negative imaginary part, insuring exponential decay as x 2 → -∞. Computing these qualifying roots analytically is not an easy matter, because here the speed v is still an unknown. Although it is actually possible to do so [START_REF] Fu | An explicit expression for the surface-impedance matrix of a generally anisotropic incompressible elastic material in a state of plane strain[END_REF][START_REF] Destrade | The speed of interface waves polarized in a symmetry plane[END_REF], this paper follows a different route and makes extensive use of the "fundamental equations" derived in [START_REF] Destrade | Elastic interface acoustic waves in twinned crystals[END_REF][START_REF] Destrade | Explicit secular equation for Scholte waves over a monoclinic crystal[END_REF][START_REF] Destrade | On interface waves in misoriented pre-stressed incompressible elastic solids[END_REF]. They read

ξ(0) • M (n) ξ(0) = 0, where M (n) := 0 I I 0 N n , ( 27
)
and n is any positive or negative integer. [START_REF] Darinskii | Supersonic Love waves in strong piezoelectrics of symmetry mm2[END_REF] The layer/substrate structure

The continuity of all fields at the layer/substrate interface x 2 = 0 imposes the boundary conditions:

U 3 (0) = Û3 (0), ϕ(0) = φ(0), t 23 (0) = t23 (0), d 2 (0) = d2 (0). ( 28
)
Now the dispersion equations are derived explicitly. In the special cases of X and Y -cuts, the dispersion equation is exact: if it is satisfied, then the Love wave exist. In the other cases, the dispersion equation is rationalized : it has spurious roots, not corresponding to a true Love wave so that a subsequent validity check is necessary.

4.1 Special cases θ = 0, 90

•
For a X-cut or a Y -cut of the substrate, the analysis simplifies considerably and a direct treatment is possible, leading to an exact dispersion equation, and not requiring the use of the fundamental equations. When m = 1, n = 0, or m = 0, n = 1, the parameters c 45 , e 14 , ǫ 12 vanish according to [START_REF] Destrade | Elastic interface acoustic waves in twinned crystals[END_REF]. Then in [START_REF] Collet | Explicit secular equations for piezoacoustic surface waves: Shear-Horizontal modes[END_REF], R ≡ 0 also, and the quartic (24) becomes the following biquadratic [START_REF] Darinskii | Supersonic Love waves in strong piezoelectrics of symmetry mm2[END_REF]:

q 4 -Sq 2 + P = 0, ( 29
)
where the non-dimensional quantities S and P are given by (30)

-
The relevant roots q 1 and q 2 come in one of the two following forms, either

(a) q 1 = -iβ 1 , q 2 = -iβ 2 , or (b) q 1 = -α -iβ, q 2 = α -iβ, ( 31
)
where β i > 0, β > 0. In either case, q 1 + q 2 has no real part and a negative imaginary part, and q 1 q 2 is a negative real number. Explicitly,

q 1 + q 2 = -i 2 √ P -S, q 1 q 2 = - √ P . (32)
The associated eigenvectors ζ 1 and ζ 2 follow from (for instance) the third column of the adjoints to the matrices N -q 1 I and N -q 2 I, as

ζ i = [a i , b i ] t , (i = 1, 2)
where

a i = q 2 i + ǫ 11 ǫ 22 , e 24 ǫ 22 q 2 i + e 15 ǫ 22 t , b i = q i T a i , (33)
(here T is given by [START_REF] Collet | Explicit secular equations for piezoacoustic surface waves: Shear-Horizontal modes[END_REF], with components evaluated at 0 • or 90 • ). Then the general solution to the equations of motion is of the form

ξ(kx 2 ) = γ 1 e ikq 1 x 2 ζ 1 + γ 2 e ikq 2 x 2 ζ 2 , ( 34
)
where γ 1 , γ 2 are disposable constants. At x 2 = 0, it can be split as

[U 3 (0), ϕ(0)] t = Aγ, [t 23 (0), d 2 (0)] t = Bγ, (35)
where

A := [a 1 |a 2 ], B := [b 1 |b 2 ], γ := [γ 1 , γ 2 ] t . (36)
Now the boundary conditions ( 28) and [START_REF] Mozhaev | Sectors of nonexistence of surface acoustic waves in potassium niobate[END_REF] give the link

Bγ = [ t23 (0), d2 (0)] t = Diag(-iĉ, -iǫ)[ Û3 (0), φ(0)] t = -i Diag(ĉ, ǫ)Aγ, (37)
from which the dispersion equation follows as

|iBA -1 -Diag(ĉ, ǫ)| = 0. ( 38
)
It is written in this form to take advantage of the many desirable properties of the surface impedance tensor iBA -1 [START_REF] Abbudi | On the existence of interfacial (Stoneley) waves in bonded piezoelectric half-spaces[END_REF]; thus, this matrix is Hermitian in the subsonic range (defined below), and of the compact form:

iBA -1 = i q 1 + q 2 ρv 2 -c 55 -c 44 q 1 q 2 -e 15 + e 24 q 1 q 2 -e 15 + e 24 q 1 q 2 ǫ 11 -

ǫ 22 q 1 q 2 . ( 39
)
Moreover, the eigenvalues of the aggregate impedance tensor in (38) are real monotonic decreasing functions of v for any fixed kh, so that the wave speed of each mode is obtained unambiguously from the roots of (38), see Shuvalov and Every [START_REF] Shuvalov | Some properties of surface acoustic waves in anisotropic-coated solids, studied by the impedance method[END_REF].

Using the identities q 2 1 + q 2 2 = S, q 2 1 q 2 2 = P and the connections (30) and (32), the exact, explicit, dispersion equation is finally derived as ρv 2

√ P -S + ǫ = 0, (40)
which is a fully explicit equation, because ĉ, ǫ are defined in [START_REF] Mozhaev | Sectors of nonexistence of surface acoustic waves in potassium niobate[END_REF] and S, P in (30). The dispersion equation is valid in the subsonic range that is, as long as the speed is less than the limiting speed v L , the smallest speed at which the biquadratic (29) ceases to have two roots q 1 , q 2 with a positive imaginary part. When q 1 , q 2 are of the form (a) in (31), then v L is root to: P = 0; when they are of the form (b), then v L is root to: 2 √ P = S. In either case, the wave becomes homogeneous at v = v L because then the roots of the biquadratic are real ; they are: ±q, where q = √ S L or q = S L /2, according to each case (here S L is S given by (30

) at v = v L ).
The vanishing of the wave away from the interface can no longer be insured then, but the continuity of the fields at x 2 = 0 can. The conclusion is that the boundary conditions (28) and [START_REF] Mozhaev | Sectors of nonexistence of surface acoustic waves in potassium niobate[END_REF] This equation has an infinity of roots in kh (corresponding to the intersections of the graph of tan with the graph of coth). Each root (kh) L say, is a cut-off parameter for each dispersion mode, at which the Love wave ceases to exist. Note that the dispersion equation ( 40) is consistent with the secular equation of a Bleustein-Gulyaev wave traveling in the substrate alone: as h tends to zero, ĉ → 0, ǫ → ∞ so that it reduces to: ρv 2 -c 55 = c 44 √ P , which, once squared, coincides with the quadratic (56), obtained in the next Section. It is also consistent with the dispersion equation for a purely elastic Love wave. Indeed, by taking e ij → 0 and ǫ ij → 0 in (30) and ( 40), the equation of Lardat et al. [START_REF] Lardat | Theory and performance of acoustical dispersive surface wave delay lines[END_REF] is recovered:

tan ρv 2 -ĉ44 ĉ44 kh = c 44 c 55 -ρv 2 c 44 ĉ44 ρv 2 -ĉ44 ĉ44 . ( 43
)
Finally, it is consistent with the dispersion equation of Love surface waves in an isotropic dielectric layer over a 6mm piezoelectric substrate [START_REF] Kessenikh | Love surface waves in piezoelectrics[END_REF], by the corresponding specialization.

Rotated cut

Combining the boundary conditions (28) with the results for t23 (0) and d2 (0) of ( 10) gives the following form for ξ(0):

ξ(0) = U 3 (0)[1, α, -iĉ, -iǫα] t , (44)
where α := ϕ(0)/U 3 (0) is complex: α = α 1 + iα 2 , say. Then the fundamental equations ( 27) read

    1 α iĉ iǫα          M (n) 11 M (n) 12 M (n) 13 M (n) 14 M (n) 12 M (n) 22 M (n) 23 M (n) 24 M (n) 13 M (n) 23 M (n) 33 M (n) 34 M (n) 14 M (n) 24 M (n) 34 M (n) 44          1 α -iĉ -iǫα     = 0, (45) or [
M (n) 12 + ĉǫM (n) 34 ](2α 1 ) + [ǫM (n) 14 -ĉM (n) 23 ](2α 2 ) + [M (n) 22 + ǫ2 M (n) 44 ](α 2 1 + α 2 2 ) = -[M (n) 11 + ĉ2 M (n) 33 ]. ( 46
)
Writing them for n = -1, 1, 2, and re-arranging the three resulting equations, leads to the following non-homogeneous system of linear equations,

[k 1 |k 2 |k 3 ] p = -k 4 , (47)
where p := [2α 1 , 2α 2 , α 2 1 + α 2 2 ] t and k 1 , k 2 , k 3 , k 4 are the vectors with components:

M (n) 12 + ĉǫM (n) 34 , ǫM (n) 14 -ĉM (n) 23 , M (n) 22 + ǫ2 M (n) 44 , M (n) 11 + ĉ2 M (n)
33 , (48) (n = -1, 1, 2) respectively. Cramer's rule gives the unique solution to the system as

2α 1 = -∆ 1 /∆, 2α 2 = -∆ 2 /∆, α 2 1 + α 2 2 = -∆ 3 /∆, (49)
where

∆ := det[k 1 |k 2 |k 3 ], ∆ 1 := det[k 4 |k 2 |k 3 ], ∆ 2 := det[k 1 |k 4 |k 3 ]
, and

∆ 3 := det[k 1 |k 2 |k 4 ]
. The dispersion equation follows then from the compatibility of the equalities (49):

∆ 2 1 + ∆ 2 2 + 4∆ 3 ∆ = 0. (50)
When (and if) this dispersion relation yields a real positive wave speed v for a given wave number k, it remains to be checked whether that speed corresponds to a valid solution. Proceed as follows. First recall that the exact boundary condition is of the form (38), where now the a i , b i (i = 1, 2) are defined by

a i := q 2 i + 2 ǫ 12 ǫ 22 q i + ǫ 11 ǫ 22 , e 24 ǫ 22 q 2 i + 2 e 14 ǫ 22 q i + e 15 ǫ 22 t , b i := (q i T + R)a i . (51)
The computation of the corresponding surface impedance tensor iBA -1 is long but perfectly possible analytically; its components depend on q 1 , q 2 through the sum q 1 + q 2 and the product q 1 q 2 . Now, having found a speed from (50), compute numerically the roots of the quartic [START_REF] Abbudi | On the existence of interfacial (Stoneley) waves in bonded piezoelectric half-spaces[END_REF]. Select q 1 and q 2 , the roots with negative imaginary parts (if there are no such roots, then v is not valid.) Then compute q 1 + q 2 , q 1 q 2 , and iBA -1 , and check whether the exact boundary condition (38) is satisfied. Finally, it is also possible to determine exactly the limiting speed v L at which the imaginary part of an attenuation factor first vanishes. Fu [START_REF] Fu | An explicit expression for the surface-impedance matrix of a generally anisotropic incompressible elastic material in a state of plane strain[END_REF] shows that v L is the smallest root of D = 0, where D is the discriminant of the cubic resolvent associated with the quartic [START_REF] Abbudi | On the existence of interfacial (Stoneley) waves in bonded piezoelectric half-spaces[END_REF]. Thus rewrite the quartic [START_REF] Abbudi | On the existence of interfacial (Stoneley) waves in bonded piezoelectric half-spaces[END_REF] in its canonical form,

p 4 + rp 2 + sp + t = 0, (52)
say, using the substitution q

= p -(1/2)(ǫ 12 /ǫ 22 )(c D 16 /c D 44 ). Then, D = -4(12t + r 2 ) 3 + (-72tr + 2r 3 + 27s 2 ) 2 . ( 53
)
Here the equation D = 0 turns out to be a sextic in the squared wave speed, according to Maple.

5 Bleustein-Gulyaev wave as h → 0

Shuvalov and Every [START_REF] Shuvalov | Some properties of surface acoustic waves in anisotropic-coated solids, studied by the impedance method[END_REF] show that a great variety of asymptotic behaviors arises for an interface wave on a coated half-space, depending on whether the layer or the substrate is "fast" or "slow", or "dense" or "light", depending on which mode is under consideration, and depending on several other factors.

When the thickness of the layer vanishes here, the asymptotic behavior of the Love wave in the layer/substrate structure is that of a shearhorizontal surface wave propagation over the piezoelectric substrate alone (the Bleustein-Gulyaev wave [START_REF] Bleustein | A new surface wave in piezoelectric materials[END_REF][START_REF] Yu | Electroacoustic surface waves in piezoelectric materials[END_REF]). For such a wave, with metalized boundary conditions, the vector ξ(0) takes the form

ξ(0) = [U 3 (0), 0, 0, d 3 (0)] t = U 3 (0)[1, 0, 0, α] t , (54)

M

(2) 14

M (2) 44      1 α + α αα   =   0 0 0   , (55)
a homogeneous linear system of three equations. Its solution is non-trivial when the determinant of the 3 × 3 matrix on the left hand-side is zero. The resulting secular equation is a cubic in ρv 2 . At θ = 0 and θ = 90 • , the secular equation factorizes into the product of a term linear in ρv 2 and a term quadratic in ρv 2 . In particular, at θ = 0 the quadratic is:

(ρv 2 -c55 ) 2 + c2 44 (ρv 2 -c55 )ǫ 11 -ẽ2 15 c44 ǫ22 + ẽ2 24 = 0. ( 56
)
Now for KbNO 3 , the material parameters of interest are [START_REF] Zgonik | Material constants of KNbO 3 relevant for electro-and acousto-optics[END_REF]: c44 = 7.43, c55 = 2.5 (10 10 N/m 2 ), ẽ24 = 11.7, ẽ15 = 5.16 (C/m 2 ), ǫ11 = 34ǫ 0 , ǫ22 = 780ǫ 0 (ǫ 0 = 8.854 × 10 -12 F/m), and ρ = 4630 kg/m 3 . Using these values, the corresponding parameters in the rotated coordinate system follow from [START_REF] Destrade | Elastic interface acoustic waves in twinned crystals[END_REF], and in turn, T , R, Q follow from [START_REF] Collet | Explicit secular equations for piezoacoustic surface waves: Shear-Horizontal modes[END_REF], N from ( 22)-( 23), and M (n) from ( 27). Then the cubic secular equation is solved for v for any value of the cut angle θ. Out of the three possible roots, only one may correspond to the Bleustein-Gulyaev wave (Ref. [START_REF] Collet | Explicit secular equations for piezoacoustic surface waves: Shear-Horizontal modes[END_REF] explains how the adequate root is selected.) It turns out that the wave exists for all angles, with a speed v BG (say) increasing from 2895.35 m/s at θ = 0 to 4450.85 m/s at θ = 90 • . Fig. 2 shows the dependence in θ, and is in agreement with the plot obtained by Nakamura and Oshiki [START_REF] Nakamura | Theoretical analysis of horizontal shear mode piezoelectric surface acoustic waves in potassium niobate[END_REF] and by Mozhaev and Weihnacht [START_REF] Mozhaev | Incredible negative values of effective electromechanical coupling coefficient for surface acoustic waves in piezoelectrics[END_REF].

Fig. 2 also displays the speed of the bulk shear wave in a layer made of germanium, for which [START_REF]Springer Handbook of Condensed Matter and Materials Data[END_REF]: ĉ44 = 67.1 × 10 10 N/m 2 , ǫ11 = 16.6ǫ 0 , and ρ = 5330 kg/m 3 ; here v = 3550.31 m/s. The angle at which v BG = v is θ 0 = 50.0817 • . Finally, the variation of the limiting speed v L in KbNO 3 with the angle of cut is shown as well, and that plot is also in agreement with Mozhaev and Weihnacht [START_REF] Mozhaev | Incredible negative values of effective electromechanical coupling coefficient for surface acoustic waves in piezoelectrics[END_REF].

6 Dispersion curves 6.1 Special case θ = 0 • At θ = 0 • and kh = 0, the interface wave travels in the substrate alone, with the Bleustein-Gulyaev wave speed of 2895.35 m/s. The limiting speed (found here as root of 2 √ P = S) is v L = 3939.33 m/s. The speed of the fundamental mode starts at the Bleustein-Gulyaev wave speed at kh = 0, increases to a maximum speed of about 3857.18 m/s at kh = 0.517, and then decreases toward the shear bulk speed of the layer, v = 3550.31 m/s. In the narrow range where kh is smaller than 7.08436 ×10 -2 , the speed v of the fundamental mode wave is smaller than v and thus ĉ is given by (11) 2 ; otherwise it is given by (11) 1 . The fundamental mode exists for all values of kh.

The speeds of the subsequent modes start from the limiting speed v L at (kh) L and tend to v in a monotone decreasing manner. The cut-off parameter (kh) L for the first mode, second mode , and third mode is: 6.24, 12.75, and 19.27, respectively. Fig. 3 shows the dispersion curves for the fundamental mode and for the first mode. Qualitatively, the plots echo those of Kielczynski et al. [START_REF] Kielczynski | Shear-horizontal surface waves on piezoelectric ceramics with depolarized surface layer[END_REF] who considered a structure made of a 6mm substrate covered with a "depolarized" layer.

6.2 Special case θ = 90 • At θ = 90 • and kh = 0 the limiting speed is root of P = 0; here v L = 4508.73 m/s. The speed of the fundamental mode starts at the Bleustein-Gulyaev wave speed of 4450.85 m/s at kh = 0, increases to a maximum speed of about 4480.69 m/s at kh = 0.0581, and then decreases toward the shear bulk speed of the layer, v = 3550.31 m/s. Here too, the fundamental mode exists for all values of kh.

The speeds of the subsequent modes start from the limiting speed v L at (kh) L and decrease toward v. The cut-off parameter (kh) L for the first mode, second mode , and third mode is: 4.05, 8.06, and 12.06, respectively. Fig. 3 shows the dispersion curves for the fundamental mode and for the first and second modes.

Special case θ = θ 0

As an example of a rotated cut, consider the case where the speed of the Bleustein-Gulyaev wave in the substrate is equal to the shear wave speed in the layer; this occurs at θ = θ 0 = 50.08 • , see Fig. 2. For this cut the limiting speed is v L = 3813.36 m/s. Starting from v BG = v at kh = 0, the speed of the fundamental mode increases rapidly with kh. At kh = 0.1131, v = v L and the wave ceases to exists. This state of affairs continues until kh reaches 0.9410, after which the wave exists and its speed decreases toward v. Hence a forbidden band of frequencies emerges for the fundamental mode, in clear contrast with the situation for non-rotated cuts. Note that the dispersion equation (50) actually gives roots below v L in that range, which must nevertheless be discarded as they do not satisfy the exact boundary condition (38).

Here the first mode starts at the cut-off parameter (kh) L of: 8.946 with v L and then decreases toward v. Fig. 4 provides a zoom into the dispersion curves of the fundamental mode around the forbidden band and of the first mode up to kh = 12.

In this example, the layered structure supports a shear-horizontal wave which in the long and short wavelength ranges travels with the speed of the layer's bulk shear wave; in the intermediate range, the wave either does not exist, or travels at a greater speed.

Fig. 1 :

 1 Fig. 1: Geometry of the layered structure

where α := d 3

 3 (0)/U 3 (0) is complex. The fundamental equations (27), written for n = -1, 1, 2 (say), can be arranged as

Fig. 2 :

 2 Fig. 2: Solid curves: Variations of the Bleustein-Gulyaev wave and limiting wave speeds with cut angle in a homogeneous KbNO 3 substrate. Dashed line: Speed of a (bulk) shear wave in germanium.

Fig. 3 :Fig. 4 :

 34 Fig. 3: Dispersion curves at θ = 0 • , θ 0 , 90 • .

  S = c 44 ǫ 11 + c 55 ǫ 22 + 2e 15 e 24 -ǫ 22 (ρv 2 ) c 44 ǫ 22 + e 2

			,
		24
	P =	c 55 ǫ 11 + e 2 15 -ǫ 11 (ρv 2 ) 24 c 44 ǫ 22 + e 2	.

  lead to the following cut-off equation,

	c 44 +	e 2 24 ǫ 22 e 24 q 2 + ǫ 11 ǫ 22 -e 15 c 44 ǫ 11 + e 24 e 15 ǫ 22	ǫ	ǫ 11 ǫ 22 ǫ 22 q 2 + ĉ q 2 + e 24 ǫ 22 e 15	= 0.	(42)


Delphine Lejri

Arnaud Can

Nicole Schiper

Ludovic Leclercq

Accounting for traffic speed dynamics when calculating COPERT and PHEM pollutant emissions at the urban scale

Keywords: Traffic microsimulation, network emission modeling, mean speed definition, vehicle kinematics representation, comparison of COPERT and PHEM

come

Introduction

Road traffic emissions have dramatic local and global effects. Pollutants such as Nitrogen Oxides (NO x ) and Particulate Matter (PM) have known detrimental impacts on human health, including respiratory and cardiovascular diseases [START_REF] Shaughnessy | Health effects of ambient levels of respirable particulate matter (PM) on healthy, young-adult population[END_REF], while Carbon Dioxide (CO 2 ) emissions greatly contribute to global warming. Accurate trafficrelated emission estimations are thus crucial to assess their evolutions over several years, and to quantify the environmental impact of sustainable transportation policies such as lowemission zones and traffic regulation strategies. Emission estimations should be based on models developed for the given scale of interest, which can range from very local (e.g., for the assessment of certain road traffic facilities) to global investigations (e.g., for compiling inventories), and therefore rely on both relevant representations of road traffic and vehicle emissions. A detailed review of vehicle emission models can be found in [START_REF] Smit | Validation of road vehicle and traffic emission models -A review and meta-analysis[END_REF], while [START_REF] Fallah Shorshani | Modelling chain for the effect of road traffic on air and water quality: Techniques, current status and future prospects[END_REF] provide a review of the complete modeling chain (traffic, emission, dispersion and storm water).

The main current research efforts focus on the urban scale, because urban road traffic causes the vehicle kinematics that generate the highest emissions, namely rapid speed variations and congestions, and they are the most difficult to take into account [START_REF] Ma | Effects of driving style on the fuel consumption of city buses under different road conditions and vehicle masses[END_REF] [START_REF] Ahn | A field evaluation case study of the environmental and energy impacts of traffic calming[END_REF], [START_REF] De | Environmental effects of driving behaviour and congestion related to passenger cars[END_REF] [START_REF] Zhang | Vehicle emissions in congestion: Comparison of work zone, rush hour and free-flow conditions[END_REF] [START_REF] Qu | Multivariate analysis between driving condition and vehicle emission for light duty gasoline vehicles during rush hours[END_REF]. Models based on vehicle representations have been developed to overcome these problems. Microscopic traffic simulation provides detailed vehicular kinematics, namely an estimate of the 1s-evolution of speed and acceleration for each vehicle on the network. Thus it can capture congestion effects at the finest spatial and temporal resolution. The estimated traffic data are provided for an instantaneous emission model which uses speed profiles [START_REF] Chen | Microscopic Traffic-Emission Simulation and Case Study for Evaluation of Traffic Control Strategies[END_REF] or derived indicators to estimate emissions [START_REF] Frey | Vehicle-Specific Emissions Modeling Based upon on-Road Measurements[END_REF]Jiménezpalacios, 1999). Such modeling has been advantageously used for assessing traffic regulation strategies such as traffic signal synchronization and speed reduction [START_REF] Madireddy | Assessment of the impact of speed limit reduction and traffic signal coordination on vehicle emissions using an integrated approach[END_REF]. However, biases introduced by this coupling have been highlighted, especially concerning simplified vehicle trajectories (Vieira da [START_REF] Vieira Da Rocha | Are vehicle trajectories simulated by dynamic traffic models relevant for estimating fuel consumption?[END_REF] and calibration processes [START_REF] Jie | Calibration of a microscopic simulation model for emission calculation[END_REF] [START_REF] Xu | Estimating Project-Level Vehicle Emissions with VISSIM and MOVES-Matrix[END_REF] [START_REF] Lu | Transportation Research Board Annual Meeting and Publication in the Transportation Research Record Washington[END_REF]. Integrating emissions in the optimization function can help to reduce these biases [START_REF] Vieira Da Rocha | Does traffic-related calibration of car-following models provide accurate estimations of vehicle emissions?[END_REF].

Microscopic emission models entail long computation times, which confine them to local investigations. Aggregated models have been developed for broader scale assessments. These models rely on simplified kinematics variables, such as vehicle mean speed (Samaras and Geivanidis, 2005), which is the traffic parameter that most influences emissions [START_REF] Hansen | The influence of driving patterns on petrol passenger car emissions[END_REF] [START_REF] Ericsson | Independent driving pattern factors and their influence on fuel-use and exhaust emission factors[END_REF] [START_REF] Joumard | Infuence of driving cycles on unit emissions from passenger cars[END_REF] [START_REF] André | Driving speeds in Europe for pollutant emissions estimation[END_REF] [START_REF] André | Analysis and modelling of the pollutant emissions from European cars regarding the driving characteristics and test cycles[END_REF]. Models based on mean driving speeds are less successful with congested situations (R. [START_REF] Smit | Do air pollution emissions and fuel consumption models for roadways include the effects of congestion in the roadway traffic flow?[END_REF]. Indeed, these models implicitly take varying levels of congestion into account, depending on the average speed and the road type chosen (and the driving cycle associated). But they have not been developed to produce accurate local emission predictions. Several alternatives have been proposed to overcome this drawback. (Robin Smit et al., 2008) reconstructed a speed distribution from the available mean speed on road sections. [START_REF] Pitsiava-Latinopoulou | Development of a system of environmental and energy consumption data for urban road traffic, Pilot application in Thessaloniki, Greece[END_REF][START_REF] Samaras | COPERT Micro : a tool to calculate the vehicle emissions in urban areas[END_REF] proposed using aggregated models at spatial scales smaller than the scale to which they were initially dedicated, up to the road segment scale. However, this calls into question the temporal and spatial scales at which the modeling choices made remain relevant. [START_REF] Samaras | Improving fuel consumption and CO2emissions calculations in urban areas by coupling a dynamic micro traffic model with an instantaneous emissions model[END_REF] proposed an evaluation for a 1.6 km long urban corridor based on micro-scale data (models and measurements). The authors showed that average speed is strongly correlated with congestion. Therefore, the aggregate emissions model (COPERT 5) is congestion sensitive and more relevant for depicting loaded rather than fluid situations and their effects on EURO 5 diesel vehicle fuel consumption.

In essence, road sections and short time scales (traditionally 6 min periods, with electromagnetic loop measurements) are relevant for observing traffic on an urban network. This scale provides information on vehicle dynamics that characterize congestion through aggregate traffic variables and raises questions regarding the emission model that should be interfaced with these traffic data at this scale of interest. Therefore, it seems crucial to highlight emission-modeling representations that are compatible with the outputs of the traffic microsimulation.

The key objectives of this article are: (i) to evaluate the impact of traffic information processing on the emissions obtained with an aggregated model, and (ii) to compare the emissions over the network obtained with an aggregated emission model and an instantaneous one. Section 2 presents the simulation framework and the models tested. Section 3 presents alternatives to usual aggregated models. Sections 4 and 5 compare the emission estimations for a case study. Finally, section 6 concludes on modeling choices for urban emission assessment.

Material

Coupling traffic variable definitions with emission models is tested by comparing estimated fuel consumptions and NOx emissions over a simulation network. The tested emission models are COPERT and PHEM, which rely on average travel speeds and instantaneous individual speeds, respectively. The versions of COPERT and PHEM used are COPERT IV and PHEM v11, respectively (see section 2.2). The traffic data required are provided by the dynamic traffic platform Symuvia, which can supply input data from individual vehicle trajectories to aggregated speeds for each vehicle category at given locations on the network and at any temporal resolution.

Traffic simulation

Traffic model and simulation

The network studied is a 3 km 2 zone covering part of the municipalities of Le Perreux-sur-Marne and Neuilly-Plaisance in the Paris region. The traffic network, displayed in Figure 1, was selected for its wide range of traffic conditions. Two structural arterials, "Boulevard d'Alsace-Lorraine" and "rue Pasteur", both oriented in the east-west direction, cross the site and feed it with dense traffic. The rest of the site is composed of a few main corridors (shown in blue in Figure 1), and a large number of residential streets crossed by light volumes of traffic limited in speed to 30 km/h. The traffic microsimulation was implemented in the Symuvia platform1 , which gives access to the position, speed and acceleration of each vehicle on the network with a 1s-resolution. Vehicle routing choices are governed by a dynamic traffic assignment model, which guides each vehicle in the network on the route that minimizes its travel time to its initially assigned destination. Vehicle movements at the microscopic scale are governed by a set of rules, including car-following modeling (Leclercq, 2007a(Leclercq, , 2007b)), lane-changes [START_REF] Laval | Microscopic modeling of the relaxation phenomenon using a macroscopic lane-changing model[END_REF] and specific movements at intersections [START_REF] Chevallier | A macroscopic theory for unsignalized intersections[END_REF]. The platform also copes with the cohabitation on the network of vehicles with different kinematics, including passenger cars, buses and heavy-duty vehicles. The question of using the platform outputs for pollutant emission estimations was addressed in (Vieira da [START_REF] Vieira Da Rocha | Are vehicle trajectories simulated by dynamic traffic models relevant for estimating fuel consumption?[END_REF]. The simulation consists of 2.5 hours representing the morning rush hour. The Origin-Destination matrix was calibrated with hourly traffic flow rates measured on typical weekdays at 24 locations of the network (see Figure 1). The calibration results are depicted in

Traffic outputs

The most refined available traffic outputs provided by Symuvia are the vehicle speeds and accelerations at each time-step (1 s), for all the vehicles on the network. This traffic model has the particular advantage of avoiding the unrealistic accelerations that can be found in traffic microsimulation (Leclercq, 2007b). Indeed, the traffic simulator used -Symuviaproduces simplified vehicle trajectories. However, it allows conforming to realistic speed and acceleration values, depending on the type of vehicle. Thus, each vehicle type is associated with a maximum acceleration per speed range. For example, between 20 and 30 km/h, light vehicles can accelerate by a maximum of 1 m/s 2 , while heavy vehicles cannot accelerate by more than 0.6 m/s 2 . Above 30km/h, these limits decrease by 0.5 m/s 2 , 0.3 m/s 2 respectively.

Figure 3 shows the speed-acceleration pairs experimented during the simulation under congested conditions for light and heavy vehicles. These vehicle trajectories can be used directly to calculate emissions with PHEM (including a default gearshift model), or aggregated into traffic variables in order to correspond to the required COPERT inputs (see section 2.2). In both cases, the traffic outputs will be provided for each road segment with a 6 min-resolution, and will be interfaced at that scale with emission models. When a traffic microsimulation is available, we can characterize the spatial mean speeds most relevant for evaluating pollutant emissions. However, it is also possible to reproduce the case where the level of traffic information is less accurate and where the average speeds are characterized by electromagnetic loops. In this case, based on punctual mean speeds, emissions can be assessed as can the deviations from a more comprehensive traffic description observed in the same conditions. In this paper, various speed definitions are compared for qualifying vehicle kinematics on each road segment, thus leading to an increasing level of detail:

1. The operational (or default) definition assumes that the average speed is the speed limit V limit , which is the first operational information available. 2. The speed experimented at one specific location on the road segment V punctual corresponds to the local measure performed for instance by electromagnetic loops. Three virtual positions are tested, corresponding to the following locations: 25%, 50% and 75% of each segment length, from the beginning of the road segment.

Associating the speed on a road segment to V punctual amounts to assuming that vehicle speeds are homogeneous along the segment. 3. The speed characterizing the vehicle kinematics on the whole road segment V spatial can be determined using Edie's definition [START_REF] Edie | Discussion of traffic stream measurements and definitions[END_REF], in which the spatial speed is the ratio between the total travel distance and the total spent time. This speed definition is the most accurate and compatible with the emission estimations, but unfortunately it relies on data not available on a real network.

These three speed definitions can differ significantly, in particular under congestion, as shown in Figure 4. As expected, the speed limit V limit overestimates the actual speeds. The wider the vehicle mean speed distribution on the road segment, the larger the errors. Punctual loops also result in speed overestimations, a long acknowledged bias determined by the location of the loop. Indeed, at the urban scale, loops close to signals are subject to much more congestion. The impact of these different speed definitions on the resulting emission errors is investigated in section 4.

Emission models

These traffic variables are used by the emission models to estimate fuel consumptions (FC) and NOx emissions in the network. The two types of emission model investigated are the modal model PHEMv11 and the aggregated model COPERT IV. Only hot exhaust emissions are considered. Indeed, although cold-start and evaporative emissions can amount to a nonnegligible proportion of total emissions, it can be assumed that both the potentially initiated traffic strategies and the vehicle kinematics reproduced affect them in the same way as they affect hot exhaust emissions.

Vehicle emissions strongly depend on the vehicle's characteristics, namely the vehicle category (passenger cars, buses, etc.) and, within each category, fuel type, age, technology, and emission standard (Euro 4, Euro 5, etc.) of the vehicle. The vehicle fleet is one of the parameters set throughout the study. This fleet is the French urban fleet for the year 2015 obtained from the IFSTTAR fleet updated in 2013. This passenger car fleet is composed of 30% EURO 5 diesel vehicles and 24% EURO 4 diesel vehicles, whereas the light commercial vehicle fleet is composed of 45% EURO 5 diesel vehicles and 30% EURO 4 diesel vehicles. The vehicle classes are divided as follows to fit the traffic microsimulation: 82% passenger cars, 14% light commercial vehicles, 3% heavy duty vehicles and 1% urban buses. The interpretation of this car fleet by both modeling approaches differs, which can be a source of discrepancies between the emissions estimated.

COPERT IV

COPERT IV has been widely used in most European Countries for compiling national emission inventories [START_REF] Ntziachristos | COPERT: A European Road Transport Emission Inventory Model[END_REF], but it is also increasingly used for emission modeling at the street level [START_REF] Borge | Comparison of road traffic emission models in Madrid (Spain)[END_REF]. Concretely, this requires relying on the average speeds and travel distances provided by a traffic model or measurements at each period, for each road segment and vehicle category, and then applying the COPERT speed continuous functions over each of these road segments. However, this use at spatial scales lower than the driving cycles is subject to question. Indeed, the actual road segment speed distributions might differ from those of the driving cycles, and lose representativeness for very small samples or specific traffic conditions (e.g. in the vicinity of intersections), yielding inaccurate emission estimates. COPERT relies on mean driving speed and travel distance (total distance in km travelled per vehicle for a given period) to predict the related exhaust emissions. The total emissions are calculated as the product of the travel distance and the unitary emission factors, according to formula (1). Unitary emission factors consist of speed continuous functions that have been constructed over driving cycles of about 6mn-length, which are representative of the traffic conditions encountered [START_REF] André | Analysis and modelling of the pollutant emissions from European cars regarding the driving characteristics and test cycles[END_REF]. These unitary emission functions are defined for each pollutant k and each vehicle class 𝑐 (e.g. passenger cars, light duty vehicles, buses and heavy duty vehicles). The unitary emission function of a specific class is obtained by using a weighted average of the vehicle technologies that compose the class. Consequently, the emissions E !,! (g) are calculated as:

E !,! = D ! . F !,! V !
(1) with 𝑭 !,! being the unitary emission factor (g/km) of pollutant k and vehicle class c , 𝐷 ! the total travel distance (km) and 𝑉 ! the mean speed for vehicle class c .

Thus, within a vehicle class, the model considers the same mean vehicle (obtained for a fixed car fleet) whatever the period and the road segment size. The model then omits the dispersion in emissions due to the variability in the actual car fleet observed for small periods and road segments, i.e. when the number of observed vehicles is small.

PHEM v11

PHEM (Passenger Car and Heavy Duty Emission Model) calculates the fuel consumption and emissions of vehicles with a 1s-time resolution, based on their longitudinal dynamics and on engine emission maps [START_REF] Hausberger | Emission factors for heavy-duty vehicles and validation by tunnel measurements[END_REF][START_REF] Hausberger | Emission factors for heavy-duty vehicles and validation by tunnel measurements[END_REF]. The model first provides an estimate of the 1s-engine power of a vehicle, on the basis of its speed time-series and the road gradient. The engine speed is estimated using transmission ratios and a gear shift model. Engine emission maps then allow estimating the time evolution of fuel consumptions and airborne pollutant and particulate matter emissions. The model also includes transient correction functions, and a cold start tool. Cold start emissions will however be disregarded to enable comparison with COPERT. PHEM is thus appropriate for coupling with dynamic traffic platforms intended to provide vehicle trajectories. Such coupling has been performed on several occasions to test the impact on emissions of road traffic strategies that modify vehicle kinematics behavior [START_REF] Zallinger | EVALUATION OF A COUPLED MICRO-SCOPIC TRAFFIC SIMULATOR AND INSTANTANEOUS EMISSION MODEL[END_REF]. However, the inadequacy between the high traffic data resolution required and the available dynamic traffic model outputs, which are much less refined or subject to inaccuracies, can lead to significant discrepancies (Vieira da [START_REF] Vieira Da Rocha | Does traffic-related calibration of car-following models provide accurate estimations of vehicle emissions?[END_REF].

In PHEM, each individual vehicle of the simulation is associated with one specific vehicle that composes the car fleet, with its own vehicle technology. Thus, in theory the model copes with the small-scale variations of the actual car fleets observed; consequently, periods with similar vehicle kinematics conditions can result in different emissions. The fit between the actual and the modeled emission variability is, however, a key factor.

3 Applying COPERT at the road segment scale COPERT emission factors are continuous and non-linear functions of vehicle average speeds. Therefore, their application at the road segment scale and with small periods first requires evaluating the sensitivity of the model to mean speed variations. Second, it requires adaptations to handle the mean speeds below 10 km/h that are frequent at small time and space scales, and cope with the vehicle trajectories provided by the dynamic traffic simulation. Adaptations to COPERT emission functions are proposed in this section to overcome these issues where appropriate.

Sensitivity of COPERT to mean speed

This section quantifies the impact of variations in the estimated mean speeds on the emission factors estimated within COPERT, thanks to formula (1). Mean speeds are relatively low under urban driving conditions, and emission factors are highly variable within this speed range [START_REF] Ntziachristos | COPERT: A European Road Transport Emission Inventory Model[END_REF]. If we hypothesize that the travel distance is associated with small variations whose effects can be neglected, then the relative errors on emissions are only due to the bias on mean speed and can be expressed as:

∆! ! = !! !! ! ! ∆! ! (2)
These relative emission errors in terms of speed errors are depicted for FC and NOx emissions in Figure 5. The figure highlights the low accuracy of speed estimates needed for actual speeds around 70 km/h, which can be explained by the known low variability of emissions in this speed range. Conversely, the need for accurate speed estimates increases at low speeds.

For example, at 30 km/h, a 2 km/h error on the mean speed estimates will lead to a 3.1% error on the FC estimate and a 3.9% error on the NOx emission estimate. At the same time, a 10km/h error will lead to a 15.7% error on the FC estimate and a 19.3% error on NOx emissions. These percentages correspond to errors on the mean speed values commonly reached through both modeling and in situ measurements. Thus, at the urban scale, such levels of bias are normally expected, especially at small time and spatial scales.

Adapting the COPERT model

Coping with low mean driving speed within COPERT

In Figure 6, the FC and NOx emission factors F k (in g/km) are drawn in terms of the mean speed. The COPERT model handles only mean speeds higher than 10 km/h, although mean speeds below 10 km/h can be encountered at the road segment scale for some of the 6minperiods. As emissions are definitely not insignificant at that speed range, the COPERT emission curves are extended by either maintaining the EF value at 10 km/h (straight extension), or drawing out the EF curve fitted with a 4 th degree polynomial (polynomial extension). The consequences of these choices are discussed in section 4.

Instantaneous mean driving speed

The dynamic traffic simulation makes it possible to describe the vehicle kinematics associated with a road segment and a period in a more detailed way than a mean driving speed. The distribution of time spent and the distribution of the travel distances per speed class are, for instance, available, possibly refining the emissions calculated with COPERT. However, the COPERT emission curves are not designed to accept such fine input data. The COPERT model is adapted in this section to make it compatible with high-frequency traffic data inputs. These adaptations consist in: (i) treating idling time separately, as it is a specific part of the speed class distribution, and (ii) rewriting the emission factor functions to reveal the dependency on instantaneous speeds rather than on spatial mean speeds. The particular case of idling time can represent from 23% to 42% of time spent on network.

The PHEM model provided the values of idling emissions for each pollutant and vehicle category. For instance, for passenger cars, idling FC is set at 0.17 g/s whereas NOx emissions are set at 0.0011 g/s.

When considering an average speed as a distribution of instantaneous speeds around this average, a mathematical formulation of the emission curves as a function of the instantaneous speeds is determined using a convolution. However, using this new expression, emissions are always calculated at the scale of a vehicle flow and not at the individual scale. In accordance with the spirit of the COPERT methodology, emission factors are defined for each vehicle category by integrating the fleet distribution within that category. COPERT emissions functions are therefore fitted and transformed into suitable functions with instantaneous speeds. We consider that each mean speed V stands for a range of experimental instantaneous speed distributed in the interval [𝑉 -𝛼; 𝑉 + 𝛼], which follows a uniform distribution 𝑑 ! , with 𝜎 ! = 𝛼 ! /3. Thus, the emission associated with a mean speed V can be seen as a mixture of emission levels experimented for the different instantaneous speed classes. The emission functions 𝑭 !,! depending on mean speed 𝑉 are rewritten as follows:

𝑭 !,! 𝑉 = 𝑮 !,! 𝑤 . 𝒅 ! 𝑑𝑤 !!! !!! (3)
This formula helps to define the instantaneous emission functions 𝑮 !,! of interest.

Considering that 𝑮 !,! can be adjusted with a 4 th degree polynomial, the first result is that the functions 𝑭 !,! maintain the same polynomial form.

With G !,! w = p ! !,! w ! + p ! !,! w ! + p ! !,! w ! + p ! !,! w + p ! !,!
, equation (3) becomes:

𝑭 !,! 𝑉 = (𝑝 ! !,! 𝑤 ! + 𝑝 ! !,! 𝑤 ! + 𝑝 ! !,! 𝑤 ! + 𝑝 ! !,! 𝑤 + 𝑝 ! !,! ). 1 2𝛼 𝑑𝑤 !!! !!!
= 1 2𝛼 𝑝 ! !,! 𝑤 ! 5 + 𝑝 ! !,! 𝑤 ! 4 + 𝑝 ! !,! 𝑤 ! 3 + 𝑝 ! !,! 𝑤 ! 2 + 𝑝 ! !,! 𝑤 !!! !!! = 𝑝 ! !,! 𝑉 ! + 𝑝 ! !,! 𝑉 ! + (𝑝 ! !,! + 6𝜎 ! 𝑝 ! !,! ) 𝑉 ! + (𝑝 ! !,! + 3𝜎 ! 𝑝 ! !,! ) 𝑉 + (𝑝 ! !,! + 𝜎 ! 𝑝 ! !,! + 9 5 𝜎 ! 𝑝 ! !,! )
Furthermore, considering 3 rd order polynomials, the expression of functions 𝑭 !,! is the same with a uniform and a normal distribution. With 4 th order polynomials, the new curve is formulated with both distributions as a fourth order polynomial but the coefficients differ. In order to provide a clear presentation of the concept, here we explain only the formulation for a 4 th order polynomial and a uniform distribution.

The functions 𝑭 !,! are fitted on existing COPERT curves for each vehicle class c and each pollutant 𝑘. The fitting process provides numerical values of polynomial coefficients and α. The latter is estimated between 8 and 13. This is in accordance with the speed standard deviation of urban driving cycles.

The 𝑮 !,! functions are then derived from the coefficients obtained. These functions, depicted for FC and NOx emissions in Figure 8, resemble the initial functions, but are slightly lower in the speed range of interest. The adapted COPERT emission functions will be further coupled to simulated traffic data in terms of time spent per speed classes. This innovative way of implementing emission functions is derived from equation ( 1):

𝐸 !,! = T ! ! ! . 𝑮 !,! 𝑉 ! (4)
with 𝑮 !,! the "instantaneous" unitary emission factor (g/s) of pollutant k and vehicle class 𝑐 𝑇 ! the total time spent by vehicle class 𝑐 (s) 𝑉 ! the speed characterizing the i th speed class (km/h) Here, 2km/h speed classes were selected and 1s vehicle speeds were reformulated with the total time spent in each speed class. The corresponding network emissions are analyzed in section 4.2.

Sensitivity of COPERT to the speed definitions for emissions calculated at the network scale

Sensitivity to speed definition

The FC and NOx emissions are calculated at each road segment every 6 min by formula (1), and then summed to observe the impact of speed definition on the emissions calculated over the network. The calculated emissions obtained with the spatial mean speed serve as references, as this speed definition is the most accurate that can be encountered. Figure 9 represents the network emissions for each 6 min period, obtained with COPERT (straight extension for the lowest speeds) and the various speed definitions. Degraded speed definitions (speed limit and punctual loops) lead globally to underestimated emissions. The discrepancies between the calculated emissions also depend on the period, and reach their maximum value during congestion (grey area). At the network scale, the positions of the virtual loops do not have a significant impact on the estimated emissions. Thus, only virtual loops positioned in the middle of the road segment (the 50% loop) are kept in the rest of this article. Using the loop detectors to estimate mean speeds introduces discrepancies on the global emissions at the 6 min period scale that range from -9.7% to -13.4% for FC, and from -13.5% to -17.1% for NOx, compared with the global emissions obtained with the spatial mean speed definition. Using the speed limit to estimate mean speeds introduces an emission bias at the 6 min period scale that ranges from -19.8% to -25.3% for FC, and from -30.7% to -36.0% for NOx. represents the spatial distribution of the discrepancies on emissions calculated at the road segment scale for a specific congested period, according to the speed definition. The local errors are higher than the global ones, especially under congestion. For NOx emissions, the mean relative error using the speed limit ranges locally from -67.2% to +4.6%, and from -64.5% to 73.0%, using the punctual mean speed. The punctual mean speed definition yields an almost null bias for 47.8% of the road segments for this period. Indeed, almost a third of these road segments do not significantly contribute to the network emissions. This can be explained by the heterogeneity of the network (in terms of geometry and/or traffic). More specifically, Figure 10 (b) represents the spatial distribution in terms of percentage of total emissions. Using a speed limit as the mean speed for emission estimation leads to a nearly null local error in 28% of the road segments, with represents around 3% of the total emissions. On the contrary, when it comes to small sets of road segments associated with huge local gaps (around -60%), these can contribute to 13% of the total emissions.

The spatial analysis reveals that locally the bias can be dramatically higher. However, it also shows that the necessary accuracy on mean speed can be considered in terms of emission stakes and is not homogeneous over the road segments. For further analyses, the mean speed definition chosen will be the spatial mean speed.

Sensitivity to the COPERT adaptations

Reference scenario

The COPERT adaptations proposed in section 3.2 are compared in this section for emissions based on the vehicle trajectories provided by the dynamic traffic simulation. Figure 11 shows the discrepancies in emissions calculated over the network for the three COPERT adaptations for free flow and congested periods (grey area).

Figure 11 COPERT network FC (left) and NOx emissions (right) evaluated thanks using three adapted emission functions (see section 3.2).

Considering the mean speed approach, using a polynomial extension of COPERT emission functions implicitly results in integrating surplus emissions associated with very low speeds. This is confirmed at the network level: FC (and NOx emission) is therefore lightly overestimated, around 1.7% (2.6% for NOx) and at most 3% (4.6% for NOx), compared with the straight extension of emission functions that serve as references. It is interesting to note that these two strategies have little impact on total emissions/consumption.

On the other hand, using adapted COPERT emission functions introduces a slight underestimation in free flow conditions and an overestimation compared with the straight extension of emission curves. This approach, treating traffic information as time spent by instantaneous speed classes, results in discrepancies at the network scale that range from -1.1% to 12.8% for FC (from -1.5% to 15.9% for NOx emissions), that is to say 4.8% (5.4% for NOx) for the whole period of 2.5 hours.

Increased demand scenario

A congested scenario is implemented, with a 30% increase in demand, to highlight the discrepancy between the models under congestion. In these conditions, the gaps associated with the use of polynomial extended COPERT emission functions (vs straight extended EF), are slightly accentuated and reach a maximum of 7% for FC and 8.8% for NOx emissions. Thus, the discrepancies between mean speed and speed distribution approaches are also increased (from -0.8% to 43.2% for FC and from -1.3% to 45.2% for NOx emission). The total gap for the whole period reaches a 15.3% discrepancy for the FC estimate and a 15.2% discrepancy for NOx emissions.

The gap between the two approaches is revealed in Figure 13 (speed distribution versus mean speed) by comparing the NOx emission contribution of each speed class and both scenarios. In Figure 13.a, the 2km/h bins are defined on instantaneous speed and integrated in the calculation as time spent in the speed class considered. In Figure 13.b, the 2km/h bins are defined for a 6 min mean speed and integrated in the calculation as distance travelled at the speed class considered. These distributions show the major contribution of emission associated with idling time, which is treated individually in the first case but mixed in the 6 min mean fleet speed. This effect is even more significant for the increased demand scenario, where the vehicles under congestion are stopped for 37.6% of the time, versus 26.2% for the reference scenario. The use of speed distributions provided by the dynamic traffic simulation leads to higher FC and NOx emissions, especially during congested periods (grey area). This traffic data treatment is certainly associated with better integration of low speeds and speed variability. These traffic conditions also correspond to the most emissive conditions. This specific trafficemission coupling seems to evaluate real emission conditions more precisely at the network scale, although this is not validated experimentally here. In the absence of validation, the emission levels could be compared to emissions evaluated using an instantaneous emission model such as PHEM.

Comparison of emissions models

This section is devoted to emission calculations at the road segment scale, using the PHEM and COPERT models. The general objective of this paper is to use a traffic microsimulation to evaluate pollutant emissions and to analyze the differences found between the implementation of different methodologies. The hypothesis tested is that the aggregate COPERT model, based on average speed, is not capable of integrating all the speed variability, which is higher during congestion. On the contrary, we assume that the PHEM model, which relies on vehicle kinematics, is able to integrate these typical phenomena. The differences observed are therefore interpreted in relation to traffic dynamics, although the absolute differences between the models may be due to other origins. This is why the analysis will not conclude on the best performance of one model compared to another, but rather aim to assess their sensitivity to traffic dynamics.

PHEM emissions

As the microsimulation provides very detailed information for each vehicle (1s trajectories), the temptation is great to use an emission model able to maintain this level of detail. Despite the huge amount of data to be implemented (associated with approximately 2 hours calculation for one run), the traffic data is compatible with the input data of the PHEM modal emissions model . For the following analysis, the emission evaluated with PHEM will be represented as the mean value of the replications for each 6 min time period.

PHEM and COPERT emissions

Finally, three modeling chains, starting from the same traffic microsimulation, were compared at the network scale in terms of fuel consumption and NOx emissions. Individual traffic trajectories were coupled with PHEM, whereas COPERT was implemented with aggregated traffic variables (spatial mean speed) and speed distribution (see section 3.2). All the environmental assessments were produced for a given fleet transcribed in the emissions models. At the driving cycle scale, PHEM overestimates overall consumption compared to COPERT in congested conditions, but the trend is less marked for NOx emissions. The issue of model consistency encountered by emission modelers is not addressed in this paper. The authors deliberately chose not to artificially modify the parameters of the emission models to make them consistent. The analysis below does not focus on the nominal values proposed by the tools but on their evolution according to traffic conditions. To do this, the emissions normalized by the average value over the period were plotted for both scenarios without involving the inherent differences due to the model versions in use (see Figures 16 & 18). As shown in Figure 15, the PHEM modal emission model obtains higher emission levels, especially for congestion periods. The gap between COPERT emissions (with spatial mean speed) compared to PHEM emissions is quantified: the relative deviations reach -34.6% for fuel consumption (-24.2% for NOx emissions). Using the adapted COPERT emission functions reduced the maximum relative discrepancies to -30.1% for fuel consumption (-16.5% for NOx emissions). This means that the result is highly dependent on the modeling chain. However, to overcome version-related deviations, the emission values were compared by model between fluid and congested periods for the two scenarios and between both scenarios (see Table 1). Figure 1 shows that all the models follow the overall traffic pattern over the period considered. However, it reveals that in a free-flowing situation, the average-speed approach tends to overestimate NOx emissions/consumption, whereas in congested conditions, it tends to underestimate them because it does not include the effects of congestion.

These differences have been quantified comparing the emission on (E on ) and off (E off ) peak for both scenarios and drawing the comparison between scenario 1 and 2, distinguished by the congestion rate:

δ = 𝐸 !" -𝐸 !"" 𝐸 !"!#$ x 100 Δ = 𝐸 !"!! -𝐸 !"!! 𝐸 !"!! + 𝐸 !"!! x 100
These indicators should be read as follows: the higher the value, the more the methodology integrates the traffic dynamics. The same conclusion with enhanced effects can be observed with the increased demand scenario Figure 17 and 18). Comparing COPERT to PHEM emissions, the relative deviations reach -46.8% for fuel consumption (-41.5% for NOx emissions). Using adapted COPERT emission functions reduces the maximum relative discrepancies to -32.1% for fuel consumption (-24.4% for NOx emissions).

Once again, this comparison says nothing about the more realistic emission level and these results cannot be seen as a recommendation for using one modeling chain rather than another. The purpose is to exhibit the potential deviations induced by different treatments of the same traffic microsimulation.

Discussion

In a context where transportation managers are urged to take the environmental impact of road traffic into account, the challenge is to provide an accurate dynamical estimation of the network emissions. Indeed, to qualify the existing situation or traffic regulation measures, officials are interested in accurate environmental evaluations. To do this, evaluation methods, even at the city scale, cannot ignore traffic dynamics, that is to say the congestion phenomenon. This point was stressed in particular in this article: traffic microsimulation provides a fine representation of vehicle kinematics, which is relevant for emission calculations. Indeed, when comparing emission calculations, each discrepancy observed, even at the network scale, depended on traffic conditions and was considerably greater during congestion.

When implementing emission calculations from vehicle trajectories, an instantaneous model, such as PHEM, seems to be more efficient. The strategy is definitely time-consuming but still affordable in terms of the amount of data processed. Nevertheless, this article specified the imprecision due to the fact that we do not know the real fleet composition: this can induce a non-negligible deviation, even at the network scale. With the simulated traffic data, the absolute global relative error can reach 5.0% for fuel consumption and 9.2% for NOx emissions.

Regarding aggregating emission models, such as COPERT, the first step is to define the useful aggregated traffic variables. This article specifically underlined the need for caution regarding the mean speed definition. The gap associated with a degraded definition was quantified. In the case of missing information, the road segment limit speed can occasionally be used to evaluate emissions, leading to a critical bias (-25% for FC, -30% for NOx).

Finally, another approach was proposed, improving traffic representation: traffic conditions were no longer synthesized with a mean speed but with a speed distribution. A theoretical reflection was proposed to adapt the COPERT emission functions to this input information. Consequently, the evaluation led to a higher emission value. The performance indicators were more or less doubled compared to a classical implementation and revealed that the integration of traffic dynamics was successful. This modeling chain appeared efficient at the urban scale in congested conditions.

This article focused on comparing emission models but cannot conclude on favoring one modeling chain over another. Its contribution lies in the description of the requisite traffic dynamics: congested periods should be characterized carefully. However, the issue of validation remains to be addressed. It would be interesting to compare these results to measurements, such as ambient concentrations, to obtain an idea of the range of real emissions over the network.

Figure 1 .

 1 Figure 1. Simulation network implemented in Symuvia. Left: road network and Origin/Destination traffic demands. Right: Mean weekday traffic counts for passenger cars and heavy duty vehicles

Figure 2 :

 2 Figure 2 : measured (red) and estimated (blue) traffic flow rates (veh/h) over time (hours) at 12 locations.

Figure 3

 3 Figure 3 Speed acceleration distributed under congested conditions for light (left) and heavy (right) vehicles

Figure 4 :

 4 Figure 4 : Left: extract of the vehicle trajectories provided by Symuvia. Right: Mean speed evolution according to its definition on the arterial "Boulevard d'Alsace-Lorraine ".

Figure 5

 5 Figure 5 Relative errors on fuel consumption and NOx emissions related to mean speed variations from 2 to 10 km/h

Figure 6

 6 Figure 6 Extended emission functions for low mean speed (left: FC function; right: NOx function)

Figure 7

 7 Figure 7 Travel distance per 2km/h speed classes over the network under freeflow (left) and congested (right) conditions

Figure 8

 8 Figure 8 Adapted COPERT emission functions: COPERT emission functions with straight extension for low speed, with the fitted COPERT emission functions permitting polynomial extension and the emission functions adapted to instantaneous speed (left: FC estimate; right: NOx estimate)

Figure 9 6

 9 Figure 9 6 min fuel consumption (left) and NOx emissions (right) over the network during the morning peak various speed definitions.

  Figure 10 (a)

  Figure 10 (a) Distribution of local relative errors on NOx for a congested period associated with a degraded speed definition: Speed limit (left); Punctual mean speed (right) -in terms of % of road segments. (b) Distribution of local relative errors on NOx for a congested period associated to a degraded definition of speed: Speed limit (left); Punctual mean speed (right) -in terms of % of the total emissions

Figure 12

 12 Figure12COPERT network FC (left) and NOx emissions (right) evaluated using three adapted emission functions for the "increased demand" scenario.

Figure 13

 13 Figure 13 Contribution to NOx emission by speed classes in terms of instantaneous speed or 6 min mean speed for a 6 mincongested period.

Figure 14

 14 Figure 14 Dynamic network PHEM emissions resulting from 15 replications.

Figure 14

 14 Figure 14 represents the time evolution of the fuel consumption (and NOx emissions) over the network obtained with PHEM. With no information on vehicle technology except for vehicle classes, the model stochastically associates a technology with a particular vehicle of the class considered (Passenger cars, Light Duty Vehicles, Buses or Heavy Duty Vehicle), conforming to the input fleet composition described in section 2.2. The boxplot shown in the figure is associated with PHEM emission values over 15 replications. The first observation is that the stochastic fleet definition not only impacts local emissions but also network emissions. The global relative gaps over time periods range from -5.0% to 3.7% for fuel consumption, and from -8.0% to 9.2% for NOx emissions.

Figure 15

 15 Figure 15 Comparison of network FC (left) and NOx emissions (right) with three modeling chains for the reference scenario

Figure 16

 16 Figure 16 Comparison of normalized network FC (left) and NOx emissions (right) with three modeling chains for the reference scenario.

Figure 17

 17 Figure 17 Comparison of network FC (left) and NOx emissions (right) with three modeling chains for the increased demand scenario.

Figure 18

 18 Figure 18 Comparison of normalized network FC (left) and NOx emissions (right) with three modeling chains for the increased demand scenario.

Table 1

 1 Performance indicators (%) reflecting the integration of traffic dynamics.

	NOx	PHEM	COPERT -speed	COPERT -mean
			classes	speed
	δ1	13.8	10.8	6.7
	δ2	21.0	15.5	6.3
	Δ	13.9	9.9	5.6

http://www.licit-lyon.eu/themes/realisations/plateformes/symuvia/

Acknowledgments

This work was performed in the framework of the Trafipollu research project funded by the National Research Agency under contract no. ANR-12-VBDU-0002-01.


L'effet photoacoustique (PA) permet de générer des ondes acoustiques à partir d'un objet éclairé par une onde lumineuse pulsée. L'imagerie photoacoustique, basée sur cet effet, présente ainsi les avantages des deux modalités dont elle est issue, notamment le contraste optique et la résolution des ultrasons. Cette nouvelle modalité, très prometteuse, s'applique parfaitement à l'imagerie de tissus biologiques et en particulier du système vasculaire. Dans ce contexte, nous avons réalisé un système expérimental d'imagerie PA dont nous présentons les caractéristiques ainsi que les résultats préliminaires obtenus sur des fantômes bi-modalités. Nous présentons les méthodes d'acquisition et de formation des images PA acquises sur ce système expérimental à l'aide d'une sonde cMUT (capacitive micromachined ultrasonic transducer). Nous montrons que la bande passante théoriquement infinie en réception est particulièrement adaptée aux signaux PA large bande. Plusieurs fantômes bimodalités, dont nous maitrisons les propriétés optiques et acoustiques ont été imagés. Nous décrivons la sensibilité du système en fonction de la nature des absorbeurs optiques choisis et de l'excitation laser multispectrale.

Introduction

L'effet photoacoustique (PA) est la génération d'ondes acoustiques par un objet éclairé par une onde lumineuse pulsée. Lorsque cet objet est soumis à des impulsions lumineuses, les absorbeurs optiques qu'il contient absorbent l'énergie apportée par chaque impulsion. Cela entraine une augmentation locale et instantanée de la température, donnant lieu à une expansion thermique, à l'origine d'une onde de pression acoustique [START_REF] Bell | Upon the production and reproduction of sound by light[END_REF][START_REF] Li | Photoacoustic tomography and sensing in biomedicine[END_REF].

Cet effet trouve des applications dans le domaine médical par l'illumination de tissus biologiques à l'aide d'un laser pulsé. Les ondes acoustiques générées au sein des tissus se propagent ensuite jusqu'à la surface où elles peuvent être détectées [START_REF] Li | Photoacoustic tomography and sensing in biomedicine[END_REF]. De ce fait, l'imagerie photoacoustique présente les avantages des deux modalités dont elle est issue, en s'affranchissant de leurs principales limites. Notamment, le contraste est celui de l'optique et la résolution celle des ultrasons. Cette modalité d'imagerie s'applique en particulier au système vasculaire qui possède un bon contraste d'absorption optique par rapport aux tissus environnants dans le proche infrarouge. De plus, l'utilisation d'une illumination multispectrale et la connaissance a priori des spectres d'absorption des absorbeurs principaux permet d'améliorer la sensibilité du système en ciblant l'absorbeur optique visé. En effet, un absorbeur sera particulièrement visible en PA aux longueurs d'ondes proches de son maximum d'absorption, et invisible à celles où il n'absorbe pas.

Le signal PA théorique d'un absorbeur sphérique a la forme d'un "N", présentant un signal triangulaire positif puis son équivalent négatif. Il s'agit du résultat de deux ondes : une onde de compression divergente partant de l'absorbeur et une onde de compression convergente vers le centre de ce dernier, devenant une onde de raréfaction qui arrive au niveau du capteur avec un certain retard [START_REF] Diebold | Photoacoustic monopole radiation in one, two, three dimensions[END_REF]. Le contenu fréquentiel d'un signal PA augmente lorsque la taille des absorbeurs est réduite. Par conséquent, on considère qu'un signal PA émit par des tissus biologiques est très large bande, de 1 MHz jusqu'à 100 MHz [START_REF] Li | Photoacoustic tomography and sensing in biomedicine[END_REF][START_REF] Diebold | Photoacoustic monopole radiation in one, two, three dimensions[END_REF].

L'acquisition de ces signaux peut être réalisée au moyen d'une sonde ultrasonore (US) semblable à celles utilisées en échographie. Toutefois, les sondes US classiques ont une bande passante limitée tant en transmission qu'en réception. La technologie cMUT (capacitive micromachined ultrasound transducer, [START_REF] Ergun | Capacitive micromachined ultrasonic transducers : fabrication technology[END_REF]) apparue il y a quelques années peut permettre de lever cette limitation car la bande passante en réception est théoriquement infinie [START_REF] Jin | Characterization of One-Dimensional Capacitive Micromachined Ultrasonic Immersion Transducer Array[END_REF].

Dans ce papier, nous présentons les caractéristiques de notre système expérimental et les résultats préliminaires in vitro obtenus sur des fantômes bi-modalités réalisés au laboratoire. Les acquisitions sont effectuées au moyen d'une sonde piezoélectrique et d'une sonde cMUT. Nous comparons ainsi les résultats obtenus par les deux sondes afin de mettre en avant l'intérêt de la technologie cMUT pour l'acquisition de signaux PA.

Sondes ultrasonores

Nous utilisons ici deux sondes linéaires ultrasonores dotées de technologies différentes. La première est une sonde linéaire classique, de type piezoélectrique (PZT) (L14-5W/60, Prosonic), la seconde est une sonde linéaire cMUT (Vermon). Les principales caractéristiques des deux sondes sont présentées dans le Tableau 1 (données constructeurs).

Fantômes

Plusieurs fantômes bi-modalités réalisés au laboratoire ont été imagés. Il s'agit de blocs de cryogel, un matériau constitué de 10% de PVA (Alcool Polyvinylique) et d'eau, dans lesquels des absorbeurs optiques ont été placés. Les propriétés optiques et acoustiques de ce milieu peuvent être ajustées selon le nombre de cycles de congélation qu'il subit [START_REF] Duboeuf | Investigation of PVA cryogel Young's modulus stability with time, controlled by a simple reliable technique[END_REF]. Avec 4 cycles, ce milieu présente un coefficient de diffusion réduit d'environ µ s = 6 cm -1 et un coefficient d'absorption négligeable (µ a << 1 cm -1 ). La vitesse de propagation du son est d'environ 1520 m.s -1 .

Les différents absorbeurs optiques utilisés, illustrés dans la Figure 2 3 Méthode

Acquisitions

Pour chacun des deux premiers fantômes, des acquisitions ont été réalisées dans la même configuration et avec les deux sondes ultrasonores. Ces acquisitions ont été effectuées à 1064 nm à différents niveaux d'énergie, jusqu'à 160 mJ/impulsion. Une image de fond, réalisée sans illumination laser, a également été acquise afin d'évaluer le bruit du système. Chaque acquisition dure quelques secondes et comporte une quarantaine d'images à partir desquelles un moyennage est réalisé.

Le dernier fantôme, présentant deux absorbeurs différents, a été imagé avec une illumination à 532 nm et 1064 nm, à 40 mJ/impulsion, à l'aide de la sonde cMUT.

Pour chacun des fantômes, une image échographique classique (mode B) a également été acquise, dans la même disposition, afin de superposer l'image PA à l'image mode B.

Reconstruction

Les images acquises sur les fantômes ont été reconstruites grâce à la méthode k-Wave [START_REF] Treeby | k-Wave : MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields[END_REF]. Les échos initialement rec ¸us par la sonde linéaire sont exprimés dans le domaine de Fourier. Les différents échantillonnages utilisés sont définis par :

                   k t ∈ [0; π/dt] k x ∈ [0; cπ/dx] k = k 2 t + k 2 x ( 1
)
avec dt l'échantillonnage temporel, dx l'échantillonnage spatial latéral. La focalisation des échos rec ¸us est donc réalisée par ré-échantillonnage du spectre. Ainsi, pour chaque ligne spatiale du spectre, une interpolation est réalisée de l'axe k t sur l'axe k. De cette manière, les divers échos sont focalisés et il a été montré que la reconstruction du signal PA est de meilleure qualité qu'avec des reconstructions basées sur des approches de retard temporel [START_REF] Treeby | k-Wave : MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields[END_REF]. Un autre avantage de cette reconstruction est sa rapidité de calcul et son implémentation temps réel.

Quantification

Deux critères sont évalués afin de comparer les acquisitions réalisées au moyen des sondes PZT et cMUT. Il s'agit du rapport signal-sur-bruit (SNR) et du rapport contraste-sur-bruit (CNR), permettant de juger de la qualité des images obtenues. Le calcul de ces deux paramètres nécessite de définir une région de l'image contenant du signal (S) et une région similaire caractérisant le bruit (N). Ces régions, que l'on appelle régions d'intérêt ou ROI, sont constituées de rectangles de 6x6 mm. La Figure 3 montre un exemple de deux ROI définies sur une des acquisitions effectuées sur l'inclusion sphérique. Les deux rapports sont calculés selon les expressions suivantes :

S NR = 20 log 10       S σ N       (2)
CNR = 20 log 10       S -N σ N       (3)
avec (.) la moyenne et σ(.) l'écart type de la ROI considérée.

Résultats

Discussion

Cette étude expérimentale montre la possibilité et l'intérêt d'utiliser une sonde cMUT pour faire de l'imagerie PA. Nous avons vérifié que celle-ci présente effectivement une bande passante en réception plus grande qu'une sonde PZT classique. Les faibles niveaux d'amplitude et le bruit électronique du système. La qualité des images obtenues se trouve ainsi améliorée, par rapport à la sonde PZT, comme le montrent les courbes du SNR, CNR et du maximum du signal, tous en faveur de la sonde cMUT.

Nous montrons également l'importance du choix de la longueur d'onde d'excitation laser, qui doit être sélectionnée en fonction de l'absorbeur optique visé, pour plus de spécificité. Grâce à la différence d'absorption des deux absorbeurs utilisés et une excitation multispectrale adaptée, nous pouvons différencier les deux encres.

Conclusion

Cette étude permet de mettre en avant l'intérêt de la technologie cMUT pour l'imagerie PA. Si la comparaison entre les sondes PZT et cMUT permet de se positionner en faveur de cette dernière, il faut poursuivre la caractérisation de cette sonde pour en connaître plus finement les performances.

Nous avons également montré que l'utilisation d'une excitation laser multispectrale permet de différencier et de localiser différents absorbeurs et pourrait ainsi permettre d'obtenir des informations fonctionnelles à partir des tissus biologiques. Une telle application soulève une autre question relevant du problème inverse optique, à savoir l'identification un absorbeur optique à partir du signal PA rec ¸u.

Une étape supplémentaire consisterait à augmenter le nombre de longueurs d'onde d'excitation laser, ce qui, pour des raisons matérielles, se ferait au détriment de l'énergie d'excitation laser.

  expérimentale La plateforme expérimentale d'imagerie PA se présente en deux parties : l'émission des impulsions laser et la réception des ondes US générées. La partie optique est constituée d'un laser pulsé Nd :YAG muni d'un tripleur de fréquence (Quanta-ray INDI, Spectra-Physics) délivrant des impulsions lumineuses de l'ordre de 5 ns à une fréquence de répétition de 10 Hz. Des acquisitions ont été réalisées à 532 nm et 1064 nm jusqu'à 160 mJ/impulsion. Le diamètre du faisceau incident est d'environ 1 cm. L'acquisition du signal PA est réalisée à l'aide d'un échographe clinique de recherche (SonixRP, Ultrasonix) couplé à un module d'acquisition (SonixDAQ, Ultrasonix). Ce module complémentaire permet l'acquisition en parallèle de 128 lignes radiofréquences (RF). Il est synchronisé sur le système laser, permettant ainsi, avant la formation de voie, d'acquérir une image au moment où débute l'impulsion laser. Cependant, un décalage de 1 µs est présent entre le début du signal de synchronisation et le début de l'acquisition des signaux RF.

Figure 1 -

 1 Figure 1 -Schéma du montage expérimental.

Tableau 1 -

 1 Caractéristiques des sondes utilisées. L14-5W/60 cMUT Fréquence centrale (MHz)

Figure 2 -

 2 Figure 2 -Photographies des fantômes bi-modalités imagés. (a) Gaine en plastique noir moulée dans du cryogel, (b) inclusion sphérique de cryogel teinté à l'encre de Chine, le tout dans du cryogel, (c) tubes remplis d'encre noire ou vermillon, moulés dans du cryogel.

Figure 3 -

 3 Figure 3 -Exemple de régions d'intérêt (ROI).

4. 1

 1 Comparaison des sondes cMUT et PZT La fréquence centrale et la bande passante des deux sondes sont calculées dans le cadre de cette application, à partir des acquisitions réalisées sur le fantôme contenant la gaine. Pour la sonde PZT nous obtenons une fréquence centrale autour de 4,5 MHz et une bande passante de 80% tandis que pour la cMUT la fréquence centrale est à 3,5 MHz et la bande passante est à 140%. La Figure 4 présente les images PA réalisées sur le fantôme contenant l'inclusion sphérique teintée, avec les deux sondes US. Les images PA sont représentées, ainsi que la superposition de ces signaux sur les images échographiques classiques. Elles correspondent aux acquisitions réalisées au maximum d'énergie, à savoir 160 mJ/impulsion. Pour ce même fantôme, nous avons calculé, à chaque niveau d'énergie et pour les deux sondes, le SNR, le CNR, ainsi que le maximum d'amplitude du signal. Les résultats sont présentés sur la Figure 5.

Figure 5 -

 5 Figure 5 -Evolutions du SNR, du CNR et du maximum d'amplitude du signal en fonction des niveaux d'énergie pour les deux sondes.

4. 2

 2 Illumination multispectrale Afin de tester la sonde cMUT avec une illumination multispectrale, nous avons acquis des images à deux longueurs d'onde différentes de notre fantôme contenant deux encres. La Figure 6 présente ces images PA superposées aux images US. Il est alors possible de comparer les amplitudes maximales du signal PA obtenu pour les deux encres et aux deux longueurs d'onde d'excitation. Ces valeurs sont données dans le Tableau 3.

Figure 4 -

 4 Figure 4 -Images PA et PA/US de l'inclusion sphérique teintée. (a) Image PA obtenue par la sonde PZT, (b) image PA superposée à l'image mode B, sonde PZT, (c) image PA obtenue par la sonde cMUT, (d) image PA superposée à l'image mode B, sonde cMUT.

Figure 6 -

 6 Figure 6 -Images PA et PA/US des tubes remplis d'encre noire (haut) et vermillon (bas). (a) et (b) Excitation à 1064 nm, (c) et (d) excitation à 532 nm.

  Tableau 3 -Amplitudes des signaux PA (u.a.).

		Encre noire Encre vermillon
	532 nm	100	54,8
	1064 nm	59,1	12,2

Remerciements

Ce travail a été réalisé grâce au soutien financier des LABEX CELYA (ANR-10-LABX-0060) et PRIMES (ANR-11-labx-0063) de l'Université de Lyon, dans le cadre du programme "Investissements d'Avenir" (ANR-11-IDEX-0007) géré par l'Agence Nationale de la Recherche (ANR). La sonde cMUT a été financée dans le cadre de l'ANR-11 TecSan-008-01 BBMUT de l'Université de Lyon.