David Lubicz

Damien Robert

A generalisation of Miller's algorithm and applications to pairing computations on abelian varieties

In this paper, we use the theory of theta functions to generalize to all abelian varieties the usual Miller's algorithm to compute a function associated to a principal divisor. We also explain how to use the Frobenius morphism on abelian varieties defined over a finite field in order to shorten the loop of the Weil and Tate pairings algorithms. This extend preceding results about ate and twisted ate pairings to all abelian varieties. Then building upon the two preceding ingredients, we obtain a variant of optimal pairings on abelian varieties. Finally, by introducing new addition formulas, we explain how to compute optimal pairings on Kummer varieties. We compare in term of performance the resulting algorithms to the algorithms already known in the genus one and two case.

Introduction

The computation of Weil and Tate pairings has important applications in arithmetic and cryptography. Almost all the known algorithms to compute these pairings on elliptic curves rely on Miller's algorithm [START_REF] Victor | The Weil pairing, and its efficient calculation[END_REF] to evaluate at a certain point a function associated to a principal divisor.

The improvements over the initial version of Miller's algorithm have followed two main approaches:

-making the basic loop of Miller's algorithm quicker with efficient arithmetic; -reducing the number of loops of the algorithm by using endomorphisms of the curve.

For a curve defined over a finite field, using the absolute Frobenius endomorphism, this last line of ideas has led to the definition of eta-pairings [START_REF] Paulo | Efficient pairing computation on supersingular abelian varieties[END_REF], ate-pairings [START_REF] Hess | The eta pairing revisited[END_REF][START_REF] Granger | Ate pairing on hyperelliptic curves[END_REF] and optimal pairings [START_REF] Vercauteren | Optimal pairings[END_REF].

The paper [START_REF] Lubicz | Efficient pairing computation with theta functions[END_REF] describes a new algorithm based on the theory of theta functions to compute Weil and Tate pairings which apply to all abelian varieties.

It is a natural question to ask whether the known optimizations of the classical Miller's algorithm on elliptic curves can be adapted to the algorithms presented in [START_REF] Lubicz | Efficient pairing computation with theta functions[END_REF]. In this paper, we focus on the optimisations which consist in reducing the number of loops in pairing computation algorithms by using non trivial endomorphisms of the abelian varieties. We won't deal here with the conversion between Weierstrass coordinates of an elliptic curve or Mumford coordinates of a Jacobian of an hyperelliptic curve and theta functions. These conversion formulas come from the well known Thomae's formula, and are described in more details in [START_REF] Gaudry | The arithmetic of characteristic 2 Kummer surfaces and of elliptic Kummer lines[END_REF][START_REF] Gaudry | Fast genus 2 arithmetic based on Theta functions[END_REF][START_REF] Wamelen | Equations for the Jacobian of a hyperelliptic curve[END_REF][START_REF] Cosset | An algorithm for computing (ℓ, ℓ)-isogenies in polynomial time on jacobians of hyperelliptic curves of genus 2[END_REF]. We note that such a formula is specific to a model of a curve. Although in the case of the Weierstrass model of an elliptic curve it is immediately provided by the definition of the group law, for other models it may require some computations to obtain it (see for instance [START_REF] Arène | Faster computation of the Tate pairing[END_REF][START_REF] Galbraith | Hyperelliptic pairings[END_REF]). Using the classical Riemann relations for theta functions, we generalize Miller's algorithm to all abelian varieties. In our context, we compute the function defined up to a constant factor by a divisor D on an abelian variety with Chern class 0 linearly equivalent to 0. Our method rely on the projective embedding of a principally polarized abelian variety provided by a power of its theta divisor and therefore is not restricted to Jacobians of curves.

Classical pairing computation algorithms rely on

Once we know how to compute the Miller's functions, it is easy to apply the optimisations of the ate and optimal ate pairings to our setting. However, the formula that we obtain to compute the Weil and Tate pairings using the classical Miller's algorithm with theta functions is slower than the algorithms presented in [START_REF] Lubicz | Efficient pairing computation with theta functions[END_REF], which use a different (and faster) subset of the Riemann relations that we call differential additions. Our second contribution is to extend the improvements of the ate and optimal ate pairings to this modified Miller's loop based on differential additions by studying the action of the absolute Frobenius endomorphism on points given by the theta coordinates.

To get even faster formulas, one can compute with a Kummer variety, which is the quotient of an abelian variety by the (-1) automorphism acting on it, by using theta functions of level 2 which are invariants by the action of (-1). There is no addition law on a Kummer variety, since we can not distinguish a point from its opposite. Still, there is an action of N on the points of a Kummer variety that one can compute with the differential additions. By carrying the pairing to this structure, we were led to introduce symmetric pairings in [START_REF] Lubicz | Efficient pairing computation with theta functions[END_REF] by using differential additions on level 2 theta functions. Our third contribution is then to extend all the results previously discussed to this setting; this can be seen as a generalization of [START_REF] Steven | Computing pairings using x-coordinates only[END_REF] to higher dimension. For this, we introduce new addition laws deduced from Riemann relations that can be used in the case of level 2 theta functions. Finally, we give an overview of the theoretical complexity of our algorithms.

In order to avoid too much formalism, we have chosen to present all the formulas of this paper using the classical analytic theory of theta functions. Nonetheless, it should be understood that most of our algorithms apply in general to abelian varieties defined over any field of characteristic not equal to 2. To see this one can invoke the Lefschetz's principle or use Mumford's theory of algebraic theta functions. There are part of the paper which are more specific to abelian varieties over finite field because we use the Frobenius endormophism in order to speed up the computations. In this case, we clearly state this hypothesis at the beginning of the section.

The paper is organised as follows: in Section 2, we recall some notations and well known results about theta functions. In Section 3, we describe different operations which can be computed on any abelian variety using Riemann relations. Section 4 is devoted to a generalisation of Miller's algorithm. In Section 5, we recall the standard definition of the Weil and Tate pairings on abelian varieties, and explain how to compute them using the results from the previous section. We give two ways of computing these pairings: the usual Miller's algorithm applied to theta functions and the utilisation of differential additions as in [START_REF] Lubicz | Efficient pairing computation with theta functions[END_REF]. In Section 6, we explain how to extend this to compute the ate pairing on an abelian variety defined over a finite field, while in Section 7, the case of the optimal ate pairing is treated. Finally in Section 8, we give an overview of the complexity.

Some notations and basic facts

In this section, in order to fix the notations, we recall some well known facts about analytic theta functions (see for instance [START_REF] Mumford | Tata lectures on theta I[END_REF][START_REF] Ichi Igusa | Theta functions[END_REF]). Let H g be the g dimensional Siegel upper-half space which is the set of g × g symmetric matrices Ω whose imaginary part is positive definite. For Ω ∈ H g , we denote by Λ Ω = ΩZ g + Z g the lattice of C g defined by Ω. Any abelian variety A of dimension g over C with a principal polarisation is analytically isomorphic to C g /Λ Ω for a certain Ω ∈ H g . In the rest of this paper, we denote by π : C g → C g /Λ Ω = A the canonical projection. The classical theory of theta functions gives a lot of functions on C g that are pseudo-periodic with respect to Λ Ω and can be used as a projective coordinate system for A . More precisely, for a, b ∈ Q g , the theta function with rational characteristics (a, b) is an analytic function on C g × H g given by:

θ [ a b ] (z, Ω) = n∈Z g exp πi t (n + a).Ω.(n + a) + 2πi t (n + a).(z + b) . (1)
For all m, n ∈ Z g , we have:

θ [ a b ] (z+Ωm+n) = exp(2πi( t a.n-t b.m)-πi t m.Ω.m-2πi t m.z)θ [ a b ] (z, Ω). (2)
We say that a function f on

C g is Λ Ω -quasi-periodic of level n ∈ N if for all z ∈ C g and m ∈ Z g , we have:f (z + m) = f (z), f (z + Ω.m) = exp(-πin t m.Ω.m - 2πin t z.m)f (z).
For any n ∈ N * , the set H Ω,n of Λ Ω -quasi-periodic functions of level n is a finite dimensional C-vector space whose basis can be given by the theta functions with characteristics:

(θ 0 b/n (z, n -1 .Ω)) b∈[0,...,n-1] g . If n = k 2 , then an alternative basis of H Ω,n is (θ a/k b/k (kz, Ω)) a,b∈[0,...,k-1] g . A theorem of
Lefschetz tells that if n ≥ 3, the functions in H Ω,n give a projective embedding of A in P n g -1 , the projective space over C of dimension n g -1. For n = 2, the functions in H Ω,2 do not give a projective embedding of A . Indeed, it is easy to check that for all f ∈ H Ω,2 , we have f (-z) = f (z). Under some well known general conditions [START_REF] Koizumi | Theta relations and projective normality of Abelian varieties[END_REF]Corollary 4.5.2], the image of the embedding defined by H Ω,2 in P n 2 -1 is the Kummer variety associated to A , which is the quotient of A by the automorphism -1.

Once we have chosen a level n ∈ N and a matrix period Ω, for the rest of this paper, we adopt the following conventions: we let Z(n) = (Z/nZ) g and for a point z P ∈ C g and i ∈ Z(n) we put θ i (z P ) = θ 0 i/n (z P , Ω/n). For ℓ, n ∈ N, such that ℓ divides n we will implicitly consider Z(ℓ) as a subgroup of Z(n) via the morphism x → (n/ℓ).x.

We denote by Θn the theta divisor of level n on A which is the divisor of zero of θ [ 0 0 ] (z, n -1 .Ω). There is an isogeny

ϕ n : A → Â = Pic 0 A , defined by x → τ *
x Θn -Θn where τ x is the morphism of translation by x on A . The kernel of ϕ n is A [n]. For n = 1 we let Θ1 = Θ. If L n is the line bundle corresponding to the polarization ϕ n (or the divisor Θn), we have L n = L n 1 . We denote by K(A ) the function field of A and if f ∈ K(A ), we denote by (f ) the divisor of the function f . Let Z 0 (A ) be the group of 0-cycles of A that is the free commutative group over the set of closed points of A . If D = i∈I n i (P i ) is an element of Z 0 (A ), we let Supp(D) be the reduced zero dimensional variety ∪ i∈I P i . If f ∈ K(A ) has no poles nor zeroes on Supp(D), we put f (D) = i∈I f (P i ) ni .

We recall the following theorem from [START_REF] Lubicz | Efficient pairing computation with theta functions[END_REF]Theorem 1] which is a version of the usual Riemann addition formula for theta functions: Theorem 1. Let i, j, k, l ∈ Z(2n). We suppose that i + j, i + k and i + l ∈ Z(n). Let Ẑ(2) be the dual group of Z(2). For all χ ∈ Ẑ(2) and z 1 , z 2 , z 3 , z 4 ∈ C g , let

L 1 = η∈Z(2) χ(η)θ i+j+η (z 1 + z 2 )θ i-j+η (z 1 -z 2 ), L 2 = η∈Z(2) χ(η)θ k+l+η (z 3 + z 4 )θ k-l+η (z 3 -z 4 ), L 3 = η∈Z(2) χ(η)θ i+k+η (z 1 + z 4 )θ i-k+η (z 1 -z 4 ), L 4 = η∈Z(2) χ(η)θ j+l+η (z 3 + z 2 )θ j-l+η (z 3 -z 2 ), then we have L 1 .L 2 = L 3 .L 4 . (3)
In this section, we explain how to compute certain operations on the set of geometric points of A using Riemann equations. These operations will be used in the algorithms of the next sections, but they may be interesting for other applications, for instance for the purpose of computing on Kummer varieties.

Normal additions

Let A be an abelian variety over C with period matrix Ω. For n a positive integer, we represent A as a closed subvariety of P Z(n) by the way of level n theta functions and we suppose that 2 divides n. Such an embedding is uniquely determined once we have chosen a numbering of a basis of H Ω,n : by default we take {θ i , i ∈ Z(n)}. With this convention, by a theorem of Mumford, if 4|n, the resulting embedding of A in P n g -1 is defined by the equations of Theorem 1 by taking z 2 = z 3 = z 4 = 0.

Let K be a number field. In the rest of this paper, we suppose that this embedding is defined over K or, say in another way, that the projective point 0 = (θ i (0)) i∈Z(n) corresponding to the neutral element of A is defined over K. Now if 4|n, from the knowledge of P = (P i ) i∈Z(n) and Q = (Q i ) i∈Z(n) , one can compute the (projective) point P + Q = ((P + Q) i ) i∈Z(n) using Theorem 1 with z 3 = z 4 = 0. We just have to check that the L 2 factor of equation ( 1) does not vanish too often which is a consequence of [START_REF] Lubicz | Efficient pairing computation with theta functions[END_REF]Proposition 3]. We write

