Computing Theta Functions with Julia

TL;DR

The paper introduces Theta.jl, a Julia package for accurate and efficient numerical computation of Riemann theta functions and their derivatives, optimized for repeated evaluations in small dimensions.

Contribution

It provides a new Julia package that enables precise and fast computation of theta functions with characteristics, including derivatives, tailored for small-dimensional problems.

Findings

Efficient implementation of theta function evaluations in Julia.

Application to the Schottky problem in genus five.

Optimized for multiple evaluations with the same Riemann matrix.

Abstract

We present a new package Theta.jl for computing with the Riemann theta function. It is implemented in Julia and offers accurate numerical evaluation of theta functions with characteristics and their derivatives of arbitrary order. Our package is optimized for multiple evaluations of theta functions for the same Riemann matrix, in small dimensions. As an application, we report on experimental approaches to the Schottky problem in genus five.

Equations

$$\theta :{\mathbb{C}}^g\times {\mathbb{H}}_g\to \mathbb{C},\theta (z,\tau )=\sum _{n\in {\mathbb{Z}}^g}\mathbf{e}\left(\frac{1}{2}{n}^t\tau n+n^{t}z\right)$$

$$\theta [m](z,\tau )=\theta \left[\begin{array}{c}\epsilon \\ \delta \end{array}\right](z,\tau )=\sum _{n\in {\mathbb{Z}}^g}\mathbf{e}\left(\frac{1}{2}{\left(n+\frac{\epsilon }{2}\right)}^t\tau \left(n+\frac{\epsilon }{2}\right)+{\left(n+\frac{\epsilon }{2}\right)}^t\left(z+\frac{\delta }{2}\right)\right)$$

$$\theta [m](\tau )=\theta [m](0,\tau ).$$

$$D\left(k^{(1)},\dots ,k^{(N)}\right)\theta (z,\tau )=\sum _{i_1,\dots ,i_N=1}^g{k}_{i_1}^{(1)}\cdots k_{i_N}^{(N)}\frac{{\partial }^N\theta (z,\tau )}{\partial z_{i_1}\cdots \partial z_{i_N}}.$$

$$ϵ=(2\pi {)}^N\frac{g}{2}{\left(\frac{2}{\rho }\right)}^g\sum _{j=0}^N(\genfrac{}{}{0ex}{}{N}{j})\frac{1}{{\pi }^{j/2}}\Vert T^{-1}{\Vert }^j{\sqrt{g}}^{N-j}\Gamma \left(\frac{g+j}{2},{\left(R-\frac{\rho }{2}\right)}^2\right).$$

$$\begin{array}{rl}& e^{\pi y^t{Y}^{-1}y}(2\pi i{)}^N\sum _{n\in C_R}\left(k^{(1)}\cdot \left(n-\eta \right)\right)\cdots \left(k^{(N)}\cdot \left(n-\eta \right)\right)\\ & \times \mathbf{e}\left(\frac{1}{2}{\left(n-\eta \right)}^{t}X\left(n-\eta \right)+{\left(n-\eta \right)}^t\left(x+\frac{\delta }{2}\right)\right)e^{-\Vert v(n+\frac{\epsilon }{2}){\Vert }^2},\end{array}$$

$$C_R=\{n\in {\mathbb{Z}}^g|\pi (n-c{)}^{t}Y(n-c)<R^2,|c_j|<1,\forall j=1,\dots ,g\}.$$

$$S_R=\{n\in {\mathbb{Z}}^g|\Vert v(n)\Vert <R\}.$$

$$\displaystyle S_{R,\varepsilon}=\left\{n\in\mathbb{Z}^{g}\,\bigg{|}\,\left\|v\left(n+\frac{\varepsilon}{2}\right)\right\|<R\right\}\,.$$

$$U_R=\{n\in {\mathbb{Z}}^g|\pi (n-c{)}^{t}Y(n-c)<R^2,|c_j|<1/2,\forall j=1,\dots ,g\},$$

$$D(k^{(1)},\dots ,k^{(N)})\theta (z,\tau )=(2\pi i{)}^N\sum _{n\in {\mathbb{Z}}^g}(k^{(1)}\cdot n)\cdots (k^{(N)}\cdot n)\mathbf{e}\left(\frac{1}{2}{n}^t\tau n+n^{t}z\right)$$

$$ϵ$$

