Benoît Grébert

Tiphaine Jézéquel

Laurent Thomann

STABILITY OF LARGE PERIODIC SOLUTIONS OF KLEIN-GORDON NEAR A HOMOCLINIC ORBIT by

Keywords: 2000 Mathematics Subject Classification. -37K45, 35Q55, 35Bxx Klein-Gordon equation, wave equation, homoclinic orbit

We consider the Klein-Gordon equation (KG) on a Riemannian surface M

, which is globally well-posed in the energy space. Viewed as a first order Hamiltonian system in the variables (u, v ≡ ∂tu), the associated flow lets invariant the two dimensional space of (u, v) independent of x. It turns out that in this invariant space, there is a homoclinic orbit to the origin, and a family of periodic solutions inside the loops of the homoclinic orbit. In this paper we study the stability of these periodic orbits under the (KG) flow, i.e. when turning on the nonlinear interaction with the non stationary modes. By a shadowing method, we prove that around the periodic orbits, solutions stay close to them during a time of order (ln η) 2 , where η is the distance between the periodic orbit considered and the homoclinic orbit.

Introduction, statement of the main results

1.1. General introduction. -Denote by M a compact Riemannian manifold without boundary of dimension 1,2 or 3 and denote by ∆ = ∆ M the Laplace-Beltrami operator on M . In this paper we study the following nonlinear Klein-Gordon (KG) equation

(KG) ∂ 2 t u -∆u -m 2 u + u 2p+1 = 0, (t, x) ∈ R × M, u(0, x) = u 0 (x), ∂ t u(0, x) = u 1 (x),
where p ≥ 1 is an integer, and (u 0 , u 1 ) ∈ H 1 (M ) × L 2 (M ) are real-valued.

When M is the torus, this is the usual Klein-Gordon equation with periodic boundary conditions. For a more general manifold M , equation (KG) can model a vibrating membrane, for example the earth during an earthquake.

Let us recall that there exists a Hilbert basis of L 2 (M ) composed with eigenfunctions (e n ) n≥0 of ∆. Moreover, there exists a sequence 0 = λ 0 < λ 1 ≤ • • • ≤ λ n ≤ . . . so that -∆e n = λ 2 n e n , n ≥ 0. We make following assumptions Assumption 1. -The parameter m satisfies 0 < m < λ 1 .

Assumption 2. -The manifold M and the integer p satisfy either:

• M is any compact manifold without boundary of dimension 1 or 2 and p ≥ 1

• M is any compact manifold without boundary of dimension 3 and p = 1.

Moreover, up to a rescaling, we can assume that Vol M = 1.

Let us recall the well posedness result proved in [START_REF] Grébert | Dynamics of Klein-Gordon on a compact surface near a homoclinic orbit[END_REF] Proposition 1.1 ([3], Theorem 2.2). -Under Assumptions (1) and (2), the equation (KG) is globally well posed in the energy space H 1 (M ) × L 2 (M ).

Under Assumptions 1 and 2, the stationary solutions of (KG) (solutions which only depend on the space variable) are exactly the constants u = 0, u = m 1/p and u = -m 1/p . The origin is an equilibrium with an unstable direction. In fact, the eigenvalues of -∆-m 2 are the (λ 2 k -m 2 ) k∈N . Since 0 < m < λ 1 , the case k = 0 only, gives the hyperbolic directions, corresponding to the solution exp(mt) for t > 0 (resp. exp(-mt) for t < 0). It turns out that (KG) admits a homoclinic orbit to the origin which is independent of x. In the previous work [START_REF] Grébert | Dynamics of Klein-Gordon on a compact surface near a homoclinic orbit[END_REF] we have proved that (KG) admits a family of heteroclinic connections to the center manifold which are close to this homoclinic orbit. The expected picture is that we have a tube of heteroclinic connections surrounding the homoclinic orbit (but the statement in [START_REF] Grébert | Dynamics of Klein-Gordon on a compact surface near a homoclinic orbit[END_REF] is not so precise, only a large family of heteroclinic orbits is constructed).

In many physical applications, the standard Klein-Gordon equation comes with a positive mass, namely

(1.1) ∂ 2 t v -∆v + m 2 v + f (v) = 0, (t, x) ∈ R × M,
in which case the linear part of the equation has no hyperbolic direction. In fact if u is a solution of (KG) then v = u -m 1/p is a solution of (1.1) for some explicit nonlinearity f . So, when we consider (KG) instead of (1.1), we center the phase space at an unstable equilibrium (u = 0) instead of an elliptic one (v = 0). This is detailed in [3, Section 1.2].

As we pointed out in [START_REF] Grébert | Dynamics of Klein-Gordon on a compact surface near a homoclinic orbit[END_REF], the dynamics around the elliptic points u = ±m 1/p can be partially described by the KAM theory or the Birkhoff normal form theory. The KAM theory gives the existence of a large family of finite dimensional invariant tori close to the equilibrium (see [START_REF] Wayne | Periodic and quasi-periodic solutions of nonlinear wave equations via KAM theory[END_REF][START_REF] Pöschel | Quasi-periodic solutions for a nonlinear wave equation[END_REF][START_REF] Kuksin | Nearly integrable infinite-dimensional Hamiltonian systems[END_REF] or the book [START_REF] Kuksin | Nearly integrable infinite-dimensional Hamiltonian systems[END_REF]). The Birkhoff normal form approach gives the stability during polynomial times for any initial condition close to the equilibrium (see [START_REF] Bambusi | Birkhoff normal form for PDEs with tame modulus[END_REF]).

In this work our point of view is different: (KG) also admits large periodic orbits inside the loops of the homoclinic orbit (see Figure 1) and we are interested in the stability of these large periodic orbits. Observe that they are not close to the origin since they turn around u = m 1/p . By a shadowing method, we prove that around the periodic orbits, solutions stay close to them during a time of order (ln η) 2 , where η is the distance between the periodic orbit considered and the homoclinic orbit. Actually, thanks to an energy method, it is easy to get a control for times of order ln η, which is the typical timescale in the presence of a hyperbolic point. In our context, ln η is the timescale needed to achieve one loop, and our contribution consists in proving that we can follow the solution for ln η loops. The two main ingredients used in the proof are

• The Hamiltonian is negative (1) on the trajectory, which implies a confinement of the solution;

• The trajectory is close to the homoclinic orbit, which gives the pattern of the solution. Combining these two facts, we conclude with a bootstrap argument.

1.2. Hamiltonian structure of (KG). -As in [START_REF] Grébert | Dynamics of Klein-Gordon on a compact surface near a homoclinic orbit[END_REF], we define the scalar product on L 2 (M ) by

f, g = 1 Vol M M f g
, where Vol M denotes the volume of M , we assume that e n L 2 = 1 and we set e 0 = 1.

Denote by

v = ∂ t u and introduce (1.2) H = 1 2 M |∇ x u| 2 + v 2 -m 2 u 2 + 1 2p + 2 M u 2p+2 .
Then, the system (KG) is equivalent to

(1.3) u = δH δv , v = - δH δu .
We write

u(t, x) = +∞ n=0 a n (t)e n (x), v(t, x) = +∞ n=0 b n (t)e n (x),
where

(a n ) n∈N ∈ h 1 (N, R) := x = (x n ) n∈N | x 2 h 1 = n∈N (1 + λ 2 n )|x n | 2 < +∞ , (b n ) n∈N ∈ 2 (N, R) := x = (x n ) n∈N | x 2 2 = n∈N |x n | 2 < +∞ ,
in such a way that to the continuous phase space

X := H 1 × L 2 corresponds the discrete one h 1 × 2 .
We endow this space with the natural norm and distance (1.4)

X X = u h 1 + v 2 , for X = (u, v) dist X (X, Y ) = X -Y X , for X, Y ∈ X .
In the coordinates (a n , b n ) n≥0 , the Hamiltonian in (1.2) reads

(1.5) H = 1 2 +∞ n=0 (λ 2 n -m 2 )a 2 n + b 2 n + 1 2p + 2 M +∞ k=0 a k e k (x) 2p+2 dx, 1.
Our Hamiltonian H is defined so that H(0) = 0. and the system (1.3) becomes

(1.6)      ȧn = b n , n ≥ 0 ḃn = -(λ 2 n -m 2 )a n - M +∞ k=0
a k e k (x) 2p+1 e n (x)dx, n ≥ 0.

1.3. Space-stationary solutions: homoclinic orbit and family of periodic orbits. -The space-stationary solutions of (KG) exactly correspond to the solutions of (1.6) satisfying a n = b n = 0 for n ≥ 1. In this case, the equation on (a 0 , b 0 ) reads

(1.7) ȧ0 = b 0 ḃ0 = m 2 a 0 -a 2p+1 0 ,
and this system possesses a homoclinic solution to 0, which we will denote in the sequel by h : t → (α(t), β(t)). We denote by K 0 ⊂ h 1 × 2 the curve which is described by a 0 (t) = h(t), b 0 (t) = h (t) and a n (t) = b n (t) = 0 for n ≥ 1 (see Figure 1). Indeed we can explicitly compute (1.8)

b 0 a 0 η K 0 K η m 1 p
h(t) = m 1/p (p + 1) 1/(2p) cosh(pmt) 1/p , h (t) = -m 1/p+1 (p + 1) 1/(2p) sinh(pmt) cosh(pmt) 1/p+1 .
For η > 0 denote by K η the trajectory of (1.6) given by the initial conditions a 0 (0) = η, b 0 (0) = 0, a n (0) = b n (0) = 0 for n ≥ 1 (see Figure 1). This orbit remains in the plane

{(a, b) ∈ h 1 × 2 | a n = b n = 0 for n ≥ 1}
and is periodic with a period of order ln 1 η (see below for the proof): most of that time is dedicated to cover the very principle and the very end of the loop, i.e. when (a 0 , b 0 ) is close to (η, 0).

Our main result states that for η small enough these periodic solutions are stable during a long time of order (ln η) 2 in the sense that if one starts η 3 -close to K η then one remains η 2 -close to K η during a time of order (ln η) 2 . More precisely, recall the definition (1.4) of dist X , then Theorem 1.2. -Let M , p and m fixed and satisfying Assumptions 1 and 2. There exists η 0 > 0, 0 < c < 1 and C > 0 such that if 0 < η < η 0 and if the initial datum (a(0

), b(0)) satisfies dist X (a(0), b(0)), K η ≤ η 3
then the solution of (1.6) satisfies (i) There exists

T η > 0 satisfying c ln 1 η ≤ T η ≤ 1 c ln 1 η such that dist X (a(T η ), b(T η )), (a(0), b(0)) ≤ Cη 2 (ii) For all |t| ≤ c ln 1 η 2 dist X a(t), b(t) , K η ≤ Cη 2 .
This result shows that a trajectory issued from a tube of thickness η 3 centered on K η describes in the mode 0 several loops around (m 1 p , 0) and remains very small in the other modes. Indeed the time T η is the time necessary for (a 0 , b 0 ) to complete one loop around (m 1 p , 0), and the previous result shows that the trajectory completes at least c 2 ln 1 η loops around (m 1 p , 0), staying inside the tube of thickness η 2 centered on K η .

