Max Wiebicke
email: max.wiebicke@tu-dresden.com

Edward Andò

Václav Šmilauer

Ivo Herle

Gioacchino Viggiani

A benchmark strategy for the experimental measurement of contact fabric

Keywords: image analysis, fabric, interparticle contacts, DEM, x-ray CT

The mechanics of granular materials can be better understood by experimental measurement of fabric and its evolution under load. X-ray tomography is a tool that is increasingly used to acquire threedimensional images and thus, enables such measurements.

Our previous study on the metrology of interparticle contacts revealed that the most common approaches either fail to accurately measure contact fabric or introduce a strong bias. Methods to improve these measurements (i.e., the detection and orientation of contacts) were proposed and validated.

This work develops a strategy to benchmark image analysis tools that can be used for the determination of contact fabric from tomographic images. The discrete element method is used to create and load a reference specimen for which the fabric and its evolution is precisely known. Chosen states of this synthetic specimen are turned into realistic images taking into account inherent image properties, such as the partial volume effect, blur and noise. The application of the image analysis tools on these images validates the findings of the metrological study and highlights the importance of addressing the identified shortcomings, i.e., the sys-

1 Introduction
Historically, soils are most commonly treated in a continuum framework in either theoretical modeling or the interpretation of experimental results. For most of the time, soil mechanicians were only able to observe and measure soil behaviour on the boundaries of a studied specimen. It was impossible to look into a specimen while loading it. It was acknowledged, however, that the micro-mechanics, i.e., the behaviour at the grainscale, drives the overall response of any granular material. Model granular materials, e.g., coins, provided first insights [START_REF] Wiendieck | Zur Struktur körniger Medien[END_REF]. In these and following experiments, two-dimensional (planar) grains were loaded and photographs were acquired at specific stages of the loading [START_REF] Calvetti | Experimental micromechanical analysis of a 2d granular material: relation between structure evolution and loading path[END_REF]. Using photogrammetry [START_REF] Desrues | Strain localization in sand: an overview of the experimental results obtained in Grenoble using stereophotogrammetry[END_REF], insight into the micromechanical behaviour was gained for either 2D model materials or granular materials in specific loading conditions that allowed photos to be acquired, e.g., in a biaxial loading apparatus.

Experiments on real granular materials were restricted to post-mortem analysis [START_REF] Oda | Initial fabrics and their relations to mechanical properties of granular material[END_REF][START_REF] Oda | Stress-induced anisotropy in granular masses[END_REF]: similar, yet not identical, specimens were loaded to specific stages and subsequently frozen in this state, e.g., using epoxy resin. This allowed for thin cuts to be taken and analysed under a microscope.

On the simulation front, the development of the discrete element method (DEM) has allowed micro-mechanical investigations of the particle-scale behaviour of granular materials [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF]. In these simulations, the response of a granular material is mainly controlled by the contact models which describe the mechanical interactions between particles. Applications of the discrete element method targeting soil behaviour have produced numerous valuable results and a better understanding of the micro-mechanisms linked to the overall behaviour [START_REF] Dobry | Discrete modelling of stressstrain behaviour of granular media at small and large strains[END_REF][START_REF] O'sullivan | Micromechanics of granular material response during load reversals: Combined DEM and experimental study[END_REF]. Despite continuously improving contact models and ever evolving DEM frameworks (e.g., [START_REF] Kawamoto | Level set discrete element method for three-dimensional computations with triaxial case study[END_REF]), DEM remains a numerical tool that needs further experimental validation. Its contact models are phenomenological, grain shapes are restricted (for an exception, see e.g., [START_REF] Kawamoto | Level set discrete element method for three-dimensional computations with triaxial case study[END_REF]) and surface properties of natural materials cannot be directly modeled.

With the rise of x-ray computed tomography (CT) in the last decade(s), full-field measurements of specimens finally became possible [START_REF] Desrues | Void ratio evolution inside shear bands in triaxial sand specimens studied by computed tomography[END_REF]. Specific set-ups in x-ray CT chambers allow for an acquisition of tomographies during a macroscopic loading of the specimen [START_REF] Hall | Discrete and continuum analysis of localised deformation in sand using X-ray µCT and volumetric digital image correlation[END_REF]. Different methods have been used and improved to measure the kinematics of the granular materials during loading, such as the digital image correlation [START_REF] Butterfield | A stero-photogrammetric method for measuring displacement fields[END_REF][START_REF] Hall | Discrete and continuum analysis of localised deformation in sand using X-ray µCT and volumetric digital image correlation[END_REF][START_REF] Andò | Experimental micro-mechanics of granular media studied by x-ray tomography: recent results and challenges[END_REF] or an individual particle tracking approach [START_REF] Andò | Grain-scale experimental investigation of localised deformation in sand: a discrete particle tracking approach[END_REF]. These experiments and tools allowed new insights into the mechanics at the grain-scale, in particular in the formation of shear bands [START_REF] Desrues | Void ratio evolution inside shear bands in triaxial sand specimens studied by computed tomography[END_REF][START_REF] Alshibli | Microscopic evaluation of strain distribution in granular materials during shear[END_REF][START_REF] Andò | Experimental micro-mechanics of granular media studied by x-ray tomography: recent results and challenges[END_REF]. Only recently, attempts to measure the granular fabric and its evolution at specific loading stages have been conducted [START_REF] Fonseca | Quantifying the evolution of soil fabric during shearing using directional parameters[END_REF][START_REF] Fonseca | Image-based investigation into the primary fabric of stress-transmitting particles in sand[END_REF][START_REF] Druckrey | 3D characterization of sand particle-to-particle contact and morphology[END_REF]. There are several inherent features in x-ray tomographic images, that influence the analysis of such images, e.g., the detection of contacts and the determination of their orientation. Regarding the measurement of fabric the most problematic and apparent of these inherent features are blur, noise and the partial volume effect (PVE) [START_REF] Cnudde | High-resolution X-ray computed tomography in geosciences: A review of the current technology and applications[END_REF][START_REF] Weis | Analyzing X-Ray tomographies of granular packings[END_REF]. Despite its significance, a quantitative analysis of the accuracy for the tools that are used to determine the fabric, i.e., contact normal orientations, principal particle orientations or void orientations, has not been undertaken until recently [START_REF] Wiebicke | On the metrology of interparticle contacts in sand from x-ray tomography images[END_REF]. This study on the metrology of individual contacts focused only on two particles being in contact instead of big granular assemblies. The accuracy of the contact detection as well as the determination of their orientation were quantitatively demonstrated. Advanced techniques to improve the measurements were proposed.

The aim of this contribution is to extend the findings of [START_REF] Wiebicke | On the metrology of interparticle contacts in sand from x-ray tomography images[END_REF] on individual contacts to a granular assembly that evolves under a macroscopic loading in order to assess what amplitude of fabric changes can be measured with realistic granular rearrangements. One main concern is the detection of small fabric changes that can occur e.g., upon unloading and reloading cycles. As these measurements and their quantification are difficult on real experiments e.g., using x-ray tomography, we start to study the problem from synthetic, but realistic1 images. The advantage of these images is the controllability of the inherent features as well as the perfect knowledge of the fabric. For this purpose, the discrete element method is used to create and load a numerical, mechanically representative specimen. Realistic images of these DEM specimens are then created using Kalisphera [START_REF] Tengattini | Kalisphera: an analytical tool to reproduce the partial volume effect of spheres imaged in 3D[END_REF] which generates analytically-correct 3D images of spheres.

Sphere packings are the object of many studies on the description of granular media [START_REF] Fox | The Pursuit of Perfect Packing[END_REF]. One of the most studied properties is the structure of these packings, e.g., a random loose packing and a jamming transition [START_REF] Song | A phase diagram for jammed matter[END_REF]. Random packings have been studied from frictionless to frictional spheres, monodisperse to polydisperse particle sizes under compression and shear to predict these structural measures [START_REF] Bi | Jamming by shear[END_REF][START_REF] Baule | Fundamental challenges in packing problems: From spherical to nonspherical particles[END_REF]. We decided to limit our analysis to the most common contact fabric measures for natural granular materials in geomechanics, which mostly deals with frictional, polydisperse and irregularly shaped grains: the coordination number and the fabric tensor.

The main object of this work is to present a benchmark test for the extraction of contact fabric from tomographies of granular materials. Our approach is somehow similar to [START_REF] Delaney | Combining tomographic imaging and DEM simulations to investigate the structure of experimental sphere packings[END_REF], where DEM simulation were conducted based on tomographic data from sphere packings. Their aim, however, was the investigation of the mechanics in contrast to ours which is the accuracy of the fabric evaluation from tomographic images. Although we employ sphere packings, the methodology is not restricted to spherical, but aimed at more general shapes. Thus, the image analysis tools never assume a particular particle shape to support the extraction of information from the tomographies. Another similar approach was used to determine constriction sizes of pores from DEM simulations [START_REF] Shire | Measurement of constriction size distributions using three grainscale methods[END_REF]. Binary images were created from the numerical sphere packings and processed by watershed segmentation. The results were compared with the constriction sizes determined directly on the DEM specimen.