$$\displaystyle\quad\quad\times\mathbf{e}\left(\frac{1}{2}(n-[Y^{-1}y])^{t}X(n-[Y^{-1}y])+(n-[Y^{-1}y])^{t}x\right)e^{-\|v(n)\|^{2}}\bigg{|}$$

$$ϵ$$

$$\le (2\pi {)}^N\frac{g}{2}{\left(\frac{2}{\rho }\right)}^g\sum _{j=0}^N(\genfrac{}{}{0ex}{}{N}{j})\frac{1}{{\pi }^{j/2}}\Vert T^{-1}{\Vert }^j\Vert Y^{-1}y{\Vert }^{N-j}\Gamma \left(\frac{g+j}{2},{\left(R-\frac{\rho }{2}\right)}^2\right).$$

$$x^6{y}^2-4x^4{y}^2-2x^3{y}^3-2x^{4}y+2x^{3}y+4x^2{y}^2+3x{y}^3+y^4+4x^{2}y+2x{y}^2+x^2-4xy-2y^2-2x+1.$$

$${\Theta }_{\tau }=\left\{z\in A_{\tau }|\theta (z,\tau )=0\right\}.$$

$$J(C)=H^0(C,{\omega }_C{)}^{\vee }/H_1(C,\mathbb{Z}),$$

$$H^0(C,{\omega }_C)\times H^1(C,\mathbb{Z})⟶\mathbb{C},(\omega ,\alpha )\mapsto {\int }_{\alpha }\omega .$$

$$z=\frac{\epsilon }{2}+\tau \frac{\delta }{2},\text{ for }m=\left[\begin{array}{c}\epsilon \\ \delta \end{array}\right]\in (\mathbb{Z}/2\mathbb{Z}{)}^{2g}.$$

$$\theta \left(\frac{\epsilon }{2}+\tau \frac{\delta }{2},\tau \right)=0\text{ if and only if }\theta \left[\begin{array}{c}\epsilon \\ \delta \end{array}\right](0,\tau )=0.$$

$${\mathcal{J}}_g\cap {\theta }_{\text{null}}\subseteq {\theta }_{\text{null}}^3$$

$$m=\left[\begin{array}{ccccc}1& 0& 0& 1& 0\\ 1& 0& 1& 1& 0\end{array}\right]$$

$$\left(\begin{array}{ccccc}-2.79665+5.29764i& -9.57825-9.04671i& 7.36305+2.28697i& 7.58338+5.34729i& 6.15667-1.90199i\\ -9.57825-9.04671i& 18.9738+8.34582i& -23.1027-3.10545i& -9.31944-0.822821i& 0.524289-3.64991i\\ 7.36305+2.28697i& -23.1027-3.10545i& 16.8441-1.15986i& 13.9363-4.56541i& -3.32248+4.10698i\\ 7.58338+5.34729i& -9.31944-0.822821i& 13.9363-4.56541i& 2.89309+1.21773i& 3.86617-0.546202i\\ 6.15667-1.90199i& 0.524289-3.64991i& -3.32248+4.10698i& 3.86617-0.546202i& -12.9726-1.928i\end{array}\right)$$

$$47.946229109152995$$

$$-15.491689246713147$$

$$-9.512858919129267$$

$$-2.7271385943272036\times 10^{-15}$$

$$-5.698014266322794\times 10^{-15}$$

Full text

1 Introduction

The Riemann theta function is the holomorphic function

[TABLE]

where $\mathbf{e}(x)=e^{2\pi ix}$ and $\mathbb{H}_{g}$ is the Siegel upper-half space, which consists of all complex symmetric $g\times g$ matrices with positive definite imaginary part. Theta functions occupy a central role throughout mathematics, appearing in fields as diverse as algebraic geometry [BL04, Mum07], number theory [Mum07, EZ85], differential geometry [AcSS20], integrable systems [KS13, Seg08], discrete mathematics [RSD17], cryptography [Gau07] and statistics [AA19].

We present a new package Theta.jl for numerical computations of theta functions, programmed in Julia [BEKS17]. Our package is specialized for multiple evaluations of theta functions for the same Riemann matrix $\tau\in\mathbb{H}_{g}$ and different $z$, for small values of the genus $g$. Our implementation is based on the algorithm from [DHB*+*04], which we extend to support computations of theta functions with characteristics and derivatives of arbitrary order. Our package is designed as an alternative to existing packages such as algcurves [DHB*+*04] in Maple, abelfunctions [SD16] in Sage and [FJK19] in Matlab, with additional functionalities and optimizations.

