David Coeurjolly
email: david.coeurjolly@liris.cnrs.fr

Valentin Blot
email: valentin.blot@ens-lyon.fr

Marie-Andrée Jacob-Da Col

Quasi-Affine Transformation in 3-D: Theory and Algorithms

Keywords: quasi-affine transform, periodic tiling, arithmetic, image transformation

In many applications and in many fields, algorithms can considerably be speed up if the underlying arithmetical computations are considered carefully. In this article, we present a theoretical analysis of affine transformations in dimension 3. More precisely, we investigate the arithmetical paving induced by the transformation to design fast algorithms.

Introduction

In many computer vision and image processing applications, we are facing new constraints due to the image sizes both in dimension with 3-D and 3-D+t medical acquisition devices, and in resolution with VHR (Very High Resolution) satellite images. This article deals with high performance image transformations using quasi-affine transforms (QATs for short), which can be viewed as a discrete version of general affine transformations. QAT can approximate rotations and scalings, and in some specific cases, QAT may also be one-to-one and onto mappings from Z n to Z n , leading to exact computations. In dimension 2, the QAT appeared in several articles [START_REF] Nehlig | Affine Texture Mapping and Antialiasing Using Integer Arithmetic[END_REF][START_REF] Jacob | Transformation of digital images by discrete affine applications[END_REF][START_REF] Nehlig | Applications quasi affines: pavages par images réciproques[END_REF][START_REF] Jacob | Applications quasi-affines[END_REF][START_REF] Col | Applications quasi-affines et pavages du plan discret[END_REF]. In higher dimension, theoretical results have been demonstrated [START_REF] Col | Quasi-linear transformations and discrete tilings[END_REF][START_REF] Blot | Quasi-affine transform in higher dimension[END_REF][START_REF] Blot | Quasi-affine transform in higher dimension[END_REF]. More precisely, authors have demonstrated the arithmetical and periodic structures embedded in n-dimensional QAT, leading to generic transformation algorithms. To implement these generic algorithms, several elements have to be fixed when considering a specific dimension. In this paper, we detail the computation of the minimal periods in dimension 3 leading to efficient transformation algorithms. Due to the space limitation, the proofs are available in the technical report [START_REF] Blot | Quasi-affine transform in higher dimension[END_REF]). In Section 2, we first present some definitions and results on n-D QAT. Section 3 focuses on the parameter computation in 3-D. Finally, Section 4 evaluate all the algorithms compared to a classical backward-mapping technique [START_REF] Foley | Computer Graphics: Principles and Practice[END_REF].

In this section, we present definitions and results obtained for QAT in higher dimension. These results have been independently proved by [START_REF] Col | Quasi-linear transformations and discrete tilings[END_REF] and [START_REF] Blot | Quasi-affine transform in higher dimension[END_REF][START_REF] Blot | Quasi-affine transform in higher dimension[END_REF]. In the following, we consider notations introduced in [START_REF] Blot | Quasi-affine transform in higher dimension[END_REF][START_REF] Blot | Quasi-affine transform in higher dimension[END_REF]. Let n denote the dimension of the considered space, V i the i th coordinate of vector V , and M i,j the (i, j) th coefficient of matrix M . We use the notation gcd(a, b, . . .) for the greatest common divisor of an arbitrary number of arguments, and lcm(a, b, . . .) for their least common multiple. Let a b denote the quotient of the euclidean division of a by b, that is the integer q ∈ Z such that a = bq + r satisfying 0 ≤ r < |b| regardless of the sign of b 4 . We also consider the straightforward generalization of these operators to n-dimensional vectors (e.g. V b is a vector where each component is the quotient of the division by b).

Definition 1. A quasi-affine transformation is a triple (ω, M, V ) ∈ Z×M n (Z)×
Z n (we assume that det(M ) = 0). The associated application is :

Z n -→ Z n X -→ M X + V ω Definition 2.
The inverse of a QAT (ω, M, V ) is the QAT:

(det(M ), ω com(M ) t , -com(M ) t V ) , (1)
where M t denotes the transposed matrix and com(M ) the co-factor matrix of M (Remind that M com(M ) t = com(M ) t M = det(M )I n .).

The associated affine application of the inverse of a QAT is therefore the inverse of the affine application associated to the QAT. However, due to the nested floor function, the composition f • f -1 is not the identity function in the general case. Let us recall the well-known Bezout Identity:

∀(a, b) ∈ Z 2 , ∃(u, v) ∈ Z 2 /au + bv = gcd(a, b) .
In Section 3, we have to consider a generalized form of the Bezout identity in dimension 3:

Proposition 1 ([7]). ∀(a, b, c) ∈ Z 3 , ∃(u, v, w) ∈ Z 3 /au+bv +cw = gcd(a, b, c) .
We present now several results and definitions that have been presented for n-dimensional QAT. All these results are given in [START_REF] Blot | Quasi-affine transform in higher dimension[END_REF][START_REF] Blot | Quasi-affine transform in higher dimension[END_REF] but we present here the main theorems, which will be used in the rest of the paper. First, the key feature of n-D QAT is that it contains a periodic paving structure. Definition 3 (Tile). Let f be a QAT. For Y ∈ Z n , we denote:

P Y = {X ∈ Z n /f (X) = Y } , (2)
P Y is called order 1 tile of index Y of f .
∀X ∈ P Y , ∃X ∈ P Z / M X + V ω = M X + V ω . (3)
Definition 5. P Y and P Z are said geometrically equivalent if:

∃v ∈ Z n /P Y = T v P Z , (4)
where T v denotes the translation of vector v.

The following theorem exhibits a relation between the notions of arithmetically and geometrically equivalent tiles.

Theorem 1 ( [START_REF] Blot | Quasi-affine transform in higher dimension[END_REF][START_REF] Blot | Quasi-affine transform in higher dimension[END_REF]). If P Y ≡ P Z , then P Y and P Z are geometrically equivalent. Since P Y ≡ P Z , there exists X ∈ P Y and X ∈ P Z such that:

M X + V ω = M X + V ω .
Then v = X -X is the translation vector:

P Y = T v P Z .
In Figure 1 we show some tiles of the QAT defined by (6, 3 1 -1 3 , 0 0 ) (a point of Z 2 is represented by a unit square whose bottom-left corner corresponds to the represented point). For each point X in a tile we provide its corresponding remainder M X+V ω . Tiles P (2,1) and P (0,0) are arithmetically equivalent, therefore they are also geometrically equivalent. It should also be noted that tiles P (1,0) and P (1,1) are geometrically equivalent but they are not arithmetically equivalent .

Definition 6 (Periodicity notations, [START_REF] Blot | Quasi-affine transform in higher dimension[END_REF][START_REF] Blot | Quasi-affine transform in higher dimension[END_REF]). For all 0 ≤ i < n, We define the set A i as follows:

A i = {α ∈ N * /∃(β j ) 0≤j<i ∈ Z i ,
∀(y 0 , . . . , y n-1 ) ∈ Z n , P y0,...,yi+α,...,yn-1 ≡ P y0+β0,...,yi-1+βi-1,yi,...,yn-1 } Furthermore, let us consider α i = min(A i ). We define

{β i j } 0≤j<i ∈ Z i and U i ∈ Z n such that ∀(y 0 , . . . , y n-1 ) ∈ Z n , P y0,...,yi+αi,...,yn-1 = T Ui P y0+β i 0 ,...,yi-1+β i i-1 ,yi,...,yn-1 .
The key theorem for the QAT periodic structure can now be presented: Theorem 2 (Periodicity, [START_REF] Blot | Quasi-affine transform in higher dimension[END_REF][START_REF] Blot | Quasi-affine transform in higher dimension[END_REF]). The set of QAT tiles is n-periodic, in other words

∀0 ≤ i < n, A i = ∅
Let us suppose that quantities α i , β i j and U i are given. The following theorem allows us to obtain P Y as the translation of a canonical tile P Y 0 . Theorem 3 ( [START_REF] Blot | Quasi-affine transform in higher dimension[END_REF][START_REF] Blot | Quasi-affine transform in higher dimension[END_REF]). ∀(y 0 , . . . , y n-1 ) ∈ Z n , we have P y0,...,yn-1 = T W P y 0 0 ,...,y 0

n-1 with W = n-1 i=0 w i U i and ∀n > i ≥ 0,        w i = yi+ n-1 j=i+1 wj β j i αi y 0 i = yi+ n-1 j=i+1 wj β j i αi .
In [START_REF] Blot | Quasi-affine transform in higher dimension[END_REF][START_REF] Blot | Quasi-affine transform in higher dimension[END_REF], we have proved that canonical tiles P Y 0 are associated to grid points of a special tile called super-tile.

Definition 7 (Super-tile, [START_REF] Blot | Quasi-affine transform in higher dimension[END_REF][START_REF] Blot | Quasi-affine transform in higher dimension[END_REF]). A super-tile of a QAT is the set P such that P = 0≤Y 0 <(α0,...,αn-1)

P Y 0
Theorem 4 ( [START_REF] Blot | Quasi-affine transform in higher dimension[END_REF][START_REF] Blot | Quasi-affine transform in higher dimension[END_REF]). P is the tile P (0,...,0) of the QAT defined by: arithmetically distinct tiles (the tiles with same color are arithmetically equivalent). In this example, for all i, j ∈ N, P (i+5,j) ≡ P (i,j) and P (i+2,j-3) ≡ P (i,j) .

  ω lcm 0≤i<n (α i ),    θ 0 • • • 0 . . . . . . . . . 0 • • • θ n-1    M,    θ 0 • • • 0 . . . . . . . . . 0 • • • θ n-1    V    , with ∀0 ≤ i < n -1, θ i = lcm 0≤j<n-1 (αj ) αi .
(i+2,j-3) (i+5,j) P P P (i,j) Fig. 2. Periodicity of thes tiles of a 2D QAT .

The set P (i,j) |i = 0, 1, 2, 3, 4, j = 0, 1, 2 contains exactly once all distinct tiles so it is a super-tile of the QAT. In [START_REF] Blot | Quasi-affine transform in higher dimension[END_REF][START_REF] Blot | Quasi-affine transform in higher dimension[END_REF], we have demonstrated that if we consider the Hermite Normal Form of the QAT matrix M , then efficient tile construction can be designed. In the following, let M H = T denote the Hermite Normal Form (with det(H) = ±1 and T upper triangular). Note that this form always exists for nonsingular integer square matrix.

Theorem 5 (Tile Construction, [START_REF] Blot | Quasi-affine transform in higher dimension[END_REF][START_REF] Blot | Quasi-affine transform in higher dimension[END_REF]). ∀Y ∈ Z n , let M H = T be the Hermite Normal Form of the QAT matrix M ,

P Y = {HX / ∀n > i ≥ 0, A i (X i+1 , . . . , X n-1 ) ≤ X i < B i (X i+1 , . . . , X n-1 )} With A i (X i+1 , . . . , X n-1 ) = - -ωY i + n-1 j=i+1 T i,j X j + V i T i,i , B i (X i+1 , . . . , X n-1 ) = - -ω(Y i + 1) + n-1 j=i+1 T i,j X j + V i T i,i .
In Algorithm 1, we give the generic algorithm applying a contracting QAT f to an image A (see Fig. 3). The principle is that we give to each pixel Y of image B the average color of the tile P Y in image A. If f is a dilating QAT, we obtain a similar algorithm in which we replace f with f -1 , and then we give the color of each pixel Y of image A to each pixel of P Y in image B (see Fig. 3 for an illustration in 2-D). In both algorithms, some elements cannot be computed in arbitrary dimension n. Indeed, even if there exist algorithms to compute the Hermite Normal Form of an arbitrary square integer matrix [START_REF] Storjohann | Asymptotically fast computation of hermite normal forms of integer matrices[END_REF], there is no generic algorithm to obtain the minimal periodicities {α i } (see discussion in Sect. 5). In the next section, we focus on the minimal periodicity computation in dimension 3.

QAT in Dimension 3

In dimension 3, we consider the following framework: we first define the Hermite Normal Form, the minimal periods and then we detail the transformation algorithm.

Hermite Normal Form and Tile Construction

Let us consider a QAT (ω, M, V

) with M =   a 0 b 0 c 0 d 0 e 0 f 0 g 0 h 0 i 0   and V =   j 0 k 0 l 0   .
In [START_REF] Blot | Quasi-affine transform in higher dimension[END_REF], we present explicit formulas to compute the Hermite Normal Form in 3-D. In the following, we define

H = H 1 H 2 H 3 H 4 and M H = T =   a b c 0 d e 0 0 f   .
Thanks to Hermite decomposition, we have a > 0, d > 0 and f > 0. To construct the tile of index (i, j, k) and thanks to Theorem 5, we have: Let α h = a h and U = HY . Then α h > 0, P i+α h ,j,k ≡ P i,j,k and ∀(i, j, k) ∈ Z 3 , P i+α h ,j,k = T U P i,j,k .

A 2 = - -ωk + l 0 f , B 2 = - -ω(k + 1) + l 0 f A 1 (z) = - -ωj + k 0 + ez d , B 1 (z) = - -ω(j + 1) + k 0 + ez d A 0 (y, z) = - -ωi + j 0 + by + cz a , B 0 (y, z) = - -ω(i + 1) + j 0 + by + cz a Algorithm 2: Tile construction in 3-D A2 ← - -ωk+l 0 f ; B2 ← - -ω(k+1)+l 0 f ; for z ← A2 to B2 -1 do A1 ← - -ωj+k 0 +ez d ; B1 ← - -ω(j+1)+k 0 +ez d ; for y ← A1 to B1 -1 do A0 ← - -ωi+j 0 +by+cz a ; B0 ← - -ω(i+1)+j 0 +by+cz a ; for x ← A0 to B0 -1 do H   x y z   ∈ P i,j,k ; 3
Proof. The proof is detailed in [START_REF] Blot | Quasi-affine transform in higher dimension[END_REF].

Theorem 7. The period α h is a minimal horizontal period, i.e. α h = α 0 .

Proof. The proof is given in [START_REF] Blot | Quasi-affine transform in higher dimension[END_REF].

Concerning the vertical period, let:

