P Gaudry

D Kohel

B Smith

Counting Points on Genus 2 Curves with Real Multiplication

We present an accelerated Schoof-type point-counting algorithm for curves of genus 2 equipped with an efficiently computable real multiplication endomorphism. Our new algorithm reduces the complexity of genus 2 point counting over a finite field Fq of large characteristic from O(log 8 q) to O(log 5 q). Using our algorithm we compute a 256-bit prime-order Jacobian, suitable for cryptographic applications, and also the order of a 1024-bit Jacobian.

Introduction

Cryptosystems based on curves of genus 2 offer per-bit security and efficiency comparable with elliptic curve cryptosystems. However, many of the computational problems related to creating secure instances of genus 2 cryptosystems are considerably more difficult than their elliptic curve analogues. Point counting-or, from a cryptographic point of view, computing the cardinality of a cryptographic group-offers a good example of this disparity, at least for curves defined over large prime fields. Indeed, while computing the order of a cryptographicsized elliptic curve with the Schoof-Elkies-Atkin algorithm is now routine, computing the order of a comparable genus 2 Jacobian requires a significant computational effort [START_REF] Gaudry | Construction of secure random curves of genus 2 over prime fields[END_REF][START_REF] Gaudry | Genus 2 point counting over prime fields[END_REF].

In this article we describe a number of improvements to the classical Schoof-Pila algorithm for genus 2 curves with explicit and efficient real multiplication (RM). For explicit RM curves over F p , we reduce the complexity of Schoof-Pila from O(log 8 p) to O(log 5 p). We applied a first implementation of our algorithms to find prime-order Jacobians over 128-bit fields (comparable to prime-order elliptic curves over 256-bit fields, and therefore suitable for contemporary cryptographic applications). Going further, we were able to compute the order of an RM Jacobian defined over a 512-bit prime field, far beyond the cryptographic range. (For comparison, the previous record computation in genus 2 was over a 128-bit field.)

While these RM curves are special, they are not "too special": Every ordinary genus 2 Jacobian over a finite field has RM; our special requirement is that this RM be known in advance and be efficiently computable. The moduli of curves with RM by a fixed ring form 2dimensional subvarieties (Humbert surfaces) in the 3-dimensional moduli space of all genus 2 curves. We can generate random curves with the specified RM by choosing random points on an explicit model of the corresponding Humbert surface [START_REF] Gruenewald | Computing Humbert surfaces and applications[END_REF]. In comparison with elliptic curves, for which the moduli space is one-dimensional, this still gives an additional degree of freedom in the random curve selection. To generate random curves with efficiently computable RM, we choose random curves from some known one and two-parameter families (see §4).

Curves with efficiently computable RM have an additional benefit in cryptography: the efficient endomorphism can be used to accelerate scalar multiplication on the Jacobian, yielding faster encryption and decryption [START_REF] Kohel | Efficiently computable endomorphisms for hyperelliptic curves[END_REF][START_REF] Park | Speeding up point multiplication on hyperelliptic curves with efficientlycomputable endomorphisms[END_REF][START_REF] Takashima | A new type of fast endomorphisms on Jacobians of hyperelliptic curves and their cryptographic application[END_REF]. The RM formulae are also compatible with fast arithmetic based on theta functions [START_REF] Gaudry | Fast genus 2 arithmetic based on Theta functions[END_REF].

Let C be a curve of genus 2 over a finite field F q , of odd characteristic, defined by an affine model y 2 = f (x), where f is a squarefree polynomial of degree 5 or 6 over F q . Let J C be the Jacobian of C; we assume J C is ordinary and absolutely simple. Points on J C correspond to degree-0 divisor classes on C; we use the Mumford representation for divisor classes together with the usual Cantor-style composition and reduction algorithms for divisor class arithmetic [START_REF] Galbraith | Efficient hyperelliptic arithmetic using balanced representation for divisors[END_REF][START_REF] Cantor | Computing in the Jacobian of a hyperelliptic curve[END_REF]. Multiplication by ℓ on J C is denoted by [ℓ], and its kernel by J C [ℓ]. More generally, if φ is an endomorphism of J C then J C [φ] = ker(φ), and if S is a set of endomorphisms then J C [S] denotes the intersection of ker(φ) for φ in S.

The Characteristic Polynomial of Frobenius

We let π denote the Frobenius endomorphism of J C , with Rosati dual π † (so ππ † = [q]). The characteristic polynomial of π has the form

χ(T ) = T 4 -s 1 T 3 + (s 2 + 2q) T 2 -qs 1 T + q 2 , (1)
where s 1 and s 2 are integers, and s 2 is a translation of the standard definition. The polynomial χ(T ) determines the cardinality of J C (F q k ) for all k: in particular, #J C (F q ) = χ(1). We refer to the determination of χ(T ) as the point counting problem.

The polynomial χ(T ) is a Weil polynomial: all of its complex roots lie on the circle |z| = √ q.

This implies the Weil bounds

|s 1 | ≤ 4 √ q and |s 2 | ≤ 4q. (2)
However, the possible values of (s 1 , s 2 ) do not fill the whole rectangle specified by the Weil bounds. Rück [START_REF] Rück | Abelian surfaces and jacobian varieties over finite fields[END_REF]Theorem 1.1] shows that in fact s 1 and s 2 satisfy s 2 1 -4s 2 ≥ 0 and s 2 + 4q ≥ 2|s 1 |, so the possible values of (s 1 , s 2 ) are in the following domain:

s 1 / √ q s 2 /q 2.

The Classical Schoof-Pila Algorithm for Genus Curves

The objective of point counting is to compute χ(T ), or equivalently the tuple of integers (s 1 , s 2 ). When the characteristic of F q is large, the conventional approach is to apply the Schoof-Pila algorithm as far as is practical, before passing to a baby-step giant-step algorithm if necessary (see §2.5). The strategy of Schoof's algorithm and its generalizations is to compute the polynomials χ ℓ (T ) = χ(T ) mod (ℓ) for sufficiently many primes (or prime powers) ℓ to reconstruct χ(T ) using the Chinese Remainder Theorem.

Since χ ℓ (T ) is the characteristic polynomial of π restricted to

J C [ℓ] (see [16, Proposition 2.1]), we have χ ℓ (π)(D) = 0 for all D in J C [ℓ].
Conversely, to compute χ ℓ (T ) we let D be a generic element of J C [ℓ] (as in §2.3 below), compute the three points (π 2 + [q]) 2 (D), (π 2 + [q])π(D), and π 2 (D), and then search for the coefficients (s 1 , s2 ) of χ ℓ (T ) in (Z/ℓZ) 2 , for which the linear relation

(π 2 + [q]) 2 (D) -[s 1 ] (π 2 + [q])π(D) + [s 2 ] π 2 (D) = 0 (3)
holds. If the minimal polynomial of π on J C [ℓ] is a proper divisor of χ ℓ (T )-which occurs for at most a finite number of ℓ dividing disc(χ)-then the polynomial so determined is not unique, but χ ℓ (T ) can be determined by deducing the correct multiplicities of its factors.

Once we have computed χ ℓ (T ) for sufficiently many ℓ, we reconstruct χ(T ) using the Chinese Remainder Theorem. The Weil and Rück bounds together with a weak version of the prime number theorem tell us how many ℓ are required: Pila notes in [16, §1] that the set of O(log q) primes ℓ < 21 log q will suffice. We analyse the complexity of the classical Schoof-Pila algorithm in §2.4.

Endomorphisms and Generic Kernel Elements

We now recall how to contruct an effective version of a generic ℓ-torsion element. We present it in a slightly more general setting, so that we can use this ingredient in the subsequent RM-specific algorithm. Therefore, we show how to compute with a generic element of the kernel of some endomorphism φ of J C , whereas φ is just [ℓ] in the classical algorithm. Definition 1. Fix an embedding P → D P of C in J C . We say that an endomorphism φ of J C is explicit if we can effectively compute polynomials d 0 , d 1 , d 2 , e 0 , e 1 , and e 2 such that if P = (x P , y P ) is a generic point of C, then the Mumford representation of φ(D P ) is given by

φ(D P ) = x 2 + d 1 (x P ) d 2 (x P ) x + d 0 (x P ) d 2 (x P )
, y -y P e 1 (x P ) e 2 (x P )

x + e 0 (x P ) e 2 (x P ) .

The d 0 , d 1 , d 2 , e 0 , e 1 , and e 2 are called the φ-division polynomials.

If φ is an explicit endomorphism, then we can use (4) (extending Z-linearly) to evaluate φ(D) for general divisor classes D in J C . In the case φ = [ℓ], the [ℓ]-division polynomials are the ℓ-division polynomials of Cantor [START_REF] Cantor | On the analogue of the division polynomials for hyperelliptic curves[END_REF]. The φ-division polynomials depend on the choice of embedding P → D P ; we will make this choice explicit when computing the φ-division polynomials for each of our families in §4.