P + Q = NormalAdd(P, Q).
We illustrate this with n = 4 and g = 1 in Algorithm 1. Let E be an elliptic curve defined by Ω ∈ H 1 ; a point z ∈ E will be represented by the projective coordinates (θ i (z)) i∈Z(4) = (ϑ 0 i/4 (z, Ω/4)) i∈Z(4) . Let (θ i (0)) i∈Z(4) = (a, b, c, d). (In all the examples that we give, we assume that we are in a generic setting so that the formulas are well defined. Otherwise we can always choose a different subset of Riemann relations, at least in level 4).

Differential additions

Denote by A the pullback of A via the natural projection κ : A n g → P n g -1 . In the following, we adopt the following convention: if P = (P i ) i∈Z(n) is a point of A , we denote by P = ( P i ) i∈Z(n) an affine lift of P that is a point P of A n g such that κ( P ) = P . We introduce the tilde notation in the formulas where we compute with affine points: it means that we want to distinguish two points lying on the same line cutting the origin of A n g . We also choose once and for all an affine lift of the theta null point 0 = (θ i (0)) i∈Z(n) defined over K.

Next, for all n such that 2|n, from the knowledge of P = (

P i ) i∈Z(n) , Q = ( Q i ) i∈Z(n) and ( P -Q) = (( P -Q) i ) i∈Z(n) , the formula of Theorem 1, defines a unique P + Q = (( P + Q) i ) i∈Z(n)
which is an affine lift of P + Q. Following [START_REF] Lubicz | Efficient pairing computation with theta functions[END_REF],

we write P + Q = DiffAdd( P , Q, P -Q).

Algorithm 1: Normal addition algorithm in genus 1 level 4 input : The points x = (x0, x1, x2, x3) and y = (y0, y1, y2, y3) on E. output : The projective coordinates (θi(x + y)) of the point t = x + y.

1 return t0 = (x 2 0 + x 2 2 )(y 2 0 + y 2 2 ) (a 2 + c 2 ) + (x 2 0 -x 2 2 )(y 2 0 -y 2 2 ) (a 2 -c 2 ) t1 = (x0x1 + x2x3)(y0y1 + y2y3) (ab + cd) + (x0x1 -x2x3)(y0y1 -y2y3) (ab -cd) t2 = (x 2 1 + x 2 3 )(y 2 1 + y 2 3 ) (a 2 + c 2 ) + (x 2 1 -x 2 3 )(y 2 1 -y 2 3 ) (a 2 -c 2 ) t3 = (x0x3 + x2x1)(y0y3 + y2y1) (ac + bc) + (x0x3 -x2x1)(y0y3 -y2y1) (ac -bc)
Chaining the algorithm DiffAdd in a classical Montgomery ladder [4, Algorithm 9.5 p. 148] yields an algorithm that takes as inputs

Q = ( Q i ) i∈Z(n) , P + Q = (( P + Q) i ) i∈Z(n) , P = ( P i ) i∈Z(n) , 0 = ( 0 i ) i∈Z(n)
and an integer ℓ and outputs Q + ℓP . We write Q + ℓP = ScalarMult(ℓ, P + Q, P , Q, 0). It is proved in [START_REF] Lubicz | Efficient pairing computation with theta functions[END_REF] that the output of ScalarMult does not depend of the addition chain used to compute it.

There is a natural action by multiplication of K * on the coordinate of a point in A n g that we denote by λ * P for λ ∈ K * and P ∈ A n g (K).

Lemma 1. Let P, Q ∈ A (K) and let P , Q, P + Q be affine lifts of P , Q and

P + Q. Let R = ScalarMult(ℓ, P + Q, P , Q, 0). Let α, β, γ, δ ∈ K * , we have [20, Remark 3] ScalarMult(ℓ, α * P + Q, β * P , γ * Q, δ * 0) = (α ℓ β ℓ(ℓ-1) /γ ℓ-1 δ ℓ(ℓ-1) ) * R, (4)
ScalarMult(ℓ, α * P , α * P , δ * 0, δ * 0) = α ℓ 2 δ ℓ 2 -1 * ScalarMult(ℓ, P , P , 0, 0). ( 5)

We illustrate differential additions with our previous g = 1, n = 4 example in Algorithm 2.

Normal additions in level 2

We consider the case n = 2 and we make the hypothesis that level 2 theta functions give a projectively normal embedding of the Kummer variety K A into P Z (2) . This condition is equivalent following [16, output : The affines coordinates P + Q = (t0, t1, t2, t3).

1 return t0 = (x 2 0 + x 2 2 )(y 2 0 + y 2 2 ) z0(a 2 + c 2 ) + (x 2 0 -x 2 2 )(y 2 0 -y 2 2 ) z0(a 2 -c 2 ) t1 = (x 2 1 + x 2 3 )(y 2 0 + y 2 2 ) z1(a 2 + c 2 ) + (x 2 1 -x 2 3 )(y 2 0 -y 2 2 ) z1(a 2 -c 2 ) t2 = (x 2 0 + x 2 2 )(y 2 0 + y 2 2 ) z2(a 2 + c 2 ) - (x 2 0 -x 2 2 )(y 2 0 -y 2 2 ) z2(a 2 -c 2 ) t3 = (x 2 1 + x 2 3 )(y 2 0 + y 2 2 ) z3(a 2 + c 2 ) - (x 2 1 -x 2 3 )(y 2 0 -y 2 2 ) z3(a 2 -c 2 )
to the fact that for all k, l ∈ Z(2) such that t k.l = 0, θ k,l (0) = 0. As for n = 2, we can not distinguish a point Q from its opposite -Q, we can not expect to have a NormalAdd(P, Q) algorithm which returns a point since the result P ± Q is not determined from the input data. Nonetheless, there exists a NormalAdd 2 (P, Q) algorithm whose output uniquely determines the set of points {P + Q, P -Q}.

From the knowledge of P = (θ i (z P )) i∈Z(2) and Q = (θ i (z Q )) i∈Z(2) , it is explained in [20, Section 5.2] that using Riemann's equations (3) with z 3 = z 4 = 0, we can recover for all i, j ∈ Z(2)

κ ij = θ i (z P + z Q )θ j (z P -z Q ) + θ j (z P + z Q )θ i (z P -z Q ), (6)
with the formula

κ ij = χ∈ Ẑ(2)   η∈Z(2) χ(η)θ i+j+η (0)θ i-j+η (0) -1 η∈Z(2) χ(η)θ i+k+η (z P )θ i-k+η (z P ) η∈Z(2) χ(η)θ j+l+η (z Q )θ j-l+η (z Q )   . ( 7
)
We can suppose that θ 0 (z P + z Q )θ 0 (z P -z Q ) = 0, if necessary by replacing the index 0 by another one. For i ∈ Z(2), if we let P i (X) = X 2 -2 κi0 κ00 X + κii κ00 , then the roots of P i (X) are

θi(z P +z Q ) θ0(z P +z Q ) , θi(z P -z Q ) θ0(z P -z Q ) ∈ K. If P, Q ∈ K A (K) are 2-torsion points, P + Q = P -Q ∈ K A (K) so each P i (X) has a double root.
Otherwise, we can suppose that there exists α ∈ Z(2), α = 0 such that the

matrix M = θ 0 (z P + z Q ) θ 0 (z P -z Q ) θ α (z P + z Q ) θ α (z P -z Q ) is invertible. Consider the algebra A = K[X]/(P α (X)) = K[x]
and where x is the image of X is A via the canonical projection. Then we can represent the set {θ α (z P + z Q ), θ α (z P -z Q )} by the set {x, 2 κα0 κ00 -x}. Note that it implies the choice of a normalisation for θ 0 (z P + z Q ): in fact, in the case that x = θ α (z P + z Q ) we have chosen θ 0 (z P + z Q ) = 1 and in the case that x = θ α (z P -z Q ), we have chosen θ 0 (z P -z Q ) = 1. These choices are made possible by the hypothesis θ 0 (z

P + z Q )θ 0 (z P -z Q ) = 0.
We can then compute the couple

{γ i x + δ i , 2 κi0 κ00 -γ i x -δ i } in A representing the set {θ i (z P + z Q ), θ i (z P -z Q )}
by solving the linear system with coefficients in A:

1 κ00 καα x κα0 κ00 -x θ i (z P -z Q ) θ i (z P + z Q ) = κ i0 κ iα . ( 8
)
The algorithm NormalAdd 2 is then defined as follows: it takes as input 0 = (θ i (0)) i∈Z( 2) 2) and outputs the polynomial defining the algebra A and the

