Sébastien Andrieux

Pierre Muller

Manish Kaushal

Nadia Sofía

Macias Vera

Robin Bollache

Clément Honorez

Alain Cagna

Wiebke Drenckhan

Microfluidic Thin Film Pressure Balance for the Study of Complex Thin Films

Investigations of free-standing liquid films enjoy an increasing popularity due to their relevance for many fundamental and applied scientific problems. They constitute soap bubbles and foams, serve as membranes for gas transport or as model membranes in biophysics. More generally, they provide a convenient tool for the investigation of numerous fundamental questions related to interface-and confinement-driven effects in soft matter science. Several approaches and devices have been developed in the past to characterise reliably the thinning and stability of such films, which were commonly created from low-viscosity, aqueous solutions/dispersions. With an increasing interest in the investigation of films made from strongly viscoelastic and complex fluids that may also solidify, the development of a new generation of devices is required to manage reliably the constraints imposed by these formulations. We therefore propose here a microfluidic chip design which allows for the reliable creation, control and characterisation of free-standing films of complex fluids. We provide all technical details and we demonstrate the device functioning for a larger range of systems via a selection of illustrative examples, including films of polymer melts and gelling hydrogels.

Introduction

Free-standing liquid films are at the core of numerous scientific and applied problems. Since they are thermodynamically metastable, their rupture is commonly delayed by the addition of surface-active agents (low molecular weight surfactants, lipids, amphiphilic polymers, proteins, particles, etc.), which leads to stabilising interactions on the film surface ("interfacial rheology" [1][2][3][START_REF] Pagureva | Colloids Surfaces A Physicochem[END_REF][5][6] ) as well as across the film ("disjoining pressure" 7 ). Thanks to this increased stability, extended liquid films of surface areas in the range of mm 2 -m 2 and thicknesses of only 10-1000 nm can be generated. Such free-standing liquid films arise most frequently in the form of soap bubbles and in foams, the latter being essentially threedimensional arrangements of free-standing films 8 . Due to their large surface-to-volume ratio, liquid films are also increasingly used as membranes for gas separation [9][10][11] . Apart from practical applications, freestanding liquid films provide excellent model systems for fundamental investigations concerning interface and/or confinement-driven phenomena in soft matter science [12][13][14][15][16][17][18][19][START_REF]Foam and Foam Films, Theory, Experiment, Application[END_REF] .

Due to the above-mentioned importance of free-standing liquid films it is necessary to avail of advanced and easy-to-handle tools to characterise their properties reliably in terms of thinning dynamics and stability. Let us assume that a difference is applied between the pressure pg of the gas surrounding the film and the pressure pl of the liquid in the film, as sketched in Figure 1a. This pressure difference drives film thinning ("drainage") and is commonly called the capillary pressure: pc = pg -pl. In many cases, it arises naturally due to gravity. The main four stages of the evolution of such a film are then:

(1) The drainage of an initially thick film (typically with a film thickness h > 100 nm) is fully governed by hydrodynamics. The viscoelastic monolayer of stabilising agents at the film surface imposes boundary conditions on the flow 1,17,[START_REF] Champougny | [END_REF] . (2) The film becomes thin enough (h < 100 nm) so that the two monolayers of stabilising agents start to interact via DLVO forces (electrostatics, van der Waals, steric…) across the film 15 , creating a disjoining pressure  (Figure 1b), which can either promote ( < 0) or hinder ( > 0) drainage.

The confinement may also create structural forces in the bulk of the film, leading, for example, to stratified drainage [22][23][24][25][26] .

(3) When  = pc, an equilibrium is reached, and film drainage stops. Different capillary pressures pc lead to different film thicknesses, allowing therefore to establish disjoining pressure isotherms (h) to quantify the DLVO interactions in the film 15 . (4) After a certain time and/or beyond a certain disjoining pressure , the film ruptures 27,28 .

The initial drainage stage is often investigated via gravity-driven drainage 1,17,[START_REF] Champougny | [END_REF] of vertical films 1,[START_REF] Champougny | [END_REF]29 or on individual bubbles 30 . However, this renders an explicit control of the capillary pressure difficult and other effects become important (marginal regeneration, stability of the film under its own weight, etc.). Therefore, a whole class of devices -so-called "Thin Film Pressure Balances" (TFPB) 15 -has been developed to investigate horizontal films under an explicitly applied capillary pressure. The main functioning of all TFPBs is identical: a controlled capillary pressure is applied to drain the film down to its equilibrium state (and potentially rupture). Optical interference techniques with white or monochromatic light are used to determine the film thickness (profile or point measurement) during drainage as well as at equilibrium. However, the devices differ in (i) how the film is held, and (ii) how the capillary pressure is applied.

Concerning (i), the historical Scheludko cell (Figure 1c.I) 31-33 contains a millimetric film at ambient gas pressure whose liquid pressure is controlled via a liquid channel. The maximum capillary pressure that can be reached is of the order of the entry pressure of the gas into the liquid channel, i.e. inversely proportional to the channel width. This leads to characteristic capillary pressures of the order of 10² Pa, which is often too low to investigate relevant phenomena. Moreover, the asymmetric feeding of the film can lead to non-negligible dynamic effects. Therefore, Mysels proposed in 1964 the use of a porous glass ring 34,35 with very fine pores (hence a large entry pressure) to hold the film. This approach was later improved by Exerowa and Scheludko 36 (Figure 1c.III) via the use of a glass disk, whose conical hole holds the film solidly in place without the formation of a large meniscus. This Exerowa-Scheludko cell has since then become one of the standard techniques, allowing to work with very homogeneous film feeding at capillary pressures up to many thousands of Pa 37,38 . However, the use of porous discs has also disadvantages, such as adsorption phenomena of the surfactants due to the high surface to volume ratio of the porous disc, clogging, and difficult cleaning. This is why Radke 39 proposed in 2001 a compromise between both film holders in form of a "bike wheel" (Figure 1c.II) feeding the film with 24 evenly distributed lithography-produced channels of less than 100 µm thickness. The same approach was successfully applied in 2016 by Beltramo et al. 40 Concerning (ii), many different approaches exist to apply a well-controlled capillary pressure in the chosen film holder. The simplest approach is to work at ambient gas pressure and to control the liquid pressure via a hydrostatic pressure feed. More sophisticated approaches control the liquid volume with a syringe pump and measure the resulting pressure, or use directly a pressure controller 25,[41][42][START_REF] Khristov | Colloids Surfaces A Physicochem[END_REF] . Other techniques apply a constant hydrostatic pressure to the liquid and control explicitly the gas pressure 22,29,[START_REF] Wang | Minerals Engineering[END_REF] .

The described approaches have been proven very useful in advancing our understanding of free-standing films generated from solutions/dispersions of reasonably low viscosities [12][13][14][15][16][17][18][19][START_REF]Foam and Foam Films, Theory, Experiment, Application[END_REF] . However, TFPBs are increasingly used for the investigation of strongly viscoelastic, complex fluids, including fluids undergoing solidification/gelation. This puts new constraints on the technique, requiring a highly flexible and rapid management of high pressures and the ease to clean the device after use. Other constraints are added naturally, such as the need to replace the gas phase on either side of the film or to replace the gas by a liquid phase to investigate liquid/liquid films.

In order to find a solution for such a versatile device, we developed a microfluidic TFPB ("TFPB"), whose core is a microfluidic Radke-style bike wheel with short feeding channels entirely integrated into a threedimensional millifluidic chip. The geometrical flexibility in the design of this part allows to adapt the microfluidic channel dimensions (width and length) to find the appropriate compromise between two conflicting parameters: maximum entry pressure and minimum flow resistance. One of the main advantages of this chip is that it is composed of several layers, containing different ensembles of channels and windows. These layers can be separated and re-assembled easily for cleaning and repeated measurements. The chip geometry and layer assembly are described in detail in Section 2.1. Both, gas and liquid pressure are explicitly set using computer-controlled pressure controllers, as detailed in Section 2.2. In contrast to existing devices, this provides a highly responsive and flexible control of the applied capillary pressure with a variable pressure offset. Section 2.3 presents the optical system used to visualise and quantify the film thickness, which allows for the study of film drainage.

The goal of this article is to provide to the community a detailed description of the developed TFPB and show the range of fluids accessible for study using a selection of illustrative examples, including a simple surfactant solution (Section 3.3.1) and a polymer melt (Section 3.3.2). We also show for the first time the use of a TFPB for solidifying hydrogel films (Section 3.3.3).

Design and functioning of the microfluidic Thin Film Pressure Balance

The main innovation of the microfluidic thin film pressure balance (µTFPB) is the microfluidic chip. It is integrated into the overall set-up as sketched in Figure 2a andb. The set-up can be divided into three main parts: the microfluidic chip (Section 2.1), the fluid control (Section 2.2) and the optical system (Section 2.3).

Figure 2: a) General setup of the microfluidic thin film pressure balance: a multi-channel pressure controller applies a well-defined capillary pressure to a free-standing film maintained by a bike-wheel type cone in the centre of the microfluidic chip. The film is visualised (and its thickness analysed) using interference of white light reflected from the film and detected by a CCD camera. b) Three-dimensional view of the chip and the optical system.

Microfluidic chip

The microfluidic chip constructs a three-dimensional channel network to create, control and analyse a millimetric, horizontal liquid film created in its centre. The chip is composed of four plates (5x5 cm, numbered from A to D, in Figure 3a), each one being independently milled in COC (Cyclic Olefin Copolymer, commercial name: TOPAS 8007X4, 4 mm-tick wafers, from Microfluidic Chip Shop) with a micromilling machine (SLS Micro Mill from Minitech). Erreur ! Source du renvoi introuvable. in the supplementary materials summarises the tools used to mill the different plates constituting the chip, and the parameters used for milling. The channel geometry is designed using Inventor® (computer-aided design software from Autodesk). The software Esprit® (computer-aided manufacturing software from DP Technology) is used to translate the design into the final data format used for the milling process.