Remark 1.3. -Numerical simulations confirm this stability result and seem to show that it holds true for larger timescales. Moreover, it is likely that there exist initial conditions as in Theorem 1.2 and which satisfy (ii) for all times. Such a result is out of reach with our method and could possibly be attacked with a KAM method.

1.4. On the linear case.-This stability problem is related to the more general question of existence of periodic or quasi-periodic solutions to (KG). This question is not evident even in the linear case. Actually, we can consider the linearized equation of (KG) around the periodic orbit K η as a first approximation to understand the possible dynamics: take a 0 (t) as a periodic solution of the system (1.7) inside the loops (see Figure 1) and search for w = k≥1 a k e k as a (small) solution of the following linear wave equation

(1.9) ∂ 2 t w -∆w -m 2 w + V (t)w = 0, (t, x) ∈ R × M where t → V (t) = (2p + 1)a 2p 0 (t
) is a time depending potential. In the discrete variables, the equation (1.9) reads

(1.10) ȧn = b n ḃn = -(λ 2 n -m 2 )a n -V (t)a n for n ≥ 1.
The existence of small periodic (or quasi-periodic) solutions of (1.9) is related to the reducibility of this equation to the autonomous case and to the Floquet theory. In the case where V is a small potential the question can be considered with the help of the KAM theory (see [START_REF] Eliasson | On reducibility of Schrödinger equations with quasiperiodic in time potentials[END_REF] in the case of a linear Schrödinger equation on a torus or [START_REF] Grébert | KAM for the Quantum Harmonic Oscillator[END_REF] in the case of a linear Schrödinger equation on the line with a harmonic potential). To the best of our knowledge, there is no result in the case of a time dependent linear perturbation which is not small. With the same techniques, we can apply our result to the linear equation (1.7),(1.10), and this shows that there is at least some stability around the solutions K η .

Remark 1.4. -In our previous paper [START_REF] Grébert | Dynamics of Klein-Gordon on a compact surface near a homoclinic orbit[END_REF], we also consider the linearised system (1.10) but in the case where t → a 0 (t) describes the homoclinic orbit, and in which case the linear flow is bounded. We refer to [START_REF] Grébert | Dynamics of Klein-Gordon on a compact surface near a homoclinic orbit[END_REF]Appendix] where we crucially use that +∞ 0 a 0 (t)dt < +∞.

1.5. Comments on the proof. -The method used in the proof of Theorem 1.2 relies on a bootstrap (or continuation) argument, is rather elementary and can be used in various situations. It can work when one tries to control perturbations of a one degree of freedom system, even when this perturbation has an infinite number of directions as in our case. The underlying idea is to give a precise description of one particular mode (here (a 0 , b 0 )) while the others are shown to remain small during the considered timescale. With these ingredients, we typically obtain a control of the solution during times of order ln η. In order to improve the control up to order (ln η) 2 we crucially use the fact that our periodic solution K η is close to the homoclinic orbit. This enables us to define a modified functional I(t) which controls the Sobolev norm of the exterior modes and which is almost constant (see Lemma 2.5).

Notations. -In this paper c, C > 0 denote constants the value of which may change from line to line. These constants will always be universal, or depend on the fixed quantities m and p. We denote by N the set of the non negative integers, and a k e k , so that in particular, u reads u = a 0 e 0 + U .

N * = N\{0}. We set X = H 1 (M ) × L 2 (M ).
In the sequel we will use the following decomposition of the energy H (see (1.5))

(2.2)

H = 1 2 (b 2 0 -m 2 a 2 0 ) + J + 1 2(p + 1) M a 0 + U 2p+2 dx = 1 2 b 2 0 + 1 2 f (a 0 ) + J + r, where (2.3) f (x) := -m 2 x 2 + 1 p + 1 x 2p+2 , (2.4) J := 1 2 +∞ k=1 (λ 2 k -m 2 )a 2 k + b 2 k ,

and

(2.5)

r(t) = 1 2(p + 1) M (a 0 + U ) 2p+2 -a 2p+2 0 dx.
Note that J ∼ U 2 X . In Subsection 2.4, we prove that J remains small. As a consequence r is also small since it is a quadratic quantity in U . Then, as long as J is small, the dynamics of the solution is governed by the space stationary dynamics (see Subsection 1.3) given by the Hamiltonian

H 0 = 1 2 b 2 0 + f (a 0 ).
Moreover, we introduce the notations q k for k ∈ N, defined by (2.6)

q k = M (a 0 + U ) 2p+1 -a 2p+1 0 -(2p + 1)a 2p 0 U e k dx,
so that the systems (1.6) read

(2.7) ȧ0 = b 0 ḃ0 = -f (a 0 ) -q 0 , (2.8) ȧk = b k , k ≥ 1 ḃk = -(λ 2 k -m 2 )a k -(2p + 1)a 2p 0 a k -q k , k ≥ 1.
We finally fix a real δ ∈]0, m 1 p [, and denote by δ the unique real of ]m 1 p , (p + 1) • When M has dimension 1 or 2, for all 2 ≤ q < +∞ U L q (M ) ≤ C q J 1/2 .

1 2p m 1 p [ such that f (δ ) = f (δ).
• When M has dimension 3, for all 2 ≤ q ≤ 6

U L q (M ) ≤ C q J 1/2 .
Proof. -By Sobolev, in each of the previous cases, there exists C q > 0 so that for all U ∈ H 1 (M ) we have U L q (M ) ≤ C q U H 1 (M ) and the result follows.

Lemma 2.2. -Let (a, b) ∈ X such that dist X (a, b), K η ≤ η 3 then (2.9) H(a, b) = - 1 2 m 2 η 2 + O(η 3 ).
Remark 2.3. -In the following, we prove in Lemma 2.6 a reciprocal result in the case J small.

Proof. -Let (ã, b) ∈ K η such that (a, b) -(ã, b) X ≤ η 3 .
On the one hand, since we have (η, 0) +∞ n=1 (0, 0) ∈ K η , then

H(ã, b) = H (η, 0) +∞ n=1 (0, 0) = - 1 2 m 2 η 2 .
On the other hand denoting by ũ =

+∞ k=0 ãk e k and J = 1 2 +∞ k=1 (λ 2 k -m 2 )ã 2 k + b2 k H(a, b) -H(ã, b) ≤ ≤ 1 2 | b2 -b 2 0 | + m 2 2 |ã 2 0 -a 2 0 | + |J -J| + 1 2(p + 1) M |u 2p+2 -ũ2p+2 |dx. ≤ C (a, b) -(ã, b) X (a, b) X + (ã, b) X + C u -ũ L 2p+2 (M ) u 2p+1 L 2p+2 (M ) + ũ 2p+1 L 2p+2 (M ) ≤ C (a, b) -(ã, b) X = O(η 3 ).
2.3. Estimates for J small. -We now fix an initial condition satisfying dist X (a(0), b(0)), K η ≤ η 3 with 0 < η 1. In view of Lemma 2.2 we can choose η small enough in such a way that (2.10) H(a(0), b(0)) < 0.

By definition of K η we have a 0 (0) > 0 and without loss of generality we can also assume that b 0 (0) ≥ 0.

Lemma 2.4. -For η small enough, as long as J(t) ≤ η 5 , there exists

C 1 = C 1 (m) > 0 such that (2.11) η/2 ≤ a 0 (t) ≤ C 1 , |b 0 (t)| ≤ C 1 .
Proof. -Assume that J ≤ η 5 . First we prove that a 0 cannot vanish: Indeed if it was the case, by (2.2) we would have H(a(t), b(t)) ≥ 0 for some t ∈ R, which is in contradiction with (2.10) and the conservation of H. Next we prove that there exists C > 0 so that (2.12)

1 2(p + 1) M a 0 + U 2p+2 dx ≥ 1 4(p + 1) a 2p+2 0 -CJ p+1 .
Write the binomial expansion (2.13)

M a 0 + U 2p+2 dx = a 2p+2 0 + 2p+2 j=1 2p + 2 j a 2p+2-j 0 M U j dx,
Apply Lemma 2.1 and use the Young inequality

c 1 c 2 = (εc 1 )(ε -1 c 2 ) ≤ ε q q c q 1 + 1 rε r c r 2 , c 1 , c 2 ≥ 0, ε > 0, 1 q + 1 r = 1,
to each term of the sum with ε > 0 small enough, and get

2p+2 j=1 2p + 2 j a 2p+2-j 0 M U j dx ≤ C 2p+2 j=1 a 2p+2-j 0 J j/2 ≤ 1 2 a 2p+2 0 + CJ p+1 ,
which together with (2.13) yields (2.12). Now, for η > 0 small enough, CJ p+1 ≤ J so by (2.12) and Lemma 2.2 we obtain

(2.14) 1 2 (b 2 0 -m 2 a 2 0 ) + 1 4(p + 1) a 2p+2 0 ≤ H = - 1 2 m 2 η 2 + O(η 3 ) ≤ - 1 4 m 2 η 2 ,
from which we deduce that |b 0 | ≤ C 1 and 0 < a 0 ≤ C 1 for some constant C 1 depending on m. Furthermore, the bound (2.14) implies -m 2 a 2 0 /2 ≤ H ≤ -1 4 m 2 η 2 which completes the proof. In the following we use the additional notations I and r:

I := 1 2 +∞ k=1 λ 2 k -m 2 + (2p + 1)a 2p 0 a 2 k + b 2 k , r(t) = 1 2(p + 1) M (a 0 + U ) 2p+2 -a 2p+2 0 -(2p + 1)a 2p 0 U 2 dx. So that H reads H = 1 2 b 2 0 + 1 2 f (a 0 ) + I + r.
Observe that as long as J ≤ Cη 5 , since a 0 is bounded (see (2.11)), there exists K 0 > 0 so that for all

(a k , b k ) k≥1 ∈ h 1 × 2 , (2.15) J ≤ I ≤ K 0 J. Lemma 2.5.
-As long as J(t) ≤ η 5 , we have the a priori bound

(2.16) | İ| ≤ Ca 2p-1 0 | ȧ0 |I + CI 3/2 .
Proof. -We recall the definition (2.6) of the q k , and that then the equations on

(a k , b k ) for k ≥ 1 read ȧk = b k , ḃk = -λ 2 k -m 2 + (2p + 1)a 2p 0 a k -q k . With this system, we compute (2.17) | İ| = p(2p + 1)a 2p-1 0 ȧ0 +∞ k=1 a 2 k - +∞ k=1 b k q k ≤ p(2p + 1)a 2p-1 0 | ȧ0 | +∞ k=1 a 2 k + b 2 (N * ) q 2 (N * ) .
We have b 2 (N * ) ≤ CI 1/2 . Then by Parseval, (2.15) and the fact that

|a 0 | ≤ C q 2 (N * ) = (a 0 + U ) 2p+1 -a 2p+1 0 -(2p + 1)a 2p 0 U L 2 (M ) ≤ C U 2 L 2 (M ) + C U 2p+1 L 2 (M ) = C U 2 L 4 (M ) + C U 2p+1 L 2(2p+1) (M ) ≤ CI + CI p+ 1 2 ≤ CI, (2.18)
where in the last line we used Lemma 2.1. Then (2.18) together with (2.17) gives the result.

Lemma 2.6. -If (a, b) ∈ X such that H(a, b) = -1 2 m 2 η 2 + O(η 3 ), η ≤ a 0 and J ≤ η 3 , then dist X (a 0 , b 0 , 0, • • • ), K η ≤ Cη 2 .
More precisely, (a 0 , b 0 , 0,

• • • ) is η 2 -close to (ã 0 , b0 , 0, • • • ) ∈ K η defined by (ã 0 , b0 ) = (a 0 , sgn(b 0 ) f (η) -f (a 0 )) if a 0 ∈ [δ, δ ], (f -1 f (η) -b 2 0 , b 0 ) if a 0 / ∈ [δ, δ ].
Proof. -From the expression (2.2) of the energy, the assumption reads

b 2 0 + f (a 0 ) + J + r = -m 2 η 2 + O(η 3
). On the one hand, observe that if J ≤ η 3 , then r ≤ cη 3 : the proof is similar to the one of (2.18). On the other hand, f (η) = -m 2 η 2 + O(η 3 ). Thus we have

(2.19) b 2 0 + f (a 0 ) = f (η) + O(η 3 ). Recall that (ã 0 , b0 , 0, • • • ) belongs to K η given that (2.20) b2 0 + f (ã 0 ) = f (η). • Assume that a 0 / ∈ [δ, δ ]. From (2.19) we deduce a 0 = f -1 f (η) -b 2 0 + O(η 3 ) .
Since a 0 ≥ η, thus |f (a 0 )| ≥ c 0 η and we can apply the mean value theorem to get

|ã 0 -a 0 | ≤ Cη 3 sup |(f -1 ) | ≤ Cη 2 . • In the region a 0 ∈ [δ, δ ],

First loop : proof of Theorem 1.2 (i).

-This part is devoted to the proof of (i) of Theorem 1.2. More precisely, we prove a stronger version of it in Lemma 2.10, which will then be useful to prove (ii) of the theorem (in the next part).

Recall that in all this section we have fixed an initial condition such that

(2.21) dist X (a(0), b(0)), K η ≤ η 3
and 0 < η 1 such that (2.10) is satisfied.

Lemma 2.7. -There exists K > 1 independent of η > 0 such that: as long as b 0 (t) ≥ 0 and |t| ≤ η -2 , J satisfies J(t) ≤ KJ(0) ≤ η 5 .

Proof. -We prove it by a bootstrap argument. As long as J(t) ≤ η 5 and b 0 (t) ≥ 0, the estimate (2.16) holds true and by Gronwall we get

J(t) ≤ I(t) ≤ I(0)e |t| 0 ( ȧ0 (s)a 2p-1 0 (s)+η 5/2 )ds ≤ I(0)e (a 0 (t) 2p -a 0 (0) 2p +|t|η 5/2 ) ≤ I(0) exp 2C 2p 1 (m) + 1 , where C 1 (m) > 0 is defined in (2.4).
Recall that by (2.11), a 0 is bounded as long as J(t) ≤ η 5 , and that I(0) ≤ K 0 η 6 by (2.21) (recall the definition of K 0 in (2.15)). Therefore, under the assumptions that b 0 (t) ≥ 0 and |t| ≤ η -2 we get that as long as J(t) ≤ η 5 , we have J(t) ≤ K 1 I(0), which gives the result for η sufficiently small, since I(0) ≤ K 0 J(0). Lemma 2.8. -Let α 0 be a real with 0 < α 0 ≤ δ. We suppose that the initial conditions (a(0), b(0)) satisfy 0 < a 0 (0) < α 0 and b 0 (0) ≥ 0. Then there exist times τ (α 0 ) and τ η such that a 0 (τ (α 0 )) = α 0 , a 0 (τ η ) = δ and b 0 (t) > 0 for t ∈]0, τ η ].