To compute a generic element of J C [φ], we generalize the approach of [START_REF] Gaudry | Construction of secure random curves of genus 2 over prime fields[END_REF] (which computes generic elements of J C [ℓ]). The resulting algorithm is essentially the same as in [7, §3] (except for the parasite computation step, which we omit) with φ-division polynomials replacing ℓdivision polynomials, so we will only briefly sketch it here.

Let D = (x 2 + a 1 x + a 0 , y -(b 1 x + b 0 )) be (the Mumford representation of) a generic point of J C . We want to compute a triangular ideal I φ in F q [a 1 , a 0 , b 1 , b 0 ] vanishing on the nonzero elements of J C [φ]. The element D equals D (x 1 ,y 1 ) +D (x 2 ,y 2 ) , where (x 1 , y 1 ) and (x 2 , y 2 ) are generic points of C. To find a triangular system of relations on the a i and b i such that D is in J C [φ] we solve for x 1 , y 1 , x 2 , and y 2 in φ(D (x 1 ,y 1 ) ) = -φ(D (x 2 ,y 2 ) ), applying (4) and using resultants computed with the evaluation-interpolation technique of [7, §3.1]. We then resymmetrize as in [7, §3.2] to express the result in terms of the a i and b i . We can now compute with a "generic" element (x 2 + a 1 x + a 0 , y -(b 1 x + b 0 )) of J C [φ] by reducing the coefficients modulo I φ after each operation.

Following the complexity analysis of [7, §3.5], we can compute a triangular representation for I φ in O(δ 2 M(δ) log δ + M(δ 2 ) log δ) field operations, where δ is the maximum among the degrees of the φ-division polynomials, and M(d) is the number of operations required to multiply polynomials of degree d over F q . Using asymptotically fast multiplication algorithms, we can therefore compute I φ in O(δ 3 ) field operations. The degree of I φ is in O(δ 2 ); with this triangular representation, each multiplication modulo I φ costs O(δ 2 ) field operations.

Complexity of Classical Schoof-Pila Point Counting

Proposition 1. The complexity of the classical Schoof-Pila algorithm for a curve of genus 2 over F q is in O((log q) 8 ).

Proof. To determine χ(T ), we need to compute χ ℓ (T ) for O(log q) primes ℓ in O(log q). To compute χ ℓ (T ), we must first compute the ℓ-division polynomials, which have degrees in O(ℓ 2 ). We then compute the kernel ideal I ℓ ; according to the previous subsection, the total cost is in O(ℓ 6 ) field operations. The cost of checking (3) against a generic element of J C [ℓ] decomposes into the cost of computing Frobenius images of the generic element in O(ℓ 4 log q) and of finding the matching pair (s 1 , s2 ) in O(ℓ 5 ) field operations. So the total complexity for computing χ ℓ (T ) is in O(ℓ 4 (ℓ 2 + log q)) field operations. In terms of bit operations, for each ℓ bounded by O(log q), we compute χ ℓ (T ) in time O((log q) 7 ), and the result follows from the addition of these costs for all the different ℓ's. ⊓ ⊔

Baby-Step Giant-Step Algorithms

In practice, computing χ ℓ (T ) with classical Schoof-Pila becomes impractical for large values of ℓ. The usual approach is to carry out the Schoof-Pila algorithm to the extent possible, obtaining congruences for s 1 and s 2 modulo some integer M , before completing the calculation using a generic group algorithm such as baby-step giant-step (BSGS). Our BSGS algorithm of choice is the low-memory parallelized variant of the Matsuo-Chao-Tsuji algorithm [START_REF] Gaudry | A low-memory parallel version of Matsuo, Chao, and Tsuji's algorithm[END_REF][START_REF] Matsuo | An improved baby step giant step algorithm for point counting of hyperelliptic curves over finite fields[END_REF]. The bounds in (2) imply that the search space of candidates for (s 1 , s 2 ) is in O(q 3/2 ), and a pure BSGS approach finds (s 1 , s 2 ) in time and space O(q 3/4 ). However, when we apply BSGS after a partial Schoof-Pila computation, we obtain a congruence for (s 1 , s 2 ) modulo M . If M < 8q, then the size of the search space is reduced to O(q 3/2 /M 2 ), and the complexity for finding (s 1 , s 2 ) is reduced to O(q 3/4 /M ). For larger M , the value of s 1 is fully determined, and the problem is reduced to a one-dimensional search space of size O(q/M ) for which the complexity becomes O( q/M ).

Point Counting in Genus 2 with Real Multiplication

By assumption, J C is ordinary and simple, so χ(T ) is an irreducible polynomial defining a quartic CM-field with real quadratic subfield Q( √ ∆). We say that J C (and C) has real multiplication (RM) by Q( √ ∆ ). For a randomly selected curve, ∆ is in O(q); but in the sequel we consider families of curves with RM by Q( √ ∆) for small ∆ (= 5 or 8), admitting an explicit (in the sense of Definition 1) endomorphism φ such that

Z[φ] = Q( √ ∆) ∩ End(J C ) (5)
(that is, Z[φ] is the full real subring of End(J C )), and

disc (Z[φ]) = ∆.
We presume that the trace Tr(φ) and norm N(φ), such that φ 2 -Tr(φ)φ + N(φ) = 0, are known. We also suppose that φ is efficient, in the following sense:

Definition 2. We say that an explicit endomorphism φ is efficiently computable if the cost of evaluating φ at points of J C (F q ) requires only O(1) field operations (comparable to a few group operations in J C ). In practice, this means that the φ-division polynomials have small degree.

The existence of an efficiently computable φ and knowledge of ∆ allows us to make significant improvements to each stage of the Schoof-Pila algorithm. Briefly: in §3.2 we use φ to simplify the testing procedure for each ℓ; in §3.3 we show that when ℓ splits in Z[φ], we can use φ to obtain a radical reduction in complexity for computing χ ℓ (T ); and in §3. [START_REF] Cohen | Heuristics on class groups of number fields[END_REF] we show that knowing an effective φ allows us to use many fewer primes ℓ.

The RM Characteristic Polynomial

Let ψ = π + π † ; we consider Z[ψ], a subring of the real quadratic subring of End(J C ). The characteristic polynomial of ψ is the real Weil polynomial

ξ(T ) = T 2 -s 1 T + s 2 ; (6) the discriminant of Z[ψ] is ∆ 0 = s 2 1 -4s 2 . The analogue for (s 1 , ∆ 0 ) of Rück's bounds is (|s 1 | -4 √ q) 2 ≥ ∆ 0 = s 2 1 -4s 2 ≥ 0. (7)
Equation ( 5) implies that Z[ψ] is contained in Z[φ], so there exist integers m and n such that

ψ = m + nφ. (8)
Both s 1 and s 2 are determined by m and n: we have

s 1 = Tr(ψ) = 2m + nTr(φ) and s 2 = N(ψ) = (s 2 1 -n 2 ∆)/4. (9)
In fact n is the conductor of

Z[ψ] in Z[φ] up to sign: |n| = [Z[φ] : Z[ψ]],
and hence

∆ 0 = disc(Z[ψ]) = s 2 1 -4s 2 = n 2 ∆.
The square root of the bounds in [START_REF] Gaudry | Construction of secure random curves of genus 2 over prime fields[END_REF] gives bounds on s 1 and n:

4 √ q -|s 1 | ≥ ∆ 0 = |n| √ ∆ ≥ 0;
In particular, |s 1 | ≤ 4 √ q and |n| ≤ 4 q/∆. Applying the relation in [START_REF] Gaudry | Genus 2 point counting over prime fields[END_REF], we have the bounds

|m| ≤ 2(|Tr(φ)| + √ ∆) q/∆ and |n| ≤ 4 q/∆. (10)
Both |m| and |n| are in O( √ q).

An Efficiently Computable RM Relation

We can use our efficiently computable endomorphism φ to replace the relation of ( 3) with a more efficiently computable alternative. Multiplying [START_REF] Gaudry | A low-memory parallel version of Matsuo, Chao, and Tsuji's algorithm[END_REF] through by π, we have

ψπ = π 2 + [q] = mπ + nφπ.
We can therefore compute m = m mod ℓ and n = n mod ℓ by letting D be a generic ℓ-torsion point, computing the three points (π 2 + [q])(D), π(D), and φπ(D),

and then searching for the m and n in Z/ℓZ such that

(π 2 + [q])(D) -[ m]π(D) -[n]φπ(D) = 0 (11)
holds; we can find such an m and n in O(ℓ) group operations. Solving [START_REF] Kohel | Efficiently computable endomorphisms for hyperelliptic curves[END_REF] rather than (3) has several advantages. First, computing (π 2 + [q])(D), π(D), and φπ(D) requires only two applications of Frobenius, instead of the four required to compute (π 2 + [q]) 2 (D), (π 2 + [q])π(D), and π 2 (D) (and Frobenius applications are costly in practice). Moreover, either s 2 needs to be determined in O(q), or else the value of n in (3) leaves a sign ambiguity for each prime ℓ, because only n 2 mod ℓ can be deduced from (s 1 , s2 ). In contrast, [START_REF] Kohel | Efficiently computable endomorphisms for hyperelliptic curves[END_REF] determines n directly.