d v = d gcd(d, ω) , ω v = ω gcd(d, ω) , a v = a gcd(a, bω v , ω) , φ = bω v gcd(a, bω v , ω) , ω v = ω gcd(a, bω v , ω) , α v = gcd(a v , ω v ), u 1 and v 1 are such that : a v u 1 + ω v v 1 = gcd(a v , ω v )(= α v ), β 0 = -φv 1 , Y =   -φu 1 ω v α v 0   Theorem 8 (Vertical Periodicity). Let α v = d v α v , U = HY .
Then α v > 0, P i,j+αv,k ≡ P i+β0,j,k and ∀(i, j, k) ∈ Z 3 , P i,j+αv,k = T U P i+β0,j,k .

Proof. The proof is given in [START_REF] Blot | Quasi-affine transform in higher dimension[END_REF].

Theorem 9. The period α v is a minimal vertical period, i.e.

α v = α 1 .
Proof. The proof is detailled in [START_REF] Blot | Quasi-affine transform in higher dimension[END_REF].

For the last period, let us consider:

f d = f gcd(ω, f ) , ω d = ω gcd(ω, f ) , d d = d gcd(d, eω d , ω) , φ = eω d gcd(d, eω d , ω) , ω d = ω gcd(d, eω d , ω)
, u 1 and v 1 are such that :

d d u 1 + ω d v 1 = gcd(d d , ω d ), ψ = cω d gcd(d d , ω d ) -bφu 1 , a d = a gcd(a, ψ, ω, ω d b gcd(d d ,ω d ) )
, ψ = ψ gcd(a, ψ, ω,

ω d b gcd(d d ,ω d ) )
, Proof. The proof is detailed in [START_REF] Blot | Quasi-affine transform in higher dimension[END_REF].

ω d = ω gcd(a, ψ, ω, ω d b gcd(d d ,ω d ) ) , χ = ω d b gcd(d d ,ω d ) gcd(a, ψ, ω, ω d b gcd(d d ,ω d ) ) , α d = gcd(a d , χ, ω d ), α d = α d gcd(d d , ω d ), u 2 , v 2 and w 2 are such that : a d u 2 + χv 2 + ω d w 2 = gcd(a d , χ, ω d )(= α d ), k = -ψ v 2 , β 0 = -ψ w 2 , β 1 = -φv 1 α d -k d d gcd(d d , ω d ) , Y =    -ψ u 2 -φu 1 α d + k
Theorem 11. The period α d is a minimal depth period, i.e. α d = α 2 .

Proof. The proof is detailed in [START_REF] Blot | Quasi-affine transform in higher dimension[END_REF].

Based on these periods, we can construct the super-tile and all the initial period tiles. To design a transformation algorithm, for each point X ∈ P, we need to determine the tile index Y to which X belongs to. Since X ∈ P Y ⇔ M X+V ω = Y , Algorithm 3 details the initial period tile construction with scanning points in P. The computational cost of Alg. 3 exactly corresponds to the number of tiles in the initial period.

Proposition 2. The number of tiles of the initial period is ω d ω d ω d .

In the Proposition statement, we do not give the closed formula. However, ω d ω d ω d is equal to ω 3 divided by a product of three gcd().

Algorithm 3: Super-tile and initial period tiles construction in 3-D.

A 2 ← - l 0 f ; for z ← A 2 to A 2 + ωα 2 f -1 do A 1 ← - k 0 +ez d ; for y ← A 1 to A 1 + ωα 1 d -1 do A 0 ← - j 0 +by+cz a ; for x ← A 0 to A 0 + ωα 0 a -1 do Y ←        T     x y z     +V ω        ; H   x y z   ∈ P Y ;
Proof. The proof is detailed in [START_REF] Blot | Quasi-affine transform in higher dimension[END_REF].

Using Theorems 4 and 5, we have

θ 0 = lcm(α 0 , α 1 , α 2 ) α 0 , θ 1 = lcm(α 0 , α 1 , α 2 ) α 1 , θ 1 = lcm(α 0 , α 1 , α 2 ) α 2 ,
and

P =    H   x y z   /A 2 ≤ z < B 2 , A 1 (z) ≤ y < B 1 (z) and A 0 (y, z) ≤ x < B 0 (y, z)    , with A 2 = -l0 f , A 1 (z) = -k0+ez d , A 0 (y, z) = -j0+by+cz a , B 2 = A 2 + ωα2 f , B 1 (z) = A 1 (z) + ωα1
d , and B 0 (y, z) = A 0 (y, z) + ωα0 a (see [START_REF] Blot | Quasi-affine transform in higher dimension[END_REF] for details).

Figure 4 shows the 16 distinct tiles of the QAT (11,

  6 -2 3 2 -2 -3 4 4 -4   ,   0 0 0   ) in Z 3
and its correponding super-tile.

QAT Algorithm in 3-D

To obtain the overall QAT algorithm, we need to find both the initial period tile index and the translation vector associated to a given tile P i,j,k . Hence, thanks to Theorem 3, we have

∀(i, j, k) ∈ Z 3 , P i,j,k = T W P i 0 ,j 0 ,k 0 with W = w 0 U 0 + w 1 U 1 + w 2 U 2 and 0 ≤ k 0 = k α 2 < α 2 , w 2 = k α 2 , 0 ≤ j 0 = j + w 2 β 2 1 α 1 < α 1 , w 1 = j + w 2 β 2 1 α 1 , 0 ≤ i 0 = i + w 1 β 1 0 + w 2 β 2 0 α 0 < α 0 , w 0 = i + w 1 β 1 0 + w 2 β 2 0 α 0 .
Fig. 4. All distinct tiles of a QAT in Z 3 and its super-tile.

Experiments

The algorithms were implemented in both 2D and 3D, with different refinements in order to be able to compare the implementations. For details on the 2-D algorithms, refer to [START_REF] Blot | Quasi-affine transform in higher dimension[END_REF]. The backward mapping (B. M. for short) implementation let Algorithm 4: QAT Algorithm in 3-D.

Input: A QAT (ω, M, V ) and an image g

: [0, t0] × [0, t1] × [0, t2] → Z Output: an image h : [min0, max0] × [min1, max1] × [min2, max2] → Z Compute mini and maxi quantities from ti; if f dilating then f ← f -1 ;
Compute the Hermite Normal Form of the matrix M ; Compute the minimal periodicities {α0, α1, α2} and vectors {U 0, U 1, U 2}; Use Algorithm 3 to compute the canonical tiles in the super-tile P;

for i ← 0 to t0 -1 do for j ← 0 to t1 -1 do for k ← 0 to t2 -1 do Compute W, i 0 , j 0 , k 0 ; h(T W P i 0 ,j 0 ,k 0 ) ← g(i, j, k) ;
else Compute the Hermite Normal Form of the matrix M ; Compute the minimal periodicities {α0, α1, α2} and vectors {U 0, U 1, U 2}; Use Algorithm 3 to compute the canonical tiles in the super-tile P; for i ← min0 to max0 -1 do for j ← min1 to max1 -1 do for k ← min2 to max2 -1 do Compute W, i 0 , j 0 , k 0 ; h(i, j, k) ← g(T W P i 0 ,j 0 ,k 0 ) ; us compare the tile periodicity method with the widely used backward mapping method [START_REF] Foley | Computer Graphics: Principles and Practice[END_REF]. The simple implementation does not use tiles periodicity and uses algorithm 2 for each tile. The periodicity implementation uses the periodicity and the algorithm 4. The noMultiply implementation additionally uses a method presented in [START_REF] Jacob | Applications quasi-affines[END_REF] Figure 5 illustrates the results in dimension 2. As expected, when comparing B.M. and Periodicity, results are similar for both contracting and isometry QATs. Differences appear when dilating QAT is considered. Indeed, since a unique color is associated to a tile in the Periodicity algorithm, the transformed image contains sharp edges (Fig 5-(l)) On the other hand, the interpolation process in the B.M. algorithm makes the image blurred. To compare the time efficiency (Table 1), we have considered two quantities: the total number of elementary operations of the main loop5 and the overall computational time in seconds. Table 1 and Figure 6 present the results in dimension 3. For the sake of clarity, we have only considered an input binary image but the transformation algorithms can be applied to 3-D color images. We have also performed a peak signal-to-noise ratio (PSNR for short, given in decibel dB) computation between the input image in 2-D and 3-D, and the result of the composition f.f -1 . This test has been designed to evaluate the propagation of the error through the transformations with a signal processing tool. As presented in Table 1, the distortion induced by the proposed method is always smaller than the one induced by the backward mapping technique. Note that to have relevant measurements, we have used a density 3-D volume for the 3-D test (see Fig. 6-(h)).