Gas-and liquid-tight sealing between the different plates of the chip is ensured by commercial rubber joints (1.6 mm thick, black rings in Figure 3c,d) positioned in specifically designed channels (width: 2.00 mm, depth: 1.28 mm). The plates are held together tightly by four bolts and nuts. The two outer parts of the chip (A and D) are milled to obtain an inclined surface with a 5° angle to avoid any parasite reflection from the chip that could disturb image analysis. 3 mm wide holes are also drilled into these plates above and below the place where the film forms to avoid any loss of intensity and leave a clear path for the light. The inclined, drilled surfaces, one of which can be seen on plate D in Figure 4a, are sealed by gluing a microscope square cover slip (see Figure 3c,d). Alignment issues of the plates may arise, yielding disturbing shades and reflections upon observation of the film. In such a case, it may be useful to add guiding bolts into holes drilled at the exact bolt diameter using the micromilling machine.

The film is formed between the parts B and C, at the intersection of two cones (see Figure 3b andf), so that the film has a diameter of 1 mm. The position of the bike wheel is pointed at by the green arrow on plate B in Figure 3a. The liquid arrives at the hole between the two cones from the side and enters the chip from the bottom of plate A, from the bottom. The gas is injected in the inlet on plate D and circulates on both sides of the film via connecting channels between the different parts of the chip (see the blue arrows in Figure 3b), ensuring an equal pressure on both sides of the film.

The liquid film forms within a bike wheel-type cell (see Figure 3f), whose geometry we modified with respect to the one proposed by Radke 39 to optimise for the study of viscous liquids at high capillary pressures. The film is fed homogeneously via six channels with an entry width of 259 µm and a height of 20 µm. The small channel height ensures an entry pressure of the order of 10000 Pa. Moreover, the channels are only 1 mm long to reduce the pressure drop required to push the liquid. The channels are fed homogeneously by a much larger and deeper ring (1 mm wide and 500 µm deep). All channels and tubing leading to the bike wheel are at least 1 mm wide. Hence the main pressure drops occurring when pushing the liquid through the chip arises in the channels of the bike wheel. This part of the chip can be easily adjusted for different liquids in order to find the appropriate compromise between maximising the entry pressure and minimising the hydrodynamic resistance of the flow. Another important advantage of this plate configuration in contrast to previous designs is that it can be opened easily in a manner that all channels are fully accessible for cleaning. The chip can be easily re-assembled for repeated experiments as long as liquids are used for the experiments and the cleaning liquid do not damage the COC.

Fluid control

Both the liquid and the gas pressure are controlled via a multi-channel pressure pump OB1 Mk3 from Elveflow with a pressure range of 0 to 2000 Pa (Figure 2a). In contrast to previous device designs, this provides a highly responsive application of the capillary pressure with a tuneable pressure offset. The latter controls the pressure difference with the atmospheric pressure, which can be an important parameter for pressure-sensitive applications (gas mixtures, etc.). If desired, a pressure controller with a higher maximal pressure can be used (at the expense of precision in the low-pressure range). We tested that the chip remains perfectly sealed up to 10 5 Pa by observing an absence of bubble formation when the chip was held under water. Beyond this range, the capillary pressure overcomes the entry pressure, leading to the entry of gas into the liquid channel. A full characterisation of films needs to combine pressure controllers of different ranges to access a wide pressure range with the required precision, especially in the low-pressure range.

The liquid is stored in a bottle with a GL45 cap linked on one side to the pressure controller and on the other side to the microfluidic chip. The liquid is pushed into the chip by applying a pressure on the bottle. A similar bottle is also inserted between the gas outlet of the pressure controller and the gas inlet of the microfluidic chip. The use of this bottle is three-fold. (1) It protects the pressure controller from being flooded by the liquid in case of a sudden pressure increase at the liquid inlet. (2) It allows for saturating the atmosphere within the chip with the solution by depositing the liquid at the bottom of the bottle [START_REF] Stubenrauch | [END_REF] . The chip is flooded with humid gas before plugging the liquid channel connector in. (3) It homogenises pressure fluctuations.

Both the liquid and gas pressures are set using the Elveflow Smart Interface software (from Elveflow), providing straightforward control over the capillary pressure and hence the film drainage and equilibrium thickness. However, one must pay attention to the liquid pressure pl, which is not equal to that set by the pressure controller. One needs to add the hydrostatic pressure difference pH = gH, where H is the height difference between the film and the surface of the solution in the bottle.  is the density difference between the solution and the gas, and g is the acceleration of gravity. To avoid this correction, it is convenient to place the bottle on a lab jack to align the level of liquid in the bottle with the centre of the microfluidic chip, i.e. H = 0.

For solutions with low viscosity, one may simply use the hydrostatic pressure to fill the chip and form the film without having to apply an additional pressure with the pressure controller. Changing the height of the liquid bottle suffices to form a soap film similar to the classical TFPB. Yet, for highly viscous solutions, the hydrostatic pressure and the pressure applied by the pressure controller may not be enough to push the liquid into the chip. One may also use a syringe pump instead to push the liquid and form the film. This requires, however, the integration of an additional pressure sensor to monitor the liquid pressure pl.

Optical System and film thickness measurement

The analysis of the film thickness is done via interferometry (Figure 2b). The procedure described here is fairly standard for a thin film pressure balance and has been validated over the last 50 years of research 46 . In our device, the film is homogeneously exposed to white light from a halogen lamp (15 V, 150 W). The lamp illuminates the chip via an optical fibre (shown in Figure 2b) directed towards a beam splitter which redirects the light towards the chip, as schematised by the yellow arrow in Figure 2a. A 4x objective focuses the light to the area where the film forms. The reflected light passes again through the objective and the beam splitter before being collected by a 12 bit CCD camera (UEye UI-3580LE-C-HQ from iDS Imaging Development Systems GmbH) with a CMOS colour sensor and a resolution of 4.2 Mpx. The ensemble is constructed using Thorlabs elements combined with 3D-printed pieces to hold the light source and the microfluidic chip. Alignment of the film with the optical path is done using a combination of a mechanical translation and rotation stage. A monochromatic filter can be used in front of the camera to obtain monochromatic images for thickness measurements. However, thanks to the quality of modern cameras it is possible to profit from the RGB splitting of unfiltered images instead without loss in precision [24][25][26] . Moreover, it has the advantage of availing of fully coloured images for general visual inspection of the film in parallel to the possibility of quantitative data analysis.

For this purpose, the images are treated with the open-source software ImageJ. First, we split the channels using the RGB splitting tool. This procedure provides three grey-scale images. All three colours can be used for treatment and even be combined for better precision [24][25][26] . We treat here only the red channel, since it provides the best precision due to its longer wavelength. This choice needs to be reconsidered, however, if strongly absorbing liquids are used. We use λ = 600 nm for the wavelength since it corresponds to the maximal sensitivity of the camera for the red channel.

Interference in liquid films is periodic with /2nl with the first constructive maximum in reflected intensity 𝐼 r max arising at h = /4nl (113 nm for water), where nl is the refractive index of the liquid. Here we are interested in quantitative measurements only for film thicknesses below this value.

We select a circular region of interest (ROI) on the film with an area of 1024 pixels, in which we measure the average grey value 𝐼 r . This operation is repeated on all the images recorded for the film using an ImageJ macro. We calculate the film thickness h at each time t for the same ROI using the Scheludko equation 46

ℎ = λ 2𝜋𝑛 l arcsin √ ∆ 1+ 4𝑅 (1-𝑅) 2 (1-∆)
, with λ the wavelength, 𝑅 the reflexion coefficient given by 𝑅 = [(𝑛 l -𝑛 g ) (𝑛 l + 𝑛 g ⁄ )]², where 𝑛 g is the refractive index of the gas. ∆ is the relative intensity which is defined using the maximal reflected intensity 𝐼 r max and the minimal reflected intensity 𝐼 r min ∆= 𝐼 r -𝐼 r min 𝐼 r max -𝐼 r min .

The maximal reflected intensity 𝐼 r max is obtained during the early stages of film thinning at h = /4nl, when the film passes through fully constructive interference, while the minimal reflected intensity 𝐼 r min is measured when the film is broken. 𝐼 r min is a collection of residual intensity resulting from the environment and different reflections within the set-up etc. Since this is a constant contribution, it can be subtracted on all images.

For reliable image analysis, the camera settings (in particular the exposure time) is chosen such that the obtained grey values cover a maximally wide range for best precision. Moreover, the gray values should be in a range between 50 and 200 where the captor sensitivity is highly linear.

Examples

In this section we demonstrate examples of thin films studied with the µTFPB. The first example (Section 3.3.1) is a film made of an SDS (sodium dodecylsulfate) solution at a concentration of 0.24 mM (i.e. 3 times the critical micellar concentration, without any added salt). The viscosity of the SDS solution is that of water (10 -3 Pa s). The second example (Section 3.2) is a film made from a polymer melt (DBP-732 from Gelest) consisting of comb polymers with a polydimethylsiloxane backbone with side chains of poly(ethylene glycol)/poly(propylene glycol) copolymers). DBP-732 has a significantly higher viscosity (2.7 Pa s 47 ) compared to the SDS solution, which makes the study of films very difficult with a classical TFPB, as discussed in Section 1. Finally, in Section 3.3 we show that the µTFPB allows us to monitor the drainage and stability of a film during solidification, using the example of a physically gelling alginate hydrogel. The examples chosen are illustrative of the range of complex fluids which can be studied with this µTFPB and do not aim to shed light on the hydrodynamics and physical chemistry involved in the drainage of these systems.