Exploiting Split Primes in

Q( √ ∆) Let Z[φ] ⊂ End(J C ) be an RM order in Q(φ) ∼ = Q( √ ∆).
Asymptotically, half of all primes ℓ split: (ℓ) = p 1 p 2 in Z[φ], where p 1 + p 2 = (1) (and this carries over to prime powers ℓ). This factorization gives a decomposition of the ℓ-torsion

J C [ℓ] = J C [p 1 ] ⊕ J C [p 2 ].
In particular, any ℓ-torsion point D can be uniquely expressed as a sum

D = D 1 + D 2 where D i is in J C [p i ].
According to the Cohen-Lenstra heuristics [START_REF] Cohen | Heuristics on class groups of number fields[END_REF], more than 75% of RM fields have class number 1; in each of the explicit RM families in §4, the order Z[φ] has class number 1. All ideals are principal in such an order, so we may find a generator for each of the ideals p i . Furthermore, the following lemma shows that we can find a generator which is not too large. Lemma 1. If p is a principal ideal of norm ℓ in a real quadratic order Z[φ], then there exists an effectively computable generator of p with coefficients in O( √ ℓ).

Proof. Let α be a generator of p, and ε a fundamental unit of

Z[φ]. Let γ → γ 1 and γ → γ 2 be the two embeddings of Z[φ] in R, indexed so that |α 1 | ≥ |α 2 | and |ε 1 | > 1 (replacing ε with ε -1 if necessary). Then R = log(|ε 1 |) is the regulator of Z[φ]. Set β = ε -k α, where k = [log(|α 1 / √ ℓ|)/R]; then β = a + bφ is a new generator for p such that - 1 2 ≤ log(|β i / √ ℓ|) R ≤ 1 2 •
From the preceding bounds, ⊓ ⊔

|β
We can now state the main theorem for RM point counting.

Theorem 1. There exists an algorithm for the point counting problem in a family of genus 2 curves with efficiently computable RM of class number 1, whose complexity is in O((log q) 5 ).

Proof. Let J C be a Jacobian in a family with efficiently computable RM by Z[φ]. Suppose that ℓ is prime, (ℓ) = p 1 p 2 in Z[φ], and that the p i are principal. By Lemma 2 we can compute representative α-division polynomials for

J C for each p in {p 1 , p 2 } in time O(ℓ), hence generic points D i in J C [p i ].
We recall that ( 11) is the homomorphic image under π of the equation

ψ(D) -[ m](D) -[n]φ(D) = 0.
When applied to D i in J C [p i ], both ψ and φ act as elements of Z[φ]/p i ∼ = Z/ℓZ. Moreover xi = φ mod p i is known, and it remains to determine ȳi = ψ mod p i by means of the discrete logarithm

ψ(D i ) = [ȳ i ](D i ) = [ m + nx i ](D i )
in the cyclic group D i ∼ = Z/ℓZ. The application of π transports this discrete logarithm problem to that of solving for ȳi in

D ′′ i = [ȳ i ]D ′ , where D ′ i = π(D i ) and D ′′ i = (π 2 + [q])(D i ).
By the CRT, from (ȳ 1 , ȳ2 ) in (Z/ℓZ) 2 we recover ȳ in Z[φ]/(ℓ), from which we solve for ( m, n) in (Z/ℓZ) 2 such that

ȳ = m + nφ ∈ Z[φ]/(ℓ).
The values of (s 1 , s2 ) are then recovered from [START_REF] Gaudry | Genus 2 point counting over prime fields[END_REF].

The ring Z[φ] is fixed, so as log q goes to infinity we find that 50% of all primes ℓ split in Z[φ] by the Chebotarev density theorem. It therefore suffices to consider split primes in O(log q). In comparison with the conventional algorithm presented in §2.2, we reduce from computation modulo the ideal for J C [ℓ] of degree in O(ℓ 4 ), to computation modulo the ideals for J C [p i ] of degree in O(ℓ 2 ). This means a reduction from O(ℓ 4 (ℓ 2 + log q)) to O(ℓ 2 (ℓ + log q)) field operations for the determination of each χ ℓ (T ), giving the stated reduction in total complexity from O((log q) 8 ) to O((log q) 5 ).

⊓ ⊔ Remark 1. Computing (m, n) instead of (s 1 , s 2 ) allows us to reduce the number of primes ℓ to be considered by about a half, since by [START_REF] Gruenewald | Computing Humbert surfaces and applications[END_REF] their product needs to be in O( √ q) instead of O(q). While this changes only the constant in the asymptotic complexity of the algorithm, it yields a significant improvement in practice.

Remark 2. If Z[φ]

does not have class number 1, and if (ℓ) = p 1 p 2 where the p i are not principal, then we may use a small complementary ideal (c) = c 1 c 2 such that c i p i are principal in order to apply Lemma 2 to a larger proportion of small ideals. Moreover, if ( m, n) is known modulo c, this can be used to reduce the discrete log problem modulo ℓ. Again, since a fixed positive density 1/2h of primes are both split and principal, where h is the class number of Z[φ], this does not affect the asymptotic complexity. Moreover, the first occurrence of a nontrivial class group is for ∆ = 65, beyond the current range for which an explicit RM construction is currently known.

Shrinking the BSGS Search Space

In the context of the conventional Schoof-Pila algorithm, we need to find s 1 in O( √ q) and s 2

in O(q). However, [START_REF] Gaudry | A low-memory parallel version of Matsuo, Chao, and Tsuji's algorithm[END_REF], and the effective form of (11) (valid for all points D of J C ), replaces the determination of (s 1 , s 2 ) with the tuple (m, n) of integers in O( √ q). As a result, the search space is reduced from O(q 3/2 ) to O(q). Thus the BSGS strategy can find (m, n) (which determines (s 1 , s 2 )) in time and space O( √ q), compared with O(q 3/4 ) when searching directly for (s 1 , s 2 ).

As in the general case, if one knows (m, n) modulo an integer M , then the area of the search rectangle is reduced by a factor of M 2 , so we find the tuple (m, n) in O( √ q/M ) group operations. Contrary to the general case of §2.5, since m and n have the same order of magnitude, the speed-up is always by a factor of M .

Examples of Families of Curves with Explicit RM

We now exhibit some families of curves and efficient RM endomorphisms that can be used as sources of inputs to our algorithm.

Correspondences and Endomorphisms

To give a concrete representation for endomorphisms of J C , we use correspondences: that is, divisors on the surface C × C. Suppose that R is a curve on C × C, and let π 1 : R → C and π 2 : R → C be the restrictions to R of the natural projections from C × C onto its first and second factors. We have a pullback homomorphism (π 1 ) * : Pic(C) → Pic(R), defined by (π 1 ) * P ∈C(Fq)

n P P = P ∈C(Fq) n P Q∈π -1 1 (P ) Q ,
where the preimages Q are counted with the appropriate multiplicities. (A standard moving lemma shows that we can always choose divisor class representatives so that each π -1 (P ) is zero-dimensional.) We also have a pushforward homomorphism (π 2 ) * : Pic(R) → Pic(C), defined by (π 2 ) * Q∈R(Fq)

n Q Q = Q∈R(Fq) n Q π 2 (Q) .
Note that (π 1 ) * maps Pic n (C) into Pic (n deg π 1 ) (R) and (π 2 ) * maps Pic n (R) into Pic n (C) for all n. Hence (π 2 ) * • (π 1 ) * maps Pic 0 (C) into Pic 0 (C), so we have an induced endomorphism

φ = (π 2 ) * • (π 1 ) * : J C → J C .
We write x 1 , y 1 and x 2 , y 2 for the coordinates on the first and second factors of C × C, respectively (so π i (x 1 , y 1 , x 2 , y 2 ) = (x i , y i )). In our examples, the correspondence R will be defined by two equations:

R = V (A(x 1 , x 2 ), B(x 1 , y 1 , x 2 , y 2 )) .
On the level of divisors, the image of a generic point P = (x P , y P ) of C (that is, a generic prime divisor) under the endomorphism φ is given by φ : (x P , y P ) -→ V (A(x P , x), B(x P , y P , x, y)) .

Using the relations y 2 P = f (x P ) and y 2 = f (x) (and the fact that correspondences cut out by principal ideals induce the zero homomorphism), we can easily replace A and B with Cantorreducible generators to derive the Mumford representation of φ(P ), and thus the φ-division polynomials.

A 1-dimensional Family with RM by

Z[ 1 + √ 5 /2]
Let t be a free parameter, and suppose that q is not a power of 5. Let C T be the family of curves of genus 2 over F q considered by Tautz, Top, and Verberkmoes in [START_REF] Tautz | Explicit hyperelliptic curves with real multiplication and permutation polynomials[END_REF]Example 3.5], defined by

C T : y 2 = x 5 -5x 3 + 5x + t. Let τ 5 = ζ 5 + ζ -1 5
, where ζ 5 is a 5th root of unity in F q . Let φ T be the endomorphism induced by the (constant) family of correspondences