All image-based calculations are done with the opensource Software spam (The Software for the Practical Analysis of Materials) [START_REF] Andò | and the spam contributors. spam: The software for the practical analysis of materials[END_REF]. The scripts to reproduce the results presented in this work are available in an open repository [START_REF] Wiebicke | Benchmark analysis of synthetical imagessource code and example dem data[END_REF].

Contact-based measurements

In the following section, the image analysis steps that are necessary to extract the contact fabric from x-ray tomographies are presented and the challenges of such analyses are highlighted. A simplified workflow of a general image analysis is presented in Figure 1.

The first step of producing an x-ray CT is the acquisition of radiographies of the specimen at various angles. These acquired radiographies serve as the basis for the reconstruction of a three-dimensional grey-scale image, the tomography.

In order to determine contact orientations, contacts have to be detected first. Therefore, the grey-scale image is transformed into a binary image where only two phases are distinguished -solid and void. This can be achieved using a variety of methods (global, local, adaptive thresholds, etc.) as described in [START_REF] Iassonov | Segmentation of X-ray computed tomography images of porous materials: A crucial step for characterization and quantitative analysis of pore structures[END_REF]. The most common and easiest method is to apply a global threshold, which attributes all grey values up to a certain threshold to either one of the two phases. The calibration can be based on physical quantities (such as the actual volume of the solid phase if known) or image statistics (such as Otsu's threshold [START_REF] Otsu | A Threshold Selection Method from Gray-Level Histograms[END_REF]). The next essential step in identifying the individual grains is the segmentation of the binary image. The most common segmentation techniques (in experimental geomechanics) are watershed algorithms, e.g., a standard topological watershed [START_REF] Meyer | Morphological segmentation[END_REF] or advanced watersheds, such as the random walker [START_REF] Grady | Random walks for image segmentation[END_REF]. Once segmented, the image can then be labelled by assigning a unique number to each individual grain. Based on the segmented and labelled image, the contacts between apparently touching particles can be detected and the contact planes can be fitted (e.g., via principal component analysis) through the contact area. The contact orientation can be identified as the normal to that plane.

The accuracy of the crucial steps is analysed in [START_REF] Wiebicke | On the metrology of interparticle contacts in sand from x-ray tomography images[END_REF]. The main findings are that the standard image analysis tools, briefly described above, are not sufficient to produce reliable results and that these techniques have either to be adapted or replaced by more advanced methods to detect contacts and determine their orientation as reliable as possible.

Major results of the study on the metrology

To methodically study the accuracy of the contact detection and orientation, the investigation started with analysing synthetic images of spheres, where particle shape and image properties such as blur and noise can be controlled, and followed up with high resolution to-mographies of manufactured spheres (to control the shape) and a natural granular material (Hostun sand).

The main source of uncertainty in contact detection is the binarisation of the grey-scale image and, particularly, the PVE [START_REF] Weis | Analyzing X-Ray tomographies of granular packings[END_REF]. The approach that was chosen to determine the accuracy of contact detection is described in detail in [START_REF] Wiebicke | On the metrology of interparticle contacts in sand from x-ray tomography images[END_REF]. Two particles, that are initially in contact, were incrementally translated from each other to check up to which distance they still appear to be in contact in the image. The results of the analysis with the most common approach (the global threshold) are shown in Figure 2 with solid lines. A perfect detection would appear as 1 at 0 distance and zero for distances > 0, i.e., with a probability of finding contacts of one at zero distance and zero at higher distances. An obvious over-detection can be observed in Figure 2. The main reason for the systematical over-detection of contacts is the PVE in regions where two particles are close to each other. In order to minimise the over-detection, a local thresholding is proposed: contacts that are detected with the global threshold are checked again with a higher, local threshold that is applied directly on the grey-values in the apparent contact region. It is found that this approach works well for rounded particles, less for more angular particles, e.g., Hostun sand grains. The results of chosen local thresholds are plotted in Figure 2 in dotted lines.

The same images, i.e., synthetic images of spheres as well as high resolution images of manufactured spheres and Hostun sand grains, were taken for the study on the accuracy of the determination of contact orientations. The contact orientations were determined by fitting a contact plane on the contact zone using principal component analysis. The two different watershed methods mentioned above, a topological watershed and the random walker [START_REF] Grady | Random walks for image segmentation[END_REF] (as implemented in [START_REF] Van Der Walt | and the scikit-image contributors. scikit-image: image processing in Python[END_REF]), were used for the segmentation of the binary image. The topological watershed results in a labelled image in which the contact region can be identified from the labels in direct contact. The random walker assigns a probability to each voxel of either belonging to one or the other particle in contact which allows to determine the contact surface from the probability field on a subvoxel basis, whereas the former method can only work on the voxel basis. The error is defined as the angle between a reference and the determined orientation.

Spheres were chosen for this and the former analysis because the reference contact orientation is known from the branch vector, that is either imposed in the synthetic images or precisely measured from the center of mass of the two apparently touching particles [START_REF] Andò | Experimental investigation of microstructural changes in deforming granular media using x-ray tomography[END_REF]. The reference orientation of the contacts in Hostun sand was defined as the orientation that can be determined at the Topological Watershed Random Walker Fig. 3 Mean error of contact orientation for a diameter of 27 pixel -from the analysis in [START_REF] Wiebicke | On the metrology of interparticle contacts in sand from x-ray tomography images[END_REF] smallest pixel size (1 µm)2 . The results of this part of the study are shown in Figure 3 for a relative resolution d 50 /pixel size of 27 pixel. The results are displayed for all materials in the analysis: the synthetic images of spheres, the high resolution tomograhies of manufactured spheres and the high resolution tomographies of Hostun sand. It can be concluded that the random walker yields more accurate orientations for all materials tested when the focus lies on individual contacts, i.e., only two particles touching each other.

Binarisation Segmentation and labelling

Grey-scale image Binary image Labelled image

Image analysis tools

All the image analysis tools mentioned above that were either tested or proposed to overcome problems in the description of inter-particle contacts were included in the open source analysis software spam [START_REF] Andò | and the spam contributors. spam: The software for the practical analysis of materials[END_REF]. It has evolved to cover needs of data analysis from 3D x-ray tomography work and correlated random fields with mechanical applications. The contact detection and orientation tools from [START_REF] Wiebicke | On the metrology of interparticle contacts in sand from x-ray tomography images[END_REF] are included in the label toolkit of spam.

The complete image analysis that is applied to obtain the results in this work is summarised in minimal python scripts using spam. In order to enable the reproduction of this image analysis, these scripts are available online in a free repository [START_REF] Wiebicke | Benchmark analysis of synthetical imagessource code and example dem data[END_REF]. The source data from the discrete element simulations (positions and radii of the spheres at every step of the analysis) can be found in the repository as well, and should be easy to obtain with any 3D DEM code.

Testing and validating the measurement tools

The study on the metrology focused completely on the detection and orientation of individual contacts, i.e., only two particles being in contact. In this work we develop a benchmark test in order to assess whether the image analysis tools work on images of assemblies of thousands of particles and how accurate the results are. In order to assess the accuracy, the reference, i.e., the true fabric, has to be known. The discrete element method is used to create assemblies of spherical particles and to obtain particle displacements due to externally imposed mechanical loads. Realistic grey-scale images of these numerical assemblies at different loading stages are then created using Kalisphera [START_REF] Tengattini | Kalisphera: an analytical tool to reproduce the partial volume effect of spheres imaged in 3D[END_REF] in the same way as in [START_REF] Wiebicke | On the metrology of interparticle contacts in sand from x-ray tomography images[END_REF] for individual particles in contact.

It has to be noted that due to the spherical shape of the particles branch vectors could be used to precisely determine the contact orientations. We employ the chosen image analysis approach because it has to be extendible to round (not necessarily spherical) particle shapes for which the equality of branch vector and contact orientation does not hold in general. The accuracy of the determination of the contact orientations of round and spherical particles, however, should be similar as the study on the metrology showed that the accuracy follows the contact topology.

As mentioned in the introduction, there exist various methods to detect contacts. Three different methods of estimating contact numbers were compared in [START_REF] Aste | Local and global relations between the number of contacts and density in monodisperse sphere packs[END_REF] on monodisperse sphere packings. The best results were obtained using a contact number scaling function presented in [START_REF] Aste | Geometrical structure of disordered sphere packings[END_REF], which however, is assuming spherical particles shapes. This approach was extended to packings of ellipsoids in [START_REF] Schaller | Tomographic analysis of jammed ellipsoid packings[END_REF] and [START_REF] Schaller | Local origin of global contact numbers in frictional ellipsoid packings[END_REF]. As our research is generally dealing with irregular shapes and polydisperse packings, our approaches cannot profit from the mentioned results.

DEM simulations as a reference

The open source, discrete element software WooDEM [START_REF] Šmilauer | Woo Documentation[END_REF] is used to create the numerical specimen and to load it mechanically. In order to directly consider the findings from [START_REF] Wiebicke | On the metrology of interparticle contacts in sand from x-ray tomography images[END_REF], the particles are chosen to have a spherical shape. The particle size distribution (PSD) of the specimen is linear and ranges from (3.2 to 4.8) • 10 -4 m and is thus in the range of the materials that were investigated in [START_REF] Wiebicke | On the metrology of interparticle contacts in sand from x-ray tomography images[END_REF] and that will be used in real experiments in the µ-CT: manufactured spheres with d = 4 • 10 -4 m and Hostun sand with d 50 = 3.38 • 10 -4 m. Periodic boundary conditions and a standard triaxial loading (axisymmetric) with a final unloading step are chosen for the simulation. The simulation itself is only needed to create mechanically meaningful states of the granular assembly for the subsequent image analysis and not to probe the mechanics. Thus, a simple, linear (Cundall) contact model [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF] is considered as sufficient. The parameters of the contact model are chosen to be: E = 50 • 10 6 kPa, k t /k n = 0.4 and tan ϕ = 0.4 with E being Young's modulus, k n , k t the normal and tangent stiffness respectively and ϕ the friction angle. The initial package of 5522 spheres based on the given PSD is created and then isotropically compressed to 50 kPa. At pre-defined axial strains the fabric of the assembly, i.e., the position of the particles and the contact orientations, is extracted to determine the evolution of the fabric upon loading. The macroscopic response, i.e., the stress-strain response calculated at the boundaries of the specimen, is plotted in Figure 4.

Turning DEM specimen into images

The chosen, exported fabric states are turned into realistic grey-scale images using Kalisphera [START_REF] Tengattini | Kalisphera: an analytical tool to reproduce the partial volume effect of spheres imaged in 3D[END_REF]. The pixel size, that is needed to convert the DEM data (positions and radii of spheres given in [m]) to actual, pixelised data for the the image (given in [pixel]), is set to 15 • 10 -6 m, which corresponds to a common setting of a tomography of a triaxial experiment on sand (e.g., [START_REF] Andò | Grain-scale experimental investigation of localised deformation in sand: a discrete particle tracking approach[END_REF]). A 32-bit floating point image of each sphere is created, where the values 0.0 and 1.0 correspond to a pure void and a pure solid voxel, respectively. The individual images are then subsequently added to an image of the big assembly. In some cases the overlap of particles on contacts leads to values higher than 1.0, which are flattened to 1.0. In order to apply blur and noise the values of the image are compressed to lie in the range of (0.25;0.75) corresponding to void and solid voxels, respectively. A Gaussian blur operator is applied to the image and random noise based on a Gaussian distribution is applied onto the blurred image. Representative levels of blur and noise for the µ-CT in Laboratoire 3SR in Grenoble were calibrated on chosen images in [START_REF] Wiebicke | On the metrology of interparticle contacts in sand from x-ray tomography images[END_REF]: the level of blur was found to be ≈ 0.8 and the standard deviation of the Gaussian distribution for random noise was chosen to 0.03. The final 32-bit floating point images are turned to 8-bit integer images (i.e., values between 0 and 255), because the image analysis on real tomographic images is usually carried out on this data type to save memory. The whole procedure is summarised in Figure 5.

Analysis of the reference and computed fabric

The advanced image analysis techniques, that were described above and in more detail in [START_REF] Wiebicke | On the metrology of interparticle contacts in sand from x-ray tomography images[END_REF], are applied to the resulting images. In order to reduce the overdetection of contacts, a local threshold of 180 (8bit unsigned integer images)3 is applied to the globally detected contacts. A ITK inter-pixel watershed as implemented in [START_REF] Andò | and the spam contributors. spam: The software for the practical analysis of materials[END_REF] is applied to the binary image for the initial segmentation and labelling. The contact orientations are determined using a principal component analysis on the contacting voxels in the initially labelled image. Furthermore, the random walker [START_REF] Grady | Random walks for image segmentation[END_REF] is then employed on contacting particles for the determination of the contact orientations from the resulting probability map.

A scalar value describing the fabric of the granular assembly is the coordination number, i.e., the number of contacts per particle. The average coordination number C n is given by:

C n = 2 • C N p (1)
with C and and N p being the number of contacts and the number of particles, respectively. The coordination number is directly determined from the amount of detected contacts and should thus reflect the findings from [START_REF] Wiebicke | On the metrology of interparticle contacts in sand from x-ray tomography images[END_REF], namely the systematic over-detection of contacts. The evolution of the average coordination number with the macroscopic loading is plotted in Figure 6. As to be expected, the number of contacts from the images is dramatically over-detected using the standard global approach with a mean deviation of 50%. The local refinement of the detected contacts leads to a coordination number that is close to the reference from the simulation with a mean deviation of 3%. Even though a local refinement is employed, some contacts are still over-detected. This is because the local threshold was chosen in order not to loose real contacts, see the dotted orange line in Figure 2 and thus, some apparent contacts are always over-detected. Contact orientations can be plotted in different ways. Herein, we choose Lambert azimuthal equal area projections for the representation of orientations. In order to use this projection, each orientation (vector) is flipped such that it has a positive vertical axis. A simplified sketch of how these orientations are plotted and binned (when considering thousands of orientations) is shown in Figure 7.

The resulting plots of orientations at different, significant stages of the loading are shown in Figure 8: the reference orientations, that are directly exported from the DEM simulation are given in the upper row. The middle row shows the orientations that are calculated directly on the initial segmentation by ITK. The lower row depicts the orientations that are determined from the random walker segmentation of two contacting particles. In order to use similar scales for every loading state and to facilitate the comparison, the plots are normalised with the average number of orientations per bin.

The distribution of contact orientations in the DEM simulation (reference) evolve from a rather random configuration in the initial and the pre-stress peak stage to a clear vertically aligned contact fabric with ongoing loading.

The orientations from the initial segmentation and labelling do not resemble the reference distribution. As found in [START_REF] Jaquet | Estimation of Separating Planes between Touching 3D Objects Using Power Watershed[END_REF][START_REF] Viggiani | Identifying and following particle-to-particle contacts in real granular media: An experimental challenge[END_REF], orientations that are calculated on the basis of pixel watersheds are strongly biased towards favourable orientations. In this specific case of the ITK watershed the orientations tend towards ≈45 • with respect to the Cartesian axes. This was found to be the reason for the high error in the determination of contact orientations [START_REF] Wiebicke | On the metrology of interparticle contacts in sand from x-ray tomography images[END_REF], see Figure 3.

The contact orientations determined based on the random walker segmentation agree with the reference distribution. The reason for the better agreement is the probability map from the random walker, which allows interpolating the contact zone within the pixel instead of following pixel boundaries.

To quantitatively describe the errors of the image analysis approaches, the individual orientations must be related to their corresponding reference from the DEM simulation. As each contact is identified by the labels of the contacting spheres and the labels in the DEM simulations are not equivalent to the labels in the images, a topographic mapping from DEM to image is established. The error of each orientation is determined by computing the angle between the orientation measured in the image and its corresponding reference in the simulation. It has to be noted that not all orientations could be compared due to the problematic detection of contacts in the image which was described before and plotted in Figure 6. A histogram of the errors of all orientations in state 2 (as highlighted in Figure 8) is plotted in Figure 9. The errors are distributed similarly in all states. As expected and found in [START_REF] Wiebicke | On the metrology of interparticle contacts in sand from x-ray tomography images[END_REF], the random walker yields more accurate results with a statistical mean error of ≈ 1.3 • than the topological watershed with ≈ 4.2 • . The results from the random walker agrees well to the error of 1.24 • determined in the more rigorous study and plotted in 3. The error of the topological watershed is lower than 5.8 • , which was determined in [START_REF] Wiebicke | On the metrology of interparticle contacts in sand from x-ray tomography images[END_REF]. This might be due to larger contact areas which arise from the inherent overlap of the spheres in the numerical simulation and the different distribution of orientations.

In order to investigate the evolution of the individual orientations, a different, scalar value has to be chosen for the comparison. The mean inclination of all orientations is determined at each loading stage of the simulation, because the main changes happen to the vertical component for this particular loading. The inclination of each individual orientation is calculated by: θ = arccos z

x 2 + y 2 + z 2 (2)
The evolution of the mean inclination with respect to the axial strain is plotted in Figure 10. The same trend as qualitatively shown in Figure 8, i.e., an increase in the vertical orientations, can be observed qualitatively using this representation. The evolution of both -the reference and the orientations from the images (ITK and random walker) -follows the same trend with some minor deviations. Upon unloading, the mean inclination partly recovers, i.e., increases, for the reference and the image analysis results.

Although, the individual orientations are strongly biased and do not resemble the reference distribution, the mean inclination of all orientations qualitatively follows and even quantitatively agrees with the reference evolution with deviations not exceeding 0.7 • .

Fabric tensors can be used in order to statistically describe the bulk of orientations at a specific state of a granular material [START_REF] Kanatani | Distribution of directional data and fabric tensors[END_REF]. Most commonly, a fabric tensor of 2nd order is used:

N = 1 C C α=1 o α ⊗ o α (3)
with C being the number of contacts and o an individual orientation. Following [START_REF] Kanatani | Distribution of directional data and fabric tensors[END_REF], a deviatoric fabric tensor D is calculated by:

D = 15 2 (N - 1 3 I) (4)
To quantitatively describe anisotropy, a scalar value can be calculated from D [START_REF] Gu | Anisotropy of elasticity and fabric of granular soils[END_REF]:

a = 3 2 D : D (5)
A fabric tensor is computed at each stage of the simulation for both, the DEM fabric and the orientations determined by image analysis. The evolution of the scalar anisotropy factor a is plotted in Figure 11. As expected from the preparation of the synthetic specimen and the projection plots in Figure 8, the fabric is initially isotropic (a ≈ 0) with an increasing anisotropy due to the vertical alignment of contacts into the loading direction. The anisotropy then increases as the contact orientations realign with the induced principal stress direction, i.e., the vertical axis.

The results of the random walker and of the ITK watershed are shown in Figure 11. Although the topological watershed yields substantially higher errors and a strong bias when determining individual orientations (see Figures 8 and9), the fabric anisotropy agrees to the one determined by the random walker. Both evolutions follow the reference fabric qualitatively well and show only minor quantitative deviations.

The reason for the agreement of the anisotropy from the ITK watershed with the random walker as well as the reference is the degree of approximation of the distribution. Even though the orientations from the ITK watershed are strongly biased towards 45 • to the general axes, see Figure 8, the distribution is symmetric to/about the vertical axis. The same holds for the orientations from the random walker and the reference. The orientations from ITK concentrate along the vertical direction with ongoing loading like the orientations from the reference and the random walker. Thus, the anisotropy changes only with the vertical alignment, as the radial distribution remains almost isotropic. We employ a typical 2nd order approximation which is not sufficient to capture the bias from the ITK watershed, because it is radially symmetric. A higher order approximation should be able to show the strong differences of the distributions.

Results from a cyclic simulation

Another question that arises is the sensitivity of the contact fabric and its measurability upon small macroscopic changes, such as unload-reload cycles. In order to investigate this, a DEM simulation similar to the one presented in Figure 4 is carried out. The constitutive parameters are identical and the specimen was prepared similarly. The specimen is unloaded and reloaded seven times at different states throughout the loading. Three intermediate states are recorded at every cycle to create the synthetic images and carry out the image analysis that is described above. The macroscopic response of the specimen is shown in Figure 12.

The distribution of contact orientations is again described and compared using the scalar anisotropy factor, see equation [START_REF] Oda | Stress-induced anisotropy in granular masses[END_REF]. Figure 13 shows the evolution of the anisotropy of the reference as well as both image analysis approaches, the ITK watershed and the random walker. Both -the anisotropy from the ITK watershed and the random walker -agree with the reference. On the right side of Figure 13 a zoom-in of the first three cycles is displayed. Both approaches are able to capture even these small changes of the contact anisotropy.

Summary and outlook

The behaviour of granular materials is known to be governed by the mechanics at the grain-scale. According to [START_REF] Calvetti | Experimental micromechanical analysis of a 2d granular material: relation between structure evolution and loading path[END_REF], the micro-mechanical behaviour may be described by (a) the fabric, (b) the kinematics and (c) the inter-granular forces. A number of contributions have introduced and improved the measurement of continuous and discrete kinematics (b). The measurement of inter-granular forces (c) is still not possible in assemblies of natural granular materials, but recent developments show promising results on small assemblies of synthetic materials [START_REF] Hall | Can intergranular force transmission be identified in sand? First results of spatially-resolved neutron and X-ray diffraction[END_REF][START_REF] Hurley | Extracting inter-particle forces in opaque granular materials: Beyond photoelasticity[END_REF][START_REF] Hurley | Force measurements in stiff, 3D, opaque granular materials[END_REF][START_REF] Imseeh | 3D finite element modelling of force transmission and particle fracture of sand[END_REF]. Some studies reported on the measurement of granular fabric [START_REF] Imseeh | 3D experimental quantification of fabric and fabric evolution of sheared granular materials using synchrotron microcomputed tomography[END_REF][START_REF] Fonseca | Image-based investigation into the primary fabric of stress-transmitting particles in sand[END_REF], but without studying the metrology of the tools that were used for the measurements. In the present paper, the results of the study on individual contacts [START_REF] Wiebicke | On the metrology of interparticle contacts in sand from x-ray tomography images[END_REF] have been related and applied to an evolving assembly of thousands of particles.

To achieve that, a discrete element simulation is carried out to provide a reference, i.e., a sequence of perfectly known fabrics at several states of a macroscopic loading. Images of these states are then created using Kalisphera [START_REF] Tengattini | Kalisphera: an analytical tool to reproduce the partial volume effect of spheres imaged in 3D[END_REF], a tool that allows creating images of spheres taking into account the PVE. The other inherent defects, i.e., blur and noise, are applied to obtain realistic and representative synthetic images for the image analysis. The results are compared in terms of detected contacts, their orientations and a second order fabric tensor, that is commonly used to statistically express orientational fabric.

In our previous study [START_REF] Wiebicke | On the metrology of interparticle contacts in sand from x-ray tomography images[END_REF], it was found that contacts are systematically over-detected using standard global thresholding approaches, i.e., particles, that are close to each other but do not touch, are detected to be in contact due to the PVE. This is confirmed for the image of the synthetic assembly by running the global contact detection: ≈ 50% more contacts are detected in the image than actually present in the reference. The local thresholding approach, that was proposed in [START_REF] Wiebicke | On the metrology of interparticle contacts in sand from x-ray tomography images[END_REF], is applied to refine the detected contacts and is found to get much closer to the reference, with an over-detection of only ≈ 3%.

The contact orientations calculated on the basis of two different watersheds, i.e., a ITK watershed, which is used for the initial segmentation, and the random walker segmentation as a refinement, are compared either using plots of all individual orientations or a mean inclination (angle with the vertical axis) per analysed state. The orientations from the ITK watershed are strongly biased and only weakly show a tendency to the vertical axis with ongoing loading. On the contrary, the orientations based on the random walker segmentation agree with the reference distribution and show a similar trend/evolution throughout the loading. However, this representation is not suitable for a quantitative comparison and the mean inclination of all orientations is thus chosen as a scalar descriptor. Surprisingly, the evolution of the mean inclination obtained from the ITK watershed agrees to the reference and shows a slight deviation of max. 0.7 • . The mean inclination from random walker images follows a similar trend as for the reference; quantitatively, a deviation of ≈0.5 • is found.

The orientations can also be statistically captured by a second-order fabric tensor. Comparing the evolution of the anisotropy of the fabric tensors shows a qualitative agreement of both watershed methods with the reference fabric. Quantitatively, the fabrics from the images agree very well with each other, but show a slight deviation from the reference.

The findings of [START_REF] Wiebicke | On the metrology of interparticle contacts in sand from x-ray tomography images[END_REF] can thus be extended to granular assemblies of spherical particles: a local contact detection is necessary to reduce the over-detection and the random walker segmentation yields qualitatively as well as quantitatively acceptable results. The ITK watershed introduces a strong bias on individual contact orientations, but performs as well as the advanced and more time-consuming random walker segmentation when focused on the evolution of anisotropy for the chosen loading and the order of approximation of the distribution. The latter observation, however, arises due to additional effects, such as the symmetric loading and fabric evolution, as highlighted in the analysis of the results. For different non-symmetric loading situations or a higher order of approximation, the high error and bias of the orientations from the ITK watershed should likely be reflected in the anisotropy as well. Thus, it is still highly recommended to use the random walker segmentation (or other advanced algorithms) to determine contact orientations.

Based on this virtual simulation we expect to be able to capture cyclic behaviour just as well in an x-ray tomography of spherical, and round particles. This will allow inroads into understanding the abrupt stiffness changes observed upon unloading and reloading.

This analysis is restricted to the measurement of contact fabric. A complete description of the geometric fabric of a granular material has also to include other parts, such as particle [START_REF] Oda | Initial fabrics and their relations to mechanical properties of granular material[END_REF] and void fabric [START_REF] Oda | Stress-induced anisotropy in granular masses[END_REF] or even combined fabric measures such as presented in [START_REF] Matthew | Stressinduced anisotropy in granular materials: fabric, stiffness, and permeability[END_REF]. Some of these measures might be extracted more accurately from tomographic images and will be part of future studies.

Fig. 1 Fig. 2

 12 Fig. 1 Simplified sketch of the image analysis steps -from a grey-scale image (left) to a binary image (centre) to a segmented and labelled image (right).

Fig. 4 Fig. 5

 45 Fig. 4 Macroscopic response of a triaxial specimen from the discrete element simulation. The chosen loading states for the creation of the synthetic images are indicated in the upper left figure with grey boxes. Four states are given by their ID to reference them specifically in the analysis.

Fig. 6

 6 Fig.6Plot of the evolution of the average coordination number with the macroscopic loading. Comparison of the reference from the simulation and the results of the global and local contact detection of the image.

Fig. 7

 7 Fig. 7 Visualisation of the plotting of individual orientations and the binning to simplify plots of many orientations. A Lambert azimuthal equal area projection is used in this work.

Fig. 8

 8 Fig.8Binned Lambert azimuthal equal area plots of the contact orientations at different stages of the loading (the labels relate to the macroscopic response in Figure4) -upper row: orientations directly from the discrete element simulation. middle row: orientations determined from the ITK watershed. lower row: orientations determined with the random walker segmentation. The colors correspond to the normalised density of orientations in the bins -red depicts a high and blue a low density. State 2 is highlighted for a more detailed analysis.

Fig. 9 Fig. 10

 910 Fig.9Histogram of the errors of the individual orientations determined by the ITK watershed and the random walker for state 2 -highlighted in Figure8. The errors represent the angles between the orientations from the respective image analysis technique and the reference from the DEM simulation.

Fig. 11

 11 Fig. 11 Evolution of the scalar anisotropy factor with the loading. Comparison of the image analysis results with the reference behaviour from the DEM simulations.

Fig. 12 Fig. 13

 1213 Fig. 12 Macroscopic response of the synthetic specimen to unload / reload cycles.

In this case realistic images are images that have similar properties as x-ray tomographies.

Note: Due to the angular and non-convex shape of the Hostun sand grains, the description of the orientation however is not objective, as highlighted in[START_REF] Wiebicke | On the metrology of interparticle contacts in sand from x-ray tomography images[END_REF]: when observing the same grains in different spatial resolutions, contacts might artificially be created (due to the PVE) or multiple contacts points at which two grains are touching might be merged to appear as a single contact.

This value corresponds to 0.7 for 32bit floating point images, which was found to yield acceptable results in[START_REF] Wiebicke | On the metrology of interparticle contacts in sand from x-ray tomography images[END_REF].

Acknowledgements We express our thanks to Félix Bertoni for implementing Kalisphera in C++. Laboratoire 3SR is part of the LabEx Tec 21 (Investissements dAvenir -grant agreement n • ANR-11-LABX-0030). We thank the Center for Information Services and High Performance Computing (ZIH) at TU Dresden for generous allocations of computing resources.

Compliance with Ethical Standards

Conflict of Interest: The authors declare that they have no conflict of interest.

Funding: The research leading to these results has received funding from the German Research Foundation (DFG) n • 254872581 and from the European Research Council under the European Union's Seventh Framework Program FP7-ERC-IDEAS Advanced Grant Agreement n • 290963 (SOMEF).


Y Abdelaziz

C Koutschan

J-M Maillard

On Christol's conjecture

Keywords: AMS Classification scheme numbers: 33C05, 33C20, 33F10, 68W30 Christol's conjecture, diagonals of rational functions, Shimura curves, creative telescoping, telescopers, D-finite series, globally bounded series

We show that the unresolved examples of Christol's conjecture 3 F 2 ([2/9, 5/9, 8/9], [2/3, 1], x) and 3 F 2 ([1/9, 4/9, 7/9], [1/3, 1], x), are indeed diagonals of rational functions.

We also show that other 3 F 2 and 4 F 3 unresolved examples of Christol's conjecture are diagonals of rational functions. Finally we give two arguments that show that it is likely that the 3 F 2 ([1/9, 4/9, 5/9], [1/3, 1], 27 • x) function is a diagonal of a rational function.

Introduction

There is a plethora of multiple integrals in physics: Feynman integrals, lattice Green functions, the summands of the magnetic susceptibility of the 2D Ising model [START_REF] Bostan | Zenine High order Fuchsian equations for the square lattice Ising model: χ(5)[END_REF][START_REF] Boukraa | High order Fuchsian equations for the square lattice Ising model: χ(6)[END_REF], that have very specific mathematical properties. These functions are D-finite, i.e., solutions of linear differential operators with polynomial coefficients, and have series expansions with integer coefficients. Furthermore, it was also shown that the linear differential operators annihilating the summands of the magnetic susceptibility of the Ising model χ(n), verify the specific property of being Fuchsian † operators: the critical exponents of all their singularities are given by rational numbers, and their Wronskians are N -th roots of rational functions [START_REF] Bostan | Zenine High order Fuchsian equations for the square lattice Ising model: χ(5)[END_REF][START_REF] Boukraa | High order Fuchsian equations for the square lattice Ising model: χ(6)[END_REF]. It was also shown that the χ(n) functions are solutions of globally nilpotent operators [START_REF] Bostan | Zenine Globally nilpotent differential operators and the square Ising model[END_REF], and that they "come from geometry" being G-operators [START_REF] André | G-functions and geometry Aspects of Mathematics E[END_REF].

The unifying scheme behind these seemingly sparse properties is that these functions are diagonals of rational functions [START_REF] Bostan | Ising n-fold integrals as diagonals of rational functions and integrality of series expansions[END_REF][START_REF] Bostan | Ising n-fold integrals as diagonals of rational functions and integrality of series expansions: integrality versus modularity Preprint[END_REF]. It was shown for example in [START_REF] Bostan | Ising n-fold integrals as diagonals of rational functions and integrality of series expansions: integrality versus modularity Preprint[END_REF], that if summands of the magnetic susceptibility χ(n) for any n have an integer coefficient series expansion reducing to algebraic series modulo any prime, it is because they are diagonals of rational functions for any integer n. In fact many problems in mathematical physics involving n-fold integrals, could be interpreted in terms of diagonals of algebraic or rational ¶ functions †.

Gilles Christol has shown in [START_REF] Christol | Diagonales de fractions rationnelles et équations de Picard-Fuchs[END_REF] that for every rational function, its diagonal f (x) has the following properties:

• it is globally bounded: there exist integers c and

d in N * , such that d f (c x) ∈ Z[[x]] and f (x) has a radius of convergence that is non-zero in C. • D-finite: there exists a differential operator L ∈ Z[x][ d dx ], with L = 0, such that L(f ) = 0.
Christol conjectured in [START_REF] Christol | Fonctions hypergéométriques bornées[END_REF] that every series that verifies these two properties is the diagonal of a rational function. In this paper [START_REF] Christol | Fonctions hypergéométriques bornées[END_REF] Christol came up with an unresolved example to his conjecture [START_REF] Christol | Fonctions hypergéométriques bornées[END_REF], and a longer list was generated by Christol and his co-authors in 2012 in [START_REF] Bostan | Ising n-fold integrals as diagonals of rational functions and integrality of series expansions: integrality versus modularity Preprint[END_REF]. In this paper we show that two of the unresolved examples of the conjecture given in [START_REF] Bostan | Ising n-fold integrals as diagonals of rational functions and integrality of series expansions: integrality versus modularity Preprint[END_REF] on page 58, namely the 3 F 2 [2/9, 5/9, 8/9], [2/3, 1], 3 6 • x and 3 F 2 [1/9, 4/9, 7/9], [1/3, 1], 3 6 • x are indeed diagonals of rational functions and provide a generalization of this result.

Recalls on diagonals of rational functions and on Christol's conjecture

Definition of the diagonal of a rational function

The diagonal of a rational function in

n variables R(x 1 , . . . , x n ) = P(x 1 , . . . , x n )/Q(x 1 , . . . , x n ), where P, Q ∈ Q[x 1 , . . . ,
x n ] such that Q(0, . . . , 0) = 0, is defined through its multi-Taylor expansion around (0, . . . , 0):

R x 1 , . . . , x n = ∞ m1 = 0 • • • ∞ mn = 0 R m1, ..., mn • x m1 1 . . . x mn n , (1)
as the series in one variable x:

Diag R x 1 , . . . , x n = ∞ m = 0 R m, m, ..., m • x m . (2)

Hadamard product of algebraic functions and Christol's conjecture

Recall that the Hadamard product of two series f

(x) = ∞ n=0 α n • x n and g(x) = ∞ n=0 β n • x n is given by: f (x) g(x) = ∞ n=0 α n • β n • x n . ( 3
)
Hypergeometric series of the form

p F p-1 ([a 1 , . . . , a p ], [b 1 . . . , b p-1 ], x) of height h = h(a 1 , . . . , a p , b 1 . . . , b p-1 )
, where the height h is given by:

h = #{1 ≤ j ≤ p | b j ∈ Z} -#{1 ≤ j ≤ p | a j ∈ Z} (4)
with b p = 1, can be written ‡ as the Hadamard product of h globally bounded series of height 1, were shown to verify Christol's conjecture. For example, the globally bounded hypergeometric series 3 F 2 ([1/3, 1/3, 1/3], [START_REF] Bostan | Re: "Trivialement diagonale[END_REF][START_REF] Bostan | Re: "Trivialement diagonale[END_REF], x) has height 3, and it can be written as the Hadamard product of three algebraic functions †:

3 F 2 ([ 1 3 , 1 3 , 1 3 ], [1, 1], x) = (1 -x) -1/3 (1 -x) -1/3 (1 -x) -1/3 , (5)
and can thus be written as the diagonal of the algebraic function in three variables:

(1 -x) -1/3 • (1 -y) -1/3 • (1 -z) -1/3 . (6)
Unlike the case of 3 F 2 ([1/3, 1/3, 1/3], [1, 1], x) , the hypergeometric functions 3 F 2 ([2/9, 5/9, 8/9], [2/3, 1], x) and 3 F 2 ([1/9, 4/9, 7/9], [1/3, 1], x), while being globally bounded functions [START_REF] Christol | Globally bounded solutions of differential equations Analytic number theory[END_REF], were constructed in a way that avoids them being written as "simple" Hadamard products of algebraic functions. Note that a p F p-1 hypergeometric function is globally bounded without restrictions if all its parameters in the denominator are integers, while a p F p-1 hypergeometric function can be shown to be globally bounded in general, by looking at Landau functions as explained in the work of Christol [START_REF] Christol | Fonctions hypergéométriques bornées[END_REF]. Furthermore, Beukers and Heckman have shown in [START_REF] Beukers | Monodromy for the hypergeometric function nFn-1[END_REF], that p F p-1 hypergeometric functions of height zero that are globally bounded are algebraic functions.

Unresolved examples to the conjecture

Generalized hypergeometric functions with regular singularities p F p-1 are a simple and natural testing ground for Christol's conjecture.

All

2 F 1 ([a, b], [c],
x) hypergeometric series that are globally bounded are diagonals of rational functions. There are two cases that fall into this category:

• Either the parameter c is an integer, and the 2 F 1 function can be written as the Hadamard product of two 1 F 0 functions, which are algebraic functions, and thus are diagonals of rational functions by Furstenberg's [START_REF] Furstenberg | Algebraic functions over finite fields[END_REF] theorem ¶.

• The parameter c is not an integer, and in this case the 2 F 1 function is a diagonal of a rational function if and only if it is an algebraic function † † (and consequently its series is globally bounded).

Moving on to 3 F 2 hypergeometric functions, one can ask when is a 3 F 2 hypergeometric function a diagonal of a rational function?

• When the parameters d and e in 3 F 2 ([a, b, c], [d, e], x) are integers, because, in this case, the 3 F 2 can be written as the Hadamard product of three 1 F 0 algebraic functions, and is thus the diagonal of a rational function, by the closure of diagonals under the Hadamard product.

• When the parameters d and

e in 3 F 2 ([a, b, c], [d, e],
x) are rational numbers but not integers, because in this case the 3 F 2 is algebraic, and is thus a diagonal by Furstenberg's theorem. † Diagonals are closed under the Hadamard product: if two series are diagonals of rational functions, their Hadamard product is also a diagonal of a rational function.

¶ Furstenberg's theorem states that any algebraic function is the diagonal of a rational function in two variables. † † The only 2 F 1 hypergeometric functions that are globally bounded with c ∈ Q are the algebraic ones: they are the ones appearing in the list of Schwarz [START_REF] Schwarz | Ueber diejenigen Fälle, in welchen die Gaussische hypergeometrische Reihe eine algebraische Function ihres vierten Elementes darstellt[END_REF].

Hence the interesting case occurs when only one of the two parameters d or e is rational, and the other is an integer. But even in this case, a lot of the 3 F 2 functions are easily seen to be diagonals of a rational function. Suppose that a 3

• 2 F 1 ([a, b], [e], x) 1 F 0 ([c], x) • 2 F 1 ([a, c], [e], x) 1 F 0 ([b], x) • 2 F 1 ([b, c], [e], x) 1 F 0 ([a], x) • 2 F 1 ([a, b], [1], x) 2 F 1 ([c, 1], [e],
x)

• 2 F 1 ([a, c], [1], x) 2 F 1 ([b, 1], [e],
x)

• 2 F 1 ([b, c], [1], x) 2 F 1 ([a, 1], [e],
x)

Now 1 F 0 ([c], x) and 2 F 1 ([a, b], [1], x)
∈ Q \ Z. Now 2 F 1 ([c, 1], [e],
x) cannot be an algebraic functions for e ∈ Q by Goursat [START_REF] Goursat | Leccons sur les séries hypérgéométriques et sur quelques fonctions qui s'y rattachent[END_REF]. Hence if one of

2 F 1 ([a, b], [e], x), 2 F 1 ([b, c], [e], x), or 2 F 1 ([a, c], [e], x) is an algebraic function, then 3 F 2 ([a, b, c], [1, e],
x) is the diagonal of a rational function. Now taking the two examples given in [START_REF] Bostan | Ising n-fold integrals as diagonals of rational functions and integrality of series expansions: integrality versus modularity Preprint[END_REF] or by G. Christol in [START_REF] Christol | Fonctions hypergéométriques et diagonales de fonctions rationnelles[END_REF] that we are looking at here, we see that neither 3 F 2 ([2/9, 5/9, 8/9], [2/3, 1], x), nor 3 F 2 ([1/9, 4/9, 7/9], [2/3, 1], x), can be obtained as diagonals of rational functions through Hadamard products † since the three 2 F 1 hypergeometric series are not globally bounded ¶:

2 F 1 [ 2 9 , 5 9 ], [ 2 3 ], x , 2 F 1 [ 2 9 , 8 9 ], [ 2 3 ], x , 2 F 1 [ 5 9 , 8 9 ], [ 2 3
], x ,

and nor are the 2 F 1 hypergeometric series:

2 F 1 [ 1 9 , 4 9 ], [ 1 3 ], x , 2 F 1 [ 4 9 , 7 9 ], [ 1 3 ], x , 2 F 1 [ 1 9 , 7 9 ], [ 1 3 ], x . (8)

The main results

The globally bounded 3 F 2 hypergeometric series

3 F 2 [ 2 9 , 5 9 , 8 9 ], [ 2 3 , 1], 27 • x , 3 F 2 [ 1 9 , 4 9 , 7 9 ], [ 1 3 , 1], 27 • x (9)
are ‡ respectively the diagonals of the two algebraic functions

3 F 2 [ 2 9 , 5 9 , 8 9 ], [ 2 3 , 1], 27 • x = Diag (1 -x -y) 1/3 1 -x -y -z , ( 10
) † † Instead of 2 F 1 ([c, 1], [e],
x), or one could take any one of the three permuted versions:

2 F 1 ([b, 1], [e], x) , etc.
Appendix F p.58 of [START_REF] Bostan | Ising n-fold integrals as diagonals of rational functions and integrality of series expansions: integrality versus modularity Preprint[END_REF], also see [START_REF] Christol | Fonctions hypergéométriques et diagonales de fonctions rationnelles[END_REF] p.19. † See [START_REF] Rivoal | Roques Hadamard products of algebraic functions[END_REF] for a proof that 3 F 2 ([1/9, 4/9, 7/9], [2/3, 1], x) cannot be written as a Hadamard product. ¶ One can see this experimentally by taking the series expansion of any of the Gauss hypergeometric functions: the prime numbers in the denominators of the coefficients grow continuously. ‡ The operators annihilating the two hypergeometric functions (9) are adjoint of each other. and

3 F 2 [ 1 9 , 4 9 , 7 9 ], [ 1 3 , 1], 27 • x = Diag (1 -x -y) 2/3 1 -x -y -z . (11)
These two hypergeometric series † † ( 9) can be recast into series with integer coefficients

3 F 2 [ 2 9 , 5 9 , 8 9 ], [ 2 3 , 1], 3 6 • x = 1 + 120x + 47124x 2 + 23483460x 3 + • • • , (12)
and

3 F 2 [ 1 9 , 4 9 , 7 9 ], [ 1 3 , 1], 3 6 • x) = 1 + 84x + 32760x 2 + 16302000x 3 + • • • (13)
Now Denef and Lipshitz in [START_REF] Denef | Algebraic power series and diagonals[END_REF] show that any power series in Q [[x 1 , . . . , x n ]], algebraic over Q(x 1 , . . . , x n ), is the diagonal of a rational function in 2n variables, and they give an algorithm to build this rational function. This means that we can construct the rational functions, whose corresponding diagonals are the We recall the algorithm of Denef and Lipshitz and apply it to the algebraic function (1-x-y) 1/3 /(1-x-y -z) in the first subsection below, and then we give the rational function and a generalization of the result in the second subsection. Finally, we give a second proof of the general result using binomial sums.

