Vincent Despr
email: vincent.despre@loria.fr

Jean-Marc Schlenker
email: jean-marc.schlenker@uni.lu

Monique Teillaud
email: monique.teillaud@inria.fr.

FLIPPING GEOMETRIC TRIANGULATIONS ON HYPERBOLIC SURFACES

We consider geometric triangulations of surfaces, i.e., triangulations whose edges can be realized by disjoint locally geodesic segments. We prove that the flip graph of geometric triangulations with fixed vertices of a flat torus or a closed hyperbolic surface is connected. We give upper bounds on the number of edge flips that are necessary to transform any geometric triangulation on such a surface into a Delaunay triangulation.

Introduction

In this paper, we investigate triangulations of two categories of surfaces: flat tori, i.e., surfaces of genus 1 with a locally Euclidean metric, and hyperbolic surfaces, i.e., surfaces of genus at least 2 with a locally hyperbolic metric (these surfaces will be introduced more formally in Section 2.1).

Triangulations of surfaces can be considered in a purely topological manner: a triangulation of a surface is a graph whose vertices, edges and faces partition the surface and whose faces have three (non-necessarily distinct) vertices. However, when the surface is equipped with a Euclidean or hyperbolic structure, it is possible to consider geometric triangulations, i.e., triangulations whose edges can be realized as interior disjoint locally geodesic segments (Definition 2.1). Note that a geometric triangulation can still have loops and multiple edges, but no contractible loop and no contractible cycle formed of two edges. We will prove that any Delaunay triangulation (Definition 2.3) of the considered surfaces is geometric (Proposition 3.3).

The flip graph of triangulations of the Euclidean plane has been well studied. It is known to be connected; moreover the number of edge flips that are needed to transform any given triangulation with n vertices in the plane into the Delaunay triangulation has complexity Θ(n 2 ) [HNU]. We are interested in generalizations on this result to surfaces. Flips in triangulations of surfaces will be defined precisely later (Definition 2.4), for now we can just think of them as similar to edge flips in triangulations of the Euclidean plane. We emphasize that geodesics only locally minimize the length, so the edges of a geometric triangulation are generally not shortest paths. We will prove that the number of geometric triangulations on a set of points can be infinite, whereas the flip graph of "shortest path" triangulations is small but not connected in most situations [CGH + ].

Definition 1.1. Let (M 2 , h) be either a torus (T 2 , h) equipped with a Euclidean structure h or a closed oriented surface (S, h) equipped with a hyperbolic structure h. Let V ⊂ M 2 be a set of n points. The geometric flip graph F M2,h,V of (M 2 , h, V ) is the graph whose vertices are the geometric triangulations of (M 2 , h) with vertex set V and where two vertices are connected by an edge if and only if the corresponding triangulations are related by a flip.

Our results are mainly interesting in the hyperbolic setting, which is richer than the flat setting. However, to help the readers' intuition, we also present them for flat tori, where they are slightly simpler to prove and might even be considered as folklore. The geometric flip graph is known to be connected for the special case of flat surfaces with conical singularities and triangulations whose vertices are these singularities [Tah].

The main results of this paper are:

• The geometric flip graph of (M 2 , h, V ) is connected (Theorems 4.4 and 4.6).

• The Delaunay triangulation can be reached from any geometric triangulation by a path in the geometric flip graph F M2,h,V whose length is bounded by n 2 times a quantity measuring the quality of the input triangulation (Theorems 5.1 and 5.4). If an initial triangulation of the surface only having one vertex is given, then the Delaunay triangulation can thus be computed incrementally by inserting points one by one in a very standard way: for each new point, the triangle containing it is split into three, then the Delaunay property is restored by propagating flips. This approach, based on flips, can handle triangulations of a surface with loops and multiarcs, which is not the case for the approach based on Bowyer's incremental algorithm [CT, BTV]. The work presented here can hardly be compared with broad results on computing Delaunay triangulations on very general manifolds [BDG].

Background and notation

2.1. Surfaces. In this section, we first recall a few notions, then we illustrate them for the two classes of surfaces (flat tori and hyperbolic surfaces) that we are interested in.

Let M 2 be a closed oriented surface, i.e., a compact connected oriented 2-manifold without boundary. There is a unique simply connected surface M 2 , called the universal cover of M 2 , equipped with a projection ρ : M 2 → M 2 that is a local diffeomorphism. There is a natural action on M 2 of the fundamental group π 1 (M 2 ) of M 2 so that for all p ∈ M 2 , ρ -1 (p) is an orbit under the action of π 1 (M 2 ). We will denote as p a lift of p, i.e., one of the elements of the orbit ρ -1 (p). A fundamental domain in M 2 for the action of π 1 (M 2 ) on M 2 is a connected subset Ω of M 2 that intersects each orbit in exactly one point, or, equivalently, such that the restriction of ρ to Ω is a bijection from Ω to M 2 [Mas]. The genus g of M 2 is its number of handles. In this paper, we consider surfaces with constant curvature (0 or -1). The value of the curvature is given by Gauss-Bonnet Theorem and thus only depends on the genus: a surface of genus 0 only admits spherical structures (not considered here); a flat torus is a surface of genus 1 and admits Euclidean structures; a surface of genus 2 and above admits only hyperbolic structures (see below).

From now on, M 2 will denote either a flat torus or a closed hyperbolic surface. Flat tori. We denote by T 2 the topological torus, that is, the product T 2 = S 1 × S 1 of two copies of the circle. Flat tori are obtained by taking the quotient of the Euclidean plane by an Abelian group generated by two independent translations. There are in fact many different Euclidean structures on T 2 ; if one considers Euclidean structures up to homothety -which is sufficient for our purposes here -a Euclidean structure is uniquely determined by a vector u in the upper half-plane R × R >0 : to such a vector u is associated the Euclidean structure (T 2 , h u ) ∼ R 2 /(Ze 1 +Zu) , where e 1 = (1, 0) and u = (u x , u y ) ∈ R 2 is linearly independent from e 1 . The orbit of a point of the plane is a lattice. The area A h of the surface is |u y |. The plane R 2 , equipped with the Euclidean metric, is then isometric to the universal cover of the corresponding quotient surface. Hyperbolic surfaces. We now consider a closed oriented surface S (a compact oriented surface without boundary) of genus g ≥ 2. Such a surface does not admit any Euclidean structure, but it admits many hyperbolic structures, corresponding to metrics of constant curvature -1, locally modeled on the hyperbolic plane H 2 . Given a hyperbolic structure h on S, the surface (S, h) is isometric to the quotient H 2 /G, where G is a (non-Abelian) discrete subgroup of the isometry group P SL(2, R) of H 2 isomorphic to the fundamental group π 1 (S). The universal cover S is isometric to the hyperbolic plane H 2 .

For completeness, we recall below some properties of the hyperbolic plane.

2.2. The Poincaré disk model of the hyperbolic plane. In the Poincaré disk model [Ber], the hyperbolic plane is represented as the open unit disk D2 of R 2 . The points on the unit circle represent points at infinity. The geodesic lines consist of circular arcs contained in the disk D 2 and that are orthogonal to its boundary (Figure 1 (left)). The model is conformal, i.e., the Euclidean angles measured in the plane are equal to the hyperbolic angles. We won't need the exact expression of the hyperbolic metric here. However, the notion of hyperbolic circle is relevant to us. Three non-collinear points in the hyperbolic plane H 2 determine a circle, which is the restriction to the Poincaré disk of a Euclidean circle or line. If C is a Euclidean circle or line and φ : D 2 → D 2 is an isometry of the hyperbolic plane, then φ(C ∩ D 2 ) is still the intersection with D 2 of a Euclidean circle or a line.

A key difference with the Euclidean case is that the "circle" defined by 3 non-collinear points in H 2 is generally not compact (i.e., it is not included in the Poincaré disk). The compact circles are sets of points at constant (hyperbolic) distance from a point. Non-compact circles are either horocycles or hypercycles, i.e., connected components of the set of points at constant (hyperbolic) distance from a hyperbolic line (Figure 1 (right)) [Gar].1 Therefore, the relatively elementary tools that can be used for flat tori must be refined for hyperbolic surfaces. Still, some basic properties of circles still hold for non-compact circles. A non-compact circle splits the hyperbolic plane into two connected regions. We will call disk the region of the corresponding Euclidean disk that lies in the Poincaré disk. When a non-compact circle is determined by the three vertices of a triangle, its associated disk is convex (in the hyperbolic sense) and contains the whole triangle. Triangulations of hyperbolic spaces have been studied [BDT] and implemented in cgal in 2D [BIT]. Note that that previous work was not considering non-compact circles as circles.

2.3. Triangulations on surfaces. Let (M 2 , h) be either a torus (T 2 , h) equipped with a Euclidean structure h or a closed surface (S, h) equipped with a hyperbolic structure h. Let V ⊂ M 2 be a finite subset of points, and let T be a triangulation of M 2 with vertex set V .

Recall that given two distinct points v, w ∈ M 2 , any homotopy class of paths on M 2 with endpoints v and w contains a unique locally geodesic segment. We can recall the following simple notion of geometric triangulation.

Definition 2.1. A triangulation T on M 2 is said to be geometric for h if it can be realized with interior disjoint locally geodesic segments as edges.

If T is a triangulation of M 2 , its inverse image 2 ρ -1 (T ) is the (infinite) triangulation of M 2 with vertices, edges and faces that are connected components of the lifted images by ρ -1 of the vertices, edges and faces of T .

