C Kalla

C Klein

ON THE NUMERICAL EVALUATION OF ALGEBRO-GEOMETRIC SOLUTIONS TO INTEGRABLE EQUATIONS

Physically meaningful periodic solutions to certain integrable partial differential equations are given in terms of multi-dimensional theta functions associated to real Riemann surfaces. Typical analytical problems in the numerical evaluation of these solutions are studied. In the case of hyperelliptic surfaces ecient algorithms exist even for almost degenerate surfaces. This allows the numerical study of solitonic limits. For general real Riemann surfaces, the choice of a homology basis adapted to the antiholomorphic involution is important for a convenient formulation of the solutions and smoothness conditions. Since existing algorithms for algebraic curves produce a homology basis not related to automorphisms of the curve, we study symplectic transformations to an adapted basis and give explicit formulae for M-curves. As examples we discuss solutions of the Davey-Stewartson and the multi-component nonlinear Schrödinger equations.

Introduction

The importance of Riemann surfaces for the construction of almost periodic solutions to various integrable partial dierential equations (PDEs) was realized at the beginning of the 1970s by Novikov, Dubrovin and Its, Matveev. The latter found the Its-Matveev formula for the Korteweg-de Vries (KdV) equation in terms of multi-dimensional theta functions on hyperelliptic Riemann surfaces. Similar formulae were later obtained for other integrable PDEs as nonlinear Schrödinger (NLS) and sine-Gordon equations. For the history of the topic the reader is referred to the reviews [START_REF] Belokolos | Algebro-geometric approach to nonlinear integrable equations[END_REF] and [START_REF] Dubrovin | Theta functions and non-linear equations[END_REF]. Krichever [START_REF] Krichever | The averaging method for two-dimensional integrable equations[END_REF] showed that thetafunctional solutions to the Kadomtsev-Petviashvili equation can be obtained on arbitrary Riemann surfaces. The problems of real-valuedness and smoothness of these solutions were solved by Dubrovin and Natanzon in [START_REF] Dubrovin | Real theta-function solutions of the Kadomtsev-Petviashvili equation[END_REF].

Novikov criticized the practical relevance of theta functions since no numerical algorithms existed at the time to actually compute the found solutions. He suggested an eective treatment of theta functions (see, for instance, [START_REF] Dubrovin | Theta functions and non-linear equations[END_REF]) by a suitable parametrization of the characteristic quantities of a Riemann surface, i.e., the periods of holomorphic and certain meromorphic dierentials on the given surface. This program is limited to genera smaller than 4 since so-called Schottky relations exist for higher genus between the components of the period matrix of a Riemann surface. The task to nd such relations is known as the Schottky problem. This led to the famous Novikov conjecture for the Schottky problem that a Riemann matrix (a symmetric matrix with negative denite real part) is the matrix of B-periods of the normalized holomorphic dierentials of a Riemann surface Date : June 18, 2011. We thank D. Korotkin and V. Shramchenko for useful discussions and hints. This work has been supported in part by the project FroM-PDE funded by the European Research Council through the Advanced Investigator Grant Scheme, the Conseil Régional de Bourgogne via a FABER grant and the ANR via the program ANR-09-BLAN-0117-01.

if and only if Krichever's formula with this matrix yields a solution to the KP equation. The conjecture was nally proven by Shiota [START_REF] Shiota | Characterization of Jacobian varieties in terms of soliton equations[END_REF].

First plots of KP solutions appeared in [START_REF] Mumford | Tata Lectures on Theta. I and II[END_REF] and via Schottky uniformizations in [START_REF] Bordag | Periodic Multiphase Solutions of the Kadomtsev-Petviashvili-equation[END_REF]. Since all compact Riemann surfaces can be dened via non-singular plane algebraic curves of the form (1.1)

F (x, y) := N n=1 M m=1
a mn x m y n = 0, x, y ∈ C, with constant complex coecients a nm , Deconinck and van Hoeij developed an approach to the symbolic-numerical treatment of algebraic curves. This approach is distributed as the algcurves package with Maple, see [START_REF] Deconinck | Computing Riemann matrices of algebraic curves[END_REF][START_REF] Deconinck | Computing Riemann Theta Functions[END_REF][START_REF] Deconinck | Computing with plane algebraic curves and Riemann surfaces: the algorithms of the Maple package algcurves[END_REF]. A purely numerical approach to real hyperelliptic Riemann surfaces was given in [START_REF] Frauendiener | Hyperelliptic theta functions and spectral methods[END_REF][START_REF] Frauendiener | Hyperelliptic theta functions and spectral methods: KdV and KP solutions[END_REF], and for general Riemann surfaces in [START_REF] Frauendiener | Algebraic curves and Riemann surfaces in Matlab[END_REF]. For a review on computational approaches to Riemann surfaces the reader is referred to [START_REF]Computational Approach to Riemann Surfaces[END_REF].

In this paper we want to address typical analytical problems appearing in the numerical study of theta-functional solutions to integrable PDEs, and present the state of the art of the eld by considering concrete examples. The case of hyperelliptic Riemann surfaces (N = 2 in (1.1)) is the most accessible, since equation (1.1) can be solved explicitly for y, and since a basis for dierentials and homology can be given a priori. Families of hyperelliptic curves can be conveniently parametrized by their branch points. The codes [START_REF] Frauendiener | Hyperelliptic theta functions and spectral methods[END_REF][START_REF] Frauendiener | Hyperelliptic theta functions and spectral methods: KdV and KP solutions[END_REF] are able to treat eectively numerically collisions of branch points, a limit in which certain periods of the corresponding hyperelliptic surface diverge. If the limiting Riemann surface has genus 0, the theta series breaks down to a nite sum which gives for an appropriate choice of the characteristic well known solitonic solutions to the studied equation.

For solutions dened on general real algebraic curves, i.e., curves (1.1) with all a nm real, an important point in applications are reality and smoothness conditions. These are conveniently formulated for a homology basis for which the A-cycles are invariant under the action of the anti-holomorphic involution. However, the existing algorithms for the computational treatment of algebraic curves produce a basis of the homology that is in general not related to possible automorphisms of the curve. To implement the reality and smoothness requirements, a transformation to the basis for which the conditions are formulated has to be constructed. We study the necessary symplectic transformations and give explicit relations for so-called M-curves, curves with the maximal number of real ovals.

To illustrate these concepts, we study for the rst time numerically theta-functional solutions to integrable equations from the family of NLS equations, namely, the multicomponent nonlinear Schrödinger equation

(1.2) i ∂ψ j ∂t + ∂ 2 ψ j ∂x 2 + 2 n k=1 s k |ψ k | 2 ψ j = 0, j = 1, . . . , n,
denoted by n-NLS s , where s = (s 1 , . . . , s n ), s k = ±1, and the (2 + 1)-dimensional Davey-Stewartson (DS) equations,

i ψ t + ψ xx -α 2 ψ yy + 2 (Φ + ρ |ψ| 2 ) ψ = 0, Φ xx + α 2 Φ yy + 2ρ |ψ| 2 xx = 0, (1.3)
where α = i or α = 1 and where ρ = ±1. Both equations (1.2) and (1.3) reduce to the NLS equation under certain conditions: the former obviously in the case n = 1, the latter if ψ is independent of the variable y and satises certain boundary conditions, for instance that Φ + ρ |ψ| 2 tends to zero when x tends to innity. Integrability of the NLS equation was shown by Zakharov and Shabat [START_REF] Zakharov | Exact theory of two-dimensional self-focusing and one-dimensional selfmodulation of waves in nonlinear media[END_REF] and algebrogeometric solutions were given by Its [START_REF] Its | Inversion of hyperelliptic integrals and integration of nonlinear dierential equations[END_REF]. The multi-component nonlinear Schrödinger equation (1.2) in the case n = 2, s = (1, 1), is called the vector NLS or Manakov system. Manakov [START_REF] Manakov | On the theory of two-dimensional stationary self-focusing of electromagnetic waves[END_REF] rst examined this equation as an asymptotic model for the propagation of the electric eld in a waveguide. Its integrability was shown for n = 2 by Zakharov and Schulman in [START_REF] Zakharov | To the integrability of the system of two coupled nonlinear Schrödinger equations[END_REF] and for the general case in [START_REF] Radhakrishnan | Integrability and singularity structure of coupled nonlinear Schrödinger equations[END_REF]. Algebro-geometric solutions to the 2-NLS s equation with s = (1, 1) were presented in [START_REF] Elgin | Eective integration of the nonlinear vector Schrödinger equation[END_REF], and for the general case in [START_REF] Kalla | New degeneration of Fay's identity and its application to integrable systems[END_REF]. The DS equation (1.3) was introduced in [START_REF] Davey | On three-dimensional packets of surface waves[END_REF] to describe the evolution of a three-dimensional wave packet on water of nite depth. Its integrability was shown in [START_REF] Anker | On the Soliton Solutions of the Davey-Stewartson Equation for Long Waves[END_REF], and solutions in terms of multi-dimensional theta functions on general Riemann surfaces were given in [START_REF] Malanyuk | Finite-gap solutions of the Davey-Stewartson equations[END_REF][START_REF] Kalla | New degeneration of Fay's identity and its application to integrable systems[END_REF].

To ensure the correct numerical implementation of the formulae of [START_REF] Kalla | New degeneration of Fay's identity and its application to integrable systems[END_REF], we check for each point in the spacetime whether certain identities for theta functions are satised. Since these identities are not used in the code, they provide a strong test for the computed quantities. Numerically the identities are never exactly satised, but to high precision. The code reports a warning if the residual of the test relations is larger than 10 -6 which is well below plotting accuracy. Typically the conditions are satised to machine precision 1 . In addition we compute the solutions on a numerical grid and numerically dierentiate them. We check in this way for low genus that the solutions to n-NLS s and DS in terms of multi-dimensional theta functions satisfy the respective equations to better than 10 -6 . These two completely independent tests ensure that the presented plots are showing the correct solutions to better than plotting accuracy.

The paper is organized as follows: in Section 2 we recall some facts from the theory of multi-dimensional theta functions and the theory of real Riemann surfaces, necessary to give almost periodic solutions to the n-NLS s and DS equations. In Section 3 we consider the hyperelliptic case and study concrete examples of low genus, also in almost degenerate situations. In Section 4 we consider examples of non-hyperelliptic real Riemann surfaces and discuss symplectic transformations needed to obtain smooth solutions. We add some concluding remarks in Section 5.

Theta functions and real Riemann surfaces

In this section we recall basic facts on Riemann surfaces, in particular real surfaces, and multi-dimensional theta functions dened on them. Almost periodic solutions to the n-NLS s and the DS equations in terms of theta functions will be given following [START_REF] Kalla | New degeneration of Fay's identity and its application to integrable systems[END_REF].

2.1. Theta functions. Let R g be a compact Riemann surface of genus g > 0. Denote by (A, B) := (A 1 , . . . , A g , B 1 , . . . , B g ) a canonical homology basis, and by (ω 1 , . . . , ω g ) the basis of holomorphic dierentials normalized via (2.1)

A k ω j = 2iπδ kj , k, j = 1, . . . , g.
The matrix B with entries B kj = B k ω j of B-periods of the normalized holomorphic differentials ω j , j = 1, . . . , g, is symmetric and has a negative denite real part. The theta 1

We work with double precision, i.e., a precision of 10 -16 ; due to rounding errors this is typically reduced to 10 -12 to 10 -14 .

function with (half integer) characteristic

δ = [δ 1 , δ 2 ] is dened by (2.2) Θ B [δ](z) = m∈Z g exp 1 2 B(m + δ 1 ), m + δ 1 + m + δ 1 , z + 2iπδ 2 ,
for any z ∈ C g ; here δ 1 , δ 2 ∈ 0, 1 2 g are the vectors of the characteristic δ; ., . denotes the scalar product u, v = i u i v i for any u, v ∈ C g . The theta function Θ[δ](z) is even if the characteristic δ is even, i.e., if 4 δ 1 , δ 2 is even, and odd if the characteristic δ is odd, i.e., if 4 δ 1 , δ 2 is odd. An even characteristic is called non-singular if Θ[δ](0) = 0, and an odd characteristic is called non-singular if the gradient ∇Θ[δ](0) is non-zero. The theta function with characteristic is related to the theta function with zero characteristic (the Riemann theta function denoted by Θ) as follows (2.3)

