Jean Kieffer

Sign choices in the AGM for genus two theta constants

Keywords: Theta functions, Genus 2, Algorithms, Borchardt mean Subject classification (2010): 11Y35, 11Y16, 11F41, 11F27

Existing algorithms to compute genus 2 theta constants in quasi-linear time use Borchardt sequences, an analogue of the arithmetic-geometric mean for four complex numbers. In this paper, we show that these Borchardt sequences are only given by good choices of square roots, as in the genus 1 case. This removes the sign indeterminacies when computing genus 2 theta constants without relying on numerical integration.

Résumé

Les algorithmes existants pour le calcul de thêta-constantes en genre 2 en temps quasilinéaire utilisent des suites de Borchardt, un analogue de la moyenne arithmético-géométrique pour quatre nombres complexes. Dans cet article, nous montrons que ces suites de Borchardt sont constituées uniquement de bons choix de signes, comme c'est le cas en genre 1. Ce résultat permet de lever les indéterminations de signes lors du calcul de thêta-constantes en genre 2 sans recours à l'intégration numérique.

Introduction

Denote by H g the Siegel half space of principally polarized abelian varieties of dimension g, consisting of all matrices τ ∈ M g (C) such that τ is symmetric and Im(τ ) is positive definite; for instance, H 1 is the usual upper half plane. The theta constants are the holomorphic functions on H g defined by

θ a,b (τ ) = m∈Z g exp iπ m + a 2 t τ m + a 2 + m + a 2 t b , (1)
where a and b run through {0, 1} g (by convention, vectors in formula (1) are written vertically). Theta constants have a fundamental importance in the theory of Siegel modular forms, as every scalar-valued Siegel modular function of any weight on H g has an expression in terms of quotients of theta constants [15, Thm. 9 p. 222]. Moreover, for 1 ≤ g ≤ 3, then the stronger result that every Siegel modular form is a polynomial in the theta constants holds [START_REF] Igusa | On the graded ring of theta-constants[END_REF][START_REF] Igusa | On the graded ring of theta-constants (II)[END_REF][START_REF] Freitag | On the variety associated to the ring of theta constants in genus 3[END_REF]. In numerical algorithms manipulating modular forms, the following operations are therefore very common: first, given (quotients of) theta constants at a given τ ∈ H g , compute τ ; second, given τ ∈ H g , compute the theta constants θ a,b (τ ). For instance, these operations are important building blocks in algorithms computing modular polynomials [START_REF] Enge | Computing modular polynomials in quasi-linear time[END_REF][START_REF] Milio | A quasi-linear time algorithm for computing modular polynomials in dimension 2[END_REF][START_REF] Milio | Modular polynomials on Hilbert surfaces[END_REF] or Hilbert class polynomials [START_REF] Enge | The complexity of class polynomial computation via floating point approximations[END_REF][START_REF] Enge | Computing class polynomials for abelian surfaces[END_REF][START_REF] Streng | Computing Igusa class polynomials[END_REF] via complex approximations.

The arithmetic-geometric mean (AGM) [START_REF] Borchardt | Theorie des arithmetisch-geometrisches Mittels aus vier Elementen[END_REF][START_REF] Cox | The arithmetic-geometric mean of Gauss[END_REF][START_REF] Bost | Moyenne arithmético-géométrique et périodes de courbes de genre 1 et 2[END_REF][START_REF] Jarvis | Higher genus arithmetic-geometric means[END_REF] gives an algorithm to find τ given its theta constants. This algorithm is quasi-linear in terms of the required precision. In order to compute theta constants in quasi-linear time as well, a well-studied strategy is to combine the AGM with Newton iterations. This strategy was first described in [START_REF] Dupont | Fast evaluation of modular functions using Newton iterations and the AGM[END_REF] in the genus 1 case, in [START_REF] Dupont | Moyenne arithmético-géométrique, suites de Borchardt et applications[END_REF] in the genus 2 case, and later extended to theta functions, in opposition to theta constants, in [START_REF] Labrande | Computing Jacobi's θ in quasi-linear time[END_REF][START_REF] Labrande | Computing theta functions in quasi-linear time in genus 2 and above[END_REF]. These references also outline extensions to higher genus.

The genus 1 case. Let us detail the genus 1 case to convey the general idea. After reducing the argument τ ∈ H 1 using Gauss's algorithm [25, §6.1], we can assume that τ belongs to the classical fundamental domain under the action of SL 2 (Z), denoted by F 1 .

First assume that theta quotients at τ ∈ F 1 are given. Then the sequence B(τ ) = θ 2 0,0 (2 n τ ) θ 2 0,0 (τ ) , θ 2 0,1 (2 n τ ) θ 2 0,0 (τ ) n≥0 is an AGM sequence, meaning that each term is obtained from the previous one by means of the transformation

(x, y) → x + y 2 , √ x √ y
for some choice of the square roots. This is a consequence of the duplication formula [23, p. 221], the correct square roots being the theta quotients themselves. In the algorithm, the sign ambiguity is easily removed using the fact that √ x and √ y should lie in a common open quarter plane [5, Thm. 2]: we say that the sequence B(τ ) is given by good sign choices. It converges quadratically to 1/θ 2 0,0 (τ ), as the series expansion (1) shows. It turns out that the sequence B(-1/τ ) is also an AGM sequence with good sign choices [START_REF] Dupont | Fast evaluation of modular functions using Newton iterations and the AGM[END_REF]Prop. 7]. Its first term can be computed from theta quotients at τ using the transformation formulas for theta constants under SL 2 (Z). The limit of B(-1/τ ) is 1/θ 2 0,0 (-1/τ ). Finally, we can recover τ using the formula

θ 2 0,0 -1 τ = -iτ θ 2 0,0 (τ ). (2)
Since AGM sequences with good sign choices converge quadratically, this gives an algorithm to invert theta functions on F 1 with quasi-linear complexity in the output precision, at least for fixed τ . This method was already known to Gauss [10, X.1, pp. 184-206], and we recommend [3, §3C] for a historical exposition of Gauss's works on the AGM and elliptic functions.

In order to compute theta functions at a given τ ∈ F 1 , the most efficient known method is to build a Newton scheme [START_REF] Dupont | Fast evaluation of modular functions using Newton iterations and the AGM[END_REF], using the AGM method to invert theta constants. This yields a quasi-linear algorithm to compute genus 1 theta constants, whose complexity can be made uniform in τ ∈ F 1 [START_REF] Dupont | Fast evaluation of modular functions using Newton iterations and the AGM[END_REF]Thm. 5].

The genus 2 case. A similar strategy can be applied to theta functions in genus 2, using Borchardt sequences, a generalization of AGM sequences for four complex numbers [START_REF] Borchardt | Theorie des arithmetisch-geometrisches Mittels aus vier Elementen[END_REF][START_REF] Bost | Moyenne arithmético-géométrique et périodes de courbes de genre 1 et 2[END_REF][START_REF] Jarvis | Higher genus arithmetic-geometric means[END_REF]. Let us refer to §2 for the definition of Borchardt sequences, the numbering of genus 2 theta constants, and the definition of the matrices γ k ∈ Sp 4 (Z) for 0 ≤ k ≤ 3. The Borchardt sequences we consider are the sequences B(γ k τ ) for 0 ≤ k ≤ 3, where

B(τ ) = θ 2 0 (2 n τ ) θ 2 0 (τ ) , θ 2 1 (2 n τ ) θ 2 0 (τ ) , θ 2 2 (2 n τ ) θ 2 0 (τ ) , θ 2 3 (2 n τ ) θ 2 0 (τ ) n≥0 for every τ ∈ H 2 .
Their first terms are given by different combinations of theta quotients at τ (see Corollary 3.3). It is known that for a given τ , all but a finite number of sign choices in these Borchardt sequences are good, and the other sign choices can be determined using certified computations of hyperelliptic integrals at relatively low precision: see the discussion before Prop. 3.3 in [START_REF] Labrande | Computing theta functions in quasi-linear time in genus 2 and above[END_REF], and [START_REF] Molin | Computing period matrices and the Abel-Jacobi map of superelliptic curves[END_REF] for an algorithm that provides this input. However, the required precision and the cost of the numerical integration algorithms depend heavily on τ .

Actually, when τ belongs to the usual fundamental domain F 2 under the action of Sp 4 (Z), practical experiments suggest that all sign choices are good in the genus 2 algorithm as well [4, Conj. 9.1], [START_REF] Enge | Computing class polynomials for abelian surfaces[END_REF]Conj. 9]. The goal of this paper is to prove this fact. More precisely, we define in §2 a subset F ′ ⊂ H 2 containing F 2 , and prove the following result.

Theorem 1.1. For every τ ∈ F ′ , every 0 ≤ k ≤ 3 and every n ≥ 0, the theta constants θ j (2 n γ k τ ) for 0 ≤ j ≤ 3 are contained in a common open quarter plane.

Dupont [4, Prop. 9.1] proved this result in the particular case of γ 0 = I 4 .

As a consequence, we can invert genus 2 theta constants in quasi-linear time by using only Borchardt sequences with good sign choices. On the practical side, this result reduces the effort needed to invert genus 2 theta constants with controlled precision losses; see for instance [4, §7.4.2] for an analysis of precision losses when computing limits of Borchardt sequences. On the theoretical side, we hope that our result can be a first step towards removing other heuristic assumptions when computing genus 2 theta constants (in particular, the assumption [4, §10.2] that the function used in the Newton scheme is analytic with invertible Jacobian), and obtaining algorithms with uniform complexity in τ ∈ F 2 .

This document is organized as follows. In Section 2, we introduce our notational conventions. In Section 3, we use the action of the symplectic group to bring the matrices 2 n γ k τ ∈ H 2 closer to the cusp at infinity: this is critical to obtain accurate information from the series expansion [START_REF] Borchardt | Theorie des arithmetisch-geometrisches Mittels aus vier Elementen[END_REF]. We give estimates on genus 2 theta constants in Section 4, and we finish the proof of the main theorem in Section 5.

Acknowledgement. The author would like to thank Aurel Page and the anonymous referees for their careful reading and helpful suggestions to improve the exposition.

Theta constants and Borchardt sequences

We define a Borchardt sequence to be a sequence of complex numbers (s

(n) b ) b∈(Z/2Z) 2 , n≥0
with the following property: for every n ≥ 0, there exist t

(n) b for b ∈ (Z/2Z) 2 such that t (n) b is a square root of s (n) b , and
s (n+1) b = 1 4 b1+b2=b t (n) b1 t (n) b2
for each b ∈ (Z/2Z) 2 .

The duplication formula [23, p. 221] states that for every τ ∈ H 2 , the sequence

B(τ ) = θ 2 0,b (2 n τ ) b∈{0,1} 2 ,n≥0
is a Borchardt sequence; the choice of square roots at each step is given by the theta constants θ 0,b (2 n τ ) themselves. By the series expansion (1), we have

θ 0,b (2 n τ ) = m∈Z 2 exp -2 n πm t Im(τ )m exp iπ 2 n m t Re(τ )m + m t b .
When n tends to infinity, all the terms except m = 0 converge rapidly to zero, because Im(τ ) is positive definite. Therefore the Borchardt sequence B(τ ) converges to (1, 1, 1, 1). We say that a set of complex numbers is in good position when it is included in an open quarter plane seen from the origin, i.e. a set of the form {r exp(i(α 0 + α)) | r > 0 and 0 < α < π/2} for some α 0 ∈ R. The property of being in good position is invariant by nonzero complex scaling. A Borchardt sequence is given by good sign choices if for every n ≥ 0, the complex numbers t Let us now detail the algorithm to recover τ ∈ H 2 from its theta quotients. We first introduce the matrices γ k ∈ Sp 4 (Z) alluded to in the introduction. Let

S 1 = 1 0 0 0 , S 2 = 0 0 0 1 , S 3 = 0 1 1 0 ,
and define the matrix

γ k ∈ Sp 4 (Z) for 0 ≤ k ≤ 3 by γ 0 = I 4 , and γ k = -I 2 -S k S k -I + S 2 k for 1 ≤ k ≤ 3.
For convenience, we also introduce a numbering of theta constants [4, §6.2]:

θ (a0,a1),(b0,b1) =: θ j where j = b 0 + 2b 1 + 4a 0 + 8a 1 ∈ 0, 15 .
Assuming that the choices of square roots in the sequences B(γ k τ ) can be determined, we can compute τ ∈ F 2 from its theta quotients as follows.

Algorithm 2.1 ([4, §9.2.3]).

Input: The projective vector of squares of theta constants θ 2 j (τ ) for j ∈ 0, 15 , for some τ ∈ H 2 . Output: The matrix τ .