As an application, we study numerical approaches to the Schottky problem in genus five. The Schottky problem asks to recognize Jacobians of curves amongst principally polarized abelian varieties, and is one of the central questions in algebraic geometry since the 19th century [Gru12]. The first nontrivial case of the Schottky problem is in genus four, which is completely solved [Igu81]. For a recent approach linking computations and tropical geometry see [CKS19]. In this paper, we describe computational approaches for studying the Schottky problem in genus five, using our new package. In particular, we use Theta.jl to compute the equations in [FGS, Acc83] which give a weak solution to the Schottky problem in genus five. We also use our package for computations on the genus five Schottky problem for Jacobians with a vanishing theta null, which is described in our companion paper [AC].

Acknowledgments: We are grateful to Bernd Sturmfels for suggesting to study the Schottky problem in genus five, and for his continuous encouragement. We thank Paul Breiding, Gavril Farkas, Jörg Frauendiener, Sam Grushevsky, Christian Klein, Riccardo Salvati Manni, Andrey Soldatenkov, Sasha Timme and Sandro Verra for useful comments and discussions. We thank the anonymous referees for their suggestions, which improved the quality of the paper and of the code. This project was initated at the Max Planck Institute for Mathematics in the Sciences in Leipzig, which both authors would like to thank for the hospitality and support at various stages of this work.

2 Theta functions

We recall here the basic definitions about theta functions with characteristics. For a more detailed account we refer to [BL04],[Mum07],[Igu81]. A characteristic is an element $m\in(\mathbb{Z}/{2}\mathbb{Z})^{2g}$, which we represent as a vector $m=\begin{bmatrix}\varepsilon\\ \delta\end{bmatrix}$ where $\varepsilon,\delta\in\{0,1\}^{g}$. The Riemann theta function with characteristic $m$ is defined as

[TABLE]

and it is a holomorphic function $\theta[m]\colon\mathbb{C}^{g}\times\mathbb{H}_{g}\to\mathbb{C}$. The Riemann theta function in (1) is a special case of (2), where the characteristic is the all-zero vector. The sign of a characteristic $m$ is defined as $e(m)=(-1)^{\varepsilon^{t}\delta}$, and we call a characteristic even or odd if the sign is 1 or $-1$ respectively. As a function of $z$, $\theta[m](z,\tau)$ is even (respectively odd) if and only if the characteristic $m$ is even (respectively odd). There are $2^{g-1}(2^{g}+1)$ even theta characteristics and $2^{g-1}(2^{g}-1)$ odd theta characteristics.

The theta constants are the functions on $\mathbb{H}_{g}$ obtained by evaluating the theta functions with characteristics at $z=0$,

[TABLE]

Theta constants corresponding to odd characteristics vanish identically.

3 Numerically approximating theta functions

We describe in this section the algorithm that we use to compute theta functions in Theta.jl. In our implementation, we modify the algorithm from [DHB*+*04], generalizing it for theta functions with characteristics and derivatives of arbitrary order.

In this section, we separate $z\in\mathbb{C}^{g}$ and $\tau\in\mathbb{H}_{g}$ into real and imaginary parts, by writing $z=x+iy$, $\tau=X+iY$, where $x,y\in\mathbb{R}^{g}$ and $X,Y$ are real symmetric $g\times g$ matrices. We also denote by $Y=T^{t}T$ the Cholesky decomposition of $Y$, where $T$ is upper-triangular. For any real vector $V\in\mathbb{R}^{g}$, we use $[V]$ to denote the vector whose entries are the entries of $V$ rounded to the closest integers, and we denote $[\![V]\!]=V-[V]$.

We set $v(n)=\sqrt{\pi}T(n+[\![Y^{-1}y]\!])$ and we define the lattice $\Lambda=\{v(n)\,|\,n\in\mathbb{Z}^{g}\}$, letting $\rho$ be the length of the shortest nonzero vector in $\Lambda$. We denote by $\Gamma(z,x)=\int_{x}^{\infty}t^{z-1}e^{-t}dt$ the incomplete Gamma function.

3.1 Theta functions with characteristics

In [DHB*+*04], Deconinck et al. derive numerical approximations of the theta function and its first and second derivatives. We extend their results for computing theta functions with characteristics and derivatives of arbitrary order.