From diagonals of algebraic functions to diagonals of rational functions: Denef and Lipshitz

We explain a method which, for a given algebraic power series in n variables, constructs a rational function in 2n variables whose diagonal equals the diagonal of the given algebraic series. Moreover, the partial diagonal of that 2n-variable rational function, with respect to the pairs of variables (x 1 , x n+1 ), . . . , (x n-1 , x 2n ), yields the original n-variable algebraic power series. The method is described in the paper by Denef and Lipshitz [START_REF] Denef | Algebraic power series and diagonals[END_REF] in the proof of their Theorem 6.2. As a running example we use the three-variable algebraic function

f (x, y, z) = (1 -x -y) 1/3 1 -x -y -z , (14)
whose multi-Taylor series expansion at 0 is actually a power series in the three variables x, y, z:

f (x, y, z) = 1 + 2 3 x + 2 3 y + z + 10 9 xy + 5 3 xz + 5 3 yz + 40 9 xyz + . . . ( 15
)
Note that the minimal polynomial of f is given by

p(x, y, z, f ) = (x + y + z -1) • f
However, by considering f = f -1, i.e. by removing the constant term of f , we can achieve an étale extension. The minimal polynomial then reads p(x, y, z, f ) = (x + y + z -1)

• (f + 1) 3 + 1 -x -y. (17)
Indeed, ∂ p ∂f (0, 0, 0, 0) = -3 = 0. According to the proof of Theorem 6.2 (i) in [START_REF] Denef | Algebraic power series and diagonals[END_REF], the rational function