3.1. Low viscosity foam films (0.24 mM SDS solution) 3). We used the refractive index of water, i.e. nl = 1.33, to obtain the film thickness. One notices a rapid film thinning down to an equilibrium thickness of 24  2 nm. One also observes the well-known stratification phenomenon, with the strata being marked by the grey dotted lines in Figure 4.I.e, typical of surfactant films above the cmc. Such stratification originates from the threedimensional packing of the electrically charged micelles within the film. The micelles leave the film layer by layer, yielding regions of discrete thicknesses depending on the number of micellar layers in the film 15,16,48 . The thickness of each step corresponds to the characteristic distance between micelles fixed by its electrostatic double layer, which here corresponds to 10 ± 3 nm. The equilibrium thickness of ca. 24 nm is reached after 180 s of drainage, which is in agreement with the literature 16,49 . This experiment thus shows that the µTFPB is indeed as efficient as the previous thin film pressure balances for the study of films of low viscosity. The refractive index used to calculate the film thickness was nl = 1.446, which was given by the manufacturer (Gelest). Similar to films from an aqueous surfactant solution, one observes the temporary coexistence of regions of different, well-defined thicknesses (i.e., regions of different colours and grey levels), which progressively decrease with time. Here the strata do not correspond to micelles, but to the characteristic dimensions of the macromolecules, as explained in a previous work 47 . The thickness of each stratum is equal to 15 ± 2 nm. Note that the equilibrium thickness of 25 nm is reached after 20 min, which is a very slow drainage compared to aqueous soap films (the equilibrium thickness was reached after 3 min for the SDS solution). This slow drainage is attributed to the high viscosity (2.7 Pa s) of the DBP-732 melt. Such highly viscous liquids can be treated without any problem with the µTFPB.

Gelling foam films

We show in this section qualitative examples of the drainage of polymer films (with an initial viscosity of 0.42 Pa s) stabilised by Glucopon 600 CSUP (an alkyl polyglycoside from Cognis, now BASF) undergoing physical cross-linking. The alginate solution contains 1.0 wt % alginate (coined as "low viscosity alginate" by the supplier, Alfa Aesar). Calcium ions (Ca 2+ ) are used as cross-linkers and are brought into the solution via calcium chloride dihydrate (from Sigma Aldrich). However, to control cross-linking, we protect the calcium ions with a chelating agent, namely Ethylene glycol-O,O'-bis(2-aminoethyl)-N,N,N',N'-tetraacetic acid, 97% (EGTA). At a neutral pH, the EGTA binds the Ca 2+ ions, which are not available for cross-linking. Upon acidification, below a pH of ca. 4, the EGTA protonates itself and releases the calcium ions, which become available to cross-link the alginates 50 . The concentration of Ca +2 ions and EGTA was the same for all the alginate films, namely 0.1 M. We use 1 wt % D-(+)-glucono-delta-lactone (GDL, from Sigma Aldrich) to acidify the solution and initiate cross-linking 51 .

Figure 5.I.a-d shows the drainage over time of a non-gelling 1 wt % alginate/0.1 M Ca 2+ /EGTA film stabilised with 0.5 wt % Glucopon 600 CSUP. The film does not gel, as no GDL was added to acidify the solution. One sees that at early times, i.e. after 4 min (Figure 5.I.a), the film shows concentric bright and colourful rings centred in the middle of the film, the colour sequence indicating that the film is thickest in its centre. This is a well-known phenomenon, coined as "dimple" 52 . After 7.2 min, one sees in Figure 5.I.b that darker spots of "common black films" 15,[START_REF]Foam and Foam Films, Theory, Experiment, Application[END_REF] appear at the edge of the film, which correspond to much thinner areas. As the film keeps draining, the darker regions percolate and the bright and colourful area on the film declines, as observed in Figure 5.I.c at 8.5 min. An even darker region appears upon further drainage, as seen on the right of the film in Figure 5.I.d, at 11.3 min, corresponding to the formation of a "Newton black film" (hydrated surfactant bilayer) accompanied by the formation of small bright spots known as Rayleigh instabilities 53 . This instability is now reasonably well understood and results from the fact that the solution is expelled from this zone faster than it can be drained away by the remaining film. This leads to the formation of a bump around the film, which, beyond a critical height, breaks into smaller "drops", where the film is locally thick [24][25][26]53 . Here we added 1 wt % of GDL to initiate cross-linking before film formation. Figure 5.II.a shows once again a drainage from the border of the film, leaving a dimple in its centre. However, after 14 min (Figure 5.II.b), we can observe a very different behaviour compared to the non-gelling film. The thickness of the overall film strops decreasing and the central area is not fully homogeneous but contains some zones which are significantly thicker than the surrounding film, as can be inferred from the colours. The precise origin of this phenomenon is not clear to us yet, but we suggest that it may result from a phase separation between polymer-rich globules and a polymer-poor film occurring as a result of confinement/drainage within the film. A more extensive and quantitative study is required to test this hypothesis. After 40 min (Figure 5.II.c) and 56 min (Figure 5.II.d), the area of the flat film keeps increasing, whereas the area of the protruding globules keeps decreasing. The white arrow in Figure 5.III.e shows that the edge of the flat film is not only very well defined but also bounded by an even thinner zone. Such a film profile is hardly imaginable for a purely liquid film and must originate from the gelation of the film. Moreover, we could observe that the film was indeed gelled as we opened the chip to clean it. This is, to the best of our knowledge, the first example of the observation of drainage of a gelling film with the help of a thin film pressure balance. This experiment was made possible by the microfluidic chip which allows for a high entry pressures and an easy cleaning. Further qualitative studies are underway to understand the thinning behaviour observed for the gelling film, but a detailed discussion of the physical and physical-chemical processes at play during the drainage of gelling films is beyond the scope of the paper at hand.

Conclusions

Using a set of illustrative examples, we have introduced a new microfluidic Thin Film Pressure Balance (µTFPB) which allows to investigate the drainage and stability of free-standing liquid films in a reliable manner via the application of a well-defined capillary pressure. The main advantage of this design in comparison to previously used TFPBs is that it is suitable for a wide range of pressures and liquids, including highly viscoelastic or solidifying liquids. It therefore provides a convenient tool to investigate the evolution of solidifying films, for example, in order to understand the influence of interfaces/confinement on gelation and to explore the role of competing processes, such as drainage and gelation. The device can be fabricated using a micromilling machine and can be easily interfaced with the adequate flow/pressure control and optical techniques for the film analysis to fit the system studied and accuracy aimed for. The device design being very flexible, it can be adapted to a multitude of complementary investigations to study, for example, gas exchange or liquid/liquid films.

While the main focus of this article is on the demonstration of the device, the presented preliminary investigation of gelling free-standing films is the first of this kind. It is expected to open the pathway to a new field dedicated to study of gelation under soft confinement and help shed light on the mechanisms responsible for pore opening in polymer foams.

Figure 1 :

 1 Figure 1: (a) Sketch of a free-standing film of thickness h stabilised by surface active agents (not drawn to scale) at a capillary pressure pc = pg -pl. The grey arrows show the liquid flow out of the film. (b) The different contributions of the disjoining pressure  (without bulk contributions), redrawn from 15 . (c) The three main approaches of holding and feeding the film: (I) Scheludko cell, (II) Radke's bike wheel, and (III) Exerowa-Scheludko cell (porous disc).

Figure 3 :

 3 Figure 3: Schematic representation of the microfluidic chip a) divided in its four components, with the green arrow showing the position of the bike wheel on part B and b,c) its assembly from two different orientations, showing where the gas (blue arrows) and the liquid (red arrows) flow. b) and c) are rotated by 90°. Picture of the assembled microfluidic chip: d) top view, showing the gas inlet; e) bottom view, showing the liquid inlet 1 . f) Photograph of the bike wheel, with the red arrows showing the motion of liquid into the wheel to form the film.

Figure 4 .

 4 Figure 4.I.a-d shows the drainage over time of a 0.24 mM SDS solution (viscosity 10 -3 Pa s) subjected to a capillary pressure of 100 Pa monitored with the µTFPB. The different colours and grey levels correspond to different film thicknesses.Figure 4.I.a-d shows the evolution of the film thickness h over time in the region depicted by the yellow circle pointed by the white arrows in the corresponding images using the Scheludko equation (Section 2.3). We used the refractive index of water, i.e. nl = 1.33, to obtain the film thickness. One notices a rapid film thinning down to an equilibrium thickness of 24  2 nm. One also observes the well-known stratification phenomenon, with the strata being marked by the grey dotted lines in Figure4.I.e, typical of surfactant films above the cmc. Such stratification originates from the threedimensional packing of the electrically charged micelles within the film. The micelles leave the film layer by layer, yielding regions of discrete thicknesses depending on the number of micellar layers in the film15,16,48 . The thickness of each step corresponds to the characteristic distance between micelles fixed by its electrostatic double layer, which here corresponds to 10 ± 3 nm. The equilibrium thickness of ca. 24 nm is reached after 180 s of drainage, which is in agreement with the literature16,49 . This experiment thus shows that the µTFPB is indeed as efficient as the previous thin film pressure balances for the study of films of low viscosity.

Figure 4 .

 4 Figure 4.I.a-d shows the drainage over time of a 0.24 mM SDS solution (viscosity 10 -3 Pa s) subjected to a capillary pressure of 100 Pa monitored with the µTFPB. The different colours and grey levels correspond to different film thicknesses.Figure 4.I.a-d shows the evolution of the film thickness h over time in the region depicted by the yellow circle pointed by the white arrows in the corresponding images using the Scheludko equation (Section 2.3). We used the refractive index of water, i.e. nl = 1.33, to obtain the film thickness. One notices a rapid film thinning down to an equilibrium thickness of 24  2 nm. One also observes the well-known stratification phenomenon, with the strata being marked by the grey dotted lines in Figure4.I.e, typical of surfactant films above the cmc. Such stratification originates from the threedimensional packing of the electrically charged micelles within the film. The micelles leave the film layer by layer, yielding regions of discrete thicknesses depending on the number of micellar layers in the film15,16,48 . The thickness of each step corresponds to the characteristic distance between micelles fixed by its electrostatic double layer, which here corresponds to 10 ± 3 nm. The equilibrium thickness of ca. 24 nm is reached after 180 s of drainage, which is in agreement with the literature16,49 . This experiment thus shows that the µTFPB is indeed as efficient as the previous thin film pressure balances for the study of films of low viscosity.

Figure 4 :

 4 Figure 4: Photographs and film thickness analysis of films thinning over time observed using the µTFPB using two different solutions with very different viscosities (10 -3 vs 2.7 Pa s). I.a-d Photographs of an SDS film (at 0.24 mM) thinning under a capillary pressure of 100 Pa; and the corresponding film thickness over time (I.e) for the region displayed by the yellow rings and shown by the white arrows. II.a-d Pictures of a DBP-732 film thinning under a capillary pressure of 100 Pa; and the corresponding film thickness over time (II.e) for the region depicted by the yellow rings and shown by the white arrows.