We denote the $N$-th order derivative of the theta function along the vectors $k^{(1)},\ldots,k^{(N)}$ as

[TABLE]

By the quasi-periodicity of the theta function, it suffices to consider inputs $z$ of the form $z=a+\tau b$, for $a,b\in[0,1)^{g}$.

Theorem 3.1**.**

Fix $\tau\in\mathbb{H}_{g}$, $\epsilon>0$. Let $k^{(1)},\ldots,k^{(N)}\in\mathbb{C}^{g}$ be unit vectors, and let $R$ be the greater of $\frac{1}{2}\sqrt{g+2N+\sqrt{g^{2}+8N}}+\frac{\rho}{2}$ and the real positive solution of $R$ in

[TABLE]

For $z$ of the form $z=a+\tau b$, for $a,b\in[0,1)^{g}$, and ${\footnotesize\begin{bmatrix}\varepsilon\\ \delta\end{bmatrix}}\in\{0,1\}^{2g}$, the $N$-th derivative $D(k^{(1)},\ldots,k^{(N)})\theta{\footnotesize\begin{bmatrix}\varepsilon\\ \delta\end{bmatrix}}(z,\tau)$ of the theta function with characteristic is approximated by

[TABLE]

with an absolute error $\epsilon$ on the product of $(2\pi i)^{N}$ with the sum, where $\eta=[Y^{-1}y]-\frac{\varepsilon}{2}$ and

[TABLE]

Proof of Theorem 3.1:   We first consider the case $N=0$ without derivatives. Then the result for characteristics $\varepsilon=\delta=0$ is proven in [DHB*+*04, Theorem 2], where they replace the deformed ellipsoid $C_{R}$ in (7) with the ellipsoid

[TABLE]

For arbitrary characteristics $\varepsilon,\delta$, we see from (2) that we can compute the corresponding theta function in a similar way as the usual theta function, by translating $z$ to $z+\frac{\delta}{2}$, and translating the lattice points in the sum from $n$ to $n+\frac{\varepsilon}{2}$. Note that this only changes the real part of $z$, while the imaginary part stays the same. Hence the approximation in Theorem 3.1 holds for theta functions with characteristics, if we take the sum over the ellipsoid

[TABLE]

To obtain a uniform approximation for any $z\in\mathbb{C}^{g}$ and any characteristic, we take the union of the ellipsoids $S_{R,\varepsilon}$ from (9) as $z$ and $\varepsilon$ vary. Since $v(n+\frac{\varepsilon}{2})=\sqrt{\pi}T(n+[\![Y^{-1}y]\!]+\frac{\varepsilon}{2})$, and the entries of $[\![Y^{-1}y]\!]+\frac{\varepsilon}{2}$ have absolute value at most 1, it follows that the deformed ellipsoid $C_{R}$ from (7) is the union of the ellipsoids $S_{R,\varepsilon}$.

To prove the result in the case of derivatives of order $N$, it will be enough to prove it for the case of zero characteristics, and then follow the same strategy as above. More precisely, we are going to prove the same statement as in 3.1, where $\varepsilon=\delta=0$ and $C_{R}$ is replaced by

[TABLE]

To do so, we write the derivative as

[TABLE]

and then the error in the approximation is

[TABLE]

Since the $k^{(i)}$ have norm one, using the triangle inequality and the Cauchy-Schwartz inequality we can bound this by

[TABLE]

Using again the triangle inequality and the binomial expansion, we get to the bound

[TABLE]

We then apply [DHB*+*04, Lemma 2] to get the bound

[TABLE]

For inputs $z$ of the form $z=a+\tau b$, we can write $z$ as $z=a+(X+iY)b=(a+Xb)+iYb=x+iy$. Then $\|Y^{-1}y\|=\|b\|\leq\sqrt{g}$. Substituting this into the expression for $\epsilon$, the result follows. $\Box$

Remark 3.2**.**

The $R$ appearing in 3.1 is computed numerically.

4 Computing theta functions in Julia

4.1 Interface

Our Julia package Theta.jl is available at the following website, which has instructions and a link to more detailed documentation.

https://github.com/chualynn/Theta.jl

We describe the basic interface of the package here. Starting with a matrix $\tau\in\mathbb{H}_{g}$, we first construct a RiemannMatrix from it. This is a type in Theta.jl which contains information needed to compute the theta function with input $\tau$. As an example, we start with a genus 5 curve defined by the singular model