Moreover, we have the estimates

τ (α 0 ) ≤ τ η ≤ C ln 1 η, and if b 0 (0) = 0, c ln 1 η ≤ τ η .
Proof. -• Firstly, we can go back to the particular case b 0 (0) = 0 : for that, we prove that there exists a time t ∈]0, η -2 [ such that b 0 (-t) = 0. We proceed by contradiction : suppose that b 0 (-t) > 0 for all t ∈ [0, 4b 0 (0) m 2 η ]. We know from Lemma 2.7, that for all those t, J(-t) ≤ η 5 . Thus q 0 (-t) ≤ Cη 5 , and by (2.11), we have that a 0 (-t) > η/2. Hence,

ḃ0 (-t) = -f (a 0 (-t)) -q 0 (t) > -f ( η 2 ) + O(η 5 ) = m 2 η 2 + O(η 3 ) > m 2 η 3 .
Then for all t in [0, 4b 0 (0)

m 2 η ], b 0 (-t) < b 0 (0) -tm 2 η 3 holds.
In particular, we obtain b 0 (-4b 0 (0) m 2 η ) < 0 : the contradiction. • So for all initial conditions satisfying the assumption of the lemma, there exists a negative time -t such that b 0 (-t) = 0 and J(-t) ≤ η 5 .

In the following we only consider the case a 0 (0) = η and b 0 (0) = 0.

Observe that the upper bound we will get for τ η in this case will be an upper bound in the general case.

• Let ϕ be the function defined by

η 5 ϕ(t) = -2(J(t) -J(0)) -2(r(t) -r(0)).
From the preservation of energy, see (2.2) and (2.9), we get that ϕ satisfies

(2.22) b 2 0 + f (a 0 ) = f (η) + η 5 ϕ(t).
Given that b 0 (0) = 0 and ḃ0 (0) > 0 for η sufficiently small, b 0 (t) > 0 for some times, and as long as

b 0 (t) ≥ 0, b 0 (t) reads (2.23) ȧ0 (t) = b 0 (t) = f (η) -f (a 0 (t)) + η 5 ϕ(t).
In particular, as long as b 0 (t) > 0, by integration and the change of variables α = a 0 (t) we obtain (2.24)

t = t 0 ȧ0 (s)ds f (η) -f (a 0 (s)) + η 5 ϕ(s) = a 0 (t) η dα f (η) -f (α) + η 5 ϕ(a -1 0 (α))
.

• Let us prove that, as long as J(t) ≤ η 5 • We now show that there exist c, C > 0 so that as long J(t) ≤ η 5 , a 0 (t) ≤ δ and b 0 (t) ≥ 0,

(2.26) c(a 2 0 -η 2 ) ≤ f (η) -f (a 0 ) + η 5 ϕ(t) ≤ C(a 2 0 -η 2 ). Observe that for all η ≤ α ≤ δ, (m 2 -δ 2p )(α 2 -η 2 ) ≤ f (η) -f (α) ≤ m 2 (α 2 -η 2 )
holds, and recall that δ < m 1 p . It is then enough to show that η 5 ϕ(t) is a negligible amount when η small. Recall the form (2.7) of the system :

ä0 = ḃ0 = m 2 a 0 -a 2p+1 0 -q 0 = -f (a 0 ) -q 0 (t).
On one hand, with the same calculations as in the proof of Lemma 2.5, we obtain that as long as J(t) ≤ η 5 , we have |q 0 (t)| ≤ η 5 . On the other hand, by a study of the variations of f , we get that for η ≤ α ≤ δ, if η small enough then f (α) ≤ -cη. We deduce that as long as J(t) ≤ η 5 and a 0 ≤ δ, we have ä0 ≥ cη and by Taylor, |a 0 (t) -a 0 (0)| ≥ cηt 2 . This in turn implies that when moreover a -1 0 is well-defined (for b 0 (t) > 0)

|a -1 0 (α)| = |a -1 0 (α) -a -1 0 (η)| ≤ C η 1/2 |α -η| 1/2
holds. Then by (2.25), we get that η 5 ϕ(t) ≤ Cη 5 η -1 |a 0 -η| ≤ Cη 3 (a 2 0 -η 2 ), which proves (2.26).

• Firstly, from (2.26) with (2.23), we get that as long as J(t) ≤ η 5 and a 0 (t) ≤ δ 2 , we have b 0 (t) > 0. Secondly, with (2.26) together with (2.24), we can compute an estimate of t in term of a 0 (t). We will see then that a 0 (t := τ η ) = δ holds for a time τ η smaller than η -2 , so that by Lemma 2.7 we are still in the regime J(t) ≤ η 5 . Precisely, observe that

α 0 η dα α 2 -η 2 = ln   α 0 η + α 2 0 η 2 -1   .
Then (2.26) together with (2.24) gives the estimate claimed in the lemma, and in particular τ η ≤ η -2 for η sufficiently small.

Consider the time τ η given in Lemma 2.8. We now show that b 0 admits zeros on ]0, +∞[ and that the smallest one is close to τ η . Lemma 2.9. -For any initial condition (a(0), b(0)) such that b 0 (0) ≥ 0, there exists a time T 1 η such that b 0 (T 1 η ) = 0, and b 0 (t) > 0 for all t ∈]0, T 1 η [. Moreover, if a 0 (0) < δ, then there exist C such that 0 < T 1 η -τ η ≤ C, where τ η is given by Lemma 2.8 and if a 0 (0) ≥ δ, then

T 1 η ≤ C. Proof. -Recall that δ ∈]0, m 1
p [, and δ is the unique real of ]m 1 p , (p + 1)