Θ[δ](z) = Θ(z + 2iπδ 2 + Bδ 1 ) exp 1 2 Bδ 1 , δ 1 + z + 2iπδ 2 , δ 1 .
Denote by Λ the lattice Λ = {2iπN + BM, N, M ∈ Z g } generated by the A and Bperiods of the normalized holomorphic dierentials ω j , j = 1, . . . , g. The complex torus J(R g ) = C g /Λ is called the Jacobian of the Riemann surface R g . The theta function (2.2) has the following quasi-periodicity property with respect to the lattice Λ:

Θ[δ](z + 2iπN + BM) (2.4) = Θ[δ](z) exp -1 2 BM, M -z, M + 2iπ( δ 1 , N -δ 2 , M ) .
For the formulation of solutions to physically relevant integrable equations in terms of multi-dimensional theta functions, there is typically a preferred homology basis in which the solution takes a simple form. Let (A, B) and ( Ã, B) be arbitrary canonical homology basis dened on R g , represented here by 2g-dimensional vectors. Under the change of homology basis (2.5)

A B C D Ã B = A B ,
where

A B C D ∈ Sp(2g, Z
) is a symplectic matrix, the theta function (2.2) transforms as (2.6)

Θ B [δ](z) = κ det K exp 1 2 zt K-1 B z Θ B[ δ](z),
where K = 2iπA + B B and

B = 2iπ (2iπ C + D B) K-1 , (2.7)
z = (2iπ) -1 K z, (2.8
)

δ 1 δ 2 = A -B -C D δ1 δ2 + 1 2 Diag B A t D C t , (2.9)
for any z ∈ C g , where Diag(.) denotes the column vector of the diagonal entries of the matrix. Here κ is a constant independent of z and B (the exact value of κ is not needed for our purposes).

The Abel map R g -→ J(R g ) is dened by (2.10)

p p 0 := p p 0 ω,
for any p ∈ R g , where p 0 ∈ R g is the base point of the application, and where ω = (ω 1 , . . . , ω g ) t is the vector of the normalized holomorphic dierentials. Now let k a denote a local parameter near a ∈ R g and consider the following expansion of the normalized holomorphic dierentials ω j , j = 1, . . . , g, (2.11) ω j (p) = (V a,j + W a,j k a (p) + . . .) dk a (p), for any point p ∈ R g lying in a neighbourhood of a, where V a,j , W a,j ∈ C. Let us denote by D a (resp. D a ) the operator of the directional derivative along the vector V a = (V a,1 , . . . , V a,g ) t (resp. W a ). According to [START_REF] Mumford | Tata Lectures on Theta. I and II[END_REF] and [START_REF] Kalla | New degeneration of Fay's identity and its application to integrable systems[END_REF], the theta function satises the following identities derived from Fay's identity [START_REF] Fay | Theta functions on Riemann surfaces[END_REF]:

(2.12)

D a D b ln Θ(z) = q 1 + q 2 Θ(z + b a ) Θ(z - b a ) Θ(z) 2 ,
(2.13)

D a ln Θ(z + b a ) Θ(z) + D 2 a ln Θ(z + b a ) Θ(z) + D a ln Θ(z + b a ) Θ(z) -K 1 2 + 2D 2 a ln Θ(z) + K 2 = 0,
for any z ∈ C g and any distinct points a, b ∈ R g ; here the scalars q i , K i , i = 1, 2 depend on the points a, b and are given by (2.14)

q 1 (a, b) = D a D b ln Θ[δ]( b a ), (2.15)
q 2 (a, b) = D a Θ[δ](0) D b Θ[δ](0) Θ[δ]( b a ) 2 , (2.16
)

K 1 (a, b) = 1 2 D a Θ[δ](0) D a Θ[δ](0) + D a ln Θ[δ]( b a ) , (2.17) K 2 (a, b) = -D a ln Θ( b a ) -D 2 a ln Θ( b a ) Θ(0) -D a ln Θ( b a ) -K 1 (a, b) 2 ,
where δ is a non-singular odd characteristic.

2.2. Real Riemann surfaces. A Riemann surface R g is called real if it admits an antiholomorphic involution τ : R g → R g , τ 2 = id. The connected components of the set of xed points of the anti-involution τ are called real ovals of τ . We denote by R g (R) the set of xed points. According to Harnack's inequality [START_REF] Harnack | Ueber die Vieltheiligkeit der ebenen algebraischen Curven[END_REF], the number k of real ovals of a real Riemann surface of genus g cannot exceed g + 1: 0 ≤ k ≤ g + 1. Curves with the maximal number k = g + 1 of real ovals are called M-curves. The complement R g \R g (R) has either one or two connected components. The curve R g is called a dividing curve if R g \ R g (R) has two components, and

R g is called non-dividing if R g \ R g (R)
is connected (notice that an M-curve is always a dividing curve).

µ 2 = 2g+2 i=1 (λ -λ i ),
where the branch points λ i ∈ R are ordered such that λ 1 < . . . < λ 2g+2 . On such a curve, we can dene two anti-holomorphic involutions τ 1 and τ 2 , given respectively by τ 1 (λ, µ) = (λ, µ) and τ 2 (λ, µ) = (λ, -µ). Projections of real ovals of τ 1 on the λ-plane coincide with the intervals [λ 2g+2 , λ 1 ], . . . , [λ 2g , λ 2g+1 ], whereas projections of real ovals of τ 2 on the λ-plane coincide with the intervals [λ 1 , λ 2 ], . . . , [λ 2g+1 , λ 2g+2 ]. Hence the curve (2.18) is an M-curve with respect to both anti-involutions τ 1 and τ 2 .

Let (A, B) be a basis of the homology group H 1 (R g ). According to Proposition 2.2 in Vinnikov's paper [START_REF] Vinnikov | Self-adjoint determinantal representations of real plane curves[END_REF], there exists a canonical homology basis (that we call for simplicity `Vinnikov basis' in the following) such that (2.19)

τ A τ B = I g 0 H -I g A B ,
where I g is the g × g unit matrix, and H is a block diagonal g × g matrix, dened as follows:

1) if R g (R) = ∅, H =              0 1 1 0 . . . 0 1 1 0 0 . . . 0              if R g is dividing, H =          1 . . . 1 0 . . . 0          if R g is non-dividing; rank(H) = g + 1 -k in both cases.
2) if R g (R) = ∅, (i.e. the curve does not have real oval), then

H =        0 1 1 0 . . . 0 1 1 0        or H =          0 1 1 0 . . . 0 1 1 0 0          ; rank(H) = g if g is even, rank(H) = g -1 if g is odd. Now let us choose the canonical homology basis in H 1 (R g ) satisfying (2.19
), take a, b ∈ R g and assume that τ a = a and τ b = b. Denote by a contour connecting the points a and b which does not intersect the canonical homology basis. Then the action of τ on the generators (A, B, ) of the relative homology group H 1 (R g , {a, b}) is given by

