F D'acapito
email: dacapito@esrf.fr

C Maurizio

M C Paul

Th S Lee

W Blanc

B Dussardier

Role of CaO addition in the local order around Erbium in SiO2-GeO2-P2O5 fiber preforms

Keywords: Optic fibers, EXAFS, Erbium, Nanoparticles

The development of materials for optical signal processing represents a major issue in present technology. In this contribution we present a study on Er-doped fiber preforms where particular attention is devoted on how the addition of CaO in the glass modifies the local environment of the rare earth. The results from photoluminescence and Extended X-ray Absorption Fine Structure (EXAFS) are compared and a clear link between the width of the emission line at 1.5 µm and the amorphous/crystalline local structure around the Er3+ ion is evidenced.

1. Introduction Knowledge of the environment of rare earth ions in glasses is essential for the understanding of the spectral properties of these glasses and the design of new glass compositions for photonics applications. Among the rare earths, erbium is especially interesting because the 4I13/2 -4I15/2 transition at ≈1.54 µm coincides with the lowest attenuation window of silica glass fiber; hence, erbium doped glasses can be used to amplify attenuated signals in optical fiber telecommunication systems. While Erbium-Doped Fiber Amplifier (EDFA) was developed 20 years ago, extensive studies are still carried out to improve its properties. In particular, linear dimensions should be reduced although the severe limitations imposed by the poor rare earth solubility into silica matrix. Moreover, spectral bandwidth should be increased to improve the Wavelength Division Multiplexing (WDM) applications. To solve these problems, various approaches have been proposed namely Er-Tm co-doping [1] or Er incorporation in nanoparticles (NP). However, to be compatible with optical waveguide applications, Rayleigh scattering must be kept at a minimum. For this purpose, Tick proposed that the particles size must be less than 15 nm with a narrow distribution size and high density [2]. There are reports on thin films containing semiconductor and metal nanoparticles acting as sensitizers, absorbing the incident light with a high crosssection and exciting erbium co-doped ions in the nanoparticle vicinity through energy transfer [3][4][5][6][7][8]. Here we propose for the first time to incorporate Er into oxide nanoparticles in optical fibers to improve the spectroscopic properties of erbium ions. To our knowledge only few studies on nano-structured silica fibers were dedicated to metal ions properties [9]. In our samples, oxide nano-structures are prepared by adjusting the composition of usual modifiers in the core of the fiber. A glass with an immiscibility gap is obtained and two phases are then formed with high and low silica content. The nanoparticles are preserved when the perform is drawn into a fiber. An effective method to study the incorporation of Er in a matrix is by Extended X-ray Absorption Fine Structure (EXAFS) [10] as this technique provides quantitative information on the geometry and chemical composition of the environment of the element under study. The Er site has been successfully studied by this technique in a variety of matrices like crystalline Si [11,12] amorphous Si [13], Si nanoparticles [14], and silicate glass [15][16][17]. In the present study we have used EXAFS to describe the site of Er and to establish the crystalline/amorphous nature of its environment.

Results

2.1. Preform preparation and characterization Fiber preforms were prepared by Modified Chemical Vapor Deposition (MCVD) associated to solution doping. Samples of composition SiO2-GeO2-P2O5-Er2O3 with and without CaO were investigated to evidence the role of Ca in the incorporation site for Er in the glass. The multiplicity of RE sites in the host matrix are known to enhance the inhomogeneous broadening of the emission/absorption lines [18] and P-codoped silicates are reported to bind in an ordered way to Er3+ ions [19]. When preforms were drawn into fibers, NP were not directly evidenced into fibers through electron microscopy analysis; however, high attenuation losses were measured, presumably due to Rayleigh scattering induced by NP. Further analyses on this point are in progress. The preforms were cut in disks of 5 mm in diameter where the Er-doped zone is limited to a central core of 1 mm in diameter. Two samples were prepared with an Er2O3 content of 0.1 mol% and ≈1 mol% of P2O5 . Sample A has 0.2 mol% of CaO whereas sample B contains no CaO. Slices of samples A and B were cut from preforms (thickness 1 mm), then successively thinned by mechanical means down to ≈10 µm and by ionic beam down to few tens of nanometers. They were characterized by electron microscopies such as Transmission Electron Microscopy (TEM) and Energy Dispersive X-Ray Micro-analysis (EDX) to analyze the chemical contrast and map the chemical composition. Fig. 1 shows a TEM image of the sample A, where a polydisperse assembly of spherical nanoparticles with a mean diameter of 50 nm can be seen. Smaller particles of 10 nm are visible and the size of the biggest ones was ≈200 nm (not shown). The formation of these particles is linked to the addition of CaO, as presented elsewhere [START_REF] Dussardier | Actes des 25 ernes journées nationales de l'optique guidée[END_REF]. High content of modifiers elements ([Ca] and [P] ≈30%) are present in these nano-structures. Moreover, erbium ions are found to be located only inside the nanoparticles. Spectroscopic characterizations at room temperature on the emission line associated to the 4I13/2 → 4I15/2 transition at 1.54 µm were made on these samples at LPCML (Lyon, France) [START_REF] Peretti | Actes des 26 emes journées nationales de l'optique e guidée[END_REF]. The results are shown in Fig. 2 where we evidence the fact that the emission spectrum of samples A is broader than that of sample B. 2.2. EXAFS measurements EXAFS measurements at the Er-LIII edge (E = 8358 eV) were carried out at the GILDA-CRG beamline [START_REF] Acapito | [END_REF] at the European Synchrotron Radiation Facility. The monochromator was equipped with a pair of Si(3 1 1) crystals and was run in dynamical focusing mode [23]. The rejection of higher order harmonics was achieved by using a pair of Pd-coated mirrors with a cutoff energy of 20 keV. The second mirror was also used as focusing element in the vertical direction and the setup was optimized to obtain the minimum spot size on the sample; in this way 130 µm × 90 µm Full Width at Half Maximum (FWHM) were obtained. The absorption coefficient was measured in fluorescence mode by selecting the Er-La emission line using an high purity Ge detector. The position of each sample respect to the X-ray beam was optimized by carrying out horizontal and vertical scans and the Er doped zone resulted to be 500 µm FWHM wide. Spectra were collected at room temperature. EXAFS data are shown in Fig. 3 whereas the related Fourier Transforms are shown in Fig. 4. Sample A exhibits a single frequency signal whereas sample B presents higher frequency oscillations that are in good agreement with those observed in the spectrum of ErPO4. Data were analyzed using the GNXAS [24,25] and Feff [26,27] codes. The theoretical paths were calculated using the muffin-tin approximation for the scattering potential and the Hedin-Lundquist form for its energy-dependent part. Atomic clusters were derived from the Er2Si2O7 [28] and ErPO4 [29] crystal structures. Sample A was quantitatively analyzed by using a model consisting in an Er3+ ion linked to SiO4 tetrahedra as described in [12,17]. Here we note that, due to the similar backscattering amplitude and phase of Si and P it is not possible to distinguish these two atoms in the second coordination shell. Sample B was reproduced with the ErPO4 structure up to the 5th shell. The results are collected in Table 1.

Discussion

In sample A the rare earth is linked to O atoms in the first coordination shell and to Si(P) atoms in the second shell in a way similar to that already observed in silicate glasses [17,30], phosphate glasses [31,[START_REF] Acapito | Proceedings of the xafs13 Conference, AIP Conf. Proc[END_REF] and in the oxide formed at purpose in crystalline Si [12]. Also the structural parameters (about seven O atoms at 2.26 A and Si(P) atoms at 3.6 A corresponding to a Er-O-Si(P) bond angle of ≈140°) are in good agreement with the cited literature. Si(P) atoms are visible as they belong to the same Si(P)O4 tetrahedron as the first shell O atoms but no further coordination shells are detected. This permits to state that an amorphous environment is realized around Er. On the other hand sample B presents a completely different EXAFS signal that is well comparable with the spectrum of ErPO4. This means that Er in this case is inserted in a locally well ordered phase of about a few coordination shells (around 4-5 A around the absorber). The fact that TEM on this sample reveals a uniform sample is not in contradiction with this result; it just means that this phase is not spatially extended to form nm-sized nanoparticles (in our TEM analyses the spatial resolution limited to few nm) but the ordering is extremely local, i.e. it is limited to only a few shells around the Rare Earth ion. From these considerations, the broadening of the emission spectrum observed in sample A can be attributed to an inhomogeneous broadening due to erbium located in a more disordered environment compare to sample B. Here we see that the cumulated effects of Ca and P within the Er-doped nanoparticles both amorphize the material structure around Er3+ ions and increase the fluorescence inhomogeneous broadening. Further investigations will be necessary to understand the combined role of these two ions in the amorphization of the local structure.

Conclusion

In this contribution we have investigated the structure around Er3+ ions in fiber preforms of composition SiO2 -GeO2-P2O5-Er2O3 with and without CaO. The addition of this component to the glass leads to the formation of calco-phospho-silicate nanoparticles. In the case without CaO the environment of Er is locally crystalline and reproduces well that observed in ErPO4 . Upon addition of Ca, amorphous nanoparticles are formed. While they contain Phosphorous, the presence of Calcium leads to an amorphization of the Erbium environment, like that observed in silicate or phosphate glasses. Correspondingly, the width of the luminescence line at 1.54 µm increases due to the inhomogeneous broadening introduced by the amorphous environment. This feature is particularly interesting in the production of materials for WDM applications. Fig. 1. TEM image from sample A (with Ca). The average particle size is 50 nm. Smaller (10 nm) particles are visible. Few bigger (≈200 nm) particles were observed, not shown on this image. Fig. 2. Room temperature optical luminescence spectra from samples A (with Ca, continuous line) and B (without Ca, dashed line), respectively. Samples were excited at 980 nm with a laser diode the system spectral resolution was better than 1 nm. The fluorescence was measured at right angle from excitation beam. The suffix a is relative to O whereas b is phosphorous in the case of samples B and ErPO 4 and it could be either P or Si for sample A; suffix c is Er. The number of first and second neighbors for sample A is 6.5 ± 0.2 whereas it is fixed to the crystallographic values of ErPO 4 (N1 = 4O, N2 = 4O, N3 = 2P, N4 = 4Er, N5 = 4P) in the other samples. Errors on the last digit are indicated in brackets.

Fig. 3 .

 3 Fig. 3. EXAFS spectra of the samples.

Fig. 4 .

 4 Fig. 4. Fourier Transforms of the EXAFS spectra. Transforms were carried out ˚ in the interval k ⊆ (2.5,. . ., 8.0) A -1 with a k 2 weight. The best fit curves are presented in dots.

Table 1

 1 Results of the quantitative EXAFS analysis: bond lengths R and Debye-Waller factors σ 2 for the first three shells around Er

Acknowledgements

This work was realized in the frame of a P2R project (Programme de Recherche en Réseau, CNRS (France) and DST (India)). We acknowledge the technical assistance of Fabrizio La Manna and Fabio d'Anca on the GILDA beamline; and Bernard Jacquier, Anne-Marie Jurdyc and Romain Peretti from LPCML (Lyon, France) for the spectroscopic measurements. Post-doctoral positions of M.C. Paul and Th. Lee Sebastian were funded by DST (Indian Department of Science and Technology, Boycast fellowship) and the MENRT (French Ministry of Research), respectively. GILDA is a project jointly financed by CNR and INFN institutes.


domain: math info

text:
Pierrick Dartois

Luciano Maino

Giacomo Pope

Damien Robert

An Algorithmic Approach to (2, 2)-isogenies in the Theta Model and Applications to Isogeny-based Cryptography

In this paper, we describe an algorithm to compute chains of (2, 2)isogenies between products of elliptic curves in the theta model. The description of the algorithm is split into various subroutines to allow for a precise field operation counting. We present a constant time implementation of our algorithm in Rust and an alternative implementation in SageMath. Our work in SageMath runs ten times faster than a comparable implementation of an isogeny chain using the Richelot correspondence. The Rust implementation runs up to forty times faster than the equivalent isogeny in SageMath and has been designed to be portable for future research in higher-dimensional isogeny-based cryptography.

Introduction

The devastating attacks on SIDH [START_REF] Castryck | An efficient key recovery attack on SIDH[END_REF][START_REF] Maino | A direct key recovery attack on SIDH[END_REF][START_REF] Robert | Breaking SIDH in polynomial time[END_REF] have highlighted the relevance of studying higher-dimension abelian varieties in isogeny-based cryptography. Following the attacks, it soon became evident that these new tools would have had applications beyond cryptanalysis. For instance, Robert leveraged these techniques both to give a representation of isogenies in polylogarithmic time [32] and to compute the endomorphism ring of ordinary elliptic curves in quantum polynomial time [START_REF] Robert | Some applications of higher dimensional isogenies to elliptic curves (overview of results)[END_REF].

On a more cryptographic side, the attacks have been used to design new protocols. Basso, Maino and Pope utilise these cryptanalytic techniques to construct a trapdoor mechanism, and using standard transformations, this trapdoor is used to derive a public-key encryption protocol named FESTA [START_REF] Basso | FESTA: Fast Encryption from Supersingular Torsion Attacks[END_REF]. Subsequently, two additional protocols employing similar ideas to FESTA have appeared [START_REF] Moriya | IS-CUBE: An isogeny-based compact KEM using a boxed SIDH diagram[END_REF][START_REF] Nakagawa | QFESTA: Efficient Algorithms and Parameters for FESTA using Quaternion Algebras[END_REF]. One of the main building blocks underlying these protocols is the computation of chains of (2, 2)isogenies between products of two elliptic curves. However, the cryptographic application of these isogeny chains extends beyond FESTA-based applications. For instance, they are at the core of computing the group action in SCALLOP-HD [START_REF] Chen | SCALLOP-HD: group action from 2-dimensional isogenies[END_REF], as well as in novel constructions of isogeny-based weak verifiable delay functions [START_REF] Decru | Towards a Quantum-Resistant Weak Verifiable Delay Function[END_REF] and verifiable random functions [START_REF] Leroux | Verifiable random function from the Deuring correspondence and higher dimensional isogenies[END_REF]. Therefore, improving algorithms to compute chains of (2, 2)-isogenies between elliptic products is of paramount importance to the progress of higher-dimensional isogeny-based protocols.

Prior to this work, the only method to compute (2, 2)-isogenies between elliptic products relied on ad-hoc procedures for gluing and splitting, and the use of the Richelot correspondence to compute isogenies between Jacobians of genus two hyperelliptic curves [START_REF] Smith | Explicit endomorphisms and correspondences[END_REF][START_REF] Oudompheng | A Note on Reimplementing the Castryck-Decru Attack and Lessons Learned for SageMath[END_REF]. This method can be considered satisfactory for cryptanalytic purposes, but it is definitely not efficient enough for constructive applications. Indeed, for the proof-of-concept implementations of [START_REF] Basso | FESTA: Fast Encryption from Supersingular Torsion Attacks[END_REF][START_REF] Nakagawa | QFESTA: Efficient Algorithms and Parameters for FESTA using Quaternion Algebras[END_REF], the two-dimensional isogenies are the bottleneck of the protocol. Richelot isogenies describe [START_REF] Basso | FESTA: Fast Encryption from Supersingular Torsion Attacks[END_REF][START_REF] Basso | FESTA: Fast Encryption from Supersingular Torsion Attacks[END_REF]-isogenies between Jacobians of genus two hyperelliptic curves in the Mumford model. Here, kernel elements are divisors, represented by a pair of univariate polynomials. The arithmetic of the group elements, as well as isogeny codomain computation and evaluation, require working in a univariate polynomial ring above the base field. This model makes doubling and evaluation of points expensive and the implementation of the isogeny chain itself is significantly more complicated than the more familiar isogeny chains between two elliptic curves, which use Vélu's formulae. A natural question is then to ask whether it could be possible to use different models that are more amenable to simple and efficient algorithmic descriptions.

In the literature, another model used to compute isogenies is already known: the theta model. Despite being suitable for isogenies between elliptic curves, the theta model has mainly been employed to compute isogenies in higher dimension due to the lack of alternatives. For instance, Cosset and Robert describe an algorithm for (ℓ, ℓ)-isogenies in the theta model for odd primes ℓ [START_REF] Cosset | Computing (ℓ, ℓ)-isogenies in polynomial time on Jacobians of genus 2 curves[END_REF], later improved in [START_REF] Lubicz | Fast change of level and applications to isogenies[END_REF]. The case ℓ = 2 has been briefly treated in [START_REF] Robert | Fonctions thêta et applications à la cryptographie[END_REF]Proposition 6.3.5] and [31, Remarks 2.10.3, 2.10.7, 2.10.14] but never formalised. The isogeny formulae in the theta model can be seen as a natural generalisation of Vélu's formulae in higher dimension. In contrast, the Richelot correspondence is a very special relationship on hyperelliptic curves which does not seem to naturally generalise to arbitrary dimension. Nevertheless, the theta model has not received much attention as many have deemed the algorithms as impractical for efficient computation.