1 2p m 1 p [ such that f (δ ) = f (δ).
In the proof we consider the case a 0 (0) < δ. To prove the lemma, we proceed in two steps, here is a summary of them :

(i) We show that there exists a time

τ 1 η > τ η such that b 0 (τ 1 η ) = b 0 (τ η ), a 0 (τ 1 η ) -δ ≤ Cη 5 and b 0 (t) > b 0 (τ η ) for t ∈]τ η , τ 1 η [. Moreover, τ 1 η -τ η ≤ C holds. (ii) Then we show the existence of T 1 η > τ 1 η such that b 0 (T 1 η ) = 0, |a 0 (T 1 η )-h(0)| ≤ Cη 2 , and b 0 (t) > 0 for all t ∈]τ 1 η , T 1 η [. Moreover, T 1 η -τ 1 η ≤ C holds. If a 0 (0) ∈]δ, δ [,
then in step (i) the proof is the same by considering a 0 (0) instead of a 0 (τ η ). If a 0 (0) > δ , then step (ii) is sufficient, considering here a 0 (0) instead of a 0 (τ 1 η ).

Step (i). From Lemma 2.8 we know that a 0 (τ η ) = δ. We then obtain that b 0 (τ η ), b 0 (τ η ) ≥ c 0 > 0 (c 0 depends only on δ which is fixed). In particular, b 0 (t) ≥ b 0 (τ η ) for small times t ≥ τ η . Let us prove this step, by contradiction : suppose that for all t in ]τ η , τ η + (h(0) -δ)/c 0 ], we have b 0 (t) > b 0 (τ η ). Then for all those t, a 0 (t) > δ + c 0 • (t -τ η ). Then we have a 0 (τ η + (h(0) -δ)/c 0 ) > h(0), and thus there exists t ∈]τ η , τ η + (h(0) -δ)/c 0 ] such that f (a 0 (t)) = f (h(0)) = 0. Hence, for this t, the assumption b 0 (t) > b 0 (τ η ) is in contradiction with the energy preservation (2.22) for η sufficiently small, given that t is in the regime |ϕ(t)| ≤ C. Finally, there exists τ 1 η such that 0 < τ 1 η -τ η ≤ C and b 0 (τ 1 η ) = b 0 (τ η ). To get the approximation of a 0 (τ 1 η ), on the one hand observe that necessarily there exists a time t ∈]τ η , τ 1 η [ such that b 0 (t) = 0. From (2.7) we thus have |f (a 0 (t))| ≤ Cη 5 , which means that |a 0 (t)m 1 p | ≤ C η 5 . On the other hand, from the preservation of the energy (2.22), given that we are still in the regime |ϕ(t)| ≤ C, we obtain that |f (a 0 (τ 1 η )) -f (δ)| ≤ Cη 5 . Then, given that a 0 (τ 1 η ) > a 0 (t), we get that for η sufficiently small, |a 0 (τ 1 η ) -δ | ≤ Cη 5 .

Step (ii). Let us prove this step by contradiction : suppose that for all t in ]τ 1 η , τ 1 η + 2b 0 (τ η )/f (δ)], we have b 0 (t) > 0. Then for all those t, a 0 (t) > a 0 (τ 1 η ) and thus b 0 (t) < -f (δ ) + O(η 5 ) < 1 2 f (δ) for η small.

Then we have b 0 (τ 1 η + 2b 0 (τ η )/f (δ)) < 0, and get the contradiction. Thus there exists T 1 η such that b 0 (T 1 η ) = 0 and 0 < T 1 η -τ 1 η ≤ C. From the energy preservation (2.22) and given that a 0 (T 1 η ) > a 0 (τ 1 η ) > m Proof. -Let us explain the strategy of the proof. We follow the periodic orbit with initial (a(0), b(0)).

When we are δ-close to 0 we use the arguments of Lemma 2.8, and when a 0 ≥ δ, we use the arguments of Lemma 2.9. Notice that by Lemma 2.6 and the precise definition of T η we know that (a 0 (0), b 0 (0)) and (a 0 (T η ), b 0 (T η )) are η 2 close to each other.

2.5. Many loops : proof of Theorem 1.2 (ii). -We would like to take (a(T η ), b(T η )) as an initial condition in Lemma 2.10 to iterate the process. But by Lemmas 2.10 and 2.6, we only know that dist (a 0 (T η ), b 0 (T η )), K η ≤ η 2 , J(T η ) ≤ Kη 6 . This is a priori too weak to apply Lemma 2.10. But we can choose to apply the Lemma 2.10 with η such that (a 0 (T η ), b 0 (T η )) ∈ K η , which gives the same estimates.

We see then that we can iterate the argument as long as J ≤ η 5 , namely N times (each time corresponds to a loop) with (2K) N η 6 = η 5 i.e. N = c ln 1 η .

This gives that for |t| ≤ c ln 1 η 2 , J ≤ η 5 and we can iteratively apply Lemma 2.10.

Figure 1 .

 1 Figure 1. Phase portrait for the space-stationary set a n = b n = 0 for n ≥ 1.

2 .

 2 Proof of the main result 2.1. Notations and strategy of the proof. -For u = +∞ k=0 a k e k we define (2.1) U := +∞ k=1

2. 2 .

 2 Preliminary results. -In the sequel we will need Lemma 2.1. -Let U ∈ H 1 (M ). Then

  we define g(b 0 ) = b 2 0 . Then we have |g (b 0 )| ≥ c and we can perform a similar argument to get | b0 -b 0 | ≤ Cη 3 .

( 2 .

 2 25) |ϕ(t)| ≤ C min (1, t 2 ). Since J ≤ η 5 and by (2.15) together with (2.18), we have |ϕ(t)| ≤ C. Then, since b 0 (0) = 0 and ȧ0 (0) = b 0 (0) = 0, from (2.22) we infer that φ(0) = 0. To obtain (2.25), it is therefore enough to prove that | φ| ≤ C. This is done by computing J and r with the relations (2.8), and with Cauchy-Schwarz we can check that | J| ≤ Cη 5 and |r| ≤ Cη 5 .

1 p

 1 and f (h(0)) > 0, we get the estimate on a 0 (T 1 η ).

Lemma 2 . 5 and c 0 ln 1 η 1 p 1 p

 25111 10. -There exist some constants c 0 , C 0 , C > 0 such that for any initial condition (a(0), b(0)) satisfying dist (a 0 (0), b 0 (0)), K η ) ≤ η 3 , J(0) ≤ η 5 , there exists a time T η such that(a 0 (T η ), b 0 (T η )) -(a 0 (0), b 0 (0)) ≤ Cη 2 , J(T η ) ≤ 2KJ(0) ≤ η ≤ T η ≤ C 0 ln 1 η .Precisely, we can define uniquely T η as the time of first return to the set(a, b) | a 0 = a 0 (0), sgn(b 0 ) = sgn(b 0 (0)) , if a 0 (0) ∈ [δ, δ ] or (a, b) | b 0 = b 0 (0), sgn(a 0 -m ) = sgn(a 0 (0) -m ) , if a 0 (0) / ∈ [δ, δ ].

B.G. was supported in part by the grant ANR-10-BLAN-DynPDE. B.G. and L.T. were supported in part by the grant "HANDDY" ANR-10-JCJC 0109 and by the grant "ANAÉ" ANR-13-BS01-0010-03.


Marie Piraud

Laurent Sanchez-Palencia

Bart Van Tiggelen

Anderson Localization of Matter Waves in 3D Anisotropic Disordered Potentials

Keywords: numbers: 03, 75, -b, 05, 60, Gg

Using a cut-off free formulation of the coherent transport theory, we show that the interference terms at the origin of localization strongly affect the transport anisotropy. In contrast to the common hypothesis, we then find that the anisotropies of incoherent and coherent diffusion are significantly different, in particular at criticality. There, we show that the coherent transport anisotropy is mainly determined by the properties of the disorder-averaged effective scattering medium while the incoherent transport contributions become irrelevant.

I. INTRODUCTION

The propagation of a coherent wave in a disordered medium is strongly affected by interference of the various multiple-scattering paths, which may suppress or even cancel diffusion. This phenomenon, known as Anderson localization (AL), is a widely studied problem at the frontier of condensed-matter and wave physics. Anderson localization has now been observed in a variety of systems with electromagnetic waves [1][2][3], acoustic waves [4], and ultracold matterwaves [5][6][7][8][9][10][START_REF] Semeghini | Measurement of the mobility edge for 3D Anderson localization[END_REF]. Basic knowledge of AL relies on the discrete Anderson model [START_REF] Anderson | [END_REF][13][14] and the one-parameter scaling theory [15], which in particular predicts a universal localization transition in dimension d > 2. Recent progress on control of disordered systems [16][17][18][19][20] triggered a renewed interest for a refined, microscopic understanding of AL, which cannot be described in the framework of the universal scaling theory.

An example where a microscopic theory is vital is that of localization in anisotropic media, which are relevant to MOSFETs [21], liquid crystals [22,23], phosphides [24], or ultracold atoms [8][9][10][START_REF] Semeghini | Measurement of the mobility edge for 3D Anderson localization[END_REF][25][26][27], for instance. These systems may be cast in two classes of anisotropy, either resulting from isotropic scatterers embedded in an anisotropic underlying medium (mass anisotropy) [28] or resulting from anisotropic disorder in an isotropic medium (disorder anisotropy). For these systems, basic knowledge presently relies on the extension of the onshell self-consistent theory to anisotropic media [29]. The latter predicts that the interference terms have the same anisotropy as the incoherent propagation terms, so that the anisotropies of localization and classical diffusion are equal. However, this prediction relies on the introduction of an elliptic cut-off whose anisotropy is arbitrarily chosen as the inverse of the incoherent transport mean free path.

In this paper, using a cut-off free formulation of the coherent transport theory, we show that interference terms at the origin of AL strongly affect the transport anisotropy. The solution of the transport equa-tions shows that the anisotropies of incoherent and coherent diffusion are significantly different, especially near the mobility edge (ME). This invalidates the use of the elliptic cut-off used so far. Our work indicates that the anisotropy in coherent diffusion is mainly determined by the properties of the average effective medium. Hence, contrary to the usual hypothesis, the interference terms significantly compensate for the large anisotropy predicted by incoherent Boltzmann transport theory. Our results provide insight into the theory of AL in anisotropic disorder and should significantly affect the interpretation of experimental data.

II. QUANTUM TRANSPORT

A. Formalism

The coherent transport equations for a wave in an anisotropic disordered medium are obtained from the exact Bethe-Salpeter equation [START_REF] Rammer | Quantum Transport Theory[END_REF] by generalizing the approach of Ref. [START_REF] Yedjour | [END_REF] to the anisotropic case. Briefly speaking (see Appendix A for details), the density propagation of a matterwave of energy E and wavevector k in a disordered medium is described in the long time and large distance limit by the static current vertex function j(E, k), which obeys the closed equation

j(E, k) = j 0 (E, k) + |G(E, k)| 2 dk ′ (2π) d j(E, k ′ ) U k,k ′ (E) (1) where j 0 (E, k) = [k + ∇ k ReΣ(E, k)] ImG(E, k) 2 - 1 2 ReG(E, k)ImG(E, k)∇ k ImΣ(E, k), G(E, k
) is the disorder average of the single-particle Green function, Σ(E, k) is the associated self-energy, and U k,k ′ (E) = U k,k ′ (E, ω = 0, q = 0) is the static irreducible vertex function (ω and q are the conjugate variables of time and space, respectively). The components of the diffu-sion tensor D(E) are then given by the Kubo formula

D uv (E) = m 1 πN (E) dk (2π) d k u j v (E, k), (2)
with

N (E) = dk (2π) d -ImG(E,k) π
the density of states per unit volume.

B. Interference term

In the weak localization regime, we retain the following form for the irreducible vertex function

U k,k ′ (E) = U B k,k ′ + U MC k,k ′ (E), (3)
where