, P = (θ i (z P )) i∈Z(2) and Q = (θ i (z Q )) i∈Z(
γ i x + δ i for i ∈ Z(2); defining the set {(θ i (z P + z Q )) i∈Z(2) , (θ i (z P -z Q )) i∈Z(2) }.
It is clear that the output of NormalAdd 2 determines the set {P + Q, P -Q} and that this algorithm only requires to compute a fixed number of operations in the base field.

We illustrate (the generic version of) this algorithm with g = 1 and n = 2 in Algorithm 3. Let E be an elliptic curve defined by Ω ∈ H 1 ; a point z of the Kummer line associated to E is represented by the projective coordinates

(θ i (z)) i∈Z(2) = (ϑ 0 i/2 (z, Ω/2)) i∈Z(2) . We let (a, b) = (θ 0 (0), θ 1 (0)) = (ϑ [ 0 0 ] (0, Ω/2), ϑ 0 1/2 (0, Ω/2). We define A = 2(a 2 + b 2 ) and B = 2(a 2 -b 2 ). If (a ′ , b ′ ) = (ϑ [ 0 0 ] (0, Ω), ϑ 1/2 0 (0, Ω))
is the coordinate of the dual theta null point, by the duplication formula we have that A = 4(a ′ ) 2 and B = 4(b ′ ) 2 .

Algorithm 3: Normal addition algorithm in genus 1 level 2 input : The points P = (x0, x1), Q = (y0, y1) on E output : The set {P + Q, P -Q}.

1 κ00 = (x 2 0 + x 2 1 )(y 2 0 + y 2 1 )/A + (x 2 0 -x 2 1 )(y 2 0 -y 2 1 )/B; 2 κ11 = (x 2 0 + x 2 1 )(y 2 0 + y 2 1 )/A -(x 2 0 -x 2 1 )(y 2 0 -y 2 1 )/B; 3 κ01 = x0x1y0y1/ab; 4 return {(κ00, X), (κ00, κ10 -X)} where X is a root of X 2 -κ10X + κ00κ11.

Differential additions in level 2

Differential additions are easier to handle: the algorithm outlined in Section 3.2 also works in level 2 with general points.

We illustrate this algorithm with our previous g = 1 and n = 2 example in Algorithm 4. Since this is the addition formula that will be used the most in the pairings algorithm, we give a factored form which may be more convenient for the computations. Note that when using a differential addition to compute ScalarMult, we can always choose affine lifts so that a = x 0 = y 0 = 1 in the notations of Algorithm 4. For the pairing algorithms that we will use, it is easy to see that we can also always multiply the formulas by the same (rational) value. This mean that we can replace the values (A, B) by (1, B/A).

Algorithm 4: Differential Addition Algorithm in genus 1 level 2 input : The points P = (x0, x1), Q = (y0, y1) on E, and P -Q = (z0, z1) with z0z1 = 0. output : The point P + Q = (t0, t1).

1 t ′ 0 = (x 2 0 + x 2 1 )(y 2 0 + y 2 2 )/A; 2 t ′ 1 = (x 2 0 -x 2 1 )(y 2 0 -y 2 1 )/B; 3 t0 = (t ′ 0 + t ′ 1 )/z0; 4 t1 = (t ′ 0 -t ′ 1 )/z1; 5 return (t0, t1)
In genus 2, the corresponding formulas are described in [START_REF] Gaudry | Fast genus 2 arithmetic based on Theta functions[END_REF]. To handle the non generic case, one can use NormalAdd 2 to compute {P + Q, P -Q}; since we know P -Q, it is easy to recover P + Q.

Compatible additions in level 2

Using the NormalAdd 2 algorithm, it is also possible to recover P + R from the knowledge of P, Q, R, P + Q, Q + R on a Kummer variety when Q is not a point of two torsion. We call this operation compatible addition and we note P + R = CompatAdd(P, Q, R, P + Q, Q + R). The idea is the following: NormalAdd 2 (P, R) gives the set {P + R, P -R} and we want to be able to identify P + R in this set. We can suppose that 2P = 0 and 2R = 0 because if this is not the case we have P + R = P -R in the Kummer variety. We remark that NormalAdd

2 (P + R, P + Q) determines {2P + Q + R, Q -R} and NormalAdd 2 (P -R, P + Q) gives the set {2P + Q -R, Q + R}.
As we know Q + R and as 2P = 0 and 2Q = 0 by hypothesis, we have a way to distinguish P + R from P -R.

We now describe an algorithm to compute the compatible addition. Let P α be the defining polynomial of the algebra A = K[X]/(P α ) = K[x] given by NormalAdd 2 (P, R). We want to find a root of P α from the knowledge of P + Q and Q + R without computing a square root in K. We use for all i ∈ Z(2), θ i (z P ±R ) = γ i x + δ i as a convenient notation. Then we define κ(x) ij as the elements of the algebra A given by equation ( 7) where we have replaced z P by z P +Q and z Q by z P ±R . We can suppose that κ(x) 00 = 0 if necessary by changing this index and we consider the polynomials

Q i (Y ) = Y 2 -2 κ(x)i0 κ(x)00 Y + κ(x)ii κ(x)00 for i ∈ Z(2)
. By representing an elements of A as a couple a elements of K via the natural evaluation morphism at the roots of P α (X), the roots of the Q i (Y ) are given by the couples

{ θi(z 2P +Q+R ) θ0(z 2P +Q+R ) , θi(z Q-R ) θ0(z Q-R ) } ,{ θi(z 2P +Q-R ) θ0(z 2P +Q+R ) , θi(z Q+R ) θ0(z Q+R ) }. As consequence, if we evaluate Q i (Y ) at θi(z Q+R ) θ0(z Q+R )
, we obtain a non invertible element of A which is non nul since 2P = 0. As θi(z Q+R ) θ0(z Q+R ) is a coordinate of P + Q, we are done. In the following, we use the notation P + R = CompatAdd(P, Q, R, P + Q, Q + R) for the preceding algorithm. It is clear from our discussion that this algorithm only requires a fixed number of operations in the base field.

We illustrate Compatible additions once more with our g = 1 and n = 2 example in Algorithm 5.

Algorithm 5: Compatible additions in genus 1 level 2 input : x, y, Y = x + z, X = y + z. output : Z = x + y. 1 Computing x ± y; 2 α = (x 2 0 + x 2 1 )(y 2 0 + y 2 1 )/A; 3 β = (x 2 0 -x 2 1 )(y 2 0 -y 2 1 )/B; 4 κ00 = (α + β), κ11 = (α -β); 5 κ10 := x0x1y0y1/ab; 6 Computing (x + z) ± (y + z); 7 α ′ = (Y 2 0 + Y 2 1 )(X 2 0 + X 2 1 )/A; 8 β ′ = (Y 2 0 -Y 2 1 )(X 2 0 -X 2 1 )/B; 9 κ ′ 00 = α ′ + β ′ , κ ′ 11 = α ′ -β ′ ; 10 κ ′ 10 = Y1Y2X1X2/ab; 11 return x + y = [κ00(κ10κ ′ 00 -κ ′ 10 κ00), κ10(κ10κ ′ 00 -κ ′ 10 κ00) + κ00(κ11κ ′ 00 -κ ′ 11 κ00)]
3.6 Three way additions (in level 2 or more)

The three way addition ThreeAdd is an important ingredient to compute the Miller functions f λ,µ,n,P of Lemma 4 below. From P , Q, R, P + Q, P + R, Q + R on a Kummer variety given in the coordinate system provided by level 2 theta functions, it allows to compute P + Q + R. It is also useful when the level n is greater than two: given any affine lifts P , Q, R, P

+ Q, P + R, Q + R in A , it allows to compute a "compatible lift" P + Q + R ∈ A .
In the following, we say that a property is true for R ∈ A (C) (resp. z ∈ C g ) a general point of A (resp. of C g ) if it is verified for all R (resp. z) taken in a Zariski dense subset of A (resp. in π -1 (U ) where U is a Zariski dense subset of A ). For instance, if f ∈ K(A ), we say that we can evaluate f at a general point of A , if we have an algorithm to compute f (R) for all R ∈ U , U a Zariski dense subset of A . Then we can state the proposition: Proposition 1. Suppose that 2 | n and let z 1 , z 2 ∈ C g . For a general point z ∈ C g , suppose that we are given the (affine) theta coordinates of level n of z 1 , z 2 , z, z 1 + z 2 , z 1 + z, z 2 + z. Then one can recover the affine coordinates

(θ i (z 1 + z 2 + z) i∈Z(n) ) of z 1 + z 2 + z. Furthermore, if 4 | n, then the Proposition holds for any z ∈ C g . Proof. For i, j, k, l ∈ Z(2n) such that i + j, i + k and i + l ∈ Z(n), χ ∈ Ẑ(2), let: L 1 (z, i, j, χ) = η∈Z(2) χ(η)θ i+j+η (z + z 1 + z 2 )θ i-j+η (z 1 ), L 2 (z, k, l, χ) = η∈Z(2) χ(η)θ k+l+η (z 2 )θ k-l+η (z), L 3 (z, i, k, χ) = η∈Z(2) χ(η)θ i+k+η (0)θ i-k+η (z 2 + z), L 4 (z, j, l, χ) = η∈Z(2) χ(η)θ j+l+η (z 1 + z)θ j-l+η (z 1 + z 2 )
Then, by Theorem 1, we have

L 1 (z, i, j, χ).L 2 (z, k, l, χ) = L 3 (z, i, k, χ).L 4 (z, j, l, χ). (9)
As 2 | n, we know that the linear system associated to Θ n has no base point. Thus, if we fix i,j with i + j ∈ Z(n), we can choose k, l such that i + k and i + l ∈ Z(n) and θ k+l (z 2 ) = 0. Suppose that L 2 (z, k, l, χ) = 0 for z a general point of C g . Then, L 2 (z, k, l, χ) is a non trivial degree 1 relation between theta functions. But such a relation does not exist since it is known that (θ i (z)) i∈Z(n) form a basis of the vector space of L Ω -quasi-periodic functions of weight n [25, Proposition 1.3 p. 122]. So we can assume (for a general point z of C g ) that for all χ, L 2 (z, k, l, χ) = 0.

We can then compute L 1 (z, i, j, χ) for a general z ∈ C g and for all χ ∈ Ẑ(2). By summing over the characters, we can thus compute for all i, j ∈ Z(2n) such that i + j ∈ Z(n) the products θ i+j (z + z 1 + z 2 )θ i-j (z 1 ). Using again the fact that the linear system associated to Θn has no base point, we obtain that we can compute

(θ i (z + z 1 + z 2 )) i∈Z(n) for z a general point of C g . If 4 | n, to show that we can compute L 1 (z, i, j, χ) for any z ∈ C g , we have to show that L 2 (z, k, l, χ) = 0 for some k, l ∈ Z(n) such that i + j + k + l ∈ 2Z(n).
By the duplication formula [15, Theorem 2, p. 139-141], with a slight abuse of notations, we get that

L 2 (z, k, l, χ) = 2 g θ χ 2 k+l 2n (z 2 + z, 2 Ω n )θ χ 2 k-l 2n (z 2 -z, 2 Ω n ). But by [22, Result a) p. 340], for any a, b ∈ Z and z ∈ C g , there exists c ∈ Z such that θ a 2 b 2n + c (z, 2 Ω n ) = 0. As for any k, l ∈ Z(n), we can find k ′ , l ′ ∈ Z(4) such that L 2 (z, k + k ′ + l ′ , l + k ′ -l ′ ,
χ) = 0 and we are done.

Definition 1. Let P , Q, R, P + Q, P + R, Q + R in A be given by their theta functions of level n. Let P + Q + R ∈ A be the point computed using the relations of the proof of Proposition 1.

We note

P + Q + R = ThreeAdd( P , Q, R, Q + R, P + R, P + Q).
From the way the three addition is computed, we immediately get the multiplicative action of Section 3.2.

Lemma 2. Let T = ThreeAdd( P , Q, R, Q + R, P + R, P + Q). Then ThreeAdd(α * P , β * Q, γ * R, A * Q + R, B * P + R, C * P + Q) = ABC αβγ * T .
We illustrate Three way additions with our g = 1, n = 2 example in Algorithm 6.

Algorithm 6: Three way addition in genus 1 level 2

input : The points x, y, z, X = y + z, Y = x + z, Z = x + y on E. output : T = x + y + z. 1 return T0 = (aX0 + bX1)(Y0Z0 + Y1Z1) x0(y0z0 + y1z1) + (aX0 -bX1)(Y0Z0 -Y1Z1) x0(y0z0 -y1z1) T1 = (aX0 + bX1)(Y0Z0 + Y1Z1) x1(y0z0 + y1z1) - (aX0 -bX1)(Y0Z0 -Y1Z1) x1(y0z0 -y1z1)

A generalisation of Miller's algorithm

In this section, we first recall the classical Miller's algorithm to compute the function of an elliptic curve defined up to a constant factor by a principal divisor, then we explain how to generalize this algorithm to all abelian varieties. Let E be an elliptic curve and denote by 0 the zero element of the group law of E . Let P be a point of E of order ℓ. Then by [27, Corollary 3.5 p. 67] the divisor ℓ(P ) -ℓ(0) is principal and we denote by f an element of K(E ) defined up to a constant factor by (f ) = ℓ(P ) -ℓ(0). More generally, let f λ,P ∈ K(E ) be defined up to a constant factor by

(f λ,P ) = λ(P ) -([λ]P ) + (λ -1)(0).
Miller's algorithm [START_REF] Victor | The Weil pairing, and its efficient calculation[END_REF] to compute f is based on the following remark. For λ > 0 an integer, we have:

f λ+µ,P = f λ,P f µ,P f λ,µ,P ,
where f λ,µ,P is a function associated to the divisor

([λ + µ]P ) -([λ]P ) -([µ]P ) + (0),
assuming a suitable normalisation of the functions. If E is given by a Weierstrass equation, f λ,µ,P can be computed easily from the usual cord and tangent definition of the group law on E . As f 1,P is a constant function, we obtain an efficient square and multiply algorithm to compute f ℓ,P . Now, let A be a principally polarised abelian variety defined by a period matrix Ω. We recall that we denote by Θ the theta divisor of A defined by Ω which is the zero divisor of the associated Riemann theta function and by Θn the theta divisor of level n (as in Section 2). Denote by 0 the zero element of A and let P be a point of order ℓ. For λ > 0 an integer, by an easy induction using the theorem of the square [24, Corollary 4, p. 67], we see that the divisor

λτ * P Θn -τ * λP Θn + (λ -1)Θ n (10)
is principal. We denote by f λ,n,P or more simply by f λ,P if n = 1 the function defined up to a constant factor by this divisor. In the following, we will use the same notation to denote f λ,n,P or its pullback by the projection π : C g → A .

As noted, f λ,n,P is defined only up to a constant factor. This will not present a problem since we will only evaluate this function on degree zero cycles. In practice, one usually take a unique representative in the class of functions defined up to a constant factor by imposing a normalizing condition at a point (for instance at 0). Also the definition of f λ,n,P depends on the divisor Θn. The divisor of any theta function of level n is linearly equivalent to Θn, but the corresponding function f λ,n,P is be different. However, we will see in Section 5 that the definition of the Weil and Tate pairings depends only on the equivalence class of Θn. Lemma 3. Let z P ∈ C g be such that P = π(z P ). For all i ∈ Z(n), we have up to a constant factor:

f λ,n,P (z) = θ i (z) θ i (z + λz P ) θ i (z + z P ) θ i (z) λ , (11)
where for i ∈ Z(n), θ i (z) = θ 0 i/n (z, Ω/n) are the theta functions of level n.

Proof. Denote by g λ,n,P the function with domain C g given by the right hand of [START_REF] Gaudry | The arithmetic of characteristic 2 Kummer surfaces and of elliptic Kummer lines[END_REF]. First, using the quasi-periodicity properties of theta functions with characteristics (2), we obtain that for all µ, ν ∈ Z g and for all z ∈ C g we have

g λ,n,P (z + µ + νΩ) = g λ,n,P (z).
As a consequence, g λ,n,P descends to a well defined function of A . As for all i ∈ Z(n), the zero divisor of θ i (z) is linearly equivalent to Θn, it is clear that the zero divisor of g λ,n,P is (10) (up to replacing Θn by an equivalent divisor).

Keeping the notations of 3, we deduce by an immediate computation using formula (11) that: Lemma 4. For all λ, µ positive integers, we have up to a constant factor f λ.µ,n,P = f λ µ,n,P .f λ,n,µ.P .

We also have the following relation

f λ+µ,n,P = f λ,n,P f µ,n,P f λ,µ,n,P , (13)
where f λ,µ,n,P is a function associated to the divisor τ * (λ+µ)P Θn -τ * λP Θn -τ * µP Θn + Θn. The function f λ,µ,n,P is uniquely defined if we impose that it is normalized on a point. From the definition, we have by using Lemma 3 Lemma 5. Let z P ∈ C g be such that P = π(z P ). Let i ∈ Z(n). For all λ, µ positive integers, we have up to a constant factor

f λ,µ,n,P (z) = θ i (z + λz P )θ i (z + µz P ) θ i (z + (λ + µ)z P )θ i (z) (14)
(for all z where f λ,µ,n,P is defined).

We can now explain how the various addition algorithms presented in Section 3 allow us to compute a normalized version of the function f λ,µ,n,P (z) on a general point.

Proposition 2. Suppose that 4 | n and let λ, µ ∈ N. Suppose that we are given θ i (λz P ), θ i (µz P ) and θ i (0) for i ∈ Z(n). Suppose that we can evaluate the functions θ i (z), θ i (z + λz P ) and θ i (z + µz P ) for all i ∈ Z(n) at a point z ∈ C g . Then we can evaluate the projective coordinates (θ i (z + (λ + µ)z P )) i∈Z(n) .

If 2 | n, the Proposition also holds for a general point z ∈ C g .

Proof. If we had the affine coordinates (θ i ((λ+µ)z P )) i∈Z(n) then by Proposition 1 one could recover the affine coordinates (θ i (z + (λ + µ)z P )) i∈Z(n) using the three way additions. But by Lemma 2, if we can only compute the projective coordinates (θ i ((λ + µ)z P )) i∈Z(n) , then the three way addition gives us the projective coordinates of (θ i (z + (λ + µ)z P )) i∈Z(n) . If 4|n, one can then use NormalAdd to compute (θ i (λ + µ)z P ) i∈Z(n) from (θ i (λz P )) i∈Z(n) and (θ i (µz P )) i∈Z(n) . In the case that n = 2, one need to use the CompatAdd algorithm instead to recover (θ i (λ + µ)z P ) i∈Z(n) from the knowledge of (θ i (λz

P )) i∈Z(n) , (θ i (µz P )) i∈Z(n) , (θ i (z)) i∈Z(n) , (θ i (z + λz P )) i∈Z(n) , (θ i (z + µz P )) i∈Z(n) for z a general point of C g .
When the level n is divisible by 4, we can use this proposition to evaluate a normalized function f λ,µ,n,P .

Corollary 1. Let 4 | n, Q ∈ A , and let R ∈ A be a point such that neither R or Q + R is a pole or zero of the divisor of f λ,µ,n,P . Then from the knowledge of λP and µP , we can evaluate f λ,µ,n,P (Q + R)/f λ,µ,n,P (R).

Proof. We fix affine lifts λP and µP of λP and µP . We compute (λ + µ)P using NormalAdd algorithm and chose an affine lift (λ + µ)P . For a point X = Q+R, R, we compute X + λP , X + µP using NormalAdd and choose any affine lifts X, X + λP and X + µP . Using the three way add, we get an affine lift X + λP + µP of X + λP + µP .

For T ∈ {λP, µP, X, (λ + µ)P, λP + X, µP + X, X + (λ + µ)P }, let z T ∈ C g be such that T = π(z T ) and let α T ∈ K be such that for all i ∈ Z(n), θ i (z T ) = α T * T i (where * is the multiplicative action described in Section 3.1).

By Lemma 2, we have

α X+(λ+µ)P = α X+λP α X+µP α X α (λ+µ)P α λP α µP .
In particular, for i ∈ Z(n), the quotient ( X+λP )i( X+µP )i ( X+(λ+µ)P )i( X)i (if defined) does not depend on choice of an affine lift for X, X + λP and X + µP , but only on the choises of λP , µP and (λ + µ)P . By applying that with X = Q + R et X = R we obtain that the function f λ,µ,n,P evaluated at the cycle (Q + R) -(R) is given by

( Q + R + λP ) i ( Q + R + µP ) i ( Q + R + (λ + µ)P ) i ( Q + R) i ( R + (λ + µ)P ) i ( R) i ( R + λP ) i ( R + µP ) i , for i ∈ Z(n) such that this fraction is defined.
If we go back to the definition of f λ,n,P given by Lemma 3, the method of [START_REF] Lubicz | Efficient pairing computation with theta functions[END_REF] provides us with another way to compute it without going through the f λ,µ,n,P functions. We assume in the Proposition that f λ,n,P is well defined on the cycle (Q) -(0) (as we will see in Section 5 this is usually the case), we leave to the reader the easy adaptation to make to evaluate it on the cycle (Q + R) -(R).

Proposition 3. Let P , Q and 0 be affine lifts of P ,Q and 0. Let P + Q be a lift of P + Q (if 4 | n, we can compute it with a normal addition, otherwise we have to assume it is given). Note Q + λP = ScalarMult(λ, P + Q, P , Q, 0) and λP = ScalarMult(λ, P , P , 0, 0). We have

f λ,n,P (z Q ) f λ,n,P (0) = Q i . λP i ( Q + λP ) i . 0 i ( P + Q) i . 0 i Q i . P i λ (15)
Proof. Let α, β, γ, δ ∈ K * be such that α * (

P + Q) = (θ i (z P + z Q )), β * P = (θ i (z P )), γ * Q = (θ i (z q )) and δ * 0 = (θ i (0)).
By definition, we have

f λ,n,P (z Q ) f λ,n,P (0) = θ i (z Q )θ i (λz P ) θ i (z Q + λz P )θ i (0) θ i (z Q + z P )θ i (0) θ i (z Q )θ i (z P ) λ (16)
Now using ( 4) and ( 5), we see that the right hand of ( 15) is equal to the right hand of ( 16) up to the factor

γβ λ 2 δ λ(λ-1) α λ β λ(λ-1) γ λ-1 δ λ 2 -1 δ αδ βγ λ = 1.
Since almost all known variations of pairing computation algorithms use the Miller's functions f λ,n,P and f λ,µ,n,P , we see that we can extend them to all abelian varieties with Corollary 1, at least if the level is divisible by 4. In the following, we make explicit how to transcribe the Weil, Tate, ate and optimal ate pairings, explain some optimizations and work out the case of level two.

The Weil and Tate pairings

In this section, we recall the definition of the Weil and Tate pairings in the general context of abelian varieties. There are several definitions of the Weil pairing leading to different formulas with their own interest in regard to algorithmic applications. Most of the proofs of the equivalence between these definitions rely on Weil's reciprocity theorem. We explain that a generalisation due to Lang of the Weil's reciprocity allows to adapt the usual proofs with minor modifications.

The Weil pairing

Let f : A → B be a separable isogeny with kernel L between two abelian varieties defined over k. Then we have the following diagram:

0 L A B 0 0 Â B L 0 f f
The Kernel L is the Cartier dual of L, so that we have a non degenerate pairing e f : L × L → K * .

We can give an explicit description of this pairing.

If Q ∈ L(K), Q defines a divisor D Q on B modulo linear equivalence. Then f (Q) = 0 corresponds to f * D Q so there is a function g Q on A such that f * D Q = P ∈L(K) τ * P D Q = (g Q ).
Then for all P ∈ L(K), as (g Q ) is invariant by translation by P , g Q (x)/g Q (x + P ) is a constant function. Its definition does not depend on the choice of D Q and g Q and we have e f (P, Q) = g Q (x)/g Q (x + P ).

Applying this to the isogeny [ℓ] : A → A , we recover the Weil pairing e W : A [ℓ] × Â [ℓ] → µ ℓ where µ ℓ is the set of ℓ th -roots of unity in K. We suppose that A has a principal polarisation Θ. Composing with the polarization ϕ : A → Â associated to the divisor Θ, we get the Weil pairing as e

W : A [ℓ] × A [ℓ] → µ ℓ .
When A is an elliptic curve, it is well known that the Weil pairing can be computed as

e W (P, Q) = f ℓ,P ((Q) -(0)) f ℓ,Q ((P ) -(0)) . ( 17
)
This result can be proved with Weil's reciprocity theorem. It can be generalized to the case where A is the Jacobian of a curve (see for instance [START_REF] Cohen | Handbook of elliptic and hyperelliptic curve cryptography[END_REF]), which is the usual setting in cryptography, because the points and the group law on A have a convenient representation in term of divisors on the curve. By using Lang's reciprocity theorem [START_REF] Lang | Reciprocity and correspondences[END_REF]Theorem 4], it is possible to obtain similar results in the general context of an abelian variety A with a principal polarisation ϕ : A → Â with minor adaptations of the proofs. To explain this, we denote by S : Z 0 (A ) → A (K), the morphism given by