This paper is dedicated to disprove this claim of impracticality. We aim to demystify the hard algebraic geometry underpinning the theta model and make it accessible to cryptographers who want to employ isogenies between higher-dimensional abelian varieties within their protocols. The end result of our work is a set of concrete algorithms which describe the necessary pieces for computing isogenies between elliptic products; written to be particularly amenable to efficient and optimised implementations which are not all that different in appearance to the one-dimensional isogenies many are more familiar with.

Contributions

This paper has been written with the aim of being modular, using an algorithmic approach. All the formulae in the paper are mainly derived from the duplication formula. As a result, a reader uniquely interested in the computational results can assume the validity of the work in Section 2 and follow along the subsequent sections, which contain the explicit algorithms. From the duplication formula, we first re-obtain the addition formulae that have already been described in [START_REF] Gaudry | Fast genus 2 arithmetic based on theta functions[END_REF] and also give a precise operation counting in the base field.

The algorithm to compute chains of (2, 2)-isogenies between elliptic products is split into various subroutines. Each subroutine is carefully described in algorithmic boxes; this allows for a precise field operation counting. The main advantage of this approach is that both reducible and irreducible abelian surfaces can be described in the same way. However, some extra care will be devoted to the splitting and gluing case.

The gluing case is the most delicate one, where zero-coordinates must be carefully handled during both arithmetic and isogeny computations. To recover the theta-model representation of an elliptic product, we build upon the formulae described in [START_REF] Barbulescu | Models of Kummer lines and Galois representation[END_REF] which constructs theta structures on elliptic curves in the Montgomery model. We efficiently compute the theta-model representation of an elliptic product using only the dimension one representation of the theta structures and a few additional multiplications to recover the product structure. A summary of the costs of the algorithms described in this paper is shown in Table 1.

Table 1. Table of base field costs of doubling, codomain computation and evaluation for both generic and gluing isogenies in the theta model. We denote by M, S, I the costs of multiplication, squaring and inversion of an element in the base field and ignore the cost of additions. The generic codomain optimisation comes from reusing certain field elements also required for doubling. These are naturally available when computing long isogeny chains for a given kernel, and so the optimised cost is the expected cost for an isogeny between elliptic products.

Isogeny Type Doubling Codomain Evaluation

Generic Optimised

Generic 8S + 6M 8S + 29M + 1I 8S + 9M + 1I 4S + 3M Gluing 12S + 12M 8S + 13M + 1I 8S + 5M + 1I
Note that unlike the case of the Richelot chain, which requires both a (2, 2)-gluing and (2, 2)-splitting isogeny, computing the elliptic product at the end of a chain of isogenies in the theta model is a case of simply converting from one model to another, which can be done efficiently.

Finally, we offer both a constant time implementation of an isogeny between elliptic products in the programming language Rust, as well as an alternative implementation for the computer algebra system SageMath [START_REF]The Sage Developers: SageMath, the Sage Mathematics Software System[END_REF]. Both are available at the following GitHub repository: https://github.com/ThetaIsogenies/two-isogenies.

The Rust implementation has been written with cryptographic applications in mind and so has been built to run in constant time, with the appropriate finite field arithmetic and no secret-dependent conditional branching. It should also be easily portable to other projects in the future. The SageMath implementation has been designed to be a drop-in replacement for the work of [START_REF] Oudompheng | A Note on Reimplementing the Castryck-Decru Attack and Lessons Learned for SageMath[END_REF]. As a result, all the protocols whose implementation relies on this work or the proof-of-concept of [START_REF] Basso | FESTA: Fast Encryption from Supersingular Torsion Attacks[END_REF] can be upgraded to (2, 2)-isogenies in the theta model with minimal effort. As a use case, we show the benefit of these algorithms in FESTA in Section 5.3.

We give explicit timings of our implementations in Table 2. In SageMath, our implementation achieves a ten times speed up for the codomain computation and more than twenty times speed up for evaluation time compared to [START_REF] Oudompheng | A Note on Reimplementing the Castryck-Decru Attack and Lessons Learned for SageMath[END_REF]. For smaller characteristic, the Rust code runs approximately forty times faster than the same algorithms written in SageMath, and more than two times as fast for very large characteristic. Concretely, on an Intel Core i7-9750H CPU with a clock-speed of 2.6 GHz with turbo-boost disabled, we compute an isogeny chain of length n = 208 between elliptic products over F p 2 with a 254 bit characteristic in only 2.85 ms.

Roadmap In Section 2, we give a concise summary of the algebraic theory of theta functions. The most important parts of this section are the duplication formula and change-of-basis algorithm; the reader willing to accept these two main building blocks can skip this section entirely. In Section 3, we derive addition formulae from the duplication formula. The isogeny formulae are described in Section 4. We discuss our implementation results in Section 5 and draw some conclusions in Section 6.

Notation Throughout the paper, M, S, I will represent the cost of multiplication, squaring and inversion of an element in the base field, respectively. In Section 2, we will introduce the Hadamard transform H; in dimension two,

H(x, y, z, w) =     1 1 1 1 1 -1 1 -1 1 1 -1 -1 1 -1 -1 1     •     x y z w    
We also define ( θ00 (P ) : θ10 (P ) : θ01 (P ) : θ11 (P )) = H (θ 00 (P ), θ 10 (P ), θ 01 (P ), θ 11 (P )) to be the dual coordinates of P , and the ⋆ operator:

(x, y, z, w) ⋆ (x ′ , y ′ , z ′ , w ′ ) = (xx ′ , yy ′ , zz ′ , ww ′ ).
Another useful operator we will introduce in Section 3 is the squaring operator S: S(x, y, z, w) = (x 2 , y 2 , z 2 , w 2 ).

When computing the cost of inverting k elements, we will use batched inversions. Batched inversions allow us to invert k elements at a cost of 3(k -1) multiplications and only one inversion [20, §10.3.1]. We refer to our implementation for an explicit description of the algorithm.

Preliminaries

We assume the reader has some familiarity with (N, N )-isogenies between principally polarised abelian surfaces; we refer to [18, §2] for a gentle introduction to the subject. Before giving an explicit description of the algorithm used to compute chains of (2, 2)-isogenies between products of two elliptic curves, we provide a concise and selfcontained summary of the algebraic theory of theta functions. Using theta functions, it is possible to perform arithmetic on principally polarised abelian varieties. The reader willing to assume the validity of the duplication formula and our method to perform a change of basis can skip this section entirely and use Algorithm 2 as a black box. For all the other readers, in what follows, we utilise the language of Mumford's theory to provide a summary of the algebraic theory of theta functions [START_REF] Mumford | On the Equations Defining Abelian Varieties[END_REF][START_REF] Mumford | On the Equations Defining Abelian Varieties[END_REF][START_REF] Mumford | On the Equations Defining Abelian Varieties[END_REF]. We will first briefly recapitulate Mumford's results, then recall the Duplication formula, which will allow us to describe isogeny formulae between principally polarised abelian surfaces, and finally provide a formula for the change of basis in the case of a product of two elliptic curves, which is the case analysed in the paper.

Mumford's Theory