U B k,k ′ = C(k -k ′ )
is the vertex associated with incoherent (Boltzmann) transport, with C(k) the power spectrum (fluctuations in the potential energy) [START_REF] Here | use C(q, ω) ≡ drdt C(r, t) exp[-i(q • r -ωt)[END_REF], and U MC k,k ′ is the leading interference contribution expressed by the maximally-crossed diagrams [START_REF] Yedjour | [END_REF],

U MC k,k ′ (E) = 2 πN (E) ImΣ(E, k-k ′ 2 ) 2 (k + k ′ ) • D(E) • (k + k ′ ) . (4)
Note that the diffusion tensor D(E) is included selfconsistently in Eqs. ( 1) and ( 2), an approximation which gives quantitative estimates in good agreement with numerical calculations [START_REF] Kroha | [END_REF]. In standard (on-shell) selfconsistent theory, the pole of Eq. ( 4) necessitates the introduction of an ultraviolet cut-off in the integral in the right-hand side of Eq. ( 1). Here, by keeping the full k dependence of the vertex function j(E, k), this artificial divergence does not appear [START_REF] Yedjour | [END_REF].

C. Disorder correlation function

For simplicity, we assume that the disorder V (r) is Gaussian-distributed. Its statistical properties are fully characterized by the two-point correlation function C(r) = V (0)V (r), where we have set the zero of energies such that V = 0. In order to carry out a number of integrations analytically, it is convenient to use a correlation function with uniaxial symmetry around, say, axis z, and to adopt

C(r) = V 2 R exp{-(x 2 + y 2 + z 2 /ξ 2 )/σ 2 R }
, where V R is the amplitude of the disorder, σ R its typical correlation length, and ξ its anisotropy. It is the simplest model of anisotropic disorder and it is relevant to ultracold-atom experiments, see e.g. Refs. [8,10] where ξ varies from 6 to 20. When ξ > 1, this correlation function corresponds to a correlated disorder made of 'grains', of typical transverse size σ R and elongated in the z direction.

III. LOCALIZATION REGIME A. Exact numerical computation

The first step in the calculation consists in determining the average Green function, which enters into Eq. ( 1). It reads

G(E, k) = [E -2 k 2 /2m -Σ(E, k)] -1
, where Σ is given by the solution of the Dyson equation. In order to study the fundamental effects due to disorder anisotropy we compute Σ in the self-consistent Born approximation (SCBA), which reads

Σ(E, k) = dk ′ (2π) d C(k -k ′ ) E -2 k ′2 /2m -Σ(E, k ′ ) . (5)
This approximation reproduces qualitatively singleparticle features such as the density of states for Gaussian disorder. It is therefore well suited for understanding basic effects of geometrical disorder anisotropy as long as the disorder is not too large [34]. The major benefit of the SCBA method compared to the on-shell first-order perturbation theory is that it properly accounts for the momentum decay of the self-energy Σ(E, k).

We now numerically solve Eq. ( 1) for the current vertex function j by iteration (see Appendix B for details) and insert the result into Eq. ( 2) to determine the diffusion tensor. Due to parity and uniaxial symmetries, the diffusion tensor admits x, y, and z as eigendirections and is isotropic in the (x, y) plane. We call D ⊥ the corresponding diffusion coefficient and D z the coefficient along the z direction. In Fig. 1(a) we show the results for the diffusion coefficients (black points) for ξ = 3 and V R = 0.5E σ with E σ = 2 /mσ 2 R and in Fig. 1(b) the anisotropy D z /D ⊥ . For comparison, the corresponding values obtained in the Boltzmann approximation, i.e. if only U B is retained in Eq. ( 3), are also shown (blue points). We find that, as expected, both diffusion tensors are anisotropic with larger values in the z direction, and that both decrease with decreasing energy. The Boltzmann diffusion tensor D B (E) "trivially" vanishes at the point where the density of states vanishes (hereafter called band edge), E be ≃ -0.35E σ . As expected in the framework of AL, the quantum diffusion tensor D(E) vanishes well before D B (E). The energy dependence of D(E) (black points) points toward a vanishing point at E c > E be , for both D ⊥ and D z . It corresponds to the ME, the states with energy E be < E < E c being localized. The iterative calculation suffers from the usual critical slowing down phenomenon, which prevents us from performing exact calculations very close to the ME. In order to find the ME, we use the extrapolating functions

D z (E) ∝ (E -E c )[1 + A z (E -E c ) 2 ] and D ⊥ (E) ∝ (E -E c
), which fit the data properly. This extrapolation is represented by the solid black line in Fig. 1(a), and gives E c ≃ -0.1E σ .

B. Approximation

In order to avoid the critical slowing down in solving Eq. ( 1), it is worth simplifying the problem. It is vital to avoid any on-shell description since they would introduce artificial divergencies at short length scales. To do so, we assume here that the k-dependence of j(E, k) is equal to the one of j 0 (E, k), which amounts to assuming that the momentum distribution of the atoms at fixed energy does not change when interference terms are included. We have checked that this is indeed a fairly good approximation, as for isotropic models of disorder [START_REF] Yedjour | [END_REF]. Then, assuming j(E, k) = D(E)D -1 0 (E) j 0 (E, k) and

j B (E, k) = D B (E)D -1
0 (E) j 0 (E, k), we get the following closed equation for the diffusion tensor

D u (E) ≃ D u B (E) 1 + 2 2 D u (E) m 2 π 2 N (E) 2 D u 0 (E) 2 (6)
× dk (2π) d dk ′ (2π) d k u |G(E, k)| 2 j u,0 (E, k ′ ) ImΣ(E, k-k ′ 2 ) 2 (k + k ′ ) • D(E) • (k + k ′ ) ,
where u ∈ {x, y, z} are the three directions. Notice that there is a great gain in computational time as we have simplified an iterative problem for a vector of three-dimensional functions j(E, k) (respectively twodimensional in the uniaxial case we consider here) to an iterative problem for a tensor with only three (respectively two) parameters. This allows us to compute the diffusion tensor for energies down to the mobility edge E c . The results of this approximation are shown in Fig. 1. We find that the eigenvalues of the full diffusion tensor [Fig. 1(a)] and the anisotropy [Fig. 1(b)] are fairly well reproduced within the approximation. In particular, we recover that the eigenvalues of the diffusion tensor all vanish at the same critical point. This may be understood by a dimensionality reduction argument. Should localization occur in one or two directions only, the remaining (transverse) dynamics would be reduced to a dimension lower than three. It would then systematically lead to localization also in the transverse directions. The same conclusion can also be found from Eq. (6). Assume that the diffusion coefficient vanishes in the eigendirection u, but not in the others. Then the integral in Eq. ( 6) would be finite, so that D u (E) = D u B (E). This would be clearly inconsistent since D u (E) = 0 by hypothesis while D u B (E) = 0. Using the approximation, we locate the ME at E c = -0.09E σ for V R = 0.5E σ and ξ = 3, which is in good agreement with the value we found by extrapolating the exact numerical results, down to 10%.

C. Mobility edge

The ME, measured with respect to the band edge, E c -E be , is shown in Fig. 2(a) versus the disorder amplitude V R for various geometrical anisotropies. For low values of V R , the ME is very close to, but above, the band edge, well below the average value of the disorder, V = 0. When V R increases, the ME moves away from the band edge. Moreover, increasing the disorder anisotropy increases the ME at fixed V R and therefore tends to favor localization. This is intuitive since a strongly elongated disorder (ξ ≫ 1) tends to reduce the relative importance of the z direction. We also define the parameter kl(E) ≡ mTr [D 0 (E)] / , which coincides with the Ioffe-Regel parameter in the case of isotropic short-range disorder [35], since l(E) is then the scattering mean-free path. When computed at the ME, we find that kl is of order unity, it weakly depends on the disorder anisotropy, and it shows a slight increase with the disorder amplitude similar to that found for isotropic disorder [START_REF] Yedjour | [END_REF]. It thus provides a generalization of the Ioffe-Regel criterion to the anisotropic case.

(E c -E be )/E σ R V R /E σ R ξ=1 ξ=2 ξ=3 0.8 1 1.2 1.4 1.6 1.8 2 0 0.2 0.4 0.6 0.8 1 D z (E c )/D ⊥ (E c ) V R /E σ R a) b) App, ξ=2 App, ξ=3 0, ξ=2 0, ξ=3 B, ξ=2 B, ξ=3

D. Transport anisotropies

We now discuss the transport anisotropies, D z B /D ⊥ B and D z /D ⊥ , which are shown in Fig. 1(b) versus the matterwave energy E. We find that both decrease with energy, similarly as found in the on-shell approximation [25]. The Boltzmann diffusion tensor becomes isotropic only at very low energy when it approaches the band edge. This is consistent with the fact that our model has a nondiverging infrared limit, lim k→0 C(k) = const. Then, for low energy, all integrals are dominated by small wavevectors k σ -1 R , which are almost insensitive to the details of the correlation function, in particular to its anisotropy. Most important, we find that the anisotropy of the quantum diffusion tensor D(E) is significantly smaller than that of the Boltzmann diffusion tensor D B (E). It shows that the localization corrections strongly affect the transport anisotropy, leading to a significant reduction for our model of disorder. For instance in Fig. 1(b), when E → E c , the coherent transport anisotropy D z /D ⊥ strongly decreases and becomes of the order of 1, while the anisotropy of the incoherent transport D z B /D ⊥ B is of the order of 1.75. This effect is the main outcome of our approach, which is free of any cutoff of arbitrary anisotropy, unlike previous work [25,29].

Figure 2(b) shows the anisotropies of the diffusion tensors at the ME, where the deviations are more pronounced. Isotropic diffusion tensors are found only in the white-noise limit V R ≪ E σR [START_REF]The limit VR/Eσ R → 0 is obtained by taking σR → 0 while V 2 R σ 3 R is fixed[END_REF]. Conversely, for larger values of V R , finite k components of the disorder contribute. Then, we find that the self-energy Σ becomes only slightly anisotropic and thus, that the anisotropy of D 0 is weakly affected. In contrast, the anisotropy of D B is significantly affected. This is due to the second term in Eq. ( 1), which couples j(E, k) to the disorder power spectrum C(k) in the Boltzmann approximation. Most importantly, we find that the anisotropy of the quantum diffusion tensor D is significantly smaller than that of D B . It confirms the strong renormalization of the transport anisotropy by localization corrections. To understand this behavior, let us rely on Eq. ( 6), which relates D to D B and properties of the disorder-average effective scattering medium. Conversely, away from the ME, the localization corrections [second term in the bracket of Eq. ( 6)] vanish and D ≃ D B . When approaching the ME localization corrections increase, which tend to lower the anisotropy of D. More precisely, at criticality, D ≃ 0 while D B = 0. Hence, D is solely determined by the cancellation of the term inside the brackets in Eq. ( 6) and D B becomes irrelevant. We can roughly estimate the anisotropy at criticality by neglecting the width of the pole in Eq. ( 6) and estimating all integrals in the range where their integrands are maximal. This yields

D z (E c )/D ⊥ (E c ) ∼ D z 0 (E c )/D ⊥ 0 (E c
) for short-range correlations, which is consistent with the weak anisotropy of the quantum diffusion tensor at criticality found numerically in Fig. 2(c).

IV. CONCLUSION