n i (P i ) → n i P i . If Z = n i (P i ) ∈ Z 0 (A ), we let ϕ(Z) = n i (ϕ(P i )) ∈ Z 0 ( Â ).
The cycle ϕ(Z) defines a line bundle on A associated to the divisor

D Z = n i (τ * Pi Θ -Θ).
If S(Z) = 0, this line bundle is linearly equivalent to 0 by the theorem of the square. This means that D Z is the divisor of a function, defined up to a constant factor, that we denote by f Z (the constant factor will play no role in the following since we only consider evaluations of f Z on degree zero cycles). With these notations, we have

Proposition 4 (Lang reciprocity). Let Z 1 , Z 2 ∈ Z 0 (A ) be such that S(Z i ) = 0 for i = 1, 2. Suppose that Supp(Z 1 ) ∩ (f Z2 ) = Supp(Z 2 ) ∩ (f Z1 ) = ∅, then we have f Z1 (Z 2 ) = f Z2 (Z 1 ). Proof. Let Z ′ 1 = ϕ(Z 1 ) ∈ Z 0 ( Â ).
Via the canonical isomorphism A →  , Z 2 ∈ Z 0 (A ) defines as above a divisor on  and because S(Z 2 ) = 0 this divisor D ′ Z2 is linearly equivalent to 0. Denote by f ′ Z2 ∈ K(  ) a function defined up to a constant factor by D ′ Z2 . Let D P be a divisor associated to a Pointcarré line bundle of  ×A . Then by [START_REF] Lang | Reciprocity and correspondences[END_REF]Theorem 4] applied to the divisorial correspondance given by D P , we have

f Z1 (Z 2 ) = f ′ Z2 (Z ′ 1 )
. In fact it is clear that f Z1 is nothing but D P (Z ′ 1 ) with the notations of [START_REF] Lang | Reciprocity and correspondences[END_REF] and because a Pointcarré bundle parametrising line bundles of A via the first projection over  is a Pointcarré bundle parametrising line bundles of  via the second projection over

A = Â , we have f ′ Z2 = t D P (Z 2 ). To finish the proof, it remains to show that f ′ Z2 (Z ′ 1 ) = f Z2 (Z 1 ) but this is an immediate consequence of the fact that ϕ * (f ′ Z2 ) = f Z2 .
In order to show the usefullness of the proposition, we prove the following theorem by adapting the proof of the same result in the case of elliptic curves given in [START_REF] Enge | Bilinear pairings on elliptic curves[END_REF].

Theorem 2. Let P, Q ∈ A [ℓ].
Let D P and D Q be two cycles equivalent to (P ) -(0) and (Q) -(0). The Weil pairing is given by

e W (P, Q) = f ℓD P (D Q ) f ℓD Q (D P ) . ( 18
)
Proof. Let P 0 (resp. Q 0 ) be a point such that P = ℓP 0 (resp.

ℓQ 0 = Q). For X = P, Q, we have [ℓ] * D X = R∈A [ℓ] (X 0 + R) -(R) and S(ℓD X ) = S([ℓ] * D X ) = 0.
Thus, we can set g X = f [ℓ] * D X . For X = P, Q, it is clear by comparing the divisors that, up to a constant factor, we have g ℓ X = [ℓ] * f D X . Let Z P = (ℓ -1)(P 0 ) + (P 0 -P ) -ℓ(0) ∈ Z 0 (A ) and let h P = f Z P as S(Z P ) = 0. Let H P = R∈A [ℓ] h P (x + R), then by comparing the divisors of the functions we obtain that (up to constant factor) H P = g ℓ P . By applying Lemma 4, we have h P (D Q ) = g Q (Z P ). This gives :

R∈A [ℓ] h P (Q 0 + R) R∈A [ℓ] h P (R) = g ℓ Q (P 0 ) g ℓ Q (0) g Q (P 0 -P ) g Q (P 0 ) . ( 19
)
As e W (P, Q) = g Q (P 0 -P )/g Q (P 0 ) and taking into acount that H P = g ℓ P , we finally obtain:

e W (P, Q) = g ℓ P ((Q 0 ) -(0)) g ℓ Q ((P 0 ) -(0)) = f ℓ D P (D Q ) f ℓ D Q (D P ) ,
where, for the second equality, we use the fact that g ℓ X = [ℓ] * f D X for X = P, Q. It is also straightforward to show that the above formula for the Weil pairing depends only on the class of the cycles (P ) -(0) and (Q) -(0) modulo equivalence and on the Chern class of the Theta divisor Θ. For more details, see [START_REF] Lang | Reciprocity and correspondences[END_REF]Section 6]. It is also easy to prove the bilinearity and non degeneracy as in [START_REF] Enge | Bilinear pairings on elliptic curves[END_REF].

Remark 1. In the elliptic case, one usually take the principal polarization coming from the neutral point (that is the point at infinity). Hence, one cannot evaluate f ℓ,P at the cycle (Q) -(0) since f ℓ,P has a pole there. However, in our case we are taking a polarization coming from the Theta divisor, which usually does not contain 0. For the Jacobian of an hyperelliptic curve, the Theta divisor corresponds to degenerate divisors translated by a theta characteristic corresponding of a choice of odd roots of the Weierstrass function (see [26]), so it contains a point of two torsion that is usually different from 0.