(a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (l) (m)

Conclusion and Future Works

In this paper, we have first re-demonstrated an existing result in dimension 2 with our formalism and provide a generalization in dimension 3 of fast image transformations with QAT. Compared to a classical image transformation technique, we have also illustrated the strength of such arithmetical methods to speed up transformations in higher dimensions. As detailed above and based on the generic n -D transformation algorithm proposed in [START_REF] Col | Quasi-linear transformations and discrete tilings[END_REF][START_REF] Blot | Quasi-affine transform in higher dimension[END_REF][START_REF] Blot | Quasi-affine transform in higher dimension[END_REF], we had to perform specific computation to obtain the minimal periods in dimension 2 and 3. A very challenging future work is to define a framework to compute these minimal periods when a specific dimension is considered. In order to achieve this goal and instead of having explicit formulas, we plan to investigate an algorithmic solution based on the n -D canonical tile counting algorithm proposed in [START_REF] Col | Quasi-linear transformations and discrete tilings[END_REF]. Furthermore, in dimension 2 and 3, comparisons to other transformation algorithms, such as Fourier based technique, are also of interest.

Fig. 1 .

 1 Fig. 1. Example of tiles and remainders

Figure 2

 2 Figure 2 illustrates tiles of the QAT (84, 12 -11 18 36 , 0 0 ) in Z 2 with 15

Algorithm 1 :Fig. 3 .

 13 Fig. 3. Illustration in dimension 2 of the QAT algorithm when f is contracting (a) and dilating (b). In both cases, we use the canonical tiles contained in the super-tile to speed-up the transformation.

. 2 MinimalTheorem 6 (

 26 Periodicity and Super-tile Construction In dimension 3, we need to compute the periodicity along each dimension. Let us first denote a h = a gcd(a,ω) , ω h = ω gcd(a,ω) Horizontal Periodicity).

Theorem 10 (

 10 Depth Periodicity). Let α d = α d f d , U = HY . Then α d > 0 P i,j,k+α d ≡ P i+β0,j+β1,k and ∀(i, j, k) ∈ Z 3 , P i,j,k+α d = T U P i+β0,j+β1,k

  which uses a handling of remains instead of computing a matrix product in 3. The experiments are performed on an Intel c Centrino c Duo T2080 (2 x 1.73 GHz) in monothread and we give on one hand the time of computation and on the other hand the number of elementary instructions. The QATs used are the following : In 2-D: dilating application The pictures are of size : 200 x 171 in 2-D and 10 x 10 x 10 in 3-D (simple cube).

Fig. 5 .

 5 Fig. 5. Results in dimension 2: (a -d) Contracting QAT (B.M. (a -b) and Periodicity (c -d)); (e -h) Isometry (B.M. and Periodicity); (i -l) Dilating (B.M. and Periodicity). (m) illustrates the tile structure of the dilating QAT on the same square as in (l).

Fig. 6 .

 6 Fig. 6. Results in dimension 3: (a -b) Contracting (B.M. and Periodicity), (c-d) Isometry and (e-f ) Dilating (B.M. and Periodicity). (g) illustrates the tile structure of the dilating QAT. (h) is the 32 × 32 × 32 density volume used in the f.f -1 composition test (from the TC18 data-set http://www.tc18.org).

Table 1 .

 1 Comparative evaluation in 2-D and 3-D. The last table presents the PSNR evaluation of the composition f • f -1 .

			2-D -instructions (time in sec.)
		B.M.		simple	Periodicity	noMultiply
	contracting 1 607 774 (0.036) 64 536 315 (0.06) 29 578 702 (0.036) 27 679 044 (0.036)
	isometry 63 058 160 (0.112) 57 619 374 (0.064) 39 682 795 (0.056) 35 875 892 (0.044)
	dilating 391 622 017 (0.404) 185 956 768 (0.12) 87 490 567 (0.084) 83 472 387 (0.078)
			3-D -instructions (time in sec.)
		B.M.		simple	Periodicity
	contracting	15 864 982 (0.02)	47 303 861 (0.052)	12 865 125 (0.012)
	isometry	750 102 224 (0.416)	51 121 827 (0.068)	15 234 007 (0.016)
	dilating 170 072 035 547 (79.637) 2 479 676 409 (1.384) 7 760 893 011 (0.632)
			2-D (PSNR in dB) 3-D (PSNR in dB)
			B.M. Periodicity B.M. Periodicity
		contracting 24.764	23.4823	17.8026	17.0304
		isometry	27.7619	25.8052	19.4115	15.9481
		dilating	31.2331	30.8375	20.4435	16.7862

obtained with the valgrind profiling tool.


Esteban Paduro
email: paduro@gmail.com

Mario Sigalotti

Approximate Controllability of the Two Trapped Ions System

We prove the approximate controllability of a bilinear Schrödinger equation modelling a two trapped ions system. A new spectral decoupling technique is introduced, which allows to analyze the controllability of the infinite-dimensional system through finitedimensional considerations.

Introduction

In this paper we study the controllability of a system modelling a two trapped ions system driven by two laser beams.

The Model

The two trapped ions system models a pair of identical charged particles confined to a small region of the space. Both ions are stabilized by the same spatial oscillation, given by a harmonic oscillator with of frequency ω. The controls are monochromatic lasers of frequencies Ω and Ω ± ω.

The state of the system is represented by a function φ = (φ gg , φ ge , φ eg , φ ee ) ∈ (L 2 (R)) 4 , which represents the different possible configurations of electrons in the ions (e for excited state, g for ground state). In the Lamb-Dicke limit, the system can be written (see [START_REF] Rouchon | Quantum systems and control[END_REF]) as the following two trapped ions system

i d dt φ gg = (u 1 + u 1r a † + u 1b a)φ eg + (u 2 + u 2r a † + u 2b a)φ ge , i d dt φ eg = (u * 1 + u * 1r a + u * 1b a † )φ gg + (u 2 + u 2r a † + u 2b a)φ ee , i d dt φ ge = (u 1 + u 1r a † + u 1b a)φ ee + (u * 2 + u * 2r a + u * 2b a † )φ gg , i d dt φ ee = (u * 1 + u * 1r a + u * 1b a † )φ ge + (u * 2 + u * 2r a + u * 2b a † )φ eg , (1.1)

Existing Literature

For a single trapped ion counterpart of (1.1), exact controllability has been proved on some finite-dimensional subspaces of (L 2 (R))

2 in [START_REF] Law | Arbitrary control of a quantum electromagnetic field[END_REF] and [START_REF] Kneer | Preparation of arbitrary entangled quantum states of a trapped ion[END_REF]. In [START_REF] Yuan | Controllability of the coupled spin-1 2 harmonic oscillator system[END_REF] an approximate controllability result on (L 2 (R))

2 is obtained using spectral properties of coupling operators, allowing the decoupling of the internal dynamics. A similar approach is proposed in a more general framework in [START_REF] Bliss | Quantum control of infinite-dimensional many-body systems[END_REF]. A different approach, based on adiabatic evolution, is proposed in [START_REF] Adami | Controllability of the Schrödinger Equation via Intersection of Eigenvalues[END_REF]. The system modelling the single trapped ion before Lamb-Dicke limit has been studied in [START_REF] Ervedoza | Approximate controllability for a system of Schrödinger equations modeling a single trapped ion[END_REF], obtaining approximate controllability in (L 2 (R)) 2 and also in with respect to some stronger norms. Such results are based on a set of explicit estimates of the approximation error with respect to the Lamb-Dicke limit. The approximate controllability of a different model for the interaction of between a bosonic mode and a two-level system is proved in [START_REF] Vittorio Boscain | On the control of spin-boson systems[END_REF].

For several trapped ions in the Lamb-Dicke limit, in [START_REF] Anthony M Bloch | Finite Controllability of Infinite-Dimensional Quantum Systems[END_REF] (see also [START_REF] Rangan | Control of trapped-ion quantum states with optical pulses[END_REF]) an approximate controllability result is obtained, based on the analysis on a sequence of nested finite-dimensional system, which can be decoupled from the rest of system. In the recent paper [START_REF] Keyl | Controlling Several Atoms in a Cavity[END_REF] the (approximate) controllability of the system is established by considering different families of controlled dynamics. Many other works studying the two trapped ions system deal with the construction of quantum gates [START_REF] Cirac | Quantum Computations with Cold Trapped Ions[END_REF][START_REF] Childs | Universal quantum computation with two-level trapped ions[END_REF][START_REF] Jonathan | Fast quantum gates for cold trapped ions[END_REF] but up to our knowledge, they do not present general controllability results for system (1.1).

Main Results

Theorem 1.1 (Approximate Control of the two trapped ions system). Let ε > 0, M > 0 and (φ 0 , φ T ) ∈ (L 2 (R)) 4 , with φ 0 = φ T . Then, there exists T 0 > 0 such that for all T > T 0 we can find

u 1 , u 1r , u 1b , u 2 , u 2r , u 2b ∈ L ∞ ((0, T ), C) with L ∞ -norm smaller than M such that the solution φ(t, x) of system (1.1) with initial data φ(0) = φ 0 satisfies φ(T ) -φ T (L 2 (R)) 4 < ε.
For the proof we introduce a new finite-dimensional approximation of system (1.1) with good control properties. In a second step we study spectral properties of the coupling operators of the original system; this allows us to construct the control laws of the original system based on those of the finite-dimensional approximation. Finally, using estimates for the approximation error, we obtain the approximate controllability of system (1.1).

The proof also shows that not all vector fields are necessary for the approximate controllability of the system. In particular, it is sufficient to exploit red shift only (respectively, blue shift only), that is, the system is still approximately controllable if one sets to zero the controls u 1b , u 2b (respectively, the controls u 1r , u 2r ). See the appendix for more details on this issue.

Outline of the Work

The paper is organized as follows. In Section 2 we develop the technical tools required for the proof. Section 3 is devoted to the proof of Theorem 1.1. In Section 3.1 we present a modal representation for the system. Section 3.2 discusses the modal approximation of the problem and some auxiliary controllability properties. In Section 3.3 we present the decoupled modal approximation and the proof of the main theorem.

A Decoupling Technique for Oscillatory Systems

In this section we study a decoupling technique for oscillatory systems. This technique is based on the spectral decomposition of skew-hermitian operators and non-resonance conditions that make the decoupling possible. This is the main technical tool used to pass from the controllability result for the finite-dimensional approximation to the approximate controllability of the infinite-dimensional system.

Preliminary Definitions

Let X be an infinite-dimensional separable Hilbert space and {φ j } j∈N an orthonormal basis of X. Let U : X → X be a skew-hermitian operator with compact resolvent. Denote by Σ(U ) = {ω j } Γ j=1 ⊂ [0, +∞) the sequence of distinct moduli of the eigenvalues of U , with Γ ∈ N ∪ {∞}. We say that the elements of Σ(U ) are the frequencies associated with the operator U .

For each frequency w ∈ Σ(U ), we define the subspaces A w ⊂ X as the space generated by eigenfunctions associated with the eigenvalues with modulus w, that is,

A w = span{x | U x = λx, |λ| = w},
and for sets of frequencies B ⊂ Σ(U ) we set

A B = w∈B A w . Definition 2.1 (Resonant Class). For ν = 0 we define the set of Q-resonant frequencies with ν as R(ν) = w ∈ Σ(U ) | w ν ∈ Q \ {0} .
We find also useful to set R(0) = {0} and to introduce, for every m ∈ N and every ν ∈ R,

R m (ν) = R(ν) ∩ {ω j } j<m .
For every m ∈ N, we define an equivalence relation in {ω 1 , . . . , ω m-1 } by ω j ∼ ω k if and only if R m (ω j ) = R m (ω k ) and we denote by N (m) the number of equivalence classes of such a relation. For each class C j = {0}, j = 1, . . . , N (m), we choose ν j > 0 such that w ν j ∈ N for each w ∈ C j (ν j exists because C j is finite). We also set v j = 0 if C j = {0}. By definition of ∼, the elements of {ν j | j = 1, . . . , N (m), ν j = 0} are Q-linearly independent. Definition 2.2 (Decoupled Decomposition). For every m ∈ N with m ≤ Γ the decoupled decomposition of U at order m is the family of operators U 1 , . . . , U N (m) , U dec , U ρ : X → X given by     

U j = U Π[A R m (ν j ) ], j = 1, . . . , N (m), U dec = U Π[A ωm ], U ρ = U -N (m) j=1 U j -U dec = U I -N (m) j=1 Π[A R m (ν j ) ] -Π[A ωm ] , (2.1)
where Π[C] denotes the orthogonal projection over the space C. By definition

U = U 1 + • • • + U N (m) + U dec + U ρ . (2.2)
Some basic properties of this decomposition are summarized in the following lemma.

Lemma 2.3. Let U be a skew-hermitian operator with compact resolvent and let U 1 , . . . , U N (m) , U dec , U ρ be the decoupled decomposition at order m defined as above. Then

(a) -im (U j ) ⊂ A R m (ν j ) , j = 1, . . . , N (m), -im (U dec ) ⊂ A ωm , -im (U ρ ) ⊂ N (m) j=1 A R m (ν j ) ⊕ Aω m ⊥ . (b) Each pair of operators in the decomposition commute. Moreover, if U, V ∈ {U 1 , . . . , U N (m) , U dec , U ρ } and U = V then U V = V U = 0.
Proof. (a) It is enough to notice that the spaces A R m (ν i ) are direct sum of eigenspaces of U . The space A R m (ν j ) is then invariant under U . (b) It follows from the fact that eigenspaces associated with different eigenvalues are orthogonal, and therefore the image of each operator is contained in the kernel of the other.

Approximate Decoupling of Skew-Hermitian Operators

In this section we develop the main technical tool used in this paper, namely, an approximation method for the operators exp (U j t) by means of exp(U τ ). This technique is based on the geometric fact classically known as the "irrational windings of the torus", which says that the integer multiples of an irrational number modulo 1 are dense in interval the [0, 1). In [START_REF] Childs | Universal quantum computation with two-level trapped ions[END_REF][START_REF] Gulde | Implementation of the Deutsch-Jozsa algorithm on an ion-trap quantum computer[END_REF][START_REF] Yuan | Controllability of the coupled spin-1 2 harmonic oscillator system[END_REF] similar decoupling techniques are used for control purposes. Our formulation gives an abstract framework which could be applied to other problems.

For the proof we need a n-dimensional formulation of the mentioned result. The proof of this result is classical (see, e.g., [START_REF] Katok | Introduction to the modern theory of dynamical systems[END_REF]Prop. 1.4.1]).

Lemma 2.4 (Irrational Winding of the Torus). Let T n = R n /Z n be the n-dimensional torus with the usual distance. Let ϕ the automorphism: x → x + ω (mod 1), where ω ∈ R n , x ∈ T n . Then the orbits of ϕ are everywhere dense if and only if

k • ω ∈ Z with k ∈ Z n ⇒ k = 0.
The main result of the section is the following. Theorem 2.5 (Approximate Decoupling of Skew-Hermitian Operators). Let X be a separable Hilbert space and U : X → X a skew-hermitian operator with compact resolvent. Let Σ(U ), R(•) and N (•) be defined as in Section 2.1. Suppose that there exists m ∈ N such that ω m = 0 and

ω j ωm ∈ (R \ Q) ∪ {0} for each j < m. Let Y be a subspace of X such that j<m A ω j ⊂ Y and j>m A ω j ⊂ Y ⊥ . (2.3)
Then, given t ∈ R, ℓ ∈ 1, . . . , N (m), and ε > 0 there exist t ∈ R and a unitary operator Σ : X → X such that exp( t U ) admits the decomposition

exp( t U ) = exp( t U ℓ )Σ = Σ exp( t U ℓ ), (2.4)
where U ℓ is defined by (2.1) and the operator Σ satisfies

(Σ -I)| Y L(Y,X) < ε.
Proof. Let ℓ ∈ {1, . . . , N (m)}, t ∈ R, and ε > 0 be fixed. We split the proof in 4 steps.

Step 1 Consider the decoupled decomposition (2.1) of the operator U , with m satisfying the hypotheses of the theorem. By Lemma 2.3 the terms in the decomposition commute and then we have

exp(tU ) = exp(tU 1 ) • • • exp(tU N (m) ) exp(tU dec ) exp(tU ρ ). (2.5)
By definition of ν ℓ , either

ν ℓ = 0 or w ν ℓ ∈ N with w ∈ R m (ν ℓ ). Let νℓ = ν ℓ if ν ℓ = 0 while, if ν ℓ = 0, choose νℓ = 1. Define t by t = t + 2π νℓ s, s ∈ N, (2.6)
where s will be chosen later. Then we have

w t = w t + 2π w νℓ s ≡ w t (mod 2π), ∀w ∈ R m (ν ℓ ). (Notice that w = 0 if ν ℓ = 0.) By definition of U ℓ , it follows that exp( tU ℓ ) = exp( tU ℓ ), (2.7)
independently of the choice of s ∈ N.

Step 2 By definition of the Q-resonant frequency classes, the elements of

{ν j } j∈{1,...,ℓ-1,ℓ+1,...,N (m)},ν j =0 ∪ {ω m } are Q-linearly independent. Let N (m) be equal to N (m) if ν j = 0 for every j ∈ {1, . . . , ℓ -1, ℓ + 1, . . . , N (m)} and to N (m) -1 otherwise. Define ν1 , . . . , ν N (m) in such a way that {ν j } j∈{1,...,ℓ-1,ℓ+1,...,N (m)},ν j =0 ∪ {ω m } = {ν j } N (m) j=1 .
Consider the application

F : R → T N (m) , s → s 2π νℓ ν1 , . . . , ν N (m) .
By Lemma 2.4 we obtain that F (N) is dense in the torus T N (m) . Since F (s) can be taken as close to -t ν1 , . . . , ν N (m) as desired, it is possible, for every c > 0 (to be chosen later suitably small) to select s ∈ N in such a way that

d T 1 s 2π νℓ νj , -tν j < cε, for j = 1, . . . , N (m).
By definition of t, the latter system of inequalities can be rewritten as

d T 1 ( t νj , 0) < cε, for j = 1, . . . , N (m).
Equivalently, there exist δ 1 , . . . , δ N (m) such that |δ j | < cε, δ j = t νj (mod 2π), for j = 1, . . . , N (m).

(2.8)

Step 3 Let h ∈ {1, . . . , ℓ -1, ℓ + 1, . . . , N (m)}. By construction, U h only has finitely many eigenvalues, so we can write its spectral decomposition as

tU h = iθ 1 t Π θ 1 + • • • + iθ n t Π θn , (2.9)
where iθ 1 , . . . , iθ n ∈ iR are the eigenvalues of U h | A R m (ν h ) and Π θ k the projection operator onto the eigenspace associated with iθ k . Using the representation (2.9) we can compute the exponential of the operator U h as follows

exp( tU h ) = e iθ 1 t Π θ 1 + • • • + e iθn t Π θn + Π[A ⊥ R m (ν h ) ].
Fix k = 1, . . . , n. Notice that either ν h = 0 = θ k or there exists j ∈ {1, . . . , N (m)} such that ν h = νj and |θ k | νj is in N. By equation (2.8) it then follows that

d T 1 (θ k t, 0) = d T 1 θ k νj νj t, 0 ≤ |θ k | νj cε.
Hence,

e tU h -I L(X) ≤ n k=1 e i t θ k -1 ≤ n k=1 d T 1 (θ k t, 0) ≤ cε n k=1 |θ k | νj ,
with e tU h = I if ν h = 0. Moreover, reasoning similarly for U dec , whose only possible nonzero eigenvalues are ±iω m , we get

e tU dec -I L(X) < 2cε.
Finally, for every c ′ > 0 there exists a choice of t as in (2.6) such that

e tU h -I L(X) < c ′ ε, for h ∈ {1, . . . , ℓ -1, ℓ + 1, . . . , N (m), dec}.
(2.10)

Step 4 Define Σ =   N (m) j=1,j =ℓ exp( t U j )   exp( t U dec ) exp( t U ρ ).
Then (2.4) follows from the commutativity of the terms in the decomposition (2.5) and from the identity (2.7).

By hypothesis, Y ⊂ Ker (U ρ ) and therefore exp(τ U ρ )φ = φ for φ ∈ Y . Then estimate (2.10) with c ′ = 1/N (m) implies that Σ satisfies

(Σ -I)| Y L(Y,X) ≤ j =ℓ e tU j -I L(Y ) + e tU dec -I L(X) < N (m) ε N (m) = ε.
Based on the unitarity of the evolution, we next provide an estimate for the tracking error made by iterated approximations. Lemma 2.6 (Iterated Approximations). Let X be a Hilbert space and Y a subspace of X. Take x 0 ∈ Y , x 0 = 1 and let Υ 1 , . . . , Υ N be a sequence of unitary operators such that

x k = Υ k x k-1 ∈ Y for k = 1, . . . , N . Let Υ k be an approximation of Υ k of the form Υ k = Σ k Υ k ,
where Σ k is an unitary operator on X that satisfies

(Σ k -I)| Y L(Y,X) < ε k , for some ε k > 0.
Then the approximate trajectory x0 = x 0 , xk = Υ k xk-1 satisfies

x n -xn < n k=1 ε k ,
for n = 1, . . . , N.

Proof. We prove the lemma by induction. For x1 , we have:

x1 = Σ 1 Υ 1 x 0 = Υ 1 x 0 + (Σ 1 -I)Υ 1 x 0 = x 1 + δ 1
where δ 1 X < ε 1 because Υ 1 x 0 ∈ Y . Now suppose that for some n ≤ N -1 we known that xn = x n + δ n , with δ n < n k=1 ε k . We have to prove that the same is valid for n + 1.

xn+1 = Σ n+1 Υ n+1 x n = Υ n+1 x n + (Σ n+1 -I)Υ n+1 x n + Σ n+1 Υ n+1 δn = x n+1 + δn+1
where δ n+1 X < ε n+1 + δ n X < n+1 k=1 ε k . Therefore we obtain

x n -xn X < n k=1 ε k , n ≤ N.
This concludes the proof of Lemma 2.6.

Control of the Two Trapped Ions System

In this section we present the decoupled modal decomposition of the two trapped ions system and we use it to prove the approximate controllability of the latter.

Galerkin Representation

For convenience, let us rewrite here the two trapped ions system introduced in (1.1), that is,

i d dt φ gg = (u 1 + u 1r a † + u 1b a)φ eg + (u 2 + u 2r a † + u 2b a)φ ge , i d dt φ eg = (u * 1 + u * 1r a + u * 1b a † )φ gg + (u 2 + u 2r a † + u 2b a)φ ee , i d dt φ ge = (u 1 + u 1r a † + u 1b a)φ ee + (u * 2 + u * 2r a + u * 2b a † )φ gg , i d dt φ ee = (u * 1 + u * 1r a + u * 1b a † )φ ge + (u * 2 + u * 2r a + u * 2b a † )φ eg , and recall that a = 1 √ 2 (x + ∂ x ), a † = 1 √ 2 (x -∂ x ).
Choose as orthonormal basis of (L 2 (R)) 4 the family {φ j } ∞ j=1 defined by

φ 4j+1 =     |j 0 0 0     , φ 4j+2 =     0 |j 0 0     , φ 4j+3 =     0 0 |j 0     , φ 4j+4 =     0 0 0 |j     , (3.1)
where |j denotes the j-th Hermite function, j = 0, 1, 2, . . . The orthonormal basis identifies a system of coordinates in the space of wavefunctions ψ ∈ (L 2 (R)) 4 through the relation

ψ = ∞ n=0     c n gg     |n 0 0 0     + c n eg     0 |n 0 0     + c n ge     0 0 |n 0     + c n ee     0 0 0 |n         . (3.2)
It is useful to split the controls in their real and imaginary parts in order to write a system with real-valued controls:

u 1 = v 1 +iw 1 , u 2 = v 2 +iw 2 , u 1r = v 1r +iw 1r , u 2r = v 2r +iw 2r , u 1b = v 1b +iw 1b , u 2b = v 2b +iw 2b . (3.3) For j, k ∈ N, define the skew-adjoint operators E j,k , F j,k : (L 2 (R)) 4 → (L 2 (R))
4 by their actions on the basis {φ j } ∞ j=1 as follows

E j,k φ j = iφ k , E j,k φ k = iφ j , F j,k φ j = -φ k , F j,k φ k = φ j , E j,k φ ℓ = 0, F j,k φ ℓ = 0, for ℓ / ∈ {j, k}. (3.4)
We can rewrite system (1.1) as

d dt φ = v 1 V 1 + w 1 W 1 + v 1r V 1r + w 1r W 1r + v 1b V 1b + w 1b W 1b (3.5) +v 2 V 2 + w 2 W 2 + v 2r V 2r + w 2r W 2r + v 2b V 2b + w 2b W 2b φ, where V 1 = - ∞ n=0 (E 4n+1,4n+2 + E 4n+3,4n+4 ) , W 1 = ∞ n=0 (F 4n+1,4n+2 + F 4n+3,4n+4 ) , V 1r = - ∞ n=0 √ n + 1 (E 4n+2,4n+5 + E 4n+4,4n+7 ) , W 1r = ∞ n=0 √ n + 1 (F 4n+2,4n+5 + F 4n+4,4n+7 ) , V 1b = - ∞ n=0 √ n + 1 (E 4n+1,4n+6 + E 4n+3,4n+8 ) , W 1b = - ∞ n=0 √ n + 1 (F 4n+1,4n+6 + F 4n+3,4n+8 ) , V 2 = - ∞ n=0 (E 4n+1,4n+3 + E 4n+2,4n+4 ) , W 2 = ∞ n=0 (F 4n+1,4n+3 + F 4n+2,4n+4 ) , V 2r = - ∞ n=0 √ n + 1 (E 4n+1,4n+7 + E 4n+2,4n+8 ) , W 2r = ∞ n=0 √ n + 1 (F 4n+1,4n+7 + F 4n+2,4n+8 ) , V 2b = - ∞ n=0 √ n + 1 (E 4n+3,4n+5 + E 4n+4,4n+6 ) , W 2b = - ∞ n=0 √ n + 1 (F 4n+3,4n+5 + F 4n+4,4n+6
) .

(3.6)

Control of the Modal Approximations

Let us study the controllability of modal (or Galerkin) approximations of system (3.6), obtained by truncating the high energy levels of the system in order to obtain a finite-dimensional reduction.

For every n ∈ N, let Y n = span{φ j | j = 1, . . . , n} ⊂ (L 2 (R)) 4 . The modal approximation of order n of the two trapped ions system is the control system in Y 4n given by

d dt φ = v 1 V (4n) 1 + w 1 W (4n) 1 + v 1r V (4n) 1r + w 1r W (4n) 1r + v 1b V (4n) 1b + w 1b W (4n) 1b +v 2 V (4n) 2 + w 2 W (4n) 2 + v 2r V (4n) 2r
+ w 2r W (4n) 2r

+ v 2b V (4n) 2b + w 2b W (4n) 2b φ, (3.7)
with coupling operators defined by

Z (4n) γ = Π[Y 4n ]Z γ , Z (4n) γ⋆ = Π[Y 4n ]Z γ⋆ , Z = V, W, γ = 1, 2, ⋆ = b, r. (3.8)
By a useful abuse of notation, the operators Z

(4n) γ , Z (4n)
γ⋆ are identified in (3.7) with operators in L(Y 4n ), while equation (3.8) actually defines operators in L((L 2 (R)) 4 ).

Up to the permutation in the coordinates of Y 4n defined by

P (c 0 gg , c 0 eg , c 0 ge , c 0 ee , . . . , c n-1 gg , c n-1 eg , c n-1 ge , c n-1 ee ) = c 0 gg , . . . , c n-1 gg , c 0 eg , . . . , c n-1 eg , c 0 ge , . . . , c n-1 ge , c 0 ee , . . . , c n-1 ee , (3.9)
the coupling operators have the following matrix representations P V

(4n) 1

P -1 = -i     I n I n I n I n     , P W (4n) 1 P -1 =     I n I n I n I n     , P V (4n) 1r P -1 = -i     D T D D T D     , P W (4n) 1r P -1 =     D T D D T D     , P V (4n) 1b P -1 = -i     D D T D D T     , P W (4n) 1b P -1 =     D D T D D T     , P V (4n) 2 P -1 = -i     I n I n I n I n     , P W (4n) 2 P -1 =     I n I n I n I n     , P V (4n) 2r P -1 = -i     D D D T D T     , P W (4n) 2r P -1 =     D D D T D T     , P V (4n) 2b P -1 = -i     D T D T D D     , P W (4n) 2b P -1 =     D T D T D D     , (3.10)
where the matrix D is defined as Remark 3.2. The boundedness on the controls implies that the infimum of the controllability time is positive (at it happens for the single ion case, see [START_REF] Law | Arbitrary control of a quantum electromagnetic field[END_REF]).

D =        0 √ 1 0 √ 2 . . . . . . 0 √ n -1 0        . ( 3
The proof is based on the classical Chow-Rashevskii theorem, recalled here below. (For a proof, see, e.g., [START_REF] Jurdjevic | Geometric control theory[END_REF].) Recall that, given a set F of smooth vector fields on a manifold M , Lie x F denotes the evaluation at a point x ∈ M of the Lie algebra LieF generated by F. Theorem 3.3. Assume that a driftless control affine system ẋ = m j=1 u j X j (x), (u 1 , . . . , u m ) ∈ U ⊂ R m , defined on a finite-dimensional compact manifold M satisfies Lie x {X 1 , . . . , X m } Figure 1: Transition diagram between energy levels of the system by considering a finite number of Fock states in the oscillator. The lines of the same type represent the action of the same control = T x M for all x in M and that the convex hull of U contains the origin in its interior. Then there exists T 0 > 0 such that for any T ≥ T 0 and any open connected set Ω in M , any two points of Ω can be connected by a trajectory x : [0, T ] → M of the control system with support in Ω.

The proof of Proposition 3.1, which can be found in the appendix, consists then in the computation of the iterated Lie brackets of the vector fields appearing in system (3.7). We show in the appendix, moreover, that it is possible to set some of the controls appearing in (3.7) identically equal to zero and still recover the controllability of the modal approximation (for every n ≥ 3).

The Decoupled Modal Approximation

The inconvenient of using the modal approximation (3.7) to study system (1.1) is that whenever we apply the controls associated with the blue lasers (u 1b and u 2b ) or the red lasers (u 1r and u 2r ), some population is transmitted from the phonon level |n to the phonon level |n + 1 (see Figure 1). The truncation defining the modal approximation does not keep track of this transfer, hence the control laws computed on the model approximations do not work for system (1.1). To overcome this issue we use a different truncation, the decoupled modal approximation of system (1.1), which depends on the spectrum of the coupling operators. Let us compare this approximation with the one proposed in [START_REF] Boscain | Multi-input Schrödinger equation: controllability, tracking, and application to the quantum angular momentum[END_REF]: the goal of both approximations is to guarantee that the admissible motions of the approximate system are also admissible, in an approximate way, in the full system. However, the technical arguments leading to the two approximations are different: in [START_REF] Boscain | Multi-input Schrödinger equation: controllability, tracking, and application to the quantum angular momentum[END_REF] (see also [START_REF] Boscain | A weak spectral condition for the controllability of the bilinear Schrödinger equation with application to the control of a rotating planar molecule[END_REF][START_REF] Chambrion | Periodic excitations of bilinear quantum systems[END_REF][START_REF] Chambrion | Controllability of the discrete-spectrum Schrödinger equation driven by an external field[END_REF]) the approximation is based on the dynamical and spectral analysis of the drift Hamiltonian (exploiting the non-resonances of its spectrum), while in the case considered here the drift is null and the non-resonances are exploited separately for different values of the control parameters and then composed together.

The following proposition lists the main spectral properties of the operators defined in (3.6). Proposition 3.4 (Spectral Properties of the Coupling Operators). Let U be one of the operators V 1 , W 1 , V 2 , W 2 . Then the spectrum of U satisfies a) the eigenvalues of U are i and -i;

b) for each n ∈ N the space Y 4n is invariant under U . If U is one of the operators V 1r , W 1r , V 1b , W 1b , V 2r , W 2r , V 2b , W 2b then the spectrum of U satisfies a ′ ) the eigenvalues of U are {±i √ j} ∞ j=0 ; b ′ ) span{x | U x = λx, |λ| < √ n} ⊂ Y 4n , n ∈ N; c ′ ) span{x | U x = λx, |λ| > √ n} ⊂ Y 4n ⊥ , n ∈ N.
The proof follows directly from the definition of the operators given in (3.6), noticing that if