3. 2 .

 2 High viscosity foam films (DBP-732 polymer melt)

Figure 4 .

 4 Figure 4.II.a-d shows the drainage over time of the polymer melt DBP-732 subjected to a capillary pressure of 100 Pa monitored with the µTFPB.Figure4.II.e shows the evolution of the film thickness with time in the region depicted by the yellow circle pointed by the white arrows in the corresponding images. The refractive index used to calculate the film thickness was nl = 1.446, which was given by the manufacturer (Gelest). Similar to films from an aqueous surfactant solution, one observes the temporary coexistence of regions of different, well-defined thicknesses (i.e., regions of different colours and grey levels), which progressively decrease with time. Here the strata do not correspond to micelles, but to the characteristic dimensions of the macromolecules, as explained in a previous work47 . The thickness of each stratum is equal to 15 ± 2 nm. Note that the equilibrium thickness of 25 nm is reached after 20 min, which is a very slow drainage compared to aqueous soap films (the equilibrium thickness was reached after 3 min for the

Figure 4 .

 4 Figure 4.II.a-d shows the drainage over time of the polymer melt DBP-732 subjected to a capillary pressure of 100 Pa monitored with the µTFPB.Figure4.II.e shows the evolution of the film thickness with time in the region depicted by the yellow circle pointed by the white arrows in the corresponding images. The refractive index used to calculate the film thickness was nl = 1.446, which was given by the manufacturer (Gelest). Similar to films from an aqueous surfactant solution, one observes the temporary coexistence of regions of different, well-defined thicknesses (i.e., regions of different colours and grey levels), which progressively decrease with time. Here the strata do not correspond to micelles, but to the characteristic dimensions of the macromolecules, as explained in a previous work47 . The thickness of each stratum is equal to 15 ± 2 nm. Note that the equilibrium thickness of 25 nm is reached after 20 min, which is a very slow drainage compared to aqueous soap films (the equilibrium thickness was reached after 3 min for the

Figure 5 :

 5 Figure 5: I.a-d) Photographs of a draining 1.0 wt % alginate film stabilised with 0.5 wt % Glucopon 600 CSUP without addition of an acidifying agent, i.e. the film does not undergo gelation and its bulk viscosity remains constant (0.42 Pa s) throughout drainage. II.a-d) Pictures of a draining 1 wt % alginate film stabilised with 0.5 wt % Glucopon 600 CSUP with addition of 1.0 wt % GDL to induce acidification which triggers cross-linking, i.e. gelation. The applied capillary pressure was 300 Pa.

Figure 5 .

 5 Figure5.II.a-d shows the drainage of a similar (yet gelling) 1 wt % alginate film with 0.1 M Ca 2+ /EGTA stabilised with 0.5 wt % Glucopon 600 CSUP. Here we added 1 wt % of GDL to initiate cross-linking before film formation. Figure5.II.a shows once again a drainage from the border of the film, leaving a dimple in its centre. However, after 14 min (Figure5.II.b), we can observe a very different behaviour compared to the non-gelling film. The thickness of the overall film strops decreasing and the central area is not fully homogeneous but contains some zones which are significantly thicker than the surrounding film, as can be inferred from the colours. The precise origin of this phenomenon is not clear to us yet, but we suggest that it may result from a phase separation between polymer-rich globules and a polymer-poor film occurring as a result of confinement/drainage within the film. A more extensive and quantitative study is required to test this hypothesis. After 40 min (Figure5.II.c) and 56 min (Figure5.II.d), the area of the flat film keeps increasing, whereas the area of the protruding globules keeps decreasing. The white arrow in Figure5.III.e shows that the edge of the flat film is not only very well defined but also bounded by an even thinner zone. Such a film profile is hardly imaginable for a purely liquid film and must originate from the gelation of the film. Moreover, we could observe that the film was indeed gelled as we opened the chip to clean it. This is, to the best of our knowledge, the first example of the observation of drainage of a gelling film with the help of a thin film pressure balance. This experiment was made possible by the microfluidic chip which allows for a high entry pressures and an easy cleaning. Further qualitative studies are underway to understand the thinning behaviour observed for the gelling film, but a detailed discussion of the physical and physical-chemical processes at play during the drainage of gelling films is beyond the scope of the paper at hand.

Some extra inlets which are sealed are visible in Figure3d,e). These are remnants of former chip designs, but we left the holes for potential alternative inlets: the chip being close to the objective, some connectors may need to be moved to avoid overload issues.

Acknowledgment

We thank Cosima Stubenrauch, Patrick Kekicheff, Tetiana Orlova and Thibaut Gaillard for fruitful discussions. We also thank Aurélie Hourlier-Fargette for precious advice on alignment issues. We acknowledge funding from the Institute Carnot MICA (ThiFilAn project). We also acknowledge funding from the IdEx Unistra framework (Chair W. Drenckhan), from the European Research Council (ERC-METAFOAM 819511) and from the Institute of Advanced Studies at the University of Strasbourg (USIAS, Fellowship C. Stubenrauch).

Authors' contributions

Conflict of interest

The authors have no conflict of interest to declare.


Andreas Enge
email: enge@lix.polytechnique.fr

COMPUTING MODULAR POLYNOMIALS IN QUASI-LINEAR TIME

Keywords:

We analyse and compare the complexity of several algorithms for computing modular polynomials. Under the assumption that rounding errors do not influence the correctness of the result, which appears to be satisfied in practice, we show that an algorithm relying on floating point evaluation of modular functions and on interpolation has a complexity that is up to logarithmic factors linear in the size of the computed polynomials. In particular, it obtains the classical modular polynomial Φ ℓ of prime level ℓ in time

where M (ℓ) is the time needed to multiply two ℓ-bit numbers.

Besides treating modular polynomials for Γ 0 (ℓ), which are an important ingredient in many algorithms dealing with isogenies of elliptic curves, the algorithm is easily adapted to more general situations. Composite levels are handled just as easily as prime levels, as well as polynomials between a modular function and its transform of prime level, such as the Schläfli polynomials and their generalisations.

Our distributed implementation of the algorithm confirms the theoretical analysis by computing modular equations of record level around 10000 in less than two weeks on ten processors.

Definitions and main result

Modular polynomials, in their broadest sense, are bivariate polynomials with a pair of modular functions as zero. Given any two modular functions f and g for arbitrary congruence subgroups, the function fields C(f ) and C(g) are finite extensions of C(j), so that there are two polynomials relating f resp. g to j. Taking the resultant of these polynomials with respect to j, one sees that a polynomial relationship between f and g exists. In practice, one is rather interested in the minimal polynomial of f over C(g), say, that will be called the modular polynomial of f with respect to g. If the functions satisfy some conditions on the rationality and integrality of their q-expansion coefficients, the modular polynomial has rational integral coefficients.

Different modular polynomials parameterise moduli spaces related to elliptic curves. Let Γ = Sl 2 (Z)/{±1} = PSl 2 (Z) be the full modular group, and C Γ = C(j) the field of modular functions invariant under Γ; j itself parameterises isomorphism classes of elliptic curves. Of special interest for applications is the congruence subgroup

Γ 0 (ℓ) = 1 0 0 ℓ -1 Γ 1 0 0 ℓ ∩ Γ = a b c d ∈ Γ : ℓ|b
for ℓ prime; its modular polynomial parameterises isomorphism classes of elliptic curves together with an isogeny of degree ℓ. This is heavily used in Atkin's and Elkies's improvements to Schoof's algorithm for counting points on elliptic curves [START_REF] Elkies | Elliptic and modular curves over finite fields and related computational issues[END_REF]. Other applications include the computation of the endomorphism ring of an elliptic curve [START_REF] Kohel | Endomorphism rings of elliptic curves over finite fields[END_REF][START_REF] Fouquet | Isogeny volcanoes and the SEA algorithm[END_REF] or of an isogeny between two elliptic curves [START_REF] Steven | Constructing isogenies between elliptic curves over finite fields[END_REF].

More general modular polynomials occur in modern complex multiplication constructions for elliptic curves [START_REF] Enge | Constructing elliptic curves over finite fields using double eta-quotients[END_REF][START_REF]Further investigations of the generalised Weber functions[END_REF][START_REF] Enge | Comparing invariants for class fields of imaginary quadratic fields[END_REF]; they are discussed in Section 4. When counting points on an elliptic curve by the Schoof-Elkies-Atkin algorithm, the modular polynomials are usually supposed to be available after a precomputation phase. This makes sense in situations where the algorithm is carried out many times for curves of about the same size, as is typical in cryptography. When establishing point counting records [START_REF]SEA in genus 1: 2500 decimal digits[END_REF], however, one quickly realises that the computation of the modular polynomials itself becomes a bottleneck that may even dominate from a complexity theoretic point of view.

In this article, after surveying in Section 2 the state of the art as it appears in the literature, I present an algorithm based on evaluation and interpolation that computes modular polynomials in time essentially linear in the bit size of the polynomial. As it is based on floating point computations, one has to assume that rounding errors do not disturb the output, a heuristic that is supported by the implementation described in Section 5. A precise analysis of rounding errors to obtain a rigorously proved result appears to be out of reach. Note, however, that the correctness of the output may be checked probabilistically. For instance, one may instantiate the variable j occurring in the modular polynomial by the j-invariant of an elliptic curve over a finite field and verify that the specialised polynomial has the expected splitting behaviour.

Under the assumption that it is sufficient to use a floating point precision of O(n) for a polynomial whose largest coefficient has n digits, the following holds:

Result 1.1 (heuristic). Let Γ ′ ⊆ Γ be a congruence subgroup, f a modular function for Γ ′ and Φ(X) ∈ C(j)[X] the characteristic polynomial of f with respect to the function field extension C Γ ′ /C(j). Assume that Φ has coefficients in Z[j], and that the largest integer coefficient occurring in it has n digits. Suppose that a system of representatives of Γ ′ \Γ is known and that f can be evaluated at precision O(n) in time O(log n M (n)), where M (n) ⊆ O n log 1+ε is the time needed to multiply two numbers with n digits. Then there is an algorithm that computes Φ in time