R T = V x 2 1 + x 2 2 -τ 5 x 1 x 2 + τ 2 5 -4, y 1 -y 2 ⊂ C T × C T .
(Note that R T and φ T are defined over F q (τ 5 ), which is equal to F q if and only if q ≡ ±2 mod 5.) The family C T has an unique point P ∞ at infinity, which we can use to define an embedding P = (x P , y P ) -→ D P := [(P ) -(P ∞ )] ↔ (x -x P , y -y P ) of C T in J C T . With respect to this embedding, the φ T -division polynomials are

d 2 = 1, d 1 = -τ 5 x, d 0 = x 2 + τ 2 5 -4, e 2 = 1, e 1 = 0, e 0 = 1. Proposition 2. The minimal polynomial of φ T is T 2 +T -1: that is, φ T acts as multiplication by -(1 + √ 5)/2 on J C T . A prime ℓ splits into two principal ideals in Z[φ T ] if and only if ℓ ≡ ±1 mod 5.
Proof. The first claim is proven in [19, §3.5]. More directly, if P and Q are generic points of C T , then on the level of divisors we find

(φ 2 T + φ T )((P ) -(Q)) = (P ) -(Q) + div y -y(P ) y -y(Q) . Hence Z[φ T ] is isomorphic to the ring of integers of Q( √ 5
). The primes ℓ splitting in Q( √ 5) are precisely those congruent to ±1 modulo 5; and Q( √ 5) has class number 1, so the primes over ℓ are principal.

⊓ ⊔

The Igusa invariants of C T , viewed as a point in weighted projective space, are (140 : 550 : 20(32t 2 -3) : 25(896t 2 -3109) : 64(t 2 -4) 2 ); in particular, C T has a one-dimensional image in the moduli space of curves of genus 2. The Jacobian of the curve with the same defining equation over Q(t) is absolutely simple (cf. [START_REF] Kohel | Efficiently computable endomorphisms for hyperelliptic curves[END_REF]Remark 15]).

A 2-dimensional Family with RM by

Z[ 1 + √ 5 /2]
Let s and t be free parameters, and consider the family of curves C H : y 2 = F H (x), where

F H (x) = sx 5 -(2s + t)x 4 + (s 2 + 3s + 2t -1)x 3 -(3s + t -3)x 2 + (s -3)x + 1.
This family is essentially due to Humbert; it is equal to the family of Mestre [13, §2.1] with (U, T ) = (s, t), and the family of Wilson [20, Proposition 3.4.1] with (A, B) = (s, -t -3s + 3). The family has a full 2-dimensional image in the moduli space of genus 2 curves. Let R H be the family of correspondences on C H × C H defined by

R H = V x 2 1 x 2 2 + s(s -1)x 1 x 2 -s 2 (x 1 -x 2 ) + s 2 , y 1 -y 2 ;
let φ H be the induced endomorphism. The family C H has a unique point P ∞ at infinity, which we can use to define an embedding

P = (x P , y P ) -→ D P := [(P ) -(P ∞ )] ↔ (x -x P , y -y P ) of C H in J C H .
With respect to this embedding, the φ H -division polynomials are For an example with ∆ = 8, we present a twisted and reparametrized version of a construction due to Mestre [START_REF] Mestre | Couples de jacobiennes isogènes de courbes hyperelliptiques de genre arbitraire[END_REF]. Let s and t be free parameters, let v(s) and n(s) be the rational functions

d 2 = x 2 , d 1 = (s 2 -s)x + s 2 , d 0 = -s 2 x + s 2 , e 2 = 1, e 1 = 0, e 0 = 1. Proposition 3. The minimal polynomial of φ H is T 2 +T -1: that is, φ H acts as multipliction by -(1 + √ 5)/2 on J C H . A prime ℓ splits into two principal ideals in Z[φ H ] if and only if ℓ ≡ ±1 mod 5.
v = v(s) := s 2 + 2 s 2 -2 and n = n(s) := 4s(s 4 + 4) (s 2 -2) 3 ,
and let C M be the family of curves defined by

C M : y 2 = F M (x) := (vx -1)(x -v)(x 4 -tx 2 + vt -1).
The family of correspondences on C M × C M defined by

R M = V x 2 1 x 2 2 -v 2 (x 2 1 + x 2 2 ) + 1, y 1 y 2 -n(x 2 1 + x 2 2 -t)(x 1 x 2 -v(x 1 + x 2 ) + 1) induces an endomorphism φ M of J C M .
The family C M has two points at infinity, P + ∞ and P - ∞ . which are generically only defined over a quadratic extension of F q (s, t). Let D ∞ = (P + ∞ ) + (P - ∞ ) denote the divisor at infinity. We can use the rational Weierstrass point P v = (v, 0) on C M to define an embedding P = (x P , y P ) -→ D P := [(P )

+ (P v ) -D ∞ ] ↔ (x -x P )(x -v), y - y P x P -v (x -v)
of C M in J C M (the appropriate composition and reduction algorithms for divisor class arithmetic on genus 2 curves with an even-degree model, such as J C M , appear in [START_REF] Galbraith | Efficient hyperelliptic arithmetic using balanced representation for divisors[END_REF].) With respect to this embedding, the φ M -division polynomials are

d 2 = x 2 -v 2 , e 2 = (x 2 -v 2 )F M (x), d 1 = 0, e 1 = n(x -v)(x 4 -tx 2 + tv 2 -1), d 0 = -v 2 x 2 + 1, e 0 = n(vx -1)(x 4 -tx 2 + tv 2 -1).
Proposition 4. The minimal polynomial of φ M is T 2 -2: that is, φ M acts as multiplication by √ 2 on J C M . A prime ℓ splits into two principal ideals in Z[φ M ] if and only if ℓ ≡ ±1 mod 8.

Proof. Let P and Q be generic points of C M . An elementary but lengthy calculation shows that on the level of divisors,

φ 2 M ((P ) -(Q)) = 2(P ) -2(Q) + div x + x(P ) x + x(Q) , so φ 2 M ([D]) = 2[D] for all [D] in Pic 0 (C M ). Hence φ 2 M = [2],
2 -tT + v 2 t -1 in F q (v, t).
Mestre's curves C ′ and C are equal (over F q (v, a 1 , a 2 )) to our C M and its quadratic twist by A = 2(v 2 -1)(v 2 + 1) 2 = (2n) 2 , respectively. We may specialize the proofs in [START_REF] Mestre | Couples de jacobiennes isogènes de courbes hyperelliptiques de genre arbitraire[END_REF] to show that C M has a two-dimensional image in the moduli space of curves of genus 2, and that the Jacobian of the curve with the same defining equation over Q(s, t) is absolutely simple. Constructions of curves with RM by Z[ √ 2] are further investigated in Bending's thesis [START_REF] Bending | Curves of genus 2 with √ 2 multiplication[END_REF].

Remark 4. The algorithms described here should be readily adaptable to work with Kummer surfaces instead of Jacobians. In the notation of [START_REF] Gaudry | Fast genus 2 arithmetic based on Theta functions[END_REF], the Kummers with parameters (a, b, c, d)

satisfying b 2 = a 2 -c 2 -d 2 have RM by Z[ √ 2
], which can be made explicit as follows: the doubling algorithm decomposes into two identical steps, since (A : B : C : D) = (a : b : c : d), and the components after one step are the coordinates of a Kummer point. The step therefore defines an efficiently computable endomorphism which squares to give multiplication by 2.

Numerical Experiments

We implemented our algorithm in C++ using the NTL library. For non-critical steps, including computations in quadratic fields, we used Magma for simplicity. With this implementation, the determination of χ(T ) for a curve over a 128-bit prime field takes approximately 3 hours on one core of a Core2 processor at 2.83 GHz. This provides a proof of concept rather than an optimized implementation.

Cryptographic Curve Generation

When looking for a cryptographic curve, we used an early-abort strategy, where we switch to another curve as soon as either the order of the Jacobian order or its twist can not be prime.

Using our adapted version of Schoof algorithm, we guarantee that the group orders are not divisible by any prime that splits in the real field up to the CRT bound used.

In fact, any prime that divides the group order of a curve having RM by the maximal order of Q( √ ∆) must either be a split (or ramified) prime, or divide it with multiplicity 2. As a consequence, the early abort strategy works much better than in the classical Schoof algorithm, because, it suffices to test half the number of primes up to our CRT bound.

We ran a search for a secure curve over a prime field of 128 bits, using a CRT bound of 131. Our series of computations frequently aborted early, and resulted in 245 curves for which χ(T ) was fully determined, and for which neither the group order nor its twist was divisible by a prime less than 131. Considering these twists, this provided 490 group orders, of which 27 were prime, and therefore suitable for cryptographic use. We give here the data for one of these curves, that was furthermore twist-secure: both the Jacobian and the twist Jacobian order are prime.

Let C/F q , where q = 2 128 + 573, be the curve in the family C T of §4.2 specialized to t = 75146620714142230387068843744286456025. The characteristic polynomial χ(T ) is determined by (s 1 , s 2 ) = (-26279773936397091867, -90827064182152428161138708787412643439), giving prime group orders for the Jacobian: 115792089237316195432513528685912298808995809621534164533135283195301868637471, and for its twist: 115792089237316195414628441331463517678650820031857370801365706066289379517451.