In summary, we have shown that localization significantly affects transport anisotropy, especially near the ME. Our results show that the anisotropy of the coherent diffusion tensor is mainly determined by the properties of the disorder-averaged effective scattering medium while the incoherent transport contributions are significantly compensated by interference when approaching the ME. More precisely, we find a small but finite anisotropy of transport near the ME. Those results are obtained in the framework of the approximation leading to Eq. ( 6), which is satisfactory for the parameters used in this work. However, it slightly underestimates the anisotropy. Relaxing this condition may slightly affect the coherent transport anisotropy. Our results indicate that a significant deviation from isotropy could be found in the case of strongly anisotropic scattering properties, a case that should probably be investigated by going beyond the SCBA for computing properties of the scattering medium. Our results provide insights into the theory of Anderson localization in anisotropic media and are directly relevant to experiments on ultracold atoms in three-dimensional anisotropic disorder [8][9][10][START_REF] Semeghini | Measurement of the mobility edge for 3D Anderson localization[END_REF]. In the future, it would be interesting to compare our predictions to direct numerical calculations using for instance an extension of the method used in Ref. [34] to anisotropic disorder.

V. ACKNOWLEDGMENTS find j(E, k) = k 2 |G(E, k)| 2 - 1 2 ∇ k ReG(E, k) (A7) + 1 2 |G(E, k)| 2 ∇ k ReΣ(E, k) +|G(E, k)| 2 dk ′ (2π) d j(E, k ′ ) U k,k ′ (E, 0) -i|G(E, k)| 2 dk ′ (2π) d ImG(E, k ′ ) ∇ q U k,k ′ (E, q)| q=0 .
Finally, we find the following expression:

j(E, k) = j 0 (E, k) (A8) + |G(E, k)| 2 dk ′ (2π) d j(E, k ′ ) U k,k ′ (E, 0) -i|G(E, k)| 2 dk ′ (2π) d ImG(E, k ′ ) ∇ q U k,k ′ (E, q)| q=0 . with j 0 (E, k) = [k + ∇ k ReΣ(E, k)] ImG(E, k ′ ) 2 - 1 2
ReG(E, k)ImG(E, k)∇ k ImΣ(E, k). In this paper, we have neglected the last term in Eq. (A8) as ∇ q U k,k ′ (E, q) vanishes in the Boltzmann approximation, and is less divergent than U k,k ′ (E, q) itself for the maximally-crossed diagrams.

Figure 1 .

 1 Figure 1. (Color online) Quantum coherent diffusion versus Boltzmann diffusion. (a) Diffusion coefficients in the transverse plane ⊥= {x, y} and in the longitudinal direction z computed from the numerical solution of Eq. (1) for our model of disorder with ξ = 3 and VR = 0.5Eσ. Shown are the results for the incoherent Boltzmann approximation (U ≃ U B , blue points), the quantum diffusion (U = U B + U MC , black points), and the approximation underlying Eq. (6) [red points]. The black solid line shows the extrapolation D ⊥ (E) ∝ (E -Ec) and D z (E) ∝ (E -Ec)[1 + A z (E -Ec) 2 ], which gives Ec = -0.102Eσ. (b) Anisotropy of the quantum and Boltzmann diffusion tensors, D z (E)/D ⊥ (E).

Figure 2 .

 2 Figure 2. (Color online) (a) Position of the mobility edge with respect to the band edge Ec -E be as a function of the amplitude of the disorder in approximation (6), for geometrical anisotropies ξ = 1, 2 and 3 (black, red and blue data points respectively). (b) Anisotropy of the diffusion tensor at the mobility edge found within the approximation for ξ = 2 and 3 (red and blue data points respectively), compared with the anisotropy of D0 and DB at the same energy (see legend).

We thank Serguey Skipetrov for enlightening discussions. This research was supported by the European Research Council (Contract No. 256294), the Agence Nationale de la Recherche (Contract No. ANR-08-blan-0016-01), the Ministère de l'Enseignement Supérieur et de la Recherche, and the Institut Francilien de Recherche sur les Atomes Froids (IFRAF). We acknowledge GM-PCS computing facilities of the LUMAT federation and HPC resources from GENCI-IDRIS (Grants No. 2013057143 and No. 2014057143).

Appendix A: Quantum transport formalism

The two-particle Green function Φ k,k ′ (q, ω, E) describes the density propagation of a wave in a disordered medium. It is defined as Φ k,k ′ (q, ω, E) ≡

|k -with G the retarded Green operator, k ± ≡ k ± q/2 and k ′ ± ≡ k ′ ± q/2, E ± ≡ E ± ω/2, and (q, ω) the Fourier conjugates of the space and time variables [START_REF] Here | we use f (q, ω) ≡ drdt f (r, t) exp[-i(q • r -ωt)[END_REF].

It is governed by the Bethe-Salpeter equation, which can be formally written as [START_REF] Rammer | Quantum Transport Theory[END_REF]

The first term in Eq. (A1) describes uncorrelated propagation of the field and its conjugate in the disordered medium. The second term involves the vertex function U, which includes all irreducible scattering diagrams, and accounts for all correlations in the density propagation.

The associated quantum kinetic equation reads [START_REF] Rammer | Quantum Transport Theory[END_REF][START_REF] Vollhardt | [END_REF] [

where ǫ 0 (k) is the free dispersion relation (ǫ 0 (k) = 2 k 2 /2m in our case, for particles in free space) and

The conservation of the total quantum probability gives the Ward identity which relates the vertex function to the singleparticle self-energy Σ associated to the average Green function G:

The conservation of quantum probability, through the Ward identity, guarantees a diffusion pole. In the long time (ω → 0) and large distance (q → 0) limit, one can therefore write Φ as [START_REF] Rammer | Quantum Transport Theory[END_REF]39]:

A4) where D(E) is the diffusion tensor, and φ(E, k, q) is the eigenfunction of the Bethe-Salpeter equation associated with the hydrodynamic diffusion [which corresponds to eigenvalue 0 as (q, ω → 0)]. The linearization of φ at small q gives φ(E, k, q) = -ImG(E, k) -iq • j(E, k) [START_REF] Yedjour | [END_REF], where j is the static current vertex function [29,[START_REF] Mahan | Many Particle Physics[END_REF].

Fick's law, which relates the diffusive flux to the concentration gradient in the diffusive regime (ω, q → 0) reads J (E, ω, k, q) = -iD(E)qP(E, ω, k, q) in Fourier space.

The current density reads

, and the probability density E,ω,q). It gives the Kubo formula for the diffusion tensor

Expanding the quantum kinetic equation (A2) into terms linear in q gives

Developing also the Ward identity (A3) permits us to

Appendix B: Numerical method When solving Eq. ( 1) for the current vertex function j(E, k), we proceed by iterations. At stage n ≥ 1 we compute j (n) (E, k) from j (n-1) (E, k) and the corresponding D (n-1) (E) according to

and compute D (n) (E) from Eq. ( 2) We first solve the Boltzmann case (U = U B ), which we initialize with j (n=0) = j 0 . In the full calculation, we then use j (n=0) = j B , which improves the convergence. The procedure is iterated until convergence of the norm of j.

The functions j(E, k) are represented in cylindrical coordinates, which are adapted to the symmetry of our problem. We discretize the integrals on a grid in spherical coordinates, with up to 500 points in the radial direction and 128 points for the angular coordinates. The different functions appearing in the integral are interpolated using spline interpolation. We have checked that increasing the grid generates negligible changes in the final results and that our procedure gives consistent results in the isotropic case.


Guillaume Bécan
email: gbecan@irisa.fr

Mathieu Acher
email: mathieu.acher@irisa.fr

Jean-Marc Jézéquel
email: jean-marc.jezequel@irisa.fr

Thomas Menguy
email: thomas@wildmoka.com

On the Variability Secrets of an Online Video Generator

Keywords: D.6.5 [Management of Computing and Information Systems]: Security and Protection, D.2.7 [Software Engineering]: Distribution, Maintenance, and Enhancement-Restructuring, reverse engineering, and reengineering, D.2.9 [Software Engineering]: Management-Software configuration management Reverse engineering, Software product line, Configurator, Video generator, Security

We relate an original experience concerning a popular online video service that offers to generate variants of an humorous video. To further the understanding of the generator, we have reverse engineered its general behavior, architecture, as well as its variation points and its configuration space. The reverse engineering also allows us to create a new generator and online configurator that proposes 18 variation points -instead of only 3 as in the original generator. We explain why and how we have collaborated and are collaborating with the original creators of the video generator. We also highlight how our reverse engineering work represents a threat to the original service and call for further investigating variability-aware security mechanisms.

INTRODUCTION

The concepts of software product lines and variability are more and more observed in numerous domains and contexts [START_REF] Apel | Feature-Oriented Software Product Lines: Concepts and Implementation[END_REF][START_REF] Erwig | The choice calculus: A representation for software variation[END_REF][START_REF] Svahnberg | A taxonomy of variability realization techniques: Research articles[END_REF][START_REF] Thüm | A classification and survey of analysis strategies for software product lines[END_REF]. Abstractions, theories, techniques, languages, and tools have been developed and are widely used for managing and building variants of user interfaces, requirements, source code, architectures, 3D models, or design models, to name a few [START_REF] Alves | An exploratory study of information retrieval techniques in domain analysis[END_REF][START_REF] Chen | An approach to constructing feature models based on requirements clustering[END_REF][START_REF] Czarnecki | Generative Programming: Methods, Tools, and Applications[END_REF][START_REF] Galindo | A variability-based testing approach for synthesizing video sequences[END_REF][START_REF] Nadi | Mining configuration constraints: Static analyses and empirical results[END_REF][START_REF] Niu | Concept analysis for product line requirements[END_REF][START_REF] Sawyer | Using constraint programming to manage configurations in self-adaptive systems[END_REF][START_REF] Weston | A framework for constructing semantically composable feature models from natural language requirements[END_REF][START_REF] Ziadi | Towards a language-independent approach for reverse-engineering of software product lines[END_REF]. The applicability of variability techniques seems infinite -any artefact is virtually subject to customisation, configuration, or extension.

We report on an experience related to an online video service. The service offers to generate variants of an humorous video. Internet users simply have to type their name, select 3 options (or click on a button for a random selection), and a particular video is launched and visualised in the browser (see Figure 1, page ). The service is quite popular and successful: more than 1.7M of video variants have been generated in 1 week.

We had an interest in the Web site, since the service claims that a particular video is resulting from a combination among billions. Moreover the term generator is explicitly employed. Our original intuition was that the video service resembles to a software product line, i.e., generative techniques and variability are likely to be present.

We have first noticed a common property of all the video variants: 18 sequences of videos are assembled; for each sequence of a video, numerous alternatives are possible. These variants represent a product line of videos. We have then contacted the creators of the video generator with two research questions and goals in mind: (1) can we automate the extraction of variation points and configurations (video variants)? (2) can we re-engineer another generator and configurator?

The creators accepted to collaborate and provided us with an access to an offline version of the generator. However we did not have access to the source code neither in the server side nor in the client side. Like this, we can study the generator as a black box system -as it would be the case in reality -but without disturbing the deployed Web generator (e.g., with hundreds of HTTP requests). The interest for the creators was to audit their systems w.r.t. variability.

To address the two research questions and further the understanding of the generator, we have reverse engineered its general behavior, architecture, as well as its variation points and its configuration space. Overall, we have collected 350,000+ configurations and 1620 video sequences. We have automatically determined to which part of the video the 1620 sequences correspond. We have also computed some statistics to understand the configuration process. The paper describes the reverse engineering process, how we gather variability information and insights about the configurations, and points out the importance of barriers for complicating or avoiding a large-scale extraction of configurations.