[TABLE]

We compute the Riemann matrix $\tau$ of the curve using the package [BSZ19] in Sage [The19], and we type it as an input in Julia. We then construct a RiemannMatrix in Theta.jl, where we specify in the input the options to compute a Siegel transformation, an error of $10^{-12}$, and to compute derivatives up to the fourth order.

julia> R = RiemannMatrix($\tau$, siegel=true, $\epsilon$=1.0e-12, nderivs=4);

We pick some input $z$ and compute the theta function $\theta(z,\tau)$ as follows.

julia> z = [1.041+0.996im; 1.254+0.669im; 0.591+0.509im; -0.301+0.599im; 0.388+0.051im]; julia> theta(z, R) -854877.6514446283 + 2.3935081163150463e6im

We can compute derivatives of theta functions by specifying the directions using the optional argument derivs. For instance, to compute $\frac{\partial^{3}\theta}{\partial z_{3}\partial z_{4}}(z,\tau)$, we use

julia> theta(z, R, derivs=[[0,0,1,0,0], [0,0,0,1,0]]) 1.0478325534969474e8 - 3.369999441122761e8im

We can also compute derivatives of theta functions with characteristics, where we specify the characteristic using the optional argument char.

julia> theta(z, R, derivs=[[1,0,0,0,0]], char=[[0,1,0,0,1],[1,1,0,0,1]]) -2.448093122926732e7 + 3.582557740667034e7im

4.2 Algorithms

We describe here some details of the algorithms and the design choices that we made in our implementation.

4.2.1 Choice of ellipsoid

We optimize our package for multiple evaluations of theta functions at the same Riemann matrix $\tau$, and with different inputs $z$, characteristics and derivatives. We do this using the approximation in Theorem 3.1, which allows us to compute derivatives of theta functions with characteristics, for inputs $z$ of the form $z=a+\tau b$, for $a,b\in[0,1)^{g}$. In this approximation, we take the sum over the deformed ellipsoid $C_{R}$ of (7), which depends only on the order $N$ of the derivative for a fixed $\tau$. Hence it suffices to compute the ellipsoids $C_{R}$ once for each order of the derivative that we are interested in, after which we can compute theta functions for any $N$-th order derivatives. These ellipsoids are stored in the RiemannMatrix type.

4.2.2 Lattice reductions

In [DHB*+*04], the authors approximate the length $\rho$ of the shortest vector of the lattice generated by $T$ using the LLL algorithm by Lenstra, Lenstra and Lovász [LLL82]. This is a reasonable choice if $g$ is large, since the LLL algorithm gives a polynomial time approximation, but with an error that grows exponentially with $g$. In our implementation, since we focus on lattices with small dimensions, we compute the shortest vector exactly using the enumeration algorithm in [SE94]. Moreover, by computing $\rho$ exactly, we obtain a smaller ellipsoid (7) than if we use the LLL algorithm.

If we are interested in computing the theta function for a fixed $\tau$ at many values of $z$, it may be more efficient if we transform $\tau$ such that the ellipsoids in (7) contain fewer lattice points. For this purpose, we use Siegel’s algorithm, which iteratively finds a new matrix where the corresponding ellipsoid has a smaller eccentricity. In our implementation, we compute the Siegel transformation once for each Riemann matrix, and work with the Siegel-transformed matrix for all computations. We use the algorithm for Siegel reduction described in [DHB*+*04, FJK19], where we use the algorithm for HKZ reduction in [ZQW12] as a subroutine.

4.3 Comparisons with other packages

The main advantage of Theta.jl over other packages [DHB*+*04, FJK19, SD16] is that we support computations of theta functions with characteristics, as well as their derivatives, which to our knowledge is not implemented elsewhere. Moreover, we make optimizations described in Section 4.2 for faster computations with a fixed Riemann matrix of low genus.

We compare the performance of Theta.jl with the Sage package abelfunctions [SD16], by comparing the average time taken to compute the genus 5 FGSM relations of Section 5.1, as well as to compute the Hessian matrix of Section 5.3. For our experiments, we sample matrices in the Siegel upper-half space as follows. First we sample $5\times 5$ matrices $M_{X},M_{Y}$ such that the entries are random floating point numbers between $-1$ and $1$, using the random number generators in Julia and NumPy. Then we sample $\tau\in\mathbb{H}_{5}$ as $\tau=\frac{1}{2}(M_{X}+M_{X}^{t})+M_{Y}^{t}M_{Y}i$. This is implemented in Theta.jl for general dimensions $g$, in the function random_siegel(g). In each experiment, we randomly sample $1000$ such matrices, then we compute the FGSM relations and the Hessian matrix using both packages on a standard laptop. We list in the table below the average time and standard deviation.