We refer to [START_REF] Milne | Abelian Varieties (v2.00[END_REF] for background material on abelian varieties and line bundles. Let A be an abelian variety defined over a field k whose characteristic is different from 2, and let L be a line bundle on A; we further assume that L is a totally symmetric line bundle [22, §2]. We define K(L) to be the set of points x such that the pullback of L under the translation-by-x map T x is isomorphic to L, i.e.

K(L)

:= {x ∈ A | T * x L ≃ L} .
Given a point x ∈ K(L), the isomorphism between L and T * x L is in general not unique. Also, given x, y ∈ K(L), together with two isomorphisms ϕ x : L ∼ -→ T * x L and ϕ y : L ∼ -→ T * y L, it is possible to define an isomorphism between L and T * x+y L as follows:

L ϕx --→ T * x L T * x ϕy ---→ T * x (T * y L) = T * x+y L.
This observation led Mumford to define the group

G(L) := {(x, ϕ x ) | x ∈ K(L), ϕ x : L ∼ -→ T * x L} ,
where (x, ϕ x ) • (y, ϕ y ) = (x + y, T * x ϕ y • ϕ x ). We will refer to G(L) as the Mumford theta group. The Mumford theta group fits in the exact sequence

0 → k * → G(L) → K(L) → 0.
Let f : A → B be an isogeny between abelian varieties, let L be a line bundle on A and M a line bundle on B. Suppose there exists an isomorphism α : f * M ∼ -→ L. Then, one can prove that Ker(f ) ⊂ K(L). On the other hand, given K ⊂ K(L), it is not generally true that the isogeny f ′ : A → B of kernel K generates an isomorphism between L and (f ′ ) * M ′ for some line bundle M ′ on B. This is exactly true when K admits a level subgroup K of G(L); a subgroup K of G(L) such that its image under the forgetful map (x, ϕ x ) → x coincides with K.

To understand when a subgroup K ⊂ K(L) admits a level subgroup, Mumford introduces a pairing on K(L): define e L : K(L) × K(L) → k * to be the pairing such that e L (x, y) = x • ỹ • x-1 • ỹ-1 , where x is any lift of x to G(L). 5 Then, there exists a level subgroup K over K if and only if e L ≡ 1 on K × K. Also, (K(L), e L ) is a symplectic space, which implies that there exists a decomposition into two maximal isotropic groups K(L) = K(L) 1 ⊕ K(L) 2 . Let δ := (d 1 , . . . , d g ) be the elementary divisors of K(L) 1 ≃ K(L) 2 . We say that L is of type δ. Define

K(δ) := g i=1 Z ⧸ d i Z , H(δ) := K(δ) ⊕ K(δ), G(δ) := k * × H(δ),
where K(δ) := Hom K(δ), k * . As K(L), H(δ) can be equipped with a pairing deriving from the duality: for all (x, χ),

(x ′ , χ ′ ) ∈ H(δ) e δ ((x, χ), (x, χ ′ )) = χ ′ (x)χ(x ′ ) -1 .
The set G(δ) is endowed with a group structure and is called Heisenberg group. This group structure is given explicitly by: for all (α, x, χ),

(α ′ , x ′ , χ ′ ) ∈ G(δ) (α, x, χ) • (α ′ , x ′ , χ ′ ) := (αα ′ χ ′ (x), x + x ′ , χ • χ ′ ).
The Heisenberg group G(δ) has a unique irreducible representation which acts trivially with k * ; such a representation is V (δ) := g : K(δ) → k . On the other hand, G(L) acts on the space of the global sections Γ (A, L), and Γ (A, L) is an irreducible representation for G(L). We denote the action of G(L) on Γ (A, L) simply by "•".

A theta structure Θ L of type δ is an isomorphism of central extensions from G(δ) to G(L) fitting in the diagram:

0 k * G(δ) H(δ) 0 0 k * G(L) K(L) 0 Θ L Θ L
Considering the horizontal arrows as the natural maps, the isomorphism Θ L induces another isomorphism Θ L between H(δ) and K(L), which is symplectic for the pairings e δ and e L . As a result, this isomorphism defines a symplectic decomposition on

K(L) = K(L) 1 ⊕ K(L) 2 , where K(L) 1 := Θ L (K(δ)) and K(L) 2 := Θ L ( K(δ)).
Using such a decomposition, we will see below that it is possible to identify a basis (θ x ) x∈K(L)1 for the space of global sections Γ (A, L). Since Θ L restricts to the identity with k * , Γ (A, L) is an irreducible representation of G(δ), which acts trivially with k * . Also, since V (δ) is the unique irreducible representation of H(δ), there exists (up to constants) a unique isomorphism β : Γ (A, L) ∼ -→ V (δ). The space V (δ) admits a canonical basis given by the Kronecker delta functions (δ i ) i∈K(δ) on K(δ). Via β, it is possible to transfer that basis to a basis (θ i ) i∈K(L)1 on Γ (A, L); the functions θ i are called theta coordinates, and, if d 1 = . . . = d g = n, these theta coordinates are called theta coordinates of level n. Using theta coordinates, it is possible to represent abelian varieties via an embedding into the projective space [25, Ch. II, Theorem 1.3]. In particular, the abelian variety A can be completely described in the projective space via the evaluation of theta coordinates at the identity using the Riemann relations; we call the projective point θ i (0 A ) i∈K(L)1 theta-null point. Given a point P ∈ A and T ∈ K(L), we can efficiently represent P + T in theta coordinates: if T corresponds to (s, χ) ∈ H(δ),

(θ i (P + T )) i = (χ(i)θ i+s (P )) i . (1)
This remark will allow us to compute (2, 2)-isogenies using projective coordinates.

Example 1. Let (a : b : c : d) be a theta-null point of level 2 obtained from the theta structure Θ L . Implicitly, we have a symplectic basis Another fundamental ingredient is the change of theta structure given by Heisenberg group automorphisms. A Heisenberg group automorphism is an automorphism of G(δ) acting as the identity on k * . In particular, such an automorphism induces a symplectic automorphism on H(δ) with respect to its natural pairing e δ . The most fundamental example is the Hadamard Transform, which is the automorphism that swaps K(δ) and K(δ).

(S 1 , S 2 , T 1 , T 2 ) of the 2-torsion. Let i 1 = ([1], [0]) ∈ K(2, 2), i 2 = ([1], [0]) ∈ K(2, 2), χ 1 ∈ K(2, 2) such that χ 1 (i 1 ) = -1 and χ 1 (i 2 ) = 1, and χ 2 ∈ K(2, 2) such that χ 2 (i 1 ) = 1 and χ 2 (i 2 ) = -1. Then, let us define S 1 = Θ L (i 1 ), S 2 = Θ L (i 2 ), T 1 = Θ L (χ 1 ) and T 2 = Θ L (χ 2 ).
Hadamard Transform Let (θ i ) i be some theta coordinates on A. The action of the Hadamard transform on (θ i ) i is described in [START_REF] Robert | Efficient algorithms for abelian varieties and their moduli spaces[END_REF]Eq. 2.4]. The resulting theta coordinates after this transform are called the dual theta coordinates; we will denote such coordinates by ( θi ) i . In what follows, we will use the Hadamard transform on level-2 theta coordinates. For the sake of clarity, we explicitly state the action of this symplectic automorphism in dimension one and two.

First, let us fix an ordering for theta coordinates. In dimension one, there are only two theta coordinates. Whenever we write (x : y) to represent a point P in theta coordinates, we actually mean (θ 0 (P ) : θ 1 (P )). Hence, specialising [31, Eq. 2.4], we obtain ( θ0 (P ) : θ1 (P )) = (x + y : x -y). In dimension two, we represent a point P in theta coordinates by a tuple (x : y : z : w), where we fix the ordering (θ 00 (P ) : θ 10 (P ) : θ 01 (P ) : θ 11 (P )). 6 Specialising [31, Eq. 2.4], we have θ00 (P ) = x + y + z + w,

θ10 (P ) = x -y + z -w, θ01 (P ) = x + y -z -w, θ11 (P ) = x -y -z + w.
Henceforth, we will use the operator H to refer to the action of the Hadamard transform on theta coordinates. Finally, we remark that H H (θ

A i ) i = (θ A i ) i .

Duplication Formula

Let (θ A i ) i be theta coordinates of level 2 on A. Implicitly, we have a symplectic de-

composition of A[2] = K(L) = K(L) 1 ⊕ K(L)
(x i ) i ⋆ (y i ) i = (x i y i ) i .
Then, a consequence of the isogeny theorem [22, Theorem 4, p. 302] (see also [START_REF] Robert | Fonctions thêta et applications à la cryptographie[END_REF]Theorem 3.6.4]) and duplication formula [22, Equation A, p. 332] shows that:

θ A i (P + Q) i ⋆ θ A i (P -Q) i = H θB i (f (P )) i ⋆ θB i (f (Q)) i , (2)
H θ A i ( f (R)) i ⋆ θ A i ( f (S)) i = θB i (R + S) i ⋆ θB i (R -S) i , (3)
where f denotes the dual isogeny of f and ( θB i ) i are the dual theta coordinates of (θ B i ) i .

Change of Basis

We now focus on products of two elliptic curves and explain how to endow these products with a theta structure of type (2, 2). First, let us recall that every theta structure of type (2, . . . , 2) comes from a symplectic basis of the

i , f i ⟩, for i = 1, 2. Then, (a 1 a 2 : b 1 a 2 : a 1 b 2 : b 1 b 2 ) is a theta-null point for E 1 × E 2 induced by the symplectic 4-torsion basis ⟨(e 1 , 0), (0, e 2 )⟩ ⊕ ⟨(f 1 , 0), (0, f 2 )⟩.
However, in the next sections, we might need to work with theta structures associated to a different symplectic 4-torsion basis. The change-of-basis algorithm we describe here is a generalisation of the case described in [START_REF] Barbulescu | Models of Kummer lines and Galois representation[END_REF], which uses a different approach based on the ramification of the Kummer line. We briefly describe the idea below and adapt it to the case of products of two elliptic curves.

Let E be an elliptic curve, and ⟨T ′ 1 , T ′ 2 ⟩ a basis of the 4-torsion. We explain how to compute the theta-null point associated to this basis. Given a point T ∈ E [START_REF] Basso | FESTA: Fast Encryption from Supersingular Torsion Attacks[END_REF], there are two symmetric elements ±g above T in the theta group G (L(2(0 E ))). We can fix a symmetric element via a point

T ′ of 4-torsion above T . Let T 1 = 2T ′ 1 , T 2 = 2T ′ 2
, and let g 1 , g 2 be these elements associated to T ′ 1 and T ′ 2 , respectively. By definition of the compatible theta basis, the theta coordinate θ 0 must be invariant under the action of g 2 , and θ 1 = g 1 • θ 0 . The coordinate θ 0 can be computed as the trace of a global section

s ∈ Γ (E, L(2(0 E ))), provided it is not equal to zero, i.e. θ 0 = id • s + g 2 • s ̸ = 0.
Working on a Montgomery curve in Weierstrass coordinates, we have a canonical point of 4-torsion T ′ = (1 : 1) above T = (0 : 1) that induces the canonical element g of the theta group acting by g • (X, Z) = (Z, X). Indeed, translation by T is given by (X : Z) → (Z : X), and the two symmetric elements above this translation act by (X, Z) → (±Z, ±X) since they have order 2. The element g T ′ in G (L(2(0 E ))) fixed by T ′ corresponds to ±g. But since g leaves invariant any affine lift of T ′ , we have

g T ′ = g.
For a general elliptic curve, if T ′ = (x, y, z) is a point of 4-torsion and 2T ′ = T = (u, v, w), we can map T ′ to the Montgomery point (1 : 1) via the linear transformation (in the Kummer line):

M : (X : Z) → (X ′ : Z ′ ) = (zwX -zuZ : (xw -zu)Z). It follows that the action of g T ′ is given by M T U M T -1 = uz/(wx -uz) zw/(wx -uz) (-wx 2 + 2uzx)/(-zwx + uz 2 ) -uz/(wx -uz) , with M = wz -zu 0 xw -uz , U = 0 1 1 0
. This computation is the output of Algorithm 1.

Example 3. Using previous notation, on a Montgomery curve we have T ′ 2 = (-1 : 1), which acts by g 2 • (X, Z) = (-Z, -X). Taking the trace of X under this action we get θ

0 = id • X + g 2 • X = X -Z. Let T ′ 1 = (a + b : a -b) be another point of 4-torsion; its double is then (a 2 + b 2 : a 2 -b 2 ). Let x = a + b, z = a -b, u = a 2 + b 2 , w = a 2 -b 2 . We compute θ 1 = g 1 • θ 0 = g 1 • (X -Z) = (uz -wxx/z + 2ux)/(wx -uz)X + z(w + u)/(wx -uz)Z = b/aX + b/aZ.
We recover the same conversion formula between Montgomery and theta coordinates as obtained in [START_REF] Barbulescu | Models of Kummer lines and Galois representation[END_REF].

We can use the same strategy to compute the theta-null point associated to a symplectic basis of the 4-torsion on a product of elliptic curves. If