U is one of the operators V 1 , W 1 , V 2 , W 2 then each space span{φ j | j = 4n+1, . . . , 4(n+1)}, n ∈ N, is invariant under U , while if U is one of the operators V 1r , W 1r , V 1b , W 1b , V 2r , W 2r , V 2b , W 2b then span{x | U x = λx, |λ| = √ n} ⊂ span{φ j | j = 4(n -1) + 1, . . . , 4(n + 1)} for every n ∈ N.
We apply in the following the construction of Section 2.1 to the operators V 1b , W 1b , V 1r , W 1r , V 2b , W 2b , V 2r , W 2r Let m ∈ N and write ω j = √ j -1 for j ≥ 1. Let us associate with {ω j } ∞ j=1 the integer N (m) and the frequencies ν j , i = 1, . . . , N (m), as detailed in Section 2.1.

Given n ∈ N, set m = n + 1 and define the decoupled modal approximation of n of the two trapped ions system as the control system in Y 4n with 8N (m) + 4 controls given by

d dt φ = v 1 V (4n) 1 + w 1 W (4n) 1 + v 2 V (4n) 2 + w 2 W (4n) 2 + N (m) j=1 v 1r,j V (4n) 1r,j + w 1r,j W (4n) 1r,j + v 1b,j V (4n) 1b,j + w 1b,j W (4n) 1b,j + v 2r,j V (4n) 2r,j + w 2r,j W (4n) 2r,j + v 2b,j V (4n) 2b,j + w 2b,j W (4n) 2b,j φ, (3.12)
where Z (4n) γ , for Z = V, W and γ = 1, 2, are defined as in (3.8) and Let M > 0 and (φ 0 , φ T ) ∈ C 4n × C 4n , n ∈ N, with φ 0 = φ T . Then, there exists T 0 > 0 such that for any T ≥ T 0 we can find controls z γ , z γ⋆,j ∈ L ∞ ((0, T ), R), z = v, w, γ = 1, 2, ⋆ = b, r, j = 1, . . . , N (n + 1), with L ∞ -norm smaller than M such that the corresponding solution of system (3.12) with initial data φ(0) = φ 0 satisfies φ(T ) = φ T .