In our case, we don't compute the Weil pairing using the principal polarization coming from the Theta divisor Θ, but we use the polarization coming from Θ n . This mean that we will compute the n-th power of the standard Weil pairing, so we will assume that ℓ is prime to n in order to have a non-degenerate pairing in the rest of the paper. To compute this pairing, we simply replace the function f ℓ((P )-( 0)) used in the definition of the Weil pairing by the function f ℓ,n,P defined in Section 4.

The Tate pairing

In this section, we suppose for simplicity that µ ℓ ⊂ K and that A [ℓ] is rational over K, the reader can consult [START_REF] Cohen | Handbook of elliptic and hyperelliptic curve cryptography[END_REF] for the general case.

Let K be the algebraic closure of K and let G = Gal(K/K) be the absolute Galois group of K. By taking the group cohomology long exact sequence associated to the Kummer exact sequence:

0 → µ ℓ → K * → K * → 0,
and using the fact that H 1 (G, K * ) = 0 by Hilbert 90, we obtain an isomorphism

δ 1 : K * /K * ℓ → H 1 (G, µ ℓ ) = Hom(G, µ ℓ ).
In the same way, from the exact sequence

0 → A [ℓ] → A (K) → A (K) → 0,
we get a morphism

δ 2 : A (K)/[ℓ]A (K) → Hom(G, A [ℓ]).
There exists a bilinear application often referred to as the Tate pairing e T :

A [ℓ]×A (K)/[ℓ]A (K) → K * /K * ℓ such that for (P, Q) ∈ A [ℓ]×A (K)/[ℓ]A (K), e W (P, δ 2 (Q)) = δ 1 (e T (P, Q)).
It is well known in the case that A is an elliptic curve that one can compute the Tate pairing by taking any divisor D linearly equivalent to (Q) -(0) and computing e T (P, Q) = f ℓ,P (D). This fact generalizes to any abelian variety with a principal polarisation.

Theorem 3. Let P, Q ∈ A (F q k ) such that P is a point of ℓ-torsion. Let D P and D Q be two cycles equivalent to (P ) -(0) and (Q) -(0). Then we have S(ℓD P ) = 0 and let f ℓD P be the corresponding function on A . The (non reduced) Tate pairing is given by e T (P,

Q) = f ℓD P (D Q ). (20)
Proof. Let Q 0 ∈ A (K) such that ℓQ 0 = Q. Following the definition of the connection morphism δ 2 , we have

δ 2 (Q) = f where f : G → A [ℓ], σ → Q σ 0 -Q 0 is a co-cycle (in fact a morphism since A [ℓ] is rational over K) representing an element of H 1 (G, µ ℓ ).
By definition of the Weil pairing, we have

e W (Q σ 0 -Q 0 , P ) = g P (Q 0 ) g P (Q σ 0 ) . (21)
On the other side, as [ℓ] * (f P ) = c.(g P ) ℓ where c ∈ K is a constant, we have

g P (Q0) g P (0) ℓ = f P (Q) f P (0) . But then δ 1 (f P ((Q) -(0))) is represented by the co-cycle g : G → g P (Q0) g P (Q σ 0 )
. Comparing this with the preceding equation concludes the proof. Let A be an abelian variety over F q a finite field of characteristic p. In order to use the preceding theory, we have to "lift" the abelian variety A over a field of characteristic 0. For this, we denote by R = W (F q ) the ring a Witt vectors with coefficients in F q and by K the quotient field of R. Denote by F q an algebraic closure of F q and let π be the absolute Frobenius morphism.

An abelian scheme A over R, the special fiber of which is A is said to be a lift of A over R. Of course such a lift is not unique in general. For the rest of the section, we fix an embedding κ : R → C so that we can consider R as a subfield of C and we suppose that A becomes an abelian variety defined over a number field K.

Let ℓ be a prime number different from p. We denote by χ ℓ the characteristic polynomial of the Frobenius morphism acting on the ℓ-adic Tate module of A . We recall (see [START_REF] Mumford | Abelian varieties[END_REF]Theorem 4 p. 206]) that χ ℓ is a degree 2g polynomial and if α i are the roots of χ ℓ then there is a permutation σ of {1, . . . , 2g} such that for i = 1, . . . , g, α σ(i) = q/α σ(2i) .

Denote by

G 1 = A [ℓ] ∩ ker(π -1) the Eigenspace of the Frobenius morphism acting on A [ℓ]. If ℓ divides the cardinal #A (F q ), then G 1 is non trivial. In the same way, we let G 2 = A [ℓ] ∩ ker(π -[q]
) be the Eigenspace associated to the Eigenvalue q, if G 1 is non trivial, then by ( 22), G 2 is also non trivial. Denote by k the embedding degree of ℓ, that is the smallest integer such that ℓ|q k -1 (so that F q k is the smallest extension of F q containing µ ℓ ). We remark that G 2 is defined over F q k . By Hensel, we can lift A (F q ) and G 2 to A , as a set of points defined over R and W (F q k ) respectively.

Reducing the Tate pairing on A modulo p we get the Tate pairing as a non degenerate pairing on the right

e T : G 2 × A (F q )/ℓA (F q ) → F * q k /F * ℓ q k . (23)
If A (F q ) has no point of ℓ 2 -torsion, we can identify A (F q )/ℓA (F q ) with G 1 , so that in particular we have that G 2 has the same rank as G 1 .

Finally, if we assume that A (F q k ) also has no points of ℓ 2 -torsion, then by looking at the Tate pairing over F q k , we get a non degenerate bilinear pairing e T,r :

A [ℓ](F q k ) × A [ℓ](F q k ) → µ ℓ ⊂ F * q k (24)
by computing the reduced Tate pairing as e T (P, Q)

q k -1 ℓ . By definition, if P and Q are points in A [ℓ](F q k ), then the reduced Tate pairing is e W (P, π(Q ′ ) -Q ′ ) where Q ′ is any geometric point with [ℓ]Q ′ = Q.
In practice, we compute the Tate pairing by using the formula of Theorem 3. We set D Q = (Q)-(0) so that f n ℓD P = f ℓ,n,P and compute e T (P,

Q) = f ℓ,n,P (D Q ) if well defined. Otherwise, we replace D Q by the equivalent cycle (Q + R) -(R)
where R is any point in A (F q k ), such that this cycle is not in the support of the divisor associated to f ℓ,D P . Then we have

e T (P, Q) = f ℓ,n,P (Q + R) f ℓ,n,P (R) .
In order to check this directly, consider the cycle (Q + R) -(Q) -(R) + (0). This cycle corresponds to a divisor linearly equivalent to 0; let g be a function associated to it. Lang's reciprocity yields f ℓ,P ((

Q + R) -(Q) -(R) + (0)) = g(ℓ(P ) -ℓ(0)) = (g(P )/g(0)) ℓ ∈ F * ℓ q k . Hence e T (P, Q) = f ℓ,P ((Q + R) -(R)).
We remark as in Section 5.1 that since we use a non principal polarisation given by Θ n , the Tate pairing that we compute is also equal to the usual Tate pairing to the power of n. Remark 2. If A is the Jacobian of a curve C, then the points P and Q can be seen as divisors D P and D Q on C. Lichtenbaum showed in [START_REF] Lichtenbaum | Duality theorems for curves over p-adic fields[END_REF] that the Tate pairing can be computed directly by evaluating f ′ ℓD P (D Q ) where f ′ ℓD P is a function in the function field of C.

This has the advantage that the Lichtenbaum pairing only uses functions defined over the curve rather than the Jacobian. However, working directly on the curve has the drawback that twists of curves are not easy to relate explicitly to twists of the Jacobian, hence it is hard to generalize the twisted ate pairing in higher genus [START_REF] Granger | Ate pairing on hyperelliptic curves[END_REF]. Since our point of view is to consider pairing on abelian varieties, it is easy for us to extend all results of [START_REF] Hess | The eta pairing revisited[END_REF] to higher dimension.

Computing the Weil and Tate pairings

Let P ∈ A [ℓ], we are going to present two methods to compute the normalized functions f ℓ,n,P /f ℓ,n,P (0). More precisely, suppose that we are given P = (P i ) i∈Z(n) and let Q = (Q i ) i∈Z(n) two geometric point in A embedded in P Z(n) via the theta coordinates of level n. Let z P , z Q ∈ C g be such that P = (θ i (z P )) and Q = (θ i (z Q )), we want to compute f ℓ,n,P (z Q )/f ℓ,n,P (0) from the knowledge of the homogeneous coordinates of the points P , Q and 0. By Section 5.1, this is sufficient to compute the Weil and Tate pairings (for the Weil pairing one has to repeat the procedure swapping P and Q).

Lemma 4, Lemma 5 and Corollary 1 give a first method we compute Weil and Tate pairings in a similar way to Miller's algorithm when the level n is divisible by 4. Take a random point R ∈ A (K). Let Q ∈ A [ℓ]; starting from the function f 1,n,P = 1, we can use relation [START_REF] Hess | Pairing lattices. Pairing-Based Cryptography-Pairing[END_REF] and Corollary 1 to compute by a square and multiply algorithm f ℓ,n,P (Q + R) and f ℓ,n,P (R). The algorithm terminates when Q + R and R do not belongs to the poles of the functions f λ,µ,n,P used in the computation. The Weil pairing can be computed in a similar manner. As explained in the introduction, all these results are valid when K is a finite field of characteristic different from 2 and give a probabilistic algorithm with expected probability arbitrarily close to 1 to output the result as the size of the K grows to infinity.

A second method is to use Proposition 3 directly. We will see in Propositions 5 and 6 that in fact we can always compute the Weil or Tate pairing.

Proposition 5. Let P + Q, P , Q and 0 be affine lifts of P +Q, P ,Q and 0, where P is a point of ℓ-torsion. Let λ 1 P , λ 0 P be such that ScalarMult(ℓ, P + Q, P , Q, 0) = λ 1 P Q and ScalarMult(ℓ, P , P , 0, 0) = λ 0 P 0. We have

e T (P, Q) = λ 0 P λ 1 P . ( 25
)
If Q is also a point of ℓ-torsion, then we can define λ 1 Q , λ 0 Q in a similar manner, and we have

e W (P, Q) = λ 0 P λ 1 P λ 1 Q λ 0 Q . ( 26
)
Proof. Immediate by Proposition 3. See also [START_REF] Lubicz | Efficient pairing computation with theta functions[END_REF].

Remark 3. Since P is a point of ℓ-torsion, the cycle ℓ(P ) -(ℓP ) -(ℓ -1)(0) is equal to the cycle (ℓ + 1)(P ) -([ℓ + 1]P ) -ℓ(0). In particular, we can also compute the Tate (and Weil) pairing by using f ℓ+1,n,P rather than f ℓ,n,P . In the context of Proposition 5, this mean that we have ScalarMult(ℓ + 1, P + Q, P , Q, 0) = λ 1 P P + Q, ScalarMult(ℓ + 1, P , P , 0, 0) = λ 0 P P and e T (P, Q) = λ 0 P /λ 1 P .

The fact that we can always compute the Weil and Tate pairings in Proposition 5 can be explained as follow. For P ∈ A [ℓ] and i ∈ Z(n), let f ℓ,n,P,i be the function whose divisor is ℓτ * P Θn,i -τ * ℓP Θn,i + (ℓ -1)Θ n,i where we make explicit the dependency with respect to the divisor Θn,i corresponding to the coordinate and we can always find i, j ∈ Z(n) such that this expression is well defined.

It is well known (see for instance [START_REF] Granger | Ate pairing on hyperelliptic curves[END_REF]) that in the case that (0) is in the pole of the function f ℓ,n,P then the Tate pairing can be defined as f ℓ,n,P (Q)/c where c correspond to the leading coefficient of f ℓ,n,P in the completion of Θ n along A . Equation (27) can be seen as a version of this since the leading coefficient of f ℓ,n,P,i is the same as f ℓ,n,P,j and the latter can be obtained as θj (ℓz P ) θj (0) (up to a ℓ th -power) when (0) is not in the pole of Θ n,j .

But the fact that we can always compute the Tate pairing can be seen as a misfeature in term of complexity. It means that computing ScalarMult(ℓ, P + Q, P , Q, 0) = λ 1 P Q amount to computing the n g Miller's functions f ℓ,n,P,j ((Q)-(0)). It would be interesting to know if we could compute only one such function f ℓ,n,P,0 ((Q) -(0)) from ScalarMult(ℓ, P , P , 0, 0) and the coordinates of Q (and eventually P + Q).

We can always use Proposition 5 when 4 | n to compute the Weil and Tate pairings. In level n = 2, we need the point P + Q. For instance, we can work at first in level 4 and then switch to level 2 (using the duplication or the isogeny formula) to compute the pairing once we have P + Q.

In general, if we start from P and Q in level 2, we can only compute the symmetric Weil (or Tate) pairing e W (P, Q) + e W (-P, Q) as defined in [20, Section 5.2]. To do this, we compute P + Q in the algebra A from Section 3, and do all the computations in this algebra of dimension 2. In practice, it seems faster to take a square root to fix a choice of P + Q in the field of definition K and do all the remaining computation in K rather than in A.