T ′ = (T ′ 1 , T ′ 2 ) ∈ E 1 × E 2 is a point of 4-torsion, the associated element g T ′ is given by g T ′ = g T ′ 1 ⊗ g T ′ 2 . We can take θ 0 as the trace of X 1 ⊗X 2 , i.e. θ 0 = i g i •X 1 ⊗X 2 = i g i,1 •X 1 ⊗g i,2 •X 2 ,
where the g i = g i,1 ⊗ g i,2 's are the elements above K 2 fixed by the 4-torsion. The other theta coordinates are computed via the action of the elements above K 1 on θ 0 .

Algorithm 1 Action by Translation

Input: A point P ′ in the 4-torsion of the Kummer line of an elliptic curve Output: The 2 × 2 submatrix M with coefficients mij describing the action of g P ′ lifting the action by translation of P = 2P ′ . 1:

P ← [2]P ′ (▷) Cost: 3M + 2S 2: Let P ′ = (X : Z) and (U : W ) = P 3: W X, W Z, U X, U Z ← W • X, W • Z, U • X, U • Z 4: δ ← W X -U Z 5: Compute δ -1 , Z -1 via batched inversions 6: m00 ← -U Z • δ -1 7: m01 ← -W Z • δ -1 8: m10 ← U X • δ -1 -X • Z -1 9: m11 ← -m00 10: return M (▷) Total cost: 14M + 2S + 1I
In practice, in the algorithm to compute the

(2 n , 2 n )-isogeny f : E 1 ×E 2 → E ′ 1 ×E ′ 2 , we only have access to ker(f )[4] = ⟨T ′ 1 , T ′ 2 ⟩ and not to a complete symplectic torsion basis of E 1 × E 2 [4]; let T ′ 1 = (P 1 , P 2 ) and T ′ 2 = (Q 1 , Q 2 ).
To bypass this problem, we define

S ′ 1 = (0, Q 2 ) and S ′ 2 = (P 1 , 0). Then, K 1 = ⟨S 1 , S 2 ⟩ and K 2 = ⟨T 1 , T 2 ⟩ , where S i = [2]S ′ i and T i = [2]T ′ i . We use the symplectic 4-torsion basis ⟨S ′ 1 , S ′ 2 , T ′ 1 , T ′ 2 ⟩
when endowing E 1 ×E 2 with a theta structure. We summarise this procedure in Algorithm 2. The output of this algorithm is a matrix N that allows for a change of coordinates as follows.

If R = (R 1 , R 2 ) ∈ E 1 × E 2 is a point in Weierstrass coordinates for the Montgomery elliptic curves, where R i = (X i : Z i ), the image of R in theta coordinates is given by N •     X 1 • X 2 X 1 • Z 2 Z 1 • X 2 Z 1 • Z 2     .
Remark 4. In Algorithm 2, it is possible to optimise the computation of the four inversions required in the four calls to Algorithm 1 in lines 1, 2, 3 and 4 by using a unique batched inversion. We decided not to show this optimisation in Algorithm 2 for the sake of a cleaner exposition.

Addition Formulae

In this section, we derive addition formulae using the equations in Subsection 2.2. These formulae have already been described in dimension one [START_REF] Barbulescu | Models of Kummer lines and Galois representation[END_REF] and dimension two [START_REF] Gaudry | Fast genus 2 arithmetic based on theta functions[END_REF]. However, we prefer to restate them in dimension two to highlight the connection with (2, 2)-isogenies and provide explicit operation counting. In what follows, we use the same notation as in Subsection 2.2.

Let P, Q ∈ A and suppose we have

(θ A i (P -Q)) i . To compute (θ(P + Q) i ) i ,
Cost: 56M + 8S + 4I 2: G2 ← action_by_translation(P ′ 2 ) 3: H1 ← action_by_translation(Q ′ 1 ) 4: H2 ← action_by_translation(Q ′ 2 ) 5: t 00|1 ← g 00|1 • h 00|1 + g 01|1 • h 10|1 (▷) Compute the first column of G1 × H1 6: t 10|1 ← g 10|1 • h 00|1 + g 11|1 • h 10|1 7: t 00|2 ← g 00|2 • h 00|2 + g 01|2 • h 10|2 (▷) Compute the first column of G2 × H2 8: t 10|2 ← g 10|2 • h 00|2 + g 11|2 • h 10|2 9: n00 ← g 00|1 • g 00|2 + h 00|1 • h 00|2 + t 00|1 • t 00|2 + 1 (▷) Compute the trace for the first row 10: n01 ← g 00|1 • g 10|2 + h 00|1 • h 10|2 + t 00|1 • t 10|2 11: n02 ← g 10|1 • g 00|2 + h 10|1 • h 00|2 + t 10|1 • t 00|2 12: n03 ← g 10|1 • g 10|2 + h 10|1 • h 10|2 + t 10|1 • t 10|2 13: n10 ← h 00|2 • n00 + h 01|2 • n01 (▷) Compute the action of (0, Q ′ 2 ) for the second row 14: n11 ← h 10|2 • n00 + h 11|2 • n01 15: n12 ← h 00|2 • n02 + h 01|2 • n03 16: n13 ← h 10|2 • n02 + h 11|2 • n03 17: n20 ← g 00|1 • n00 + g 01|1 • n02
(▷) Compute the action of (P ′ 1 , 0) for the third row 18:

n21 ← g 00|1 • n01 + g 01|1 • n03 19: n22 ← g 10|1 • n00 + g 11|1 • n02 20: n23 ← g 10|1 • n01 + g 11|1 • n03 21: n30 ← g 00|1 • n10 + g 01|1 • n12 (▷) Compute the action of (P ′ 1 , Q ′ 2 ) for the final row 22: n31 ← g 00|1 • n11 + g 01|1 • n13 23: n32 ← g 10|1 • n10 + g 11|1 • n12 24: n33 ← g 10|1 • n11 + g 11|1 • n13 25: return N (▷) Total cost: 100M + 8S + 4I
be computed as

θB i (f (P )) i ⋆ θB i (0) i = H θ A i (P ) i ⋆ θ A i (P ) i ,
and similarly for θB

i (f (Q)) i . The quantity θB i (0) i is actually not needed, as we only need θB i (0) i ⋆ θB i (0) i if we use θB i (f (P )) i ⋆ θB i (0) i ⋆ θB i (f (Q)) i ⋆ θB i (0) i , to compute θB i (f (P )) i ⋆ θB i (f (Q)) i .
For the sake of compactness, we introduce the operator S that, on input θ A i (P ) i , returns θ A i (P ) i ⋆ θ A i (P ) i . We formalise this procedure in Algorithm 3.

Let (a : b : c : d) be a theta-null point for A and define (α : β : γ : δ) to be the dual coordinates ( θB i (0)) i of the theta-null point (θ B (0) i ) i . For simplicity, let us assume that α • β • γ • δ ̸ = 0; the (rare) case when one of the dual coordinates is zero is briefly treated in Remark 5. Equation 2 proves that (α 2 :

β 2 : γ 2 : δ 2 ) = H(a 2 : b 2 : c 2 : d 2 )
. However, since we are working projectively, we need to use the quantities (α 2 /β 2 , α 2 /γ 2 , α 2 /δ 2 ) and (a/b, a/c, a/d), which can we computed via batched inversions.

Algorithm 3 Differential addition