Z (4n) γ⋆,j = Z γ⋆ Π[A R m (ν j ) ], Z = V, W, γ = 1, 2, ⋆ = b, r. ( 3
The proposition follows easily from Proposition 3.1, since the set of controlled operators in (3.7) is contained in the one in (3.12), as it follows from (3.14).

Remark 3.6. Clearly, it is also true that, if

U ⊂ {z γ , z γ⋆ | z = v, w, γ = 1, 2, ⋆ = b, r}
is a family of controls which makes (3.7) controllable for every n ≥ 3 (i.e., (3.7) is controllable if we set to zero all the controls which are not in U ) then the corresponding decoupled modal approximation (i.e., system (3.12) where all z γ which are not in U are set to zero, together with all z γ⋆,j such that z γ⋆ is not in U ) is also controllable. Remark 3.7. It is well-know for finite-dimensional control systems that the controls in (3.5) can be taken piecewise constant. Moreover, applying Chow-Rashevskii theorem (see Theorem 3.3) to the control system having as admissible vector fields ±Z (4n) γ and ±Z (4n) γ⋆,j , with Z = V, W , γ = 1, 2, j = 1, . . . , N (p), and ⋆ = b, r, one deduces that φ 0 can be steered to φ T by the concatenation of the flows of such vector fields. Equivalently said, (3.5) is controllable by piecewise constants controls such that at each time instant at most one of them is nonzero.

Main Theorem

Theorem 3.8 (Approximate Controllability of the Two Trapped Ions System). Let ε > 0, M > 0 and φ 0 , φ T ∈ (L 2 (R)) 4 , with φ 0 = φ T . Then, there exists T 0 > 0 such that for all T > T 0 we can find piecewise constant controls u 1 , u 1r , u 1b , u 2 , u 2r , u 2b ∈ L ∞ ((0, T ), C), with norm smaller than M such that the solution of system (1.1) with initial data φ(0) = φ 0 satisfies φ(T ) -φ T (L 2 (R)) 4 < ε.

Proof. First notice that it is enough to prove the theorem for φ 0 , φ T ∈ Y 4n for any given n ∈ N . Indeed, assume that n is large enough so that there exist φ0 , φT ∈ Y 4n with

φ 0 = φ0 = φT = φ T , φ0 -φ 0 (L 2 (R)) 4 < ε 3 , φT -φ T (L 2 (R)) 4 < ε 3 .
Assuming that the theorem is true in the case of initial and final data in Y 4n implies that there exists an admissible trajectory φ(•) of (1.1) such that

φ(0) = φ0 , φ(T ) -φT (L 2 (R)) 4 < ε 3 .
The conclusion then follows from the unitarity of the flow of (1.1). Assume then that φ 0 , φ T ∈ Y 4n , n ∈ N. Let p ≥ n be a prime number and define Y = Y 4p . Consider the decoupled modal approximation (3.12) of order p. By Proposition 3.5 we know that it is controllable. Following Remark 3.7, there exists a sequence of times t 1 , . . . , t N > 0, amplitudes u 1 , . . . , u N ∈ [-M, M ], and coupling operators S 1 , . . . , S N ∈ {Z

(4p) γ , Z (4p) γ⋆,j | Z = V, W, γ = 1, 2, j = 1, . . . , N (p + 1), ⋆ = b, r} such that φ T = e t N u N S N • • • e t 1 u 1 S 1 φ 0 .
In order to mimic the control scheme of the finite-dimensional system in system (1.1), we are going to replace the controls of the decoupled modal approximation by those given by Theorem 2.5 and then study how close the final state φ(T ) is to the target state φ T .

For a given S k , k = 1, . . . , N , define

Υ k = e t k u k S k
and consider the following two cases.

• Either S k is of the type Z (4p) γ

, with Z ∈ {V, W }, γ ∈ {1, 2}. Then we drive the original system (3.5) by the corresponding coupling operator Z γ with the same control u k and for the same time t k (i.e., we set the corresponding control z γ at the value u k , all the other controls to zero, and we follow the flow of the system for a time t k ).

We are then approximating the flow Υ k by an admissible flow Υ k = e t k u k Zγ of (3.5) which can be written as Σ k Υ k with Σ k unitary and Σ k | Y = I. The latter fact follows from Proposition 3.4, point b).

• Or S k is of the type Z (4p) γ⋆,j , with Z ∈ {V, W }, γ ∈ {1, 2}, ⋆ ∈ {b, r}, j ∈ {1, . . . , N (p+1)}. Then we apply Theorem 2.5 with m = p + 1, U = Z γ⋆ , ℓ = j, t = u k t k , and taking as ε the quantity ε/N . Let us check that the hypotheses of the theorem are satisfied. On the one hand, condition (2.3) on Y is verified thanks to Proposition 3.4, points b ′ ) and c ′ ) (taking n = p). On the other hand, in order to check that ω h ωm ∈ (R \ Q) ∪ {0} for each h < m, let us assume by contradiction that there exists 1 < h < p such that

√ h √ p = a b
with a, b ∈ N relatively prime. Therefore, p = h b 2 a 2 , which implies that a 2 divides to h. Then we can write h = sa 2 , with s ∈ N, and we have p = sb 2 . However, as p is prime, this implies that b = 1 and s = p. Hence, h = pa 2 . Since h ≤ p -1, however, this leads to a contradiction. The hypotheses of Theorem 2.5 are then satisfied.

As a consequence of Theorem 2.5 we deduce that there exists t ∈ R such that

exp( tZ γ⋆ ) = Σ k exp(u k t k Z γ⋆ Π[A R p+1 (ν j ) ]) = Σ k Υ k ,
where Σ k is a unitary operator that satisfies

(Σ k -I)| Y L(Y,(L 2 (R)) 4 ) < ε N .
The strategy is then to set in the original system (3.5) the control z γ⋆ corresponding to Z γ⋆ at the value u ′ k = ±u k , all other controls to zero, and let the system flow for a time τ k > 0, where the choice of τ k and of the sign of u ′ k are such that u ′ k τ k = t. Hence Υ k = exp( tZ γ⋆ ) is an admissible flow for system (3.5).

     . (4.4)
Special types of matrix appearing in (4.4) are diag(A (n) , B (n) , A (n) , B (n) ), diag(C (n) , C (n) , D (n) , D (n) ), for A (n) , B (n) , C (n) , D (n) ∈ su(n). Taking brackets between them (with A (n) = 0 (n) or B (n) = 0 (n) and C (n) = 0 (n) or D (n) = 0 (n) ) and since [su(n), su(n)] = su(n), one deduces that, (n) , B (n) , C (n) , D (n) ) ∈ LieF (4.5)

for every A (n) , B (n) , C (n) , D (n) ∈ su(n).

Another special type of matrices appearing in (4.4) are

    0 (n) A (n) 0 (n) 0 (n) -(A (n) ) † 0 (n) 0 (n) 0 (n) 0 (n) 0 (n) 0 (n) A (n) 0 (n) 0 (n) -(A (n) ) † 0 (n)     (4.6)
with A (n) ∈ gl(n). Taking brackets of matrices in (4.5) and in (4.6), and exploiting the fact that su(n)gl(n) = gl(n), we conclude that all matrices of the type    

0 (n) A (n) 0 (n) 0 (n) -(A (n) ) † 0 (n) 0 (n) 0 (n) 0 (n) 0 (n) 0 (n) B (n) 0 (n) 0 (n) -(B (n) ) † 0 (n)     (4.7)
with A (n) , B (n) ∈ gl(n) are in LieF. Taking linear combinations of matrices of the type (4.7) with those of the type (4.5), it turns out that LieF contains all matrices of the type diag(A (2n) , B (2n) ) ∈ LieF, A (2n) , B (2n) ∈ su(2n).

We are left to prove that 0 2n C 2n -(C 2n ) † 0 2n is in LieF for C 2n ∈ gl(2n). Because of (4.4) and (4.5), LieF contains all brackets between matrices of the type diag(A (n) , B (n) , 0 (2n) ) and    

0 (n) 0 (n) C (n) 0 (n) 0 (n) 0 (n) 0 (n) C (n) -(C (n) ) † 0 (n) 0 (n) 0 (n) 0 (n) -(C (n) ) † 0 (n) 0 (n)     , for A (n) , B (n) ∈ su(n), and
C (n) ∈ gl(n). Hence     0 (n) 0 (n) A (n) 0 (n) 0 (n) 0 (n) 0 (n) B (n) -(A (n) ) † 0 (n) 0 (n) 0 (n) 0 (n) -(B (n) ) † 0 (n) 0 (n)     ∈ LieF, for all A (n) , B (n) ∈ gl(n). (4.8)
Taking brackets between matrices of the type (4.7) and (4.8), we easily deduce that    

0 (n) 0 (n) 0 (n) C (n) 0 (n) 0 (n) D (n) 0 (n) 0 (n) -(D (n) ) † 0 (n) 0 (n) -(C (n) ) † 0 (n) 0 (n) 0 (n)     ∈ LieF, for every C (n) , D (n) ∈ gl(n).
This, together with (4.8), completes the proof of the equality LieF = su(4n).

. 11 )

 11 Proposition 3.1 (Exact Controllability of the Modal Approximations of the Two Trapped Ions System). Let n ≥ 3, M > 0 and (φ 0 , φ T ) ∈ Y 4n × Y 4n , with φ 0 = φ T . Then, there exists T 0 > 0 such that for any T ≥ T 0 we can find controls in L ∞ ((0, T ), C) with L ∞ -norm smaller than M such that the solution of system (3.7) with initial data φ(0) = φ 0 satisfies φ(T ) = φ T .

. 13 )

 13 Notice that, byProposition 3.4 point b ′ ), since √ n = ω n+1 = ω m , the operators Z (4n)γ⋆,j are indeed well-defined operators in L(Y 4n ). By Proposition 3.4 point c ′ ), moreover,Z (4n) γ⋆ = N (m) j=1 Z (4n) γ⋆,j , Z = V, W, γ = 1, 2, ⋆ = b, r. (3.14) Proposition 3.5 (Exact Controllability of the Decoupled Modal Approximation).

Acknowledgments

The authors would like to thank Eduardo Cerpa and Alberto Mercado for many fruitful discussions.

The first author was supported in part by scholarship CONICYT Magíster Nacional 2012, project FONDE-CYT 1120610 and project CONICYT ACT-1106.

‡ INRIA, Team GECO, Centre de Recherche Saclay -Île-de-France, France, mario.sigalotti@inria.fr § The second author was partially supported by the European Research Council, ERC StG 2009 GeCoMethods, contract number 239748 and by the iCODE institute, research project of the Idex Paris-Saclay.

1

With the previous procedure we have obtained a control strategy for the two trapped ions system. Finally, applying Lemma 2.6 with X := (L 2 (R)) 4 , x 0 = φ 0 x j = Υ k x j-i , we obtain

Therefore system (1.1) is approximately controllable.

Remark 3.9. It follows from the argument used to prove the theorem and from Remark 3.6 that if U ⊂ {z γ , z γ⋆ | z = v, w, γ = 1, 2, ⋆ = b, r} is a family of controls which makes (3.7) controllable for every n ≥ 3 then the same controls are also sufficient to approximately control system (3.5) (or, equivalently, system (1.1) up to identification of real and imaginary parts as in (3.3)). The appendix discusses which conditions on the set U guarantee such a controllability property.

Conclusions

In this work we have studied a new decomposition method for Schrödinger systems based on spectral techniques. The proposed decomposition allows to obtain approximate controllability results by using finite-dimensional techniques. The method provides satisfactory theoretical results (sufficient conditions for approximate controllability), although it requires large times for the decoupling procedure. The method is applied to the two trapped ions model, for which we obtain a new approximate controllability result. Since the underlying approximate controllability result is based on constructive considerations, the result on the two trapped ions model actually provides a motion planning algorithm (as detailed in [START_REF] Paduro | Approximate Controllability for the Two Trapped Ions System[END_REF]).

Appendix: Controllability of modal approximations

The scope of this appendix is to present the Lie algebra computations necessary for the proof of Proposition 3.1 and, more generally, to identify subfamilies of controls using which it is possible to control the modal approximation (3.7) for every n sufficiently large.

As recalled above, the key criterion allowing the controllability analysis of (3.7) is the Chow-Rashevskii theorem (Theorem 3.3). A useful approach is to apply such a criterion to the lift of system (3.7) in SU(4n) and to exploit the structure of homogeneous space of S 8n-1 (the unit sphere in C 4n ) with respect to SU (4n). Based on such lifting procedure, Albertini and D'Alessandro proved in [START_REF] Albertini | Notions of controllability for bilinear multilevel quantum systems[END_REF] a result implying that (3.7) is controllable if and only if the Lie algebra generated by {Z γ , Z γ⋆ | Z = V, W, γ = 1, 2, ⋆ = r, b} is equal to su(4n) or contains a subalgebra conjugate to sp(2n).

The main result of the appendix, which implies Proposition 3.1 as a particular case, is the following. A preliminary step in the proof of Proposition 4.1 is the analysis of the controllability of the modal approximation of order n of the Law-Eberly system, namely, the control system in

where the controls v, w, v, w r , v b , w b , are real-valued and the coupling operators are defined by

where the matrix D is defined as in (3.11).

Proposition 4.2. Let n ≥ 2 and assume that

Then the Lie algebra generated by F EL is equal to su(2n).

Proof. The proof is contained in [START_REF] Yuan | Controllability of the coupled spin-1 2 harmonic oscillator system[END_REF] in the case ⋆ = r. The case ⋆ = b can be treated in complete analogy, since by a simple reordering of coordinates one can transform

, while preserving V (2n) and W (2n) . The same arguments as in [START_REF] Yuan | Controllability of the coupled spin-1 2 harmonic oscillator system[END_REF] then allow to conclude.