O deg X Φ deg j Φ log 2 max(deg X Φ, deg j Φ) + log n M (n) .
In particular, the classical modular polynomial Φ ℓ such that Φ ℓ (j(z), j(ℓz)) = 0 for ℓ prime is obtained in time

O ℓ 2 log 2 ℓ M (ℓ log ℓ) ⊆ O ℓ 3 log 4+ε ℓ .
If Φ is a dense polynomial and all of its coefficients (or at least a constant fraction of them) are of bit size about n, then this algorithm is linear in the size of the polynomial, except for the logarithmic factors.

The time bound of O(log n M (n)) for evaluating modular functions at precision n is motivated by an algorithm due to Dupont; it relies on Newton iterations on an expression involving the arithmetic-geometric mean and reaches this complexity for a wide class of modular functions, including those built from Dedekind's η function and in particular k, k ′ and j [7, Theorem 5 and Section 7.2]. Alternatively, one may use an algorithm relying on multievaluation of q-expansions as described in [START_REF]The complexity of class polynomial computation via floating point approximations[END_REF]Section 6.3]. Its complexity, worse by a factor of log n, is still sufficient for the running time of Result 1.1 to hold. The bound M (n) ∈ O n log 1+ε may be achieved by the classical Schönhage-Strassen algorithm [START_REF] Schönhage | Schnelle Multiplikation großer Zahlen[END_REF] or the more recent and asymptotically faster algorithm due to Fürer [START_REF] Fürer | Faster integer multiplication[END_REF].

The roots of Φ, that is, the algebraic conjugates of f , are given by the f (M ν z) with Γ ′ M ν running through a system of representatives of the residue classes Γ ′ \Γ (see [START_REF] Deuring | Die Klassenkörper der komplexen Multiplikation[END_REF]), so that

(1.1) Φ(X) = [Γ:Γ ′ ] ν=1 (X -f (M ν z)) = X [Γ:Γ ′ ] + [Γ:Γ ′ ] r=1 c r X [Γ:Γ ′ ]-r ,
where the coefficients c r of Φ are the elementary symmetric functions of the conjugates f (M ν z).

By the q-expansion principle (cf. [6, §3]), a sufficient condition for the c r to be polynomials in Z[j] is that f has rational integral q-coefficients, is holomorphic in

H = {z ∈ C : ℑz > 0} and all conjugates f az+b cz+d with a b c d ∈ Γ have integral algebraic q-coefficients.
2. Approaches to computing modular polynomials for Γ 0 (ℓ)

The case of Γ 0 (ℓ) with ℓ prime, which is the most important one for applications, is also the simplest one from a theoretical point of view: The residue classes Γ 0 (ℓ)\Γ are represented by

T ν = 1 ν 0 1
, ν = 0, . . . , ℓ -1, and S = 0 -1 1 0 .

(The literature often concentrates on Γ 0 (ℓ) = ℓ 0 0 1

-1 Γ ℓ 0 0 1 ∩ Γ. Since Γ 0 (ℓ) = S -1 Γ 0 (ℓ)S
, a modular function f for Γ 0 (ℓ) yields the function f (Sz) for Γ 0 (ℓ), and f (Sz) is actually a conjugate of f . So both result in the same modular polynomial Φ. Or phrased differently: C Γ 0 (ℓ) /C Γ is not normal, and another one of its embeddings is C Γ0(ℓ) /C Γ .) Different functions have been suggested in the literature:

• j(z/ℓ), leading to the so-called "classical" or "traditional" modular polynomials;

• w ℓ (z) 2s = η(z/ℓ) η(z) 2s
for s = 12/ gcd(12, ℓ -1), yielding the so-called "canonical" modular polynomials;

• functions that are moreover invariant under the Fricke-Atkin-Lehner involution z → -ℓ z ; in particular, functions of the form T r (η(z)η(ℓz)) η(z)η(ℓz) evaluated in -1 z with T r the r-th Hecke operator and linear combinations thereof as described in [START_REF] Volker Müller | Ein Algorithmus zur Bestimmung der Punktanzahl elliptischer Kurven über endlichen Körpern der Charakteristik größer drei[END_REF]Ch. 5.3].

All these functions lead to modular polynomials with coefficients in Z.

2.1.

Heights of modular polynomials. The (logarithmic) height of a modular polynomial with coefficients in Z is defined as the logarithm of the largest absolute value of the coefficients. It provides a lower bound on the arithmetic precision required to compute the polynomial. For the classical polynomial between j(z) and j(z/ℓ), the height is given by 6(ℓ + 1) (log ℓ + O(1)) ⊆ O(ℓ log ℓ) according to [START_REF] Cohen | On the coefficients of the transformation polynomials for the elliptic modular function[END_REF]. For alternative modular polynomials, one observes a similar growth of the coefficients with ℓ, but the constant is usually smaller. It should be possible to adapt the argumentation in [START_REF] Cohen | On the coefficients of the transformation polynomials for the elliptic modular function[END_REF] to obtain a similar bound for further classes of modular polynomials.

2.2.

Using q-expansions. The system of representatives of Γ 0 (ℓ)\Γ and thus the conjugates of the modular function f being explicitly known, one may use (1.1) to compute the modular polynomial Φ, if only a procedure for recognising its coefficients c r as polynomials in Z[j] is available. The straightforward and most popular approach is to use q-expansions: The expansion of j is of the form q -1 + • • • , so that knowing the expansion c r = c r,k q -k + • • • , one deduces that c r has leading term c r,k j k as a polynomial in j; one subtracts the q-expansion of this term and continues in the same vein. In fact, only the non-positive q-powers in the series expansions of the c r are needed. Several variants of the algorithm are described in [START_REF] Morain | Calcul du nombre de points sur une courbe elliptique dans un corps fini: aspects algorithmiques[END_REF] with a brief complexity analysis in [8, Section 3]. We shall describe only the algorithm with the best complexity and give a detailed analysis of its running time.

Let f be a modular function for Γ 0 (ℓ) with valuation v < 0 at infinity; that is, f admits a q-expansion of the form

f (z) = ∞ k=vℓ a k q k/ℓ ∈ C((q 1/ℓ )) with q 1/ℓ = e 2πiz/ℓ .
Suppose that the a k are rational integers. (For instance, f = j(z/ℓ) resp. f = w 2s ℓ are valid choices with v = -1/ℓ resp. v = -s(ℓ-1) 12ℓ .) Then the q-expansions of the conjugates f ν = f (T ν z) are of the same form with the coefficients a k being multiplied by powers of ζ ℓ = e 2πi/ℓ .

The last conjugate f ∞ = f (Sz) poses difficulties since it cannot be described generically, q -1 z = q -2πi/z being unrelated to q in a simple way, and is thus treated separately. (And it is the fact that there is only one of them that makes the case of Γ ′ = Γ 0 (ℓ) stand out so that it is comparatively easily handled via the q-expansion approach.) Denote by v ∞ the order of f ∞ at infinity. For

f = j(z/ℓ), one has f ∞ = j(ℓz) = q -ℓ + • • • and v ∞ = -ℓ, and for f = w 2s ℓ , the transformation behaviour of η yields f ∞ = √ ℓ η(ℓz) η(z) 2s = ℓ s q s(ℓ-1)/12 + • • • with v ∞ = s(ℓ-1) 12 = ℓ-1 gcd(12,ℓ-1) .
In particular, v ∞ may be positive or negative. If it is negative, the term with lowest q-exponent occurring in

(1.1) is f ∞ • ℓ-1 ν=0 f ν of exponent ℓv + v ∞ ; so the degree of Φ in j is |ℓv + v ∞ |, reached for c ℓ+1 .
For instance, the degree in j is ℓ + 1 for the classical modular polynomial Φ ℓ . If v ∞ is positive, the term with lowest q-exponent is ℓ-1 ν=0 f ν , so the degree |ℓv| of Φ in j is reached in c ℓ . For instance, the degree in j of the canonical modular polynomials is s(ℓ-1)

12 .

We first consider the complexity of determining the partial product

ℓ-1 ν=0 (X -f ν ) = X ℓ + ℓ r=1 cr X ℓ-r ∈ Z[ζ ℓ ]((q 1/ℓ ))[X].
One may proceed by computing the Newton sums s r = ℓ-1 ν=0 f r ν for 1 ≤ r ≤ ℓ and use Newton's recurrence formulae to deduce the cr . Following [START_REF] Morain | Calcul du nombre de points sur une courbe elliptique dans un corps fini: aspects algorithmiques[END_REF], one obtains with

f r = ∞ k=rvℓ a k,r q k/ℓ that (2.1) s r = ∞ k=⌈rv⌉ ℓ a ℓk,r q k ∈ Z((q)),
so that in fact the roots of unity do not interfere with the computations. One sees that the valuations of the s r at infinity are bounded below by ⌈rv⌉ (and in general, they will be equal to it). Assume that d + 1 terms of them are known, that is, s r is known up to and including the coefficient of q ⌈rv⌉+d . Then one shows by induction on Newton's relation

(2.2) cr = 1 r r k=1 (-1) k+1 s k cr-k
that the valuation of cr is also bounded below by ⌈rv⌉ and that it is also known up to q ⌈rv⌉+d . If v ∞ ≥ 0, we need only the terms of the cr with non-positive exponents; for v ∞ < 0, we loose some precision, in fact |v ∞ | terms, through the final multiplication by X -f ∞ . Thus, the value required for d is given by

d = ℓ|v| + max(0, -v ∞ ) = deg j Φ.
This is an indication that the algorithm given above should be optimal among those relying on q-expansions: Outputting polynomials with up to deg j Φ + 1 coefficients, it manipulates series with as many terms. Unfortunately, the computation of the s r involves decimating the series for f r by ℓ. Eventually, a factor of ℓ is lost in the running time since O(ℓd) terms of the f r are needed.

Once the q-coefficients of f are known, the f r and s r are thus computed by O(ℓ) multiplications of series with O(ℓd) terms and the cr by O(ℓ 2 ) multiplications of series with O(d) terms. Let M q (d) be the number of arithmetic operations in Z required to multiply two dense q-expansions with d terms. As M q is at least linear, the total complexity of the series computations becomes

O(ℓ M q (ℓd)).
For functions f related to η quotients, such as