Input: The theta coordinates (xP : yP : zP : wP ) of P , the theta coordinates (xQ : yQ :

zQ : wQ) of Q, the theta coordinates (xP -Q : yP -Q : zP -Q : wP -Q) of P -Q and ( λ1, λ2, λ3) = (α 2 /β 2 , α 2 /γ 2 , α 2 /δ 2
X f (P )f (Q) ← XP • XQ 4: Y f (P )f (Q) ← λ1 • YP • YQ 5: Z f (P )f (Q) ← λ2 • ZP • ZQ 6: W f (P )f (Q) ← λ3 • WP • WQ 7: XP Q, YP Q, ZP Q, WP Q ← H(X f (P )f (Q) , Y f (P )f (Q) , Z f (P )f (Q) , W f (P )f (Q) ) 8: xyP -Q, ztP -Q ← xP -Q • yP -Q, zP -Q • tP -Q 9: xP +Q ← XP Q • ztP -Q • yP -Q 10: yP +Q ← YP Q • ztP -Q • xP -Q 11: zP +Q ← ZP Q • xyP -Q • wP -Q 12: wP +Q ← WP Q • xyP -Q • zP -Q 13: return xP +Q, yP +Q, zP +Q, wP +Q (▷) Total cost: 8S + 17M
To obtain an algorithm to double a point P ∈ A, we proceed as before with the only difference that we only need (a/b, a/c, a/d). We provide a detailed description of the doubling in Algorithm 4.

Algorithm 4 Doubling

Input: The theta coordinates (xP : yP : zP : wP

) of P , ( λ1, λ2, λ3) = (α 2 /β 2 , α 2 /γ 2 , α 2 /δ 2 )
and (λ1, λ2, λ3) = (a/b, a/c, a/d). Output: The theta coordinates (x2P : y2P : z2P : w2P ) of 2P .

1: XP , YP , ZP , WP ← H • S(xP , yP , zP , wP ) 2:

X 2 f (P ) ← XP • XP 3: Y 2 f (P ) ← λ1 • YP • YP 4: Z 2 f (P ) ← λ2 • ZP • ZP 5: W 2 f (P ) ← λ3 • WP • WP 6: X ′ P , Y ′ P , Z ′ P , W ′ P ← H(X 2 f (P ) , Y 2 f (P ) , Z 2 f (P ) , W 2 f (P ) ) 7: return X ′ P , λ1 • Y ′ P , λ2 • Z ′ P , λ3 • W ′ P (▷) Total cost: 8S + 6M
Remark 5. In practice, we will always be in the case that α • β • γ • δ ̸ = 0. We may incur in such an exception when we are working on a product of elliptic curves but with a non product theta structure (see [35, § 16.4]): this is due to the fact we may perform a change of basis as in Subsection 2.3 to put the kernel in the "right position" -see Section 4.1.

In this case we do not have to worry since we could perform arithmetic on the elliptic curves and then convert to the theta model afterwards. However, if one wants to deal with this case in the theta model, one may first act with a symplectic automorphism ψ sending the dual theta-null point to a all-non-zero one, use the addition formulae and eventually switch back to the former theta structure acting by ψ -1 .

To compute (2 n , 2 n )-isogenies, we will only use doublings. Algorithm 4 works only

if a • b • c • d ̸ = 0. The case a • b • c • d = 0 can
happen only if the codomain B is a product of elliptic curves with non product theta structure by [35, § 16.4]. In this case, we can use the same solution as above, by using ψ = H as our symplectic transformation.

The Isogeny Formula

In this section, we explain how to derive an isogeny formula for (2, 2)-isogenies from Subsection 2.2. Ultimately, we focus on isogenies between products of two elliptic curves. However, we first show how to compute isogenies when we have an abelian surface already endowed with a theta structure compatible with the kernel of the (2, 2)isogeny we want to compute.

Let A be an abelian surface endowed with a theta structure of type (2, 2), and let ⟨S 1 , S 2 , T 1 , T 2 ⟩ be the canonical symplectic basis associated with the symplectic decomposition A[2] = K 1 ⊕ K 2 ; to be more specific,

K 1 = ⟨S 1 , S 2 ⟩ and K 2 = ⟨T 1 , T 2 ⟩.
Let us recall that a theta-null point is fixed by a symplectic basis of the 4-torsion [22, Remark 4, p. 319], and let

⟨S ′ 1 , S ′ 2 , T ′ 1 , T ′ 2 ⟩ be such a basis; in particular 2S ′ i = S i and 2T ′ i = T i . Our goal is to compute a (2, 2)-isogeny f : A → B.
If the kernel ker f does not coincide with K 2 , we can always apply a change of basis as described in Subsection 2.3. Hence, we can assume that ker f = K 2 .

Before outlining the explicit procedure, let us assume that we have

T ′′ 1 , T ′′ 2 such that ⟨T ′′ 1 , T ′′ 2 ⟩[4] = ⟨T ′ 1 , T ′ 2 ⟩, 2T ′′ i = T ′
i and their Weil pairing e 8 (T ′′ 1 , T ′′ 2 ) = 1. 7 These conditions are not restrictive since they are naturally verified for chains of (2, 2)-isogenies, which are our end goal. In particular, ⟨f (T ′′ 1 ), f (T ′′ 2 )⟩ are two of the 4-torsion points inducing the theta-null point on B laying above the two 2-torsion points in K 2 . Hence, the points ⟨f (T ′′ 1 ), f (T ′′ 2 )⟩ lay above the canonical points in K 1 for the dual coordinates. Let us remark that f (

T ′′ i ) + 2f (T ′′ i ) = -f (T ′′ i ).
So, as highlighted in Equation 1 and since we are on the Kummer, we have

θB 00 (f (T ′′ 1 )) : θB 10 (f (T ′′ 1 )) : θB 01 (f (T ′′ 1 )) : θB 11 (f (T ′′ 1 )) = θB 10 (f (T ′′ 1 )) : θB 00 (f (T ′′ 1 )) : θB 11 (f (T ′′ 1 )) : θB 01 (f (T ′′ 1 )) , and
θB 00 (f (T ′′ 2 )) : θB 10 (f (T ′′ 2 )) : θB 01 (f (T ′′ 2 )) : θB 11 (f (T ′′ 2 )) = θB 01 (f (T ′′ 2 )) : θB 11 (f (T ′′ 2 )) : θB 00 (f (T ′′ 2 )) : θB 10 (f (T ′′ 2 )) .
7 Such a subgroup ⟨T ′′ 1 , T ′′ 2 ⟩ is said to be isotropic.

Define (α : β : γ : δ) to be the dual theta-null point of B, i.e. θB 00 (0) : θB 10 (0) : θB 01 (0) : θB 11 (0) = (α : β : γ : δ).

Then, combining Equation 2 with the above observations, we have

H • S(θ A 00 (T ′′ 1 ), θ A 10 (T ′′ 1 ), θ A 01 (T ′′ 1 ), θ A 11 (T ′′ 1 )) = (xα, xβ, yγ, yδ), H • S(θ A 00 (T ′′ 2 ), θ A 10 (T ′′ 2 ), θ A 01 (T ′′ 2 ), θ A 11 (T ′′ 2 )) = (zα, wβ, zγ, wδ),
for some unknown x, y, z, t. Hence, we can recover the dual theta-null point (α : β : γ : δ) for B, and in turn its theta-null point H(α : β : γ : δ).

Remark 6 (Technical Remark). It is possible to prove that all x, y, z, t must be different from zero. If it had not been the case, we would have ended up with a theta-null point with at least two zero coordinates. This is a contradiction since it implies we have more than a zero even theta-null coordinate of level (2, 2) -see Subsection "Gluing Isogeny" for the definition of level-(2, 2) theta coordinates.

In general, the dual theta-null point (α : β : γ : δ) has all coordinates different from zero. The only exceptions can be found for certain cases of the gluing isogeny -we will discuss how to handle this case in Subsection 4.1. Once we have computed (α : β : γ : δ), we can evaluate the isogeny f at any point P using (again) Equation 2. First, we compute

(x ′ , y ′ , z ′ , w ′ ) = H • S (θ A i (P )) i . Then (θ B i (f (P ))) i = H(α -1 x ′ , β -1 y ′ , γ -1 z ′ , δ -1 w ′ ).
We give a detailed description in Algorithms 5 and 6.

Algorithm 5 Codomain

Input: Theta coordinates (x T ′′ 1 :

y T ′′ 1 : z T ′′ 1 : w T ′′ 1 ) of T ′′ 1 , theta coordinates (x T ′′ 2 : y T ′′ 2 : z T ′′ 2 : w T ′′ 2 ) of T ′′ 2 ,
where T ′′ i is a 8-torsion point laying above the K2 part of the symplectic 4-torsion basis inducing the theta-null point. Output: Dual theta-null point (1 : β : γ : δ), the inverse of the dual theta-null point (1 :

β -1 : γ -1 : δ -1 ) and the theta-null point (a

′ : b ′ : c ′ : d ′ ) on B. (▷) Case β • γ • δ ̸ = 0 1: (xα, xβ, yγ, yδ) ← H • S(x T ′′ 1 , y T ′′ 1 , z T ′′ 1 , w T ′′ 1 ) 2: (zα, wβ, zγ, wδ) ← H • S(x T ′′ 2 , y T ′′ 2 , z T ′′ 2 , w T ′′ 2 ) 3: Compute (xα) -1 , (xβ) -1 , (yγ) -1 , (yδ) -1 , (zα) -1 , (wβ) -1 , (zγ) -1 , (wδ) -1 via batched in- versions. 4: β ← xβ • (xα) -1 5: γ ← zγ • (zα) -1 6: δ ← wδ • (wβ) -1 • β 7: β -1 ← xα • (xβ) -1 8: γ -1 ← zα • (zγ) -1 9: δ -1 ← wβ • (wδ) -1 • β -1 10: (a ′ , b ′ , c ′ , d ′ ) ← H(1, β, γ, δ) 11: return (1, β, γ, δ), (1, β -1 , γ -1 , δ -1 ), (a ′ , b ′ , c ′ , d ′ ) (▷) Total cost: 8S + 29M + 1I

Algorithm 6 Evaluation

Input: Theta coordinates (xP : yP : zP : wP ) of P and the dual theta-null point (1 : β -1 : γ -1 : δ -1 ) on B. Output: Theta coordinates (x f (P ) : y f (P ) : z f (P ) :

w f (P ) ) of f (P ). (▷) Case β • γ • δ ̸ = 0 1: (XP , YP , ZP , WP ) ← H • S(xP , yP , zP , wP ) 2: (X ′ f (P ) , Y ′ f (P ) , Z ′ f (P ) , W ′ f (P ) ) ← (XP , β -1 • YP , γ -1 • ZP , δ -1 • WP ) 3: (x f (P ) , y f (P ) , z f (P ) , w f (P ) ) ← H(X ′ f (P ) , Y ′ f (P ) , Z ′ f (P ) , W ′ f (P ) ) 4: return (x f (P ) , y f (P ) , z f (P ) , w f (P ) ) (▷) Total cost: 4S + 3M
Remark 7. As we explained in Section 3, the inverse squared dual theta-null point (α 2 /β 2 , α 2 /γ 2 , α 2 /δ 2 ) are also needed for the addition and doubling formulae. If such a quantity has already been computed, we can use it to lower down the cost in Algorithm 5. In line 3, we need only to invert three elements, namely xα, zα, wβ. Then, to obtain (β -1 , γ -1 , δ -1 ), we can simply multiply component-wise (α 2 /β 2 , α 2 /γ 2 , α 2 /δ 2 ) by (β, γ, δ). The total cost in this optimised case is 8S + 9M + 1I.

Computing

(2 n , 2 n )-isogenies between Elliptic Products
Now, we specialise to the case of a (2 n , 2 n )-isogeny f :

E 1 × E 2 → E ′ 1 × E ′ 2
between elliptic products, which will be computed as a chain of (2, 2)-isogenies. Let K be the kernel of this isogeny and suppose that we have two points of order 2 n+2 on E 1 × E 2 above K forming an isotropic group. To apply the formulae we described above, we need to be sure that K [START_REF] Chen | SCALLOP-HD: group action from 2-dimensional isogenies[END_REF] is in "the right position". Given K [START_REF] Chen | SCALLOP-HD: group action from 2-dimensional isogenies[END_REF], we apply Algorithm 2 to obtain a theta-null point induced by the symplectic 4-torsion decomposition

⟨S ′ 1 , S ′ 2 ⟩ ⊕ ⟨T ′ 1 , T ′ 2 ⟩, where K[4] = ⟨T ′ 1 , T ′ 2 ⟩ If n > 2, we have that K[8] = ⟨T ′′ 1 , T ′′ 2 ⟩
is the isotropic 8-torsion above K [START_REF] Chen | SCALLOP-HD: group action from 2-dimensional isogenies[END_REF]. This means we could apply Algorithms 5 and 6 to compute the first step of the isogeny f , i.e. the isogeny f 1 : E 1 × E 2 → A 1 with kernel K [START_REF] Basso | FESTA: Fast Encryption from Supersingular Torsion Attacks[END_REF]. However, we should be careful as on the product structure, one of the coordinates of the dual theta-null point on A 1 may be equal to zero. We explain why this happens and how to bypass this obstacle in the Subsection "Gluing Isogeny" below.

After the first step, we also end up with a complete description of the theta structure on A

1 ; let A 1 [2] = K 1 ⊕ K 2 . The points f 1 (T ′′ 1 )
and f 1 (T ′′ 2 ) are two of the 4-torsion elements describing the theta-null point on A 1 . If n > 3, we can use f 1 (K) [START_REF] Decru | Towards a Quantum-Resistant Weak Verifiable Delay Function[END_REF] to describe the 8-torsion above ⟨f 1 (T ′′ 1 ), f 1 (T ′′ 2 )⟩ and iterate the process. Once we reach the second last step f n-1 : A n-2 → A n-1 , we cannot inherit the 8-torsion above the K 2 part of the A n-2 anymore. However, thanks to the assumption that we have to two points of order 2 n+2 on E 1 × E 2 above K forming an isotropic group, we can use the same exact strategy using the images of such points. We explain how to relax this condition in Subsection 4.2.

In the last step

f n : A n-1 → E ′ 1 ×E ′ 2 ,
we map onto an elliptic product. Even though we can reuse the same computational strategies when we stay in the theta model, for most of the cryptographic applications we have to explicitly recover the equations of the curves E ′

1 and E ′ 2 , and we also have to have map points onto these curves. We describe how to do so in the Subsection "Splitting Isogeny".

Gluing Isogeny

In this subsection, we focus on the first step of the isogeny chain, an isogeny originating from an elliptic product; let f : E 1 × E 2 → A be such an isogeny. Theta structures on elliptic products E 1 ×E 2 verify some additional properties with respect to level-(2, 2) theta coordinates. We refer to Dupont's PhD thesis [START_REF] Dupont | Moyenne arithmético-géométrique, suites de Borchardt et applications[END_REF] for background material. For the case at hand, we briefly recall some fundamental facts.

Theta coordinates of level (2, 2) are indexed by a pair of elements in K(2, 2). A level-(2, 2) theta coordinate U i,j is said to be even if i • j T = 0 (mod 2); otherwise it is said to be odd. Moreover, at most one of the even indices (i, j) verifies U i,j (0) = 0, and there is exactly one zero even index if and only if the theta structure is associated with a product of two elliptic curves.

Given level-2 theta coordinates (θ i (P )) i , we can compute the square of its level-(2, 2) theta coordinates as

U 2 i,j (P ) = t (-1) i•t T θ t (P )θ t+j (P ).
In [9, Proposition 6.5], Dupont shows that a theta-null point (θ i (0)) i comes from the product theta structure of two elliptic curves if and only if U 11,11 (0) = 0. This means that if we are working on a product structure, all the coordinates of the dual theta-null point (α : β : γ : δ) are non-zero since (U 00,00 (0) : U 10,00 (0) : U 01,00 (0) : U 11,00 (0)) = (α : β : γ : δ). However, when we perform a change of basis, we might move the zero even index of the level-(2, 2) theta-null point around, and potentially we might have one of the dual theta-null coordinates equal to zero. In fact, unless A is a product of elliptic curves (which would be the case if the kernel is a product kernel), then we know that one of α, β, γ, δ is zero. If it were not the case, we could compute f (P ) from P . But since we work with theta coordinates of level 2, on the product E 1 ×E 2 we are really working with the product of Kummer lines E 1 / ± 1 × E 2 / ± 1. The automorphism group by which we quotient is thus Z/2Z × Z/2Z compared to Z/2Z when working on the Kummer surface A/ ± 1 of an abelian surface with a non product principal polarisation. Thus, when going from P to f (P ), there is an ambiguity coming from an action of Z/2Z, which can only be resolved by either taking a square root, or as we will explain next, by using extra information coming from the arithmetic of E 1 × E 2 (rather than E 1 / ± 1 × E 2 / ± 1).

Let us handle the case where one of the coordinates of the dual theta-null point (α : β : γ : δ) is zero; let us first analyse the case α = 0. In Algorithm 5, we normalised everything with respect to δ. This actually simplifies the codomain computation. We explain how to do so in Algorithm 7.

As explained above, mapping points under this isogeny requires extra care. If we simply use Algorithm 6, we cannot retrieve the first coordinate of a point. To be precise, if we want to evaluate f at the point (θ E1×E2 i (P )) i , we have

H • S((θ E1×E2 i (P )) i ) = (0, β θA 10 (f (P )), γ θA 01 (f (P )), δ θA 11 (f (P ))). (4)
Multiplying by β -1 , γ -1 and δ -1 the components β θA 10 (f (P )), γ θA 01 (f (P )), δ θA 11 (f (P )), we retrieve all the dual components but θA 00 (f (P )). The component θA 00 (f (P )) can be computed using the additional information coming from the theta structure. Let T ′ 1 be the point above T 1 ∈ K 2 as in the previous Algorithm 7 Special Codomain, α = 0

Input: In the same notation as before, theta coordinates

(x T ′′ 1 : y T ′′ 1 : z T ′′ 1 : w T ′′ 1 ) of T ′′ 1 , theta coordinates (x T ′′ 2 : y T ′′ 2 : z T ′′ 2 : w T ′′ 2 ) of T ′′ 2 ,
where T ′′ i is a 8-torsion point laying above the K2 part of the symplectic 4-torsion basis inducing the theta-null point. Output: Dual theta-null point (0 : β : γ : 1), the "inverse" of the dual theta-null point (0 : β -1 : γ -1 : 1) and the theta-null point (a

′ : b ′ : c ′ : d ′ ) on A. (▷) Case α = 0 1: (0, xβ, yγ, yδ) ← H • S(x T ′′ 1 , y T ′′ 1 , z T ′′ 1 , w T ′′ 1 ) 2: (0, wβ, zγ, wδ) ← H • S(x T ′′ 2 , y T ′′ 2 , z T ′′ 2 , w T ′′ 2 ) 3: Compute (yγ) -1 , (wβ) -1 , (yδ) -1 , (wδ) -1 via batched inversions. 4: β ← wβ • (wδ) -1 5: γ ← yγ • (yδ) -1 6: β -1 ← wδ • (wβ) -1 7: γ -1 ← yδ • (yγ) -1 8: (a ′ , b ′ , c ′ , d ′ ) ← H(0, β, γ, 1) 9: return (0, β, γ, 1), (0, β -1 , γ -1 , 1), (a ′ , b ′ , c ′ , d ′ ) (▷) Total cost: 8S + 13M + 1I
subsection. Then,

H • S((θ E1×E2 i (P + T ′ 1 )) i ) = (0, β θA 00 (f (P )), γ θA 11 (f (P )), δ θA 01 (f (P ))). (5)
However, multiplying the component β θA 00 (f (P )) by β -1 is not enough since we are working up to projective factors.

Once we recover θA 10 (f (P )), θA 01 (f (P )), θA 11 (f (P )) from Equation 4, we can compute λ θA 01 (f (P )) from Equation 5for some projective factor λ: we simply multiply the last component of H • S((θ E1×E2 i We summarise everything in Algorithm 8. Note that for both Algorithm 7 and 8, the case for β, γ or δ = 0 follows almost identically, see the implementation for a concrete example of how all cases can be considered concisely. We note that Algorithm 8 requires the knowledge not only of the theta coordinates of P , but also of P + T ′ 1 . From the knowledge of the (θ i (P )) and (θ i (T ′ 1 )), we may only recover (θ i (P ± T ′ 1 )), hence extract (θ i (P + T ′ 1 )) via a square root, consistent with the fact that in a gluing isogeny we have an ambiguity for images coming from an action by Z/2Z. Luckily we can compute this addition on each elliptic curve separately, using Weierstrass coordinates, before switching to the level-2 theta coordinates on the surface

(P + T ′ 1 )) i ) by δ -1 . If θA 01 (f (P )) ̸ = 0,
E 1 × E 2 .
Splitting Isogeny In this last step of the isogeny chain, we need to compute an isogeny f : A → E ′ 1 × E ′ 2 mapping onto an elliptic product. We can compute the thetanull point of E ′ 1 × E ′ 2 and mapping points under f using Algorithms 5 and 6. However, Algorithm 8 Special Evaluation, α = 0

Input: Theta coordinates (xP : yP : zP : wP ) of P , theta coordinates (x P +T ′ 1 : y P +T ′ 1 : z P +T ′ 1 : w P +T ′ 1 ) of P + T ′ 1 and the inverse of the dual theta-null point (0 : β -1 : γ -1 : 1) on A. Output: Theta coordinates (x f (P ) : y f (P ) : z f (P ) : w f (P ) ) of f (P ).

(▷) Case α = 0 1: (0, YP , ZP , WP ) ← H • S(xP , yP , zP , wP )

2: (0, YP +T 1 , ZP +T 1 , WP +T 1 ) ← H • S(xP +T 1 : yP +T 1 : zP +T 1 : wP +T 1 ) 3: (Y ′ f (P ) , Z ′ f (P ) , W ′ f (P ) ) ← (β -1 • YP , γ -1 • ZP , WP ) 4: if Z ′ f (P ) ̸ = 0 then 5: λ -1 ← Z ′ f (P ) /WP +T 1 6: else 7: Z ′ f (P +T 1 ) ← γ -1 • ZP +T 1 8: λ -1 ← W ′ f (P ) /Z ′ P +T 1 9: X ′ f (P ) ← λ -1 • β -1 • YP +T 1 10: (x f (P ) , y f (P ) , z f (P ) , w f (P ) ) ← H(X ′ f (P ) , Y ′ f (P ) , Z ′ f (P ) , W ′ f (P ) ) 11: return (x f (P ) , y f (P ) , z f (P ) , w f (P ) ) (▷) Total cost: 8S + 5M + 1I
we still need to retrieve the explicit equations for the curves E ′ 1 and E ′ 2 . This can be done using level-(2, 2) theta coordinates U i,j .

Since the theta structure on the image surface underlies an elliptic product, we know that one of the even indices -say (i, j) -of the level-(2, 2) theta-null point is equal to zero. Also, we know that if we compute a symplectic automorphism ψ mapping (i, j) onto [START_REF] Gaudry | Fast genus 2 arithmetic based on theta functions[END_REF][START_REF] Gaudry | Fast genus 2 arithmetic based on theta functions[END_REF], the action of ψ on the theta-null point obtained via Algorithm 5 gives back a theta-null point associated with the product theta structure.

From the above, it can be seen that there are ten distinct even indices. For each of these indices, we computed a symplectic automorphism sending this index to [START_REF] Gaudry | Fast genus 2 arithmetic based on theta functions[END_REF][START_REF] Gaudry | Fast genus 2 arithmetic based on theta functions[END_REF]. For efficiency reasons, we hardcoded the action of each of the symplectic automorphisms onto theta coordinates of level 2 in the reference implementation. These symplectic automorphisms and their actions have been derived from [31, p. 28] using the following sequential steps.

Let (i, j) be the even index such that U i,j (0) = 0, and, for ease of notation, let (a 00 : a 10 : a 01 : a 11 ) be the underlying theta-null point.

1. If i = j = 00, we act by the symplectic automorphism with matrix form

    1 0 2 0 0 1 0 0 0 0 1 0 0 0 0 1     .
We then obtain the theta-null point (a 00 : √ -1 • a 10 : a 01 : √ -1 • a 11 ), which means that U 10,00 (0) = 0. 2. If j = 00 and i ̸ = 00, we act by H, which swaps the roles of i and j. We can now assume that j ̸ = 00.

3. Let A be any invertible matrix such that A • j T = 11 T . Then, the action of the symplectic automorphism with matrix

A 0 0 A T -1
maps the theta-null point (a 00 : a 10 : a 01 : a 11 ) to (a 00 :

a 10•A T : a 01•A T : a 11•A T ).
This means that U i ′ ,j ′ (0) = 0, where i ′ = i • A and j ′ = 11. We can now assume that j = 11. 4. Now, either i = 00 or i = 11. If i = 11, we are done. Otherwise, we act by the symplectic automorphism with matrix form

    1 0 2 0 0 1 0 2 0 0 1 0 0 0 0 1     .
We then obtain the theta-null point (a 00 : √ -1 • a 10 : √ -1 • a 01 : a 11 ), which means that U 11,11 (0) = 0.

Thus, we can assume we are now working on the product theta structure.

If (a : b : c : d) is theta-null point on E 1 × E 2 , from Proposition 2, it follows that (a : b) is a theta-null point for E 1 and (b : d) is a theta-null point for E 2 . Also, if f (P ) = (P 1 , P 2 ) ∈ E 1 × E 2 is
represented in theta coordinates as (x : y : z : w), we have that (x : y) is the representation of P 1 in theta coordinates for E 1 and (y : w) is the representation of P 2 in theta coordinates for E 2 . Finally, to convert from theta coordinates to the Montgomery model we can use the formulae in [START_REF] Barbulescu | Models of Kummer lines and Galois representation[END_REF].

Computing Isogenies without Extra Isotropic Information

In this subsection, we relax the condition on the two points of order 2 n+2 on E 1 × E 2 above K forming an isotropic group. This does not represent a problem when computing a (2 n , 2 n )-isogeny, except for the two last steps. We discuss two cases: when we can work with 2 n+2 -torsion, and when we cannot. 8Let us discuss the former case.

Let K = ⟨P 1 , P 2 ⟩ ⊂ E 1 × E 2 .
To apply the previous algorithm, we would like to have an isotropic ⟨P ′′ 1 , P ′′ 2 ⟩ above K such that P i = [4]P ′′ i . However, it suffices to pick up any Q ′′ 1 , Q ′′ 2 , not necessarily isotropic, as long as

P i = [4]Q ′′
i . Indeed, one can check that applying the algorithm of section 4.1 on these Q ′′ i gives a theta-null point that differs from the one given by isotropic P ′′ i by an automorphism of the theta group (see Subsection 2.1) induced by a symplectic automorphism,. Hence, it still corresponds to the correct codomain, but with a different theta structure.

In the latter case, we cannot use the 2 n+2 -torsion at all. A way to circumvent this problem is to use square roots to compute the codomains for the last two steps. Once we have the codomain, the image evaluation is unaffected. There is no way to avoid the square root computations: the theta-null point requires a theta structure of level 2, so in particular a full basis of the 2-torsion and some extra information on the 4-torsion. If we don't have the 2 n+2 -torsion at the beginning, we miss the necessary information on the 4-torsion at the penultimate step and on the 2-torsion on the last step. To reconstruct this information requires making choices, hence taking square roots.

At the penultimate step f : A → B, we have T ′ 1 and T ′ 2 of 4-torsion but not the 8torsion points T ′′ 1 and T ′′ 2 anymore. This means that on the codomain, we only have the 2-torsion determined. We have several choices of possible compatible theta structure, but we still want to use the information at hand.

Let (α : β : γ : δ) be the dual theta-null point on B. Applying Equation 2 to the theta-null point (a : b : c : d), we have

H • S(a : b : c : d) = (α 2 : β 2 : γ 2 : δ 2 ). (6)
Also, since f (T ′ 1 ) is in K 1 for the dual theta structure, we have

θB i (f (T ′ 1 )) i = (β : α : δ : γ).
Therefore, from Equation 2,

H • S( θ A i (T ′ 1 ) i ) = (αβ, αβ, γδ, γδ). (7)
Fix α = 1. From Equation 6, we can compute any square root of β 2 for β and any square root of γ 2 for γ. From Equation 7and β we can recover the correct lifting of γδ, and in turn, we can recover δ. The four choices we can make on the square roots of γ 2 and δ 2 describe different theta structures underlying the same abelian surface since they differ by the action of a symplectic automorphism [START_REF] Robert | A note on optimising 2 n -isogenies in higher dimension[END_REF]Example B.3].