[TABLE]

One major reason for the faster runtime on Theta.jl is the use of the Siegel transformation on the Riemann matrix, which is not implemented in abelfunctions. This also leads to the higher standard deviation in the computations for the latter.

5 Applications to the Schottky problem in genus five

Here we describe the main application that we had in mind when designing our package: experiments around the Schottky problem in genus five. We start with a brief account of the background of the problem, referring to [Gru12] for more details.

An abelian variety is a projective variety that has the structure of an algebraic group, and it is a fundamental object in algebraic geometry. Especially important are principally polarized abelian varieties, which can all be described in terms of Riemann matrices. For every $\tau\in\mathbb{H}_{g}$, we define the corresponding principally polarized abelian variety (ppav) as the quotient $A_{\tau}=\mathbb{C}^{g}/\Lambda_{\tau}$, where $\Lambda_{\tau}=\mathbb{Z}^{g}\oplus\tau\mathbb{Z}^{g}$ is a sublattice of $\mathbb{C}^{g}$. The polarization on $A_{\tau}$ is given by the theta divisor

[TABLE]

Two ppavs $A_{\tau}$ and $A_{\tau^{\prime}}$ are isomorphic if and only if the corresponding Riemann matrices are related via an action of the symplectic group $\Gamma_{g}=\operatorname{Sp}(2g,\mathbb{Z})$. Hence, the quotient $\mathcal{A}_{g}=\mathbb{H}_{g}/\operatorname{Sp}(2g,\mathbb{Z})$ is the moduli space of principally polarized abelian varieties of dimension $g$. This is a quasi-projective variety of dimension $\dim\mathcal{A}_{g}=\dim\mathbb{H}_{g}=\frac{g(g+1)}{2}$, and we can look at the theta constants $\theta[m](0,\tau)$ as homogeneous coordinates on (a finite cover of) $\mathcal{A}_{g}$.

Perhaps the most important example of abelian varieties are Jacobians of Riemann surfaces. For a Riemann surface $C$ of genus $g$, its Jacobian is defined as the quotient

[TABLE]

where the lattice $H^{1}(C,\mathbb{Z})$ is embedded in $H^{0}(C,\omega_{C})^{\vee}$ via the integration pairing

[TABLE]

The Jacobian is a principally polarized abelian variety, and the corresponding Riemann matrix $\tau\in\mathcal{A}_{g}$ can be obtained by computing bases of $H^{0}(C,\omega_{C})$ and $H^{1}(C,\mathbb{Z})$, as well as the integration pairing. This is implemented numerically in the packages abelfunctions [SD16] and RiemannSurfaces [BSZ19] in Sage, and algcurves [DHB*+*04] in Maple.

The Schottky locus $\mathcal{J}_{g}$ is the closure of the set of Jacobian varieties in $\mathcal{A}_{g}$, and the Schottky problem asks for a characterization of $\mathcal{J}_{g}$ inside $\mathcal{A}_{g}$. It is one of the most celebrated questions in algebraic geometry, dating from the 19th century. There are many possible interpretations of and solutions to the Schottky problem. Here we focus on the most classical one, which asks for equations in the theta constants $\theta[m](0,\tau)$ that vanish exactly on the Schottky locus. In this form, the Schottky problem is completely solved only in genus $4$, with an explicit equation given by Schottky and Igusa [Igu81]. A computational implementation and analysis of this solution was presented in [CKS19].

The weak Schottky problem asks for explicit equations that characterize Jacobians up to extra irreducible components. A solution to this problem was given in genus $5$ by Accola [Acc83], and in a recent breakthrough, by Farkas, Grushevsky and Salvati Manni in all genera [FGS]. In the rest of this section, we discuss briefly these two solutions, together with related algorithms that we implemented in Theta.jl. We also present a computational solution of a weak Schottky problem for genus five Jacobians with a theta null, from our companion paper [AC].

5.1 Farkas, Grushevsky and Salvati Manni’s solution