(2.20)   τ A τ B τ   =   I g 0 0 H -I g 0 N t M t 1     A B   ,
where the vectors N, M ∈ Z g are related by (see [START_REF] Kalla | New degeneration of Fay's identity and its application to integrable systems[END_REF]) (2.21)

2 N + HM = 0.
2.3. Theta-functional solutions of the n-NLS s equation. Algebro-geometric data associated to smooth theta-functional solutions of the n-NLS s equation (1.2) consist of {R g , τ, f, z a }, where R g is a compact Riemann surface of genus g > 0 dividing with respect to an anti-holomorphic involution τ , and admitting a real meromorphic function f of degree n + 1; here z a ∈ R is a non critical value of f such that the ber f -1 (z a ) = {a 1 , . . . , a n+1 } over z a belongs to the set R g (R). Local parameters k a j near a j are dened by k a j (p) = f (p) -z a , for any point p ∈ R g lying in a neighbourhood of a j .

Denote by (A, B, j ) the generators of the relative homology group H 1 (R g , {a n+1 , a j }). Let d ∈ R g and θ ∈ R. Then the following functions ψ j , j = 1, . . . , n, give smooth solutions of the n-NLS s equation (1.2), see [START_REF] Kalla | New degeneration of Fay's identity and its application to integrable systems[END_REF],

(2.22)

ψ j (x, t) = |A j | e iθ Θ(Z -d + r j ) Θ(Z -d) exp {-i (E j x -F j t)} ,
where

|A j | = |q 2 (a n+1 , a j )| 1/2 exp 1 2 d, M j .
The vector M j ∈ Z g is dened by the action of τ on the relative homology group H 1 (R g , {a n+1 , a j }) (see (2.20)). Moreover, r j = j ω, and the vector Z reads

Z = i V a n+1 x + i W a n+1 t,
where vectors V a n+1 and W a n+1 are dened in (2.11). The scalars E j , F j are given by (2.23)

E j = K 1 (a n+1 , a j ), F j = K 2 (a n+1 , a j ) -2 n k=1 q 1 (a n+1 , a k ),
and scalars q i , K i , i = 1, 2 are dened in (2.14)-(2.17). According to [START_REF] Kalla | New degeneration of Fay's identity and its application to integrable systems[END_REF], necessary conditions for the functions ψ j in (2.22) to solve the n-NLS s equation are the identities (2.12) and (2.13) with (a, b) := (a n+1 , a j ).

The signs s 1 , . . . , s n in (1.2) are given by (2.24)

s j = exp {iπ(1 + α j )} ,
where α j ∈ Z denote certain intersection indices on R g dened as follows: let ãn+1 , ãj ∈ R g (R) lie in a neighbourhood of a n+1 and a j respectively such that f (ã n+1 ) = f (ã j ). Denote by ˜ j an oriented contour connecting ãn+1 and ãj . Then (2.25)

α j = (τ ˜ j -˜ j ) • j
is the intersection index of the closed contour τ ˜ j -˜ j and the contour j ; this index is computed in the relative homology group H 1 (R g , {a n+1 , a j }).

In particular, it was shown in [START_REF] Kalla | New degeneration of Fay's identity and its application to integrable systems[END_REF] that solutions of the focusing n-NLS s equation, i.e., for s = (1, . . . , 1), are obtained when the branch points of the meromorphic function f are pairwise conjugate.

Theta-functional solutions of the DS equations. Now let us introduce smooth solutions of the DS equations. In characteristic coordinates

ξ = 1 2 (x -iα y), η = 1 2 (x + iα y), α = i or 1,
the DS equations (1.3) take the form

iψ t + 1 2 (∂ ξξ + ∂ ηη )ψ + 2 φ ψ = 0, ∂ ξ ∂ η φ + ρ 1 2 (∂ ξξ + ∂ ηη )|ψ| 2 = 0, (2.26)
where φ := Φ + ρ |ψ| 2 , ρ = ±1. Recall that DS1 ρ denotes the Davey-Stewartson equation when α = i (in this case ξ and η are both real), and DS2 ρ when α = 1 (in this case ξ and η are pairwise conjugate).

In both cases, for the DS1 ρ and DS2 ρ equations, the solutions have the form [START_REF] Malanyuk | Finite-gap solutions of the Davey-Stewartson equations[END_REF][START_REF] Kalla | New degeneration of Fay's identity and its application to integrable systems[END_REF]:

(2.27)

ψ(ξ, η, t) = |A| e iθ Θ(Z -d + r) Θ(Z -d) exp -i G 1 ξ + G 2 η -G 3 t 2 ,
(2.28)

φ(ξ, η, t) = 1 2 (ln Θ(Z -d)) ξξ + 1 2 (ln Θ(Z -d)) ηη + h 4 .
Here r = b a ω for some distinct points a, b ∈ R g , and the vector Z is dened as (2.29)

Z = i κ 1 V a ξ -κ 2 V b η + (κ 2 1 W a -κ 2 2 W b ) t 2 .
Moreover, the scalars G 1 , G 2 and G 3 read (2.30)

G 1 = κ 1 K 1 (a, b), G 2 = κ 2 K 1 (b, a), (2.31)
G 3 = κ 2 1 K 2 (a, b) + κ 2 2 K 2 (b, a) + h,
where the scalars K 1 , K 2 are dened in (2.16) and (2.17). As shown in [START_REF] Kalla | New degeneration of Fay's identity and its application to integrable systems[END_REF], necessary conditions for the functions ψ (2.27) and φ (2.28) to solve the DS equations are the identities (2.12) and (2.13).

Algebro-geometric data associated to smooth solutions (2.27), (2.28) of the DS1 ρ equation consist of {R g , τ, a, b, k a , k b }, where R g is a compact Riemann surface of genus g > 0, dividing with respect to an anti-holomorphic involution τ , a, b are two distinct points in R g (R), and k a , k b denote local parameters near a and b respectively which satisfy k a (τ p) = k a (p) for any p lying in a neighbourhood of a, and k b (τ p) = k b (p) for any p lying in a neighbourhood of b. The remaining quantities satisfy the conditions:

d ∈ R g , θ, h ∈ R, κ 2 ∈ R \ {0},
(2.32)

κ 1 = -ρ κ2 1 κ 2 q 2 (a, b) exp 1 2 BM, M + r + d, M ,
for some κ1 ∈ R, where M ∈ Z g is dened in (2.20). The scalar |A| is given by

|A| = |κ 1 κ 2 q 2 (a, b)| exp { d, M } ,
where the quantity q 2 is dened in (2.15). Algebro-geometric data associated to smooth solutions (2.27), (2.28) of the DS2 ρ equation consist of {R g , τ, a, b, k a , k b }, where R g is a compact Riemann surface of genus g > 0 with an anti-holomorphic involution τ , a, b are two distinct points such that τ a = b, and k a , k b denote local parameters near a and b respectively which satisfy k b (τ p) = k a (p) for any point p lying in a neighbourhood of a. Moreover,

d ∈ iR g , θ, h ∈ R, κ 1 , κ 2 ∈ C \ {0} satisfy κ 1 = κ 2 ,
and the scalar |A| is given by

|A| = |κ 1 | |q 2 (a, b)| 1/2 .
Smooth solutions of the DS2 + equation are obtained when the curve R g is an M-curve with respect to τ , whereas solutions to DS2 -are smooth if the associated Riemann surface does not have real oval with respect to τ , and if there is no pseudo-real function of degree g -1 on it (i.e., function which satises f (τ p) = -f (p) -1 ), see [START_REF] Malanyuk | Finite-gap solutions of the Davey-Stewartson equations[END_REF].

Remark 2.1. The symmetric structure of the DS equations (2.26) with respect to ξ and η implies that a solution ψ = Ψ(ξ, η, t) to DS1 + leads to a solution Ψ(-ξ, η, t) of DS1 -.

Hyperelliptic case

Here we consider concrete examples for the solutions, in terms of multi-dimensional theta functions, to DS and n-NLS s on hyperelliptic Riemann surfaces. We rst review the numerical methods to visualize the solutions and discuss how the accuracy is tested.

3.1. Computation on real hyperelliptic curves. The simplest example of algebraic curves are hyperelliptic curves,

µ 2 = 2g+2 i=1 (λ -λ i ), without branching at innity 2g+1 i=1 (λ -λ i )
, with branching at innity , where g is the genus of the Riemann surface, and where we have for the branch points λ i ∈ C the relations λ i = λ j for i = j. If the number of nite branch points is odd, the curve is branched at innity. Recall that all Riemann surfaces of genus 2 are hyperelliptic, and that the involution σ which interchanges the sheets, σ(λ, µ) = (λ, -µ), is an automorphism on any hyperelliptic curve. A vector of holomorphic dierentials for these surfaces is given by (1, λ, . . . , λ g-1 ) t dλ/µ. For a real hyperelliptic curve, the branch points are either real or pairwise conjugate. As we saw in Example 2.1, if all branch points λ i are real and ordered such that λ 1 < . . . < λ 2g+2 , the hyperelliptic curve is an M-curve with respect to both anti-holomorphic involutions τ 1 and τ 2 dened in the example. The other case of interest in the context of smooth solutions to n-NLS s and DS are real curves without real branch point. For the involution τ 1 , a curve given by µ 2 = g+1 i=1 (λ -λ i )(λ -λ i ), with λ i ∈ C\R, i = 1, . . . , g+1, in this case is dividing (two points whose projections onto C have respectively a positive and a negative imaginary part cannot be connected by a contour which does not cross a real oval), whereas a curve given by µ

2 = -g+1 i=1 (λ -λ i )(λ -λ i )
has no real oval, and vice versa for the involution τ 2 .

In the following, we will only consider real hyperelliptic curves without branching at innity and write the dening equation in the form

µ 2 = (λ-ξ)(λ-η) g i=1 (λ-E i )(λ-F i ).
It is possible to introduce a convenient homology basis on the related surfaces, see Fig. 1 for the case η = ξ. The simple form of the algebraic relation between µ and λ for hyperelliptic curves makes the generation of very ecient numerical codes possible, see, for instance, [START_REF] Frauendiener | Hyperelliptic theta functions and spectral methods[END_REF][START_REF] Frauendiener | Hyperelliptic theta functions and spectral methods: KdV and KP solutions[END_REF] for details. These codes allow the treatment of almost degenerate Riemann surfaces, i.e., the case where the branch points almost collide pairwise, where the distance of the branch points is of the order of machine precision: |E i -F i | ∼ 10 -14 . The homology basis Fig. 1 is adapted to this kind of degeneration.

N .. E 1 F N+1 E g F g E N+1 F 1 E N F
The Abel map b a ω between two points a and b is computed in the following way: the sheet identied at the point a = (λ(a), µ(a)) (where we take for µ the root computed by Matlab) is labeled sheet 1, and at the point (λ(a), -µ(a)), sheet 2. Then the ramication point whose projection to the λ-sphere has the minimal distance to λ(a) is determined. For simplicity we assume always that this is the point ξ in Fig. 1 (for another branch point, this leads to the addition of half-periods, see e.g. [START_REF] Belokolos | Algebro-geometric approach to nonlinear integrable equations[END_REF]). This means we compute b a ω as

b a ω = b ξ ω - a ξ ω.
The choice of a branch point as the base point of the Abel map has the advantage that a change of sheet of a point a just implies a change of sign of the integral:

(λ(a),µ(a)) ξ ω = - (λ(a),-µ(a)) ξ ω.
To compute the integral a ξ ω, one has to analytically continue µ on the connecting line between λ(a) and ξ onto the λ-sphere. Whereas the root µ is not supposed to have any branching on the considered path, the square root in Matlab is branched on the negative real axis. To analytically continue µ on the path [λ(a), ξ], we compute the Matlab root at some λ j ∈ [λ(a), ξ], j = 0, . . . , N c and analytically continue it starting from µ(a) by demanding that |µ(λ j+1 ) -µ(λ j )| < |µ(λ j+1 ) + µ(λ j )|. The so dened sheets will be denoted here and in the following by numbers, i.e., a point on sheet 1 with projection λ(a) into the base is denoted by (λ(a)) (1) .

Thus the computation of the Abel map is reduced to the computation of line integrals on the connecting line between λ(a) and ξ in the complex λ-plane. For the numerical computation of such integrals we use Clenshaw-Curtis integration (see, for instance, [START_REF] Trefethen | Spectral Methods in Matlab[END_REF]): to compute an integral 1 -1 h(x) dx, this algorithm samples the integrand on the N c + 1 Chebyshev collocation points x j = cos(jπ/N c ), j = 0, . . . , N c . The integral is approximated as the sum: 1 -1 h(x) dx ∼ Nc j=0 w j h(x j ) (see [START_REF] Trefethen | Spectral Methods in Matlab[END_REF] on how to obtain the weights w j ). It can be shown that the convergence of the integral is exponential for analytic functions h as the ones considered here. To compute the Abel map, one uses the transformation λ → λ(a)(1 + x)/2 + ξ(1 -x)/2, to the Clenshaw-Curtis integration variable. The same procedure is then carried out for the integral from ξ to b.

The theta functions are approximated as in [START_REF] Frauendiener | Hyperelliptic theta functions and spectral methods[END_REF] as a sum,

(3.1) Θ B [δ](z) ∼ N θ m 1 =-N θ . . . N θ mg=-N θ exp 1 2 B(m + δ 1 ), m + δ 1 + m + δ 1 , z + 2iπδ 2 .
The periodicity properties of the theta function (2.4) make it possible to write z = z 0 + 2iπN + BM for some N, M ∈ Z g , where

z 0 = 2iπα + Bβ with α i , β i ∈ ] -1 2 , 1 2
] for i = 1, . . . , g. The value of N θ is determined by the condition that all terms in (2.2) with |m i | > N θ are smaller than machine precision, which is controlled by the largest eigenvalue of the real part of the Riemann matrix (the one with minimal absolute value since the real part is negative denite), see [START_REF] Frauendiener | Hyperelliptic theta functions and spectral methods[END_REF][START_REF] Frauendiener | Algebraic curves and Riemann surfaces in Matlab[END_REF].

To control the accuracy of the numerical solutions, we use essentially two approaches. First we check the theta identity (2.13), which is the underlying reason for the studied functions being solutions to n-NLS s and DS, at each point in the spacetime. This test requires the computation of theta derivatives not needed in the solution (which slightly reduces the eciency of the code since additional quantities are computed), but provides an immediate check whether the solution satises (2.13) with the required accuracy. Since this identity is not implemented in the code, it provides a strong test. This ensures that all quantities entering the solution are computed with the necessary precision. In addition, the solutions are computed on Chebyshev collocation points (see, for instance, [START_REF] Trefethen | Spectral Methods in Matlab[END_REF]) for each of the physical variables. This can be used for an expansion of the computed solution in terms of Chebyshev polynomials, a so-called spectral method having in practice exponential convergence for analytic functions as the ones considered here. Since the derivatives of the Chebyshev polynomials can be expressed linearly in terms of Chebyshev polynomials, a derivative acts on the space of polynomials via a so called dierentiation matrix. With these standard Chebyshev dierentiation matrices (see [START_REF] Trefethen | Spectral Methods in Matlab[END_REF]), the solution can be numerically dierentiated. The computed derivatives allow to check with which numerical precision the PDE is satised by a numerical solution. With these two independent tests, we ensure that the shown solutions are correct to much better than plotting accuracy (the code reports a warning if the above tests are not satised to better than 10 -6 ).

3.2. Solutions to the DS equations. The elliptic solutions are the well known travelling wave solutions and will not be discussed here. The simplest examples we will consider for the DS solutions are given on hyperelliptic curves of genus 2. As we saw in Section 2.4, for DS1 ρ reality and smoothness conditions imply that the branch points of the curve are either all real (M-curve) or all pairwise conjugate (dividing curve). The points a and b must project to real points on the λ-sphere and must be stable under the anti-holomorphic involution τ (we use here τ = τ 1 , as dened in Example 2.1, except for DS2 -). For DS2 ρ , we have τ a = b where the projection of a onto the λ-sphere is the conjugate of the projection of b. For DS2 + the curve must have only real branch points (M-curve), whereas for DS2 -it must have no real oval.

For DS we will mainly give plots for xed time since for low genus, the solution is essentially travelling in one direction. For higher genus, we show a more interesting time dependence in Fig. 9.

We rst consider the defocusing variants, DS1 + and DS2 + on M-curves. In genus 2 we study the family of curves with the branch points -2, -1, 0, , 2, 2 + for = 1 and = 10 -10 . In the former case the solutions will be periodic in the (x, y)-plane, in the latter almost solitonic since the Riemann surface is almost degenerate (in the limit → 0 the surface degenerates to a surface of genus 0; the resulting solutions are discussed in more detail in [START_REF] Kalla | Breathers and solitons of generalized nonlinear Schrödinger equations as degenerations of algebro-geometric solutions[END_REF]). To obtain non-trivial solutions in the solitonic limit, we use