At the last step, we only have T 1 and T 2 . As a result, we can only recover the squares (α 2 : β 2 : γ 2 : δ 2 ) of the dual theta-null point (α : β : γ : δ). We can fix α = 1 and compute β, γ, δ via three square roots. Once again, we can check that these 8 choices all come from a valid theta structure [START_REF] Robert | A note on optimising 2 n -isogenies in higher dimension[END_REF]Example B.3].

To sum up, if the 2 n+2 -torsion is available, we need no square root. If the 2 n+1torsion is available, we need two square roots. If only the 2 n -torsion is available, we need 2 + 3 = 5 square roots.

Implementation

We have implemented the computation of an isogeny between elliptic products in the theta model using both the programming language Rust and the computer algebra system SageMath. The SageMath implementation has been designed to follow the API of isogenies between elliptic curves and is intended to be a tool in both experimentation and in constructing proof-of-concept implementations of isogeny-based cryptographic primitives. For those who have previously relied on the SageMath implementation of [START_REF] Oudompheng | A Note on Reimplementing the Castryck-Decru Attack and Lessons Learned for SageMath[END_REF], the function EllipticProductIsogeny(kernel, n) has been designed to be a drop-in replacement for the (2 n , 2 n )-isogeny computed using via the Richelot correspondence and the algorithms presented in [START_REF] Smith | Explicit endomorphisms and correspondences[END_REF].