f = w 2s ℓ or f = j(z/ℓ) = (w 24 2 (z/ℓ)+16) 3 w 24 2 (z/ℓ)
, the effort for computing their q-expansions is negligible, since it corresponds to a constant number of arithmetic operations with series starting from the easily written down series expansion of η.

To write the c r as polynomials in j, one has to compute the non-positive parts of the q-expansions of the j k for 1 ≤ k ≤ d. This can be done in time O(d M q (d)), which is negligible compared to the previous steps because d ∈ O(ℓ). Identifying the c r as polynomials then corresponds to solving triangular systems of linear equations, requiring altogether O ℓd 2 operations with integers.

The total complexity thus becomes

O ℓ M q (ℓ 2 ) integer operations.

Let n ∈ O(ℓ log ℓ) as in Section 2.1 be a bound on the height of the modular polynomial. Then the bit complexity of the algorithm is given by

O ℓ 3 log ℓ M (n) ⊆ O ℓ 4 log 3+ε ℓ
In [3, Appendix A] it is argued that in fact the running time is higher, on grounds that the coefficients of the powers of j grow faster than O(ℓ log ℓ). However, this objection can be dismissed by carrying out all computations modulo a sufficiently large prime of bit size n ∈ O(ℓ log ℓ). Alternatively, and probably preferably in practice, one may work modulo small primes and use the Chinese remainder theorem. Both approaches yield the desired complexity.

Remark. If Φ ℓ is sought only modulo some prime p, then the complete algorithm can be carried out modulo p as soon as p is larger than ℓ to make the divisions in (2.2) possible. The bit complexity becomes O ℓ 3 log ℓ log p (log log p) 1+ε .

2.3. Charles-Lauter. In [START_REF] Charles | Computing modular polynomials[END_REF], the authors describe an algorithm to compute the classical modular polynomial Φ ℓ directly modulo a prime p without recourse to q-expansions of modular functions. Instead, they rely on the moduli interpretation of Φ ℓ , which characterises pairs of ℓ-isogenous elliptic curves. So the algorithm manipulates only elliptic curves and explicit isogenies over extension fields of F p .

Its basic building block is the computation of the instantiated polynomial Φ = Φ ℓ (X, j) with j ∈ F p 2 the j-invariant of a supersingular elliptic curve E. The roots of Φ are the j-invariants of the ℓ + 1 curves that are ℓ-isogenous to E, and that may be obtained via Vlu's formulae [START_REF] Vélu | Isogénies entre courbes elliptiques[END_REF] once the complete ℓ-torsion of E is known. These are ℓ 2 points defined over an extension of F p 2 of degree O(ℓ) in the supersingular case; in fact, one expects a degree of Θ(ℓ) virtually all of the time. Thus writing down the ℓ-torsion points requires a time of Ω(ℓ 3 log p).

After having repeated the procedure for deg j Φ ℓ = Θ(ℓ) different values of j, one obtains the coefficients of Φ ℓ modulo p by interpolation. (This is analogous to the floating point approach of Section 3, and more details are given there.) The complexity of the algorithm is thus at least Ω ℓ 4 log p , and an upper bound of O ℓ 4+ε log 2+ε p + log 4+ε p is proved in [3, Theorem 3.2] under the generalised Riemann hypothesis.

Hence, this algorithm is slower by a factor of order ℓ than the one presented in Section 2.2 relying on q-expansions. This is due to the costly determination of isogenies via their kernels: While the isogenies are defined over F p 2 , their kernels lie in an extension of degree O(ℓ). The q-expansions of the conjugates of j(z/ℓ), on the other hand, provide a synthetic description of the curves that are ℓ-isogenous to the one with j-invariant j(z), and using them it is sufficient to carry out all computations in F p in order to obtain Φ ℓ mod p.

Evaluation-interpolation

The approach for calculating modular polynomials that is described in this section is in fact neither new nor particularly involved. It has been suggested to me by R. Schertz during our work on the class invariants of [START_REF] Enge | Constructing elliptic curves over finite fields using double eta-quotients[END_REF], and later R. Dupont pointed out to me that it can already be found in [START_REF] Borwein | Pi and the AGM[END_REF]Chapter 4.5]. However, it has not been noticed before that the algorithm allows to lower the exponent of the computational complexity by 1 and thus to reach an essentially optimal complexity if fast polynomial arithmetic and fast techniques for evaluating modular functions as described in [START_REF] Dupont | Fast evaluation of modular functions using Newton iterations and the AGM[END_REF][START_REF]The complexity of class polynomial computation via floating point approximations[END_REF] are used.

The basic idea of the evaluation-interpolation approach is to specialise and to compute the identity (1.1) between modular functions in several complex floating point arguments (the evaluation phase); and then to interpolate the coefficients in order to recognise them as polynomials in j. Precisely, (1.1) can be rewritten as

(3.1) Φ(X, j(z)) = [Γ:Γ ′ ] ν=1 (X -f (M ν z)) = X [Γ:Γ ′ ] + [Γ:Γ ′ ] r=1 c r (z)X [Γ:Γ ′ ]-r
for all z in the upper complex half plane.

Evaluating the conjugates f (M ν z) of f in a number of complex arguments z k with ℑz k > 0, multiplying out the left hand side as a polynomial in C[X] and separating the coefficients according to powers of X yields the values c r (z k ). Writing c r = deg j Φ s=0 c r,s j deg j Φ-s , the c r (z k ) are actually the values of the polynomial c r in j(z k ), so that the coefficients c r,s ∈ C can be retrieved by interpolation as soon as the c r (z k ) are known for deg j Φ+1 values of z k . The final step is to round the c r,s to rational integers. The degree of the modular polynomial in j or an upper bound on it may be known beforehand; if it is not, then the algorithm can be repeated with increasing guesses for deg j Φ until rounding to integers succeeds (and doubling the guess at each time results in the same asymptotic complexity as taking the correct value from the beginning).

Let E(n) be the bit complexity for evaluating the modular functions f or j at a precision of n bits (here, f is considered to be fixed, while n tends to infinity with ℓ). Then the evaluation phase requires (ℓ + 1)(

deg j Φ + 1) evaluations of f at a cost of O ℓ deg j Φ E(n)
and the reconstruction of deg j Φ + 1 polynomials of degree ℓ + 1 from their roots. Multiplying complex polynomials by the FFT, this reconstruction step takes

O deg j Φ ℓ log 2 ℓ M (n) + log n M (n) ,
where the (eventually negligible) term log n M (n) accounts for the computation of a primitive root of unity of sufficiently high order to carry out all the FFTs, and as usual M (n) is the time needed to multiply two n-bit numbers.

The interpolation phase consists of ℓ interpolations of polynomials of degree deg j Φ. Employing again fast algorithms such as [START_REF] Von | Modern computer algebra[END_REF]Algorithm 10.11], this takes

O ℓ deg j Φ log 2 (deg j Φ) M (n) ,
once the roots of unity are available.

It is shown in [START_REF] Dupont | Fast evaluation of modular functions using Newton iterations and the AGM[END_REF] that among others, Dedekind's η function can be evaluated at precision n in time

E(n) ∈ O (log n M (n)) uniformly for the argument in F = z ∈ H : |ℜz| ≤ 1 2 , |z| ≥ 1 .
(The restriction to F is not crucial since we may freely choose our interpolation points.) So in particular, all functions built from η such as j and w ℓ satisfy E(n) ∈ O (log n M (n)). (Alternatively, one may use multievaluation as described in [10, Section 6.3] with a slightly worse amortised cost of O log 2 n M (n) per value, which still allows to reach the final complexity of Result 1.1.)

In total, the steps add up to a running time of

O (ℓ deg j Φ log 2 max(ℓ, deg j Φ) + log n)M (n) ⊆ O (deg X Φ deg j Φ log 2 max(deg X Φ, deg j Φ) + log n)M (n) ,
and the latter formulation is valid for arbitrary congruence subgroups Γ ′ in the place of Γ 0 (ℓ).

Here, n must be at least as large as the logarithmic height of Φ, and maybe bigger to account for rounding errors. In the case of the classical modular polynomial Φ ℓ the bound O(ℓ log ℓ) of Section 2.1 yields a complexity of O ℓ 2 log 2 ℓM (ℓ log ℓ) .

This proves Result 1.1.

Further kinds of modular polynomials

4.1. Modular functions of composite level. Besides its better complexity compared to the approach using q-expansions, the evaluation-interpolation algorithm has the advantage of a great flexibility, which makes it easily adaptable to a variety of different modular polynomials. In fact, the proof of Result 1.1 in Section 3 does not use special properties of Γ 0 (ℓ) and is valid for any congruence subgroup.

An application is given by the computation of the polynomial relationship between j and a modular function f of composite level N . At first sight, these polynomials are not of great interest. In the point counting or more generally isogeny computation context, it is more efficient to express an isogeny of composite degree as a composition of prime degree isogenies. But this kind of modular polynomials occurs naturally in the context of [START_REF] Enge | Constructing elliptic curves over finite fields using double eta-quotients[END_REF], in which elliptic curves with complex multiplication are obtained via modular functions f of composite level. Levels N = pq or N = p 2 with p and q prime are examined in [START_REF] Enge | Constructing elliptic curves over finite fields using double eta-quotients[END_REF], but more general settings can be devised in a straightforward way. The polynomial relationship between f and j is used to derive from a special value of f the j-invariant of the corresponding elliptic curve.

If the q-expansion approach were to be pursued to compute this kind of modular polynomials, it would be necessary to somehow obtain the q-expansions of the conjugates f (M ν z) for a system of representatives (M ν ) of Γ ′ \Γ. This is straightforward only for translations; for all other matrices, it is necessary to take the transformation behaviour of f under unimodular matrices into account, which requires ad hoc computations for each particular function. (This is illustrated by j(Sz/ℓ) = j(ℓz) and w 2s ℓ (Sz) = ℓ s η(ℓz) η(z)

2s

, which are not derived from j(z/ℓ) resp. w ℓ (z) by the same generic transformation.) In the case of Γ 0 (ℓ), only the matrix S asks for special treatment; in the case of composite level, many more special matrices appear. For instance, the full article [START_REF]Modular curves of composite level[END_REF] is concerned with deriving the conjugates of only the functions η(z/p1)η(z/p2) η(z/(p1p2))η(z) s for p 1 , p 2 prime and s = 24/ gcd(24, (p 1 -1)(p 2 -1)).