Definition 2.2. The diameter ∆(T ) of T is the smallest diameter of a fundamental domain that is the union of lifts of the triangles of T (with geodesic edges) in M 2 .

The diameter ∆(T ) is not smaller than the diameter of (S, h). It is unclear how to compute ∆(T ) algorithmically and the problem looks difficult. However bounds are easy to obtain: ∆(T ) is at least equal to the maximum of the diameters of the triangles of ρ -1 (T ) in M 2 and is at most the sum of the diameters of these triangles.

Definition 2.3. We say that a triangulation T of M 2 is a Delaunay triangulation if for each face f of T and any face f of ρ -1 (T ), there exists an open disk in M 2 inscribing f that is empty, i.e., that contains no vertex of ρ -1 (T ).

We will see in Section 3 that any Delaunay triangulation of M 2 is geometric. Remark that, even for a hyperbolic surface, every empty disk in the universal cover H 2 is compact. Indeed, any non-compact disk contains at least one disk of any diameter, so, at least one disk of diameter ∆(T ), thus it contains a fundamental domain (actually, infinitely many fundamental domains) and cannot be empty.

Let us now give a natural definition for flips in triangulations of surfaces. It is based on the usual notion of flips in the Euclidean plane.

Definition 2.4. Let T be a triangulation of M 2 . Let (v 1 , v 2 , v 3 ) and (v 2 , v 1 , v 4 ) be two adjacent triangles in T , sharing the edge e = (v 1 , v 2 ). Let us lift the quadrilateral (v 1 , v 2 , v 3 , v 4 ) to M 2 so that ( v 1 , v 2 , v 3 ) and ( v 2 , v 1 , v 4 ) form two adjacent triangles of ρ -1 (T ) sharing the edge e = ( v 1 , v 2 ).
Flipping e in T consists in replacing the diagonal e in the quadrilateral

( v 1 , v 2 , v 3 , v 4 ) (which lies in M 2 , i.e., R 2 or H 2 ) by the other diagonal ( v 3 , v 4 ), then projecting the two new triangles ( v 3 , v 4 , v 2 ) and ( v 4 , v 3 , v 1 ) to M 2 by ρ.
We say that the flip of T along e is Delaunay if the triangulation is locally Delaunay in the quadrilateral after the flip, i.e., the disk inscribing

( v 3 , v 4 , v 2 ) does not contain v 1 (and the disk inscribing ( v 4 , v 3 , v 1 ) does not contain v 2 ).
An edge e is said to be Delaunay flippable if the flip along e is Delaunay.

Note that even if T is geometric, the triangulation after a flip is not necessarily geometric. We will prove later (Lemma 4.1) that a Delaunay flip transforms a geometric triangulation into a geometric triangulation. Degenerate sets of points on a surface. Let us quickly examine here the case of degenerate sets of points, i.e., sets of points V on M 2 such that the infinite Delaunay triangulation of ρ -1 (V ) is not unique, i.e., at least two adjacent triangles in the possible Delaunay triangulations of ρ -1 (V ) in M 2 have cocircular vertices. In such a case, any triangulation of the subset C of ρ -1 (V ) consisting of c cocircular points is a Delaunay triangulation. Any of these triangulations can be transformed in any other by O(c) flips [HNU]. From now on, we can thus assume that the set of points V on the surfaces that we consider is always non-degenerate. Triangulations and polyhedral surfaces. The Euclidean plane can be identified with the plane (z = 1) in R 3 , while the Poincaré model of the hyperbolic plane can be identified with the unit disk in that plane. We can now use the stereographic projection σ : S 2 \ {s 0 } → R 2 to send the unit sphere S 2 to this plane (z = 1), where s 0 = (0, 0, -1) is the pole. In this projection, each point p = s 0 on the sphere is sent to the unique intersection with the plane (z = 1) of the line going through s 0 and p. The inverse image of the plane (z = 1) is S 2 \ {s 0 }, while the inverse image of the disk containing the Poincaré model of the hyperbolic plane is a disk, which is the set of points of S 2 above a horizontal plane.

Let T be a triangulation of the Euclidean or the hyperbolic plane -for instance, T could be the inverse image ρ -1 (T ) of a triangulation T of a surface (M 2 , h), in which case T has infinitely many vertices. We associate to T a polyhedral surface Σ in R 3 , constructed as follows. The construction is similar to the classic duality originally presented with a paraboloid in the case of (finite) triangulations in a Euclidean space [ES]. It can also be seen as a simpler version, sufficient for our purpose, of the construction presented for triangulations in hyperbolic spaces using the space of spheres [BDT].

• The vertices of Σ are the inverse images on § 2 by σ of the vertices of T .

• The edges of Σ are line segments in R 3 corresponding to the edges of T and the faces of Σ are triangles in R 3 corresponding to the faces of T . Note that Σ is not necessarily convex. We can make the following well-known remarks. Let t 1 and t 2 be two triangles of T sharing an edge e, and let t Σ 1 and t Σ 2 be corresponding faces of the polyhedral surface Σ, sharing the edge e Σ . Then Σ is concave at e Σ if and only if e is Delaunay flippable. Flipping e in the triangulation T in the plane corresponds to replacing the two faces t Σ 1 and t Σ 2 of Σ by the two other faces of the tetrahedron formed by their vertices. That tetrahedron lies between Σ and S 2 . We obtain a new edge e Σ at which the new polyhedral surface Σ is convex, and which is strictly closer to S 2 than Σ. By an abuse of language, we will say that Σ contains Σ, which we will denote as Σ ⊂ Σ .

As a consequence, Σ is convex if and only if T is Delaunay. There is a direct corollary of this statement: Given a (non-degenerate, see above) discrete set V of points in R 2 or H 2 , there is a unique Delaunay triangulation with this set of vertices.

However we are going to see in the next two sections that there can be infinitely many geometric (non-Delaunay) triangulations on a surface, with the same given finite vertex set.

Geometric triangulations of surfaces

We consider now Dehn twists, which are usually considered as acting on the space of metrics on a surface [CB], but are defined here equivalently, for simplicity, as acting on triangulations of a closed oriented surface (M 2 , h) equipped with a fixed Euclidean or hyperbolic structure (figures in this section illustrate the flat case, but the results are proved for both flat and hyperbolic cases). Let T be a triangulation of (M 2 , h), with vertex set V , and let c be an oriented homotopically nontrivial simple closed curve on M 2 \ V . We define a new triangulation τ c (T ) of M 2 by performing a Dehn twist along c: whenever an edge e of T intersects c at a point p, we orient e so that the unit vectors of the tangent plane along e and c form a positively oriented basis (see Figure 2 (left)), and then replace e by the oriented path following e until p, then following c until it comes back to p, then following e until its endpoint (see Figure 2 (right)). This defines a map τ c from the space of triangulations of T 2 with vertex set V to itself. Note that, even if T is a geometric triangulation, τ c (T ) is not necessarily geometric. If we denote by -c the curve c with the opposite orientation, then one easily checks that τ -c = τ Lemma 3.1. There exists a geometric triangulation T of (M 2 , h) and a simple closed curve c ⊂ M 2 such that for all k ∈ Z, τ k c (T ) is geometric. Proof. We choose a simple closed geodesic c on (M 2 , h) and ε > 0. We denote by c -, c + the two geodesics at distance ε from c on the positive and negative sides of c. The value of ε must be sufficiently small so that the region between c -and c + is an annulus drawn on M 2 . We then choose a geometric triangulation T of (M 2 , h) with no vertex in the open annulus bounded by c -and c + and containing c, such that each edge crossing c intersects c exactly once, and has one endpoint on c -and another on c + .