The Rust implementation has been designed with constructive cryptographic implementations in mind, and in particular, it has been written to be constant time. 9 The finite field arithmetic and certain elliptic curve functions have been adapted from the crrl library [START_REF] Pornin | crrl: Rust library for cryptographic research[END_REF] maintained by Thomas Pornin as well as other ongoing collaborations. An effort has been made to ensure the code is (reasonably) flexible so that without too much tweaking, this work can be ported to other Rust projects. As an example of this flexibility, we show timings of isogenies of various lengths between elliptic products over three distinct base fields.

Both the SageMath and Rust implementations are made available via the following GitHub repository: https://github.com/ThetaIsogenies/two-isogenies

Performance

In this subsection, we include the performance of our algorithm for three distinct isogeny chains between elliptic products over a range of base fields. We include the timings for both the constant-time Rust implementation as well as the proof-of-concept SageMath implementation, together with a comparison to previous work on isogenies between elliptic products in the Mumford model [START_REF] Oudompheng | A Note on Reimplementing the Castryck-Decru Attack and Lessons Learned for SageMath[END_REF] using the optimisations introduced in the implementation of [START_REF] Basso | FESTA: Fast Encryption from Supersingular Torsion Attacks[END_REF]. This triplet of comparisons has a twofold advantage. Firstly, the Rust implementation we present is the first (to our knowledge) constant time implementation of dimension two isogenies between elliptic products. By including the timings of both our Rust implementation and the SageMath implementation, we hope that researchers can estimate a performance gain if they were to write efficient and cryptographically minded implementations following the proof-of-concept scripts which currently exist in the higher-dimensional isogeny-based cryptography literature.

Secondly, our SageMath implementation allows an honest comparison of the isogenies in the theta model to the Richelot isogenies in the Mumford model. We compare against the implementation of [START_REF] Oudompheng | A Note on Reimplementing the Castryck-Decru Attack and Lessons Learned for SageMath[END_REF] together with the additional optimisations introduced for the proof-of-concept of [START_REF] Basso | FESTA: Fast Encryption from Supersingular Torsion Attacks[END_REF] which offered more than a two times speed up by optimising both the arithmetic on Jacobians as well as the isogenies themselves.

The run-times displayed in Table 2 were captured on a Intel Core i7-9750H CPU with a clock-speed of 2.6 GHz with turbo-boost disabled for stable measurements. The Rust code was compiled with the Rust compiler version 1.75.0-nightly with the flag -C target-cpu=native to allow the compiler to use CPU specific opcodes (specifically, mulx for the finite field arithmetic). The arithmetic is written using Rust, rather than optimised assembly for each base field; the inclusion of which would allow dramatically faster results, especially for base fields with large characteristic. This form of optimisation is better suited to particular protocols, and we would expect to see this in optimised implementations of isogeny-based cryptographic primitives.

Comparing our SageMath implementation to the isogeny chain in the Mumford model, we find that the codomain computation is consistently faster by a factor of ten, while the image computation is more than twenty times faster. For smaller characteristic, the Rust implementation is approximately forty times faster than the same algorithm written in SageMath, but this gap closes significantly for larger primes. For example, the FESTA sized parameters run only 2.5 times faster than the SageMath code. Note that the Rust implementation has been written to run in constant time and so the underlying arithmetic between these two implementations is incomparable. 10 We note here that an alternative and faster implementation of (2, 2)-isogenies in the Mumford model is available in [START_REF] Kunzweiler | Efficient computation of (2 n , 2 n )-isogenies[END_REF]. In this work, Kunzweiler uses Jacobians of hyperelliptic curves in specific models which allows (2, 2)-isogeny chains to be computed particularly efficiently. In the initial treatment of this work, isogenies between elliptic products were not considered, leading to FESTA [START_REF] Basso | FESTA: Fast Encryption from Supersingular Torsion Attacks[END_REF] and other projects to rely on [START_REF] Oudompheng | A Note on Reimplementing the Castryck-Decru Attack and Lessons Learned for SageMath[END_REF]. However, Kunzweiler's work can be adapted to the case of isogenies between elliptic products. Additionally, Kunzweiler also has an unpublished SageMath implementation of (2, 2)-isogenies using Kummer surfaces in the Mumford model rather than in the Jacobian model that she kindly provided us. 11 Comparing our results to the computations of Kummer surfaces in the Mumford model is a fairer comparison as we work with level-2 theta coordinates, which are also on Kummer surfaces. Comparing against this implementation, the codomain computation in the theta model is around four times faster than in the Mumford model, and evaluations are around four times faster. We give detailed comparison timings in Table 3.

Although theta coordinates are faster, working in the Mumford model is interesting if the level-2 theta coordinates are not rational, so using theta coordinates would require 10 It is not surprising to see this gap close though, as we expect for very large characteristic that the SageMath overhead becomes negligible compared to the cost of the arithmetic. As such, the comparisons of the two run-times boil down to comparing the Rust finite field arithmetic against the SageMath calls to the optimised arithmetic of the C libraries it is built upon. 11 Kunzweiler's isogenies between elliptic products using both Jacobians and Kummer surfaces are now available via GitHub [START_REF] Kunzweiler | Efficient Computation of (2 n , 2 n )-isogenies[END_REF] Table 3. Comparison of the SageMath running times for a (2 n , 2 n )-isogeny between elliptic products in the theta model against Kunzweiler's implementation in the Mumford model using both Jacobians and Kummer surfaces [START_REF] Kunzweiler | Efficient Computation of (2 n , 2 n )-isogenies[END_REF]. to work on a field extension. Since our domain is a product of elliptic curves, the theta coordinates are rational when each elliptic curve is described by rational theta coordinates. By [START_REF] Barbulescu | Models of Kummer lines and Galois representation[END_REF], this is the case for an elliptic curve E whenever

E[4] = K ′ 1 ⊕ K ′ 2
and the K ′ i are rational subgroups (equivalently they are generated by a 4-torsion point T ′ i which is rational on the Kummer line E/ ± 1). In particular, a supersingular elliptic curve E/F p 2 always has rational level-2 theta coordinates (assuming p odd). However, if E is a supersingular elliptic curve defined over F p , the theta coordinates are rational only if E lies at the top of the 2-isogeny volcano.

Comparison with Dimension One

In Table 4, we provide a timing comparison using SageMath between a 2 n -isogeny in dimension one using the efficient formulae of [START_REF] Renes | Computing isogenies between montgomery curves using the action of (0,0)[END_REF] to our formulae in dimension two over the same base field. The dimension two isogeny has degree 2 2n so is expected to be slower. Our timings show a consistent factor-two slow down both for the codomain and image computations in dimension two compared to dimension one. This is essentially the best we could hope given the degrees, and actually better than expected.

The dominating costs of a 2 n -isogeny are the intermediate doublings and images. In the following we consider the more costly doublings and images in dimension two which arise from avoiding inversion when computing the codomain. First of all, we have around twice as many doublings and images in dimension two than in dimension one beause the kernel is of rank 2. The cost of doubling in dimension one is 4M + 2S compared to 8M + 8S in dimension two, and an image is 4M compared to 4M + 4S in dimension two. Thus, while images are twice slower, doublings are around 2.5× slower, and the intermediate codomain computations are also slower. Furthermore, a lot of doublings are done on the first step of the chain to get the first kernel, so on the elliptic product.

While it might seem at first glance that these doublings would only occur a twofold slowdown, in practice, for the gluing images, we need to compute these points in affine (x, y) coordinates rather than x-only coordinates to allow access to addition laws 12So, all in all, we should expect a slowdown around 4× for perfect implementations. Our benchmarks show a slowdown slightly less than 2×, making two-dimensional iso-genies perform better than expected (by contrast, the 2× slowdown for images is consistent with the theory). This is probably due to SageMath overhead and the fact that the dimension one implementation has been designed to allow arbitrary degree rather than only chains of two isogenies and is missing some optimisations. A final caveat is that in dimension one, it is faster to split the 2 n -isogeny using the fast 4-isogenies from [6] rather than using 2-isogenies -we did not do that in our comparison because we do not have efficient 4-isogeny formulae in dimension two yet. Still, taking into account the degrees of the respective isogenies, this shows that our dimension two formulae are quite competitive with the best dimension one formulae.

Implementation details

In this subsection, we explain two optimisations we applied in the implementation. The first one is a direct consequence of Remark 7, where we describe how to lower the complexity of the codomain computation by reusing some constants. The second optimisation consists in the application of optimal strategies [10] to our case.

Reduce, Reuse, Recycle A simple and obvious optimisation is to reuse as many computations as possible throughout the isogeny chain. As mentioned in Remark 7, for each step on the isogeny chain, we can precompute six field elements at a cost of 4S + 21M + 1I. Knowledge of these values allows the doubling of any theta point on the corresponding theta structure to have a cost of 8S + 6M, but it also allows the following codomain computation cost to be lowered to 8S + 9M + 1I.

In a similar way, by computing the codomain together with all the images we require at each step, we can precompute two additional field elements which can be reused for each evaluation, bringing the cost of an image to half that of doubling: 4S + 3M.

For the gluing isogeny, the basis change is determined from the kernel and so cannot be precomputed. However, the last step at the end of the isogeny chain requires to find a symplectic transformation that maps the zero even index to the position [START_REF] Gaudry | Fast genus 2 arithmetic based on theta functions[END_REF][START_REF] Gaudry | Fast genus 2 arithmetic based on theta functions[END_REF]. As there are only ten even indices, we can precompute ten symplectic transforms which map any given zero even index to [START_REF] Gaudry | Fast genus 2 arithmetic based on theta functions[END_REF][START_REF] Gaudry | Fast genus 2 arithmetic based on theta functions[END_REF]. Computing the basis change is then only a matter of finding the the current zero index and from this, selecting the precomputed matrix and applying the transformation.

On Inversions At each step of the isogeny chain, we compute one inversion for the intermediate codomain. This inversion allows us to reduce the cost of doublings on this codomain from 8M + 8S to 6M + 8S and the cost of images from 4M + 4S to 3M + 4S. However, at the end of the isogeny chain, there remain fewer doublings and images to compute, so it would be more efficient to skip this inversion and occur the higher cost. The precise cutoff point would depend on the relative cost of the inversion compared to a multiplication and the length of the chain. This optimisation has not yet been implemented in our code.

On Square Roots As explained in Section 4.2, when we do not have the 2 n+2 -torsion available, we need to compute some square roots at the end of the chain (five square roots in total). This only changes the computation cost of the last two codomains, and do not affect the images computations. The longer the isogeny chain, the less impactful these square roots will be. Benchmarking our SageMath implementation, we observed that the impact of these five square roots is completely negligible for the chains we consider.

Optimal Strategies

As is now standard with computing long isogeny chains, we can reduce the complexity of isogeny chains from a quadratic number of edges in the graph of doublings and evaluations to quasi-linear following the "optimal strategies" introduced in [START_REF] Feo | Towards quantum-resistant cryptosystems from supersingular elliptic curve isogenies[END_REF]. Essentially, the saving comes from reducing the total number of doublings when computing the kernel for each step in the chain by pushing through intermediate points encountered in the repeated doubling. For isogenies in the theta model, the cost of images is half that of doubling, and so shifting the cost in this way is particularly useful in optimisations.

Although this strategy was first discussed in dimension one for the case of isogenies between elliptic curves, using it in dimension two is a natural generalisation. For the dimension one case, the strategy is computed from balancing the costs of doubling and evaluating the kernel generator through the chain. In dimension two, the (2, 2)isogeny is generated by a pair of elements which means twice the number of evaluations, but as the pair of elements must also be doubled to obtain the kernel for each step, essentially nothing changes. The cost weighting for the optimal strategies is a ratio between doublings and evaluations, which means we can naively use an identical method as described in [START_REF] Feo | Towards quantum-resistant cryptosystems from supersingular elliptic curve isogenies[END_REF] to compute a strategy for our isogeny chain.

Implementing the strategy with the weighting of doublings and images at a cost of (2 : 1), we find an approximate ten times speed up in comparison to an implementation with no strategy. Concretely, for the Rust implementation of the isogeny chain of length n = 208, we see a speed up from 107 ms to 11.4 ms.

However, unlike the isogeny chains between elliptic curves, the isogeny chain between elliptic products in our implementation does not have the same costs for every step. For steps in the chain between generic theta structures, the cost weighting is indeed (2 : 1). However, for the first gluing isogeny, doubling an element on the product structure has a cost of 12S + 12M while the cost for the image is much more expensive.