In a sense, the second method, by computing differential additions on affine lifts, can be thought of as using the definition of the Weil pairing as the commutator pairing of the theta group [START_REF] Mumford | On the equations defining abelian varieties[END_REF] (or in the analytic version, as the symplectic pairing induced by the Riemann form [START_REF] Mumford | Tata lectures on theta I[END_REF]). But if we unravel the first version using the classical Miller algorithm with theta functions, by looking at Lemma 3, we see that Proposition 2 actually is just another way to compute an affine scalar multiplication. Proposition 6. Let P + Q, P , Q and 0 be affine lifts of P + Q, P ,Q and 0, where P is a point of ℓ-torsion. Use the following algorithm to compute lifts ℓP and Q + ℓP by a double and add method:

Input: Affine lifts λP and Q + λP . Double: (at each step.) Compute 2λP and Q + 2λP with two differential additions. Add: (only if the current bit of ℓ is one.) From the points 2λP , P , Q + P and Q + 2λP use a compatible addition to compute the projective point (2λ + 1)P (if 4 | n one can of course use a normal addition directly); and take an arbitrary lift.

From the affine lifts 2λP , Q, P , Q + P , (2λ + 1)P and Q + 2λP , do a three way addition to compute an affine lift Q + (2λ + 1)P .

Let λ 1 P , λ 0 P be such that Q + ℓP = λ 1 P Q and ℓP = λ 0 P 0, then we have

e T (P, Q) = λ 0 P λ 1 P . ( 28
)
And a similar result holds for the Weil pairing. In particular, we can always compute the Tate and Weil pairing when 4 | n.

Proof. This is a direct application of Proposition 2 and Lemma 5. Indeed, by Corollary 1, the result does not depend on the choice of lift of (2λ + 1)P done after each addition step.

Q, we have first to compute the coordinates of the point P + Q. This can be done in two ways: either we work with a system of coordinates where we can perform group law addition, like theta functions of level 4 or Mumford coordinates if A is the Jacobian of an hyperelliptic curve. After computing P + Q in this system of coordinates, we can convert P , Q and P + Q to theta coordinates of level 2 (if working with theta functions of level four, we can convert to level 2 by using the duplication formula, or by using the isogeny formula, which is essentially free since the level two coordinates of the isogenous points are a subset of the level four coordinates of the original points). If we only have the points P and Q in theta coordinates of level 2, then we can use the results of Section 3.3 to compute P ± Q either by computing with a degree 2 an algebra over the field of definition of the points or taking a square root (if q is congruent to 3 modulo 4 this can be done by an exponentiation). In this case, we can only compute the symmetric pairing. In all cases, to be able to work with theta coordinates of level 2, we need the theta null point of level 2 to be rational.

The main loop of the different algorithms rest in ScalarMult to compute the functions f λ,n,P . Each step of the loop will then consists in a doubling, and two differential additions with the same point Q as difference. In the optimal ate pairings, there is some additional computations required to obtain the functions f λ,µ,n,P , but there is at most k -1 such functions to compute, while the length of the loop is in log(r 1/ϕ(k) ), so we can safely ignore their contribution to the running time.

This justify that we focus on the complexity of one step in the evaluation of the Miller function f λ,µ,n,P (Q), according to whether P and Q is defined over the "big field" F q k or the base field F q . We denote by M a multiplication in F q k , S a square in F q k , m a multiplication by a "constant" in F q k coming from the coordinate of P , Q or P + Q. The corresponding operations over F q are denoted by M , S and m respectively, and we denote by m 0 a multiplication by a constant depending only on the abelian variety (that comes from the theta null point). Lastly, we let M, m and m 0 be the multiplication between an element of F q k and an element of F q (respectively a constant in F q depending only on the coordinates of P , Q or P + Q, and a constant in F q depending only on the theta null point).

We first focus on the case where P and Q are defined over the big field. This is not the case in cryptography where we take a point in F q to speed up the computations, but this is unavoidable in situation where one has to compute the Weil pairing between points of ℓ-torsion (for instance, to get a symplectic basis of the ℓ-torsion).

In dimension 1, one step of the "Miller loop" using ScalarMult as in Proposition 5 is given in Algorithm 8 (for the Weil pairing, we will need two such loops). We first note that we can always choose the lifts P , Q and P + Q so that x P = 1, x Q = 1, x P +Q = 1. We see that a step takes 5M + 2m + 7S + 2m 0 in dimension 1. We note that we can get a 1S + 1m 0 -1M trade-off by computing

Z n = β(x 2 n -z 2 n ) as Z n = A B (x 2 n -z 2 n ) 2 .
In dimension 2, as in [START_REF] Lubicz | Efficient pairing computation with theta functions[END_REF] we obtain that the cost of one step is 11M + 6m + 13S + 6m 0 . Similar to the dimension one case, there is a possible trade-off of 3S + 3m 0 -3M.

+ P = (x ′ n+1 , z ′ n+1 ). 2 Output 2nQ = (x2n, z2n); (2n + 1)Q = (x2n+1, z2n+1); (2n + 1)Q + P = (x ′ 2n+1 , z ′ 2n+1 ). 1. α = (x 2 n + z 2 n ); β = A B (x 2 n -z 2 n ). 2. Xn = α 2 ; Xn+1 = α(x 2 n+1 + z 2 n+1 ); X ′ n+1 = α(x ′2 n+1 + z ′2 n+1 ); 3. Zn = β(x 2 n -z 2 n ); Zn+1 = β(x 2 n+1 -z 2 n+1 ); Z ′ n+1 = β(x ′2 n+1 + z ′2 n+1 ); 4. x2n = Xn + Zn; x2n+1 = (Xn+1 + Zn+1)/xP ; x ′ 2n+1 = (X ′ n+1 + Z ′ n+1 )/xQ; 5. z2n = a b (Xn -Zn); z2n+1 = (Xn+1 -Zn+1)/zp; z ′ 2n+1 = (X ′ n+1 -Z ′ n+1 )/zQ; 6. Output (x2n, z2n); (x2n+1, z2n+1); (x ′ 2n+1 , z ′ 2n+1 ).
We note that the pairing algorithm relying on theta functions, despite being a very generic algorithm available for all abelian varieties is actually pretty fast in small dimension. As a comparison, just doubling a point using the fastest known arithmetic for Mumford projective coordinates on a Jacobian of a curve of genus 2 (which is a necessary step for the Miller algorithm) already takes 33M + 7S + 1m 0 [START_REF] Lange | Formulae for arithmetic on genus 2 hyperelliptic curves[END_REF]! For elliptic curve, just the Doubling step in the Miller loop of a point on an Edwards curve takes 9M + 7S + 2m 0 (since there is no denominator elimination to compute the Weil pairing). Also, as already noted, if the Hamming weight of ℓ is small, the algorithm from Proposition 6 could be faster than the one from Proposition 5. Now for pairings used in cryptography, one can usually choose which type of points to compute with. For the Tate pairing, this means that P is defined over F q while Q will is defined over F q k . Moreover, when the embedding degree k is even one can do denominator elimination, and in genus 2 one can take for Q a degenerate divisor. Table 2 shows the comparison between the algorithm from Proposition 5 with the usual Miller algorithm, with or without denominator elimination. More precisely, in genus 1 the cost for one step of the Miller loop using theta coordinates is 1m + 2S + 2M + 3M + 1m + 5S + 2m 0 (one can also use the 1S+1m 0 -1M trade-off). In genus 2, the cost is 3m+4S+4M+7M +3m+9S+6m 0 (or with the 3S + 3m 0 -3M trade-off). Here, the algorithm from Proposition 6 is not interesting because it only reduces operations in the smaller field for Doubling, while adding operations in the bigger field for each Addition. As a comparison, a pairing step with Edwards coordinates using denominator elimination costs 1M+1S+1M+6M +5S +2m 0 for each Doubling, and 1M+1S+1M+12M +1m 0 for each Addition. We see that for the Tate pairing, our algorithm performs poorly because while we do less operations in the small field, we compute more in the big field. In genus 1, Table 2 shows that this is because we do not have denominator elimination.

For the ate (and optimal ate) pairing, the Miller loop is shortened, but P lives in F q k while Q lives in F q . Since most operations take place in the big field, we expect our algorithm to be competitive, since it does less computations overall when all points are in the big field. For instance in dimension 1, one step costs 5M + 1m + 7S + 1m + 2m 0 , and in dimension 2 11M + 3m + 13S + 3m + 6m 0 (with the possible 1S + 1m 0 -1M and 3S + 3m 0 -3M trade-offs respectively). In general, the cost of doing more computations in the bigger field offsets the reduced loop, so one usually use ate pairings in presence of twists. If there is a twist of degree d | k, then there is a twist that will send P to a point in F q e (where e = k/d) and Q in F q k . Computing the pairing entirely on the twist then costs about the same as the Tate pairing, but with operations in F q replaced by operations in F q e (see [START_REF] Costello | Faster pairing computations on curves with high-degree twists[END_REF]), but with still the same loop length gain as the ate pairing. Depending on the size of e and k, our algorithm may be competitive in this case (and it could extend the range where the ate pairing is faster than the Tate pairing).

In dimension 2, according to [START_REF] Granger | Ate pairing on hyperelliptic curves[END_REF] the ate pairing using affine Mumford coordinates costs 1I + 29M + 5S + 7M for an addition, and 1I + 29M + 9S + 7M for a doubling, where I denotes the cost of an affine inversion in F q k . Even when using degenerate divisors, the cost is still of 1I + 27M + 3S + 4M for an addition and 1I + 27M + 7S + 4M for a doubling, so our formulas are way faster. In dimension 3, our formula for the ate pairing (without any optimization as settings some projective coordinates to one) will be in 32M + 24S + 8M + 8m + 16m 0 , which is faster than only a doubling step using degenerates divisor using affine Mumford coordinates.

In genus 1, we expect the optimal ate pairing to gain a factor of ϕ(k) in the loop length where k is the embedding degree. We would prefer k to be uneven, since this will increase the size of ϕ(k). But the denominator elimination trick only work in the even case (even if there are some adaptations available in the uneven case). In dimension g > 1, we expect to have ℓ > q, so the reduction µ mod ℓ in Section 6 is equal to q. In this case the ate pairing is already reduced [START_REF] Granger | Ate pairing on hyperelliptic curves[END_REF] so there is no need for the final exponentiation. This means that one can't use a denominator elimination. But since in our algorithm we don't use denominator elimination anyway, these cases are actually favorable for the algorithms presented in this paper. 1. P ∈ A[ℓ](F q ), Q ∈ A[ℓ](F q k ) (counting only operations in F q k ).

Conclusion

The main purpose of the paper is to give an algorithm to compute all known pairings (in particular the optimal ate pairing) on any abelian variety represented by theta functions (of even level). For efficiency reason, a particular focus was given on the level 2 case, which correspond to Kummer varieties rather than abelian varieties which lead to some difficulties. As seen in Section 8, while generic, the algorithm is surprisingly fast in lower dimension. To be truly competitive with the best pairing algorithms used in cryptography, it would be interesting to know if a denominator elimination is possible. It would also be worthwhile to investigate the case of degenerate divisors to speed up the pairing computation.

According to a recent paper [START_REF] Joppe W Bos | Two is greater than one[END_REF], the arithmetic in dimension 2 can be faster than in dimension 1 (for a 128-bit security level). To achieve this speed, they use the representation given by level 2 theta functions. The ability to compute optimal pairings with such functions as explained in this paper, definitively show that the dimension 2 case is worth studying for pairings applications!

  Miller's algorithm to compute the function f of a genus g curve C defined up to a constant factor by a divisor D ∈ Pic 0 (C) linearly equivalent to 0. For D a , D b ∈ Pic 0 (C), denote by f Da,D b the function given (up to a constant factor) by the divisor D a + D b -(D a ⊕ D b ) where D a ⊕ D b is the reduced divisor associated to D a + D b . Miller's algorithm is based on a double and add loop that iterates a formula which gives the function f Da,D b .

Algorithm 2 :

 2 Corollary 4.5.2 and Remark (2)] Differential addition algorithm in genus 1 level 4 input : The points P = (x0, x1, x2, x3), Q = (y0, y1, y2, y3) and ( P -Q)(z0, z1, z2, z3) on E.

Algorithm 8 :

 8 One step of the differential addition Miller loop 1 Input nQ = (xn, zn); (n + 1)Q = (xn+1, zn+1), (n + 1)Q

Compared to the method of Proposition 5, which uses three differential additions at each step (whatever the current bit of ℓ is), the method of Proposition 6 only uses two differential addition at each "doubling" step. However, it requires a compatible addition and a three way addition at each "addition" step. Still, it may be worthwhile to use when the hamming weight of ℓ is small, especially combining with a NAF method.

Ate and symmetric ate pairings

In this section, we give a generalisation of ate pairings to all abelian varieties. We keep the notations of Section 5.2.