In a recent preprint [FGS], H. Farkas, Grushevsky and Salvati Manni give a solution to the weak Schottky problem in arbitrary genus. More precisely, for every genus $g\geq 4$ they give $\binom{g-2}{2}=\frac{(g-2)(g-3)}{2}$ explicit homogeneous equations of degree $2^{3\cdot 2^{g-4}+1}$ in the theta constants, such that their zero locus contains the Schottky locus as an irreducible component.

In the case of genus $5$, this gives $3$ equations of degree $128$. We implement them in the function fgsm() in Theta.jl. Using the same example matrix $\tau$ from Section 4.1, the function fgsm($\tau$) gives us the output $7.850462293418876\mbox{e-}16$. This is expected since $\tau$ is the Jacobian of a genus 5 curve.

5.2 Accola’s equations in genus 5

A solution to the weak Schottky problem in genus $5$ was given already by Accola [Acc83] in 1983, in the form of eight equations of degree $32$ in the theta constants whose zero locus contains the Schottky locus as an irreducible component. We implement these equations in the function accola() in Theta.jl. Again using the example $\tau$ from Section 4.1, the function accola($\tau$) gives us the output $3.062334813867916\mbox{e-}9$, which is expected since $\tau$ is in the Schottky locus.

5.3 Schottky problem for Jacobians with a vanishing theta null

We describe here a variant of the Schottky problem focusing on two-torsion points on Jacobians, referring to our companion article [AC] for a more complete account. A two-torsion point on an abelian variety $A_{\tau}$ is a point $z\in A_{\tau}$ such that $2z=0$. These can be written as

[TABLE]

Hence two-torsion points correspond to characteristics, and we say that such a point is even or odd if the corresponding characteristic is. Observe that

[TABLE]

Thus the two-torsion points in $\Theta_{\tau}$ correspond to the characteristics $m$ such that the theta constants $\theta[m](0,\tau)$ vanish. For this reason, we say that $A_{\tau}$ has a vanishing theta null if it has an even two-torsion point in the theta divisor. The abelian varieties with this property have been intensely studied and they form a divisor $\theta_{\text{null}}$ in $\mathcal{A}_{g}$. The Jacobians with a vanishing theta null lie in the locus $\mathcal{J}_{g}\cap\theta_{\text{null}}$ and they correspond to Riemann surfaces with an effective even theta characteristic. The Schottky problem in this case becomes that of recognizing $\mathcal{J}_{g}\cap\theta_{\text{null}}$ inside $\theta_{\text{null}}$.

The first observation is that a vanishing theta null is automatically a singular point of the theta divisor, because the partial derivatives $\frac{\partial\theta[m]}{\partial z_{i}}$ are odd. Hence one is led to study the local structure of $\Theta_{\tau}$ around the singular point, and the first natural invariant is the rank of the quadric tangent cone, which corresponds to the rank of the Hessian matrix of $\theta$ at the theta null. In particular, if a Jacobian has a vanishing theta null, then the quadric tangent cone has rank at most three. Hence

[TABLE]

where we denote by $\theta^{3}_{\text{null}}$ the locus of abelian varieties with a vanishing theta null whose quadric tangent cone has rank at most three. Conversely, Grushevsky and Salvati Manni proved in [GM08] that this inclusion is actually an equality in genus $4$, confirming a conjecture of H. Farkas. In the same paper, they ask whether $\mathcal{J}_{g}\cap\theta_{\text{null}}$ is an irreducible component of $\theta^{3}_{\text{null}}$ in higher genera, which would imply a solution to the weak Schottky problem for Jacobians with a vanishing theta null. The main result of our companion paper [AC] is an affirmative answer in genus $5$.

Theorem 5.1**.**

[AC]** In genus five, the locus $\mathcal{J}_{5}\cap\theta_{\text{null}}$ is an irreducible component of $\theta^{3}_{\text{null}}$.

We observe that the containment $\tau\in\theta^{3}_{\text{null}}$ can be checked explicitly. Indeed, the condition of having an even two-torsion point in the theta divisor can be checked by evaluating the finitely many theta constants $\theta[m](0,\tau)$, and then numerically computing the rank of the Hessian matrix. We present such a computation here, which is also in [AC]. From the example in Section 4.1, we use the function schottky_null($\tau$) in Theta.jl. The output gives the even characteristic

[TABLE]

where the theta constant vanishes. The output also gives the corresponding Hessian matrix

[TABLE]

The Hessian has the eigenvalues

[TABLE]

so it has numerical rank 3 as expected.