To compute the image of a point P ∈ E 1 ×E 2 one must first compute the shift P +T ′ 1 for a cost of 10S + 32M to projectively add a pair of points. Then, for each of these two points on the product, there is a cost of 4M to compute the corresponding theta point from elliptic curve coordinates, and an additional 16M required perform the matrix multiplication for the basis change to ensure a compatible representation. Altogether, this precomputation costs 10S + 72M. Given the theta point corresponding to the pair of points on the product structure, there is still then the final cost of 8S + 5M + 1I for the special image itself. Furthermore, for this to be implemented in constant time, both branches depending on whether a coordinate is zero or not must be evaluated, raising the practical cost to 8S + 10M + 1I.

All things considered, a gluing image costs 18S+82M+1I, making it approximately seven times the cost of the doubling for this first step and fourteen times the cost of a regular image. Visualising the graph of doublings and images as in [START_REF] Feo | Towards quantum-resistant cryptosystems from supersingular elliptic curve isogenies[END_REF]Figure 2], this means we must weigh the cost of moving down the left most branch with the product doubling and the first step right from the leftmost branch with this high-cost gluing image.

Taking this into account, an optimised strategy for the isogeny between elliptic products for our formula must be tweaked from the original case to find the right balance of doublings and expensive images from this left branch. Applying this modification, we are able to find the "proper" optimised strategy, which further reduces the run-time of the isogeny chain computation by approximately 2%. 13 For the same chain as above, we see a computation time improve from 11.4ms to 11.2ms. For an explicit description for computing the optimised strategy with a different costs on the left-most branch, see the implementation.

An Application: FESTA

As an explicit, cryptographic example of the new isogeny formula, we can take our implementation and use it to compute the isogeny between elliptic products which is required within the decryption algorithm of the isogeny-based public key encryption protocol FESTA-128 [START_REF] Basso | FESTA: Fast Encryption from Supersingular Torsion Attacks[END_REF]. Concretely, this requires computing an isogeny of length n = 632, where the base field has a characteristic with log p = 1293 bits, and the evaluation of a pair of points on the elliptic product L 1 = (R 1 , R 2 ) and L 2 = (S 1 , S 2 ), L i ∈ E 1 ×E 2 .

A direct swap from the isogeny chain derived from the Richelot correspondence used in the FESTA proof-of-concept would require using Section 4.2 to compute the final two steps without the eight-torsion above the kernel. An implementation of this is available in SageMath, but for the purpose of FESTA, we instead propose to tweak the 128-bit parameter set to instead allow for the additional torsion information to be known, allowing the isogeny chain to be computed as fast as possible while only including an additional 2-bits in the masked torsion data.

We find that the our SageMath implementation of the codomain computation has a ten times speed up compared to the proof-of-concept code accompanying [START_REF] Basso | FESTA: Fast Encryption from Supersingular Torsion Attacks[END_REF], and evaluating the pair of points is now thirty times faster. As a hint to what approximate 13 As an aside, in the original discussion of the optimised costings, it is shown that a 2-3% improvement is gained by moving from a balanced to optimised strategy. Seeing a similar saving from the naive (2 : 1) weighted optimisation to one carefully handling the cost of the gluing step is then within our expectations.

running times may be for FESTA, computing the codomain and both images using our Rust implementation takes only 563ms, a 2.5 times speed up over the SageMath implementation. Note that these computation are precisely that of the final row of Table 2. We again note that further optimisations of the finite field arithmetic could offer substantial speed ups, as seen in the optimised assembly implementations for large characteristic implementations of SIDH [12, Table 2.1] and the efficient algorithms of [START_REF] Longa | Efficient Algorithms for Large Prime Characteristic Fields and Their Application to Bilinear Pairings[END_REF]. Concretely, we propose the following parameters for FESTA-128 using the same notations as in the original paper b := 632, In practice, we only tweak the cofactor f in the prime p = 2 b d 1 (d A,1 d A,2 ) sf d 2 f -1 to be f = 0 (mod 4). 14Within the context of the proof-of-concept decryption times (in SageMath), the novel algorithms we present here offer a four times speed up in decryption, with runtimes for FESTA-128 being reduced from 20.7s to only 5.4s. 15 We see that through computing the dimension two isogeny in the theta model, the time spent computing a (2 n , 2 n )-isogeny shrinks from 70% of the run-time to only 25%, with the remaining computation time spent in dimension one, computing various discrete logs and Weil pairings to complete the decryption routine.

Conclusions

In this paper, we have described and implemented formulae to compute (2 n , 2 n )isogenies between elliptic products in the theta model. The main goal was to provide a comprehensive and self-contained treatment of the theta model, specialising to the two-dimensional case. The extension of this work to dimension four is therefore a natural consequence, which will greatly benefit isogeny-based cryptography (see [START_REF] Dartois | SQISignHD: New Dimensions in Cryptography[END_REF] for instance), as well optimising the dimension two case for other degrees.

Our algorithm significantly outperforms the previous method in [START_REF] Oudompheng | A Note on Reimplementing the Castryck-Decru Attack and Lessons Learned for SageMath[END_REF][START_REF] Basso | FESTA: Fast Encryption from Supersingular Torsion Attacks[END_REF]: in SageMath, the codomain computation is ten times faster, while the isogeny evaluation is more than twenty times faster. The implementation in Rust has been written to run in constant time, with cryptographic implementations in mind. It runs up to forty times faster than the same algorithm written in SageMath.

We tested our algorithm on the proof-of-concept implementation in [START_REF] Basso | FESTA: Fast Encryption from Supersingular Torsion Attacks[END_REF] and showed a fourfold speed up in decryption, highlighting that the slowest part is now given by the computations in dimension one. Furthermore, our SageMath implementation has been designed to allow protocols whose implementation relies on the previous proofof-concept in [START_REF] Basso | FESTA: Fast Encryption from Supersingular Torsion Attacks[END_REF] to be easily upgradeable, allowing the theta-model code to be used in many more projects without too much work. Ultimately, the aim is to provide a new tool to facilitate research in higher-dimensional isogeny-based cryptography, allowing us to better understand the practical role of higher-dimensional isogenies in constructive applications.

If

  P = (x : y : z : t), we have S 1 = (b : a : d : c) and P + S 1 = (y : x : t : z), S 2 = (c : d : a : b) and P + S 2 = (z : t : x : y), T 1 = (a : -b : c : -d) and P + T 1 = (x : -y : z : -t), T 2 = (a : b : -c : -d) and P + T 2 = (x : y : -z : -t).

  we can actually compute the inverse of the projective factor by θA 01 (f (P ))/(λ θA 01 (f (P ))). Otherwise, we repeat the same process with the second last component of H • S((θ E1×E2 i (P + T ′ 1 )) i ). Once we have λ, we extract θA 00 (f (P )) from the second component of H•S((θ E1×E2 i (P + T ′ 1 )) i ): we multiply the second component of H • S((θ E1×E2 i (P + T ′ 1 )) i ) by λ -1 • β -1 . Finally, we obtain the image of the point P under f via H( θA 00 (f (P )), θA 10 (f (P )), θA 01 (f (P )), θA 11 (f (P ))).

d 1 :

 1 = (3 3 • 19 • 29 • 37 • 83 • 139 • 167 • 251 • 419 • 421 • 701 • 839 • 1009 • 1259 • 3061 • 3779) 2 , d 2 := 7 • (5 2 • 7 • 11 • 13 • 17 • 41 • 43 • 71 • 89 • 127 • 211 • 281 • 503 • 631• 2309 • 2521 • 2647 • 2729) 2 , d A,1 := (59 • 6299 • 6719 • 9181) 2 , d A,2 := (3023 • 3359 • 4409 • 5039 • 19531 • 22679 • 41161) 2 , m 1 := 1492184945093476592520242083925044182103921,m 2 := 25617331336429939300166693069, f := 4 • 71.

  [START_REF] Basso | FESTA: Fast Encryption from Supersingular Torsion Attacks[END_REF] -from now on, we will drop the dependence on the line bundle L and simply writeA[2] = K 1 ⊕ K 2 . Let ⟨S 1 , . . . , S g ⟩

	be the canonical basis induced by Θ	L (K(2, . . . , 2)) and ⟨T 1 , . . . , T g ⟩ the canonical basis
	induced by Θ

L ( K(2, . . . , 2)), which means K 1 = ⟨S 1 , . . . , S g ⟩ and K 2 = ⟨T 1 , . . . , T g ⟩. Now, let us consider the isogeny f : A → B, where ker(f ) = K 2 . The abelian variety B is principally polarised and in turn can be endowed with a type-(2, . . . , 2) theta structure, whose theta coordinates are denoted by (θ B i ) i . Also, let us define ⋆ to be the operator such that

  The points (P ′ 1 , P ′2 ) and(Q ′ 1 , Q ′ 2 )in the four-torsion of E1 × E2 below the kernel. Output: The 4 × 4 change of basis matrix N .

	1: G1 ← action_by_translation(P ′ 1 )	(▷) Algorithm 1:
		we can
	use Equation 2, but first we need to recover θB i (f (P ))

i and θB i (f (Q)) i , which can Algorithm 2 Change of Basis Input:

  ). Output: The theta coordinates (xP +Q : yP +Q : zP +Q : wP +Q) of P + Q. 1: XP , YP , ZP , WP ← H • S(xP , yP , zP , wP ) 2: XQ, YQ, ZQ, WQ ← H • S(xQ, yQ, zQ, wQ) 3:

Table 2 .

 2 Running times of computing the codomain and evaluating a (2 n , 2 n )-isogeny between elliptic products over the base field F p 2 . Times were recorded on a Intel Core i7-9750H CPU with a clock-speed of 2.6 GHz with turbo-boost disabled.

				Codomain			Evaluation
			Theta	Theta	Richelot	Theta	Theta	Richelot
	log p	n	Rust	SageMath SageMath [27]	Rust	SageMath SageMath [27]
	254 126 2.85 ms	108 ms	1028 ms	161 µs	5.43 ms	114 ms
	381 208 11.2 ms	201 ms	1998 ms	411 µs	8.68 ms	208 ms
	1293 632 495 ms	1225 ms	12840 ms	17.8 ms	40.8 ms	1203 ms

  Times were recorded on a Intel Core i7-9750H CPU with a clock-speed of 2.6 GHz with turbo-boost disabled.

				Codomain			Evaluation
	log p	n	Theta	Jacobian Kummer Theta Jacobian Kummer
	254 126 108 ms	760 ms	467 ms	5.43 ms 66.7 ms	18.4 ms
	381 208 201 ms	1478 ms	858 ms	8.68 ms	119 ms	31.4 ms
	1293 632 1225 ms 9196 ms 5150 ms 40.8 ms	593 ms	170 ms

Table 4 .

 4 Comparison of the running times for a 2 n -isogeny in dimension one and dimension two over the same base field. Times were recorded on a Intel Core i7-9750H CPU with a clock-speed of 2.6 GHz with turbo-boost disabled.

			Codomain	Evaluation
	log p	n	Montgomery	Theta	Montgomery Theta
	254 126	63 ms	108 ms	2.24 ms	5.43 ms
	381 208	136 ms	201 ms	4.4 ms	8.68 ms
	1293 632	727 ms	1225 ms	20 ms	40.8 ms

Technically, eL maps onto G(L). However, the first component of eL(x, y) is the identity on the abelian variety, whereas the second one is an automorphism of L, and the group of such automorphisms can be identified with k * .

The subscript ij refers to the pair([i], [j]) ∈ K(2, 2).

For instance, it is preferable not to work with the 2 n+2 -torsion when it is defined over a field extension of the base field.

The implementation assumes the kernel generators are good with respect to them generating an isogeny between elliptic products. Designing the algorithm to run in constant time with malformed input extends beyond the goals of this paper but may be necessary for protection against side-channel attacks against schemes which rely on this algorithm.

We could also use differential additions to compute [m]P, [m]P +T ′ 1 , but this would be more expensive than just doubling in the affine model.

The notation t sf refers to the square-free part of the integer t.

The speed up here is more than just a faster isogeny chain, but also some other speedups which come for free, as the previous implementation required further isomorphisms to ensure the codomain product curves were in the Montgomery model, but for the theta-model isogeny, we get this for free.

Acknowledgments. Huge thanks are given to Thomas Pornin for his advice and previous collaborations, both of which were instrumental in the design and implementation of the constant time Rust implementation. We also thank Sabrina Kunzweiler for fruitful discussion.