We note that correctness of the orders is easily verified on random points in the Jacobians.

A Kilobit Jacobian

Let q be the prime 2 512 + 1273, and consider the curve over F q from the family C T of §4.2 specialized at t = 2908566633378727243799826112991980174977453300368095776223 2569868073752702720144714779198828456042697008202708167215 32434975921085316560590832659122351278.

This value of t was randomly chosen, and carries no special structure. We computed the values of the pair (s 1 mod ℓ, n mod ℓ) for this curve for each split prime ℓ up to 419; this was enough to uniquely determine the true value of (s 1 , n) using the Chinese Remainder Theorem. The numerical data for the curve follows: The order of the Jacobian is therefore N = (1 + q) 2 -s 1 (1 + q) + s 2 = 179769313486231590772930519078902473361797697894230657273 430081157732675805502375737059489561441845417204171807809 294449627634528012273648053238189262589020748518180898888 687577372373289203253158846463934629657544938945248034686 681123456817063106485440844869387396665859422186636442258 712684177900105119005520.

The total runtime for this computation was about 80 days on a single core of a Core 2 clocked at 2.83 GHz. In practice, we use the inherent parallelism of the algorithm, running one prime ℓ on each available core. We did not compute the characteristic polynomial modulo small prime powers (as in [START_REF] Gaudry | Genus 2 point counting over prime fields[END_REF]), nor did we use BSGS to deduce the result from partial modular information as in §3.4 (indeed, we were more interested in measuring the behaviour of our algorithm for large values of ℓ). These improvements with an exponential-complexity nature bring much less than in the classical point counting algorithms, since they have to be balanced with a polynomial-time algorithm with a lower degree. For this example, we estimate that BSGS and small prime powers could have saved a factor of about 2 in the total runtime.

Degrees of Division Polynomials

For each prime ℓ splitting in Z[φ T ], we report the degree of the α-division polynomial d 2 (where α is the endomorphism of norm ℓ that was used). By Lemma

1 + β 2 |Lemma 2 .

 122 = |2a + bTr(φ)| and |β 1 -β 2 | = |b √ ∆| are bounded by 2e R/2 √ ℓ. Since Tr(φ), ∆ and R are fixed constants, |a| and |b| are in O( √ ℓ). The "effective" part of the result follows from classical algorithms for quadratic fields.⊓ ⊔ Let J C be the Jacobian of a genus 2 curve over a finite field F q with an efficiently computable RM endomorphism φ. There exists an algorithm which, given a principal ideal p of norm ℓ in Z[φ], computes an explicit generator α of p and the α-division polynomials in O(ℓ) field operations. Proof. By Lemma 1, we can compute a generator α = [a] + [b]φ with a and b in O( √ ℓ). The [a]-and [b]-division polynomials have degrees in O(ℓ), and can be determined in O(ℓ) field operations. The division polynomials for the sum α = [a] + [b]φ require one sum and one application of φ; and since φ is efficiently computable, this increases the division polynomial degrees and computing time by at most a constant factor.