θ 2 0 (γ 1 τ ) = -iz 1 θ 2 4 (τ ), θ 2 0 (γ 2 τ ) = -iz 2 θ 2 8 (τ ), θ 2 0 (γ 3 τ ) = -det(τ )θ 2 0 (τ ).
In the sequel, we use the following notational conventions. For τ ∈ H 2 , we write

τ = z 1 (τ ) z 3 (τ ) z 3 (τ ) z 2 (τ ) and x j (τ ) = Re z j (τ ) y j (τ ) = Im z j (τ ) for 1 ≤ j ≤ 3.
For 1 ≤ j ≤ 3, we also write q j (τ ) = exp(-πy j (τ )).

We denote by λ 1 (τ ) the smallest eigenvalue of Im(τ ), and define

r(τ ) = min λ 1 (τ ), y 1 (τ ) 2 , y 2 (τ ) 2 .
We often omit the argument τ to ease notation. We define F ′ to be the set of all τ ∈ H 2 such that the following conditions are satisfied:

|x j (τ )| ≤ 1 2 for each 1 ≤ j ≤ 3, 2 |y 3 (τ )| ≤ y 1 (τ ) ≤ y 2 (τ ), y 1 (τ ) ≥ √ 3 2 , |z j (τ )| ≥ 1 for j ∈ {1, 2}. (3)
The domain F ′ contains the classical fundamental domain F 2 for the action of Sp 4 (Z) on H 2 [START_REF] Klingen | Introductory lectures on Siegel modular forms[END_REF]Prop. 3 p. 33]. Assumptions similar to (3) are usual when giving analytic estimates on theta constants: for instance, the domain B in [START_REF] Streng | Computing Igusa class polynomials[END_REF] is defined by the first three inequalities of (3). Finally, for each τ ∈ H 2 , we write

ξ 4,6 (τ ) = 2 exp iπ z 1 (τ ) 4 , ξ 8,9 (τ ) = 2 exp iπ z 2 (τ ) 4 , ξ 0 (τ ) = 1 + 2 exp(iπz 1 (τ )) + 2 exp(iπz 2 (τ )), ξ 0,2 (τ ) = 1 + 2 exp(iπz 1 (τ )), ξ 0,1 (τ ) = 1 + 2 exp(iπz 2 (τ )), and
ξ 12 (τ ) = exp iπ z 1 (τ ) + z 2 (τ ) 4 exp iπ z 3 (τ ) 2 + exp -iπ z 3 (τ ) 2 . (4)
These complex numbers correspond to the first term(s) of the series defining theta constants at τ . For instance, ξ 4,6 (τ ) approximates both θ 4 (τ ) and θ 6 (τ ).

We will recall the definitions (4) before using them in the computations of §4.

3 Other expressions for theta constants at 2 n γ k τ

For every n ≥ 0, we define

η (n) 1 =     0 0 -1 0 0 1 0 0 1 0 2 n 0 0 0 0 1     , η (n) 2 =     1 0 0 0 0 0 0 -1 0 0 1 0 0 1 0 2 n     , η (n) 3 =     0 0 0 -1 0 0 -1 0 0 1 2 n 0 1 0 0 2 n     , and η (n) 4 =     0 0 1 0 0 1 0 0 -1 0 0 0 0 0 0 1     η (n) 3 . Lemma 3.1. Let n ≥ 0.
1. For every 1 ≤ k ≤ 4, the matrix η

(n) k belongs to Sp 4 (Z).
2. For every τ =

z 1 z 3 z 3 z 2 ∈ H 2 , we have τ (n) 1 := η (n) 1 (2 n γ 1 τ ) = 2 -n z 1 z 3 z 3 2 n z 2 , τ (n) 2 := η (n) 2 (2 n γ 2 τ ) = 2 n z 1 z 3 z 3 2 -n z 2 , τ (n) 3 := η (n) 3 (2 n γ 3 τ ) = 2 -n τ, and
τ (n) 4 := η (n) 4 (2 n γ 3 τ ) = -2 n /z 1 -z 3 /z 1 -z 3 /z 1 2 -n (z 2 -z 2 3 /z 1 )
.

(5)

Proof.

1. The lines of each η

(n) k define a symplectic basis of Z 4 .
2. The action of Sp 4 (Z) on H 2 extends to an action of the larger group

GSp 4 (Q) = γ ∈ GL 4 (Q) | ∃µ ∈ Q × , γ t 0 I 2 -I 2 0 γ = µ 0 I 2 -I 2 0 . The matrix 2 n γ k τ is the image of τ under -2 n I 2 -2 n S k S k -I + S 2 i ∈ GSp 4 (Q).
When we multiply this matrix by η

(n) k
on the left, we obtain

Diag(-1, -2 n , -2 n , -1) for k = 1, Diag(-2 n , -1, -1, -2 n ) for k = 2, and
Diag(-1, -1, -2 n , -2 n ) for k = 3.
We recall the transformation formulas for theta constants in genus 2. For a square matrix m, we denote by m 0 the column vector containing the diagonal of m.

γ = A B C D ∈ Sp 4 (Z). Define α β = γ t a -(CD t ) 0 b -(AB t ) 0 .
Then, for every τ ∈ H 2 , we have

θ a,b (γτ ) = κ(γ) ζ ε(γ,a,b) 8 det(Cτ + D) 1/2 θ a ′ ,b ′ (τ )
where

ζ 8 = e iπ/4 , a ′ b ′ = α β mod 2, ε(γ, a, b) = 2(Bα) t (Cβ) -(Dα) t (Bα) -(Cβ) t (Aβ) + 2((AB t ) 0 ) t (Dα -Cβ),
and κ(γ) is an eighth root of unity depending only on γ, with a sign ambiguity coming from the choice of a holomorphic square root of det(Cτ + D).

Corollary 3.3. For every τ ∈ H 2 , we have the following equalities of projective tuples:

(θ j (2 n γ 1 τ )) 0≤j≤3 = (θ 4 (τ ) : θ 0 (τ ) : θ 6 (τ ) : θ 2 (τ )) if n = 0, (θ 0 (τ (n) 1 ) : θ 4 (τ (n) 1 ) : θ 2 (τ (n) 1 ) : θ 6 (τ (n) 1 )) if n ≥ 1, (θ j (2 n γ 2 τ )) 0≤j≤3 = (θ 8 (τ ) : θ 9 (τ ) : θ 0 (τ ) : θ 1 (τ )) if n = 0, (θ 0 (τ (n) 2 ) : θ 1 (τ (n) 2 ) : θ 8 (τ (n) 2 ) : θ 9 (τ (n) 2 )) if n ≥ 1, (θ j (2 n γ 3 τ )) 0≤j≤3 = (θ 0 (τ (n) 3 ) : θ 8 (τ (n) 3 ) : θ 4 (τ (n) 3 ) : θ 12 (τ (n) 3 )) for every n ≥ 0, (θ j (2 n γ 3 τ )) 0≤j≤3 = (θ 0 (τ (n) 4 ) : θ 8 (τ (n) 4 ) : θ 1 (τ (n) 4 ) : θ 9 (τ (n) 4 ))
for every n ≥ 0,

where the τ

(n) j
are defined as in [START_REF] Dupont | Fast evaluation of modular functions using Newton iterations and the AGM[END_REF].

Proof. Apply Proposition 3.2 to the matrices η

(n) i .
When τ ∈ F ′ , the real and imaginary parts of τ

(n) k
for 1 ≤ k ≤ 3 are easy to study: for instance, from the second inequality in (3) we always have

y 3 (τ (n) k ) 2 ≤ 1 4 y 1 (τ (n) k )y 2 (τ (n) k ).
Such estimates are less obvious for the matrices τ

(n) 4 . Lemma 3.4. Let τ ∈ F ′ . Then, for every n ≥ 0, we have

y 3 (τ (n) 4 ) ≤ 3 2 n+2 y 1 (τ (n) 4 ), y 3 (τ (n) 4 ) 2 ≤ 3 7 y 1 (τ (n) 4 )y 2 (τ (n) 4
), and

x 2 (τ (n) 4 ) ≤ 9 2 n+3 .
Proof. Write z 1 for z 1 (τ ), etc. We have

y 3 (τ (n) 4 ) = Im(-z 3 /z 1 ) = 1 |z 1 | 2 (x 3 y 1 -y 3 x 1 ), so y 3 (τ (n) 4 ) ≤ 3y 1 4 |z 1 | 2 = 3 2 n+2 y 1 (τ (n) 4 ), since y 1 (τ (n) 4 ) = 2 n y 1 / |z 1 |
2 by (5). For the second inequality, we have

Im(τ (n) 4 ) = 2 -n z 1 -2 -n z 3 0 1 -t (2 -n Im τ ) 2 -n z 1 -2 -n z 3 0 1 -1 so det Im(τ (n) 4 ) = 1 |z 1 | 2 det Im τ.
Moreover det Im τ ≥ 3 4 y 2 1 , so

y 3 (τ (n) 4 ) 2 y 1 (τ (n) 4 )y 2 (τ (n) 4 ) ≤ y 3 (τ (n) 4 ) 2 y 3 (τ (n) 4 ) 2 + 3y 2 1 4|z1| 2 ≤ 1 1 + 4 3 |z 1 | 2 ≤ 3 7 .
For the last inequality, we compute

2 n x 2 (τ (n) 4 ) = x 2 - 1 |z 1 | 2 ((x 2 3 -y 2 3 )x 1 + 2x 3 y 3 y 1 )
and

1 |z 1 | 2 (x 2 3 -y 2 3 )x 1 ≤ 1 2 max{x 2 3 , y 2 3 |z 1 | 2 } ≤ 1 8 , so 2 n x 2 (τ (n) 4 ) ≤ 1 2 + 1 8 + 1 2 = 9 8 .

Bounds on theta constants