The reverse engineering work have also allowed us to reengineer a new generator and configurator. This time, users can configure in a fine-grained way the 18 variations points -instead of only 3 in the original version. 15 new variation points are now apparent, out of which 13 are configurable (i.e., exhibiting at least 2 alternatives). For the 3 variation points also present in the original configurator, users can now explicitly choose the alternatives. Thumbnails representative of the video sequences are depicted and help the users during the selection. The paper describes the re-engineering process, highlights the large degree of automation of devising a new configurator, and discusses the differences with the original configurator.

Finally, we explain why the re-engineering of a configurator and the exhibit of new variation points could be a threat for the creators of the original video generator. The main reason is that the new configurator lets users control, choose and visualize any alternative -thus dramatically limiting the surprise effect. Our experience shows that the scoping of variability (i.e., what is visible to users) is a strategic solution and, as such, mechanisms to prevent the reverse engineering of variability have to be considered. We report on our preliminary exchanges with the creators of the video generator.

This paper is both an experience report and a research-inprogress. The intended audience are (1) researchers working on reverse engineering of variability; (2) developers of generators, configurators, and variability-intensive systems. The former shall benefit from an extraction and analysis of a large configuration set -in an original setting. The latter can learn from our experience of (re-)engineering a configurator and scoping variability. We also highlight the need to develop defensive mechanisms for increasing the difficulty of reverse engineering variability. Besides, we believe that the originality of the considered artefacts (videos, Web applications), the popularity of the online service, and the intuitiveness of the application domain can be of interest to a broad audience -not necessarily experts of product lines or variability (e.g., students).

The remainder of the paper is organized as follows. Section 2 introduces the online video generator "Bref, 30 ans". Section 3 presents the reverse engineering process which allows us to understand in a fine-grained way the variability and the general behavior of the generator. Section 4 describes the re-engineering of the configurator and generator. Section 5 discusses related work. Section 6 summarizes the contributions of the paper and concludes.

"BREF, 30 ANS": AN ONLINE VIDEO GENERATOR

Bref, 30 ans: A Family of Video Variants

Canal+ (a French premium cable television channel) launched several initiatives for celebrating its 30th anniversary. Among others, the Web site http://bref30ans.canalplus.fr has been developed and provides an online service to generate variants of an humorous video. A video variant follows a story line in which the historic perception of Canal+ by the main character of Bref 1 is presented in an humoristic way. A video variant typically includes sequences related to famous programs of Canal+ while several guests (e.g., actors, presenters) are part of the story line.

Internet users simply have to type their name (step 1) and select 3 options (step 2, see Figure 1). A particular video is launched and visualised in the browser. The information of step 1 is used during the first sequence (the name appears in a blank screen). The three options of step 2 correspond to programs ("souvenirs" in Figure 1) of Canal+ that a user wants to include in the variant. Users can also choose a randomized selection, corresponding to "Je ne veux pas choisir" 1 Bref is a popular television show promoted by Canal+ (see Figure 1) that can be translated as "I do not want to choose". The service is quite popular and successful. More than 1.7M of video variants have been generated in 1 week.

Motivation and Research Questions

We had an interest in the Web site since the video service resembles to a software product line. From a research perspective, the online video generator is likely to constitute an interesting case study (i.e., "an empirical inquiry that investigates a contemporary phenomenon within its real-life context" [START_REF] Yin | Case Study Research: Design and Methods[END_REF]). From a dissemination perspective, the relative accessibility, intuitive nature, and concreteness of the case study are good qualities for a broad audience. For example, the case study could be used for presenting concepts of product lines and variability to students [START_REF] Acher | A survey on teaching of software product lines[END_REF]. Interestingly, the main character of Bref states, at the very end, that a particular video is resulting from a combination among billions; the term generator is even explicitly employed. Hence we suspect that generative techniques and variability are likely to be present.

We sought to verify this assumption and considered the following research questions: What is the variability and commonality of the generator? How is variability modeled and implemented? How many variation points, commonalities, constraints, and configurations are there? How configurations and video variants are generated on-demand?

In short, our original motivation was to understand in a fine-grained way the online video generator. Another related motivation was to determine whether we can re-engineer the generator and provide an alternative configurator, e.g., with more than 3 options. For both understanding and reengineering the generator, we adopted a reverse engineering approach. In essence, reverse engineering "consists in deriving information from the available software artifacts and translating it into abstract representations more easily understandable by humans" [START_REF] Canfora | Achievements and challenges in software reverse engineering[END_REF]. Here, we wanted to extract variability information and to translate it into an understandable form, so that we can extract insights and start the re-engineering of a new configurator.

Research Method

After some preliminary observations, we contacted the creators of the video generator for collaborating on the topic. The creators accepted and provided us with an access to an offline version of the generator. The interest for the creators was to audit their systems w.r.t. variability. Importantly, we did not have access to the original source code neither in the server side nor in the client side (e.g., the JavaScript was minimized and obfuscated, as in the original version). By doing so, we wanted to replicate the exact conditions in which is confronted a person in charge of reverse engineering (e.g., an attacker). Overall, we studied the generator as a black box system -as it would be the case in reality -but without disturbing the deployed service. For instance, we performed hundreds of HTTP requests without overwhelming the network with traffic.

REVERSE ENGINEERING VARIABILITY SECRETS

Reverse engineering scenarios usually assume an access to the source code or any form of documentation, for example, for achieving design recovery of an existing application [START_REF] Canfora | Achievements and challenges in software reverse engineering[END_REF]. It is not the case in our context. Not having access to the original source code of the application was a deliberate decision, challenging us to retrieve the variability secrets of the generator.

Preliminary Inspection

We started the investigation of the Web generator through a manual analysis of its behavior. We first sought to understand, at a high level, what are the properties of the generated videos. We manually generated some videos and we reported the following observations:

• the story line seems fixed, with some recurrent sequences, as well as the length of the video (3 minutes); • the variability not only concerns the 3 souvenirs of the configurator, but also some parts of the video (e.g., some actors are replaced by others)

In a sense, we quickly verified the intuition that the generator exhibits commonality and variability.

We then investigated, at the technical level, how the videos are actually produced. We relied on basic Web tooling support to analyze the JavaScript and the requests with the server; we learned that:

• a video variant is a combination of video sequences;

• each video is encoded in .ts format and accessible via a specific URL, based on the filename; • there is a specification for executing the video sequences in order (a playlist-like file).

At this step, we understood the general principle of the generator but ignored (or had some doubts) about the variation points (e.g., how many variation points are there in the generator?).

Large Scale Extraction of Configurations

Due to the limits of a manual inspection, we thus sought to automatically extract as many video variants as possible. We first considered the use of crawlers to emulate the user behavior when selecting options in step 2 (see Figure 1); then we could get videos. However the approach is not really scalable since it requires to download every sequence of video, for each variant.

We thus decided to reverse engineer the way the playlist is retrieved. We analyzed the HTTP and JSON requests of the client with the server. We came to the conclusion that a playlist is documented as an ordered list of 18 string-based values (an example is given below):

"sq": ["dwlcjv", "1y60t3z", "1lyfhk", "wqzv0y", "1xxivi2", "1oxnvtu", "lolbe9", "wvo06o", "1u6y5t2", "1eqb8bw", "1j9aij7", "nr7jom", "1jmv11y", "1qgn9dh", "1bv7rka", "19ykyyw", "5znrg7", "116hv1k"] Each value of the ordered list actually corresponds a subplaylist, accessible via an URL. For example, for ml7ila: #EXTINF:05.72, http://XXX/../ml7ila_med0.ts #EXTINF:05.96, http://XXX/../ml7ila_med1.ts #EXTINF:03.96, http://XXX/../ml7ila_med2.ts #EXTINF:02.12, http://XXX/../ml7ila_med3.ts From a conceptual point of view, the ordered list corresponds to a configuration (we will use this term in the remainder of the paper). Overall we automatically generated 363,281 configurations. We also retrieved the sub-playlists corresponding to each value of the configurations. We combined wget and curl to operate over URLs.

We stopped at 350,000+ configurations as we did not observe new values composing a configuration, compared to a 50,000+ configurations sample we first obtained. An attacker is likely to limit the number of requests in order to not draw attention. A defender can also apply defensive strategies to limit the requests. An interesting research direction is thus to determine the number of requests for which the reverse engineering and understanding remains effective (e.g., comprehensive). We discuss this perspective in the conclusion.

Automated Analysis of Configurations

We processed the collected set of configurations to further understand the variability of the generator. We observed that each configuration is constituted of 18 values, that is, there are 18 variation points (we will use this term in the remainder of the paper). We collected some statistics for each of the variation point (see Table 1). 2 variation points, VP2 and VP18, have only 1 alternative. We can consider them as two commonalities or mandatory features of a video variant. 16 variation points exhibit at least 2 alternatives and are thus configurable.

The number of alternatives varies for a given variation point. VP1, VP7, and VP16 correspond to "souvenirs" and the three options originally proposed in the video generator of Bref, 30 ans. These variation points have the highest number of alternatives (resp. 63, 51, and 86). There are only two alternatives for VP9. The first alternative appears in 362,903 configurations; the second alternative only 378 times, i.e., 0.1 % of the configurations contains this video sequence. More generally, we observed that the frequencies of alternatives for a given variation point: Besides we observed hard-constraints between the alternatives within the sample, e.g., there are some logical, binary exclusions. For example, some alternatives of VP16 are mutually exclusive with the least frequent alternative of VP9. Similarly there are soft constraints [START_REF] Czarnecki | Sample spaces and feature models: There and back again[END_REF] within the sample, i.e., there is a conditional probability of a configuration to contain one alternative together with another. It is tempting to interpret these hard-or soft-constraints as design decisions of the creators. However, we recall that the configuration set is only a subset of the whole configuration set. Some hard-or soft-constraints may be due to the incompleteness of the extracted configuration set. It is thus premature to draw any conclusions.

•
Finally, it should be noted that an alternative is constituted of a series of videos (e.g., in the example of ml7ila, there are 4 videos); in total, we collect 1620 video sequences. Another form of variability observed in the generator is as follows. The videos constituting an alternative may be partly assembled. For example, in the case of ml7ila, only the 3 last videos (instead of the 4 videos) can be assembled -the goal is to reduce the length of the video or propose yet another subtle variant of a video. A thorough analysis of the obfuscated client side code could reveal such kind of variability. We leave it as future work.

Overall there is a deliberate intention and strategy for surprising the users of the video generator. In the same vein, it should be noted that each configuration of our sample is unique. The rationale is certainly that users have to obtain a new video variant each time they visit the online service.

RE-ENGINEERING THE GENERATOR AND CONFIGURATOR

The result of the previous reverse engineering process is as follows: It is simply an ordered combination of video sequences.

Our objective was to demonstrate the feasibility of devising a new generator and configurator -this time exhibiting all variation points. The re-engineering of the same exact version of the generator could have been another objective.