For the Tate pairing, one usually take P ∈ G 1 and Q ∈ G 2 for efficiency reasons. However, if P ∈ G 2 (F q k ), then π(f ℓ,n,P ) = f ℓ,n,qP .

(See [START_REF] Hess | The eta pairing revisited[END_REF]Lemma 3]). The idea of the ate pairing is to switch the role of P and Q (that is take P ∈ G 2 and Q ∈ G 1 ) and use this relation in order to decrease the number of iterations in the computation of the Tate pairing. For this we take µ ≡ q mod ℓ. We remark that as ℓ|q k -1, we have ℓ|µ k -1 and we put m = µ k -1 ℓ . Let n be the level of our theta functions. We compute:

The first equality is a consequence of Lemma 4 and the second one comes from the definition of f µ k -1,n,P .

Then by a repeated use of Lemma 4 and (29) we get

We can define the (n-power of the usual) ate paring e A :

By the equations above we have,

where c = k-1 i=0 µ i q k-1-i . The reduced ate pairing is then equal to the reduced Tate pairing to the power of m, so this pairing is non degenerate if ℓ does not divide m q k -1 ℓ . A trick to reduce m is to divide it by q k -1 ∧ µ k -1. Indeed, the (reduced) hℓ Tate pairing computed on a point of ℓ-torsion is equal to its (reduced) ℓ-Tate pairing. Remark 4. If A is an elliptic curve, let t be the trace of the Frobenius morphism. If ℓ = q + 1 -t is a prime number, we can take µ = 1 -t which by the Weil bound is a O( √ q). In this case, the expected number of iteration needed to compute the ate pairing is less than half of those required for the computation of the Tate pairing so that we can expect a speed up.

If A is an abelian variety of dimension greater or equal to 2, we have µ q while ♯A (F q ) = O(q g ), so usually ℓ is greater than q. In this case we gain a g-fold speedup in the number of iteration to compute the Miller function.

Note that in the case µ = q (for instance when g 2), then f µ,n,P ((Q) -(0)) is already reduced so there is no need for the final exponentiation, and we can replace the cycle (Q) -(0) by any equivalent cycle [START_REF] Granger | Ate pairing on hyperelliptic curves[END_REF].

An algorithm to compute e A is provided by Proposition 3. Before the final exponentiation, we need to compute ScalarMult(µ, P + Q, P , Q, 0), ScalarMult(µ, P , P , 0, 0)

As a consequence, these calculations can be done at the cost of O(log µ) iterations using a fixed number of operations in the field F q k . Explicitly: Proposition 7. Let P + Q, P , Q and 0 be affine lifts of P +Q, P ,Q and 0, where

We have

Remark 5. By looking at the differential additions, we can recover the power in which the (non reduced) Tate pairing and ate pairing correspond.

Recall that µ k -1 = mℓ. Suppose that ℓP + Q = α Q and ℓP = β 0, so that α/β gives the Tate pairing. Then mℓP

If we write µ i P + Q = γ ′ i π( µ i-1 P + Q) and µ i P = δ ′ i π( µ i-1 P ), then since π commutes with differential additions we get by Lemma 1

By a trivial recursion, we have

Finally, we can apply the same trick to the symmetric Tate pairing of [START_REF] Lubicz | Efficient pairing computation with theta functions[END_REF] to obtain in the same way a symmetric ate pairing on Kummer varieties.

Twisted ate pairing

We have remarked that in the ate pairing, the length of the Miller loop is less than in a regular Tate pairing. This comes at a cost however. In the ate pairing, we need to compute f µ,n,P , where P ∈ G 2 lives in the "big field" F q k . When doing the Tate pairing, one can instead compute f ℓ,n,P , where P ∈ G 1 lives in the "small field" F q ; in this case only the evaluation f ℓ,n,P (Q) is done in F q k .

In the supersingular case, since the action of the Verschiebung π is inseparable and acts by multiplication by q on G 1 and 1 on G 2 , we can define in a similar manner the eta pairing as f µ,n,P (Q) where P ∈ G 1 and Q ∈ G 2 like the Tate pairing [START_REF] Hess | The eta pairing revisited[END_REF].

In the ordinary case, to alleviate this problem, one usually combines twists with the ate pairing. Suppose that there exists a twist of degree d | k, and let e = k/d. Then this twist correspond to a d th -root of unity ζ in the endomorphism ring. We have

Mapping a point of G 2 via this twist, we get a point rational over F q e (while the points in G 1 are defined over F q k on the twist). So one can use the twist map to compute f µ,n,P in the twist and come back to evaluate at Q. This has the advantage that in the twist, P only lives in the extension F q e . However, to apply this idea to our setting, we need a rational theta structure on both the abelian variety A and its twist, which is impossible since the theta structure rigidify the moduli stack of abelian varieties by [START_REF] Mumford | On the equations defining abelian varieties[END_REF] (at least in level 4).

One can instead compute everything in the twist, as explained in [START_REF] Costello | Faster pairing computations on curves with high-degree twists[END_REF]. In this case, we only need a rational theta structure on the twist, and P is defined over F q e while Q is defined over F q k . (If P ′ and Q ′ are the twisted points, then e A (P ′ , Q ′ ) d = e A (P, Q) d . In [START_REF] Costello | Faster pairing computations on curves with high-degree twists[END_REF], they have a finer result for elliptic curves by analysing the updates of Miller functions).

Another method is to compute the twisted ate pairing e A ′ (P, Q) = f µ e ,n,P ((Q)-(0)) [START_REF] Hess | The eta pairing revisited[END_REF], where P ∈ G 1 and Q ∈ G 2 so that the computation of f µ e ,n,P can be done in the smaller field, but at the cost of a larger loop than the ate pairing. The idea behind the twisted ate pairing is that pulling back the action of the Frobenius of F q e of the twist A ′ to A , we get that ζ i π e acts as an inseparable endomorphism ψ of degree q e on A , with ψ(P ) = [q e ]P and

Optimal and symmetric optimal pairings

We recall the following proposition from [29, Theorem 1], stated for elliptic curves but which adapt easily to abelian varieties, as in Section 6 (from where we take the notations G 1 and G 2 ). Proposition 8. Let λ = mℓ = c i q i be a multiple of ℓ (and such that ℓ ∤ m). The pairing

j>i cj q j ,ciq i ,n,Q (P )

is non degenerate when

The idea is then to find a multiplie λ such that the coefficients c i are small. More precisely, since ϕ k (q) = 0 mod ℓ, one can use LLL to find a small relation among the powers q, q 2 , . . . , q ϕ(k)-1 . The discussion in [START_REF] Vercauteren | Optimal pairings. Information Theory[END_REF]Section 3.3] shows that we can expect to find λ such that c i ≈ ℓ 1/ϕ(k) . Algorithm 7: Optimal ate

1 Take affine lifts P , Q and P + Q; 2 for i = k -1, . . . , i = 0 do 3 Compute P + ci Q and ci Q using ScalarMult; Apply the Frobenius endomorphism to obtain P + ciq i Q and ciq i Q ; 5 Compute ciq i Q + j>i cjq j Q (up to a constant, using NormalAdd) and then use the extended Riemann relations from Proposition 2 to compute P + ciq i Q + j>i cjq j Q (up to the same constant); 6 Find the constants C0 and C1 such that we have λ Q = C0 * Q and P + λ Q = C1 * P ;

Return : (C1/C0) q k -1 ℓ

We can easily compute such a pairing by using the results of Section 4. In Algorithm 7 we give the corresponding algorithm when 4 | n. When n = 2, as in Section 5 we either compute P + Q by working in A all the time, or by taking a square root to get back in K. All the steps of the algorithm are then the same, except for Step 5 where we use CompatAdd on c i q i Q, j>i c j q j Q , P + c i q i Q, P + j>i c j q j Q to compute c i q i Q + j>i c j q j Q.

The extension to pairing lattices [START_REF] Hess | Pairing lattices. Pairing-Based Cryptography-Pairing[END_REF] is also straightforward.

Performance comparisons

For the performance analysis, we only consider the case of level n = 2 since it gives the fastest representation and arithmetic. To compute the pairing between P and


Jean Vandemeulebrouck

P Roux

E Cros

The plumbing of Old Faithful Geyser revealed by hydrothermal tremor

Keywords: Hydrothermal systems, Volcano seismology, Instruments and Techniques

Key points -Sources of hydrothermal tremor are precisely localized using acoustic array processing. -Tracking of tremor sources reveals a lateral cavity that acts as a recharge reservoir. -Oscillating behavior can be explained by fluid compressibility in the recharge cavity.

Introduction

Seismic records recorded around Old Faithful Geyser (OFG) (Fig. 1a) highlighted the remarkably energetic, seismo-acoustic activity that is induced by the continuous cavitation of steam bubbles in the water column [START_REF] Kieffer | Seismicity at Old Faithful Geyser: an isolated source of geothermal noise and possible analogue of volcanic seismicity[END_REF][START_REF] Kedar | The origin of harmonic tremor at Old Faithful Geyser[END_REF]. Concurrent pressure and seismic measurements in the conduit [START_REF] Kedar | Bubble collapse as the source of tremor at Old Faithful Geyser[END_REF], demonstrated that the boiling and cavitation is focused in the upper one meter of the water column in the conduit. Time-dependent localization of this powerful acoustic source during the eruption cycle allows us to track water level and phase separation in the conduit, which in turn provides new insight into the dynamics of the geyser.

The geometry of the Old Faithful Geyser conduit was investigated to a depth of 14 m using a video camera [START_REF] Hutchinson | In-situ observations of Old Faithful Geyser[END_REF]. It consists of a fissure-like channel that is elongated in the eastwards direction. Video recordings show the fissure width varies with depth, and the narrowest part in the conduit, ~10 cm across, is likely responsible for the choked (supersonic) flow during the eruptions [START_REF] Kieffer | Geologic nozzles[END_REF]. High temperatures and adverse flow conditions below 14 m precluded the camera from passing to deeper levels.

Based on the combination of visual observations and seismic, temperature and water level data, [START_REF] Kieffer | Seismicity at Old Faithful Geyser: an isolated source of geothermal noise and possible analogue of volcanic seismicity[END_REF] recognized different stages during the OFG cycle. These stages include a ~15-min-long recharge period just after an eruption characterized by a low water level and weak seismic activity. This period is followed by a ca. 50-min-long period of increasing seismic activity associated with the rise of the water level prior to the next eruption (Fig. 2a).

In the present study, seismic data obtained in 1992 [START_REF] Kedar | The origin of harmonic tremor at Old Faithful Geyser[END_REF][START_REF] Kedar | Bubble collapse as the source of tremor at Old Faithful Geyser[END_REF] with a dense array around OFG (Fig. 1b) are re-processed using acoustic localization techniques [START_REF] Legaz | Self-potential and passive seismic monitoring of hydrothermal activity: a case study at Iodine Pool, Waimangu geothermal valley[END_REF], with the ambition to characterize the spatial and temporal patterns of the seismic sources inside the geyser conduit (Fig. 2b).

Data and Localization Method

In volcano seismology, source localization is generally performed on a series of similar and individual events. However, continuous tremor activity is traditionally encountered in hydrothermal areas with high event rate (e.g., up to 100 events/min of Old Faithful), which prevents discrete seismic events from being isolated. Instead, hydrothermal tremor signals can be analyzed using beamforming techniques [START_REF] Legaz | Self-potential and passive seismic monitoring of hydrothermal activity: a case study at Iodine Pool, Waimangu geothermal valley[END_REF][START_REF] Cros | Locating hydrothermal acoustic sources at Old Faithful Geyser using Matched Field Processing[END_REF], the goal of which is to build a probability map of the dominant acoustic source for a defined time window that might contain a large number of individual events.

A fundamental requirement for beamforming techniques [START_REF] Baggeroer | An overview of matched field methods in ocean acoustics[END_REF] is that the signal is spatially coherent across all the geophones in the array. This implies that the selected time window is long enough to accumulate the signal-to-noise ratios from bubble collapses that occur at the same location, but not too long to prevent spatial decorrelation when the dominant noise source inside the geyser conduit shifts. An estimate of the medium's seismic velocity is also required, as velocity mismatch can degrade the resolution of the source position. This often promotes the application of beamforming techniques to low-frequency signals. In contrast, higher frequencies (short wavelength) result in a better spatial resolution of beamforming techniques. The balance between the resolution at high frequencies and the robustness at low frequencies is site specific, and is determined on a case-by-case basis.