r(x, y, z, f ) = f 2 • ∂ p ∂f (xf, yf, zf, f ) p(xf, yf, zf, f ) (18)
has the property that D r(x, y, z, f ) = f (x, y, z), and hence D r(x, y, z, f ) = f (x, y, z) for r(x, y, z, f ) = r(x, y, z, f ) + 1. Here the operator D denotes a special kind of "diagonalization" with respect to the last variable: for

f (x 1 , . . . , x n , y) = a i1,...,in,j • x i1 1 • • • x in n y j , (19)
one defines

D f (x 1 , . . . , x n , y) = j=i1+•••+in a i1,...,in,j • x i1 1 • • • x in n . (20)
In our running example we obtain:

r(x, y, z, f ) = 3 f 2 • (f + 1) 2 • (xf + yf + zf -1) 3 (f + 1) 3 • (xf + yf + zf -1) 3 -xf -yf + 1 + 1. (21)
In the second step, which is explained in the proof of Theorem 6.2(ii) of [START_REF] Denef | Algebraic power series and diagonals[END_REF], one has to transform the rational function r (that has n + 1 variables) into another rational function (having 2n variables) such that its "true" (partial) diagonal gives the nvariable algebraic series f . It consists of a sequence of n -1 elementary steps, each of which is adding one more variable. In our example, we have to do the following