d = 1 2 [ 1 1 0 0 ]
t in all examples. In Fig. 2 it can be seen that these are in fact dark solitons, i.e., the solutions tend asymptotically to a non-zero constant and the solitons thus represent `shadows' on a background of light. The well known features from soliton collisions for (1+1)-dimensional integrable equations, namely, the propagation without change of shape, and the unchanged shape and phase shift after the collision, can be seen here in the (x, y)-plane.

The corresponding solutions to DS2 + can be seen in Fig. 3. We only show the square modulus of the solution here for simplicity. For the real and the imaginary part of such a solution for the DS1 --case, see Fig. 6. Because of remark 2.1 all solutions shown for DS1 + on M-curves are after the change of coordinate ξ → -ξ solutions to DS1 -. For this reason DS1 -solutions on M-curves will not be presented here.

In the same way one can study, on a genus 4 hyperelliptic curve, the formation of the dark 4-soliton for these two equations. We consider the curve with branch points -4, -3, -2, -2 + , 0, , 2, 2 + , 4, 4 + for = 1 and = 10 -10 , and use d = 1 2 [ 1 1 1 1 0 0 0 0 ] t . The DS1 + solutions for this curve can be seen in Fig. 4. The corresponding solutions to DS2 + is shown in Fig. 5. 2. Solution (2.27) to the DS1 + equation at t = 0 on a hyperelliptic curve of genus 2 with branch points -2, -1, 0, , 2, 2 + and a = (-1.9) (1) , b = (-1.1) (2) for = 1 on the left and = 10 -10 , the almost solitonic limit, on the right. (2) for = 1 on the left and = 10 -10 , the almost solitonic limit, on the right.

Solutions to the focusing variants of these equations can be obtained on hyperelliptic curves with pairwise conjugate branch points. For such solutions the solitonic limit cannot be obtained as above since the quotient of theta functions in (2.27) tends to a constant in this case. To obtain the well-known bright solitons (solutions tend to zero at spatial innity) in this way, the hyperelliptic curve has to be completely degenerated (all branch points must collide pairwise to double points) which leads to limits of the form '0/0' in the expression for the solution (2.27) which are not convenient for a numerical treatment; see [START_REF] Kalla | Breathers and solitons of generalized nonlinear Schrödinger equations as degenerations of algebro-geometric solutions[END_REF] for an analytic discussion. Therefore we only consider non-degenerate hyperelliptic curves here. To obtain smooth solutions, we use d = 0. A solution in genus 2 of the DS1 - equation is studied on the curve with the branch points -2 ± i, -1 ± i, 1 ± i in Fig. 6.

A typical example of a DS1 -solution on a hyperelliptic curve of genus 4 with branch points -2 ± i, -1 ± i, ±i, 1 ± i, 2 ± i is shown in Fig. 7. and a = (-3.9) (1) , b = (-3.1) (2) for = 1 on the left and = 10 -10 , the almost solitonic limit, on the right. Figure 5. Solution (2.27) to the DS2 + equation at t = 0 on a hyperelliptic curve of genus 4 with branch points -4, -3, -2, -2 + , 0, , 2, 2 + , 4, 4 + and a = (-1.5+2i) (1) , b = (-1.5-2i) (1) for = 1 on the left and = 10 -10 , the almost solitonic limit, on the right.

Smooth solutions to DS2 -can be obtained on Riemann surfaces without real oval for points a and b satisfying τ a = b. As mentioned above, hyperelliptic curves of the form

µ 2 = -2g+2 i=1 (λ -λ i )
with pairwise conjugate branch points have no real oval for the standard involution τ 1 as dened in Example 2.1. On the other hand, surfaces dened by the algebraic equation µ 2 = 2g+2 i=1 (λ -λ i ) have no real oval for the involution τ 2 (see Example 2.1). We will consider here τ 2 for the same curves as for DS1 -. An example for genus 2 can be seen in Fig. 8. An example for a DS2 -solution of genus 4 can be seen in Fig. 9.

3.3. Solutions to the n-NLS s equation. A straightforward way to obtain solutions (2.22) to the n-NLS s equation is given on an (n + 1)-sheeted branched covering of the complex plane, an approach that will be studied in more detail in the next section. As can be seen from the proof of Theorem 4.1 in [START_REF] Kalla | New degeneration of Fay's identity and its application to integrable systems[END_REF], the crucial point in the construction of these Figure 6. Solution (2.27) to the DS1 -equation at t = 0 on a hyperelliptic curve of genus 2 with branch points -2 ± i, -1 ± i, 1 ± i and a = (-4) (1) , b = (-3) (2) . The square modulus of the solution is shown on the left, real and imaginary parts on the right. 7. Solution (2.27) to the DS1 -equation at t = 0 on a hyperelliptic curve of genus 4 with branch points -2 ± i, -1 ± i, ±i, 1 ± i, 2 ± i and a = (-4) (1) , b = (-3) (2) . solutions is the fact that n+1 k=1 V a k = 0. This implies that it is also possible to construct theta-functional n-NLS s solutions on hyperelliptic surfaces by introducing constants γ k via n+1 k=1 γ k V a k = 0 in the following corollary of Theorem 4.1 in [START_REF] Kalla | New degeneration of Fay's identity and its application to integrable systems[END_REF]: Corollary 3.1. Let R g be a real hyperelliptic curve of genus g > 0 and denote by τ an anti-holomorphic involution. Choose the canonical homology basis which satises (2.19). Take n ≥ g and let a 1 , . . . , a n+1 ∈ R g (R) not ramication points having distinct projection λ(a j ), j = 1, . . . , n + 1, onto the λ-sphere. Denote by j an oriented contour between a n+1 and a j which does not intersect cycles of the canonical homology basis. Let d R ∈ R g , T ∈ Z g , and dene d = d R + iπ 2 (diag(H) -2 T). Take θ ∈ R and let γ g+1 , . . . , γ n ∈ R be arbitrary constants with γ n+1 = 1. Put ŝ = (sign(γ 1 ) s 1 , . . . , sign(γ n ) s n ) where s j is Figure 8. Solution to the DS2 -equation at t = 0 on a hyperelliptic curve of genus 2 with branch points -2 ± i, -1 ± i, 1 ± i and a = (-1.5 + 2i) (1) , b = (-1.5 -2i) (2) . Figure 9. Solution to the DS2 -equation for several values of t on a hyperelliptic curve of genus 4 with branch points -2 ± i, -1 ± i, ±i, 1 ± i, 2 ± i and a = (-1.5 + 2i) (1) , b = (-1.5 -2i) (2) .

given in (2.24), and the scalars γ j , j = 1, . . . , g, are dened by n+1 k=1 γ k V a k = 0. Then the following functions ψ j , j = 1, . . . , n, give solutions of the n-NLS ŝ equation (1.2) (3.2)

ψ j (x, t) = |γ j | 1/2 |A j | e iθ Θ(Z -d + r j ) Θ(Z -d) exp{-i (E j x -F j t)},
where

|A j | = |q 2 (a n+1 , a j )| 1/2 exp 1 2 d, M j . Here Z = i V a n+1 x + i W a n+1 t,
where the vectors V a n+1 and W a n+1 were introduced in (2.11), and r j = j ω. The scalars E j , F j are given by

E j = K 1 (a n+1 , a j ), F j = K 2 (a n+1 , a j ) -2 n k=1 γ k q 1 (a n+1 , a k ),
where q i , K i for i = 1, 2 are dened in (2.14)-(2.17). If R g is dividing and if d ∈ R g , functions (3.2) give smooth solutions of n-NLS ŝ.

As an example we consider, as for DS in genus 2, the family of curves with the branch points -2, -1, 0, , 2, 2 + for = 1 and = 10 -10 . In the former case the solutions will be periodic in the (x, t)-plane, in the latter almost solitonic. To obtain non-trivial solutions in the solitonic limit, we use

d = 1 2 [ 1 1 0 0 ]
t in all examples. In Fig. 10 we show the case a 1 = (-1.9) (1) , a 2 = (-1.1) (1) and a 3 = (-1.8) (1) , which leads to a solution of 2-NLS ŝ with ŝ = (-1, -1). Interchanging a 2 and a 3 in the above example, we obtain a solution to 2-NLS ŝ with ŝ = (1, -1) in Fig. 11. -1, -1) on a hyperelliptic curve of genus 2 with branch points -2, -1, 0, , 2, 2 + and a 1 = (-1.9) (1) , a 2 = (-1.1) (1) and a 3 = (-1.8) (1) for = 1 on the left and = 10 -10 , the almost solitonic limit, on the right. Solutions of 4-NLS ŝ can be studied in the same way on the hyperelliptic curve of genus 4 with branch points -4, -3, -2, -2 + , 0, , 2, 2 + , 4, 4 + . We use

d = 1 2 [ 1 1 1 1 0 0 0 0 ]
t and the points a 1 = (-3.9) (1) , a 2 = (-3.7) (1) , a 3 = (-3.5) (1) , a 4 = (-3.3) (1) and a 5 = (-3.1) (1) . Since the vectors V a j and W a j are very similar in this case, the same is true for the functions ψ j . Therefore, we will only show the square modulus of the rst component ψ 1 in Fig. 12 for ŝ = (1, -1, 1, -1) on the left. Interchanging a 4 and a 5 in this case, one gets a solution to 4-NLS ŝ with ŝ = (-1, 1, -1, -1) which can be seen on the right of Fig. 12. The almost solitonic limit = 10 -10 produces well known solitonic patterns as shown for instance for the DS equation in the previous subsection. -1) on a hyperelliptic curve of genus 2 with branch points -2, -1, 0, , 2, 2 + and a 1 = (-1.9) (1) , a 2 = (-1.8) (1) and a 3 = (-1.1) (1) for = 1 on the left and = 10 -10 , the almost solitonic limit, on the right. for a 1 = (-3.9) (1) , a 2 = (-3.7) (1) , a 3 = (-3.5) (1) , a 4 = (-3.3) (1) and a 5 = (-3.1) (1) , which leads to ŝ = (1, -1, 1, -1), on the left, and for a 1 = (-3.9) (1) , a 2 = (-3.7) (1) , a 3 = (-3.5) (1) , a 4 = (-3.1) (1) and a 5 = (-3.3) (1) , which leads to ŝ = (-1, 1, -1, -1) on the right.

Hyperelliptic solutions to the n-NLS ŝ equation with all signs ŝj satisfying ŝj = 1, can be constructed on a curve without real branch points. To obtain smooth solutions, we use d = 0. A solution of the 2-NLS ŝ equation is studied on the curve of genus 2 with the branch points -2 ± i, -1 ± i, 1 ± i in Fig. 13.

A typical example for a hyperelliptic 4-NLS ŝ solution with ŝ = (1, 1, 1, 1) can be obtained on a curve of genus 4 with branch points -2 ± i, -1 ± i, ±i, 1 ± i, 2 ± i, as shown in Fig. 14.

General real algebraic curves

The quantities entering theta-functional solutions of the DS and n-NLS s equations are related to compact Riemann surfaces. Since all compact Riemann surfaces can be dened ) on a hyperelliptic curve of genus 2 with branch points -2±i, -1±i, 1±i and a 1 = (-1.9) (1) , a 2 = (-1.8) (2) and a 3 = (-1.1) (1) . ) on a hyperelliptic curve of genus 4 with branch points -2 ± i, -1 ± i, ±i, 1 ± i, 2 ± i and a 1 = (-3.9) (1) , a 2 = (-3.7) (2) , a 3 = (-3.5) (1) , a 4 = (-3.3) (2) and a 5 = (-3.1) (1) . via compactied non-singular algebraic curves, convenient computational approaches as [START_REF] Deconinck | Computing Riemann matrices of algebraic curves[END_REF][START_REF] Deconinck | Computing Riemann Theta Functions[END_REF] and [START_REF] Frauendiener | Algebraic curves and Riemann surfaces in Matlab[END_REF] are based on algebraic curves: dierentials, homology basis and periods of the Riemann surface can be obtained in an algorithmic way. We refer the reader to the cited literature for details. The identication of the sheets of the covering dened by the algebraic curve (1.1) via the projection map (x, y) → x, is done, as in the hyperelliptic case, by analytic continuation of the roots y i , i = 1, . . . , N for some non-critical point x b on the x-sphere, along a set of contours specied in [START_REF] Frauendiener | Algebraic curves and Riemann surfaces in Matlab[END_REF]. In the context of real algebraic curves for which solutions of n-NLS s and DS are discussed here, an additional problem is to establish the action of the anti-holomorphic involution τ on points on dierent sheets. A typical problem is to nd points a ∈ R g and b ∈ R g with the same projection onto the x-sphere such that τ a = b; here τ is dened via τ a = (x(a), y(a)). To this end, the roots y i , i = 1, . . . , N, identied at x = x b , are analytically continued to the points projecting to x(a) on the x-sphere. It is then established which pairs of points in the dierent sheets satisfy τ a = b.

In contrast to the hyperelliptic curves of the previous section, it is not possible for general curves to introduce a priori a basis of the homology. Thus the cited codes use an algorithm by Tretko and Tretko [START_REF] Tretko | Combinatorial group theory, Riemann surfaces and dierential equations[END_REF] which produces a homology basis for a given branching structure of the covering which is in general not adapted to possible automorphisms of the curve. In the context of theta-functional solutions to integrable PDEs one is often interested in real curves. As discussed in [START_REF] Kalla | New degeneration of Fay's identity and its application to integrable systems[END_REF], the Vinnikov basis (i.e., the canonical homology basis which satises (2.19)) is convenient in this context. Since solutions and smoothness conditions for n-NLS s and DS equations are formulated in this basis, a symplectic transformation relating the computed basis to the Vinnikov basis needs to be worked out. This transformation is discussed in the present section and will be applied to examples of real algebraic curves. 4.1. Symplectic transformation. Let R g be a real compact Riemann surface of genus g and τ an anti-holomorphic involution dened on it. Let (ν 1 , . . . , ν g ) be a basis of holomorphic dierentials such that (4.1)

τ * ν j = ν j , j = 1, . . . , g,
where τ * is the action of τ lifted to the space of holomorphic dierentials: τ * ω(p) = ω(τ p) for any p ∈ R g . For an arbitrary canonical homology basis (A, B), let us denote by P A and P B the matrices of A and B-periods of the dierentials ν j :

(4.2)

(P A ) ij = A i ν j , (P B ) ij = B i ν j , i, j = 1, . . . , g.
In what follows (A, B) denotes the Vinnikov basis. From (4.1) and (2.19) we deduce the action of the complex conjugation on the matrices P A and P B :

(4.3)

(P A ) ij ∈ R, (4.4)
P B = -P B + HP A .
Denote by ( Ã, B) the homology basis on R g produced by the Tretko-Tretko algorithm. From the symplectic transformation (2.5) we obtain the following transformation law between the matrices P Ã, P B and P A , P B dened in (4.2): According to (4.6), the matrix A Re P Ã +B Re P B is invertible, since the matrix P A of A-periods of a basis of holomorphic dierentials is always invertible (see, for instance, [START_REF]Computational Approach to Riemann Surfaces[END_REF]). Moreover, it is well known that the Riemann matrix B = 2iπ P B (P A ) -1 has a (negative) denite real part, which is equal to -2π Im(P B ) Im((P A ) -1 ) for the real matrix P A here. Then, by (4.9) the matrix C Im P Ã + D Im P B is also invertible. (4.10)

B t = -Im P Ã C Im P Ã + D Im P B -1
(4.11)

C t = 1 2 A t H -Re P B A Re P Ã + B Re P B -1
(4.12)

D t = 1 2 B t H + Re P Ã A Re P Ã + B Re P B -1 .
(4.13)

Proof. Recall that symplectic matrices

M = A B C D ∈ Sp(2g, Z) are characterized by A t D-C t B = I g , (4.14)
A t C = C t A, (4.15)
D t B = B t D. ( 4
A t C Re P Ã + A t D Re P B = 1 2 A t HP A C t A Re P Ã + (I g + C t B) Re P B = 1 2 A t HP A C t A Re P Ã + B Re P B = 1 2 A t HP A -Re P B ,
which by (4.6) leads to (4.12). Identity (4.13) can be proved analogously.

Remark 4.1. Lemma 4.1 implies that it is sucient to know the matrices A and B (or C and D) to determine the symplectic matrix in (4.5). In practice, this means that a convenient ansatz for one of the matrices has to be found. The others then follow from the relations in Lemma 4.1.

Thus to construct these matrices one rst checks which of the matrices Re P Ã , Re P B , Im P Ã , Im P B are invertible. This way a matrix can be identied (e.g. A) in terms of which the others can be expressed. The task is thus reduced to provide an ansatz for this matrix such that the others will have entire components. We illustrate this approach at the example of the Trott curve below. Proof. Under the change of the canonical homology basis (2.5), the vector ω = (ω 1 , . . . , ω g ) t of normalized holomorphic dierentials transforms as (4.20)

ω = 2iπ K-1 ω,
where K = 2iπA+B B. After straightforward calculations, according to the transformation law (2.6) of theta functions, it can be checked that under this change of homology basis, Fay's identity [START_REF] Fay | Theta functions on Riemann surfaces[END_REF] expressed in the Vinnikov basis transforms as

(4.21) Ẽ(a, b) Ẽ(c, d) Θ B[ δ] z + a c ω Θ B[ δ] z + d b ω + Ẽ(a, c) Ẽ(d, b) Θ B[ δ] z + a b ω Θ B[ δ] z + d c ω = Ẽ(a, d) Ẽ(c, b) Θ B[ δ](z) Θ B[ δ] z + a c ω + d b ω ,
where z = (2iπ) -1 K z, for any z ∈ C g and for any a, b, c, d ∈ R g . Here we used the fact that the transformation law (2.6)-(2.9) does not change the parity of the characteristic in theta functions. Now let us compute the vectors of the characteristic δ. Inversion of the symplectic matrix in (2.5) leads to

à B = D t -B t -C t A t A B .
Moreover, since the characteristic in Fay's identity used in [START_REF] Kalla | New degeneration of Fay's identity and its application to integrable systems[END_REF]

It would be possible to compute the theta-functional solutions in the Vinnikov basis once the symplectic transformation between this basis and the basis determined by the code is known. However, since this symplectic transformation is not unique, the found Vinnikov basis leads in general to a Riemann matrix for which the theta series converges only slowly, i.e., the value N θ in (3.1) has to be chosen very large. To avoid this problem, we compute the theta function always in the typically more convenient Tretko-Tretko basis with the characteristic of the theta functions given by (4.17)-(4.19). [START_REF] Trott | Applying Groebner Basis to Three Problems in Geometry[END_REF] given by the algebraic equation (4.24) 144 (x 4 + y 4 ) -225 (x 2 + y 2 ) + 350 x 2 y 2 + 81 = 0 is an M-curve with respect to the anti-holomorphic involution τ dened by τ (x, y) = (x, y), and is of genus 3. Moreover, this curve has real branch points only (and 28 real bitangents, namely, tangents to the curve in two places). Our computed matrices of à and B-periods read 2

Trott curve. The Trott curve

P Ã =   0.0235i 0.0138i 0.0138i 0 0.0277i 0 -0.0315 0 0.0250   , P B =   -0.0315 + 0.0235i 0.0138i -0.0250 + 0.0138i 0 -0.025 + 0.0277i 0.0250 -0.0235i 0.0138i 0.0138i   . 2
For the ease of representation we only give 4 digits here, though at least 12 digits are known for these quantities.

The Trott curve being an M-curve, the vectors of the characteristic δ satisfy (4.22) and (4.23), which leads to δ = 1 2 [ 0 0 0 1 1 0 ] t . A possible choice of a symplectic transformation bringing the computed basis to the Vinnikov basis is: The computed matrix Im P Ã Im P B -1 being (within numerical precision) equal to

A =   1 0 0 0 1 0 0 0 1   , B =   -1 0 0 0 -1 0 0 0 0   , C =   1 0 0 0 1 0 0 0 0   , D =   0 0 0 0 0 0 0 0 1   .

Note that the matrices

Im P Ã Im P B -1 =   -1 0 0 0 -1 0 0 0 0   ,
and with C, D ∈ M 3 (Z), we get from (4.26) that det A = 1. Since A ∈ M 3 (Z), the condition det A = 1 implies A ∈ Gl 3 (Z). For any A ∈ Gl 3 (Z), one can see from (4.25), (4.12) and (4.13) that B, C, D ∈ M 3 (Z), and therefore that the matrices A, B, C, D give a solution of (4.6)-(4.9). The choice A = I 3 leads to the above matrices.

The Trott curve has real bers and can thus be used to construct solutions to the 3-NLS equation via the projection map f : (x, y) → x, which is a real meromorphic function of degree 4 on the curve. We consider the points on the curve stable with respect to τ and projecting to the point with x = 0.1 in the x-sphere, and choose d = 0. The corresponding solution to the 3-NLS equation can be seen in Fig. 15.

A solution to the DS1 + equation on this curve can be constructed for points a and b stable with respect to the involution τ . The solution for a = (-0.2) (1) , b = (0.2) (2) and the choice d = 0 can be seen in Fig. 16. Note that in accordance with Remark 2.1, one would obtain a solution of DS1 -for the choice a = (-0.2) (1) and b = (0.2) (1) .

Similarly, a solution to the DS2 + equation can be obtained for points a and b subject to τ a = b. For a = (0.1 + i) (1) and b = (0.1 -i) (1) we get Fig. 17 30x 4 -61x 3 y + 41y 2 x 2 -43x 2 -11y 3 x + 42xy + y 4 -11y 2 + 9 = 0 which was studied in [START_REF] Dubrovin | Matrix nite-zone operators[END_REF] and [START_REF] Vinnikov | Self-adjoint determinantal representations of real plane curves[END_REF]. It is a genus 3 curve, dividing with respect to the antiholomorphic involution τ , without real branch point. This curve admits two real ovals. In this case the matrix H has the form After some calculations, one nds that the following matrices A, B, C, D provide a solution of (4.6)-(4.9):

H =   0 1 0 1 0 0 0 0 0   .
A =   -1 2 -1 2 -1 0 0 2 -1   , B =   1 0 1 0 1 0 1 0 0   , C =   1 -1 -1 -1 1 -1 0 0 1   , D =   0 1 1 1 0 1 0 0 -1   .
From (4.17) and (4.18) one gets for the characteristic: δ = 1 2 [ 0 0 1 1 1 0 ] t . The curve (4.27) has real bers and can thus be used to construct solutions to the focusing 3-NLS equation. We consider the points on the curve with x = 2.5 and stable with respect to τ , and we choose d = 0. The corresponding solution to the focusing 3-NLS equation can be seen in Fig. 18.

A solution to the DS1 -equation can be constructed by choosing the points a = (-4) (1) and b = (-3) (2) see Fig. which leads to the characteristic: δ = 1 2 [ 0 0 1 0 1 0 ] t . To construct a solution of the DS2 -equation on the Fermat curve, we choose the points a = (-1.5 + i) (1) and b = (-1.5 -i) (3) . The resulting solution for the choice d = 0 can be seen in Fig. 20. (1) and b = (-1.5 -i) (3) at t = 0.

Conclusion

In this paper we have presented the state of the art of the numerical evaluation of solutions to integrable equations in terms of multi-dimensional theta functions associated to real Riemann surfaces by using an approach via real algebraic curves. It was shown that real hyperelliptic curves parametrized by the branch points can be treated with machine precision for a wide range of the parameters. Even almost degenerate situations where the branch points coincide pairwise can be handled as long as at least one cut stays nite. This approach to real hyperelliptic curves [START_REF] Frauendiener | Hyperelliptic theta functions and spectral methods[END_REF][START_REF] Frauendiener | Hyperelliptic theta functions and spectral methods: KdV and KP solutions[END_REF] is being generalized to arbitrary hyperelliptic curves.

As discussed in [START_REF] Frauendiener | Algebraic curves and Riemann surfaces in Matlab[END_REF], the main diculty for general algebraic curves is the correct numerical identication of the branch points. The case of degenerations for given branch points

has not yet been studied numerically, but is planned for the future. In what concerns the solutions (2.22) to n-NLS s and similar solutions to the DS and the Kadomtsev-Petviashvili equations, the main problem in the context of real Riemann surfaces is to nd the symplectic transformation leading to the homology basis introduced in [START_REF] Vinnikov | Self-adjoint determinantal representations of real plane curves[END_REF], for which the solutions of the studied equations, with regularity conditions, can be conveniently formulated. This problem has been reduced to nd a single g × g-matrix for given periods and real ovals, the latter encoded by the matrix H. For M-curves, where the matrix H vanishes, a general formula for the characteristic (4.17)-(4.19) could be given. In the general case, an algorithm along the lines indicated in the previous section to nd the transformation will be based on a suciently general ansatz for one of the matrices entering the symplectic transformation which is the subject of future work.

Example 2 . 1 .

 21 Consider the hyperelliptic curve of genus g dened by the equation

Figure 1 .

 1 Figure 1. Homology basis on real hyperelliptic curves, contours on sheet 1 are solid, contours on sheet 2 are dashed. A-cycles are the closed contours entirely on sheet 1.

Figure

  Figure 2. Solution(2.27) to the DS1 + equation at t = 0 on a hyperelliptic curve of genus 2 with branch points -2, -1, 0, , 2, 2 + and a = (-1.9)(1) , b = (-1.1)(2) for = 1 on the left and = 10 -10 , the almost solitonic limit, on the right.

Figure 3 .

 3 Figure 3. Solution(2.27) to the DS2 + equation at t = 0 on a hyperelliptic curve of genus 2 with branch points -2, -1, 0, , 2, 2 + and a = (-1.5 + 2i)[START_REF] Anker | On the Soliton Solutions of the Davey-Stewartson Equation for Long Waves[END_REF] , b = (-1.5 -2i)(2) for = 1 on the left and = 10 -10 , the almost solitonic limit, on the right.

Figure 4 .

 4 Figure 4. Solution (2.27) to the DS1 + equation at t = 0 on a hyperelliptic curve of genus 4 with branch points -4, -3, -2, -2 + , 0, , 2, 2 + , 4, 4 +and a = (-3.9)(1) , b = (-3.1)(2) for = 1 on the left and = 10 -10 , the almost solitonic limit, on the right.

Figure

  Figure 7. Solution(2.27) to the DS1 -equation at t = 0 on a hyperelliptic curve of genus 4 with branch points -2 ± i, -1 ± i, ±i, 1 ± i, 2 ± i and a = (-4)(1) , b = (-3)(2) .

Figure 10 .

 10 Figure 10. Solution (3.2) to the 2-NLS ŝ equation with ŝ = (-1,-1) on a hyperelliptic curve of genus 2 with branch points -2, -1, 0, , 2, 2 + and a 1 = (-1.9)(1) , a 2 = (-1.1)(1) and a 3 = (-1.8)(1) for = 1 on the left and = 10 -10 , the almost solitonic limit, on the right.

2 Figure 11 .

 211 Figure 11. Solution (3.2) to the 2-NLS ŝ equation with ŝ = (1, -1) on a hyperelliptic curve of genus 2 with branch points -2, -1, 0, , 2, 2 + and a 1 = (-1.9)(1) , a 2 = (-1.8)(1) and a 3 = (-1.1)(1) for = 1 on the left and = 10 -10 , the almost solitonic limit, on the right.

2 Figure 12 .

 212 Figure 12. Solution to the 4-NLS ŝ equation on a hyperelliptic curve of genus 4 with branch points -4, -3, -2, -2 + , 0, , 2, 2 + , 4, 4 + and = 1 for a 1 = (-3.9)(1) , a 2 = (-3.7)(1) , a 3 = (-3.5)(1) , a 4 = (-3.3)(1) and a 5 = (-3.1)(1) , which leads to ŝ = (1, -1, 1, -1), on the left, and for a 1 = (-3.9)(1) , a 2 = (-3.7)(1) , a 3 = (-3.5)(1) , a 4 = (-3.1)(1) and a 5 = (-3.3)(1) , which leads to ŝ = (-1, 1, -1, -1) on the right.

2 Figure 13 .

 213 Figure 13. Solution to the 2-NLS ŝ equation with ŝ = (1, 1) on a hyperelliptic curve of genus 2 with branch points -2±i, -1±i, 1±i and a 1 = (-1.9)(1) , a 2 = (-1.8)(2) and a 3 = (-1.1)(1) .

Figure 14 .

 14 Figure 14. Solution to the 4-NLS ŝ equation with ŝ = (1, 1, 1, 1) on a hyperelliptic curve of genus 4 with branch points -2 ± i, -1 ± i, ±i, 1 ± i, 2 ± i and a 1 = (-3.9)(1) , a 2 = (-3.7)(2) , a 3 = (-3.5)(1) , a 4 = (-3.3)(2) and a 5 = (-3.1)(1) .

C

  Therefore, by (4.3) one getsA Re P Ã + B Re P B = P A (4.6)A Im P Ã + B Im P B = 0, (4.7) and by (4.4) C Re P Ã + D Re P B Im P Ã + D Im P B = Im (P B ) .(4.9)

Lemma 4 . 1 .

 41 The matrices A, B, C, D ∈ M g (Z) solving (4.6)-(4.9) satisfy:A t = Im P B C Im P Ã + D Im P B -1

  .

4. 3 .

 3 Dividing curves without real branch point. We consider the curve given by the equation(4.27)

Figure 15 .

 15 Figure 15. Solution (2.22) to the 3-NLS s equation on the Trott curve for the points with x = 0.1 on the x-sphere. The sheets are identied at the points projecting to x = -1.0129, (0.9582i, -0.9582i, 0.1146i, -0.1146i). The vector of signs equals s = (1, -1, -1) from top to bottom.

Figure 16 .

 16 Figure 16. Solution to the DS1 + equation on the Trott curve for the points a = (-0.2) (1) and b = (0.2) (2) at t = 0.The period matrices computed by the code read

Figure 17 .

 17 Figure 17. Solution to the DS2 + equation on the Trott curve for the points a = (0.1 + i)(1) and b = (0.1 -i)(1) at t = 0.

19

 19

  .

4. 4 .

 4 Fermat curve. The Fermat curves (4.28)y n + x n + 1 = 0, n > 2,n even, are real curves without real oval with respect to τ . We consider here the curve with n = 4 that has genus 3. The matrix H has the form

Figure 18 .

 18 Figure 18. Solution to the 3-NLS s equation on the dividing curve (4.27) of genus 3 for the points with x = 2.5 on the x-sphere. The sheets are identied at the ber over -2.1404+0.4404i, (-12.2492+2.0113i, -5.1634+ 1.3519i, -4.5915 + 0.9380i, -1.5405 + 0.5429i). The vector of signs is s = (1, 1, 1).

Figure 19 .

 19 Figure 19. Solution to the DS1 -equation on the dividing curve (4.27) of genus 3 for the points a = (-4)(1) and b = (-3)(2) at t = 0. and we nd

Figure 20 .

 20 Figure 20. Solution to the DS2 -equation on the Fermat curve (4.28) of genus 3 for the points a = (-1.5 + i)(1) and b = (-1.5 -i)(3) at t = 0.

  Im P Ã + C t B Im P B = 0 C t A Im P Ã + (A t D -I g ) Im P B = 0 A t C Im P Ã + A t D Im P B = Im P B ,which leads to (4.10). Equality (4.11) can be checked analogously with (4.14) and (4.16). To prove (4.12), multiply equality (4.8) from the left by the matrix A t . Using (4.14) and (4.15) one gets:

	.16)
	Multiplying equality (4.7) from the left by the matrix C t , we deduce from (4.14) and (4.15)
	that:
	C t A

  Proposition 4.1. Let ( Ã, B) be the canonical homology basis obtained with the Tretko-Tretko algorithm; we denote with a tilde the quantities expressed in this basis. Under the change of homology basis (2.5), solutions of n-NLS s and DS equations given in (2.22) and (2.27), respectively, which are expressed in the basis satisfying (2.19), transform as follows: the vector d appearing in the solutions becomes (2iπ) -1 K d where K = 2iπA + B B, and the theta function Θ = Θ B with zero characteristic, transforms to the theta function Θ B[ δ] with characteristic δ = [ δ1 , δ2 ] given by diag B t HB -2 Re P à M-1 Im P t diag A t HA -2 Re P B M-1 Im P t Im P t B Re P à -Im P t à Re P B . Moreover, the vectors N, M dened in (2.20) become A t N + C t M and B t N + D t M respectively.

	(4.17)	δ1 =	1 4	Ã	,
	(4.18)	δ2 =	1 4	B	,
	where
	(4.19)	M =

  to construct solutions (2.22), (2.27) of n-NLS s and DS is zero, from (2.9) we deduce that Im P t B Re P à -Im P t à Re P B . This completes the proof since solutions (2.22) and (2.27) are derived from Fay's identity expressed in the Vinnikov basis. Remark 4.2.In the case where the spectral curve is an M-curve, i.e. H = 0, the vectors of characteristic (4.17) and (4.18) do not depend explicitly on the symplectic matrix appearing in the change of homology basis and are uniquely dened by: diag Re P à Im P t B Re P à -Im P t à Re P B -1 Im P t

						δ1 δ2	=	1 2	Diag	D t B C t A
	(note that D t B and C t A are symmetric matrices). Substituting (4.10) and (4.11) in (4.12)
	(resp. (4.13)), it can be checked that
				C t A =	1 2	A t HA -2 Re P B M-1 Im P t B	,
				D t B =	1 2	B t HB -2 Re P Ã M-1 Im P t Ã	,
	with
	(4.22)	M = δ1 = 1 2		Ã	,
		δ2 =	1 2	diag Re P

B Im P t B Re P à -Im P t à Re P B -1 Im P t B .

  A, B, C, D are not unique since the action (2.19) of the antiholomorphic involution on the basic cycles allows for permutations of A j -cycles for instance. These matrices can be computed as follows. Since the Trott curve is an M-curve, one has H = 0. Moreover, the matrix Im(P B) being invertible here, by (4.7) one gets: (4.25) B = -A Im P Ã Im P B

	-1 .
	With (4.14) and (4.15) it follows that
	(4.26)

A t D + C Im P Ã Im P B -1 = I 3 .


Vicente Palacios

L David Mech

Problems with studying wolf predation on small prey in summer via global positioning system collars

Keywords: Canis lupus, Global positioning system (GPS) collars, Predation, Telemetry, White-tailed deer, Wolf

We attempted to study predation on various-sized prey by a male and female wolf (Canis lupus) with global positioning system (GPS) collars programmed to acquire locations every 10 min in the Superior National Forest of Minnesota. During May to August 2007, we investigated 147 clusters of locations (31% of the total) and found evidence of predation on a white-tailed deer (Odocoileus virginianus) fawn and yearling, a beaver (Castor canadensis), ruffed grouse (Bonasa umbellus), and fisher (Martes pennanti) and scavenging on a road-killed deer and other carrion. However, we missed finding many prey items and discuss the problems associated with trying to conduct such a study.

Introduction

Wolves (Canis lupus) are large carnivores that usually prey on large ungulates. However, many wolves throughout the world prey primarily on juvenile ungulates or small prey in summer due to their availability and vulnerability. Wolves in the Superior National Forest, Minnesota, prey primarily on white-tailed deer (Odocoileus virginianus) fawns during summer [START_REF] Kunkel | Wolf and bear predation on whitetailed deer fawns in northeastern Minnesota[END_REF][START_REF] Nelson | Mortality of white-tailed deer in northeastern Minnesota[END_REF][START_REF] Van Ballenberghe | Ecology of the timber wolf in northeastern Minnesota[END_REF]. Roe deer (Capreolus capreolus) fawns and wild boar (Sus scrofa) piglets constitute the main prey in some areas of the Iberian Peninsula and Italy during summer [START_REF] Barja | Prey and prey-age preference by the Iberian wolf Canis lupus signatus in a multiple-prey ecosystem[END_REF][START_REF] Capitani | A comparative analysis of wolf (Canis lupus) diet in three different Italian ecosystems[END_REF]. Ungulates less than 40 kg such as blackbucks (Antilope cervicapra) represent the main prey of wolves in India [START_REF] Jethva | Foraging ecology, economics and conservation of Indian wolves in the Bhal region of Gujarat, Western India[END_REF]. Therefore, the study of wolf predation on small prey during summer represents a key need for understanding the wolf's ecology.

Radiotelemetry has been widely used for locating wolves in winter and studying predator-prey dynamics [START_REF] Peterson | The wolf as a carnivore[END_REF]. The difficulty of locating wolves in summer (especially in forested areas) and finding prey remains or evidence of predation on small prey are possible explanations for the scarcity of summer predation studies. Recently, use of global positioning system (GPS) collars and geographic information systems (GIS) to monitor free-ranging species brought a major advance in the study of predatory behavior by wolves. Using GPS collars, [START_REF] Demma | Testing global positioning system telemetry to study wolf predation on deer fawns[END_REF] demonstrated the possibility of investigating summer wolf predation on white-tailed deer fawns, and [START_REF] Sand | Summer kill rates and predation pattern in a wolf-moose system: can we rely on winter estimates?[END_REF] studied summer wolf predation on moose (Alces alces) calves. However, most wolf predation studies using GPS collars have been conducted in winter and have been most

Communicated by P. Acevedo V. Palacios (*) Instituto Cavanilles de Biodiversidad y Biología Evolutiva, Universidad de Valencia, Apdo 22085, 46071 Valencia, Spain e-mail: v_palacios_s@hotmail.com effective for large prey [START_REF] Franke | Prediction of wolf (Canis lupus) kill-sites using hidden Markov models[END_REF][START_REF] Sand | Using GPS technology and GIS cluster analyses to estimate kill rates in wolf-ungulate ecosystems[END_REF][START_REF] Webb | Statistical methods for identifying wolf kill sites using global positioning system locations[END_REF].

High rates of data acquisition seem to be one of the keys to study predation on small prey. The amount of time that wolves spend at kills (handling time) tends to be longer for large prey. Wolves usually spend >48 h handling a moose carcass [START_REF] Ballard | Ecology of an exploited wolf population in south-central Alaska[END_REF][START_REF] Hayes | Population dynamics and prey relationships of an exploited and recovering wolf population in the southern Yukon[END_REF][START_REF] Messier | Moose (Alces alces) and wolf (Canis lupus) dynamics and the natural regulation of moose populations[END_REF][START_REF] Peterson | Wolves of the Kenai Peninsula, Alaska[END_REF]. Thus, by obtaining wolf locations hourly or longer, it is still possible to find most large-prey kill sites [START_REF] Franke | Prediction of wolf (Canis lupus) kill-sites using hidden Markov models[END_REF][START_REF] Sand | Using GPS technology and GIS cluster analyses to estimate kill rates in wolf-ungulate ecosystems[END_REF][START_REF] Webb | Statistical methods for identifying wolf kill sites using global positioning system locations[END_REF]. However, finding small prey with 1-h intervals is more problematic [START_REF] Sand | Using GPS technology and GIS cluster analyses to estimate kill rates in wolf-ungulate ecosystems[END_REF]. [START_REF] Webb | Statistical methods for identifying wolf kill sites using global positioning system locations[END_REF] found that a 30-min sampling interval was too long for detecting small prey, and [START_REF] Zimmermann | Wolf movement patterns: a key to estimation of kill rate?[END_REF] recommended using GPS position intervals of ≤ 30 min. [START_REF] Demma | Testing global positioning system telemetry to study wolf predation on deer fawns[END_REF] tried 10-to 60-min intervals and found that they held promise for studying wolf predation on deer fawns. In the present study, we used GPS collars on wolves programmed to acquire locations at 10-min intervals to (1) investigate if 10-min sampling intervals are useful for detecting small prey, (2) analyze the characteristics of clusters of locations that could be used for studying wolf predation, and (3) test the usefulness of short-location intervals for this type of study.

Study area

We conducted our study during the spring and summer of 2007 in a 130-km 2 area in the Superior National Forest of northeastern Minnesota (48°N, 92°W). The area is basically flat, being part of an old peneplain eroded by weather, water, and glaciers. Topography varies from large stretches of swamps to rocky ridges, with elevation ranging from 325 to 700 m above sea level. Temperatures range from -45°C in the winter to 37°C in the summer. Average seasonal temperatures are 6°C in the spring, 18°C in the summer, 7°C in the fall, and -11°C in the winter. Snow depths (usually from about mid-November through mid-April) generally range from 50 to 75 cm on the level. Conifers predominate in the forest overstory, with the following species present: jack pine (Pinus banksiana), white pine (Pinus strobus), red pine (Pinus resinosa), black spruce (Picea mariana), white spruce (Picea glauca), balsam fir (Abies balsamea), northern white cedar (Thuja occidentalis), and tamarack (Larix laricina). However, as a result of extensive cutting and fires, much of the conifer is interspersed with large stands of paper birch (Betula papyrifera) and aspen (Populus tremuloides). In summer, the vegetation is very dense at ground level.

Wolves occurred throughout the study area at densities of 30 to 36/1,000 km 2 during the study [START_REF] Mech | Long-term research on wolves in the Superior National Forest[END_REF]. White-tailed deer occurred at densities of 12 to 15/10 km 2 (M.H. Dexter, Minnesota Department of Natural Resources, unpublished report) and constituted the major prey of wolves in the area, primarily fawns during summer [START_REF] Kunkel | Wolf and bear predation on whitetailed deer fawns in northeastern Minnesota[END_REF][START_REF] Nelson | Mortality of white-tailed deer in northeastern Minnesota[END_REF][START_REF] Van Ballenberghe | Ecology of the timber wolf in northeastern Minnesota[END_REF].

Methods

During May 2007, we live-trapped [START_REF] Mech | Current techniques in the study of elusive wilderness carnivores[END_REF]) two wolves (a female and a male) that belonged to the Perch Lake pack, aged them [START_REF] Gipson | Accuracy and precision of estimating age of gray wolves by tooth wear[END_REF], and fitted them with Tellus GPS/VHF collars (Televilt, Lindesberg, Sweden) for 3 months (89 days), May 25 until August 22, following the same process as detailed by [START_REF] Demma | Testing global positioning system telemetry to study wolf predation on deer fawns[END_REF]. Both wolves were 1-to 2-year-old nonbreeding individuals. No evidence of pups was found in the Perch Lake pack during the period of study. We programmed GPS collars for positioning at 10-min intervals, assuming an error of 5 to 30 m of the exact location 50% and 95% of the time, respectively, as it has been estimated from landscapes similar to ours [START_REF] Dussault | Influence of satellite geometry and differential correction on GPS location accuracy[END_REF][START_REF] Moen | Accuracy of GPS telemetry collar locations with differential correction[END_REF]. We remotely downloaded location data transmitted by GPS collars (attempts every 2 days) using a VHF-receiver data logger (RX-900; Televilt International) and a handheld antenna and plotted the locations on digital topographic maps of the study area (TOPO!, National Geographic Society, Hanover, PA). (Mention of brand names does not constitute endorsement by the US government.) Downloaded data included activity levels in the x and y axes (changes in collar position during the time the collar uses to obtain the GPS fix). We assumed that the highest x and y activity levels were related to greater activity and that activity levels = 0 corresponded with the wolves resting. We defined clusters as groups of consecutive locations <100 m away from the next without directionality (minimum duration = 1 h). For small clusters (distance between the most distant locations both in the latitudinal and longitudinal axes ≤ 50 m), we calculated the cluster centroids and entered the centroid positions into handheld GPS units. For larger clusters, besides centroids, we transferred other cluster locations to assure that the overall cluster area was checked. We visited all positions transferred to handheld GPS units in search of wolf signs and prey remains in a minimum radius of 100 m per position. Clusters were visited 1 to 7 days after remote downloading. We did not visit clusters when VHF signals indicated the wolves were still present.

The GPS collars were programmed to release 90 days after initialization, according to the expected battery life. We downloaded GPS data from the released collars and plotted them in Arcview 3.2, using the Animal Movement Analysis extension [START_REF] Hooge | Animal movement extension to Arcview[END_REF], to calculate wolf movements and cluster characteristics. We added lines connecting successive GPS points to estimate wolf travel paths. To minimize the possible effects of the capture, we excluded GPS locations during the first 3 days postcapture. To avoid errors due to missing locations during the remote downloading process, we conducted all the analyses from the released collars downloads.

For each cluster, we measured nine variables related to spatiotemporal characteristics and activity levels registered in the collars (Table 1). Moreover, we created a variable (daytime) related to the time of the day at which clusters occurred. Based on daylight and local time during the study, we divided the day into four periods (central daylight, 1000-1800; central night time, 2200-0600; and transition times, 0600-1000 and 1800-2200). Depending on the percentage of time spent at each interval, we assigned a value for each cluster, assigning value 1 to strictly diurnal clusters and value 5 with strictly nocturnal clusters (Table 2). Based on our field checks, we classified the clusters into three groups: resting clusters, when only wolf beds were found; feeding clusters, when prey remains or other food sources were found inside the minimum convex polygon (MCP) of the cluster; and unknowns. We used paired t test and Mann-Whitney U test to compare resting and feeding cluster characteristics.

At a broader scale, we defined events as a group of clusters and movements related to a prey/food item. Predation events were those when we found predation evidence as hematomas in the remains or when predation was the most probable option (depending on the type of prey, remains found, activity levels, and wolves' movements related with the findings), and scavenging events were those involving food resources not related to predation. We buffered prey remains and other food source positions with radii of 100 and 500 m. We considered the beginning of the event the first wolf location inside the buffer and the end of the event the last location inside the buffer. Because wolves sometimes return to food resources, we considered an event the same when locations within a buffered area occurred on consecutive days. We measured the number of consecutive days a wolf was inside the buffered areas related to the event, the number of locations related to the event inside the buffered areas, and the overall event duration considering the 100 and 500 m buffer radii. We also measured the number of days, weeks, and location a wolf was inside the buffered areas considering the whole period of study. We used paired t test and Mann-Whitney U test to determine if there were differences in these characteristics between scavenging and predation events. To explore possible relationships between event characteristics and the estimated prey weights, we used Pearson correlation coefficients.

For statistical analyses, we used SPSS (12.0) for Windows (SPSS Inc., Chicago, IL). Statistical significance in all analyses was assumed at P≤0.05.

Results

Data from the collars resulted in 12,194 locations for the female and 12,199 locations for the male, 5% of the overall We conducted 76 remote downloads, with an average of 210 locations per download (35 h), and remotely downloaded 52% of the female and 60% of the male acquired locations. We checked 147 location clusters (31% of clusters resulting from the analysis of total collar data) and found signs of wolf activity in 88 clusters (60% of the clusters checked and 19% of the overall clusters). Clusters were checked 2.75±0.17 (SE) days after the wolf left the area.

We identified 75 resting clusters, 18 feeding clusters, and 59 unknowns (Table 3). We also considered as feeding clusters unchecked clusters related to the events when the remains/food position was inside the cluster MCP. Due to the small sample of kill sites found, we did not distinguish between predation and scavenging, considering only one class for feeding behavior.

Mean linear distances from previous points (Ldpp) were greater for feeding clusters than for resting clusters (t=-5.059, P<0.001). Wolf activity as measured by the collar (Meanact) was greater for feeding clusters (t=-4.462, P<0.001). When food remains were present, wolves were inactive (Perinact) 50% of the time, whereas in resting clusters, wolves were inactive a mean of 63% (t=-2.57, P= 0.01). Feeding clusters were more nocturnal than resting clusters (Mann-Whitney U=326.0, P<0.001). Differences in the other measured variables were not significant (P> 0.1) except that duration of resting clusters tended to be longer (t=1,777, P=0.08).

We defined 12 feeding events, five of them predation and the other seven scavenging. During the 3-month study, we found remains of one white-tailed deer fawn and a beaver (Castor canadensis) related to the male wolf and a white-tailed deer yearling, a ruffed grouse (Bonasa umbellus), and a fisher (Martes pennanti) related to the female wolf, which we considered as kill sites (Figs. 1a, b; 2a, b, andc). Based on the number of close beds found in the clusters (one single wolf bed found per cluster versus multiple beds and flattened vegetation), we assumed that with the fawn, grouse, and fisher kills, the wolf was alone, whereas the female was with pack mates at the yearling deer kill, and the male was with pack mates at least at the end of the beaver kill event. We found a correlation between estimated prey weight and the number of consecutive days the wolf used the 500-m-buffered areas (Pearson correlation coefficient r=0.96, n=5, P<0.001) and the overall event duration assuming a 500-m-buffer radius (Pearson correlation coefficient r=0.95, n=5, P=0.01).

The scavenging events included a road-killed adult deer ("road kill"), a pile of unidentified rotting meat, probably a bear-hunting (Ursus americanus) bait station, two bear-bait stations, and a pile of fish in a garbage dump. Both wolves visited the road kill (Fig. 3). The female and male were together at the road kill for 12 h. Event characteristics were similar for both the female and the male, although the male remained much longer inside the buffered area during successive visits (Table 4). The other food item visited by both wolves was the rotting meat (Fig. 4). During this event, the female and the male were together for 3 h. The use of this food source was completely different for both wolves based on the variables measured (Table 4) and the movements and clusters of the event (Fig. 4).

During one bear-bait station event, the female was alone, and during the dead-fish event, the male was alone. The male was with other pack mates during the other bait station event.

There were no significant differences between scavenging and predation event characteristics in the variables measured, although the small sample for this analysis needs to be considered for properly interpreting this result.

Discussion

Programming collars to obtain locations at 10-min intervals provides information unobtainable a decade ago. However, we still found many problems in trying to study wolf predation on small prey in summer via GPS collars. The first problem is related to the specific equipment used and the field work needed for obtaining appropriate sample sizes. By terrestrial, remote downloading, we downloaded 50% to 60% of possible locations. The VHF-receiver data logger capacity was <2 days of information. Every download made after 2 days of the previous one resulted in real-time missing data. In other studies, collars were programmed for positioning at 30-min intervals [START_REF] Sand | Summer kill rates and predation pattern in a wolf-moose system: can we rely on winter estimates?[END_REF], 15-min to 6-h intervals [START_REF] Webb | Statistical methods for identifying wolf kill sites using global positioning system locations[END_REF], and 60-min intervals [START_REF] Franke | Prediction of wolf (Canis lupus) kill-sites using hidden Markov models[END_REF], and remote-data downloading was conducted weekly or biweekly [START_REF] Sand | Using GPS technology and GIS cluster analyses to estimate kill rates in wolf-ungulate ecosystems[END_REF][START_REF] Zimmermann | Wolf movement patterns: a key to estimation of kill rate?[END_REF]). [START_REF] Demma | Testing global positioning system telemetry to study wolf predation on deer fawns[END_REF] programmed one GPS collar for positioning at 10-min intervals during 1800 to 1000 h. They tried to remotely download the GPS collar data on 25 days, being 68% successful. In our study, we tried to locate the wolves almost every 2 days during the 89-day study; even so, we missed 40% to 50% of locations due to the fact that download trials were not always successful. Furthermore, the downloading in such circumstances (high fix rates and limited data logger capacity) directly reduced the time invested in checking clusters and, as a result, a smaller number of clusters checked.

We also found other problems that would be unrelated to the type of GPS collar or technology used. We classified 40% of the checked clusters as unknowns because we did not find wolf sign or predation evidence. In spite of being difficult to find small-prey remains itself, summer conditions and habitat characteristics, often swampy and with dense thicket, could affect the percentage of unknown clusters. A possible solution to increase the percentage of knowns in future studies could be the use of well-trained dogs for checking clusters. Dogs have been used for carnivore scat detection [START_REF] Mackay | Scat detection dogs[END_REF]) and for checking clusters in wolf kill-rate studies as well [START_REF] Sand | Summer kill rates and predation pattern in a wolf-moose system: can we rely on winter estimates?[END_REF].

We classified known clusters as resting and feeding clusters. Although it might be simplistic to reduce wolf behavior at clusters to resting and feeding, these results could serve to decide in future studies which clusters to check, avoiding diurnal clusters with low linear distance between consecutive locations, a high percentage of inactive time, and low activity levels if the purpose of the study is to estimate kill rates. We did not further evaluate Fig. 1 GPS-collared male wolf feeding events in the Superior National Forest during springsummer 2007. The figures represent all the MCPs of GPS location clusters and movements since the wolf arrived in the 500-m buffered area until the end of the event. a Deer fawn predation. b Beaver predation. c Dead-fish feeding. d Bear-bait station feeding. c and d events were chronologically consecutive. Based on our definition of an event, buffering fish and baiting station with 500-m radius, both events overlap in some extent. It implies that the collared wolf visited both food resources alternatively. Findings: (1) prey remains/food resources; (2) wolf beds; (3) wolf scats the accuracy of the activity levels or try to determine which behaviors could be related with different activity levels, but this type of study could be very useful in the future for understanding wild wolves' behavior.

In this study, we found that wolves fed on carrion and bear-bait stations, and killed two deer, a beaver, a fisher, and a grouse. As it has been reported, wolves can kill large prey and small prey and eat carrion and garbage, taking advantage of a variety of resources [START_REF] Peterson | The wolf as a carnivore[END_REF]. Four of the five preys we found were small and medium. Other studies reported that wolves preyed on grouse and beavers [START_REF] Zimmermann | Wolf movement patterns: a key to estimation of kill rate?[END_REF][START_REF] Sand | Summer kill rates and predation pattern in a wolf-moose system: can we rely on winter estimates?[END_REF], similar to our results. We know of no other reports of wolves killing fishers. Although only minor scavenging was reported in the Canadian Rocky Mountains [START_REF] Webb | Statistical methods for identifying wolf kill sites using global positioning system locations[END_REF]), our wolves spent considerable time scavenging. Thus, attempts to estimate kill rates in such a multiprey/ food resources system become difficult.

Other problems in determining summer kill rates include the fact that wolves of the same pack do not travel together all the time during summer. Our wolves belonged to the same pack, and their home ranges overlapped 89% (unpublished data), but they travelled separately almost throughout the study, similar to other packs near our study area [START_REF] Demma | Wolf use of summer territory in northeastern Minnesota[END_REF][START_REF] Demma | Testing global positioning system telemetry to study wolf predation on deer fawns[END_REF]). [START_REF] Demma | Wolf use of summer territory in northeastern Minnesota[END_REF] suggested that this pattern could be due to the fact that while they seek food independently, pack mates visit each other's kill remains. In our study, the most time the wolves were together was at the road kill and rotting meat. One wolf arrived at the area first, and the other joined it. Travelling separately but remaining loosely associated could increase the chance of finding as many different food sources as possible for all the pack mates. On the other hand, our results suggest that predation on small prey tends to be related to wolves travelling alone. Hence, pack member independence can be a problem for estimating kill rates when few animals per pack are radiocollared [START_REF] Webb | Statistical methods for identifying wolf kill sites using global positioning system locations[END_REF]. Moreover, as we found at the road kill and rotting meat, for the same event, different wolves can behave differently, making it more complex to build predictive models of predation patterns in such circumstances.

Despite our small sample, we found a strong correlation between the time a wolf is near a kill and the prey's weight, similar to [START_REF] Webb | Statistical methods for identifying wolf kill sites using global positioning system locations[END_REF]. However, [START_REF] Zimmermann | Wolf movement patterns: a key to estimation of kill rate?[END_REF] reported great variability in handling times, due to pack member independence and stochastic factors such as human disturbance and others. Spending more time around large prey should be normal especially for wolves travelling alone and with no disturbance.

A final problem to consider is the fact that during summer, most wolf packs have pups. We did not find evidence of pups in our Perch Lake pack during the study. The pups' feeding needs should affect handling time at kill sites, as wolves commonly carry food to the pups [START_REF] Packard | Wolf behavior: reproductive, social, and intelligent[END_REF]. Highlighted weeks reflect the event dates.

a

Based on both a 100-and 500-m buffer, the female spent 44 h near the road kill, and the male spent 27 h.

b

Based on a 100-m buffer, the female spent 17 h, and the male 109 h near the meat; based on a 500-m buffer, the female spent 24 h and the male 161 h. (3) wolf scats Thus, during summer, the movement patterns depending on the prey species, presence of pups, and the pack role of the radioed wolf could vary considerably. We conclude that it is difficult to get an accurate account of the number and type of small prey via GPS collars even with the high location acquisition rates we used. Despite obtaining positions at 10-min intervals, we found the same problems as others studying wolf predation, especially on small prey such as deer fawns. Although [START_REF] Demma | Testing global positioning system telemetry to study wolf predation on deer fawns[END_REF] located four deer fawn kill sites checking clusters of one GPS-collared wolf for a 37-day-period, they reported the need for improving the detection of kills as well. There are multiple factors (prey size, prey availability, presence of other food sources, presence of pups, pack role, stochastic factors, human disturbances, whether the study animal is alone or with other pack mates, among others) that increase the variability of predation patterns, making it difficult to estimate accurate kill rates and modeling predation patterns for wolves.

However, our results could be useful for future studies of wolf predation on deer fawns. Improving the datadownloading process, increasing the chance of finding wolf signs and prey remains, and rejecting resting clusters could be helpful to increase the number of small-prey kill sites found. This type of study still provides valuable information for the understanding of wolf predation, although extensive field work must still accompany it.

Fig. 2

 2 Fig. 2 GPS-collared female wolf feeding events in the Superior National Forest during spring-summer 2007. The figures represent all the MCPs of GPS location clusters and movements since the wolf arrived in the 500-m buffered area until the end of the event. a Yearling deer predation. b Grouse predation. c Fisher predation. d Bear-bait station feeding. Findings: (1) prey remains/ food resources; (2) wolf beds; (3) wolf scats; (4) digging

Fig. 4

 4 Fig.4GPS-collared male wolf (left) and female wolf rotting meat feeding events in the Superior National Forest during spring-summer 2007. The figures represent all the MCPs of GPS location clusters and movements since the wolf arrived in the 500-m buffered area until the end of the event. The male arrived 31 July and left 6 August. The

Table 2

 2 Times of day used to define daytime values for GPS location clusters for wolves in the Superior National Forest during springsummer 2007

	Table 1 Variables measured for each GPS location cluster for	Variable	Description
	wolves in the Superior National Forest during spring-summer 2007	Duration Time between the first and last location in the cluster Prevdist Minimum distance travelled from the previous cluster (adding linear distances between consecutive locations)
		Postdist	Minimum distance travelled to the next cluster (adding linear distances between
			consecutive locations)
		Prevtime Time since the last location of the previous cluster until the first location of the current cluster
		Posttime	Time since the last location of the current cluster until the first location of the next cluster
		Area	Cluster minimum-convex-polygon area
		Ldpp	Mean value of the linear distance from previous point
		Perinact	Percent in a cluster of locations when both activity levels (x and y axis) = 0
		Meanact	Mean value of activity levels (x+y) when x and y are distinct from 0
			Daytime	Central daylight	Transition times	Central
			value	(1000-1800)	(1800-2200 and	night time
					0600-1000)	(2200-0600)
				% of cluster	% of cluster	% of cluster
				duration	duration	duration
			1	100	0	0
			2	50 to <100	<50	<50
			3	<50	0-100	<50
			4	<50	<50	50 to <100
			5	0	0	100

Table 4

 4 Characteristics of visits to food items by both collared wolves in the Superior National Forest during spring-summer 2007 assuming a variable 500-m buffer

	Event	Wolf	Week of Study

Eur J Wildl Res (2011) 57:149-156

Acknowledgments This study was supported by the Biological Resources Division, United States Geological Survey, and the United Sates Department of Agriculture NorthCentral Research Station. We thank M.E. Nelson, D. Thompson, and A. Bennett for field assistance, and especially M. Clark for special devotion to this project, and D.J. Demma and S.M. Barber-Meyer for reviewing the manuscript and offering helpful suggestions for improvement.