The spatial coherence of the OFG data is optimal between 10 Hz and 15 Hz, corresponding to the lower range of the frequency spectrum that is excited by bubble cavitation in geothermal areas [START_REF] Kedar | Bubble collapse as the source of tremor at Old Faithful Geyser[END_REF][START_REF] Legaz | Self-potential and passive seismic monitoring of hydrothermal activity: a case study at Iodine Pool, Waimangu geothermal valley[END_REF]. Analysis of the seismic response to hammer shots [Kedar et al., 1995] reveals a mean surface velocity of ~130 m/s, slightly increasing with depth. This value was used for the velocity model, in agreement with previous studies [Kedar et al., 1995;[START_REF] Cros | Locating hydrothermal acoustic sources at Old Faithful Geyser using Matched Field Processing[END_REF]

New insights on the geometry and dynamics of Old Faithful Geyser

The 137 min of recorded signals from Old Faithful were centered on a long duration eruption and were processed to locate the time-dependent dominant noise-source position. The seismic tremor signals were truncated into consecutive 20-s-long time windows, which contain at most 40 overlapping bubble collapse events.

High-resolution beamforming techniques produce a ~2 m focal spot that is much smaller than the average acoustic wavelength (λ~10 m) but still not sufficient to provide insight regarding to the diameter of the water conduit (from ~10 cm to ~1 m), especially in the constricted areas [START_REF] Cros | Locating hydrothermal acoustic sources at Old Faithful Geyser using Matched Field Processing[END_REF]. Note that, as expected from a surface array, the size of the focal spot increases with depth from ~2 m for tremor sources at a depth of 10 m to ~4 m for tremors at 20 m. Taking the half of the focal spot as the acoustic resolution, the source location uncertainty ranges then from ~1 m close to the surface to ~2 m for the deeper tremor sources.

The seismo-acoustic localization of the dominant tremor source during the cycle provides information from depths greater than the 14 m reached by the video observations [START_REF] Hutchinson | In-situ observations of Old Faithful Geyser[END_REF], and enables the reconstruction of both the 3D structure and the dynamics of the geyser during the eruption cycle.

The distribution of the dominant tremor sources along a SW-NE vertical plane highlights two distinct structures in the geyser subsurface where water boiling is concentrated (Figs. 1b,3). The first structure has a cylindrical geometry and corresponds to the extension of the geyser conduit down to 20 m (Fig. 1b, black dots). At a depth of 10 m, the X-Y elongation of the noise source locations is consistent with a fissure-like conduit (Fig. 3a) that was observed in situ [START_REF] Hutchinson | In-situ observations of Old Faithful Geyser[END_REF]. Note that the tremor sources in the vent between 12 m and 20 m dip towards the southwest, at an angle of ~24° (Fig. 3b). This trend points toward the vent orifice when extrapolated to the surface.

Another set of seismic sources is located at a depth range of 14 m to 28 m, which are only excited during the post-eruption recharge period. In the time-domain, these signals resemble the signals in the vent but their spectrum shows a slight shift to lower frequencies, with a central frequency decreasing by 0.5 Hz. These sources reveal a previously unknown lateral cavityshaped structure (Fig, 1b, red dots), which is located between Old Faithful Geyser and Split Cone, a weakly active geyser to the southwest (Fig. 3a). The flat roof of this structure, at ∼15 m below the surface (Fig. 3b) clearly stands out as the upper horizontal boundary of the seismic tremor location, and may correspond to a lithological boundary, e.g. the base of the sinter layer.

This deep structure, composed of a sub-vertical conduit connected at depth to a lateral sub-horizontal cavity is similar to first models of geyser with a siphon-like conduit connected to a bubble trap [Mackenzie, 1811] and to horizontal bubble-trap structures in Kamchatka geysers inferred from video images [START_REF] Belousov | Video observations inside conduits of erupting geysers in Kamchatka, Russia, and their geological framework: implications to the geyser mechanism[END_REF]. The temporal link between the geyser vent and the "recharge cavity" is evident from depth variations of the dominant tremor sources shown in Figure 2b. The location of seismic sources in the cavity area during the 15 min following the eruption (Fig. 1b, red dots) indicates the recharge of this deep structure. After this stage, the tremor sources slowly rise towards the geyser conduit, following an exponential time-dependence (Fig. 2b), up to a depth of ~10 m. The temporal evolution of the tremor source shows recurrent periods of water level stability at this depth, which has already been observed [START_REF] Kedar | The origin of harmonic tremor at Old Faithful Geyser[END_REF] at Old Faithful. It was explained by the enlargement of the conduit at a depth of 10 m, which forces the water level in the conduit to be nearly constant even though water is still filling the conduit (Fig. 4). This steady hydraulic state is generally followed by an unstable period in which some variability can occur both in the location of the source and in the amplitude of the seismic signal. This stage finally leads to an eruption, when the strong seismic signal observed at the surface is due to water splashing at the surface near the geophones.

In the data record, there is an initial rising phase that occurs 45 min before eruption. This ends with the sudden return of the seismic sources from the geyser conduit to the recharge zone, with no eruption and splashing at the surface. This could indicate an aborted eruption that would be in agreement with the bimodal eruption cycle observed by [START_REF] Kieffer | Seismicity at Old Faithful Geyser: an isolated source of geothermal noise and possible analogue of volcanic seismicity[END_REF].

Periods of rising dominant tremor sources in the geyser conduit can be closely matched with an exponential function (Fig. 2b). During the initial rising cycle, the characteristic time (12.5 +/-- 2 minutes) is twice as long as in the following two cycles (6.3 +/-- 2.3 minutes), and the critical depth of 10 m is not reached. This could explain why no eruption is observed during this cycle. The amount of water filling and its rate are governed by the hydraulic boundary conditions, and these depend on the volume of water discharged during the previous eruption. This behavior, whereby the triggering of an eruption depends on the amount and rate of water filling in the geyser (Steinberg, 1999), could be an explanation for the bimodal distribution of the eruption interval.

Discussion

During the stable periods with tremor sources at a depth of ~10 m, the sources have harmonic oscillations with amplitudes of about 1 m, and periods of 1 min (Fig. 4a). Similar periodic water temperature oscillations were observed [START_REF] Hutchinson | In-situ observations of Old Faithful Geyser[END_REF] using temperature probes lowered into the conduit, and explained to result from convective instability. We suggest that this periodic instability results from free oscillations [START_REF] Berest | Vibrations and free oscillations in salt caverns[END_REF] of a column of liquid water in the conduit above a compressible steam-rich fluid region in the recharge cavity, which has the role of a bubble trap.

Indeed, water in the upper part of the conduit (Fig. 3c) can be considered as a rigid mass that moves up and down as if above a spring, with this movement obeying the classical harmonic equation. The period of oscillation (T) depends on the geometry and compressibility of the spring-equivalent region, and this can be described by T = 2π/ω, where 2 S V h ω βρ = , and S is the cross-section of the conduit, h is the height of the liquid mass of density ρ, V is the volume of the cavity, and β is the compressibility of the fluid in the cavity. For slow pressure variations, the compressibility β can be written as 1/(γp 2 ), where γ ≈ 1.3 for steam, and p 2 is the hydrostatic pressure in the cavity.

A conduit diameter of ~0.1 m (as observed in videos at depths of 7 or 10 m), a gas-filled cavity at 15 m in depth, with volume ~62 m 3 and a 4 m characteristic size that fits with the seismic structure revealed in Fig. 3b would give the observed oscillation period of T = 60 s. It may be noted that pockets of air trapped in the horizontal enlargement of the upper conduit can act as additional barometric springs. The cavity volume inferred from this study is larger than the erupted volume at OFG, which was estimated from 20 m 3 [START_REF] Kieffer | Seismicity at Old Faithful Geyser: an isolated source of geothermal noise and possible analogue of volcanic seismicity[END_REF] to 40 m 3 [START_REF] Allen | Hot Springs of the Yellowstone National Park[END_REF], in agreement with the observation that several consecutive eruptions were required to clear OFG from an introduced tracer [START_REF] Fournier | Old Faithful: A Physical Model[END_REF].

As shown in Figure 4, the oscillation pattern is observed simultaneously in the depth of the seismic source, the seismic intensity and the central frequency of the seismic records. During an oscillation, water level rise (Fig. 4a) corresponds to fluid pressure increase, which causes a decrease in bubble diameter at the top of the water column. If a constant steam-injection rate during this stability period is assumed, then these smaller bubbles can induce both (1) a larger number of bubbles, to maintain a constant gas volume, and which implies a stronger seismic intensity (Fig. 4b), and (2) an increase in the central frequency of the seismic records (Fig. 4c).

In summary, these time-varying locations of seismic hydrothermal tremor source suggest that there is a compressible fluid-filled cavity at depth in the 'plumbing' system that can produce free oscillations of a water column that is excited by hydrodynamic pressure transients. This water level oscillation in the conduit is of fundamental importance for eruption triggering, as it favors water overflowing from the vent when the water level is close to the surface. Successive overflowing episodes constitute the preplay stage [START_REF] Kieffer | Seismicity at Old Faithful Geyser: an isolated source of geothermal noise and possible analogue of volcanic seismicity[END_REF], which induces boiling by reducing the hydrostatic pressure in the conduit that triggers the eruption.

Conclusion

In conclusion, data from dense seismic arrays surrounding very energetic seismic radiators such as geysers that are processed with precise localization techniques enable detailed mapping of structures in hydrothermal systems. This in turn allows monitoring of multiphase fluid flow in the system, which improves our knowledge of geyser dynamics. In recent models of geyser dynamics, the existence of a constriction in the vertical conduit [START_REF] Kieffer | Seismicity at Old Faithful Geyser: an isolated source of geothermal noise and possible analogue of volcanic seismicity[END_REF] or of highpermeability contrast in the fractured rock medium [START_REF] Ingebritsen | Controls on geyser periodicity[END_REF] was necessary to explain the generation of the flow intermittency in geysers. In this study, the presence of a subterranean lateral cavity revealed by the location of hydrothermal tremor sources during a cycle eruption appears to play a major role in the eruptive dynamics of Old Geyser. ., 179, 11-18. Mackenzie G. P. (1811), Travels in the Island of Iceland (T. Allan, Edinburgh). Steinberg, A.S. (1999), An experimental study of geyser eruption periodicity, Power engineering, 366, 47-50. Normalized root-mean-square amplitude averaged over the whole network. The gray rectangle with higher amplitudes corresponds to water splashing at the surface during the geyser eruption. The steady period during which water level oscillations are observed (see Fig. 4) is delimited by the two vertical dashed lines. (B) Depth-versus-time location of the tremor sources after beamforming, over three successive geyser cycles. The black dots correspond to tremor sources located in the geyser conduit, and the red dots to tremor sources in the recharge cavity (see Figs. 1 and3). The dotted blue lines correspond to the exponential fit of the water level rising periods in the geyser conduit. Fig. 3. Spatial representation of the tremor sources obtained from beamforming and the proposed structure of OFG reservoir. (A) X-Y projection of the tremor sources showing the sources inside the geyser conduit (black dots) and the sources in the recharge cavity (red dots). The insert corresponds to the X-Y elliptic shape of the cluster of tremor sources in the fissure-like conduit (from 12 m to 20 m; dashed line), to be compared to the X-Y elongation observed at the surface (full line). (B) Vertical projection of the noise sources along the OFG-SC segment and proposed model of structure, with a sub-vertical conduit connected to a lateral cavity. The flat roof of the recharge cavity is clearly highlighted by the density of red dots. The upper part of the conduit and the constriction are adapted from [START_REF] Hutchinson | In-situ observations of Old Faithful Geyser[END_REF], the widening at about 10 m depth is represented as horizontal fractures. (C) Schematic diagram of an oscillating liquid-water column over a compressible bubble-trap cavity, equivalent to a mass-spring harmonic oscillator.

Res

Figure captions

Fig. 1 .

 1 Fig. 1. Seismic array and 3D representation of the hydrothermal tremor sources at Old Faithful Geyser (OFG). (A) Shaded relief map of OFG area, with UTM coordinates (m). The red square indicates the study area. (B) Elevation model (with 1-m-interval in color) and position of the 96 vertical geophones (+), with reference to OFG vent. SC refers to Split Cone geyser. The location of the dominant seismic sources at depth discriminates between two structures: the geyser conduit and its vertical downward extension (black dots) and a lateral, cavity-shape structure (red dots).

Fig. 2 .

 2 Fig. 2. Time variation of the tremor amplitude and depth of the dominant tremor source. (A)Normalized root-mean-square amplitude averaged over the whole network. The gray rectangle with higher amplitudes corresponds to water splashing at the surface during the geyser eruption. The steady period during which water level oscillations are observed (see Fig.4) is delimited by the two vertical dashed lines. (B) Depth-versus-time location of the tremor sources after beamforming, over three successive geyser cycles. The black dots correspond to tremor sources located in the geyser conduit, and the red dots to tremor sources in the recharge cavity (see Figs.1 and 3). The dotted blue lines correspond to the exponential fit of the water level rising periods in the geyser conduit.

Fig. 4 .

 4 Fig. 4. Harmonic oscillations observed over a 5-min-long interval in a steady period (see Figure 2A). The blue stars correspond to data points. The black line is a band-passed filtered interpolation. (A) Depth-versus-time of the tremor sources. (B) Root-mean square amplitude recorded by the geophones around the vent. (C) Time evolution of the central frequency change of the seismic tremor.

  Figure2