Typically, when τ ∈ H 2 is close enough to the cusp at infinity (more precisely when Im z 1 (τ ), Im z 2 (τ ), and det Im(τ ) are large), useful information on theta constants at τ can be obtained from the series expansion (1). Our computations are similar in spirit to those found in [17, pp. 116

f (k + 2) -f (k + 1) ≥ f (k + 1) -f (k) for every k ≥ 0. Let 0 < q < 1. Then ∞ k=0 q f (k) ≤ q f (0) 1 -q f (1)-f (0) . Proof. Use that f (k) ≥ f (0) + k(f (1) -f (0)) for all k. Lemma 4.2. Let k ≥ 1, and let τ ∈ H 2 such that y 3 (τ ) 2 ≤ 1 4 y 1 (τ )y 2 (τ ) and k |y 3 (τ )| ≤ y 2 (τ ). Define ξ 4,6 (τ ) = 2 exp iπ z 1 (τ ) 4 and ρ (k) 4,6 (q 1 , q 2 ) = q 2 1 1 -q 4 1 + q 1-1 k 2 1 -q 3-1 k 2 + q 1+ 1 k 2 1 -q 3+ 1 k 2 + q 7/8 1 q 1/2 2 (1 -q 3/2 2 )(1 -q 2 1 ) + q 25/8 1 q 3/2 2 (1 -q 9/2 2 )(1 -q 6 1 )
.

Then for j ∈ {4, 6}, we have

θ j (τ ) ξ 4,6 (τ ) -1 ≤ ρ (k)
4,6 q 1 (τ ), q 2 (τ ) .

Proof. Write u = 1/2 0 . Using the definition, we obtain

θ j (τ ) ξ 4,6 (τ ) -1 ≤ 1 2 q -1/4 1 m∈Z 2 m =( 0 0 ), -1 0 exp -π(m + u) t Im(τ )(m + u) .
We split this sum in two parts, according to whether the second coordinate of m is zero or not. The first part gives

q -1/4 1 m∈N+ 3 2 q m 2 1 ≤ q -1/4 1 q 9/4 1 1 -q 4 1 = q 2 1 1 -q 4
1 by Lemma 4.1. The second part is

q -1/4 1 m1∈N+ 1 2 m2≥1 q m 2 1 1 q m 2 2 2 • 2 cosh(2πy 3 m 1 m 2 ).
We use the fact that for every (m 1 , m 2 ) ∈ R 2 + ,

2y 3 m 1 m 2 ≤ y 1 2 m 2 1 + y 2 2 m 2 2 .
When m 1 = 1/2, we use the following bound instead:

2y 3 m 1 m 2 = |y 3 m 2 | ≤ y 2 m 2 k .
Therefore the total contribution of the second part is bounded by

q -1/4 1 m2≥1 q 1/4 1 q m 2 2 2 • 2 cosh π y 2 k m 2 + q -1/4 1 m1∈N+ 3 2 m2≥1 q m 2 1 1 q m 2 2 2 • 2 cosh π y 1 2 m 2 1 + y 2 2 m 2 2 ≤ q 1-1 k 2 1 -q 3-1 k 2 + q 1+ 1 k 2 1 -q 3+ 1 k 2 + q 7/8 1 q 1/2 2 (1 -q 3/2 2 )(1 -q 2 1 ) + q 25/8 1 q 3/2 2 (1 -q 9/ 2
2 )(1 -q 6 1 ) by other applications of Lemma 4.1. Define

ξ 8,9 (τ ) = 2 exp iπ z 2 (τ ) 4 , and
ρ (k) 8,9 (q 1 , q 2 ) = q 2 2 1 -q 4 2 + q 1-1 k 1 1 -q 3-1 k 1 + q 1+ 1 k 1 1 -q 3+ 1 k 1 + q 7/8 2 q 1/2 1 (1 -q 3/2 1 )(1 -q 2 2 ) + q 25/8 2 q 3/2 1 (1 -q 9/2 1 )(1 -q 6 2 )
.

Then for j ∈ {8, 9}, we have

θ j (τ ) ξ 8,9 (τ ) -1 ≤ ρ (k)
8,9 q 1 (τ ), q 2 (τ ) .

Proof. We proceed in a similar fashion as in the proof of Lemma 4.2 by switching the roles of q 1 and q 2 . Lemma 4.4. Let τ ∈ H 2 such that

y 3 (τ ) 2 ≤ 1 4 y 1 (τ )y 2 (τ ). Define ξ 0 (τ ) = 1 + 2 exp(iπz 1 (τ )) + 2 exp(iπz 2 (τ )), ξ 0,2 (τ ) = 1 + 2 exp(iπz 1 (τ )), ξ 0,1 (τ ) = 1 + 2 exp(iπz 2 (τ )), and
ρ 0 (q 1 , q 2 ) = 2q 4 1 1 -q 5 1 + 2q 4 2 1 -q 5 2 + 2q 1/2 1 q 1/2 2 (1 -q 3/2 1 )(1 -q 3/2 2 ) + 2q 3/2 1 q 3/2 2 (1 -q 9/2 1 )(1 -q 9/2 2 )
.

Then we have

|θ 0 (τ ) -ξ 0 (τ )| ≤ ρ 0 (q 1 (τ ), q 2 (τ )), |θ j (τ ) -ξ 0,2 (τ )| ≤ ρ 0 (q 1 (τ ), q 2 (τ )) + 2q 2 (τ ) for j ∈ {0, 2}, |θ j (τ ) -ξ 0,1 (τ )| ≤ ρ 0 (q 1 (τ ), q 2 (τ )) + 2q 1 (τ )
for j ∈ {0, 1}, and |θ j (τ ) -1| ≤ ρ 0 (q 1 (τ ), q 2 (τ )) + 2q 1 (τ ) + 2q 2 (τ ) for 0 ≤ j ≤ 3.

Proof. We proceed again in a similar fashion as in the proof of Lemma 4.2. The terms of ρ 0 (q 1 , q 2 ) are obtained by considering the following subsets of indices m ∈ Z 2 :

{( m1 0 ) | |m 1 | ≥ 2}, { 0 m2 | |m 2 | ≥ 2}, and
{( m1 m2 ) | |m 1 | ≥ 1, |m 2 | ≥ 1}. Lemma 4.5. Let τ ∈ H 2 such that |x 3 (τ )| ≤ 1 2 and 2 |y 3 (τ )| ≤ min{y 1 (τ ), y 2 (τ )}.
Write

ξ 12 (τ ) = exp iπ z 1 (τ ) + z 2 (τ ) 4 exp iπ z 3 (τ ) 2 + exp -iπ z 3 (τ ) 2 ,
and

ρ 12 (q 1 , q 2 ) = q 3/2 1 1 -q 7/2 1 + q 5/2 1 1 -q 9/2 1 + q 3/2 2 1 -q 7/2 2 + q 5/2 2 1 -q 9/2 2 + q 7/8 1 q 7/8 2 (1 -q 2 1 )(1 -q 2 2 ) + q 25/8 1 q 25/8 2
(1 -q 6 1 )(1 -q 6 2 )

.

Then we have θ 12 (τ ) ξ 12 (τ ) -1 ≤ ρ 12 q 1 (τ ), q 2 (τ ) .

Proof. By (1), we have

θ 12 (τ ) = 2 m1∈N+ 1 2 m2∈N+ 1 2 exp iπ(m 2 1 z 1 + m 2 2 z 2 ) • exp(2πim 1 m 2 z 3 ) + exp(-2πim 1 m 2 z 3 ) .
We leave the term corresponding to (m 1 , m 2 ) = ( 1 2 , 1 2 ) aside, and write

θ 12 (τ ) 2 exp(iπ(z 1 + z 2 )/4) -(exp(iπz 3 /2) + exp(-iπz 3 /2)) ≤ (m1,m2)∈(N+ 1 2 ) 2 (m1,m2) =( 1 2 , 1 2 ) q m 2 1 -1 4 1 q m 2 2 -1 4 2 • 2 cosh(2πm 1 m 2 y 3 ). Since |x 3 | ≤ 1 2
, the absolute value of the argument of exp(iπz 3 /2) is at most π/4. Therefore, exp(iπz 3 /2) + exp(-iπz 3 /2) ≥ exp(π |y 3 | /2).

We obtain

θ 12 (τ ) ξ 12 (τ ) -1 ≤ (m1,m2)∈(N+ 1 2 ) 2 (m1,m2) =( 1 2 , 1 2 ) q m 2 1 -1 4 1 q m 2 2 -1 4 2 • 2 cosh 2π(m 1 m 2 -1 4 )y 3 .
We separate the terms corresponding to m 2 = 1 2 . Since 2 |y 3 | ≤ y 1 , their contribution is bounded by

m1∈N+ 3 2 q m 2 1 -1 2 m1 1 + q m 2 1 + 1 2 m1-1 2 1 ≤ q 3/2 1 1 -q 7/2 1 + q 5/2 1 1 -q 9/2 1
.

Similarly, the contribution from the terms with m 1 = 1/2 is bounded by

q 3/2 2 1 -q 7/2 2 + q 5/2 2 1 -q 9/2 2
.

For the remaining terms, we use the majoration

2π(m 1 m 2 -1 4 )y 3 ≤ 2πm 1 m 2 y 3 ≤ π |y 3 | (m 2 1 + m 2 2 ) ≤ π 2 (m 2 1 y 1 + m 2 2 y 2 ).
Thus, the rest of the sum is bounded by

m1,m2∈N+ 3 2 q m 2 1 -1 4 1 q m 2 2 -1 4 2 • 2 cosh π 2 (m 2 1 y 1 + m 2 2 y 2 ) ≤ m1,m2∈N+ 3 2 q 1 2 m 2 1 -1 4 1 q 1 2 m 2 2 -1 4 2 + q 3 2 m 2 1 -1 4 1 q 3 2 m 2 2 -1 4 2 ≤ q 7/8 1 q 7/8 2 (1 -q 2 1 )(1 -q 2 2 ) + q 25/8 1 q 25/8 2 (1 -q 6 1 )(1 -q 6 2 )
.

We give another version of these estimates that we will use for τ

(n) 4 .
Lemma 4.6. Let k ≥ 2, and let τ ∈ H 2 such that

y 3 (τ ) 2 ≤ 3 7 y 1 (τ )y 2 (τ ) and k |y 3 (τ )| ≤ y 1 (τ ). Let α = 3/7. Define ρ ′(k) 0,1 (q 1 , q 2 ) = 2q 4 2 1 -q 5 2 + 2q 1 1 -q 3 1 + 2q 1-2 k 1 q 2 1 -q 3-2 k 1 + 2q 1+ 2 k 1 q 2 1 -q 3+ 2 k 1 + 2q 1-α 1 q 4(1-α) 2 (1 -q 3(1-α) 1 )(1 -q 5(1-α) 2 ) + 2q 1+α 1 q 4(1+α) 2 (1 -q 3(1+α) 1 )(1 -q 5(1+α) 2 ) and ρ ′(k) 8,9 (q 1 , q 2 ) = q 2 2 1 -q 4 2 + q 1-1 k 1 1 -q 3-1 k 1 + q 1+ 1 k 1 1 -q 3+ 1 k 1 + q 2-9 4 α 2 q 1-α 1 (1 -q 4(1-α) 2 )(1 -q 3(1-α) 1 ) + q 2+ 9 4 α 2 q 1+α 1 (1 -q 4(1+α) 2 )(1 -q 3(1+α) 1
)

.

Then we have

|θ j (τ ) -ξ 0,1 (τ )| ≤ ρ ′ 0,1 (τ ) for j ∈ {0, 1} and θ j (τ ) ξ 8,9 (τ ) -1 ≤ ρ ′ 8,9 (τ ) for j ∈ {8, 9}.
Proof. We bound the cross-product terms by

2y 3 m 1 m 2 ≤ αy 1 m 2 1 + αy 2 m 2 2 , 2y 3 m 1 m 2 ≤ 1 k y 1 m 1 if m 2 = 1 2
, and

2y 3 m 1 m 2 ≤ 2 k y 1 m 1 if m 2 = 1.
For j ∈ {0, 1}, we separate the terms with |m 2 | ≤ 1 or m 1 = 0, and obtain

|θ j (τ ) -ξ 0,1 (τ )| ≤ 2 m2≥2 q m 2 2 2 + 2 m1≥1 q m 2 1 1 + 2 m2≥1 q 2 (q m 2 1 -2 k m1 1 + q m 2 1 + 2 k m1 1 ) + 2 m1≥1 m2≥2 q m 2 1 1 q m 2 2 2 • 2 cosh α(y 1 m 2 1 + y 2 m 2 2 ) ≤ 2q 4 2 1 -q 5 2 + 2q 1 1 -q 3 1 + 2q 1-2 k 1 q 2 1 -q 3-2 k 1 + 2q 1+ 2 k 1 q 2 1 -q 3+ 2 k 1 + 2q 1-α 1 q 4(1-α) 2 (1 -q 3(1-α) 1 )(1 -q 5(1-α) 2 ) + 2q 1+α 1 q 4 (1+α) 2 (1 -q 3(1+α) 1 )(1 -q 5(1+α) 2
) .

For j ∈ {8, 9}, we separate the terms with

|m 2 | = 1 2 or m 1 = 0. We obtain θ j (τ ) ξ 8,9 (τ ) -1 ≤ q -1/4 2 m2∈N+ 3 2 q m 2 2 2 + m1≥1 q m 2 1 -1 k m1 1 + q m 2 1 + 1 k m1 1 + q -1/4 2 m2∈N+ 3 2 m1≥1 q m 2 2 2 q m 2 2 1 • 2 cosh α(y 1 m 2 1 + y 2 m 2 2 ) ≤ q 2 2 1 -q 4 2 + q 1-1 k 1 1 -q 3-1 k 1 + q 1+ 1 k 1 1 -q 3+ 1 k 1 + q 2-9 4 α 2 q 1-α 1 (1 -q 4(1-α) 2 (1 -q 3(1-α) 1 ) + q 2+ 9 4 α 2 q 1+α 1 (1 -q 4(1+α) 2 )(1 -q 3(1+α) 1
) .

Finally, when n is large, we will show that the theta constants θ j (2 n γ k τ ) for 0 ≤ j ≤ 3 are in good position using the following lemma. Recall the definition of r(τ ) and λ 1 (τ ) from §2. 1. If r(τ ) ≥ 0.4, then the θ j (τ ) for 0 ≤ j ≤ 3 are in good position.

2. If λ 1 (τ ) ≥ 0.6, then the θ j (τ ) for 0 ≤ j ≤ 3 are in good position.

Proof.

1. Write q = exp(-πr(τ )).

For 0 ≤ j ≤ 3, we have

|θ j (τ ) -1| ≤ 4q 2 + n∈Z 2 , n 2 ≥2 exp(-πλ 1 (τ ) n 2 ) ≤ 8q 2 + 4q 4 + 8q 5 + 4q 8 + 4 1 + q (1 -q) 2 q 9 . (6)
In this inequality, the first term 4q 2 comes from the four vectors n ∈ Z 2 with n = 1. Then we separate the terms n = ( n1 n2 ) such that |n 1 | ≥ 3 and |n 2 | ≥ 3; this accounts for the term 4q 9 (1 + q)/(1 -q) 2 , as in the proof of [START_REF] Dupont | Moyenne arithmético-géométrique, suites de Borchardt et applications[END_REF]Prop. 6.1]. We leave the remaining terms as they are.

If q ≤ 0.287, then the quantity on the right hand side of ( 6) is less than √ 2/2, and the θ j (τ ) are contained in a disk which is itself contained in a quarter plane. We have q ≤ 0.287 when r(τ ) ≥ 0.4.

Write

q = exp(-πλ 1 (τ )).

Then for 0 ≤ j ≤ 3, we have |θ j (τ ) -1| ≤ 4q + 4q 2 + 4q 4 + 8q 5 + 4q 8 + 4 1 + q (1 -q) 2 q 9 . This quantity is less than √ 2/2 when λ 1 (τ ) ≥ 0.6.

We conclude this section with lower bounds on r or λ 1 at γ k τ for τ ∈ F ′ and 1 ≤ k ≤ 3. Proof. We have Therefore,

Im(γ 1 τ ) = z 1 z 3 0 -1 -t Im(τ ) z 1 z 3 0 -1 -1 = 1 |z 1 | 2 y 1 α α β with α = y 1 x 3 -y 3 x 1 , so |α| ≤
λ 1 (γ 1 τ ) ≥ det Im(γ 1 τ ) Tr Im(γ 1 τ ) ≥ y 1 |z 1 | 2 1 1 + 25 16 y 2 1 |z1| 2 det Im(τ ) ≥ 9y 1 34 |z 1 | 2 .
We did not use the property that y 1 ≤ y 2 , so the same proof works for γ 2 τ . Finally, we consider γ 3 τ . We have

Im(γ 3 τ ) = 1 |det τ | 2 β 1 α α β 2 with β 1 = y 1 |z 3 | 2 + y 2 |z 1 | 2 -y 3 (z 1 z 3 + z 3 z 1 )
,

β 2 = y 1 |z 2 | 2 + y 2 |z 3 | 2 -y 3 (z 2 z 3 + z 3 z 2 ).
We compute

|det τ | 2 Tr Im(γ 3 τ ) = β 1 + β 2 ≤

Proof of the main theorem

In this final section, we prove Theorem 1.1 by separating different cases according to the value of n. If n is large enough, then Lemmas 4.7 and 4.8 are enough to conclude; if n is smaller, then we apply the theta transformation formula (Proposition 3.2) and the bounds on other theta constants given in §4.

In the proofs, we use numerical calculations, typically in order to show that a given angle α(q) is smaller than π/2 for certain values of q. Such calculations are easily certified using interval arithmetic, since the functions α(q) we consider are simple: they are either increasing or convex functions of q.

In order to help the reader visualize the estimates, we created pictures using GeoGebra [12]. Proposition 5.1. Let τ ∈ F ′ . Then for every n ≥ 0, the theta constants θ j (2 n τ ) for 0 ≤ j ≤ 3 are in good position.

Proof. For every n ≥ 0, we have

r(2 n τ ) = 2 n r(τ ) ≥ √ 3/4 ≥ 0.4,
so the result follows from Lemma 4.7.

Lemma 5.2. Let τ ∈ F ′ .

1. For every n ≥ 0 such that 2 n ≤ 8.77y 1 (τ ), the theta constants θ j (τ

(n)
1 ) for j ∈ {0, 2, 4, 6} are in good position.

2. For every n ≥ 0 such that 2 n ≤ 8.77y 2 (τ ), the theta constants θ j (τ (n)

2 ) for j ∈ {0, 1, 8, 9} are in good position.

Proof. We only prove the first statement, the second one being symmetric. We separate three cases: n = 0, n = 1, and n ≥ 2. Case 2: n = 1. We study the relative positions of ξ 0,2 and ξ 4,6 at τ

1 . Since |2 -n x 1 (τ )| ≤ 1/4, the absolute value of the argument of ξ 4,6 (τ (1) 1 ) is bounded above by π/16. Moreover,

ξ 0,2 (τ (1) 1 ) ≥ 1, arg(ξ 0,2 (τ (1) 1 )) ≤ arctan 2q 1 sin(π/4) 1 + 2q 1 cos(π/4) ,
and the arguments of ξ 0,2 and ξ 4,6 have the same sign. Therefore the angle between any two θ j (τ

1 ) for j ∈ {0, 2, 4, 6} is at most max π 16 , arctan 2q 1 sin(π/4) 1 + 2q 1 cos(π/4) +arcsin ρ

4,6 (q 1 , q 2 )+arcsin(ρ 0 (q 1 , q 2 )+2q 2 )

by Lemmas 4.2 and 4.4. This quantity is less than π/2 because q 2 (τ

1 ) ≤ exp(-π √ 3) and q 1 (τ

1 ) ≤ exp(-π √ 3/8).

Case 3: n ≥ 2. We proceed as in Case 2, but we now have

q 2 (τ (n) 1 ) ≤ exp(-2π √ 3), 8 y 3 (τ (n) 1 ) ≤ y 2 (τ (n)
1 ), and x 1 (τ

(n) 1 ) ≤ 1 8 .
Therefore the angle between the θ j (τ

(n) 1 ) for j ∈ {0, 2, 4, 6} is bounded by max π 32 , arctan 2q 1 sin(π/8) 1 + 2q 1 cos(π/8) + arcsin(ρ 0 (q 1 , exp(-2π √ 3)) + 2 exp(-2π √ 3)) + arcsin ρ (8) 4,6 (q 1 , exp(-2π √ 3)). (7)
This angle remains less that π/2 when q 1 (τ

(n)
1 ) ≤ 0.699. This is the case when 2 n ≥ 8.77y 1 (τ ).

The geometric situation in Case 3 of Lemma 5.2 can be represented as follows.

In this picture, we take q 1 = 0.699, and represent two complex numbers ξ 0,2 and ξ 4,6 with modulus one, separated by an angle of max π 32 , arctan 2q 1 sin(π/8) 1 + 2q 1 cos(π/8) ≃ 0.22.

Then we draw disks centered in ξ 0,2 and ξ 4,6 with radii ρ 0 (q 1 , exp(-2π √ 3)) and ρ (8) 4,6 (q 1 , exp(-2 √ 3)) respectively. Finally we represent the smallest angular sector seen from the origin containing these two disks. The green angle is equal to the quantity [START_REF] Enge | Computing modular polynomials in quasi-linear time[END_REF], and is indeed smaller than π/2. Proposition 5.3. Let τ ∈ F ′ . 1. For every n ≥ 0, the theta constants (θ j (2 n γ 1 τ )) 0≤j≤3 are in good position.

2. For every n ≥ 0, the theta constants (θ j (2 n γ 2 τ )) 0≤j≤3 are in good position. because y 1 (τ ) ≥ √ 3/2. By Lemma 4.7, the θ j (2 n γ 1 τ ) for 0 ≤ j ≤ 3 are in good position when 2 n r(γ 1 τ ) ≥ 0.4. This is the case when 2 n ≥ 1.96y 1 . On the other hand, Lemma 5.2 applies when 2 n ≤ 8.77y 1 . The second statement is proved in the same way.

Lemma 5.4. Let τ ∈ F ′ . Then, for every n ≥ 0 such that 2 n ≤ 1.66y 1 , the theta constants θ j (τ

(n)
3 ) for j ∈ {0, 4, 8, 12} are in good position.

Proof. Write q = q 1 (τ

(n)
3 ) for short. We separate two cases: n ≥ 1, and n = 0.

Case 1: n ≥ 1. In this case, we have

x j (τ (n) 3 ) ≤ 1/4 for each 1 ≤ j ≤ 3.
Therefore, given the expressions of ξ 0 , ξ 4,6 , ξ 8,9 and ξ 12 (see ( 4)), and by Lemmas 4.2 to 4.5,

• The angle between θ 4 (τ 4,6 (q, q).

• The angle between θ 4 (τ

(n)
3 ) (or θ 8 ) and θ 0 (τ

(n)
3 ) is bounded by π 16 + arcsin ρ

4,6 (q, q) + 2q sin(π/4) + arcsin ρ 0 (q, q).

• The angle between θ 12 (τ

(n)
3 ) and θ 4 (τ

(n)
3 ) (or θ 8 ) is bounded by 3π 16 + arcsin ρ 12 (q, q) + arcsin ρ

(2) 4,6 (q, q).

• The angle between θ 12 (τ

(n)
3 ) and θ 0 (τ

(n)
3 ) is bounded by π 4 + arcsin ρ 12 (q, q) + arcsin ρ 0 (q, q). All these quantities remain less than π/2 when q ≤ 0.151. This is the case when 2 n ≤ 1.66y 1 .

Case 2: n = 0. In this case, we have q ≤ exp(-π √ 3/2). Therefore,

• The angle between θ 4 and θ 8 is bounded by

π 4 + 2 arcsin ρ (2)
4,6 (q, q) < π 2 .

• The angle between θ 4 (or θ 8 ) and θ 0 is bounded by

π 8 + arcsin ρ (2)
4,6 (q, q) + arcsin(ρ 0 (q, q) + 4q) < π 2 .

• The angle between θ 12 and θ 4 (or θ 8 ) is bounded by 3π 8 + arcsin ρ 12 (q, q) + arcsin ρ

(2) 4,6 (q, q) < π 2 .

These estimations can be represented as follows, with similar conventions as in the picture after Lemma 5.2:

We finally study the angle between θ 12 and θ 0 . The argument of ξ 12 (τ ) is x 1 /4 + x 2 /4 + β with β = arg(exp(iπz 3 /2) + exp(-iπz 3 /2)). Up to conjugating, we may assume that y 3 ≥ 0 and x 3 ≥ 0. Then exp(iπz 3 /2) + exp(-iπz 3 /2) = exp(-iπz 3 /2)(1 + exp(iπz 3 )) so β + πx 3 2 ≥ arctan q 3 sin(πx 3 ) 1 + q 3 ≥ arctan 2x 3 q 3 1 + q 3 .

In general, we have |β| ≤ π 4 arctan q 1/2 1 + q 1/2 .

  (n) b for b ∈ (Z/2Z) 2 are in good position.

1 .

 1 For each 0 ≤ k ≤ 3, compute the first term of the sequence B(γ k τ )/θ 2 0 (γ k τ ) using the transformation formulas for theta constants under Sp 4 (Z) (see Igusa [15, Thm. 2 p. 175 and Cor. p. 176], or Corollary 3.3); 2. For each 0 ≤ k ≤ 3, compute 1/θ 2 0 (γ k τ ) as the limit of the Borchardt sequence B(γ k τ )/θ 2 0 (γ k τ ); 3. Use the input and the newly computed θ 2 0 (γ 0 τ ) = θ 2 0 (τ ) to compute all squares of theta constants at τ ; 4. Recover τ = z 1 z 3 z 3 z 2 using the relations given in [4, §6.3.1]:

  Proposition 3.2 ([15, Thm. 2 p. 175 and Cor. p. 176]). Let a, b ∈ {0, 1} 2 , and let

Lemma 4 . 3 . 4 y 1

 4341 Let k ≥ 1, and let τ ∈ H 2 such that y 3 (τ ) 2 ≤ 1 (τ )y 2 (τ ) and k |y 3 (τ )| ≤ y 1 (τ ).

Lemma 4 . 7 .

 47 Let τ ∈ H 2 .

Lemma 4 . 8 .

 48 For every τ ∈ F ′ , we haver(γ 1 τ ) ≥ 9 y 1 (τ ) 34 |z 1 (τ )| 2 , r(γ 2 τ ) ≥ 9 y 2 (τ ) 34 |z 2 (τ )| 2 , and λ 1 (γ 3 τ ) ≥ 9 44 y 2 (τ ).

Case 1 :

 1 n = 0. Then τ (n) 1 = τ . By [24, Prop. 7.7], we have |θ j (τ ) -1| ≤ 0.405 for j ∈ {0, 1, 2, 3}, andθ j (τ ) ξ 4,6 (τ ) -1 ≤ 0.348 for j ∈ {4, 6}.The absolure value of the argument of ξ 4,6 (τ ) is at most π/8. Therefore the angle between any two θ j (τ ) for j ∈ {0, 1, 2, 3, 4, 6} is at most π 8 + arcsin(0.348) + arcsin(0.405) < π 2 .

Proof. 1 34 |z 1 | 2 ≥

 12 By Lemma 4.8, we have r(γ 1 τ ) ≥ 9 y

On the other hand, 1 ≤ Re(ξ 0 (τ )) ≤ 1 + 4q, Im(ξ 0 (τ )) = 2q 1 sin(πx 1 ) + 2q 2 sin(πx 2 ).

We discuss two cases according to the signs of x 1 and x 2 :

• If x 1 and x 2 have opposite signs, then the angle between θ 12 and θ 0 is at most 3π 8 + arctan(2q) + arcsin ρ 12 (q, q) + arcsin ρ 0 (q, q).

• If x 1 and x 2 have the same sign, say positive, then

Therefore the angle between θ 12 and θ 0 is at most π 2 arctan q 1/2 1 + q 1/2 + arcsin ρ 12 (q, q) + arcsin ρ 0 (q, q). This function of q is not increasing, but it is convex.

A numerical investigation shows that both quantities remain less than π/2 when q ≤ exp( √ 3/2).

Lemma 5.5. Let τ ∈ F ′ , and let n 0 ∈ N such that 2 n0 > 1.66y 1 . Then, for every n ≥ n 0 such that 2 n ≤ 4.2y 2 (τ ), the theta constants θ j (τ

4 ) for j ∈ {0, 1, 8, 9} are in good position.

Proof. By assumption, we have y 1 (τ

4 ) ≤ 0.021. Moreover we must have n ≥ 1, so by Lemma 3.4, x 2 (τ (n) 4 ) ≤ 9/16, and

Therefore, we can apply Lemma 4.6 with k = 8/3: we have

(0.021, q 2 (τ

Let us investigate the difference between the arguments of ξ 8,9 (τ

) and ξ 0,1 (τ

4 ), which we may assume to be positive. If the argument of ξ 8,9 is the largest, then the difference is bounded by

.

If the argument of ξ 0,1 is the largest, we distinguish two cases. If x 2 (τ

On the other hand, if x 2 (τ ) ≤ 3π/8, then arg ξ 0,1 (τ

) is always greater than cos(π/16). Therefore the angle between the θ j (τ (0.021, q 2 ) + arcsin ρ ′(8/3) 0,1 (0.021, q 2 ) cos(π/ 16) .

This quantity is less than π/2 when q 2 (τ

On the other hand, if 2 n > 2.43y 2 (τ ), then we must have n ≥ 2. Moreover, y

4 ) > 2.43 (0.0033, q 2 ) + arcsin ρ ′(16/3) 0,1 (0.0033, q 2 ).

This quantity is less than π/2 when q 2 (τ (n) 4 ) ≤ 0.571, and the latter inequality holds when 2 n ≤ 4.2y 2 (τ ). Proposition 5.6. Let τ ∈ F ′ . Then, for every n ≥ 0, the theta constants θ j (2 n γ 3 τ ) for 0 ≤ j ≤ 3 are in good position.

Proof. By Lemma 4.8, we have

.

Therefore, by Lemma 4.7, the theta constants are in good position as soon as 2 n 9 44y 2 (τ ) ≥ 0.6, or 2 n ≥ 2.94y 2 (τ ).

When n is smaller, we use the transformation formulas. Lemma 5.4 applies when 2 n ≤ 1.66y 1 (τ ), and Lemma 5.5 applies when 1.66y 1 (τ ) < 2 n ≤ 4.2y 2 (τ ).


R Pohl

G Kaiser

Labeyrie

M Gaudesius

Y.-C Zhang

T Pohl

R Kaiser

G Labeyrie

Three-dimensional simulations of spatiotemporal instabilities in a magneto-optical trap

Three-dimensional simulations of spatiotemporal instabilities in a magneto-optical trap

I. Introduction

A magneto-optical trap (MOT) is nowadays extensively used to study properties of neutral atoms, or as a part of more elaborate setups for creating, e.g., quantum degenerate gases. The limit of a large atom number N , in particular, has attracted interest in various applications, notably random lasing [START_REF] Baudouin | A cold-atom random laser[END_REF], Anderson localization [START_REF] Cottier | Microscopic and macroscopic signatures of 3D Anderson localization of light[END_REF], selforganization [START_REF] Labeyrie | Optomechanical self-structuring in a cold atomic gas[END_REF], superradiance [START_REF] Araújo | Superradiance in a large and dilute cloud of cold atoms in the linear-optics regime[END_REF] and subradiance [START_REF] Guerin | Subradiance in a large cloud of cold atoms[END_REF]. This limit can, moreover, allow for investigations of nonlinear phenomena that possess similarities with pulsating stars [START_REF] Labeyrie | Self-sustained oscillations in a large magneto-optical trap[END_REF][START_REF] Cox | Theory of stellar pulsation[END_REF] and unstable plasmas [START_REF] Mendonça | Collective oscillations in ultracold atomic gas[END_REF]- [START_REF] Hansen | Statistical mechanics of dense ionized matter. III. Dynamical properties of the classical one-component plasma[END_REF].

At low N , a MOT is governed by single-atom physics, such that each atom in the cloud is independently subjected to a cooling and confining force applied by the laser beams in the presence of the magnetic field. At large N , many-atom physics become important as collective forces appear, e.g., the shadow force [START_REF] Dalibard | Laser cooling of an optically thick gas: The simplest radiation pressure trap?[END_REF] and the rescattering force [START_REF] Sesko | Behavior of neutral atoms in a spontaneous force trap[END_REF]. The shadow force is caused by an imbalance of beam intensities in the cloud due to attenuation, which occurs because of light's scattering as it traverses through the cloud. This force is compressive but is countered by the repulsive, Coulomb-like force caused by photon rescattering between the atoms. With these antagonistic forces present, large MOT clouds can exhibit spatiotemporal instabilities in the form of spontaneous oscillations.

Spatiotemporal instabilities have been studied in several MOT configurations [START_REF] Labeyrie | Self-sustained oscillations in a large magneto-optical trap[END_REF][START_REF] Sesko | Behavior of neutral atoms in a spontaneous force trap[END_REF]- [START_REF] Kim | Observation of the Hopf bifurcation in parametrically driven trapped atoms[END_REF], and many theoretical models have been explored to gain insight into the experimentally observed features [START_REF] Pohl | Self-driven nonlinear dynamics in magneto-optical traps[END_REF]- [START_REF] Moon | Nonlinear, nonequilibrium and collective dynamics in a periodically modulated cold atom system[END_REF]. In our case, the socalled balanced MOT configuration is considered, where the laser beams are independent and have the same, constant intensities before entering the cloud. In our recent works [START_REF] Gaudesius | Instability threshold in a large balanced magnetooptical trap[END_REF] and [START_REF] Gaudesius | Phase diagram of spatiotemporal instabilities in a large magneto-optical trap[END_REF], we have made experimental observations of instability thresholds and regimes, respectively, and reached qualitative agreements with results of our three-dimensional (3D) simulations. In the present work, we explain the kinetic model and the technicalities employed in these simulations. Moreover, the simulations are employed to investigate the role of different effects on the instabilities.

The remainder of this article is structured as follows. In Sec. II.A, the kinetic model employed in the 3D simulations is described, followed up by Sec. II.B, where details surrounding its implementation are covered. Then, Sec. III investigates the impact of its physical effects on the instabilities. Finally, Sec. IV provides a conclusion and a discussion on the future perspectives.

II. Simulation method and approximations

A. Kinetic model

In the development of our kinetic model, the considered conventions for the MOT magnetic field and beams are displayed in Fig. 1(a). The magnetic field, which is quadrupole in nature, is assumed to be B(r) = B -x 2 , -y 2 , z , where B > 0 is the field gradient along the z axis, and r = (x, y, z) is the atom position with the origin at the trap center. This assumption holds when |r| is much smaller than both the radius of the MOT coils and the separation between them. In accordance with the chosen magnetic field convention, the z-axis beams have right-handed circular helicity, whereas the xand y-axis beams have left-handed circular helicity.

We base our model on the hyperfine transition F = 0 → F = 1. It is evidently much simpler than in our experiments, F = 2 → F = 3 (totaling 12 Zeeman sublevels) [START_REF] Gaudesius | Instability threshold in a large balanced magnetooptical trap[END_REF][START_REF] Gaudesius | Phase diagram of spatiotemporal instabilities in a large magneto-optical trap[END_REF], but allows for a proper description of the features related to the magnetic field and light polariza-

(b) (a) x z y m = 0 m' = -1 m' = 0 m' = +1 σ + σ - π m = 0 m = 0 m = 0 σ + σ - π F' = 1 F = 0 μB(r) μB(r) m' = 0 m' = -1 m' = +1 μB(r) μB(r) ⇒ Figure 1: (a)
Sketch of a MOT, displaying the considered conventions for the MOT magnetic field and beams in the derivation of the kinetic model employed in our 3D simulations of instabilities. The coils produce a magnetic field (curved arrows) that has a positive gradient along the z axis. Accordingly, the beams have the following circular helicities (circular arrows): Right-handed for the z-axis beams, and left-handed for the xand y-axis beams. (b) Simplification of the Zeeman sub-level structure of the hyperfine transition F = 0 → F = 1, in the employed kinetic model. Each of the three Zeeman transitions m = 0 → m = -1, 0, +1 between F = 0 and F = 1 is treated as an independent 2-level system (right picture). They are induced by, respectively, only σ -, π, σ + polarized light. The MOT magnetic field leads to the Zeeman splitting of the hyperfine levels depending on the atom position [µB(r)].

tion. As depicted in Fig. 1(b), each of the three Zeeman transitions m = 0 → m = -1, 0, +1 between the hyperfine levels is treated as an independent 2-level system and driven by, respectively, σ -, π, σ + polarized light. We expect this approximation to hold only in the regime of low saturation [START_REF] Hanley | Quantitative simulation of a magneto-optical trap operating near the photon recoil limit[END_REF], s = I∞/Isat

1+ 4∆ 2 Γ 2
1, where I ∞ is the intensity of a beam before entering the cloud, I sat is the saturation intensity, ∆ = ω L -ω 0 is the detuning of the laser frequency ω L from the atomic transition m = 0 → m = 0 frequency ω 0 , and Γ is the natural linewidth.

We also neglect sub-Doppler effects, which seems to be reasonable considering that large MOT clouds have been observed to be mostly unaffected by these effects [START_REF] Kim | Direct observation of the sub-Doppler trap in a parametrically driven magneto-optical trap[END_REF].

We include the following main four physical effects in our model: (i) the mean cooling and confining force stemming from the MOT beams, hereafter referred to as the trapping force; (ii) the diffusion resulting from its fluctuations; (iii) the beam intensity attenuation caused by the light's scattering in the cloud; and (iv) the rescattering force due to photon exchange between the atoms. In the following, we explain how these effects are described.

Trapping force

To describe the trapping force, we use the standard Doppler model. It relies on an assumption, expected to be valid for s 1, that this force can be expressed as a sum of six radiation pressure forces, one for each beam:

F tr (r, v) = α=x,y,z F + α (r, v) + F - α (r, v) (1)
where v = (v x , v y , v z ) is the atom velocity, and F + α and F - α are the radiation pressure forces exerted by the beams traveling in, respectively, positive and negative directions of α = x, y, z, given by

F ± α (r, v) = F ± α,-(r, v) + F ± α,0 (r, v) + F ± α,+ (r, v) (2)
where F ± α,-, F ± α,0 , F ± α,+ are the radiation pressure forces acting on the atom's 2-level transitions that are driven by, respectively, σ -, π, σ + polarized light, given by

F ± α,q (r, v) = ± p ± α,q (r)I ± α (r, v)σ ± α,q (r, v) c α
(3) In this equation, q = -, 0, + refers to the respective σ -, π, σ + transitions, c is the vacuum light speed, and the remaining quantities are defined as follows.

The coefficient p ± α,q denotes the fraction of the positive or negative (±) α = x, ŷ, ẑ directed light that drives the σ -, π, or σ + (q = -, 0, or +) transition [START_REF] Brickman | Magneto-optical trapping of cadmium[END_REF]. It introduces anisotropy to the trapping force. There is a total of 18 fractions, as there are three transitions corresponding to each of the six MOT beams:

p ± α,q (r) =          1 2 1 ± α B 2B(r) 2 , q = + (σ + ) 1 2 1 ∓ α B 2B(r) 2 , q = -(σ -) 1 -(p ± α,+ + p ± α,-) , q = 0 (π) (4
) where α = x, y, 2z for, respectively, α = x, y, z, B(r) = B z 2 + 1 4 (x 2 + y 2 ) is the magnitude of the magnetic field B(r). The quantization axis has been chosen to be along the direction of B(r). Consequently, as shown in Fig. 1(b), the Zeeman shifts of the excited levels with m = -1, 0, +1 are given by µ q (r) = qµB(r), where q = -, 0, +, respectively, and µ is the gyromagnetic ratio.

For the final quantities in Eq. (3), I ± α is the corresponding beam intensity, which is subject to attenuation as will be covered later, and σ ± α,q is the corresponding scattering cross-section for a single 2-level transition of the atom, given by

σ ± α,q (r, v) = σ 0 1 + Itot,q(r,v) Isat + 4 (∆∓k L vα-µq(r)) 2 Γ 2
(5) where σ 0 = 6π/k 2 L is the resonant scattering cross section, where k L = ω L /c denotes the laser wavenumber; ∓k L v α is the Doppler shift for a positive or negative (±) beam; µ q is the Zeeman shift defined previously; and

I tot,q (r, v) = α=x,y,z p + α,q (r)I + α (r, v) + p - α,q (r)I - α (r, v) (6)
is the total beam intensity that a single 2-level transition receives, which is observed to be generally different for all three (σ -, π, and σ + ) due to p ± α,q . Because this intensity enters into scattering cross sections, the beam cross-saturation effect gets naturally introduced into our model. In the section on attenuation, we mention a difficulty (luckily, solvable) that is encountered by having this effect included.

Diffusion

The fluctuating part of the trapping force can be introduced via a momentum diffusion coefficient. Here we describe such diffusion processes approximately by adapting the known description for a 2-level atom in a single laser beam [START_REF] Cohen-Tannoudji | Atomic motion in laser light[END_REF] to our 3D situation involving six laser beams. Explicitly, we write the following for the momentum diffusion coefficient:

D(r, v) = D vac (r, v) + D las (r, v) (7)
where

D vac (r, v) = 2 k 2 L Γ 4
s tot (r, v) 1 + s tot (r, v) (8)
D las (r, v) = 2 k 2 L Γ 4
s tot (r, v) [1 + s tot (r, v)] 3 × 1 + 12∆ 2 -Γ 2 4∆ 2 + Γ 2 s tot (r, v) + s 2 tot (r, v) (9)
are the momentum diffusion coefficients of, respectively, the vacuum and laser field, with being the reduced Planck constant and s tot being the total saturation parameter that is a sum of the total saturation parameters s tot,-, s tot,0 , s tot,+ for the atom's 2-level transitions that are driven by, respectively, σ -, π, σ + polarized light, i.e.,

s tot (r, v) = q=-,0,+ s tot,q (r, v) (10)
where

s tot,q (r, v) = α=x,y,z s + α,q (r, v) + s - α,q (r, v) (11)
and

s ± α,q (r, v) = p ± α,q (r)I ± α (r, v)/I sat 1 + 4 (∆∓k L vα-µq(r)) 2 Γ 2 (12)
denotes the saturation parameter for a single beam and atomic transition.

Note that the diffusion is affected by the attenuation of the laser beams and, therefore, is a collective effect that depends on N [START_REF] Drewsen | Investigation of sub-Doppler cooling effects in a cesium magneto-optical trap[END_REF][START_REF] Cooper | The temperature of atoms in a magneto-optical trap[END_REF]. While our model necessarily involves an approximate treatment of diffusion and its dependence on the attenuation, its inclusion is nevertheless important for unstable balanced MOTs [START_REF] Heimonen | Synchronization of a self-sustained cold-atom oscillator[END_REF]. In Sec. III, we discuss the significance of diffusion in more detail.

Beam attenuation

The beam attenuation modifies the properties of all the main physical effects of our model, including itself (see below). Specifically for trapping [see Eq. ( 3)], the attenuation's inclusion results in an additional compression known as the shadow force, as illustrated in Fig. 2(a). Note that the attenuation in our model is fully justified considering that the experimentally observed optical depths at the instability threshold are typically on the order of 1 (along the z axis). To describe this effect, we employ the low-saturation-regime assumption, s 1, under which a given beam intensity decays exponentially as

I ± α (r, v) = I ∞ e -O ± α (r,v) (13)
where O ± α is the optical depth of the cloud for a given beam. For, e.g., the positive ẑ directed beam,

O + z (r, v) = z -∞ dz ρ(x, y, z ) × q=-,0,+ p + z,q (x, y, z )σ + z,q (x, y, z , v(x, y, z )) ( 14
)
where ρ is the cloud density, and v(x, y, z ) refers to the velocity at (x, y, z ). For the remaining five MOT beams, analogous expressions of the optical depth can readily be obtained. Note that for the negative α = x, ŷ, ẑ beam, the integral limit is from, respectively, x, y, z to +∞.

We point out that calculating intensity requires prior knowledge of I tot,q , due to it entering into scattering cross sections [Eq. ( 5)]. Although the precise determination of I tot,q is numerically demanding, this difficulty can be overcome by implementing an iterative procedure that is presented in Sec. II.B; thus, the beam cross saturation can be included. Moreover, in Sec. III, we show how attenuation impacts the instabilities.

Rescattering force

To understand the origin of the rescattering force, let us consider the basic situation in Fig. 2 has Coulomb-like character, as the rescattered intensity I R is diminished according to the inverse-square law of propagation. Importantly, the scattering cross section of atom 1, σ L , is different from the rescattering cross section of atom 2, σ R , due to, e.g., inelastic scattering at atom 1. In the well-known Wieman model for the multiplescattering regime [START_REF] Sesko | Behavior of neutral atoms in a spontaneous force trap[END_REF], one has σ R > σ L , and thus the cloud expands as N is increased; this inequality relies critically on the presence of inelastic scattering in the cloud.

In our model, the rescattering is more complex than described above, and in Fig. 3 we provide an illustration of our employed approximate description. As atom 1 is exposed to all the MOT beams, each of its σ -, π, σ + transitions produces a characteristic radiation pattern of scattered light. At atom 2, each transition of this atom rescatters the scattered light by an amount depending on the fractions that play a role similar to that of the fractions in Eq. ( 4). Under this description, the rescattering force on an atom with the position r j and the velocity v j , due to surrounding atoms with the positions r l and the velocities v l , is

atom 1 atom 2 𝑏 - - 𝑏 0 - 𝑏 + - 𝜂 - 𝜂 0 𝜂 + σ -(r 1 ) π(r 1 ) σ + (r 1 ) σ -(r 2 ) π(r 2 ) σ + (r 2 ) 𝐼 z - 𝐼 z + 𝐼 x + 𝐼 y + 𝐼 y - 𝐼 x -
F rsc (r j , v j ) = l =j P R (r l , r j , v l , v j ) c rl,j (15)
where rl,j is the unit vector corresponding to r l,j = r j -r l , which points from atom l to atom j, and P R is the power rescattered by atom j due to atom l. It is given by

P R (r l , r j , v l , v j ) = q =-,0,+ I R,q (r l , r j , v l )σ R,q (r l , r j , v l , v j ) = q =-,0,+ q =+,0,- b q q (r l , r j )I S,q (r l , r j , v l ) σ R,q (r l , r j , v l , v j ) (16)
where q and q refer to the 2-level transitions of, respectively, atom l and atom j. The sum in the brackets is the intensity I R,q rescattered by a single 2-level transition of atom j. The coefficient b q q denotes the fraction of the scattered radiation with the intensity I S,q that drives its σ -, π, or σ + (q = -, 0, or +) transition. Finally, σ R,q is the rescattering cross section for a single 2-level transition of this atom.

The coefficients are given by

b ± q (r l , r j ) =                      1 2u     n=l,j rl,j • Bn   ± 1   2 , q = + (σ + ) 1 2u     n=l,j rl,j • Bn   ∓ 1   2 , q = -(σ -) 1 -(b ± + + b ± -) , q = 0 (π) (17) b 0 q (r l , r j ) =          1 2 rl,j • Bj 2 , q = + (σ + ) 1 2 rl,j • Bj 2 , q = -(σ -) 1 -(b 0 + + b 0 -) , q = 0 (π) (18)
where Bl and Bj are the unit vectors corresponding to the MOT magnetic fields at, respectively, the atom l and atom j positions, and u = 1 + rl,j • Bl 2 is a normalization constant common to both b + q and b - q . The scattered intensity is given by

I S,q (r l , r j , v l ) = η q (r l , r j ) × P L,tot,q (r l , v l ) 4π|r l,j | 2 (19)
where P L,tot,q = α=x,y,z |F + α,q | + |F - α,q | c is the total power scattered by a single 2-level transition of atom l [refer to Eq. ( 3)], and η q is the corresponding normalized radiation pattern:

η q (r l , r j ) =            3 4 1 + rl,j • Bl 2 , q = + (σ + ) 3 4 1 + rl,j • Bl 2 , q = -(σ -) 3 2 1 -rl,j • Bl 2 , q = 0 (π) (20)
The rescattering cross section σ R,q is found by evaluating an overlap integral between (a) the emission spectrum that atom l produces when illuminated by the laser field and (b) the absorption spectrum of a single 2-level transition of atom j for the scattered field in the presence of the laser field. Denoting the emission spectrum of atom l by S tot,q and the absorption spectrum of a single 2-level transition of atom j by σ A,q , we have σ R,q (r l , r j , v l , v j ) = dω S tot,q (ω, r l , r j , v l )σ A,q (ω, r j , v j ) [START_REF] Hennequin | Stochastic dynamics of the magnetooptical trap[END_REF] where S tot,q (ω, r l , r j , v l ) = q =+,0,-w q q (r l , r j )S q (ω, r l , v l ) [START_REF] Moon | Nonlinear, nonequilibrium and collective dynamics in a periodically modulated cold atom system[END_REF] is composed of the normalized emission spectra S -, S 0 , S + for the respective σ -, π, σ + transitions of atom l, together with their spectral weight w q q (r l , r j ) = b q q (r l , r j )η q (r l , r j )

q =+,0,-b q q (r l , r j )η q (r l , r j )

where w - q + w 0 q + w + q = 1, such that S tot,q [Eq. ( 22)] is a normalized emission spectrum.

In Refs. [START_REF] Mollow | Power spectrum of light scattered by twolevel systems[END_REF] and [START_REF] Mollow | Stimulated emission and absorption near resonance for driven systems[END_REF], quite general expressions for the emission and absorption spectra for a 2-level atom have respectively been derived. Normalizing these expressions and adapting them to our approximate description, we obtain the expressions for S q (ω, r l , v l ) = S q (ω) and σ A,q (ω, r j , v j ) = σ A,q (ω) seen below. We have defined the total Rabi frequency for a single 2-level transition of the

F = 0 → F = 1 model atom, Ω tot,q (r, v) = Γ Itot,q(r,v) 2Isat
, and the corresponding detuning incorporating the Zeeman effect, ∆ q (r) = ∆ -µ q (r). Note that S q is written as a sum of two contributions. The first term, which involves a Dirac delta function, is an elastically scattered spectrum, which is centered at a single frequency specified by ∆ q . The second term is an inelastically scattered spectrum, which, on the other hand, is polychromatic. Because, in the Wieman model, the inelastic scattering is critical for the repulsion to win over the compression, we have been motivated to investigate the impact of elastically and inelastically scattered spectra on the instabilities, and in Sec. III we see the outcome. The impact of rescattering (as a whole) is also investigated.

Note that the Doppler effect can be included in ∆ q (r) in an alternative description, where each of the six MOT beams is separately scattered by the respective three atomic transitions. In such a case, 18 rescattering crosssections would have to be used, instead of the current 3 [Eq. [START_REF] Hennequin | Stochastic dynamics of the magnetooptical trap[END_REF]]. Nevertheless, neither of these descriptions is fully justified considering the complexity involved in modeling of the atom's behavior as it is coupled to the interacting field of several beams. In our simulations, we usually have k L |v α | Γ, |∆|, and thus the Doppler effect can be omitted. However, when deeply in the unstable regime, this is no longer true, and k L |v α | can become on the order of Γ.

S q (ω) = Γ 2 + 4∆ 2 q Γ 2 + 4∆ 2 q + 2Ω 2 tot,q δ(ω -∆ q ) + ΓΩ 2 tot,q 2π    (ω -∆ q ) 2 + 1 2 Ω 2 tot,q + Γ 2 Γ 2 1 2 Ω 2 tot,q + ∆ 2 q + 1 4 Γ 2 -2(ω -∆ q ) 2 2 + (ω -∆ q ) 2 Ω 2 tot,q + ∆ 2 q + 5 4 Γ 2 -(ω -∆ q ) 2 2    (24) σ A,q (ω) = σ0Γ 4 × Γ 2 + 4∆ 2 q Γ 2 + 4∆ 2 q + 2Ω 2 tot,q × (-iω + i∆ q + Γ) -iω + i2∆ q + Γ 2 + 1 2 iΩ 2 tot,q (ω -∆ q )/ i∆ q + Γ 2 (-iω + i∆ q + Γ) -iω + i2∆ q + Γ 2 -iω + Γ 2 + Ω 2 tot,q -iω + i∆ q + Γ 2 + c.c. ( 25
)
We finally discuss the similarities and differences between our kinetic model and that of Ref. [START_REF] Pohl | Self-driven nonlinear dynamics in magneto-optical traps[END_REF]. In both cases, the same main physical effects except for diffusion are used. In particular, both models include Doppler trapping force, attenuation, rescattering force, with atomic cross-sections possessing complete spatial dependence, thus resulting in all the forces being nonlocal. However, an important difference is that, unlike in Ref. [START_REF] Pohl | Self-driven nonlinear dynamics in magneto-optical traps[END_REF], we do not make any assumption of spherical symmetry of the forces. This allows us to observe, e.g., center-ofmass (COM) oscillations of the cloud, which are also seen in experiments [START_REF] Gaudesius | Phase diagram of spatiotemporal instabilities in a large magneto-optical trap[END_REF]. Moreover, we work with the more complex F = 0 → F = 1 system compared to the 2-level system in Ref. [START_REF] Pohl | Self-driven nonlinear dynamics in magneto-optical traps[END_REF]. This improvement allows us to properly describe, e.g., the anisotropy of the trapping force or, more generally, the features related to the magnetic field and light polarization.

B. Implementation

The collective system dynamics in the F = 0 → F = 1 model (Sec. II.A) are described by the following collisionless Vlasov-type kinetic equation for the atomic phasespace density f (r, v, t) [START_REF] Pohl | Self-driven nonlinear dynamics in magneto-optical traps[END_REF][START_REF] Romain | Phase-space description of the magneto-optical trap[END_REF]:

∂ ∂t f + v ∂ ∂r f + 1 M ∂ ∂v {[F tr (r, v) + F rsc,co (r, v)]f } - 1 M 2 ∂ 2 ∂v 2 [D(r, v)f ] = 0 ( 26
)
where M is the atomic mass, F tr (r, v) is the trapping force given by Eq. ( 1), D(r, v) is the momentum diffusion coefficient given by Eq. ( 7), and F rsc,co (r, v) = dr dv F sg rsc (r, r , v, v )f (r , v , t) is the continuous space analog of F rsc (r, v) in Eq. ( 15), where F sg rsc defines the rescattering force due to a single surrounding atom; the positions and velocities depend implicitly on time t. Given both the local and nonlocal spatial dependencies of the model's forces, finding direct numerical solution to Eq. ( 26) is highly demanding. As an equivalent yet more sensible way of solving the dynamics of our system, we use a numerical procedure based on a superparticle [START_REF] Birsdall | Plasma physics via computer simulation[END_REF] treatment similar to the one in Ref. [START_REF] Pohl | Self-driven nonlinear dynamics in magneto-optical traps[END_REF].

Note that a hydrodynamical description may not be directly comparable to our (kinetic) description. Hydrodynamics typically simplifies the phase space by restricting to, e.g., mean-field velocities. This is not fulfilled especially in the unstable regime. As a consequence, the applicability of such simplified descriptions (see, e.g., the MOT photon bubble model [START_REF] Mendonça | Photon bubbles in ultracold matter[END_REF]) to the present problem is not obvious.

Before getting into detail about the different elements involved in our implementation, we sketch first how it proceeds. At the initial time t 0 , a Gaussian cloud is generated, composed of superparticles, i.e., particles that represent collections of regular particles for increasing the simulation speed. Next, the beam intensity attenuation is evaluated at the superparticle positions, with help of our developed tube method [START_REF] Gaudesius | Self-oscillating clouds in magneto-optical traps[END_REF]. Forces acting on each superparticle are then computed, and, using these forces, the Leapfrog algorithm [START_REF] Kaya | The Leapfrog algorithm: Background and demonstration[END_REF] propagates the cloud in time by one time-step δt. At the new positions, the attenuation is first evaluated, then the forces are computed, and the algorithm propagates the cloud by one δt -this cycle is repeated until the simulation end time t end .

Spatiotemporal instabilities can arise in the simulations, and in Refs. [START_REF] Gaudesius | Instability threshold in a large balanced magnetooptical trap[END_REF][START_REF] Gaudesius | Phase diagram of spatiotemporal instabilities in a large magneto-optical trap[END_REF] we have demonstrated that we are successful in reproducing experimental behaviors. The 3D nature of our simulations is showcased in Fig. 4 (see also online Supplemental Material [38]), displaying clouds belonging to different regimes obtained also experimentally (see Ref. [START_REF] Gaudesius | Phase diagram of spatiotemporal instabilities in a large magneto-optical trap[END_REF] for the discussion). Note that in Ref. [START_REF] Pohl | Self-driven nonlinear dynamics in magneto-optical traps[END_REF], their kinetic model has been used to achieve simulations of balanced MOT instabilities in quasi-1D.

The exact iterative scheme we employ in the Leapfrog algorithm for updating the superparticle velocity and position in time is as follows:

1. Compute total force F tot (n) = F tot (r(t n ), v(t n-1/2 )). 2. Update velocity v(t n+1/2 ) = v(t n-1/2 ) + Ftot(n) Msup δt. 3. Update position r(t n+1 ) = r(t n ) + v(t n+1/2 )δt.
Here the integer n = 0, 1, 2, . . . marks the iteration step of fixed size δt, with n = 0 used when initiating the scheme, and M sup = εM is the superparticle mass, with ε = N/N sup being a fixed superparticle number N sup scaling (i.e., the number of regular particles represented by one superparticle) that we note also scales σ 0 , so one correspondingly has σ sup,0 = εσ 0 . We have in Refs. [START_REF] Gaudesius | Instability threshold in a large balanced magnetooptical trap[END_REF][START_REF] Gaudesius | Phase diagram of spatiotemporal instabilities in a large magneto-optical trap[END_REF] used N sup = 7 × 10 3 -a choice that was based on tests for the simulation outcome to be independent of N sup ; we found that below N sup = 10 3 the clouds are stable. Observe, in the iterative scheme, that the velocity and position are updated in an interleaved manner, i.e., they leapfrog over each other. When the scheme is initiated (n = 0), the velocity v(t -1/2 ) and the position r(t 0 ) are specified. To obtain v(t -1/2 ), we utilize the approximation v(t -1/2 ) ≈ v(t 0 ), which is assumed to be valid for a small enough δt. In Refs. [START_REF] Gaudesius | Instability threshold in a large balanced magnetooptical trap[END_REF][START_REF] Gaudesius | Phase diagram of spatiotemporal instabilities in a large magneto-optical trap[END_REF], the x, y, z components of v(t 0 ) have been picked to be random between 0 and 0.01 m/s [much less than the root-mean-square (RMS) The total force acting on a superparticle with the position r s = r(t n ) and the velocity

v s = v(t n-1/2 ) is F tot (r s , v s ) = εF tr (r s , v s ) + ε 2 F rsc (r s , v s ) + εF D (r s , v s ) ( 27
)
where F D is the stochastic force given by [START_REF] Sachs | Langevin dynamics with variable coefficients and nonconservative forces: From stationary states to numerical methods[END_REF]

F D (r s , v s ) = 2D(r s , v s ) 3δt × g(n)r s ( 28
)
where D is the momentum diffusion coefficient given by Eq. ( 7), and g(n) is a Gaussian white noise with the autocorrelation function equal to δ(n). We note that, in addition to F D , the velocity dependence in our total force makes the Leapfrog algorithm no longer time reversible. In Eq. ( 27), F tr is scaled by ε as it is written as a sum of terms that contains σ ± α,q ∝ σ 0 [see Eq. ( 3)]. F rsc is scaled by ε 2 as it is written as a sum of terms that contains a product of σ ± α,q ∝ σ 0 and σ R,q ∝ σ 0 [see Eqs. [START_REF] Terças | Driven collective instabilities in magneto-optical traps: A fluiddynamical approach[END_REF][START_REF] Felinto | Influence of the number of atoms in a ring-shaped magneto-optical trap: Observation of bifurcation[END_REF][START_REF] Hennequin | Stochastic dynamics of the magnetooptical trap[END_REF][START_REF] Hanley | Quantitative simulation of a magneto-optical trap operating near the photon recoil limit[END_REF]]. F D is scaled by ε for the following reasons. As the diffusion coefficient D describes the equilibrium between the diffusive heating and Doppler cooling processes, one can write D ∝ γT lim [START_REF] Cohen-Tannoudji | Atomic motion in laser light[END_REF], where the trap friction constant γ is scaled by ε as F tr is scaled likewise, and the limit temperature T lim is scaled by ε as T lim ∝ M according to the equipartition theorem. Moreover, F D involves a square root of D [see Eq. ( 28)]. Taking everything into account, F D is thus scaled by ε.

For evaluating the beam intensity attenuation at the superparticle positions, we employ our developed tube method, whose name is attributed to the fact the attenuation of each beam is calculated in rectangular tubesegments parallel to the beam's propagation direction. This method is implemented by extending its 2D illustration detailed in Fig. 5. In particular, we numerically generate a fixed grid of points in 3D space and calculate the intensity of each beam at the positions of the grid-points that contain the superparticles, after which the intensity at each superparticle position is found by means of interpolation. The interpolated intensity values are then used in the calculation of the forces [see Eq. ( 27)]. Finding a beam intensity at a given gridpoint position r g = (x g , y g , z g ) involves the assumption that a given superparticle at r s = (x s , y s , z s ) is represented by a Dirac delta function δ(r g -r s ), allowing us to write the cloud density ρ(r g ) = s δ(r g -r s ) = s δ(x g -x s )δ(y g -y s )δ(z g -z s ). With this assumption, the intensity for, e.g., the positive ẑ directed beam [see Eqs. [START_REF] Romain | Spatial instabilities in a cloud of cold atoms[END_REF][START_REF] Kim | Observation of the Hopf bifurcation in parametrically driven trapped atoms[END_REF]] is numerically determined from

I + z (xg, yg, zg) = I∞e -O + z (xg,yg,zg) , O + z (xg, yg, zg) ≈ ε W 2 × zs<zg |xs-xg|<W/2 |ys-yg|<W/2
q=-,0,+ p + z,q (xs, ys, zs)σ + z,q (xs, ys, zs, vs(xs, ys, zs))

where the superparticle scaling ε is taken into account, W is a fixed tube-width, and W 2 is the corresponding transverse tube-area. From this equation, we observe that the only superparticles that contribute at a given grid-point position are those inside the tube of this point (|x s -x g |, |y s -y g | < W/2) and positioned before it (z s < z g ). The remaining beam intensities are found in an analogous way. The value of W is picked according to our tests' results with uniform and Gaussian clouds, where numerically determined attenuation profiles were checked to converge with corresponding analytically calculated profiles. In the simulations of Refs. [START_REF] Gaudesius | Instability threshold in a large balanced magnetooptical trap[END_REF][START_REF] Gaudesius | Phase diagram of spatiotemporal instabilities in a large magneto-optical trap[END_REF], we have used W = 0.15σ, where σ is the RMS width of the initial Gaussian cloud (t 0 ), picked to be as close as possible to the RMS radius of the cloud after transient behavior (Fig. 3 in Ref. [START_REF] Gaudesius | Instability threshold in a large balanced magnetooptical trap[END_REF] exemplifies this behavior). We note that the choice of the value of W is a compromise. Indeed, if W is too small compared to σ, too few superparticles will participate in the determination of the attenuation, thus yielding large spatial fluctuations. In the opposite limit where W is larger than σ, the spatial dependence of the attenuation will be washed out. In this case, W is expected to limit the size of the structures that can appear in a simulated cloud (such as the ones seen in Fig. 4).

Lastly, we explain the implementation of the beam cross-saturation effect, which we recall from Sec. II.A appears because I tot,q enters into scattering cross sections [see Eqs. [START_REF] Guerin | Subradiance in a large cloud of cold atoms[END_REF][START_REF] Labeyrie | Self-sustained oscillations in a large magneto-optical trap[END_REF]]. To implement this effect, the beam intensities at the grid-point positions are calculated selfconsistently. To initiate this calculation, one must specify initial beam intensity values (discussed below), which are used to determine I tot,q . With this I tot,q , new beam intensities can be found [as in Eq. ( 29

III. Impact of different effects on the instabilities

The simulations offer a better understanding of the complex instability mechanism by varying the magnitude of the different effects. In the following, we illustrate this with preliminary tests on the impact of diffusion, attenuation, and rescattering, as well as of elastic versus inelastic scattering. In doing so, we look at how the detuning at threshold, ∆ thr , is affected. Indeed, the detuning is the most sensitive parameter for determining the state (stable or unstable) of a MOT [START_REF] Labeyrie | Self-sustained oscillations in a large magneto-optical trap[END_REF].

Let us begin by investigating how ∆ thr is affected as the stochastic force F D [Eq. [START_REF] Cohen-Tannoudji | Atomic motion in laser light[END_REF]] is scaled by a constant factor d. We use d = 0, 0.5, 2, 5 and concentrate on the case with B = 3 G/cm. We display, in Fig. 6, the outcome of this investigation. As can be seen, with no diffu- sion (d = 0), ∆ thr is the same as in the ordinary simulation (d = 1; -∆ thr /Γ = 2.99). Importantly, because the instabilities can exist without the diffusion, this tells us that it is not essential for triggering them. We note that such a conclusion stands in agreement with the results of Ref. [START_REF] Pohl | Self-driven nonlinear dynamics in magneto-optical traps[END_REF], where diffusion was neglected. Our studied instabilities are thus different from retroreflected MOT instabilities of the stochastic type, where noise is necessary [START_REF] Hennequin | Stochastic dynamics of the magnetooptical trap[END_REF]. With an increased diffusion (d > 1), ∆ thr is seen to shift closer to the resonance, with the change being ∼0.4 Γ per half a decade. A smaller -∆ thr means a reduced instability range and, therefore, a weakened instability mechanism. This may not be surprising conside-ring that, in a feedback system, the phase of the feedback is critical in determining whether the system is stable or unstable [START_REF] Åström | Feedback systems: An introduction for scientists and engineers[END_REF]; increasing the diffusion may be regarded as destroying the phase relationship, thus preventing the unstable regime from being entered. In conclusion, the diffusion is not an essential effect but, otherwise, suppresses the instability mechanism if it becomes great enough.

Next, we continue with our investigation on whether the instabilities persist when either attenuation is removed (i.e., the beam intensity is constant) or rescattering is removed, for B = 3 G/cm at different ∆ values in the previous range of simulated unstable and stable clouds [START_REF] Gaudesius | Phase diagram of spatiotemporal instabilities in a large magneto-optical trap[END_REF]. Note that while the former case is experimentally relevant (for relatively low atom number and/or weak magnetic field gradient), the latter one is not, as one cannot have large attenuation without multiple scattering. Without attenuation, we find for all the explored parameter range that the MOT is stable. When attenuation is present but rescattering is turned off, we observe small clouds with temporal fluctuations in their COM positions, but relatively stable RMS radii. The COM fluctuations are in this particular case attributed to diffusion. These tests thus seem to show that both attenuation and rescattering are necessary to reach the unstable regime. This is consistent with previous models [START_REF] Labeyrie | Self-sustained oscillations in a large magneto-optical trap[END_REF][START_REF] Pohl | Self-driven nonlinear dynamics in magneto-optical traps[END_REF] and supports the view that the shadow force (produced by attenuation) takes part in a feedback mechanism, where this force works against the cloud expansion due to the rescattering force in order to produce unstable motion in a balanced MOT. The beam attenuation alone can be noted to be critical for instabilities in a retroreflected MOT [START_REF] Hennequin | Stochastic dynamics of the magnetooptical trap[END_REF][START_REF] Di Stefano | Deterministic instabilities in the magneto-optical trap[END_REF] as well as in the case of collective phenomena that parameter-modulated MOT instabilities can exhibit [START_REF] Moon | Nonlinear, nonequilibrium and collective dynamics in a periodically modulated cold atom system[END_REF]. Moreover, Ref. [START_REF] Gaudesius | Self-oscillating clouds in magneto-optical traps[END_REF] has recently identified both attenuation and rescattering as necessary effects for a spatiotemporal instability in a misaligned MOT, using a slightly modified version of our simulation model.

Finally, we investigate the impact of elastically versus inelastically scattered light on the instabilities. Here we use the notations σ el,q and σ inel,q in denoting the parts of σ R,q [Eq. ( 21), with q = q] that result from the contribution of, respectively, elastically and inelastically scattered spectra; they satisfy σ R,q = σ el,q + σ inel,q . In Fig. 7, we display how ∆ thr is affected for B = 3 G/cm after removing either part. Importantly, the fact the instabilities still are obtained indicates that none of these parts alone is necessary. This is surprising considering the (beforementioned) Wieman model prediction that the cloud expansion relies critically on the presence of inelastic scattering in the cloud. The thresholds are seen to be shifted further away from the resonance compared to the ordinary simulation result (-∆ thr /Γ = 2.99), by ∼2.5 times (-∆ thr /Γ = 7.22 for σ el,q = 0 and -∆ thr /Γ = 8.06 for σ inel,q = 0). The shift to larger absolute values is correlated with the fact that we observe the simulated clouds to become smaller (by ∼1.5 times, with the σ inel,q = 0 case being slightly smaller in cloud size), which is consistent with the Wieman model (there is less rescattering). With this decrease, one is led to an increase in the optical depth and thus the shadow force. On the opposite hand, we find that artificially increasing σ R,q , and consequently the size, shifts the threshold closer to the resonance. This result makes perfect sense considering the second investigation's finding (on attenuation), as for a vanishing shadow force the instabilities should disappear.

!" thr =! < el;q = 0 , < inel;q 6 = 0 < el;q 6 = 0 , < inel;q 6 = 0 < el;q 6 = 0 , < inel;q = 0 Figure 7: Investigation on how the threshold detuning ∆ thr for the magnetic field gradient B = 3 G/cm is affected after removing either σ el,q or σ inel,q , being the parts of rescattering cross section that result from the contribution of, respectively, the elastically and inelastically scattered spectra. The small asterisks are the test results, and the large asterisk is the result of the ordinary simulation (from Fig. 6 in Ref. [START_REF] Gaudesius | Instability threshold in a large balanced magnetooptical trap[END_REF]).

To summarize, both attenuation and rescattering seem to be necessary for generating instabilities in a balanced MOT, and their mechanism is strengthened when the shadow force (produced by attenuation) gets larger compared to the rescattering force. Also, great enough diffusion suppresses this mechanism, and neither elastically nor inelastically scattered light alone is critical for the generation.

IV. Conclusion

In this paper, we presented a numerical tool for studying MOT instabilities in a full-blown 3D environment. It has been successfully employed to predict instability thresholds [START_REF] Gaudesius | Instability threshold in a large balanced magnetooptical trap[END_REF] and various unstable regimes [START_REF] Gaudesius | Phase diagram of spatiotemporal instabilities in a large magneto-optical trap[END_REF] in a large balanced MOT. These simulations can be used in the future to investigate features that are challenging to access experimentally, such as, e.g., velocity fields and 3D density distributions. Other lines of research include the analysis of the cloud dynamics in terms of turbulence, to compare with recent experimental observations [START_REF] Griffin | Turbulence in a magneto-optical trap[END_REF]. With some minor modifications, our model is applicable to different MOT configurations, as recently exemplified for a misaligned MOT [START_REF] Gaudesius | Self-oscillating clouds in magneto-optical traps[END_REF]. In the case of a retroreflected MOT, the added value of the spatial information could be strongly beneficial for studying the observed spatiotemporal instabilities [START_REF] Romain | Spatial instabilities in a cloud of cold atoms[END_REF]. Improving the understanding of MOT instabilities can effectively continue through investigations on how these are impacted by simulation model's physical effects, and broader perspectives can be opened up by refining the descriptions of currently included effects and/or by adding new ones (e.g., dipole forces, higher-order rescattering).

Figure 2 :

 2 Figure 2: (a) Depiction of the mechanism behind the shadow force. As the oppositely directed beams travel through the cloud, their respective intensities I +α and I - α (α = x, y, z) become attenuated, giving rise to an intensity imbalance in the cloud. This imbalance yields an additional compression, i.e., the shadow force. (b) Depiction of the mechanism behind the rescattering force with two 2-level atoms. Atom 1 first scatters the laser light of intensity I L , and atom 2 then rescatters the scattered light, thus experiencing a repulsion. The rescattered light has the intensity I R determined by the inverse-square law of propagation, such that this repulsion is Coulomb-like. The scattering cross section σ L and the rescattering cross section σ R determine respectively the power scattered and rescattered.

Figure 3 :

 3 Figure 3: Illustration of the rescattering in the F = 0 → F = 1 model. Each of the σ -, π, σ + transitions of atom 1 produces a characteristic radiation pattern of light scattered from the total laser field (respectively, η -, η 0 , η + ). The scattered light polarization seen by atom 2 depends on the quantization axis orientations of the two atoms, such that each transition of atom 2 rescatters light in specific fractions (respectively, b - -, b - 0 , b - + , when considering light from the σ -transition of atom 1).

igure 4 : 2 L Γ 3 ∆(Γ 2

 4232 Display of 3D images (upper row) and corresponding 2D images (lower row) of simulated clouds belonging to different instability regimes discussed in Ref.[START_REF] Gaudesius | Phase diagram of spatiotemporal instabilities in a large magneto-optical trap[END_REF]. Each image is a single-shot image. In the 3D images, the individual dots are superparticles. The 2D images are the cloud densities integrated along the z axis and are Gaussian filtered; the diagonals correspond to the directions of two pairs of MOT beams; the field of view is 10 × 10 cm 2 for all these images. The x and y axes have been rotated by 45 • between the 3D and 2D images. A video version of this figure is available as online Supplemental Material [38]. velocity of Rb-87 at the Doppler temperature]. To pick δt, a rule of thumb is employed, telling that for a given set of MOT parameters we use δt < 0.1/ω tr , where ω tr = κ/M is the trap frequency, with κ = +6Γ 2 I∞/Isat+4∆ 2 ) 2 being the trap spring constant [found from the linear expansion of Eq. (1) for low velocities and positions near the trap center, with fixed intensity I ∞ ]. In Refs. [23, 24], we have used the following constants (related to Rb-87 and its D2 line): M = 1.443 × 10 -25 kg, µ = 2π × 1.4 × 10 6 Hz/G, Γ = 2π × 6.07 MHz, k L = 2π 780×10 -9 m -1 , and I sat = 1.67 mW/cm 2 .

  )]. These new intensities are next used to construct intensities that are equal mixtures of new and last intensities, i.e., we construct [New intensities + Last intensities]×1/2. The constructed intensities are compared to the last intensities, and this is reiterated until convergent intensities are found. (These convergent intensities are then used in the interpolation of the intensities at the superparticle positions, after which the forces acting on the superparticles are found.) Regarding initial beam intensity values, these are picked to be equal to I ∞ ; any values can, in principle, work, due to the convergent intensities being independent of such a choice. At the second Leapfrog algorithm iteration step, the intensity calculation starts with the convergent beam intensities found in the previous iteration step, and so on for the remaining iteration steps (this increases the simulation speed).

Figure 5 :

 5 Figure 5: (a) First and (b) second part of a 2D illustration of the numerical method called the tube method. (a)The upper two drawings display the same cloud composed of superparticles (small dots), and it is imagined that each beam is segmented into rectangular tubes parallel to the beam's propagation direction. During the cloud's evolution, the tubes remain at fixed positions, with their width W being of a set size. The lower drawing displays the same cloud, with the added large dots indicating the positions of the points where the intensity of each beam is calculated first. These dots are placed at grid-point positions located through the center of the tubes. (b) In the calculation of, e.g., the +ẑ directed beam's intensity at the position of a given large dot, the superparticles that contribute are those inside the tube of this dot and positioned before it. Once the intensities of each beam at the positions of the large dots are calculated, the intensities at the positions of the superparticles are found by means of interpolation.

Figure 6 :

 6 Figure6: Investigation on how the threshold detuning ∆ thr for the magnetic field gradient B = 3 G/cm is affected as the stochastic force F D [Eq.[START_REF] Cohen-Tannoudji | Atomic motion in laser light[END_REF]] is scaled by a constant factor d. The open stars are the test results, and the solid star is the result of the ordinary simulation (from Fig.6in Ref.[START_REF] Gaudesius | Instability threshold in a large balanced magnetooptical trap[END_REF]).

Acknowledgments

Part of this work was performed in the framework of the European training network ColOpt, which is funded by the European Union (EU) Horizon 2020 programme under the Marie Sklodowska-Curie action (Grant No. 721465), and of the European project ANDLICA, ERC advanced grant agreement (Grant No. 832219). It has also been supported by the Danish National Research Foundation through a Niels Bohr Professorship to T.P. and through the Center of Excellence "CCQ" (Grant No. DNRF156).

Data availability:

The data supporting this work are available from the corresponding author upon reasonable request.