In the evaluation-interpolation approach, the only difference between Γ 0 (ℓ) and other congruence subgroups Γ ′ is the enumeration of the system of representatives (M ν ). For most interesting Γ ′ (and in particular for arbitrary Γ 0 (N )), this is easy; in any case, this step is independent of the particular function f . Then Result 1.1 applies, and one obtains an algorithm that is essentially linear in its output size. 4.2. Schläfli equations. In [START_REF] Schläfli | Beweis der Hermiteschen Verwandlungstafeln für die elliptischen Modulfunktionen[END_REF], Schläfli examines transformations of prime level ℓ of special modular functions different from j that lead to particularly simple modular equations. Weber gives a systematic treatment of them in [START_REF] Weber | Elliptische Funktionen und algebraische Zahlen[END_REF]

, § §73-74]. Let f = ζ -1 48 η((z+1)/2) η(z)
= ζ 48 w 2 (z + 1) be one of "the Weber functions", a modular function of level 48. Let ℓ be a prime not dividing this level. Then g(z) = f(z/ℓ) is the root of a monic polynomial Φ f ℓ (X) of degree ℓ + 1 with coefficients in Z[f]. The evaluation-interpolation approach allows to easily obtain this polynomial with only minimal modifications to the algorithm. The function f being modular for a subgroup of Γ(48) and ℓ being different from 2 and 3, a quick computation reveals that g is modular for Γ(48) ∩ Γ 0 (ℓ). The polynomial Φ f ℓ is thus given by an equation analogous to (1.1):

(4.1) Φ f ℓ (X) = ℓ+1 ν=1 (X -g(M ν z)) = X ℓ+1 + ℓ r=1 c r X ℓ+1-r ,
where the M ν range over a set of representatives of (Γ(48) ∩ Γ 0 (ℓ))\Γ(48). Such a set may be obtained by multiplying the standard representatives of Γ 0 (ℓ)\Γ from the left by a matrix in Γ 0 (ℓ) such that they end up in Γ(48). Precisely, a possible set of representatives is given by 1 48ν 0 1 for ν = 0, . . . , ℓ -1 (corresponding to the translations) and 1 -48k 48k -48k 1 + 48k with k = 48 -1 mod ℓ, corresponding to the inversion S. So the only modification required for the evaluation-interpolation algorithm to work is this adaptation of the matrices M ν , and of course j has to be replaced by f in the interpolation phase to recover the coefficients c r as elements of Z[f].

Hence Result 1.1 clearly applies also to the Schläfli equations, showing that they can be computed in essentially linear time with respect to the output size.

Taking specifics of the function f into account, a more efficient algorithm can be obtained, gaining a constant factor. Weber shows in [32, p. 266] that only every 24-th coefficient in f of the c r may be non-zero. Precisely, f i g k having a non-zero coefficient implies ℓi + k ≡ ℓ + 1 (mod 24). Hence the number of evaluations can be reduced by a factor of about 24. (In [START_REF] Borwein | Pi and the AGM[END_REF]Ch. 4.5], similar sparse modular equations are studied, in which one out of eight coefficients is non-zero. The Borwein's suggest in this chapter the evaluation-interpolation approach while profiting of this sparseness.) Weber shows that the Φ f ℓ are symmetric, which could also be taken into account during the interpolation phase.

Weber derives the Schläfli equations in a very compact form. He considers the new functions

A = f g r + g f r and B = (fg) s + 2 ℓ
2 s (fg) s with r, s such that 12|(ℓ -1)r, 12|(ℓ + 1)s and (ℓ-1)r 12 ≡ (ℓ+1)s 12 (mod 2). Depending on ℓ mod 24, these conditions are satisfied by some r and s such that 2r|ℓ + 1 and 2s|ℓ -1. Weber [32, p. 268] then shows that

Φ f ℓ = (fg) (ℓ+1)/2   A (ℓ+1)/(2r) -B (ℓ-1)/(2s) + (ℓ+1)/(2r)-1 α=0 (ℓ-1)/(2s)-1 β=0 c α,β A α B β   .
In this form, the polynomial has about ℓ 2 /(4rs) coefficients, which is often less than the roughly ℓ 2 /24 coefficients if it is written as an equation between f and g. However, the more compact form appears to be less suited to the evaluationinterpolation algorithm: There is no easy way, in the spirit of (1.1), of obtaining its values when interpreting it as a univariate polynomial in A(z) with coefficients in Z[B(z)]. So instead of performing O(ℓ) interpolations of polynomials of degree O(ℓ), one would apparently have to obtain all coefficients at once by solving a linear system with O(ℓ 2 ) unknowns, resulting in a complexity that is worse than that of Result 1.1.

4.3.

Generalised Schläfli equations. The Schläfli equations of the previous paragraph can obviously be generalised to further modular functions: Given a modular function f and a prime ℓ, one may want to compute the polynomial relating f (z) and f (z/ℓ).

In [START_REF] Hart | Evaluation of the Dedekind eta funktion[END_REF]Ch. 5], a theory analogous to Weber's is developped for the functions w ℓ . The derived polynomials are used to obtain symbolically minimal polynomials of special algebraic values of η quotients.

In the context of a p-adic algorithm for the computation of class polynomials and ultimately elliptic curves with complex multiplications, another application is developed in [START_REF] Bröker | Constructing elliptic curves of prescribed order[END_REF]Ch. 6.8]. Let τ be a quadratic integer. Then the singular value j(τ ) lies in the ring class field for the order O = [1, τ ] Z and is the j-invariant of an elliptic curve with complex multiplication by O. In [START_REF] Couveignes | Action of modular correspondences around CM points, Algorithmic Number Theory -ANTS-V[END_REF], the authors describe an algorithm to compute the canonical lift of j(τ ) to the p-adic numbers at arbitrary precision by some kind of Newton iteration. The main algorithmic step consists of applying an isogeny to j(τ ) (or more precisely to the elliptic curve with this j-invariant) that corresponds to a smooth principal ideal and that can thus be realised by composing isogenies of manageable prime degrees. In [START_REF] Bröker | Constructing elliptic curves of prescribed order[END_REF], the algorithm is generalised to other class invariants. Let f be a modular function for Γ 0 (N ) such that f (τ ) is a class invariant in the sense that it also lies in the ring class field of O. Let ℓ be a prime number not dividing N that splits as (ℓ) = ll in O. As in the case of j(τ ), one needs to "apply an isogeny" and compute f (l -1 ). This value is given as a root of the modular polynomial Ψ between f (z) and f (z/ℓ), specialised with f (τ ) in the place of f (z). (The modular polynomial Φ between f (z) and j(z), treated in Section 4.1, also plays a role as it helps to choose whenever Ψ has multiple roots: f (l -1 ) is also a root of Φ specialised with j(l -1 ) in the place of j(z).)

A quick computation shows that g(z) = f (z/ℓ) is modular for Γ 0 (ℓ) ∩ Γ 0 (N ) = Γ 0 (ℓN ), so that (4.1) yields the minimal polynomial of g over C Γ 0 (N ) if the M ν are chosen as a system of representatives of Γ 0 (ℓN )\Γ 0 (N ). Such a system is obtained precisely as for the Schläfli equations by multiplying the standard representatives of Γ 0 (ℓ)\Γ by matrices in Γ 0 (ℓ) so that they end up in Γ 0 (N ). For instance, a set of representatives is given by 1 N ν 0 1 for ν = 0, . . . , ℓ -1 and

1 -N k N k -N k 1 + N k with k = N -1 mod ℓ.
So at first sight, the evaluation-interpolation algorithm seems to apply.

The problem lies in recognising the coefficients c r as polynomials in f . In fact, the c r are modular for Γ 0 (N ), so unless C Γ 0 (N ) = C(f ) (otherwise said, f is a Hauptmodul for Γ 0 (N )), there is no reason that they are rational in f . Generically, one expects C Γ 0 (N ) = C(j, f ) and may hope (under suitable integrality and rationality conditions) to recover the c r as elements of Z[j, f ].

A necessary condition for f being a Hauptmodul is that the modular curve X 0 (N ) is of genus 0, which limits the possible values of N to a very short list. In general, it will be necessary to replace the minimal polynomial of g over C Γ 0 (N ) by that over C(f ); its degree will be ℓ + 1 multiplied by [C Γ 0 (N ) : C(f )]. It is not clear whether the evaluation-interpolation approach can be adapted to this context, since C(f ) is not the field of modular functions for any congruence subgroup.

Implementation

The evaluation-interpolation algorithm has been implemented in C by the author for different kinds of modular polynomials. One big advantage of the algorithm (besides its superior complexity) is that it is rather straightforward to implement, especially in the context of complex multiplication (cf. [START_REF]The complexity of class polynomial computation via floating point approximations[END_REF]), where the major building blocks such as the evaluation of modular functions and computations with polynomials over floating point numbers are already in use. The author's implementation relies on the existing libraries gmp [START_REF] Granlund | gmp -GNU multiprecision library[END_REF] with an assembly patch for AMD64 by P. Gaudry [START_REF] Gaudry | Assembly support for GMP on AMD64[END_REF] for fast basic multiprecision arithmetic, mpfr [START_REF] Hanrot | mpfr -a library for multiple-precision floating-point computations with exact rounding[END_REF] and mpc [START_REF] Enge | mpc -a library for multiprecision complex arithmetic with exact rounding[END_REF], which provide elementary operations with arbitrary precision floating point real and complex numbers with exact rounding. A library for fast polynomial operations via Karatsuba, Toom-Cook and the FFT has been written for the occasion [START_REF] Enge | mpfrcx -a library for univariate polynomials over arbitrary precision real or complex numbers[END_REF]. All timings mentioned in the following have been obtained on Opteron processors clocked at 2.4 GHz. The heights are given as logarithms in base 2, and the arithmetic precision is given in bits.

Details of the implementation.