In particular we could have tried to replicate the process of generating unique configurations that respect the frequencies of alternatives. We leave it as future work.

New Configurator

By reusing our observations on how the variability is realized in the original service and by analyzing some parts of the Javascript code, we were able to re-engineer a new configurator. The result is depicted in Figure 2. Users can configure in a fine-grained way the 16 variations pointsinstead of only 3 in the original version -and visualize the result of the two common videos of all variants. For instance, the variation point "Tu m'écoutes", corresponding to VP8 in Table 1, depicts 30 alternatives. For the 9th sequence (VP9, "Jean Jacques") two videos can be selected, etc.

15 new variation points are now apparent, out of which 13 are configurable (i.e., exhibiting at least 2 alternatives). For the 3 variation points also present in the original configurator (the "souvenirs"), users can now explicitly choose the alternatives.

For each variation point, thumbnails representative of the video sequences are depicted and help the users during the selection. We generated the thumbnails from the videos sequences using ffmpeg 2 . Overall the creation of the new configurator was largely automated thanks to the result of the reverse engineering process. The manual task consists in giving a proper name to each variation points (in the form of french sentences) and engineering the generator for playing the video sequences. It should be noted that our generator assembles and reuses videos located in the servers of the original service. 2 https://www.ffmpeg.org/

The Threats of Variability Secrets

We presented the re-engineered configurator to the original creators of the video generator. The new configurator represents an important threat and could "kill the original idea". That is, with our solution, the surprise effect is limited when getting a new video variant. Instead, users can control, choose and visualize any alternative (video sequence). The creators of the generator did not have this intention. In other words, they have taken a strategic decision for, onpurpose, delimiting the scope of variability and restraining the visibility of some variation points and alternatives.

The main lessons we can learn from this experience is that (1) the large scale extraction of configurations is key and allows us to re-engineer another competing configurator that could have annihilated the original innovation of the generator; (2) whenever the scoping of variability does matter, mechanisms should be used to prevent or limit the mining in the large of non visible alternatives.

With regards to the last point, some defensive techniques have been considered, for instance, the use of block lists for banning IP addresses that perform numerous HTTP requests in a short amount of time. Other security mechanisms have been discussed with the creators of the video generator but are out of the scope of this paper. Our experience so far shows that the problem of protecting the variability requires a careful attention and certainly the use of innovative techniques.

RELATED WORK

Variability and Video

The use of variability techniques in the video domain is not new and has already been reported [START_REF] Acher | Modeling variability from requirements to runtime[END_REF][START_REF] Galindo | A variability-based testing approach for synthesizing video sequences[END_REF][START_REF] Moisan | Run Time Adaptation of Video-Surveillance Systems: A Software Modeling Approach[END_REF]. However the underlying motivations and realizations vastly differ. In [START_REF] Acher | Modeling variability from requirements to runtime[END_REF][START_REF] Moisan | Run Time Adaptation of Video-Surveillance Systems: A Software Modeling Approach[END_REF], Moisan et al. model variability of video processing algorithms. The goal is to systematize the parametrization of highly configurable algorithms composing a processing chain. The algorithms exhibit variability at runtime and are able to adapt to evolving and variable contexts (e.g., a change of luminosity in the video).

In [START_REF] Galindo | A variability-based testing approach for synthesizing video sequences[END_REF], Galindo et al. presented an approach for generating synthetic variants of a video. The developed generator relies on a base (a video sequence) that is modified for changing the luminosity, adding elements (e.g., vehicles) in the scene, etc. The derivation process is thus quite different. Instead of assembling pre-defined video sequences, the derivation consists in altering an existing video sequence through a series of transformations. The techniques are more advanced and are based on video processing algorithms. Another important difference concerns the objective. The industrial motivation of [START_REF] Galindo | A variability-based testing approach for synthesizing video sequences[END_REF] is that companies providing or consuming algorithms have severe difficulties to evaluate or compare the quality of existing solutions. The use of variability techniques allows practitioners to generate as much as possible videos, with different characteristics, in order to test or benchmark algorithms on a rich and diverse dataset.

It is possible to apply the video processing techniques exposed in [START_REF] Galindo | A variability-based testing approach for synthesizing video sequences[END_REF] in the case of the online video generator. We envision to process variants of pre-defined video sequences and then assemble the new variants for deriving a comprehensive video.

Reverse Engineering Variability

Numerous works address the reverse engineering of variability [1, 5-7, 10, 12-14, 18-20, 23, 26, 27, 31, 33, 33, 36, 38, 39].

The artefacts considered as input are product descriptions [START_REF] Czarnecki | Sample spaces and feature models: There and back again[END_REF][START_REF] Davril | Feature model extraction from large collections of informal product descriptions[END_REF][START_REF] Haslinger | Reverse engineering feature models from programs' feature sets[END_REF][START_REF] Haslinger | On extracting feature models from sets of valid feature combinations[END_REF], requirements [START_REF] Alves | An exploratory study of information retrieval techniques in domain analysis[END_REF][START_REF] Chen | An approach to constructing feature models based on requirements clustering[END_REF][START_REF] Niu | Concept analysis for product line requirements[END_REF][START_REF] Weston | A framework for constructing semantically composable feature models from natural language requirements[END_REF], configurators [START_REF] Abbasi | Reverse engineering web configurators[END_REF], source code [START_REF] Al-Msie'deen | Mining features from the object-oriented source code of software variants by combining lexical and structural similarity[END_REF][START_REF] Nadi | Mining configuration constraints: Static analyses and empirical results[END_REF][START_REF] She | Reverse engineering feature models[END_REF][START_REF] Ziadi | Towards a language-independent approach for reverse-engineering of software product lines[END_REF], or models [START_REF] Martinez | Identifying and visualising commonality and variability in model variants[END_REF]. The objectives are usually to locate features or configuration options, to formalize logical relations, and possibly to re-engineer a system with a generative, variability-based approach.

The experience report of this paper leans on (1) an intensive and manual observation of a generator (2) the use of automated techniques for a large-scale extraction of configurations and a classification of variations points. The approach has the originality of targeting videos and Web artefacts. Abbasi et al. [START_REF] Abbasi | Reverse engineering web configurators[END_REF] propose generic techniques to reverse engineer Web configurators, but do not develop techniques to extract configurations. Our reverse engineering process leads to a fine-grained understanding of variability (see Section 3) and a re-engineering of the configurator (see Section 4).

Configuration

Techniques for automatic reasoning have been proposed to manage the configuration process, being for guiding users (e.g., [START_REF] Mazo | Recommendation heuristics for improving product line configuration processes[END_REF][START_REF] Nöhrer | Optimizing user guidance during decision-making[END_REF]) or for piloting an adaptation [START_REF] Sawyer | Using constraint programming to manage configurations in self-adaptive systems[END_REF]. Decision or feature-based configuration techniques and environment have been proposed [START_REF] Abbasi | A toolset for feature-based configuration workflows[END_REF][START_REF] Dhungana | Integrating heterogeneous variability modeling approaches with invar[END_REF][START_REF] Nöhrer | Optimizing user guidance during decision-making[END_REF][START_REF] Pleuss | Visualization techniques for application in interactive product configuration[END_REF][START_REF] Rabiser | A qualitative study on user guidance capabilities in product configuration tools[END_REF]. Web configurators are engineered along different practices [START_REF] Khalil Abbasi | The anatomy of a sales configurator: An empirical study of 111 cases[END_REF] for representing options, managing constraints, guiding users, scheduling the configuration process, etc. We can reuse and apply part of these techniques when re-engineering the configurator (see Section 4).

CONCLUSION

We reported an experience in reverse engineering and reengineering the variability secrets of an online video generator. We described our collaboration with the original creators of the video generator to audit their systems w.r.t. variability. The contributions can be summarized as follows:

• We described how we have inferred the general behavior of the generator (including its variation points) and how we have automatically extracted and classified 1620 video sequences out of 350,000+ configurations; • We presented the creation of a new configurator exposing 18 variation points instead of only 3 as in the original video generator; • We explained why letting an attacker re-engineers a new configurator is a threat for the creators, thus motivating the use of defensive mechanisms for complicating his or her task.

Education material. Describing the implementation of the video generator (or even relating our experience) can be considered for illustrating the concept of generative techniques and variability. We believe the concreteness and the relative accessibility of the case study are likely to reach a broad audience (e.g., students [START_REF] Acher | A survey on teaching of software product lines[END_REF]). At University of Rennes 1 (graduate level), we have already started to illustrate the software product line course with the video generator.

Lessons learned. A key aspect of our work is that the ability to reverse engineer a large number of configurations dramatically helps to understand the variability and behavior of the generator. It also helps us to devise a quite comprehensive configurator -at least competing with the original solution. Our conclusion -and the main lesson learned -is that defensive mechanisms should be considered and developed to prevent such large-scale extraction of configurations.

Future work. The main limitation of our current and ongoing work is that we do not yet validate our understanding and re-engineering with the creators of the generator. For instance, we ignore if the extraction of alternatives and 1620 video sequences is complete. When considering the sample of 350,000+ configurations, we did not observe new alternatives compared to 50,000+ configurations. However we may still have missed very low frequencies of alternatives. Similarly we cannot draw definitive conclusions w.r.t. hard and soft constraints between alternatives since the configuration set is by construction only a sample. Mining a larger number of configurations is a possible strategy and can increase the statistical significance of the results. Generally speaking it is hard to cover the whole configuration space due to the exponential combinations of alternatives. An even more ambitious research direction is to characterize the minimum number of configurations an attacker needs to have an accurate and complete understanding of the generator and configurator. We are currently continuing the discussions with the creators of the video generator along these lines. A comprehensive report of our future exchanges is planned but is out of the scope of this research-in-progress paper.

We conjecture that similar security issues, due to the presence of hidden variability, arise in other contexts -beyond the case of the Web video generator. Some proprietary systems or applications hide on purpose some features. Users typically have to pay for activating the features. For example, Windows 7 was subject to a hack that allows a clean install of the operating system using an upgrade disc rather than the full version upgrades [START_REF] Ligman | Regardless of what any hack says, a windows 7 upgrade is an upgrade. what you need to know[END_REF]. It calls for investigating variability-aware security mechanisms.

Figure 1 :

 1 Figure 1: Online video generator, step 2/3 of the configuration process (left-hand side): a user can select 3 options ("souvenirs") or let the system randomly choose. A variant of a video is then played (step 3).

Figure 2 :

 2 Figure 2: Re-engineering of the novel configurator (excerpt): users can now select a specific video for the 18 variation points identified during the reverse engineering of the original generator

Table 1 :

 1 Number of alternatives per variation point• an abstraction of variation points and alternatives. We can infer a variability model (e.g., a feature model); • an understanding of how video variants are generated.

	Variation Point # Alternatives
	VP1	63
	VP2	1
	VP3	9
	VP4	34
	VP5	15
	VP6	25
	VP7	51
	VP8	30
	VP9	2
	VP10	6
	VP11	6
	VP12	12
	VP13	21
	VP14	28
	VP15	6
	VP16	86
	VP17	4
	VP18	1