We realize the image by τ c of an edge e of T as a geodesic segment -there is a unique choice in the homotopy class of the path described above (Figure 3 Corollary 3.2. For any closed oriented surface (M 2 , h), there exists a finite set of points V ⊂ M 2 such that the graph of geometric triangulations with vertex set V is infinite.

We can now prove the following result:

Proposition 3.3. Any Delaunay triangulation of a closed oriented surface (M 2 , h) is geometric.

Proof. Let V be a finite set of points on M 2 , and let T be the Delaunay triangulation of (M 2 , h) with vertex set V . Realize every edge of T as a the unique geodesic segment in its homotopy class. We argue by contradiction and suppose that T is not geometric, so that there are two edges e 1 and e 2 that intersect in their interiors. We then lift e 1 and e 2 to edges e 1 and e 2 of ρ -1 (T ) whose interiors still intersect at one point.

There are at least two distincts faces f 1 and f 2 of ρ -1 (T ) such that e 1 is an edge of f 1 and e 2 is an edge of f 2 . Let C 1 and C 2 be the circles inscribing f 1 and f 2 , respectively. Since ρ -1 (T ) is Delaunay, C 1 and C 2 bound empty disks D 1 and D 2 , i.e., open disks not containing any point of ρ -1 (V ). Recall that, as mentioned in Section 2.3, empty disks are compact even in the hyperbolic case, and that e 1 ⊂ D 1 and e 2 ⊂ D 2 (edges are considered as open).

The two circles C 1 and C 2 do not intersect more than twice. Let L be the geodesic line through their two intersection points. The endpoints of e 1 are on C 1 \ D 2 and those of e 2 are on C 2 \ D 1 , so the two pairs of endpoints are on opposite sides of L. As a consequence, e 1 and e 2 are on opposite sides of L, so they cannot intersect. This leads to a contradiction.

The flip algorithm

Let us consider a closed oriented surface (M 2 , h). The flip algorithm consists in performing Delaunay flips in any order, starting from a given input geometric triangulation of M 2 , until there is no more Delaunay flippable edge.

In this section, we first define a data structure that supports this algorithm, then we prove the correctness of the algorithm. 4.1. Data structure. In both cases of a flat or hyperbolic surface, the group of isometries defining the surface is denoted as G. We assume that a fundamental domain Ω 0 is given. By definition (Section 2.1), M 2 is the union G(Ω 0 ) of the images of Ω 0 under the action of G.

To represent a triangulation on the surface, we propose a data structure generalizing the data structure previously introduced for triangulations of flat orbifolds [CT] and triangulations of the Bolza surface [IT]. The combinatorics of the triangulation is given by the set of its vertices V on the surface and the set of its triangles, where each triangle gives access to its three vertices in V and its three adjacent triangles, and each vertex gives access to one of its incident triangles. The geometry of the triangulation is given by the set V 0 of the lifts of its vertices that lie in the fundamental domain Ω 0 and one lift t 0 in M 2 of each triangle t = (v 0,t ; v 1,t ; v 2,t ) of the triangulation, chosen among the (one, two, or three) lifts of t in M 2 having at least one vertex in Ω 0 : t 0 has at least one of its vertices v i,t 0 in Ω 0 (i = 0, 1, or 2); then the other vertices of t 0 are images g i+1,t • v i+1,t 0 and g i+2,t • v i+2,t 0 of two vertices in V 0 , where g i+1,t and g i+2,t are elements of G (indices are taken modulo 3). In the data structure, each vertex v on the surface has access to its representative v 0 , and each triangle t on the surface has access to the isometries g 0,t , g 1,t , and g 2,t allowing to construct t 0 , at least one of the isometries being the identity 1 G . Note that two triangles t and t of T that are adjacent on the surface are represented by two triangles t 0 and t 0 , which are not necessarily adjacent in M 2 (Figure 5 (left)). However, there is an isometry g in G such that t 0 and g • t 0 are adjacent.

Let T be an input triangulation given as such a data structure. Figure 5 illustrates a Delaunay flip performed on two adjacent triangles t and t on the surface. The triangle t 0 is first moved so that the vertices of the edge to be flipped coincide. Then the edge is flipped. The isometries in the two triangles created by the flip are easy to compute from the isometries stored in t and t . Note that the order in which isometries are composed is crucial in the hyperbolic case, as they do not commute. We have shown that the data structure can be maintained through flips. Right: the isometries in the two triangles created by the flip.

Ω 0 t 0 g i = 1 Γ g i+1 g i+2 g i+1 g -1 j • t 0 g i+1 g i+2 g j g j+1 g j+2 1 Γ 1 Γ g j+1 g

4.2.

Correctness of the algorithm. The following statement is a key starting point.

Lemma 4.1. Let T be a geometric triangulation of (M 2 , h), and let T be obtained from T by a Delaunay flip. Then T is still geometric.

Proof. Let e be a Delaunay flippable edge and e a lift in M 2 . Denote the vertices of e by v and v . Let t 1 and t 2 be the triangles of ρ -1 (T ) incident to e. To prove that T is geometric, it is sufficient to prove that t 1 ∪ t 2 is a strictly convex quadrilateral.

Let C 1 (resp. C 2 ) be the circle through the three vertices of t 1 (resp. t 2 ). Note that C 1 and C 2 may be non-compact. Let D 1 and D 2 be the corresponding disks (as defined in Section 2.2 on case of non-compact circles). The disk D 1 (resp. D 2 ) is convex (in the Euclidean plane if M 2 is a flat torus, or in the sense of hyperbolic geometry if M 2 is a hyperbolic surface) and contains t 1 (resp. t 2 ). The fact that e is Delaunay flippable then implies that t 1 and t 2 are contained in D 1 ∩ D 2 (see Figure 6). As a consequence, the sum of angles of t 1 and t 2 at v is smaller than the interior angle at v of D 1 ∩ D 2 , which is at most π, and similarly at v . As a consequence, the quadrilateral t 1 ∪ t 2 is strictly convex at v and v . Since it is strictly convex at its other two vertices (as each of these vertices is a vertex of a triangle), it is strictly convex, and the statement follows.

The following lemma, using the diameter of the triangulation (Definition 2.2), is central in the proof of the termination of the algorithm (Theorem 4.6) for hyperbolic surfaces and in its analysis for both flat tori and hyperbolic surfaces (Section 5).

Lemma 4.2. Let T be a geometric triangulation of (M 2 , h). Then, the flip algorithm starting from T will never insert an edge longer than 2∆(T ).

Note that the length of an edge can be measured on any or its lifts in the universal covering space M 2 .

Proof. Let T k be the triangulation obtained from T = T 0 after k flips and let Σ k be the corresponding polyhedral surface of R 3 as defined in Section 2.3. Since we perform only Delaunay flips, Σ 0 ⊂ . . . ⊂ Σ k ⊂ Σ k+1 (with the abuse of language mentioned in Section 2.3).

We will prove the result by contradiction. Let us assume that T k has an edge e of length larger than 2∆(T ). Let Ω be a fundamental domain of M 2 having diameter ∆(T ), given as the union of lifts of triangles of T = T 0 (it is not clear how to compute such a fundamental domain efficiently but its existence is clear). Let v be the midpoint of e and v its lift in Ω. Let e = ( v 1 , v 2 ) be the unique lift of e whose midpoint is v. The domain Ω is strictly included in the disk D of radius ∆(T ) and centered at v, by definition of ∆(T ) (see Figure 7 (left)).

Let P D denote the plane in R 3 containing the circle on S 2 that is the boundary of σ -1 (D) (recall that σ denotes the stereographic projection, see Section 2.3), and let p denote the point σ -1 ( v) on S 2 . As p ∈ σ -1 (Ω) ⊂ σ -1 ( D), the projection p Σ0 of p onto Σ 0 lies above P D (Figure 7 (right)).

S 2 s 0 R 2 D v 1 v 2 p 1 p 2 D v p Σ k p p Σ0 Ω D v v 2 v 1 D Ω e P D Figure 7.
Illustration for the proof of Lemma 4.2 (for a hyperbolic surface). Left: notation in H 2 . Right: contradiction seen in a cutting plane in R 3 . Now, denote the edge σ -1 ( e) on S 2 as (p 1 , p 2 The points p 1 and p 2 lie outside σ -1 (D). So, the corresponding edge e Σ = [p 1 , p 2 ] of Σ k lies below the plane P C , thus the projection p

Σ k ∈ [p 1 , p 2 ] of p onto Σ k lies below P C .
From what we have shown, p Σ k is a point of Σ k that lies strictly between the pole s 0 and the point p Σ0 of Σ 0 , which contradicts the inclusion Σ 0 ⊂ Σ k .

We will now show that, for any order, the flip algorithm terminates and returns the Delaunay triangulation of the surface. The proof given for the hyperbolic case would also work for the flat case. However we propose a more elementary proof for the flat case.

Flat tori. The case of flat tori is easy, and might be considered as folklore. However, as we have not found a reference, we give the details here for completeness.

We define the weight of a triangle t of a geometric triangulation T of T 2 as the number of vertices of ρ -1 (T ) that lie in the open circumdisk of a lift of t. The weight w(T ) of T is defined as the sum of the weights of its triangles.

Lemma 4.3. The weight w(T ) of a triangulation T of a flat torus (T 2 , h) is finite. Let T be the triangulation obtained from a geometric triangulation T after performing a Delaunay flip. Then w(T ) ≤ w(T ) -2.

Proof. A circumdisk of any triangle in R 2 is compact, so, it can only contain a finite number of vertices of ρ -1 (T ). The sum w(T ) of these numbers over triangles of T is clearly finite as the number of triangles of T is finite. Let us now focus on a quadrilateral in R 2 that is a lift of the quadrilateral on T 2 whose diagonal e is flipped. Let D 1 and D 2 denote the two open circumdisks in R 2 before the flip and D 1 and D 2 denote the two open circumdisks after the flip, then Hyperbolic surfaces. To show that the flip algorith terminates in the hyperbolic case, we cannot mimic the proof presented for the flat tori since the circumcircle of a hyperbolic triangle can be non-compact (see Section 2.2) and thus can have an infinite weight. Note also that the proof cannot use a property on the angles of the Delaunay triangulation similar to what holds in the Euclidean case: in H 2 , the locus of points seeing a segment with a given angle is not a circle arc, and thus the Delaunay triangulation of a set of points in H 2 does not maximize the smallest angle of triangles. The proof relies on Lemma 4.2.

D 1 ∪ D 2 ⊂ D 1 ∪ D 2 and D 1 ∩ D 2 ⊂ D 1 ∩ D 2 (see Figure 8
Theorem 4.6. Let T be a geometric triangulation of a closed hyperbolic surface with finite vertex set V . The flip algorithm terminates and outputs the Delaunay triangulation of V .

Proof. We use the same notation as in the proof Lemma 4.2. Once an edge of T k is flipped, it can never reappear in the triangulation, as the corresponding segment in R 3 becomes interior to the polyhedral surface Σ k+1 (see Section 2.3) and further surfaces Σ k , k ≥ k + 1. In addition, all the introduced edges have length smaller than 2∆(T ) by Lemma 4.2. Moreover, there is only a finite number of edges with vertices in V that are shorter than 2∆(T ) on S, as a circle given by a center and a bounded radius is compact. So, the flip algorithm terminates. The output does not have any Delaunay flippable edge, so, it is the Delaunay triangulation.

Corollary 4.7. The geometric flip graph F S,h,V is connected.

Algorithm analysis

For a triangulation on n vertices in the Euclidean plane, counting the weights of triangulations leads to the optimal O(n 2 ) bound. However the same argument does not yield a bound even for the flat torus, since points must be counted in the universal cover.

Theorem 5.1. For any triangulation T with n vertices of a torus (T 2 , h), there is a sequence of flips of length C h • ∆(T ) 2 • n 2 connecting T to a Delaunay triangulation of (T 2 , h), where C h only depends on h.

Proof. Let e = (v 1 , v 2 ) be an edge appearing during the flip algorithm, and v 1 (resp. v 2 ) be a lift of v 1 (resp. v 2 ), such that ( v 1 , v 2 ) is a lift e of e. The point v 2 lies in a circle C of diameter 4∆(T ) centered at v 1 by Lemma 4.2. Let M be the affine transformation that maps the lattice of the lifts of v 2 to the square lattice Z 2 . M (C) is a convex set and from Pick's theorem [Tra],3 the number of points of Z 2 in M (C) is smaller than area(M (C)) + 1/2 • perimeter(M (C)) + 1, which is also a bound on the number of possible points v 2 in C and thus the number of possible edges e. The area of M (C) is 1/A h • area(C) since det(M ) = 1/A h , but there is no simple formula for its perimeter. As already mentioned in the proof of Theorem 4.6, an edge can never reappear after it was flipped. Moreover, there are n 2 /2 sets of points {v 1 , v 2 } (v 1 and v 2 may be the same point), which yields the result.

The rest of this section is devoted to computing the number of edges not longer than 2∆(T ) between two fixed points v 1 and v 2 on a hyperbolic surface (S, h). Counting the number of points in a disk of fixed radius would give an exponential bound because the area of a circle in H 2 is exponential in its radius [Mar]. Note that we only consider geodesic edges, so we only need to count homotopy classes of simple paths. The behavior of the number N l of simple closed curves smaller than a fixed length l is well understood: N l /l 6g-6 converges to a positive constant depending "continuously" on h [Mir]. However, we need a result for geodesic paths instead of geodesic closed curves, and Mirzakhani's proof is too deep and relies on too sophisticated structures to easily be generalized. So, we will only prove an upper bound on the number of paths. Such an upper bound could be derived from the theory of measured laminations of Thurston, which is also quite intricate. Fortunately, a more comprehensible proof, specific to simple closed geodesic curves on hyperbolic structures, can be found in a book published by the French Mathematical Society [FLP + , 4.III, p.61-67] [FLP]. While recalling the main steps of the proof, we show how to extend it to geodesic paths.

Let Γ = {γ, i = 1, . . . , 3g -3} be a set of 3g -3 simple disjoint closed geodesics on (S, h) not containing v 1 and v 2 that forms a pants decomposition on S, where each γ i belongs to two different pairs of pants. A set {γ i , i = 1, . . . , 3g -3} of disjoint closed annuli is defined on S, where each γ i is a tubular neighborhood of γ i containing none of v 1 , v 2 . This yields a decomposition of S into 3g -3 annuli γ i (i = 1, . . . , 3g -3) and 2g -2 pairs of "short pants" P j (j = 1, . . . , 2g -2). For i = 1, . . . , 3g -3, let us denote as ∂γ i any one of the two curves bounding the annulus γ i (this is an abuse of notation but should not introduce any confusion). In each pair of pants P j , j = 1, . . . , 2g -2, for each boundary ∂γ, an arc J γ i is drawn in P i , going from the boundary of γ to itself that separates the other two boundaries of P i and that has minimal length.

Two curves γ and γ are associated to each γ ∈ Γ in the following way (Figure 9). The annulus γ is glued with the two pairs of pants P i and P j between which it is lying, which yields a sphere with four boundaries: ∂γ i,1 and ∂γ i,2 bounding P i and ∂γ j,1 and ∂γ j,2 bounding P j . A curve γ is then defined: it coincides with J γ i in P i and J γ j in P j , it separates ∂γ i,1 and ∂γ j,1 from ∂γ i,2 and ∂γ j,2 , and it has exactly 2 crossings with γ. The curve γ is defined in the same way, separating ∂γ i,1 and ∂γ j,2 from ∂γ i,2 and ∂γ j,1 .

For each P i and m i,1 , m i,2 , m i,3 ∈ N, a model multiarc is fixed in P i , having m i,1 , m i,2 and m i,3 intersections with the three boundaries ∂γ i,1 , ∂γ i,2 and ∂γ i,3 of P i (if one exists). The model is chosen among all the possible model multiarcs as the one that has a minimal number of intersections with the three arcs J γi,j i (j = 1, 2, 3) of P i . The model multiarcs is unique, up to homeomorphisms of the pair of pants, and those homeomorphisms are rather simple to understand since they can be decomposed into three Dehn twists around curves homotopic to the three boundaries of the pair of pants.

Let now f be a path between v 1 and v 2 on S. We decompose f into three parts: (v 1 , w 1 ), (w 1 , w 2 ) and (w 2 , v 2 ) where w 1 and w 2 are the first and the last point of f on an annulus boundary. We "push" all the twists of f w into the annuli γ, γ ∈ Γ, and obtain a normal form homotopic to f , whose definition adapts the definition given in the book [FLP + ] for closed curves:

(1) It is simple.

(2) It has a minimal number m i of intersections with each γ i , i = 1, . . . , 3g -3.

(3) In each P j , j = 1, . . . , 2g -2, it is homotopic with fixed endpoints to the model multiarc that corresponds to the number of intersections with its boundaries. For P j1 (resp. P j2 ) containing v 1 (resp. v 2 ), only the intersections different from w 1 (resp. w 2 ) are counted. (4) Between v 1 and w 1 (resp. w 2 and v 2 ), it has a minimal number of intersections with the three arcs J

γ j 1 ,k j1 (k = 1, 2, 3) in P j1 containing v 1 (resp. J γ j 2 ,k j2
in P j2 containing v 2 ). ( 5) It has a minimal number t i of intersections with γ i inside γ i , for any i = 1, . . . , 3g -3. (6) It has a minimal number s i of intersections with γ i inside γ i , for any i = 1, . . . , 3g -3. The existence of a normal form is clear but its uniqueness is unclear (uniqueness is not required for the upper bound that we are looking for, but it can actually be proved by extension of the next lemma). The two forms of the path f are used to define two notions of complexity: its geodesic form is used to define its length, which can be seen as a geometric complexity, whereas its its normal coordinates m i , s i and t i can be seen as a combinatorial complexity. Lemma 5.3 shows some equivalence between the two notions of complexity. We first show that a fixed set of coordinates corresponds to a finite number of possible non-homotopic paths.

Lemma 5.2. For any set of coordinates m i , t i , s i , i = 1, . . . , 3g-3, there are at most 9(max {i=1,...,3g-3} (m i )) 2 non-homotopic normal forms.

Proof. Let f be a path, decomposed as above into (v 1 , w 1 ), f w = (w 1 , w 2 ) and (w 2 , v 2 ). The uniqueness for closed curves comes from the facts that in each pair of pants, fixing the m i , s i and t i leads to a unique homotopy class of model multiarcs [FLP + , Lemma 5, p.63]. Everything remains true but the uniqueness of the homotopy class of model multiarcs in the two (not necessarily different) pairs of pants P j1 and P j2 containing v 1 and v 2 . However, w 1 and w 2 are fixing unique models (see Figure 10). There are three possible annulus boundaries ∂γ j,i , i = 1, 2, 3 for w 1 in the pair of pants P j that contains v 1 (resp. γ j,i for w 2 ), so, at most 3 max {i} (m i ) possibilities for each of them. The choices for w 1 and w 2 are independent and the result follows.

Lemma 5.3. Let f be a geodesic path of length l, then there exists a constant c h such that the coordinates m i , t i and s i , i = 1, . . . , 3g -3 of the normal form of f are smaller than c h • l.

Proof. For any simple closed geodesic δ on S, the geodesic form of f intersects δ in a minimal number k δ of points, since they are both geodesics. If ε δ is the width of a tubular neighborhood of δ, then l ≥ ε δ (k δ -1) [START_REF] Birman | Geodesics with bounded intersection number on surfaces are sparsely distributed[END_REF]Lemma 3.1]. Each coordinate m i , t i and s i of f corresponds to the minimal number of intersections with a curve. The number m i corresponds to γ i . The number t i is actually not larger than the number of intersections of f with the geodesic curve that is homotopic to γ i (γ i is generally not geodesic), and similarly s i is not larger than the number of intersections of f with the geodesic homotopic to γ i . These curves γ i , γ i , γ i only depend on (S, h), so, we can take ε h to be the largest of all the 9g -9 widths ε γi , ε γ i , ε γ i and we obtain l ≥ ε h • max(m i , t i , s i ) and thus max(m i , t i , s i ) ≤ 1/ε h • l.

Theorem 5.4. For any hyperbolic structure h on S and any triangulation T of (S, h), there is a sequence of flips of length at most C h • ∆(T ) 6g-4 • n 2 in the geometric flip graph connecting T to a Delaunay triangulation of (S, h).

Proof. Let N v1,v2 be the number of paths from v 1 to v 2 shorter than l = 2 • ∆(T ). From the previous lemma, we obtain that the 9g -9 coordinates m i , t i , and s i of any such path f are smaller than c h • 2∆(T ). It appears that, ∀i, m 1 = t i + s i , t i = m i + s i or s i = m i + t i [FLP + , Lemma 6, p.64]. So, if we fix m i and t i there are at most 3 possible s i . Lemma 5.2 and 5.3 proves that there are 9(c h • 2∆(T )) 2 potential paths for each coordinate set. We obtain a bound for N v1,v2 : N v1,v2 ≤ 9(c h • 2∆(T )) 2 • 3(c h • 2∆(T )) 6g-6 and thus, there is a constant C h such that N v1,v2 ≤ C h • ∆(T ) 6g-4 . Since there are 1/2 • n 2 possible sets {v 1 , v 2 }, we obtain the bound on the number of edges.

Figure 1 .

 1 Figure 1. The Poincaré disk. Left: Geodesic lines (black) and compact circles (red) centered at point ω. Right: A horocycle (green). A hypercycle (blue), whose points have constant distance from the black geodesic line.

Figure 2 .

 2 Figure 2. Transformation of an edge e by the Dehn twist along c on a flat torus T 2 . Here the black parallelepiped is a fundamental domain, and the gray one, used for the construction of the image of e by τ c , is another fundamental domain, image through an element of the the group Γ of isometries.

Figure 3 .Figure 4 .

 34 Figure 3. Image of e by a Dehn twist (middle), realized as a geodesic edge (right).

Figure 5 .

 5 Figure 5. A flip. Here (hyperbolic) triangles are represented schematically with straight edges. Left: the two triangles t 0 and t 0 before the flip. Here g i = 1 G .

Figure 6 .

 6 Figure 6. The quadrilateral is convex (edges are represented schematicaly as straight line segments).

Figure 8 .

 8 Figure 8. Circumdisks D 1 and D 2 before flipping e and D 1 and D 2 after the Delaunay flip.

Figure 9 .

 9 Figure9. Two adjacent pairs of pants P i and P j .

Figure 10 .

 10 Figure 10. Three possible choices for w 1 . The two left choices correspond to the same model multiarcs, but the orderings on the upper boundary lead to nonhomotopic paths. The right choice leads to different model multiarcs.

A synthetic presentation can be found at http://en.wikipedia.org/wiki/Hypercycle_(geometry)

the notion of pull-back would be more correct but we stay with inverse image for simplicity

See also https://en.wikipedia.org/wiki/Pick's_theorem#Inequality_for_convex_sets

The authors were supported by the grant(s) ANR-17-CE40-0033 of the French National Research Agency ANR (project SoS) and INTER/ANR/16/11554412/SoS of the Luxembourg National Research fund FNR (https:// members.loria.fr/Monique.Teillaud/collab/SoS/).

1


Awatef Ouhibi

Maroua Saadaoui

Nathalie Lorrain
email: nathalie.lorrain@univ-rennes1.fr

Mohammed Guendouz

Noureddine Raouafi

Adel Moadhen
email: adel.moadhen@istmt.utm.tn

Application of Doehlert Matrix for an Optimized Preparation of a Surface-Enhanced Raman Spectroscopy (SERS) Substrate Based on Silicon Nanowires for Ultrasensitive Detection of Rhodamine 6G

Keywords: Silicon nanowires, silver nanoparticles, surface-enhanced Raman spectroscopy substrate, SERS substrate, Doehlert matrix, detection

Introduction

Fast results analysis, cost-effectiveness, and simultaneous study of several parameters are the main requirements for analytical methodologies optimization and development. Surfaceenhanced Raman spectroscopy (SERS) is well recognized as a powerful tool to analyze chemical and biological samples due to the wealth of information it provides from the vibrational DOI: 10.1177/0003702819881222 2 fingerprinting of molecules. [START_REF] Ru | Principles of Surface-Enhanced Raman Spectroscopy and Related Plasmonic Effects[END_REF] SERS effect takes place when the target is adsorbed on the metal surface. [START_REF] Ru | Principles of Surface-Enhanced Raman Spectroscopy and Related Plasmonic Effects[END_REF][START_REF] Fleischmann | Raman Spectra of Pyridine Adsorbed at a Silver Electrode[END_REF] Using optimal parameters, Raman signals can be increased by up to fifteen orders of magnitude to ensure the ability to detect very low concentration of molecules, or even the detection of a single molecule. [START_REF] Albrecht | Anomalously Intense Raman Spectra of Pyridine at a Silver Electrode[END_REF][START_REF] Jeanmaire | Blu Excited State of the Tetracyanoquinodimethane Anion Radical[END_REF] Electromagnetic and chemical effects are recognized as the most contributing factors to Raman signal enhancement. [START_REF] Ru | Principles of Surface-Enhanced Raman Spectroscopy and Related Plasmonic Effects[END_REF] These effects are respectively due to the increase of electromagnetic fields by the resonant excitation of the plasmon surface via the photon interaction and to the charge-transfer gaps between adsorbed target and metal. [START_REF] Wang | Surface Plasmon Resonance Sensors on Raman and Fluorescence Spectroscopy[END_REF][START_REF] Zou | Silver Nanoparticle Array Structures that Produce Giant Enhancements in Electromagnetic Fields[END_REF] Silver and gold nanoparticles are the two most often used SERS substrates due to their outstanding optical properties, their high enhancement factors and their low reactivity.

The large surface-to-volume ratio and the one-dimensional structure of silicon nanowires (SiNWs) dramatically modify the electrical and optical properties of the bulk silicon, [START_REF] Li | Structure and Antireflection Properties of SiNWs Arrays Form mc-Si Wafer Through Ag-Catalyzed Chemical Etching[END_REF][START_REF] Zhang | Fabrication of Ultra-Low Antireflection SiNWs Arrays from mc-Si Using One Step MACE[END_REF][START_REF] Han | Metal-Assisted Chemical Etching of Silicon and Nanotechnology Applications[END_REF][START_REF] Gonchar | Optical Properties of Silicon Nanowires Fabricated by Environment-Friendly Chemistry[END_REF][START_REF] Hutagalung | Optical and Electrical Characteristics of Silicon Nanowires Prepared by Electroless Etching[END_REF] which makes them suitable candidates for numerous applications such as diagnosis of diseases, [START_REF] Zhang | Silicon Nanowire Biosensor and its Applications in Disease Diagnostics: A Review[END_REF] sensing, [START_REF] Rahman | Enhanced Sensing of Dengue Virus DNA Detection Using O2 Plasma Treated-Silicon Nanowire Based Electrical Biosensor[END_REF][START_REF] Cao | A Facile Fluorescent Sensor Based on Silicon Nanowires for Dithionite[END_REF][START_REF] Qin | Ultrasensitive Silicon Nanowire Sensor Developed by a Special Ag Modification Process for Rapid NH3 Detection[END_REF] catalysis, [START_REF] Wang | Catalase-Like and Peroxidase-Like Catalytic Activities of Silicon Nanowire Arrays[END_REF] manufacture of batteries, [START_REF] Cohn | New Insight into the Discharge Mechanism of Silicon-Air Batteries Using Electrochemical Impedance Spectroscopy[END_REF] and micro-super capacitors and solar cells. [START_REF] Soam | Controlling the Shell Microstructure in a Low-Temperature-Grown SiNWs and Correlating it to the Performance of the SiNWs-Based Micro-Supercapacitor[END_REF][START_REF] Chen | Efficiency Enhancement of Silicon Nanowire Solar Cells by Using UV/Ozone Treatments and Micro-Grid Electrodes[END_REF] Over the past two decades, many strategies have been developed to combine a hierarchical nano-array effect of silicon nanowires decorated with a metallic surface of silver nanoparticles (AgNPs) in order to design high-performance SERS substrates. [START_REF] Gebavi | Silicon Nanowires Substrates Fabrication for Ultra-Sensitive Surface Enhanced Raman Spectroscopy Sensors[END_REF][21][START_REF] Shao | Ag-Modified Silicon Nanowires Substrate for Ultrasensitive Surface-Enhanced Raman Spectroscopy[END_REF] The optimization of SERS substrates by the traditional univariate approach wherein only one parameter is changed while the others are fixed requires ongoing testing and much time for results interpretation. Moreover, this traditional method is less effective. For these reasons, the choice of the univariate approach for optimization and development of analytical methodologies is being replaced by multivariate approaches for their advantages of identifying the different intertwining of variables, avoiding time-consuming experiments and obtaining an efficient data analysis. [START_REF] Zhang | A High-Efficiency Surface-Enhanced Raman Scattering Substrate Based on Silicon Nanowires Array Decorated with Silver Nanoparticles[END_REF] In this work, we combine the effect of an etched surface area of silicon with the surface Plasmon of AgNPs to optimize SERS substrate preparation for ultrasensitive detection of Rhodamine 6G (R6G) by varying silver nitrate concentration, silver nitrate and R6G immersion times. The interdependence of these different parameters dictated the use of a multivariate method allowing simultaneous modification of the mentioned above three parameters.

The Doehlert plan was used first to plan the experiment, then to determine the optimal conditions by studying the effect of the experimental parameters on the R6G SERS response and the interactions between them. [START_REF] Doehlert | Uniform Shell Design[END_REF][START_REF] Pinkowska | Application of the Doehlert Matrix for the Determination of the Optimal Conditions of Hydrolysis of Soybean Protein in Subcritical Water[END_REF][START_REF] Hannachi | Preparation of Manganese Sulfide (MnS) Thin Films by Chemical Bath Deposition: Application of the Experimental Design Methodology[END_REF][START_REF] Bigan | Cold Remote Plasma Modification of Wood: Optimization Process Using Experimental Design[END_REF]

Experimental

Reagents and Materials

Silver nitrate, hydrofluoric acid (40%), hydrogen peroxide (34%), nitric acid (70%), sulfuric acid (97%), and Rhodamine 6G from Sigma-Aldrich (Germany) were used for all the experiments. All solutions were prepared using ultrapure water produced by a Millipore system.

The morphological properties of the as obtained SiNWs and AgNPs/SiNWs surface were studied by scanning electron microscopy (SEM) using HITACHI SU3500 microscope.

The SERS measurements were performed at room temperature (RT) on a Micro Raman HORIBA system (LabRAM HR800). A helium-neon laser with a power of about 1 mW was used for excitation at 632.8 nm of the sample surface. Data were collected with an integration time of ten seconds.

The construction of the design experimental methodology, the modeling of the SERS responses and the drawing of the various representations were carried out using NemrodW software (v.2000, LPRAI, France).

Elaboration of SERS Substrates

The SiNWs were prepared using metal-assisted chemical etching (MACE) of n-Si (100) wafers. [START_REF] Zhang | Preparation of Large-Area Uniform Silicon Nanowires Arrays Through Metal-Assisted Chemical Etching[END_REF] Briefly, wafers with a surface of 1 cm 2 were cleaned in a acetone and ethanol solution for 15 minutes. Then, they were immersed into a solution of H2O2 and H2SO4 in a ratio 1:3 [v/v] for 10 minutes. After that, the wafers were immersed in a 5% HF aqueous solution for three minutes to obtain H-terminated silicon surfaces. The wafers were immediately immersed in a aqueous solution containing 5.10 -3 M of Ag NO3 and 4.8 M of HF for one minute in order to form AgNPs before being thoroughly washed with ultrapure water. Subsequently, the as obtained samples were immersed in an aqueous etching solution composed of 4.8 M HF and 0.4 M H2O2 for 45 minutes to prepare SiNWs. Moreover the samples were immersed in aqueous solution of HNO3 for one hour to remove the Ag catalyst. Finally, they were rinsed with ultrapure water and dried at room temperature. DOI: 10.1177/0003702819881222 4 The AgNPs were deposited on the SiNWs samples through chemical bath immersion in an aqueous solution containing Ag NO3 and HF. [START_REF] Liu | Improved Photovoltaic Performance of Silicon Nanowire/Organic Hybrid Solar Cells by Incorporating Silver Nanoparticles[END_REF] The SERS substrates were immersed in 10 -6 M of R6G solution. The samples were rinsed after each immersion step using ultrapure water and they were dried at room temperature.

Experimentation Using the Doehlert Design

The response of the SERS substrate based on SiNWs is influenced by several experimental conditions. We present here the optimization of AgNPs-SiNWs elaboration by varying silver nitrate concentration, immersion time in silver nitrate solution and R6G solution using the multivariate method to improve the SERS response. This method is interesting because it enables to simultaneously evaluate the effect of the three chosen factors (silver nitrate concentration, immersion times in silver nitrate and R6G solution) on the SERS response and to examine their interaction as well as to determine the optimal conditions. [START_REF] Lewis | Pharmaceutical Experimental Design[END_REF] The as-adopted method is based on statistical tools to limit the number of experiments in selected experimental field. [START_REF] Pinkowska | Application of the Doehlert Matrix for the Determination of the Optimal Conditions of Hydrolysis of Soybean Protein in Subcritical Water[END_REF] In our case, Doehlert matrix was chosen as the most adequate multivariate systems because the experimental points are uniformly distributed in the studied field. Furthermore, it allows adding new experimental points without changing the matrix and minimizing the number of experimental points. For a number of factors k, the number of experiments is k 2 +k+1+n0, where n0 is the number of experiments repetition of the center point. [START_REF] Doehlert | Uniform Shell Design[END_REF] Using the three factors mentioned above, the number of Doehlert experiment is 13 without the experiment repetition as can be observed in Figure 1. [START_REF] Uchinokura | Raman Scattering by Silicon[END_REF] Silver nitrate concentration has five levels, silver nitrate time has seven levels and R6G time has three levels. Taking into account two additional experiments at the center, the designed experimental methodology is then composed of fifteen experiments presented in Table I with the experimental values of the studied factors U1, U2, U3, which represent respectively the silver concentration in M, the immersion time in silver nitrate in second, the immersion time in RG6 in minute. All experiments were conducted in a random order to avoid any tendency in the measurements. [START_REF] Sena | Application of Doehlert Experimental Design in the Optimization of Experimental Variables for the Pseudozyma sp. (CCMB 306) and Pseudozyma sp. (CCMB 300) Cell Lysis[END_REF] To evaluate the influence of the operating parameters and the interaction between them on the experimental SERS response, coded variables were used.

The transformation of real variable (Ui) into coded variable (Xi) is made using the following equation.

X i =(U i -U i0 )/Δu i (1)
where Xi is the coded variable i, Ui is the real variable i, Ui0 is the value of the real variable i in the center, and Δui is the range of variation of the real variable i.

The investigated experimented domain of the studied factors Ui are summarized in Table II.

For three variables, the experimental response associated by a Doehlert matrix was represented by a quadratic polynominal model Y:

Y=b 0 + b 1 X 1 + b 2 X 2 + b 3 X 3 + b 11 X 1 2 + b 22 X 2 2 + b 33 X 2 2 + b 12 X 1 X 2 + b 13 X 1 X 3 + b 23 X 2 X 3 (2)
with bi as the principal estimation of the effect of factor Ui on the exprimental response, bii as the second estimation of the effect of factor Ui on the exprimental response, bij as the estimation of the effect of the interaction between the factor Ui and the factor Uj.

The calculation of the coefficients b (bi, bii or bij) was obtained using the least squares method:

b=(X T X) -1 X T Y (3)
where Y is the vector of experimental response, X is the model matrix, and X T is the transposed model matrix. however, the intensities of Raman peaks in the spectrum of R6G (10 -6 M)/SiNWs are negligible compared to those of SERS peaks in the R6G (10 -6 M)/AgNPs/SiNW spectrum. These modes are in agreement with the R6G vibration modes defined by Hildbrand in 1984. [START_REF] Hildebrandt | Surface-Enhanced Resonance Raman Spectroscopy of Rhodamine 6G Adsorbed on Colloidal Silver[END_REF][START_REF] Dieringer | Surface-Enhanced Raman Excitation Spectroscopy of a Single Rhodamine 6G Molecule[END_REF] We present here some of the strongest R6G vibration modes: The mode at 610 cm -1 ( noted L1) is attributed to C-C-C ring in-plane bending mode, the modes at 1360 cm -1 ( noted L2), 1509 cm -1 ( noted L3), and 1650 cm -1 ( noted L4) are attributed to aromatic C-C stretching modes. LSi correspond to silicon mode at 520 cm -1 . The appearance of a strong R6G SERS signal after the addition of the silver nanoparticles clearly demonstrates the theoretical approach regarding Raman signals enhancement via the resonant excitation of the surface Plasmon by photon interaction.

Results and Discussion

Morphology of SiNWs and AgNPs/SiNWs

The SERS effect mainly arises from the nanostructure Plasmon resonance, it is typically wavelength dependent. In order to optimize the wavelength excitation, the R6G SERS substrate (R6G (10 -6 M)/AgNPs/SiNWs) was excited by two excitation wavelengths which are 488 nm and 632.8 nm. Figure 4 shows that the majority of R6G vibration modes of the substrate excited at 632.8 nm are more enhanced than those of the substrate excited at 488 nm. For this reason, we have completed the work by the excitation wavelength 632.8 nm.

Optimization of SERS Substrate Elaboration Using the Doehlert Matrix

After the elaboration of the fifteen experiments previously presented in Figure 1, Table I, and

Table II, the samples were characterized using Raman spectroscopy. Figure 5 shows the evolution of the four R6G vibrational modes at 610 cm -1 (L1), 1360 cm -1 (L2), 1509 cm -1 (L3), and 1650 cm -1 (L4) normalized by the silicon mode at 520 cm -1 (LSi) for the fifteen samples.

These normalized modes are respectively noted (L1/LSi), (L2/LSi), (L3/LSi) and (L4/LSi).

According to Figure 5, we can conclude that the four Raman vibrational modes of the fifteen samples have the same evolution. In addition, the most intense R6G peak is located at 1360 cm -1 (L2/LSi). Therefore, L2/LSi was chosen in the Doehlert model in order to follow the evolution of the SERS response by varying the three parameters. In addition, Figure 5 clearly shows that the Raman enhancement is optimal for both samples N°2 and N°15 (rand) of the experimental design methodology (Table I).

The effects of silver nitrate concentration as well as immersion times in silver nitrate and R6G solution on SERS response were investigated using Doehlert matrix.

Taking into consideration the experimental results previously presented, the rate L2/LSi were introduced for all the fifteen samples in Nemrod software, as well as the the data defined in Table I. Then, the coefficients b (bi, bii, bij) of the polynominal model defined in Eq. 2 were determined by the treatement of all the data. They are presented in the following equation Y=0.067 + 0.024 X 1 -0.021 X 2 + 0.16 X 3 + 0.048 X 1 2 + 0.092 X 2 2 + 0.16 X 3 2 +0.017

X 1 X 2 -0.12 X 1 X 2 -0.048 X 1 X 3 (4)
The obtained coefficients are also presented in graphical form in Figure 6 in order to visualize easily the estimation of bi and bii coefficients, respectively corresponding to the principal and the second effect of the factor Ui on the experimental response, as well as the estimation of bij coefficients corresponding to the interaction effect between the factor Ui and the factor Uj.

The coefficient b1 provides information on the effect of silver nitrate concentration, whereas b2 and b3 provide information on immersion time respectively in silver nitrate and R6G solution.

The coefficients presented in the Figure 6 well show that the main influencing parameter on the studied response is the R6G immersion time (b3 is greater than b1 and b2). Silver nitrate concentration and the R6G immersion time have positive estimation effects on the experimental response (since b1 and b3 are positive). Thus, an increase on the both parameters leads to the augmentation of the observed response. However, an increase in silver nitrate incubation time decreases the response (since b2 is negative). The b12 coefficient, which corresponds to the interaction between the silver nitrate concentration and the silver nitrate immersion time, is positive as noted in Figure 6. Thus, there is an interaction between both parameters. [START_REF] Hannachi | Preparation of Manganese Sulfide (MnS) Thin Films by Chemical Bath Deposition: Application of the Experimental Design Methodology[END_REF][START_REF] Mathieu | New Efficient Methodology for Research Using Optimal Design (NEMROD) Software[END_REF] To properly interpret the results of the parameter study, the predicted contour plots and their corresponding three-dimensional representations are given in the Figures 6e to 7a Accordingly, we can conclude that the maximum response is obtained by increasing R6G immersion time.

Figures 7c andd also indicate the SERS response is enhanced by increasing silver nitrate concentration from above 5-7 mM. Moreover, the Raman response seems to be more sensitive to the variation of R6G immersion time than the variation of the silver nitrate concentration. This could be explained by a longer immersion time which would enable the adsorption of R6G molecules in organized manner on the silver nanoparticle surface.

For an R6G immersion time of 38 minutes, the contour plots and their corresponding three-dimensional are displayed in Figures 7e andf. An analysis of these figures shows that the response is maximal for a broad area of silver nitrate concentration and silver nitrate immersion time. The contour plots confirm also the interaction between both parameters.

Consequently, to obtain a maximized SERS response, the SiNWs must be immersed in a silver nitrate concentration ranging from 5 to 7 mM during 67 to 68 seconds, followed by R6G solution immersion for more than 38 minutes.

Limit of Detection and Enhancement Factor

The limit of detection (LOD) is recognized as one of the main parameters to characterize a given sensing platform. It informs about the lowest target concentration to be detected by the platform.

Based on the optimal conditions obtained by the previously work, Figure 8 presents the SERS spectra of the optimized AgNPs/SiNWs substrates (Sample 2 of the experimental design) soaked in R6G solutions of different concentrations ranging from (b) 10 -11 M to (e) 10 -6 M. As noted, the most R6G vibration modes are visible at the R6G SERS spectra of the different concentrations greater than 10 -11 M. The most prominent modes located at 610 cm -1 , 779 cm -1 , 1509 cm -1 , 1570 cm -1 , and 1650 cm -1 have undergone a decrease in their scattered intensity when the R6G concentration varies from 10 -6 to 10 -11 M. We can conclude that the limit of detection is about 10 -11 M.

Enhancement factor (EF) is also a very important parameter to demonstrate the effectiveness of SERS. In this framework, we estimate the EF by the following equation, which was used by Ngan et al.: [START_REF] Ngan | Improvement of Raman Enhancement Factor Due to the Use of Silver Nanoparticles Coated Obliquely Aligned Silicon Nanowire Arrays in SERS Measurements[END_REF] EF = I SERS /N Raman I Raman /N SERS (5) where ISERS and IRaman are the Raman intensity of a diagnostic band of the probe molecules adsorbed on the SERS and non-SERS substrate, respectively. NSERS is the number of molecules on the SERS substrate and NRaman is the number of molecules on the non-SERS substrate, which caused respectively ISERS and IRaman. To estimate the EF of the SERS substrate, the SERS spectrum of R6G at a concentration of 10 -10 M (R6G (10 -10 )/AgNPs/SiNWs) and a Raman spectrum of R6G at a concentration of 10 -6 M (R6G (10 -6 )/SiNWs) were used. The Raman band at 1360 cm -1 (L2) has been chosen to calculate the EF.

With the assumption that R6G solution is evenly distributed on the surface of the sample, the ratio NSERS/NRaman is considered as equal to the ratio of the R6G concentrations in which the non-SERS and the SERS substrates were incubated, respectively. The SERS intensity at 1360 nm is equal to 1600 and the Raman intensity is equal to 29. The EF estimated to 0.6 10 6 .

Conclusion

In this study, an enlarged surface area of vertically aligned silicon nanowires decorated with silver nanoparticles was used in order to build SERS substrates capable of detecting a remarkably lower R6G concentration under optimal conditions. The Doehlert design was adopted in order to estimate the effects of the studied parameters on the R6G SERS response and to obtain the optimal conditions for SERS substrates development. The multivariate method showed that the R6G immersion time was the main influential factor to the response. In addition, an interaction between the silver nitrate concentration and silver nitrate immersion time was

found. An optimal SERS response can be achieved by immersion the SiNWs in 5 to 7 mM DOI: 10.1177/0003702819881222 10 nitrate silver concentrations during 67-68 seconds followed by R6G solution immersion for more than 38 minutes. The sensing platform characterized by a limit of detection about 10 -11 M and an enhancement factor estimated to 0.6 10 6 .

Tables

Figure 6 SurfaceFigure 3

 63 Figure2ashows SEM images of SiNWs prepared using metal-assisted chemical etching (MACE). The SiNWs cross sectional image presented in Figure2aindicates an identical large area of nanowires having a one-dimensional structure. They are distinguishable and the majority of them are vertical with regard to the sample surface. The layer thickness of the nanostructure is about 20 µm. The synthesized AgNPs have a spherical shape as it can be seen in Figure2b, which present a scanning electron microscopy (SEM) image of the surface AgNPs/SiNWs. The nanoparticles are well dispersed on the SiNWs surface with a non-uniform particle size in the range around of 40-120 nm.

  . The zone characterized by a maximal SERS response overall the three-dimensional representations is in red as it can be seen in the various representations. A AgNO3 concentration equals 5 mM, Figures 6b and 7a depict the contour plots and their corresponding three-dimensional representations showing the SERS response with regard to the silver nitrate and R6G immersion time. Both indicate that the maximal SERS response is obtained by increasing the R6G immersion time to above 38-40 minutes. The contours are exhibited by circles centered around the middle of the studied field of silver nitrate immersion time denoting a maximal SERS response located in the middle of this field. For a silver nitrate immersion time of 67-68 seconds, Figures 6d and 7c depict the contour plots and their corresponding three-dimensional representations showing the SERS response with respect to silver nitrate concentration and immersion time in R6G solutions.

Figure 2 .

 2 Figure 2. Results obtained using SEM: (a) Cross-sectional SEM image of SiNWs; (b) surface SEM image of AgNPs/SiNWs.

Figure 3 .

 3 Figure 3. Raman spectra of SiNWs (a) ; R6G (10 -6 M)/ SiNWs (b) and R6G (10 -6 M)/AgNPs/SiNWs (c) for an excitation at 632.8 nm; L1, L2, L3 and L4 peaks correspond to the peaks of the main R6G vibration modes, LSi corresponds to a Si vibrational mode. The insert represents a zoom of the spectra a and b in the range 1200-1800 cm -1 .

Figure 4 .

 4 Figure 4. SERS spectra of R6G (10 -6 M) /AgNPs/SiNWs with two excitation wavelengths 488 nm and 632.8 nm.

Figure 5 .

 5 Figure 5. Evolution of normalized intensity of selected modes from R6G of the fifteen experiments prepared according the experimental design.

Figure 6 .

 6 Figure 6. Graphical analysis of effects on the experimental response. bi coefficients represent the estimation of the main effects of Ui ; bii coefficients represent the estimation of the second order effects of Ui ; bij coefficients represent the interactions between Ui and Uj.

Figure 7 . 18 Figure 8 .

 7188 Figure 7. Results obtained using the Doehlert matrix: (a) Contour plots of SERS responses with respect to silver nitrate time (seconds) and R6G time (minutes) at a fixed silver nitrate concentration 5 mM; (b) three-dimensional corresponding plot; (c) contour plots of SERS responses versus silver nitrate concentration (M) and R6G time (minutes) at a fixed silver nitrate time 67.5 seconds; (d) three-dimensional corresponding plot; (e) contour plots of SERS responses with respect to silver nitrate time (seconds) and silver nitrate concentration (M) at a fixed R6G time 38 minutes and F) three-dimensional corresponding plot.

Table I .

 I Doehlert experimental matrix.

	DOI: 10.1177/0003702819881222
	11	2	-0.5000 0.2887	0.8165	0.00258	82.73	36.49
	12	15	0.0000	-0.5778 0.8165	0.00505	37.04	36.49
	13	10	0.0000	0.0000	0.0000	0.00505	67.50	20.50
	14	6	0.0000	0.0000	0.0000	0.00505	67.50	20.50
	15	14	0.0000	0.0000	0.0000	0.00505	67.50	20.50
	Experiment	Coded variables		Real variables
	number
	Before	After	X1	X2	X3	U1	U2	U3
	rand	rand
	1	4	1.0000	0.0000	0.0000	0.01000	67.50	20.50
	2	7	-1.0000 0.0000	0.0000	0.00010	67.50	20.50
	3	8	0.5000	0.8660	0.0000	0.00753	113.18	20.50
	4	12	-0.5000 -0.8660 0.0000	0.00258	21.82	20.50
	5	3	0.5000	-0.8660 0.0000	0.00753	21.82	20.50
	6	11	-0.5000 0.8660	0.0000	0.00258	113.18	20.5
	7	13	0.5000	0.2887	0.8165	0.00753	82.73	36.49
	8	5	-0.5000 -0.2887 -0.8165 0.00258	52.27	4.51
	9	1	0.5000	-0.2887 -0.8165 0.00753	52.27	4.51
	10	9	0.0000	0.5778	-0.8165 0.00505	97.96	4.51
					19

Table II .

 II The experimental field.

		Factor Ui	Unit	Center (Ui0)	the range of variation (∆ui)
	U1	Silver nitrate concentration	M	0.00505	0.00495
	U2	Silver nitrate time	Second	67.50	52.72
	U3	R6G time	Minute	20.50	19.58

Acknowledgments

The authors acknowledge the financial support from the University of Tunis El Manar for the mobility grant (Bourse d'alternance) awarded to Awatef Ouhibi.


Xavier Bonnetain
email: xavier.bonnetain@inria.fr

Collisions on Feistel-MiMC and univariate GMiMC

Keywords: MiMC, MPC, symmetric cryptanalysis

published or not. The documents may come

Description of the ciphers 1.MiMC-2n/n

MiMC-2n/n [START_REF] Albrecht | MiMC: Efficient Encryption and Cryptographic Hashing with Minimal Multiplicative Complexity[END_REF] is a 2n-bit block size, n-bit key block cipher. It claimed n bits of security. Its round function is described in Figure 1, and can be written as

R i k (x L , x R ) = x R ⊕ (x L ⊕ k ⊕ c i ) 3 , x L . x L k ⊕ c i x 3 x R

GMiMC

GMiMC [START_REF] Albrecht | Feistel Structures for MPC, and More[END_REF] generalizes the MiMC-2n/n construction to generalized Feistels. Two key schedules are proposed. The univariate key schedule uses a fixed key for each round, while the multivariate key schedule uses t initial keys and updates the round keys. Their claimed security corresponds to the number of bits of the key. Four generalized feistel constructions are proposed:

GMiMC-crf. GMiMC-crf has t branches and adds a function of t -1 branches on one branch. The round function is

R i k (x 1 , . . . , x t ) = x 2 , . . . , x t , x 1 ⊕   t ⊕ j=2 x j ⊕ k ⊕ c i   3 .
GMiMC-erf. GMiMC-erf has t branches, and adds a function of one branch on all the other. The round function is

R i k (x 1 , . . . , x t ) = x 2 ⊕ (x 1 ⊕ k ⊕ c i ) 3 , . . . , x t ⊕ (x 1 ⊕ k ⊕ c i ) 3 , x 1 .
GMiMC-Nyb. GMiMC-Nyb has 2t branches, and adds a function of each odd branch to the next branch. The round function is

R i k (x 1 , . . . ,x t ) = x 2 ⊕ (x 1 ⊕ k ⊕ c ti ) 3 , x 3 , x 4 ⊕ (x 3 ⊕ k ⊕ c ti+1 ) 3 , . . . , x 2t ⊕ (x 2t-1 ⊕ k ⊕ c ti+t-1 ) 3 , x 1 .
GMiMC-mrf. GMiMC-mrf is a generalization of the previous construction with a permutation of the branches that change for each round.

Attacks

Attack on MiMC-2n/n

The attack relies on an invariant property of the round function, and can be seen as a slight generalization of a slide attack presented in [START_REF] Bonnetain | On Quantum Slide Attacks[END_REF]. The invariant property is described in Figure 2.

x k ⊕ c i x 3 y z t c i x 3 y ⊕ k x ⊕ k z ⊕ k t ⊕ k Figure 2: Illustration of Lemma 1 Lemma 1. Let R i k be the round function of MiMC-2n/n with the key k for round i. Then for all x, y, k, i, R i k (x, y) ⊕ (k, k) = R i 0 (x ⊕ k, y ⊕ k) Proof. R i 0 (x ⊕ k, y ⊕ k) = ( y ⊕ k ⊕ (x ⊕ k ⊕ c i ) 3 , x ⊕ k ) = ( y ⊕ (x ⊕ k ⊕ c i ) 3 , x ) ⊕ (k, k) = R i k (x, y) ⊕ (k, k) Theorem 1. Let E k be MiMC-2n/n with the key k. Then, for all x, y, k, E k (x, y) ⊕ (k, k) = E 0 (x ⊕ k, y ⊕ k).
Proof. By induction over the number of rounds. The base case is Lemma 1. If the property holds after i -1 rounds, then

(R i-1 k • R i-2 k • • • • R 1 k )(x, y) ⊕ (k, k) = (R i-1 0 • R i-2 0 • • • • R 1 0 )(x ⊕ k, y ⊕ k).
By Lemma 1,

(R i 0 • R i-1 0 • • • • R 1 0 )(x ⊕ k, y ⊕ k) = R i 0 ((R i-1 0 • R i-2 0 • • • • R 1 0 )(x ⊕ k, y ⊕ k)) =R i 0 ((R i-1 k • R i-2 k • • • • R 1 k )(x, y) ⊕ (k, k)) = (R i k • R i-1 k • • • • R 1 k )(x, y) ⊕ (k, k) Corollary 1. Let E k be MiMC-2n/n with the key k. Let f (x) = E k (x, x) ⊕ (x, x) and g(x) = E 0 (x, x) ⊕ (x, x). Then f (x) = g(x ⊕ k).
Proof.

g(x ⊕ k) = E 0 (x ⊕ k, x ⊕ k) ⊕ (x ⊕ k, x ⊕ k) = E k (x, x) ⊕ (k, k) ⊕ (x ⊕ k, x ⊕ k) = E k (x, x) ⊕ (x, x) = f (x)
Key recovery. The key recovery simply consists in looking for a collision between f and g from Corollary 1, which can be done in time 2 n/2 as the two functions have an n-bit input. This contradicts the claim of n bits of security of MiMC-2n/n.

Hash function.

MiMC can be used keyless as a permutation for a sponge-based hash function. As there is no key in this construction, it is unclear how Theorem 1 could be used to attack the hash function.

Attacks on GMiMC

In most cases ,the same property can be found in univariate GMiMC, that is,

E k (x 1 , . . . , x t )⊕ (k, . . . , k) = E 0 (x 1 ⊕ k, . . . , x t ⊕ k)
, which allows to apply the same attack as in the MiMC-2n/n case.

GMiMC-Nyb and GMiMC-mrf. One round of GMiMC-Nyb and GMiMC-mrf can be seen, up to a permutation of the branches, as t Feistel in parallel. Hence, the property holds.

GMiMC-erf.

The added function only depends on one input branch, hence the property also holds.

GMiMC-crf.

The function is slightly different in that case, as it depends on more than one branch. For the property to hold, we must have that ((

⊕ t j=2 x j ) ⊕ k ⊕ c i ) 3 = ( ⊕ t j=2 (x j ⊕ k) ⊕ c i ) 3 .
Hence, the property holds only if t is even.

Variants in large characteristics

MiMC and GMiMC can also be defined over a finite field of large characteristic. In that case, the property we have is E k (x 1 , . . . , x t ) + (k, . . . , k) = E 0 (x 1 + k, . . . , x t + k), and the same attack can be applied. The only exception is GMiMC-crf, where we need to have k + k = 0 for the property to hold.

Quantum attacks

The collision property corresponds to a hidden period, and as such, permits a key recovery in O(n) quantum queries. With a restriction to classical queries, these attacks happens to be in a form suitable for the offline Simon's algorithm [START_REF] Bonnetain | Quantum Attacks without Superposition Queries: the Offline Simon Algorithm[END_REF], which allows to make a key recovery in O(2 n/3 ) classical queries and quantum time.

Conclusion

We have shown that MiMC-2n/n and all the versions of univariate GMiMC except some instances of GMiMC-crf are vulnerable to a collision attack. More generally, this demonstrates that using round constants is not enough for a key schedule to secure a Feistel or generalized Feistel construction. This attack does not appear to be applicable to the other MiMC construction, MiMCn/n, nor to the hash functions based on any version of MiMC or GMiMC.

Figure 1 :

 1 Figure 1: MiMC-2n/n round function

Acknowledgements. The author would like to thank Chaoyun Li for his presentation of MiMC and for suggesting to apply the attack to GMiMC. The author would also like to thank the authors of MiMC and GMiMC for their comments and discussions on this result. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement no. 714294 -acronym QUASYModo).