r 1 (x, y, z, u 1 , v 1 ) = u 1 • r(x, y, z, u 1 ) -v 1 • r(x, y, z, v 1 ) u 1 -v 1 , (22)
r 2 (x, y, z, u 1 , u 2 , v 2 ) = u 2 • r 1 (x, y, z, u 1 , u 2 ) -v 2 • r 1 (x, y, z, u 1 , v 2 ) u 2 -v 2 ,
and obtain with r 2 the desired rational function in six variables.

Generalization of the previous result

By the algorithm of Denef and Lipshitz given in the previous section, it is possible to show that the algebraic function

(1 -x -y) a/b 1 -x -y -z , (23)
corresponds to the following rational function in six variables, by taking the diagonal with respect to (x, u), (y, v) and (z, w):

a • u 3 v • (1 -ux -uy -uz) • (1 + u) a-1 • (1 -ux -uy -uz) a-1 (1 + u) a • (1 -ux -uy -uz) a -(1 -ux -uy) b • (u -v) • (v -w) - a • v 4 • (1 -vx -vy -vz) • ((1 + v) • (1 -vx -vy -vz)) a-1 (1 + v) a • (1 -vx -vy -vz) a -(1 -vx -vy) b • (u -v)(v -w) (24)
- a • u 3 w • (1 -ux -uy -uz)((1 + u) • (1 -ux -uy -uz)) a-1 (1 + u) a • (1 -ux -uy -uz) a -(1 -ux -uy) b • (u -w) • (v -w) - aw 4 • (1 -wx -wy -wz) • (1 + w) a-1 • (1 -wx -wy -wz) a-1 (1 + w) a • (1 -wx -wy -wz) a -(1 -wx -wy) b • (u -w) • (v -w) + 1.
The diagonal of the rational function ( 24) is annihilated by the linear differential operator of order three:

b 3 x 2 (1 -27 x) • D 3 x + b 2 x ((27 a -135 b) • x -a + 3 b) • D 2 x (25) -b • ((9 a 2 -63 a b + 114 b 2 ) • x + a b -b 2 ) • D x + (a -3 b) • (a -2 b) • (a -b),
and can be expressed as the 3 F 2 hypergeometric function

3 F 2 [ 3 a -b 3 a , 2 a -b 3 a , a -b 3 a ], [ a -b a , 1], 27 • x . (26)
In particular, the two hypergeometric functions 3 F 2 ([2/9, 5/9, 8/9], [2/3, 1], 27• x) and 3 F 2 ([1/9, 4/9, 7/9], [1/3, 1], 27• x) appearing in [START_REF] Christol | Diagonales de fractions rationnelles et équations de Picard-Fuchs[END_REF], correspond respectively to the parameters (a, b) = (1, 3), and (a, b) = (2, 3) in the algebraic function [START_REF] Rivoal | Roques Hadamard products of algebraic functions[END_REF]. Other values of the parameters (a, b) are not necessarily unresolved examples of Christol's conjecture.

For example if we consider the parameter values a = 1 and b = 7 , we see that the diagonal of ( 24) is given by the globally bounded † series (27)

3 F 2 [ 2 7 , 13
being a globally bounded series, which means that it can be written as a diagonal using one of the procedures given in Section 2.3. We note that algebraic functions close to the algebraic functions appearing in [START_REF] Christol | Fonctions hypergéométriques bornées[END_REF] and [START_REF] Christol | Fonctions hypergéométriques et diagonales de fonctions rationnelles[END_REF], also give 3 F 2 or 4 F 3 hypergeometric functions as their diagonals that are unresolved examples to Christol's conjecture: 13 21 , 20 21 ], [

Diag (1 -x -2 y) 2/3 1 -x -y -z = 3 F 2 [ 1 9 , 4 9 , 7 9 ], [ 2 3 , 1], 27 • x , (30)
Diag (1 -x -2 y) 1/3 1 -x -y -z = 3 F 2 [ 2 9 , 5 9 , 8 9 ], [ 5 6 , 1], 27 • x , (31)
Diag (1 -x) 1/3 1 -x -y -z = 4 F 3 [ 2 9 , 5 9 , 8 9 , 1 2 ], [ 1 3 , 5 6 , 1], 27 • x , (32)
Diag (1 -x -y) 1/3 1 -x -z = 4 F 3 [ 2 9 , 5 9 , 8 9 , - 1 3
], [ 1 3 , 5 6 , 1], 27 • x . ( 33
) † 3 F 2 [ 2 7 ,

Proof

A computer algebra proof of this result can easily be obtained using the creative telescoping program [16]: one computes the operator (25) using the program [16], and verifies that this operator does annihilate the diagonal of (23) † †. Another longer way to do it which we provide below, is through binomial sums.

The denominator of the algebraic function (1 -x -y) a/b /(1 -x -y -z) can be expanded as a geometric series:

(1 -x -y -z) -1 = ∞ n=0 ∞ m=0 n m • (x + y) m z n-m = ∞ n=0 ∞ m=0 ∞ l=0 n m m l • x l y m-l z n-m , (34)
while the numerator can be written as the sum:

(1 -(x + y)) a/b = ∞ k=0 (-a/b) k k! • (x + y) k = ∞ k=0 k j=0 (-a/b) k k! • k j x j y k-j . ( 35
)
Multiplying these two sums (34) and ( 35) and re-indexing, we obtain:

∞ s=0 ∞ t=0 ∞ u=0 x s y t z u • s j=0 ∞ k=0 (-a/b) k k! • k j s + t + u -k s + t -k s + t -k s -j . (36)
Now taking the coefficients corresponding to the diagonal in (36), i.e. such that s = t = u = n, we get: n-j , we find that (37) can be written as

S(n) = 2n n • 2n k=0 (-a/b) k k! • 3n -k 2n -k , (38)
and by using a computer algebra tool like Mathematica or Maple to simplify this sum into a closed form, from which we can read off the hypergeometric function representation of the diagonal. More precisely, we used creative telescoping (Zeilberger's algorithm) to prove that (38) satisfies the first-order recurrence: (40) † † One also needs to note that initial conditions have to be compared.

F 2 (

 2 [a, b, c], [1, e], x) is globally bounded, with the parameters a, b, c, e ∈ Q \ Z, then there are six ways to write the 3 F 2 ([a, b, c], [1, e], x) function as the diagonal of a rational function. This corresponds to the six ways to write the 3 F 2 ([a, b, c], [1, e], x) as a Hadamard product of hypergeometric functions:

(a - 3 b

 3 -3 b n) • (a -2 b -3 b n) • (a -b -3 b n) • S(n) = b 2 • (n + 1) 2 • (a -b -b n) • S(n + 1). (39)Together with the initial condition S(0) = 1, we obtain the closed formS(n) = 3 3n • (b -a)/3b n • (2b -a)/3b n • (3b -a)/3b n (b -a)/b n • n! 2 .

  being series that are not globally bounded. Hence the hypergeometric series (27) cannot be easily written as a Hadamard product, as explained in Section 2.3.In contrast, for a = 3 and b = 4 the diagonal of (24) which is given by the globally bounded ‡ series (28)

				21	,		20 21		], [	6 7	, 1], 27 x = 1 +	260 49	x +	188190 2401	x 2 + • • •	(27)
	with the 2 F 1 series
	2 F 1 [	13 21	,		20 21	], [	6 7	], 27 x , 2 F 1 [	2 7	,	20 21	], [	6 7	], 27 x , 2 F 1 [	2 7	,	13 21	], [	6 7	], 27 x ,
	3 F 2 [	3 4	,	5 12	,	1 12	], [	1 4	, 1], 27 x = 1 +	45 16	x +	41769 1024	x 2 + • • •	(28)
	with the 2 F 1 series
												2 F 1 [	5 12	,	1 12	], [	1 4	], 27 x ,

Acknowledgments. J-M. M. and Y.A. would like to thank G. Christol for many enlightening discussions on diagonals of rational functions. We would like to thank A. Bostan for useful discussions on creative telescoping. Y.A. would like to thank A. Bostan for many explanations on Christol's conjecture including the details of 2.3 in an enlightening email correspondence [1]. J-M. M. would like to thank the School of Mathematics and Statistics of Melbourne University where part of this paper was written for hospitality. J-M. M. would like to thank A.J. Guttmann for many discussions on D-finite series. We would like to thank A.J. Guttmann for proofreading the paper. Y.A. would like to thank J. Voight for providing enlightening explanations on Shimura curves and for providing reference [20]. C.K. was supported by the Austrian Science Fund (FWF): P29467-N32. Y. A. was supported by the Austrian Science Fund (FWF): F5011-N15. Y. A. would like to thank the RICAM for hosting him on several occasions. Y.A. would like to thank Elaine and Rob for their hospitality in Linz, and for the great discussions he had with them there. Y.A. would like to thank his parents and his family for all their support and love. We thank the Research Institute for Symbolic Computation for giving us access to the RISC software packages.

Telescopers of algebraic functions versus diagonals of algebraic functions

The diagonal of a rational function and a solution of a telescoper † of a rational function are close, yet distinct notions. A telescoper annihilates an n-fold integral of a rational function over all integration cycles † †. For example the 3 F 2 ([a, b, c], [d, 1], x) can be written through the well-known integral representation as:

( (1 -y) 10/9 • (1 -x y 2 ) 1/9 • (1 -z) 7/9 , (42) admits as a solution the hypergeometric function 3 F 2 ([ 1 9 , 4 9 , 7 9 ], [ 1 3 , 1], x). Yet the diagonal of the algebraic function ( 42) is equal to zero, and it is through the algebraic function [START_REF] Christol | Fonctions hypergéométriques et diagonales de fonctions rationnelles[END_REF] that we were able to obtain it as the diagonal of an algebraic function. Other 3 F 2 unresolved examples to Christol's conjecture like [START_REF] Christol | Fonctions hypergéométriques bornées[END_REF][START_REF] Christol | Fonctions hypergéométriques et diagonales de fonctions rationnelles[END_REF]

were not obtained here as diagonals of a rational function, yet they are solution of the telescoper of an algebraic function and can thus be seen as a period of an algebraic variety over a non-evanescent cycle ‡, but not necessarily as a diagonal of an algebraic function (i.e. a period over an evanescent cycle). We give two arguments in favour of the fact that the 3 F 2 hypergeometric function ( 43) is most probably a diagonal of an algebraic function.

Diagonal: algebraic mod p

If one expects 3 F 2 hypergeometric functions like (43) to be diagonals of an algebraic function, one should find [START_REF] Bostan | Ising n-fold integrals as diagonals of rational functions and integrality of series expansions[END_REF][START_REF] Bostan | Ising n-fold integrals as diagonals of rational functions and integrality of series expansions: integrality versus modularity Preprint[END_REF] that the corresponding series expansion reduces to an algebraic series modulo any prime number p, or power of a prime number p r . In order to verify this fact on (43) we look at the series expansion of

, 1], 27

(45) † By "telescoper" of a rational function R(x, y, z) we denote the output of the creative telescoping program [16], applied to the transformed rational function R(x/y, y/z, z)/(yz), which is a differential operator that annihilates the diagonal of R. † † Diagonals correspond only to evanescent integration cycles over rational functions. ‡ To be totally rigorous, one has to consider the two certificates of the telescoping equation see if that the integral of the derivatives of these two certificates over that cycle are actually zero.

Straightforward guessing gives the infinite product formula

which is solution of

i.e. 3 F 2 [ 1 9 , 4 9 , 5 9 ], [ 1 3 , 1], 27 2 • x is an algebraic function modulo 2 satisfying:

or:

Modulo 3 we have the following expansion

where:

which is solution of

i.e. F (x) is an algebraic function modulo 3:

Unlike the situation on the 3 F 2 [ 1 9 , 4 9 , 7 9 ], [ 1 3 , 1], 27 2 • x hypergeometric series, it is less obvious how to obtain the 3 F 2 [ 1 9 , 4 9 , 5 9 ], [ 1 3 , 1], 27 2 • x as the diagonal of a rational function.

It is however possible to obtain the solution of

], [ 1 3 , 1], 27 2 • x , as a telescoper of an algebraic function, and this solution is an algebraic function modulo p. The diagonal of the product of algebraic functions

A relation between

is given by the 4 F 3 hypergeometric function H which is the Hadamard product of

This 4 F 3 hypergeometric series (55) can also be written as the Hadamard product:

So even though we did not find a rational (or algebraic) function whose diagonal is given by (43), knowing that 3 F 2 [ 1 9 , 4 9 , 7 9 ], [ 1 3 , 1], 27 • x is the diagonal of a rational function, we see that the Hadamard product of (43) with a simple algebraic function (1 -x) -7/9 is actually a diagonal of an algebraic (or rational) function. This suggests but does not prove, that 3 F 2 [ 1 9 , 4 9 , 5 9 ], [ 1 3 , 1], 27 • x could also be a diagonal of a rational function.

Conclusion

Because of the crucial role played by diagonals of rational functions in physics [START_REF] Bostan | Ising n-fold integrals as diagonals of rational functions and integrality of series expansions[END_REF][START_REF] Bostan | Ising n-fold integrals as diagonals of rational functions and integrality of series expansions: integrality versus modularity Preprint[END_REF], Christol's conjecture is an important open problem.

We have shown that the hypergeometric series 3 F 2 ([2/9, 5/9, 8/9], [2/3, 1], x) and 3 F 2 ([1/9, 4/9, 7/9], [1/3, 1], x) appearing in [START_REF] Bostan | Ising n-fold integrals as diagonals of rational functions and integrality of series expansions: integrality versus modularity Preprint[END_REF][START_REF] Christol | Fonctions hypergéométriques et diagonales de fonctions rationnelles[END_REF] are diagonals of rational functions. We did so by first finding two algebraic functions whose diagonals were given by these two hypergeometric functions, and through an algorithm outlined in the paper [START_REF] Denef | Algebraic power series and diagonals[END_REF], we were able to recover the rational functions whose diagonals are given by these two 3 F 2 hypergeometric functions.

We were also able to give a generalization of this result, and obtain other unresolved examples of Christol's conjecture as diagonals of rational functions. Furthermore, even though we were not able to write the 3 F 2 ([1/9, 4/9, 5/9], [1/3, 1], 27 • x) given by Christol in [START_REF] Christol | Fonctions hypergéométriques bornées[END_REF], as a diagonal of a rational function, we gave two arguments that suggested that it was likely to be so. More generally, we believe after writing the 3 F 2 ([2/9, 5/9, 8/9], [2/3, 1], x) and 3 F 2 ([1/9, 4/9, 7/9], [1/3, 1], x) as diagonal of rational functions, that it is likely that the other 3 F 2 unresolved examples of Christol's conjecture are diagonals of rational functions.

Appendix A. Counterexamples and links with Shimura curves

The Gauss hypergeometric function appearing on the left in (9)

, 1], 27 x (A.1) can be seen as the Hadamard product of a Gauss hypergeometric function and an algebraic function given by: [START_REF] Bostan | Ising n-fold integrals as diagonals of rational functions and integrality of series expansions[END_REF] ], 27 x happens to be a hypergeometric function corresponding to an automorphic form associated with a Shimura curve [START_REF] Abdelaziz | Heun functions and diagonals of rational functions[END_REF][START_REF] Voight | Shimura curves of genus at most two[END_REF][START_REF] Voight | Three lectures on Shimura curves[END_REF]. One has the identity:

], -1728