Proof of Proposition 4.1. First of all let us fix the following notation: in order to ease the reading of the proof, we add here below an index to each square matrix indicating its size, so that A (k) denotes a general k × k matrix and 0 (k) the k × k null matrix. We also write diag(A (k 1 ) , . . . , A (kr) ) to denote the square matrix of size k 1 + • • • + k r having A (k 1 ) , . . . , A (kr) as block-diagonal terms. Let us apply Proposition 4.2 to each of the two subsystems of (3.7) obtained by setting to zero either all controls of the type z 1 , z 1⋆ or all controls of the type z 2 , z 2⋆ . In each of the two cases, we recover two decoupled copies of the Law-Eberly modal approximation of order n, controlled simultaneously by the same controls. Proposition 4.2 then implies that for any choice of


domain: math spi

text:
Philippe Destuynder
email: philippe.destuynder@cnam.fr

Caroline Fabre
email: caroline.fabre@u-psud.fr

Singularities

Keywords: wave guides, transmission, singularities, transparent boundary conditions Mathematics Subject Classification 2000: 22E46, 53C35, 57S20

published or not. The documents may come

Singularities occuring in multimaterials with transparent boundary conditions

Introduction

The theory of transparent boundary conditions is a large and deep domain of mathematical studies. To our knowledge, the first main theoretical paper on the subject concerns simulations of solutions of a linear wave equation in the exterior of a bounded domain and is due to B. Engquist and A. Majda ( [START_REF] Engquist | Absorbing Boundary Conditions for the Numerical Simulation of Waves[END_REF]) : the authors suggest a method, based on a Fourier transform in time and the transverse direction, which leads to exact and non local transparent boundary conditions. They also perform several series developments (with different orders) in order to obtain approximate but local transparent boundary conditions. Let us also notice the works of L. Halpern ([10] [START_REF] Halpern | Well-Posedness of One-Way Wave Equations and Absorbing Boundary Conditions[END_REF]) : she (and co-authors) studied well-posedness of different boundary conditions and their numerical schemes.

When one studies solutions of linear wave equations with a unique velocity, the boundary conditions on the input and output boundary are transparent boundaries conditions : they ensure that waves should go out of the domain (there is no reflection) and they avoid singularity (see Fig. 2). In case of a transmission problem, the difference between wave'velocities involves singularities which are localized at the intersections between the boundary and the interface of transmission (see Fig. 3 and Fig. 4). We give a precise description of them with a mathematical analysis of this phenomenon : we prove that the singularities are involved by the gap across the interface on the lateral boundary of the antisymmetric part of the Neumann derivative of the solution. Let us notice that these singularities do not appear for homogeneous Neumann boundary conditions on the whole boundary of Ω.