5.1.1. Details slowing the implementation down asymptotically. Of the different techniques for evaluating modular functions described in [START_REF]The complexity of class polynomial computation via floating point approximations[END_REF], the asymptotically fast ones are not beneficial in the present context. For computing univariate class polynomials, they do not pay off at degrees below 100 000. Such high degrees are for the time unachievable for bivariate modular polynomials, since the number of coefficients would make it impossible to store the result. Hence all examples have been computed via the sparse representation of the η function as a q-series.

The interpolation phase turns out to take negligible time; thus a simple quadratic algorithm using iterated differences is employed instead of a quasilinear one. 5.1.2. Details speeding the implementation up. In order to work with real instead of complex numbers as much as possible, purely imaginary values are chosen for the z k . Then the q(z k ) are real, and so are the values j(z k ) and thus, by (3.1), the Φ(X, j(z k )). While the values of the conjugates f (M ν z k ) still have to be computed as complex numbers, roughly half of them suffice for functions f for Γ 0 (ℓ) with rational q-coefficients: f (z k ) is real, f (z k -r) is the complex conjugate of f (z k + r) for r ∈ Z, and finally f (-1/z k ) has to be real again to obtain a real polynomial in the end. When reconstructing Φ(X, j(z k )) from its roots, the complex conjugate roots can be grouped so that only real arithmetic is required.

Also the interpolation phase uses only real numbers. It is furthermore helped by choosing the z k such that the j(z k ) lie in a simple arithmetic progression, for instance, j(z k ) = 1728 + k for k ≥ 1.

5.1.3. Parallelisation. Another nice feature of the evaluation-interpolation algorithm is that it can be parallelised almost arbitrarily, in principal up to distributing at the level of each modular function evaluation. A more natural, coarser parallelisation of the evaluation phase, which needs almost no communication, is amply sufficient in practice: Each processor treats a different z k and computes an approximation Φ(X, j(z k )) by evaluating the conjugates (f (M ν z k )) ν and reconstructing the polynomial from these roots. The result is written into a file on a shared hard disk. The interpolation phase has been too fast to warrant a parallel implementation. If it were to become a bottleneck, a natural way of distributing the work would be to have each processor compute the coefficients of a different polynomial c r ∈ Z[j].

Implementing the communication via files solves a second problem. During the evaluation phase, the matrix (c r (z k )) k,r is computed row-wise; for the interpolation, it is accessed column-wise. However, the matrix is roughly as big as the final modular polynomial and for the largest computed examples (see Section 5.2.1) does not fit into main memory any more. Writing each row into a different file during the evaluation and reading all files in parallel during the interpolation can be seen as a means of transposing the matrix on disk without having to keep it in main memory. 5.1.4. Arithmetic precision. Except for the classical modular polynomials, no upper bound on the height of the polynomials is known, and even in the classical case, the bound is not completely explicit, but contains an O(1) term (see Section 2.1). In the implementation, an approximate bound is obtained by carrying out all computations first at very low precision (100 bits). Due to the numeric instability of interpolation, the resulting coefficients are wrong already in the first digit, but one obtains an idea of the height of the correct polynomial (actually, the height at low precision turns out to be larger than the real one). This estimate is multiplied by a factor (between 1.1 and 2 depending on the function and the level), determined experimentally, to deal with rounding errors.

In the context of class polynomial computations, it has been observed in [START_REF]The complexity of class polynomial computation via floating point approximations[END_REF] that virtually no rounding errors occur: Evaluating modular functions via the qdevelopment of η as well as multiplication of polynomials is rather uncritical from a numerical point of view. It suffices to add a few bits to the target precision. This should also hold for the evaluation phase in our context. In practice, it can be observed, however, that a significantly higher precision is needed than just the height of the result. This can only be due to the interpolation phase, which implicitly handles Vandermonde matrices, that are known to be ill conditioned. It is an open problem to choose the interpolation points j(z k ) so as to minimise the rounding errors.

Examples.

5.2.1. Modular polynomials for Γ 0 (ℓ). For achieving the point counting record for elliptic curves with the Schoof-Elkies-Atkin algorithm [START_REF]SEA in genus 1: 2500 decimal digits[END_REF], the polynomials have been systematically computed with the function

(5.1) f ℓ,r = T r (η(z)η(ℓz)) η(z)η(ℓz)
(evaluated in -1/z) suggested in [START_REF] Volker Müller | Ein Algorithmus zur Bestimmung der Punktanzahl elliptischer Kurven über endlichen Körpern der Charakteristik größer drei[END_REF]. Here,

T r (f ) = 1 r r-1 ν=0 f z + 24ν r + f (rz)
is the r-th Hecke operator for a prime r ≥ 5 such that 24|(r -1)(ℓ + 1), r is a quadratic residue modulo ℓ and ℓ is a quadratic residue modulo r. (The unusual factor 24 in front of ν takes care of the fact that the exponents in the q-expansion of η are not integral due to the factor q 1/24 ; it eliminates 24-th roots of unity.) It is not advisable to obtain the values of f ℓ,r via its q-expansion. First of all, the expansion is dense and thus much slower to evaluate than that of η. But more importantly, q tends to infinity when z approaches the cusp 0 of the fundamental domain for Γ 0 (ℓ), so that the q-series converges arbitrarily slowly. This is not an issue when relying on the evaluation of η only, since its known transformation behaviour under unimodular matrices allows to transform the arguments into the fundamental domain for Γ, whose only cusp is at infinity and corresponds to q = 0. So the implementation evaluates (5.1) numerically with 2r + 4 evaluations of η per value of f ℓ,r . The resulting dependence of the running time on r can be seen clearly in the following examples. The largest computed polynomial has a level of 10079. The computation takes less than two weeks on a small cluster with ten processors; as a compressed text file, the polynomial fills about 16 Gb of space. This illustrates that the limiting factor for the algorithm becomes not so much its running time, but rather the space requirements of its output, as can be expected for algorithms that have an essentially linear time complexity with respect to their output size. ℓ 10079 r 5 degree in j 672 estimated height 31865 precision 35051 height 28825 time for evaluation 10 000 000 s ≈ 120 d time for interpolation 56 000 s ≈ 16 h It would be natural to try to linearly combine functions f ℓ,r with different r to reduce the pole order at infinity, or even to find a function on X 0 (ℓ) with minimal pole order. As a result, the degree of the modular polynomial in j (and thus the number of interpolation points) and the size of its coefficients (and thus the required precision) could be reduced. A high price, however, would have to be paid by an increased running time for the evaluation. As high performance clusters are becoming increasingly available, while bandwidth and storage space appear to be the bottleneck for modular polynomial computation, this might, however, be the road to follow for constructing a manageable database of modular polynomials up to a high level. 5.2.2. Schläfli equations. As explained in Section 4.2, a trivial modification to the enumeration of the ℓ + 1 matrices M ν adapts the code to computing the Schläfli equations between f(z) and f(z/ℓ). It is then imperative to change the interpolation code so as to take advantage of the sparseness of the polynomial, in which only every 24-th coefficient is non-zero. This reduces the number of evaluations and (after a suitable transformation) the degree of the interpolated polynomials roughly by a factor of 24. The latter is also important since the interpolation phase becomes numerically more stable. 3, the minimal polynomial of g(z) = f (z/ℓ) over C(f ) for some function f for Γ 0 (N ) need not be of degree ℓ + 1 any more, in which case it is a priori unclear how to obtain it by evaluation-interpolation. But besides the functions f that are Hauptmoduln for Γ 0 (N ), a few others may be handled by this approach. Namely, let w p1,p2 = η(z/p 1 )η(z/p 2 ) η(z)η(z/(p 1 p 2 ))

with p 1 , p 2 prime and 24|(p 1 -1)(p 2 -1) be the functions of level N = p 1 p 2 suggested as class invariants in [START_REF] Enge | Constructing elliptic curves over finite fields using double eta-quotients[END_REF]. These functions are invariant under the Fricke-Atkin-Lehner involution. Now, X + 0 (N ) is of genus 0 for N = 35 and N = 39 [START_REF] Ogg | Hyperelliptic modular curves[END_REF], and apparently, the double η quotients are Hauptmoduln in this case. The degree of the modular equation of transformation level ℓ thus remains ℓ + 1 (an observation made for N = 35 in [START_REF] Bröker | Constructing elliptic curves of prescribed order[END_REF]Ch. 7.4]).

For instance, with f = w 3,13 and g(z) = f (z/ℓ), one obtains the following polynomials Φ w3,13 ℓ : Φ w3,13 2 = g 3 + f 3 -g 2 f 2 -gf + 2(g 2 f + f g 2 ) Φ w3, 13 5 = g 6 + f 6 -g 5 f 5 -gf + 35g 3 f 3 +5(g 5 f 4 + g 4 f 5 + g 5 f + gf 5 -g 4 f 3 -g 3 f 4 -g 3 f 2 g 2 f 3 + g 2 f + gf 2 ) +10(-g 5 f 2 -g 2 f 5 + g 4 f 4 + g 4 f 2 + g 2 f 4 -g 4 f -gf 4 + g 2 f 2 ) Φ w3, 13 7 = g 8 + f 8 -g 7 f 7 -gf +7(g 7 f 6 + g 6 f 7 -g 7 f 5 -g 5 f 7 + g 6 f 4 + g 4 f 6 + g 6 f 2 + g 2 f 6 ) +7(g 4 f 2 + g 2 f 4 -g 3 f -gf 3 + g 2 f + gf 2 ) +14(-g 7 f -gf 7 + g 5 f 4 + g 4 f 5 + g 5 f 3 + g 3 f 5 + g 4 f 3 + g 3 f 4 ) +21(-g 7 f 4 -g 4 f 7 + g 7 f 3 + g 3 f 7 -g 6 f 5 -g 5 f 6 + g 6 f 3 + g 3 f 6 ) +21(g 5 f 2 + g 2 f 5 + g 5 f + gf 5 -g 4 f -gf 4 -g 3 f 2 -g 2 f 3 ) +42(g 6 f 6 + g 2 f 2 ) -182g 4 f 4 Like the Schläfli equations, the polynomials are symmetric in f and g, but as can be seen from these and further examples, they are no more sparse.

3 .

 3 Generalised Schläfli equations. As argued in Section 4.