Proof.

  The first assertion is [13, Proposition 2] with n = 5; the rest of the proof is exactly as in Proposition 2. ⊓ ⊔ 4.4 A 2-dimensional Family with RM by Z[ √ 2]

  As noted above, this construction is a twisted reparametrization of a family of isogenies described by Mestre in [14, §2.1]. Let a 1 and a 2 be the roots of T

	and Z[φ M ] is isomorphic to the 2) are precisely those congruent to 2) has class number 1, so the primes over ℓ are principal. 2). The primes ℓ splitting in Q( √ ±1 modulo 8; further, Q( maximal order of Q( √ √ ⊓ ⊔
	Remark 3.

  1, deg(d 2 ) is in O(ℓ); the table below gives the ratio deg(d 2 )/ℓ, thus measuring the hidden constant in the O() notation. /ℓ 1.82 2.05 2.07 1.94 2.05 2.10 1.97 2.03 2.01 2.02 1.98 2.02 2.02 ℓ 139 149 151 179 181 191 199 211 229 239 241 251 269 deg d 2 /ℓ 2.12 2.04 1.99 2.00 2.01 2.09 2.21 1.99 2.18 2.01 2.05 2.07 2.17 ℓ 271 281 311 331 349 359 379 389 401 409 419 deg d 2 /ℓ 2.01 1.99 2.11 2.12 2.13 2.02 2.00 2.16 2.03 2.10 2.00 We have deg(d 1 ) = deg(d 2 ) + 1 and deg(d 0 ) = deg(d 2 ) + 2. All of these degrees depend only on the curve family C T , and not on the individual curve chosen.

	ℓ	11 19 29 31 41 59 61 71 79 89 101 109 131
	deg d 2


Petr Keil
email: pkeil@seznam.cz

Tomáš Herben
email: herben@site.cas.cz

James Rosindell
email: j.l.rosindell@leeds.ac.uk

David Storch
email: storch@cts.cuni.cz

Predictions of Taylor's power law, density dependence and pink noise from a neutrally modelled time series

Keywords: random walk, spectral analysis, mean-variance scaling, stochasticity, neutral theory of biodiversity, autocorrelation, 1/f noise

There has recently been increasing interest in neutral models of biodiversity and their ability to reproduce the patterns observed in nature, such as species abundance distributions. Here we investigate the ability of a neutral model to predict phenomena observed in single-population time series, a study complementary to most existing work that concentrates on snapshots in time of the whole community. We consider tests for density dependence, the dominant frequencies of population fluctuation (spectral density) and a relationship between the mean and variance of a fluctuating population (Taylor's power law). We simulated an archipelago model of a set of interconnected local communities with variable mortality rate, migration rate, speciation rate, size of local community and number of local communities. Our spectral analysis showed 'pink noise': a departure from a standard random walk dynamics in favor of the higher frequency fluctuations which is partly consistent with empirical data. We detected density dependence in local community time series but not in metacommunity time series. The slope of the Taylor's power law in the model was similar to the slopes observed in natural populations, but the fit to the power law was worse. Our observations of pink noise and density dependence can be attributed to the presence of an upper limit to community sizes and to the effect of migration which distorts temporal autocorrelation in local time series. We conclude that some of the phenomena observed in natural time series can emerge from neutral processes, as a result of random zero-sum birth, death and migration. This suggests the neutral model would be a parsimonious null model for future studies of time series data.

Introduction

The neutral theory of biodiversity (NTB; [START_REF] Hubbell | The unified neutral theory of biodiversity and biogeography[END_REF][START_REF] Chave | Neutral theory and community ecology[END_REF][START_REF] Alonso | abundance of animal and plant species[END_REF][START_REF] Leigh | Neutral theory: a historical perspective[END_REF][START_REF] Rosindell | Protracted speciation revitalizes the neutral theory of biodiversity[END_REF] has revolutionized the way we think about species coexistence over the last decade. The assumptions of the theory are controversial [START_REF] Tilman | Niche tradeoffs, neutrality, and community structure: A stochastic theory of resource competition, invasion, and community assembly[END_REF][START_REF] Leibold | Coexistence of the niche and neutral perspectives in community ecology[END_REF] and some argue that there is limited empirical support [START_REF] Mcgill | Empirical evaluation of neutral theory[END_REF], but the theory still has huge potential as a null model [START_REF] Gotelli | Null versus neutral models: what's the difference?[END_REF], serving as an essential first approximation to reality and a base for more complex theories. Even with its apparently oversimplified set of assumptions, the NTB is capable of predicting realistic speciesabundance distributions [START_REF] Hubbell | The unified neutral theory of biodiversity and biogeography[END_REF][START_REF] Etienne | A new sampling formula for neutral biodiversity[END_REF] and species-area curves (Rosindell andCornell, 2007, 2009), as well as many other ecological patterns [START_REF] Herben | Invasibility and species richness of a community: a neutral model and a survey of published data[END_REF][START_REF] Herben | Invasibility of neutral communities[END_REF]. Most of the current research is focused on the community-level predictions of a single snapshot in time at equilibrium (Chave et al., 2002;[START_REF] Mcgill | A test of the unified neutral theory of biodiversity[END_REF][START_REF] Turnbull | Niches, rather than neutrality, structure a grassland pioneer guild[END_REF][START_REF] Volkov | Density dependence explains tree species abundance and diversity in tropical forests[END_REF][START_REF] Dornelas | Coral reef diversity refutes the neutral theory of biodiversity[END_REF]. Dynamical behaviour of populations in NTB has rarely been studied. [START_REF] Azaele | Dynamical evolution of ecosystems[END_REF] used ratios of population sizes between two snapshots of time to estimate species turnover rates and extinction times of trees in Barro Colorado Island, Panama. [START_REF] Mutshinda | Species abundance dynamics under neutral assumptions: a Bayesian approach to the controversy[END_REF] produced the only study (as far as we are aware) that examines sufficiently long real-world time series from the perspective of NTB. They used the Bayesian statistical framework and a neutral model of community dynamics (zero-sum There are some characteristics of natural time-series that are often used as summary statistics and can be employed here in making comparisons with the time series produced by a neutral model. For example, natural time series have been shown

to obey some degree of stabilizing density dependence [START_REF] Woiwod | Patterns of density dependence in moths and aphids[END_REF][START_REF] Brook | Strength of evidence for density dependence in abundance time series of 1198 species[END_REF][START_REF] Sibly | On the stability of populations of mammals, birds, fish and insects[END_REF]. Natural populations also follow a specific relationship between mean population size (m) and variance of population size (s 2 ) called Taylor's Power Law (TPL) [START_REF] Taylor | Behavioural dynamics[END_REF][START_REF] Tokeshi | On the mathematical basis of the variance-mean power relationship[END_REF][START_REF] Ballantyne | The observed range for temporal mean-variance scaling exponents can be explained by reproductive correlation[END_REF] which has the form b m s v 2 . In natural populations, the exponent b of TPL is usually observed to satisfy 2 1 b [START_REF] Kendal | Taylor's ecological power law as a consequence of scale invariant exponential dispersion models[END_REF]. Natural time series can be considered as being the net result of a number of independent fluctuations each with a distinct frequency and amplitude. The amplitude (or 'spectral density') S(f) is a function of the frequency f and often satisfies [START_REF] Halley | Ecology, evolution and 1/f-noise[END_REF]. The exponent Ȗ usually falls between 0 and 2 for natural time series [START_REF] Arino | On the nature of population extremes[END_REF][START_REF] Halley | Ecology, evolution and 1/f-noise[END_REF]Inchausti andHalley, 2001, 2002), a property known as 'pink noise'.

J f f S / 1 ) ( v (
The null model for discrete population dynamics is an unbounded random walk of population size (N) and is described as

t t t e N N ln ln 1 (1)
where e t varies at random according to a given distribution with zero mean and is independent on population size [START_REF] Murdoch | Population regulation in theory and practice[END_REF][START_REF] Hanski | Random walks in a metapopulation: how much density dependence is necessary for long-term persistence?[END_REF]. The model is a discrete first-order Markov chain in which all the information that can influence the Although the term random walk is sometimes used to describe the population time series in neutral models [START_REF] Hubbell | The unified neutral theory of biodiversity and biogeography[END_REF][START_REF] Alonso | abundance of animal and plant species[END_REF], there are reasons why neutral models can behave differently from the random walk as defined above.

Firstly, there is usually a constraint on the total number of individuals (the 'zero-sum' rule) which imposes a clear upper limit to any population fluctuations [START_REF] Hubbell | The unified neutral theory of biodiversity and biogeography[END_REF][START_REF] Hubbell | The neutral theory of biodiversity an biogeography, and beyond[END_REF]. Secondly, immigration of individuals into the local community can cause a locally extinct species to reappear (zero is not necessarily absorbing; [START_REF] Hubbell | The unified neutral theory of biodiversity and biogeography[END_REF] and finally immigration can distort the temporal autocorrelation of N, (the number of individuals belonging to the species of interest in a local community). In the presence of immigration, N t+1 for a given local community does not depend solely on N t for that community (as in the autocorrelated random walk), it also depends on the N t for all the other local communities (the metacommunity).

In this paper we ask if a simple neutral model can produce time series that are similar to natural time series. In particular, we show that populations in the neutral model can produce realistic slopes of the Taylor's power law, have spectral properties similar to natural populations and seem to reveal density-dependency more frequently than expected in a random walk. Additionally, we explore the influence of the NTB parameters (size of local communities, number of local communities, mortality rate, migration and mutation rate) on these properties of population time series. Law. We used the simulation approach because it enabled us to record the species identity of every individual at each time step and across all local communities.

Furthermore, we wished to collect time series of limited lengths that are comparable to those of empirically collected time series in order to accurately mimic any sampling effects present in empirical data. This would be particularly difficult to approach analytically. Our simulations were conducted using a 'forwards in time' approach. A 'backwards in time' or coalescence approach can be much faster [START_REF] Rosindell | Coalescence models for spatial neutral ecology[END_REF], but the main strengths of this approach are only present when collecting a snap shot in time of a sample from a very large or infinite metacommunity. We are instead collecting single species time series data and sampling all individuals from a relatively small metacommunity: a scenario where coalescence is far less beneficial.

Four processes operate in the local community at each step: mortality (d; 0 < d < 1), local replacement (l; 0 < l < 1), immigration from outside of the local community (m; l m 1 ) and speciation (s; 1 0 s ) in that order (Fig. 1, Table 1).

To model mortality, d J u individuals were randomly removed from each local Etienne and Haegeman, 2010;[START_REF] Haegeman | Self-consistent approach for neutral community models with speciation[END_REF]. Pilot testing of a point mutation model [START_REF] Hubbell | The unified neutral theory of biodiversity and biogeography[END_REF] showed that the choice of a speciation model does not affect the outcome of the simulations, this is not entirely surprising because for a single species population time series, speciation can be regarded as another form of mortality: individuals belonging to our focal species that switch to a new species reduce the population of our focal species in the same way as dead individuals. Each simulation was launched with only one species that occurred in all local communities and was left running for 5000 steps which we found was always sufficient for the system to reach a steady state. The last 100 steps of the simulation were then taken and used for analyses.

A c c e p t e d m a n u s c r i p t

8

Spectral analysis

We performed 200 simulations for each combination of model parameters (Table 1).

From each simulation we randomly selected one species and from this species we used two time series: (1) time series of population sizes within a randomly selected local community and (2) time series of population sizes within the whole metacommunity. In cases where the species went extinct in one or both time series we randomly selected another one thereby conditioning the time series as being drawn from a locally extant species. For these time series we conducted spectral analysis [START_REF] Halley | Ecology, evolution and 1/f-noise[END_REF]pages 680-681 in Legendre and Legendre, 2003) and estimated the value of the exponent Ȗ in

J f f S / 1 ) ( v where S(f) is spectral density at given frequency f ( 2 / 0 n f f d d
; n is length of the time series). Values of S(f) were obtained through Fast Fourier Transform of the log-transformed time series data. The spectral exponent Ȗ was estimated as minus the regression slope of the )) ( log( f S versus ) log( f . This procedure is exactly the same as used by [START_REF] Inchausti | Investigating long-term ecological variability using the global population dynamics catabase[END_REF] and hence our estimates of Ȗ are directly comparable to theirs (Fig. 2 Database (kindly provided by Pablo Inchausti; see also Inchausti andHalley, 2001, 2002).

Testing for density dependence

Using the same data, we conducted three distinct tests for density dependence in the time series: Pollard's randomization test [START_REF] Pollard | The detection of density-dependence from a series of annual censuses[END_REF], Dennis and Taper's parametric bootstrap maximum likelihood test (PLBR; [START_REF] Dennis | Density dependence in time series observations of natural populations: Estimation and testing[END_REF]) and

Bulmer's test [START_REF] Bulmer | A statistical analysis of the 10-year cycle in Canada[END_REF]. The Pollards's test is a distribution-free likelihood ratio randomization test in which a random walk (Eq. 1) or a random walk with trend are density independent hypotheses, whereas the density dependent hypothesis is defined as a first-order linear autoregression model (stochastic Gompertz model):

t t t t e N r N N ln ln ln 1 E ) 1 ( z E (2)
where r and ȕ are model parameters. The Bulmer's test is a parametric predecessor of the Pollard's test and has been criticized for weak performance in cases where there was a trend in the data ( 0 z r

). The PLBR test is a likelihood ratio test which uses the first-order non-linear autoregression model (stochastic logistic model) as the density dependent hypothesis:

t t t t e N r N N E ln ln 1 (3)
All of these tests have been criticized for disputable performance when there is a census error affecting the estimation of population abundance in the field [START_REF] Shenk | Sampling-variance effects on detecting density dependence from temporal trends in natural populations[END_REF]Freckleton et al., 2006), but this is not the case for our simulated data. All three tests gave nearly identical results and we therefore only present the output from Pollard's randomization test [START_REF] Pollard | The detection of density-dependence from a series of annual censuses[END_REF] which is 'Pollard's P': the probability of observing the empirical data in the case of the density independent hypothesis (Eq. 1). Following the approach of [START_REF] Woiwod | Patterns of density dependence in moths and aphids[END_REF] we

Taylor's power law

To test if populations in our model obey the Taylor's Power Law (TPL) we performed 50 simulations for each combination of parameters (Table 1). For each species in the simulations we created a metapopulation matrix (MPM) of abundances in which rows represent the last 100 time steps of the simulation and columns represent the populations in different local communities. We excluded all MPMs in which the species had gone globally extinct. We studied the TPL in its temporal context [START_REF] Taylor | Comparative synoptic dynamics. I. Relationships between inter-and intra-specific spatial and temporal variance/mean population parameters[END_REF] which means that variances and means are calculated over each column (one column is one local population) in the MPM -when plotted, each point in the log-log plot therefore represents one local population. Using rows instead of columns in the MPM would yield the TPL in a spatial context, but the temporal context is more germane to this time series study. A linear regression was fitted in each of the log-log plots and the TPL slope (parameter b; Taylor, 1961;[START_REF] Tokeshi | On the mathematical basis of the variance-mean power relationship[END_REF] and the coefficients of determination (R 2 ) were calculated. As with our other results, the simulated values of TPL slopes and R 2 s were compared with empirically observed values (taken from Taylor andWoiwod, 1980, 1982). 2, Fig. 3) with higher mortality and higher migration each producing time series that were closer to pink noise ( 2J ) (Fig. 2).

Results

Spectral analysis

On the metacommunity scale, the median of Ȗ was 1.75 (2.5% and 97.5% quantiles were 1.05 and 2.07 respectively). We did not detect any influence of the model parameters on the value of Ȗ at the metacommunity scale (Table 2, Fig. 3). Our simulations produced Ȗ values falling within a narrower interval than is generally observed in nature (Fig. 2), however, they departed from 2 (the value for an unbounded random walk) in the same direction as for natural populations.

In general, the population time series emerging from our neutral model had lower Ȗ than would be expected for an autocorrelated random walk ( 2

J

). This could partly be an artifact caused by the insufficient length of time series. It was recently shown that random walks of length ~100 actually give estimates of between 1.5 and 1.8 instead of 2 because of the missing low-frequency fluctuations [START_REF] Miramontes | Estimating 1/f Į scaling exponents from short timeseries[END_REF]. Nevertheless, since empirically collected time series are rarely longer than 100 sequential observations, our simulations of time series of similar lengths are amply justified to purposefully include the effects of sampling over a limited time scale. Since we also showed that the vales of Ȗ tended to be lower at local scales when comparing time series of equal lengths, we can still conclude that some whitening of spectra at local community sales is caused by more than just the length of the sampled time series.

Testing for density dependence

A c c e p t e d m a n u s c r i p t 12

On the scale of local communities the median of Pollard's P was 0.14 (2.5% and 97.5% quantiles were 0 and 0.82 respectively) indicating a distribution skewed towards the lower values of Pollard's P. We also found that the values of Pollard's P were dependent on mortality rate and migration rate (Table 2, Fig. 4). Similarly to spectral analysis, simulations with higher mortality rate and higher migration rate produced time series with lower values of Pollard's P (Fig. 4). On the metacommunity scale the median value of Pollard's P was 0.5 (2.5% and 97.5% quantiles being 0.01 and 0.97 respectively). This indicates no clear signal of density dependence on the metacommunity scale. We did not detect any influence of the other model parameters on the values of logit-transformed Pollard's P (Table 2, Fig. 4).

Pollard's P varied with relative abundance of species in local communities but not in metacommunities (Fig. 6).

Taylor's power law (TPL)

The median slope of TPL was 1.73 (2.5% and 97.5% quantiles were 0.059 and 2.49 respectively). The coefficients of determination (R 2 ) of power law regressions varied, following skewed, non-normal distributions (Figs. 2, 5), having generally lower values to what has been observed in nature (Fig. 2). The effects of model parameters on the slope of TPL was not as distinct as in spectral analysis and density dependence analysis; we detected a weak influence of number of local communities and mortality rate on the slopes of TPL (Table 2, Fig. 5). Both the slopes and R 2 s showed negative dependence on the mean relative abundance of species (Fig. 7).

Discussion

A c c e p t e d m a n u s c r i p t 13

Spectral analysis and density dependence

The observed whitening of spectra in our model is not what we would expect under random walk dynamics. One possible reason for this is the effect of migration which can cause the spectral whitening. A population in an isolated local community (low or no immigration from outside) follows dynamics similar to an autocorrelated random walk ( 2 | J ) in which there is little or no dependence on the abundances in other local communities. However, after introducing immigration from outside of the local community, the autocorrelative nature of the process is broken and the species abundance in the whole metacommunity becomes important. The influence of single local community dynamics on the metacommunity are likely to be small (inversely proportional to the number of local communities), but in contrast the influence of the metacommunity on local community dynamics can be substantial (proportional to the migration rate). Hence, the migration causes the local time series to be whitened (nonautocorrelated in time) while the metacommunity time series remains "brownian" (autocorrelated in time). The mortality rate can also cause spectral whitening by adding random population fluctuations at high frequencies, but not at low ones subsequently lowering Ȗ. We expect that this effect of mortality rate would not be apparent in very long time series (> 400 time steps; [START_REF] Miramontes | Estimating 1/f Į scaling exponents from short timeseries[END_REF].

However, such long time series are nearly never available in ecology and we wanted to keep our time series comparably long to the real-world ones.

The spectral analysis results are mirrored by the results of the tests for density dependence. Density dependence is modeled either as stochastic Gompertz (Eq. 2) or stochastic logistic (Eq. 3) models in the tests we used (Shenk et al., 1998, see also Methods).

A c c e p t e d m a n u s c r i p t

14 It has been recognized for a long time [START_REF] Eberhardt | Correlation, regression and density dependence[END_REF]Freckleton et al., 2006) that a sequence of uncorrelated random numbers ("white noise") tends to follow a pattern consistent with density dependence. If you take white noise time series and plot changes of ln population size (

t t N N ln ln 1
) against ln N t (where N t is the number of individuals at time t) you will get a linear negative relationship, corresponding (spuriously) to the density dependent Gompertz model (Eq. 2) with 1 | E . Moreover, white noise can be very similar to the output of the stochastic logistic model (Eq. 3) with 0

E

. We have learned that the whitening of a spectra means bringing it from 2 | J down towards 0, which is white noise. We also know that migration (by breaking the autocorrelation structure) and mortality rate (by increasing spectral density at high frequencies in time series of limited length) cause spectral whitening. The detection of density dependence can also be attributed to the same mechanisms, although most ecologists would consider such density dependence to be spurious (Freckleton et al., 2006).

Ucorrelated noise in the form of census error also causes spurious detection of density dependence has long been known to population ecologists [START_REF] Shenk | Sampling-variance effects on detecting density dependence from temporal trends in natural populations[END_REF]Freckleton et al., 2006). Moreover, Freckleton et al. (2006) showed an example of a population that follows a random walk but consists of two subpopulations each of which apparently shows density dependence because there is random migration between them. That immigration can cause time series to resemble a sequence of uncorrelated random numbers was also noted by Anderson et al. (1982) and this is indeed the case for our neutral model.

There should also be real density dependence in our model. In an isolated community, the closer N t gets to zero (extinction) or to J (monodominance) the lower the probability of N t+1 differing from N t will be. When a species almost fills the local A c c e p t e d m a n u s c r i p t 15 community, its probability of increasing its abundance in the next time step is small, because it is unlikely that one of the few individuals belonging to a different species will be selected for mortality. Similarly, when the abundance of a species is low, its probability of reproducing is low because after mortality it is not likely to be selected to reproduce. The highest rate of fluctuations and lowest incidence of density dependence is then expected when a population has a relative abundance of 0.5 [START_REF] Hubbell | The unified neutral theory of biodiversity and biogeography[END_REF]. We detected this form of density dependence only in local communities at low relative abundances (Fig. 6). Species with high relative abundances were rare in our simulations and hence we were unable to clearly show the real density dependence for them (Fig. 6) even though in theory it must exist.

Taylor's power law (TPL)

The archipelago neutral model produced slopes of Taylor's Power Law with median close to the range observed in natural populations (1<b<2), but with higher incidence of low coefficients of determination. In order to explain this we need to summarize the mean-variance scaling properties of simple stochastic models. If a species undergoes an unrestricted random walk (Eq. 1), the TPL has a slope 1 b at very low population densities and gradually changes to 2 b at a certain density (Anderson et al., 1982;[START_REF] Tokeshi | On the mathematical basis of the variance-mean power relationship[END_REF][START_REF] Keeling | Simple stochastic models and their power-law type behaviour[END_REF]. The population density at which b starts to change from 1 to 2 depends on the relative magnitude of stochastic changes because of the high incidence of extinctions at low population densities which lower the variance. This effect can explain 2 b at very low densities. resembling the power-law with slope of 2. Although [START_REF] Hubbell | The unified neutral theory of biodiversity and biogeography[END_REF] reasoning focuses on the relationship between mean and variance in a set of species within one community (one point in the mean-variance plot represents one species) it can apply to the set of single-species local populations in our model: the mean relative abundances of species in local communities were rarely higher than 0.2 in our simulations and indeed, there is a much better fit of the TPL for species with low abundances (Fig. 7). The poor fit at higher relative abundances is then a result of fitting the linear TPL on the actually parabolic mean-variance relationship.

Consequently, as the mean relative abundance of species increases, the slope of the TPL decreases (Fig. 7).

Comparison with empirical data and conclusions

Our model produces qualitatively similar patterns to those observed in nature. This is perhaps surprising given the assumptions made by the neutral model, but it does again highlight its use as a powerful null model.

There is a general agreement that natural time series are not "brownian" (Fig. 1; [START_REF] Arino | On the nature of population extremes[END_REF]Inchausti andHalley, 2001, 2002) -they clearly fall into the whitened range of scaling exponents ( 2 0 J ; Fig. 2). Our simulations indicate that neutral models have the potential to explain part of this observed whitening, although this potential may be limited. Although higher migration rates can result in [START_REF] Volkov | Neutral theory and relative species abundance in ecology[END_REF][START_REF] Mutshinda | Species abundance dynamics under neutral assumptions: a Bayesian approach to the controversy[END_REF]. [START_REF] Mutshinda | Species abundance dynamics under neutral assumptions: a Bayesian approach to the controversy[END_REF] already showed that sizes of local communities fluctuate more than is expected from neutral theory. This resonates with our simulations of individual species time series which show rather high values of Ȗ (i.e. lower importance of high-frequency fluctuations)

when compared to empirical data. Recent progress in population ecology suggests that a promising explanation for these high fluctuations of natural populations might be environmental stochasticity [START_REF] Lande | Stochastic population dynamics in ecology and conservation[END_REF]Mutshinda et al., 2009) which is absent in neutral models, but not impossible to include in future and more general versions of the original theory [START_REF] Allen | Setting the absolute tempo of biodiversity dynamics[END_REF].

Realistic slopes of TPL can be produced by a number of different mechanisms [START_REF] Tokeshi | On the mathematical basis of the variance-mean power relationship[END_REF][START_REF] Keeling | Simple stochastic models and their power-law type behaviour[END_REF][START_REF] Kendal | Taylor's ecological power law as a consequence of scale invariant exponential dispersion models[END_REF][START_REF] Ballantyne | The observed range for temporal mean-variance scaling exponents can be explained by reproductive correlation[END_REF] and it is therefore difficult to infer any definitive underlying mechanisms from its shape [START_REF] Mitzenmacher | A brief history of generative models for power laws and lognormal distributions[END_REF]. Here we add neutral models to the list of possible mechanisms that can produce realistic TPL slopes, although not necessarily realistic fits.

In conclusion, our results show that realistic slopes of Taylor's power law and The neutral model encompasses these factors without being overly complex and we therefore suggest that neutral models make more powerful null models for time series data compared to alternatives such as Brownian random walks. (1980,1982). Empirical values of scaling exponent Ȗ of power spectra of 544 time series were taken from the Global Population Dynamics Database (GPDD, kindly provided by Pablo Inchausti; see also Inchausti andHalley, 2001, 2002). Box and whisker plots show medians, quartiles and outliers. Pollard's P is the probability of obtaining the data in case of density independence (Eq. 1; [START_REF] Pollard | The detection of density-dependence from a series of annual censuses[END_REF].

  to show that sizes of natural communities fluctuate more than expected under neutrality. However,[START_REF] Mutshinda | Species abundance dynamics under neutral assumptions: a Bayesian approach to the controversy[END_REF] focused only on the dynamics of local community size and they did not examine time series of individual species.

  t+1) is fully captured in the present state (t). Random walk time series are often termed "brownian" or autocorrelated noise. Random walk is synonymous for a completely density-independent process in ecology, it produces a slope of 2

  Instead of the classic neutral model used by[START_REF] Hubbell | The unified neutral theory of biodiversity and biogeography[END_REF] and[START_REF] Volkov | Neutral theory and relative species abundance in ecology[END_REF], we simulate a metacommunity consisting of an archipelago of L isolated local communities interconnected by migration at rate m (see Fig.1below[START_REF] Volkov | Patterns of relative species abundance in rainforests and coral reefs[END_REF]. Each local community contains J individuals all of which behave in an identical manner irrespective of the species to which they belong. The archipelago model enabled us to obtain data on abundances and their variances for one species at multiple local communities, which is necessary for the analysis of Taylor's Power

  time step with all individuals having an equal probability of removal. These removed individuals were then replaced by offspring from others in the local community and other local communities (immigration). Local replacement was modeled by drawing individuals to reproduce at random from the replaced by immigration from all of the other local communities. The probability of immigration of a species is given by its relative abundance in the whole metacommunity (excluding the local community in focus). A local population can occasionally give rise to a new species (a speciation event), this happens with a fixed probability at each time step for each local population of a species. When speciation occurs, individuals of the given species in the local community were randomly split into two parts, one of which forms the initial abundance of a new species (a version of random-fission model; Hubbell, 2001;

  noise'). We use the term 'whitening' of the spectra for a deviation of Ȗ from 2 towards 0.To assess effects of the neutral model parameters on Ȗ we fitted a multiple ANOVA model with Ȗ as the response variable and five explanatory variables: mortality rate, migration rate, speciation rate, number of local communities and number of individuals. Since we were dealing with simulated data we only report R 2 s and do not report any P-values because these could be modified arbitrarily by changing number of simulations. Comparisons were made with the values of Ȗ measured in natural populations taken from the Global Population Dynamics

  of Pollard's P by logit transformation. To assess effects of the neutral model parameters on Pollard's P we again fitted a multiple ANOVA model in which we used the values of logit-transformed Pollard's P as the response variable.

  Figure2. For local communities, the median value of Ȗ was 1.59 (2.5% and 97.5% quantiles were 1.26 and 2.19 respectively). The value of Ȗ was only influenced by mortality rate and migration rate (Table 2, Fig.3) with higher mortality and higher

Hubbell ( 2001 )

 2001 shows what happens with the mean-variance scaling if the stochastic population fluctuations are restricted by an upper boundary (the zero-sum rule) and are fed by immigration, although Hubbell does not use the term Taylor's Hubbell (2001) predicts hump-shaped parabolic relationship between abundance of species in a metacommunity and the variance of a local population size undergoing the zero-sum ecological drift. He also hypothesizes that most of the species in nature have low relative abundances and hence should occur in the left part (relative abundance < 0.2) of the parabola which is nearly linear in the log-log space,

  more realistic values of Ȗ, our migration rates are similar (10%) or much higher (30%) than what is usually observed in natural metacommunities (see

  density dependence can emerge from a simple neutral model, whilst spectral colours are less likely, but still possible to be reproduced by the neutral model. Still some other natural phenomena remain unexplained by our model, such as the possibility for very small Ȗ closer to white noise. Our findings suggest the interesting possibility that the properties of population time series observed in nature are in part due to simple factors such as sampling effects and dispersal limitation, which would then represent the most parsimonious explanation of observed patterns. The patterns of population fluctuations observed in nature may not reflect resource-driven species-specific processes at the metacommunity scale, with an important effect of stochastic mortality and random immigration of individuals into local communities.

Figure 1

 1 Figure legends

Figure 2

 2 Figure 2 Some aspects of real-world populations and those produced by our neutral model. Empirical slopes of Taylor's Power Law and coefficients of determination of TPL of 429 species of aphids, moths and birds were taken from Taylor & Woiwod

Figure 3 .Figure 4

 34 Figure 3 Results of the spectral analysis of population time series produced by a

Figure 5

 5 Figure 5 The slopes of temporal Taylor's power law (TPL) as produced by neutral models with various levels of model parameters. Solid horizontal lines delimit the range observed in natural populations. Shown as a box and whisker plot with median and quartiles shown, outliers are not displayed.

Figure 6

 6 Figure 6 Dependence of Pollard's P on relative abundance of species in local communities and in metacommunities. Shades of grey represent density of points in hexagonal bins.

Figure 7

 7 Figure 7 Dependence of slope and coefficient of determination (R 2 ) of Taylor's Power Law on mean relative abundance of species in local communities. Shades of grey represent density of points in hexagonal bins.

Table 2

 2 The strength of the effect of each parameter on the value of γ and the logittransformed Pollard's P. Parameter γ describes the 'colour' of the spectra of a time series. Pollard's P is the probability of obtaining the data in case of density independence (random walk). The values of R 2 were estimated using a multiple ANOVA. Figures in bold show values higher than 0.001.

		Local communities		Metacommunities
		R 2 (γ)	R 2 (logit of	R 2 (γ)	R 2 (logit of	R 2 (slope of
			Pollard's P)		Pollard's P)	TPL)
	Mortality rate	0.045	0.078	<0.001	<0.001	0.019
	Migration rate	0.11	0.19	<0.001	<0.001	0.0014
	Speciation rate	<0.001	<0.001	<0.001	<0.001	<0.001
	# of local	<0.001	0.014	<0.001	<0.001	0.048
	communities
	Size of local	<0.001	<0.001	<0.001	<0.001	0.0015
	community

Acknowledgements

We are grateful to Rampal Etienne, Stephen Cornell, William E. Kunin and Samuel Zajicek for discussions on the subject and to one anonymous referee for helpful comments. Pablo Inchausti generously provided us the GPPD data on spectral exponents. The study was supported by the EU FP7 project SCALES (project 226852), by the Grant Agency of Charles University (GAUK 106108), by the Grant Agency of the Academy of Sciences of the CR (IAA601970801) and by the Czech Ministry of Education (LC06073, MSM0021620845). James Rosindell was supported with funding from the EPSRC (grant EP/F043112/1). Petr Keil was supported with Marie Curie Early Stage Training fellowship at University of Leeds (UK).