(A.

3)

The Gauss hypergeometric function 2 F 1 [ 1 36 , 19 36 ], [ 8 9 ], x can be also expressed as:

Now the Gauss hypergeometric function 2 F 1 ([ 1 36 , 13 36 ], [ 8 9 ], x) which occurs in p.14 of [START_REF] Van Hoeij | Belyi functions for hyperbolic Hypergeometric-to-Heun transformations[END_REF], corresponds to a hypergeometric function related to a Shimura curve since it has exponent differences † (1/9, 1/2, 1/3), and these exponent differences are listed in the exhaustive list of hypergeometric functions that are associated with Shimura curves appearing in Table 1 of [START_REF] Takeuchi | Commensurability classes of arithmetic triangle groups[END_REF]. Other 3 F 2 functions that are unresolved examples to Christol's conjecture that we found to be related to 2 F 1 hypergeometric functions related to Shimura curves are given by:

], [

], 3 6 x , (A.5)

], 3 6 x , (A.6)

], [

, 1], 3

], 3 3 x , (A.7) ], [ 4 3 , 1], 27 x), and the three globally bounded 3 F 2 hypergeometric series (A.6), (A.7) and (A.8), we were not able to write the other examples given in this section as a Hadamard product involving a 2 F 1 hypergeometric function associated to a Shimura curve. In any case, since the class of potential counterexamples formulated by Christol is infinite, while the list of Shimura in Table 1 of [START_REF] Takeuchi | Commensurability classes of arithmetic triangle groups[END_REF] is finite, a list of 3 F 2 functions both related to Shimura curves and to Christol's conjecture is bound to be finite. † See [START_REF] Tu | Algebraic Transformations of hypergeometric functions arising from of Shimura curves RIMS Kôkyûroku Bessatsu[END_REF] p.10 for a definition of exponent difference.