Let us illustrate our claim with the following different situations : figures 1, 2, 3 and 4 below concern solutions of a transmission problem (see system (1.1)) in the square Ω =]0, 1[×]0, 1[ and during the time T = 1.2. The interface is localized at y = 0.2 and the velocities are c = c -if y < 0.2 and c = c + if y > 0.2 with 0 < c -≤ c + . Initial data are null and the right hand side f is localized in the region with velocity c -. Fig. 1 shows the support of the right hand side f . Numerical computations are performed in Matlab with a time step of order 0.001 and space steps (both in the x and y directions) of 1 120 = 0.0083. Figures 2, 3 and 4 are snap-shots at time T = 0.41294 or T = 0.41367 in the three following situations : in Fig. 2, one has c -= c + = 1 and the right hand side is a high frequencies time excitation. One can see that there is no singularity. In Fig. 3, one has c + = 2 and c -= 1 and the right hand side f is a low frequencies time excitation. In Fig. 4, one has c + = 2 and c -= 0.5 and the right hand side is a high frequency time excitation. In both Fig. 3 and4, one can see singularities appearing at the intersection between the interface y = 0.2 and the the boundary x = 0. This is what we explain in the following. Others numerical computations are given in order to exhibit more precisely these singularities. In a forthcoming paper, we suggest a way to avoid these singularities introducing efficient new tranparent boundary conditions. Let us notice that our method leads to a non local computation which is classical in the obtention of exact transparent boundary conditions at higher order but new for first order.

Our plan is the following one : in a first section, we present the mathematical problem and we focus on the singularity. In a second section, we study existence and regularity results of its solutions. In a third section, we focus on singularities on the boundary of the interface. In this section, we state our main theorem. All along, the paper is completed by numerical simulations which have been performed with Matlab and which illustrate our theoritical results.

The wave model used for the discussion

Let us consider a two dimensional open set Ω =]0, L[×] -a, a[ (L > 0 and a > 0) of R 2 as shown on figure 5. A point x of Ω has coordinates x = (x 1 , x 2 ). The wave velocity is denoted by c and is piecewise constant. It is c + in Ω + = Ω ∩ (x 2 > 0) and c -in Ω -= Ω ∩ (x 2 < 0) and we assume that 0 < c -< c + .

For any given functions f = f (x, t), u 0 = u 0 (x) and u 1 = u 1 (x), let us consider

                               ∂ 2 u ∂t 2 -div(c 2 ∇u) = f in Q = Ω×]0, T [, ∂u ∂ν = 0 on (Γ + ∪ Γ -)×]0, T [, ∂u ∂t + c ∂u ∂ν = 0 on (Γ e ∪ Γ s )×]0, T [, u(x, 0) = u 0 (x) ∂u ∂t (x, 0) = u 1 (x) in Ω. (1.1) We write Γ e+ = Γ e × (x 2 > 0), Γ e-= Γ e × (x 2 < 0), Σ + = Γ + ×]0, T [, Σ -= Γ -×]0, T [, Σ e = Γ e ×]0, T [, and Σ s = Γ s ×]0, T [. The interface is Γ i = Ω ∩ [x 2 =
0] and Σ i = Γ i ×]0, T [. The letter ν denotes the unit outward normal vector on the boundary Γ of Ω.

Existence and uniqueness of solutions of system (1.1) are classical results, even if the transparent boundary conditions on Γ e ∪ Γ s (which imply first order time derivative) require a slightly different strategy as the usual one. Let us summarize them in the following statement.

Proposition 1.1. Let us assume that u 0 ∈ H 1 (Ω), u 1 ∈ L 2 (Ω) and f ∈ L 2 (Q).
Then, there exists a unique solution u to the system (1.1) with :

u ∈ L ∞ (]0, T [; H 1 (Ω)) ∩ W 1,∞ (]0, T [; L 2 (Ω)).
We don't give the full proof of proposition 1.1, one can refer to [START_REF] Dautray | Analyse mathématique et calcul numérique T8. Evolution : Semi-groupe, Variationnel, collection enseignement[END_REF] in case of interest : it is based on a Galerkin method, a priori estimate and weak convergence of a subsequence solution of a finite dimensional approximation. Let us recall the variational formulation of (1.1).

• The function u is solution of

∀v ∈ H 1 (Ω), Ω ∂ 2 u ∂t 2 (x, t)v(x)dx + Ω c 2 ∇u(x, t).∇v(x)dx+ Γe∪Γs c ∂u ∂t (x, t)v(x) = Ω f (x, t)v(x)dx (1.2)
• The energy is defined by

E(t) = 1 2 Ω ∂u ∂t (x, t) 2 dx + 1 2 Ω c 2 |∇u(x, t)| 2 dx (1.3)
and one has for f = 0

E(t) + Σe∪Σs c ∂u ∂t (x, t) 2 dx = E(0).
Therefore it decreases with respect to the time variable. In the general case, there exists a constant d > 0, such that for every data

(u 0 , u 1 , f ) ∈ H 1 (Ω) × L 2 (Ω) × L 2 (Q) : E(t) + Σe∪Σs c| ∂u ∂t (x, t)| 2 dx ≤ d[E(0) + ||f || Q ] 2 ,
where || || Q denotes the L 2 (Q)-norm. These results prove that :

∂u ∂t ∈ L 2 (Σ e ∪ Σ s ). Since ∂u ∂ν = - 1 c ∂u ∂t , one has ∂u ∂ν ∈ L 2 (Σ e ∪ Σ s ). Furthermore, even if ∂u ∂t ∈ H -1 (0, T ; H 1 (Ω)), the discontinuity of c implies that in general ∂u ∂ν ∈ D ′ (0, T ; H 1/2 (Γ e )) and ∂u ∂ν ∈ D ′ (0, T ; H 1/2 (Γ s )).
The function u, unique solution of (1.1), satisfies :

                     ∂ 2 u ∂t 2 -c 2 + ∆u = f in Q + = Ω + ×]0, T [, ∂ 2 u ∂t 2 -c 2 -∆u = f in Q -= Ω -×]0, T [, c 2 + ∂u ∂ν i -c 2 - ∂u ∂ν i = 0 on Σ i , (1.4)
where ν i is one of the unit normal vector to Γ i and ∂u ∂ν i = ∇u.ν i .

Polylog singularities 7

The third equation of (1.4) implies that, whatever the smoothness of the data are, there is in general no hope to have u ∈ D ′ (0, T ; H 2 loc (Ω)) if c is not constant across the separation line Γ i between the two materials.

We now turn to the regularity study of solutions of (1.1).

On the smoothness of u in time and space

Let us study the regularity of u first with respect to the time variable and then with respect to space variables. Singularities will be discussed in the next section.

Based on derivative methods with respect to the time t, one can upgrade the results stated in Proposition 1.1 by assuming that the initial condition u 0 doesn't cross the interface Γ i .

u 0 ∈ H 2 (Ω) with supp(u 0 ) ⊂ O + ∪ O -, u 1 ∈ H 1 (Ω) and if f ∈ H 1 (]0, T [; L 2 (Ω)), then the solution u of (1.1) is such that: u ∈ W 1,∞ (]0, T [; H 1 (Ω)) ∩ W 2,∞ (]0, T [; L 2 (Ω)).
From classical inclusions (see [START_REF] Brezis | Analyse Fonctionnelle, collection mathématiques appliquées pour la maîtrise[END_REF]), this implies for instance that

u ∈ C 0 ([0, T ]; H 1 (Ω)). Proof. Let us set u = ∂u ∂t . The function u is solution of (1.1) with data ḟ = ∂f ∂t ∈ L 2 (Q), u(x, 0) = u 1 (x) ∈ H 1 (Ω) and ∂ u ∂t (x, 0) = f (x, 0) + div(c 2 ∇u 0 ).
The assumption on the initial data u 0 ensures that ∇u 0 is null on a neighborhood of the interface Γ i in Ω and thus div(c

2 ∇u 0 ) ∈ L 2 (Ω). Moreover, f ∈ C 0 ([0, T ]; L 2 (Ω)) and f (x, 0) is in L 2 (Ω). Applying Proposition 1.1, we eas- ily deduce that u ∈ W 1,∞ (]0, T [; H 1 (Ω)) ∩ W 2,∞ (]0, T [; L 2 (Ω)).
By iterating Theorem 2.1, one easily gets Theorem 2.2. With the same notations as in Theorem 2.1, and if the initial conditions are such that

(u 0 , u 1 ) ∈ H 3 (Ω) × H 2 (Ω) both with compact supports in O + ∪ O -, and if f ∈ H 2 (0, T ; L 2 (Ω)), then the solution u of (1.1) is such that: u ∈ W 2,∞ (]0, T [; H 1 (Ω)) ∩ W 3,∞ (]0, T [; L 2 (Ω)) and thus u ∈ C 1 ([0, T ], H 1 (Ω)).
Remark 2.1. Let us first recall that, if γ is a non-empty and open part of the whole boundary Γ of Ω, the restriction operator to γ maps H 1/2 (Γ) onto H 1/2 (γ) whereas the extension by zero of a function of H 1/2 (γ) is not in general a function of H 1/2 (Γ). In both cases of theorems 2.1 and 2.2, one has

u |Γ e ∪Γs ∈ C 1 ([0, T ]; H 1/2 (Γ e ∪ Γ s ))
and therefore

c ∂u ∂ν |Γe∪Γs ∈ C 0 ([0, T ]; H 1/2 (Γ e ∪ Γ s )).
The function c is discontinuous across Γ i and therefore in general (even a piecewise continuous function is not globally H 1/2 of the whole open set)

∂u ∂ν |Γe ∈ D ′ (0, T ; H 1/2 (Γ e )) and ∂u ∂ν |Γs ∈ D ′ (0, T ; H 1/2 (Γ s )). (2.1)
We get at least for 0 ≤ s < 1/2,

∂u ∂ν |Γe ∈ C 0 ([0, T ]; H s (Γ e )) and ∂u ∂ν |Γs ∈ C 0 ([0, T ]; H s (Γ s )).
Let us now turn to the smoothness with respect to the space variables and there are several cases. For x ∈ R 2 , we denote by V x a small enough open and non-empty neighborhood of x in R 2 . We write 1 O the characteristic function of a set O of R 2 and

• Γ i the interior of Γ i . Theorem 2.3. Let f ∈ H 1 (Q)
and let us assume that the initial data as in Theorem 2.1. Let u be the solution of system 1.1. One has :

(1) u1 Ω+ ∈ L ∞ (]0, T [; H 2 loc (Ω + )) and u1 Ω-∈ L ∞ (]0, T [; H 2 loc (Ω -)) (2) If x ∈ Γ + ∪ Γ -then u ∈ L ∞ (]0, T [; H 2 (V x ∩ Ω)) (3) If x ∈ (Γ e ∩Ω + )∪(Γ e ∩Ω -)∪(Γ s ∩Ω + )∪(Γ s ∩Ω -) then u ∈ L ∞ (]0, T [; H 2 (V x ∩ Ω)) (4) If x ∈ • Γ i then ∂u ∂x 1 ∈ L ∞ (]0, T [; H 1 (V x )). Proof. (1) On Ω + ×]0, T [, one has div(c 2 ∇u) = c 2 + ∆(u) = ∂ 2 u ∂t 2 ∈ L ∞ (]0, T [; L 2 (Ω + )).
From a localisation argument, it is classical to prove that u| Ω+ ∈ L ∞ (]0, T [; H 2 loc (Ω + )). Of course, the same argument can be applied in Ω -. (2) Let x ∈ Γ + for example. The boundary conditions ∂u ∂ν = 0 on Σ + allows us to apply a symmetry argument and to consider an extension ū of u across Γ + which is still solution of 1.1 in a neighborhood V x of x. The first point of this theorem leads to ū

∈ L ∞ (]0, T [; H 2 (V x )) and thus u ∈ L ∞ (]0, T [; H 2 (V x ∩ Ω)). (3) Let x ∈ Γ e ∩ Ω + . Since ∂u ∂t ∈ L ∞ (]0, T [; H 1 (Ω)), we obtain ∂u ∂ν = - 1 c ∂u ∂t ∈ L ∞ (]0, T [; H 1/2 (Γ e ∩Ω + )). After localisation, we get ∆u ∈ L 2 ((V x ∩Ω)×(0, T )) and ∂u ∂ν ∈ L ∞ (0, T ; H 1/2 (∂(V x ∩Ω)). Classically, this leads to u ∈ L ∞ (0, T ; H 2 (V x ∩ Ω)). (4) Let x ∈ • Γ i .
In this case, let us notice that we can assume that V x ⊂ Ω. There is no hope (in general) that u ∈ D ′ (0, T ; H 2 loc (V x )) excepted for c + = c -. Let us consider ρ ∈ D(Ω) with ρ = 1 on V x and let us write w 1 = ρ ∂u ∂x 1 . Since Γ i is parallel to the axis and boundary x 1 , the function w 1 is still solution of a transmission problem similar to 1.1 with

     ∂ 2 w 1 ∂t 2 -div(c 2 ∇w 1 ) ∈ L 2 (Q), ∂w 1 ∂ν ∈ L ∞ (]0, T [; H 1/2 (∂V x )),
and thus

w 1 ∈ L ∞ (]0, T [; H 1 (V x )).
Remark 2.2. One can easily prove that

∂u ∂x 1 is L ∞ (]0, T [, H 1 (]a, b[×] -a, a[)) in any rectangle ]a, b[×] -a, a[ with 0 < a < b < L.
We now turn to the main point : the description of the singularities focusing at the point (x 1 , x 2 ) = (0, 0).

Polylog-2 Singularities at the jonction of the interface and the

boundary Γ e (or Γ s )

We focus our study at the point with the coordinates (0, 0) but, of course, an analogous result is valid at the point with coordinates (L, 0). We introduce an even function

ρ with ρ = ρ(x 2 ) ∈ C ∞ (R)    ρ(x 2 ) = 1 for - a 2 < x 2 < a 2 ρ(x 2 ) = 0 for |x 2 | > 3a 4 (3.1) Let us denote W =]0, inf(a, L) 2 [×] - a 2 , a 2 [. We consider a function ρ a ∈ C ∞ (Ω) such that ρ a = 1 in W, ρ a = 0 in a neighborhood of Γ + ∪ Γ -∪ Γ s . (3.2)
Let us set :

V = {w ∈ H 1 (Ω), ρw1 Ω+ ∈ H 2 (Ω + ), ρw1 Ω-∈ H 2 (Ω -); and ∂ρw ∂x 1 ∈ H 1 (W )}. (3.3)
For a given function h, we denote by T s and T a the following symmetric and non-symmetric part of h defined by

T s (h)(x 1 , x 2 ) = 1 2 [h(x 1 , x 2 ) + h(x 1 , -x 2 )], (3.4)
and

T a (h)(x 1 , x 2 ) = 1 2 [h(x 1 , x 2 ) -h(x 1 , -x 2 )]. (3.5)
When no mistake can be made, we write h s and h a instead of T s (h) and T a (h). Of course, h = h s + h a .

We denote by Im(z) the imaginary part of the complex number z and we introduce the two following functions :

Li 2 is the polylogarithm function of order 2 defined by

Li 2 (z) = n≥1 z n n 2 for |z| ≤ 1. With z = exp[- π a (x 1 -ix 2 )],
we consider the function S defined by :

S(x 1 , x 2 ) = 2a π 2 (c 2 + + c 2 -) Im[Li 2 (z) -Li 2 (-z)] [ c - c + 1 (x2>0) + c + c - 1 (x2<0) ].
Let us notice that S ∈ H 1 (Ω) and that S is null on the interface Γ i . The graphs of the function S and its partial derivatives are given on figure 6 in the case where c + = c -= 1. Our main result in the paper is the following one : it gives a decompostion of the solution u in a regular part (u r ) and a singular one (u sg ).

Theorem 3.1. Let f ∈ H 1 (Ω×]0, T [), u 0 ∈ H 1 (Ω) and u 1 ∈ L 2 (Ω). Then, there exist two functions u r and u sg in L ∞ (]0, T [;

H 1 (Ω)) with u = u r + u sg in Q. The functions u r satisfies ∂u r ∂x 1 ∈ H -2 (]0, T [; V).

Polylog singularities 11

Furthermore, there exist a function α 0 = α 0 (t) ∈ H -2 (0, T ) and ε > 0, such that

u sg (x 1 , x 2 , t) = (c + -c -)α 0 (t)S(x 1 , x 2 ) in ]0, ε[×] -ε, ε[. On ]0, ε[×] -ε, ε[, we have ∇u sg (0, x 2 , t) = (c + -c -)[ c - c + 1 (x2>0) + c + c - 1 (x2<0) ]α 0 (t)    - 2 π Arg( 1 + e iπ x 2 a 1 -e iπ x 2 a ) 2 π Ln|cotan πx 2 2a |   (3.6) Remark 3.1.
Let us notice that S is null if c + = c -thus the existence of the singularity is due to the presence of two materials and we then get Waves in a bimaterial with interface at x 2 = 0.2 The end of the paper is devoted to the proof of Theorem 3.1.

lim x2→0 + ∂u sg ∂x 2 (0, x 2 ) = lim x2→0 - ∂u sg ∂x 2 (0, x 2 ) = +∞.
Proof of Theorem 3.1. The idea is the following one : we prove that the singular part of u on Γ e comes from the antisymmetric part of the singular function -1 c ∂u ∂t on the boundary Γ e . More precisely, the singularity is strictly connected to the gap at the origin (0, 0) of this function. In order to point out this fact, we split the solution of (1.1) in several parts and we write g = -

1 c ∂u ∂t 1 Σe∪Σs .
Let us recall that for (u 0 , u

1 ) ∈ H 1 (Ω) × L 2 (Ω) and f ∈ L 2 (Q), we have g ∈ L 2 (Σ) but g ∈ H -1 (]0, T [; H 1/2 (Γ)). In a same way, if (u 0 , u 1 ) ∈ H 2 (Ω) × H 1 (Ω) with supp(u 0 ) ⊂ K where K is a compact set in Ω + ∪ Ω -, even if f ∈ H 1 (Q), then g ∈ L 2 (]0, T [; H 1/2 (Γ))
. Since our interest is not far from the interface Γ i , we first localize the function u introducing ũ = ρu where ρ is defined in (3.1). The function ũ is solution of the following system :

                             ∂ 2 ũ ∂t 2 -div(c 2 ∇ũ) = f in Q = Ω×]0, T [, ũ = 0 on (Γ + ∪ Γ -)×]0, T [, ∂ ũ ∂t + c ∂ ũ ∂ν = 0 on (Γ e ∪ Γ s )×]0, T [, ũ(x, 0) = ρu 0 (x) ∂ ũ ∂t (x, 0) = ρu 1 (x) in Ω.
(3.7)

with f = ρf + c 2 ∇ρ.∇u + div(c 2 ∇ρu). Since ρ = 1 in a neighborhood of Γ i , we have div(c 2 ∇ρu) ∈ L 2 (Q). We denote g = - 1 c ∂ ũ ∂t 1 Σe∪Σs . We get g ∈ L 2 (Σ) and g = g = - 1 c
∂u ∂t in a neighborhood in Γ e of the point (0, 0).

We recall that g s (respectively g a ) is the symmetric (respectively the antisymmetric) part of g. We introduce u (1) ∈ L 2 (0, T ; H 1 (Ω)) and u (2) ∈ L 2 (0, T ; H 1 (Ω)) solutions of :

               div(c 2 ∇u (2) (t)) = 0 in Q = Ω×]0, T [, u (2) (t) = 0 on {Γ + ∪ Γ -}×]0, T [, ∂u (2) ∂ν (t) = ga (t) on {Γ e ∪ Γ s }×]0, T [, (3.8)
and u (1) = ũ -u (2) . We obtain that u (1) ∈ L 2 (0, T ; H 1 (Ω) is solution of

               div(c 2 ∇u (1) ) = div(c 2 ∇ũ) = ü -f in Q, u (1) = 0 on {Γ + ∪ Γ -}×]0, T [, ∂u (1) ∂ν = ḡs on (Γ e ∪ Γ s )×]0, T [, (3.9)
The following lemmas will be useful for the study of the regularity of u (1) .

Lemma 3.1. Let w ∈ H 1 (Ω) be a solution of                    div(c 2 ∇w) = h in Ω, w = 0 on Γ + ∪ Γ -, ∂w ∂ν = g 0 on Γ e ∪ Γ s , (3.10)
with h ∈ L 2 (Ω), g 0 1 Γe ∈ H 1/2 (Γ e ) and g 0 1 Γs ∈ H 1/2 (Γ s ). Then w ∈ V (see (3.3)).

Proof of Lemma 3.1. The function w 1 = ρw is solution of

                   div(c 2 ∇w 1 ) = h 1 in Ω, w 1 = 0 on Γ + ∪ Γ -, ∂w 1 ∂ν = g 1 on Γ e ∪ Γ s , (3.11)
with h 1 = ρh + c 2 ∇ρ.∇w + div(c 2 w∇ρ) and g 1 = ρg 0 + w∇ρ.ν.

Since ∇ρ = 0 in a neighborhood of Γ i , we deduce that c 2 ∇ρ ∈ C ∞ (R) and div(c 2 w∇ρ) ∈ L 2 (Q) and thus h 1 ∈ L 2 (Q). Furthemore, g 0 1 Γe ∈ H 1/2 (Γ e ) implies that g 0 ρ1 Γe ∈ H 1/2 00 (Γ e ) (and same argument is valid on Γ s ). Let us denoted by ḡ1 the extension by zero on the whole boundary of the function g 1 . We then have ḡ1 ∈ H 1/2 (Γ).

Let us set w

2 = ∂w 1 ∂x 1 = ρ ∂w ∂x 1 . We obtain                    div(c 2 ∇w 2 ) = ∂h 1 ∂x 1 in Ω, w 2 = 0 on Γ + ∪ Γ -,
w 2 = -g 1 on Γ e and w 2 = g 1 on Γ s .

(3.12)

Since g 1 1 Γe ∈ H 1/2 00 (Γ e ) (and same on Γ s ), we get that the function 1 Γs g 1 -

1 Γe g 1 ∈ H 1/2 (Γ). Introducing now G 1 ∈ H 1 (Ω) with G 1 = 1 Γs g 1 -1 Γe g 1 on Γ, and w 3 = w 2 -G 1 , we deduce that        div(c 2 ∇w 3 ) = ∂h 1 ∂x 1 -div(c 2 ∇G 1 ) in Ω, w 3 = 0 on Γ, (3.13)
with

h 3 = ∂h 1 ∂x 1 -div(c 2 ∇G 1 ) ∈ H -1 (Ω).
We easily deduce that

w 3 ∈ H 1 (Ω) thus w 2 = ρ ∂w ∂x 1 ∈ H 1 (Ω).
Assertions ρw1 Ω+ ∈ H 2 (Ω + ) and ρw1 Ω-∈ H 2 (Ω -) are easy consequences of ∆(ρw1 Ω+ ) ∈ L 2 (Ω + ) and ∂ρw ∂x 1 ∈ H 1 (Ω) ( same for Ω -). Of course, there is no hope to have w ∈ H 2 (W ). Lemma 3.1 is therefore proved.

Lemma 3.2. Let g ∈ L 2 (Γ e ) with g1 Γe + ∈ H 1/2 (Γ e+ ) and g1 Γe -∈ H 1/2 (Γ e-).
Then, we have

g s ∈ H 1/2 (Γ e ).
Let us suppose that Lemma 3.2 is proved. Since the function g satisfies g1 Γe+ ∈ H -1 (0, T, H 1/2 (Γ e+ )), (and same on Γ e-), we obtain that gs ∈ H -1 (0, T ; H 1/2 (Γ e ). Moreover, üf ∈ H -1 (0, T ; L 2 (Ω)). We can apply Lemma 3.1 to the function u (1) , and we get that u (1)

∈ H -1 (0, T ; V) if (u 0 , u 1 ) ∈ H 1 (Ω)×L 2 (Ω).
In the case where (u 0 , u

1 ) ∈ H 2 (Ω) × H 1 (Ω) with supp(u 0 ) ⊂ K where K is a compact set in Ω + ∪ Ω -, one can easely prove that u (1) ∈ L 2 (0, T ; V) since ü -f ∈ L 2 (0, T ; L 2 (Ω)) gs ∈ L 2 (0, T ; H 1/2 (Γ e )).
Proof of Lemma 3.2. Let us point out that the time variable is a parameter in Lemma 3.2 and thus by linearity, the regularity in time comes from the one of the function g. There is nothing to prove for this.

Let us introduce :

T s (g 1 , g 2 )(x 2 ) =      1 2 [g 1 (x 2 ) + g 2 (-x 2 )] if x 2 > 0 1 2 [g 2 (x 2 ) + g 1 (-x 2 )] if x 2 < 0 for g 1 ∈ L 2 (Γ e+ ) and g 2 ∈ L 2 (Γ e-).
We easily get

T ∈ L(L 2 (Γ e+ ) × L 2 (Γ e-); L 2 (Γ e )). If (g 1 , g 2 ) ∈ H 1 (Γ e+ )×H 1 (Γ e-) then T s (g 1 , g 2 ) ∈ H 1 (Γ e ) (there is no gap across x 2 = 0) hence T ∈ L(H 1 (Γ e+ ) × H 1 (Γ e-); H 1 (Γ e )).
By interpolation of order 1/2, we deduce Lemma 3.2.

We now turn to the proof of Theorem 3.1 with the study of u (2) , solution of (3.8) and we first prove the following result concerning the solution u : it defines the gap of u| Γe at the origin.

Proposition 3.1. Suppose (u 0 , u 1 ) ∈ H 1 (Ω) × L 2 (Ω). We have ũ1 Ω+ ∈ L 2 (0, T ; H 3/2 (Ω + )), ũ1 Ω-∈ L 2 (]0, T [; H 3/2 (Ω -)) and we can write ga = T a (g) = α(t) sign(x 2 ) + g ar where α ∈ H -2 (0, T ) and g ar ∈ H -2 (]0, T [; H 1/2 (Γ e )).
Proof of Proposition 3.1. Let us prove that ũ1 Ω+ ∈ L 2 (]0, T [; H 3/2 (Ω + )). We introduce ũS and ũA defined by ũA = ũ -ũS and

ũS (x 1 , x 2 , t) =            c 2 + ũ(x 1 , x 2 , t) + c 2 -ũ(x 1 , -x 2 , t) c 2 + + c 2 - if x 2 > 0 c 2 + ũ(x 1 , -x 2 , t) + c 2 -ũ(x 1 , x 2 , t) c 2 + + c 2 - if x 2 < 0
We have ũS ∈ L 2 (0, T ; H 1 (Ω) and ũS = ũ on Γ i . Let us prove that ũS ∈ H -1 (0, T ; H 2 (Ω)). We obtain in Ω + :

∂ ũS ∂x 2 = 1 c 2 + + c 2 - [c 2 + ∂ ũ ∂x 2 (x 1 , x 2 , t) -c 2 - ∂ ũ ∂x 2 (x 1 , -x 2 , t)]
thus (recall that ũ satisfies the transmission condition), we get ∂ ũS ∂ν = 0 on Γ i . An analogous calculus on Ω -proves that there is no gap of the normal derivative of ũS through the interface Γ i . Since ũS ∈ L 2 (0, T ;

H 1 (Ω)), ∆(ũ S 1 Ω+ ) ∈ H -1 (0, T ; L 2 (Ω + )), ∆(ũ S 1 Ω-) ∈ H -1 (0, T ; L 2 (Ω -))
, with no gap of the normal derivative through Γ i , we deduce that ∆ũ S ∈ H -1 (0, T ; L 2 (Ω)).

On another hand, we have on

Ω + ∂ ũS ∂x 1 = 1 c 2 + + c 2 - [c 2 + ∂ ũ ∂x 1 (x 1 , x 2 , t) + c 2 - ∂ ũ ∂x 1 (x 1 , -x 2 , t)]
Polylog singularities 17

and on Ω -

∂ ũS ∂x 1 = 1 c 2 + + c 2 - [c 2 - ∂ ũ ∂x 1 (x 1 , x 2 , t) + c 2 + ∂ ũ ∂x 1 (x 1 , -x 2 , t)],
thus ∂ ũS ∂ν is a symmetric function on Γ e and on Γ s . We can apply Lemma 3.2 and we get that

∂ ũS ∂ν ∈ H -1 (0, T ; H 1/2 (Γ e ∪ Γ s )).
We have proved that ũS is solution of

       ∆ũ S ∈ H -1 (0, T ; L 2 (Ω)) ũS = 0 on Γ + ∪ Γ - ∂ ũS ∂ν ∈ H -1 (0, T ; H 1/2 (Γ e ∪ Γ s )).
Since ũS is null in a neighborhood of Γ + and Γ -, we deduce that ũS ∈ H -1 (0, T ; H 2 (Ω)).

Let us now prove that ũA 1 Ω+ = (ũ -ũS )1 Ω+ ∈ H -1 (0, T ; H 3/2 (Ω + )). We have ũA ∈ L 2 (0, T ; H 1 (Ω)). Since ũ = ũS on Γ i , we get ũA | Γi = 0. Moreover, using that the normal direction on Γ e (respectively Γ s ) is the -x 1 's one (respectively x 1 's one), one can easily prove that u A satisfies in Ω + × (0, T )

               ∆ũ A ∈ H -1 (0, T ; L 2 (Ω + )) ∂ ũA ∂ν = ∂ ũ ∂ν - ∂ ũS ∂ν ∈ H -1 (0, T ; H 1/2 (Γ e+ ∪ Γ s+ )) u A = 0 on Γ i ∪ Γ + and therefore ũA ∈ H -1 (]0, T [; H 3/2 (Ω + )) if (u 0 , u 1 ) ∈ H 1 (Ω) × L 2 (Ω). We proved that ũ ∈ H -1 (]0, T [; H 3/2 (Ω + )) + H -1 (]0, T [; H 2 (Ω + )) and thus ũ ∈ H -1 (]0, T [; H 3/2 (Ω + )) for initial data (u 0 , u 1 ) in H 1 (Ω) × L 2 (Ω). Of course, the same is valid in Ω -.
We deduce that

∂ ũ ∂t | Γe + and ∂u ∂t | Γe -make sense in H -2 (]0, T [; H 1 (Γ e+ )
) and

H -2 (]0, T [; H 1 (Γ e-))
and thus their values at point (0, 0) exist. Furthermore,

∂ ũ ∂t ∈ H -1 (]0, T [; H 1/2 (Γ e ))
, thus these values are equal. Let us write α = ga | Γe + be the trace of the function ga on Γ e+ . We have

α(t) = ga (0+, t) = - 1 2 [ 1 c + ∂ ũ ∂t (0, 0+, t) - 1 c - ∂u ∂t (0, 0-, t)] = 1 2c + c - (c + -c -) ∂ ũ ∂t (0, 0, t).
We write ga (x 2 , t) = α(t) sign (x 2 ) + g ar (x 2 , t)

where sign denotes the sign-function defined by sign (x) = 1 if x > 0 and sign (x) = -1 if x < 0. The function g ar is an odd function with respect to x 2 , g ar 1 Γe + ∈ H -2 (0, T ; H 1 (Γ e+ )) and g ar 1 Γe -∈ H -2 (0, T ; H 1 (Γ e-)) and their values at the point (0, 0) are null therefore g ar ∈ H -2 (]0, T [; H 1/2 (Γ e )) and proposition 3.1 is proved.

Let us return to the proof of our main theorem. We write u (2) = u (3) + u (4) where u (3) (respectively u (4) ) is solution of (3.8) with ∂u (3) ∂ν = g ar (respectively ∂u (4) ∂ν = α sign (x 2 )) on Γ e .

Lemma 3.1 can be apply to u (3) and leads to u (3) ∈ V. We deduce that the singular part of the function u is involved by the lack of continuity at the origin of the odd part of the function 1 c ∂u ∂t . Let study w = u (4) .

We use the same splitting as in proposition 3.1 and we write w = w S + w A . We get w S ∈ H -1 (0, T ; H 2 (Ω)) and thus it is sufficient to study u (5) = w A . We know that ∆u (5) = 0 separately in Ω + and Ω -and that u (5) = 0 on Γ i ∪ Γ + ∪ Γ -.

Let us compute

∂u (5) ∂ν = -∂u (5) ∂x 1 on Γ e+ and Γ e-. We have on Γ e+ :

∂u (5) ∂x 1 (x 1 , x 2 ) = ∂u (4) ∂x 1 (x 1 , x 2 ) - ∂w S ∂x 1 (x 1 , x 2 ) = -α - c 2 + ∂u (4) ∂x 1 (x 1 , x 2 ) + c 2 - ∂u (4) ∂x 1 (x 1 , -x 2 ) c 2 + + c 2 - = -α - c 2 + (-α) + c 2 -α c 2 + + c 2 - = -α 2c 2 - c 2 + + c 2 - and thus
β + (t) = ∂u (5) ∂ν = c - c + (c 2 + + c 2 -) (c + -c -) ∂ ũ ∂t (0, 0, t) on Γ e+ . (3.14)
On Γ -, we get

∂u (5) ∂ν = -α - c 2 - ∂u (4) ∂ν (x 1 , x 2 ) + c 2 + ∂u (4) ∂ν (x 1 , -x 2 ) c 2 + + c 2 - = -α - c 2 -(-α) + c 2 + α c 2 + + c 2 - = -α 2c 2 + c 2 + + c 2 - and thus
β -(t) = ∂u (5) ∂ν = - c + c -(c 2 + + c 2 -) (c + -c -) ∂ ũ ∂t (0, 0, t) on Γ e-. ( 3

.15)

Polylog singularities 19

We deduce that the function u (5) is solution of

                   ∆u (5) = 0 in Ω + , u (5) = 0 on Γ + ∪ Γ i , ∂u (5
) ∂ν = β + on Γ e+ ∪ Γ s+ , and
                   ∆u (5) = 0 in Ω -, u (5) = 0 on Γ -∪ Γ i , ∂u (5
) ∂ν = β -on Γ e-∪ Γ s-, (3.16)
In order to describe this singularity involved by the function u (5) and du to the fact that the role of the time variable is here a parameter, we begin with Proposition 3.2. Let β ∈ R and let w = w(x) ∈ H 1 (Ω + ) be solution of :

               ∆w = 0 in Ω + , w = 0 on Γ + ∪ Γ i ∪ Γ s+ , ∂w ∂ν = β on Γ e+ (3.17)
then there exists ε > 0 such that the behavior near the origin is given by

w(x) = 2a π 2 β Im[Li 2 (e -2π( x 1 L +i x 2 a ) ) -Li 2 (-e -2π( x 1 L +i x 2 a ) )] in ]0, ε[ 2 Proof.
For sake of simplicity we consider another system which leads to the same singularity (the problem is a local one). Let us write Ω a = R + * ×]0, a[ and consider the solution v ∈ H 1 (Ω a ) of the system

                   ∆v = 0 in R + * ×]0, a[, v = 0 on R + * × {0, a}, ∂v ∂ν = β on {0}×]0, a[. (3.18) We write for n ∈ N * , w n (x 2 ) = 2 a sin( nπx 2 a ) and v(x 1 , x 2 ) = n≥1 a n w n (x 2 )e - nπx 1 a
the serie being convergent in the spaces H 1 (Ω a ) and H 2 (Ω a ∩ (x 2 > 0)). Let us recall that the functions w n (n > 0) represent an orthonormal basis of L 2 (]0, a[). We obtain for p ∈ N * ,

a 0 v(x 1 , x 2 )w p (x 2 )dx 2 = a p e - pπx 1 a and ∂v ∂x 1 (x 1 , x 2 ) = - π a n a n nw n (x 2 )e - nπx 1 a thus : a 0 ∂v ∂x 1 (0, x 2 )w p (x 2 )dx 2 = - pπ a a p
which leads to :

a p = aβ pπ a 0 w p (x 2 )dx 2 = a 2 β p 2 π 2 2 a [1 -(-1) p ].
We then get Let us recall that the dilogarithm function Li 2 is defined by Li 2 (z) = n≥1 z n n 2 ( |z| ≤ 1) and we refer to [START_REF] Zagier | The Dilogarithm Function[END_REF] for very surprising properties of this function. We proved that v(x Theorem 3.1 is proved. Let us remark that in case of a unique velocity (c + = c -), there is no singular part.

Conclusion -

In this paper, we have discussed the singularity which appears when one uses a transparent boundary condition for a bimaterial. This situation occurs for instance when one tries to simulate wave propagations in an infinite strip replaced by a finite one. The singularity which belongs to the Dilog family is an artefact which doesn't exist in the physical model. Therefore it seems necessary to eliminate it from the solution. The first step was obviously to make it explicit in order to be able to suggest a method which would improve the boundary condition. For instance, a numerical method could be to compute the contribution of this artificial singularity and to substract it from the global solution. In fact, this is a way for defining an upgrade transparent boundary condition for this kind of problems. This will be discussed in a forthcoming paper. Finally, let us remark that these phenomena appear in the study of the detection of cracks at the interface between two materials, as it is noticed in [START_REF] Ph | Few remarks on the use of Love waves in non destructive testing[END_REF].

Fig. 1 . 3 Fig. 2 .Fig. 4 .

 1324 Fig. 1. The right hand side's support

Polylog singularities 5 Fig. 5 .

 55 Fig. 5. The open set Ω and the notations for the theoretical discussion

Theorem 2 . 1 .

 21 Let us define by O+, respectively O -, two open subsets whose closures are subsets of Ω + and Ω -. If the initial conditions are such that

Fig. 6 .

 6 Fig. 6. Graphs of the function S and its partial derivatives : ∂S ∂x 1 and ∂S ∂x 2

  Numerical illustrations -On figures 7 and 8, the velocities satisfy c + = 2 and c -= 0.5. The graph of ∂u ∂x 2 is plotted on Γ e in Fig. 7(a), on an interior line x = 0.5 in Fig. 7(b) whereas figures 8(a) and 8(b) point out the singularity on the outgoing side at the interface. One can see that the function ∂u ∂x 2 presents singularities on the boundaries Γ e and Γ s but not on the interior line of the domain Ω.. On figures 9 (a), (b) and and 10 (a) and (b), the graph of ∂u ∂x 2 is given for same data and time than for figures 7 and 8 but in the case where c + = c -. One can obviously see that there is no singularity in this homogeneous case.

Fig. 7 . 5 Fig. 8 .

 758 Fig. 7. (a) ∂u ∂x 2 on Γe (b) ∂u ∂x 2 on x 1 = 0.5

Fig. 9 . 5 Fig. 10 .

 9510 Fig. 9. (a) ∂u ∂x 2 on Γe (b) ∂u ∂x 2 on x 1 = 0.5

  Arg takes its values in the open set ] -π, +π[.

	December 8, 2014 20:28	Singularities
														Polylog singularities 21
	∂v ∂x 1	= -	4 π	β	p≥0	1 (2p + 1)	sin(	(2p + 1)π a	x 2 )e	-	(2p + 1)πx 1 a
	and
	∂v ∂x 2	=	4 π	β	p≥0	1 (2p + 1)	cos(	(2p + 1)π a	x 2 )e	-	(2p + 1)πx 1 a
	thus (notice that	1 + z 1 -z	∈ R -)
	∂v ∂x 2	-i	∂v ∂x 1	=	4 π	β	p≥0	z 2p+1 (2p + 1)	=	2 π	β log(	1 + z 1 -z	)
	where the function One deduces that								=	2 π	β [Ln|	1 + z 1 -z	| + iArg(	1 + z 1 -z	)]
		∂v ∂x 2		=			1 π	βLn|	1 + 2e 1 -2e	--	πx 1 a cos( πx 1 a cos(	πx 2 a πx 2 a	) + e ) + e	-2 -2	πx 1 a πx 1 a	|
	On Γ i , one has
	and	∂v ∂x 2		=			2 π	βLn|	1 + e 1 -e	--	πx 1 a πx 1 a	| =	4 π	β Ln| coth	πx 1 a	|,
													lim x1→0 +	∂v ∂x 2	= +∞.
	In a same way, one has on Γ e+ ,
												lim x2→0 +	∂v ∂x 2	(0, x

1 , x 2 ) = 2a π 2 βIm[Li 2 (z) -Li 2 (-z)].

Let us compute the gradient of v. One obtains with (3.19) : 2 ) = +∞.