Dan Xu

Etienne Cantin

Florian Frank

Nicolas Quintin

Frédéric Meynadier

Philip Tuckey

Anne Amy-Klein

Olivier Lopez

Paul-Eric Pottie
email: paul-eric.pottie@obspm.fr

Two-branch fiber link for international clock networks

We present our work on realizing two-branch fiber links enabling multiple-partner clock comparisons in Europe. We report in detail on the setup connecting two long-haul links, and report for the relative frequency stability and accuracy. We report on a combined uptime of 90% for our dual-branch link during almost one month. We finally discuss the combined uncertainty contribution of the ensemble of the link architecture to a clock comparison. We show that the frequency transfer uncertainty is 2 × 10 -19 .

I. INTRODUCTION

Comparisons of clocks through fiber links at the continental scale have become a reality in Europe [START_REF] Lisdat | A clock network for geodesy and fundamental science[END_REF], [START_REF] Guéna | First international comparison of fountain primary frequency standards via a long distance optical fiber link[END_REF], [START_REF] Delva | Test of special relativity using a fiber network of optical clocks[END_REF]. One challenge is now to combine several coherent fiber links to build a clock network, with the means of comparison making a negligible contribution to the uncertainty compared to the resolution of the clocks. These prospects rely on the ability to realize a multi-branch fiber network. A second challenge is to make the fiber links highly available for on demand optical clocks comparison and for long-period of comparison of fountain cold atoms clocks, over campaign duration of several weeks and months.

For one aiming at effective optical frequency dissemination service to multi-users, the ability to build optical network and a high signal availability are the key features that are needed. Indeed the comparison arises from the multiplication of the uptime of all the elements constituting the metrology chain. For instance, a two-partite comparison based on 2 clocks, 2 combs and 2 links with uptime of 50% for each of them leads to a total uptime of only 1.6%. This means that 10 days of operation are needed to obtain 4 hours of non-continuous data. The more partners join the network, the lower the uptime will be, and it makes the use of the dissemination network almost impractical. This calculation illustrates the strong motivation and need to achieve optical link uptime as close as possible to 100%. This requires a high level of automation, remote control, and technical readiness level, but also a delicate control of many experimental aspects when connecting two links together. With that respect, the work presented here is a key step towards a new definition of the SI second based on optical frequency standards and towards novel opportunities for stringent tests of fundamental physics [START_REF] Riehle | Towards a redefinition of the second based on optical atomic clocks[END_REF], [START_REF] Hees | Searching for an Oscillating Massive Scalar Field as a Dark Matter Candidate Using Atomic Hyperfine Frequency Comparisons[END_REF].

Another issue that arises with long-duration comparison campaign is the limitation due to the accuracy and the stability of the frequency and time references available for counting the optical beat notes. Fortunately the best performances are not needed everywhere, but at least some level of performances are needed not to loose the benefit of the optical frequency dissemination. Nowadays commercial products based on ultrastable radio-frequency (RF) oscillator disciplined on Global Positioning System (GPS) are available, and some of them may comply with the needed accuracy in time and frequency. However Global Navigation Satellite System (GNSS) based solutions are not always available. So alternative solutions need to be considered.

In this paper, we report on the first implementation of a twobranch fiber link, as part of a wider concatenation of links between National Metrological Institutes (NMIs) in Europe enabling multiple clock comparison [START_REF] Often | Optical Frequency Transfer in a European Network[END_REF]. We discuss then the results and the limitations due to the RF stability and accuracy.

II. EXPERIMENTAL SET-UP

The sketch of the experiment is depicted on Fig. 1. We start our experiment with an ultra-stable laser operated at 194.4 THz, continuous wave [START_REF] Argence | Prototype of an ultra-stable optical cavity for space applications[END_REF]. A slave laser is phase-locked to the ultra-stable laser, and actively de-drifted using an optical frequency comb to measure the drift versus an H-Maser [START_REF] Targat | Experimental realization of an optical second with strontium lattice clocks[END_REF]. The light is then injected into an actively stabilized singlemode fiber of 30-m long and transferred to another experimental room of our laboratory where the links are operated [START_REF] Ma | Delivering the same optical frequency at two places: accurate cancellation of phase noise introduced by an optical fiber or other time-varying path[END_REF]. Then, we have two repeater laser stations, that are the starting point for two long-haul fiber links [START_REF] Chiodo | Cascaded optical fiber link using the internet network for remote clocks comparison[END_REF], [START_REF] Xu | Studying the fundamental limit of optical fiber links to the 10 -21 level[END_REF]. Each of these links make use of two fibers, one to transfer the optical frequency from our laboratory to the remote point, and the second fiber to evaluate the performance of the frequency transfer.

The upper link on the right of Fig. 1 is a 4-span cascaded link of 2×705 km from SYRTE to University of Strasbourg and back [START_REF] Lisdat | A clock network for geodesy and fundamental science[END_REF], [START_REF] Chiodo | Cascaded optical fiber link using the internet network for remote clocks comparison[END_REF]. The fiber pair is used for parallel data traffic using the French academic network of RENATER. Therefore the link is equipped with 40 optical add-drop multiplexers (OADMs). 16 bi-directional Erbium-doped fiber amplifiers (EDFAs) compensate partially for the losses, that are about 410 dB, i.e., about 0.29 dB/km. At University of Strasbourg, the signal transferred with this link is compared with the output signal from another long-haul link from PTB in Braunschweig [START_REF] Lisdat | A clock network for geodesy and fundamental science[END_REF]. For that purpose, repeater laser stations are operated at the remote end of these two long haul links and the optical frequency comparisons between the international links are done using two-way method between the userend outputs [START_REF] Cantin | Progress on the REFIMEVE+ Project for Optical Frequency Standard Dissemination[END_REF], which is useful for the clock comparison campaign [START_REF] Lisdat | A clock network for geodesy and fundamental science[END_REF], [START_REF] Guéna | First international comparison of fountain primary frequency standards via a long distance optical fiber link[END_REF], [START_REF] Delva | Test of special relativity using a fiber network of optical clocks[END_REF].

The lower link on the right of Fig. 1 uses a pair of 43-km dedicated fibers from SYRTE to LPL in Villetaneuse. It is a hybrid fiber link, combining active compensation on one fiber and passive compensation on the second fiber, using a two-way beatnote between the local laser and the laser transferred at LPL. Many details of the implementation are reported in [START_REF] Xu | Studying the fundamental limit of optical fiber links to the 10 -21 level[END_REF]. The link is equipped with four OADMs, so that other time and frequency services can be operated in parallel. We are using three bi-directional EDFAs. The losses are about 32 dB in total, i.e., about 0.37 dB/km. At LPL, the transferred frequency signal is compared with the signal transferred from NPL through a 769-km noise-compensated optical link, using two-way frequency comparison [START_REF] Delva | Test of special relativity using a fiber network of optical clocks[END_REF].

For both upper and lower links on Fig. 1, the residual instability is estimated from the end-to-end measurements which are made with the integrated multi-branch Michelson interferometer available in the repeater laser stations [START_REF] Chiodo | Cascaded optical fiber link using the internet network for remote clocks comparison[END_REF], [START_REF] Lopez | Cascaded multiplexed optical link on a telecommunication network for frequency dissemination[END_REF]. With this first generation of interferometers built couple of years ago, the length difference between the references arms were not well minimized [START_REF] Stefani | Tackling the limits of optical fiber links[END_REF]. The length imbalance of the reference arms of the Michelson interferometer of the second repeater laser station is about 15 cm [START_REF] Xu | Studying the fundamental limit of optical fiber links to the 10 -21 level[END_REF]. The remote stations are engineered in such a way that the unaccuracy or instability of their RF reference is not degrading the transferred signal frequency and the end-to-end signal frequency [START_REF] Lopez | Cascaded multiplexed optical link on a telecommunication network for frequency dissemination[END_REF].

For the specific purpose of a multipartite clock comparison, the frequency of the infra-red (IR) laser sent to the two links has to be precisely chosen so that the beat note frequencies recorded with the remote IR lasers are within the available bandpass filters and the capabilities of the frequency counters. Once the frequencies of the offset phase locks and acoustooptic modulators (AOMs) depicted at the upper part of the sketch Fig. 1 are chosen, the frequency offset of the lower repeater laser station can be used to tune the beat frequency with the remote IR lasers at the far end of the second link. The issue is to control the optical phase between the two repeater laser stations, that is the frequency fluctuations induced in the short 5-m link between these stations.

We implemented the preliminary solution depicted in Fig. 1 for compensating these frequency fluctuations. An end-user output of the upper repeater laser station is used to feed the input port of the lower station. The offset frequency of the second repeater laser station is shifted to a well chosen frequency by an AOM. Its RF frequency is 37 MHz. We use the Faraday mirror at the input of the second repeater laser station to reflect back the laser light to the photodiode at the output of the first repeater laser station. On this photodiode, we measure the round-trip delay fluctuations at the carrier RF frequency of 74 MHz. This signal is tracked and divided by 74, amplified and counted with a dead-time free phase-frequency counter (K-K Messtechnik FXE) with 1-s gate time. It is then subtracted from the clock comparison data performed at the remote end of the link. This solution is a variant of a two-way technique, with one laser light, one AOM and one Faraday mirror. This experimental arrangement is very easy and fast to implement and enables us to connect the two stations without using any servo loops.

III. EXPERIMENTAL RESULTS

We now present the data acquired within the frame of the 2017 campaign of comparisons between NPL, SYRTE and PTB, for which we used only one laser source as a seed for two links. We focus on the 5-m short link connecting the two stations. The acquired data represent the free-running round-trip frequency noise of the link to be compensated in post-processing. Figure 2 shows the histogram for 12 days of continuous data set and a slight asymmetry is observed. One can denote that the tail of the distribution is longer towards smaller frequencies. For a more quantitative analysis, the mean value is 7.9 × 10 -19 and the median is 1.7 × 10 -17 . By calculating the rolling mean with a time window of 3600 s, we found that the mean of the relative frequency fluctuation over 3600 s samples varies typically from ±2 × 10 -17 with a daily period, and at worse from -4 × 10 -17 to +6 × 10 -17 . These fluctuations are two orders of magnitude higher than the residual noise of the long-haul link. These free running fluctuations are fortunately removed in postprocessing.

We also present the Allan deviation and Modified Allan deviation of the same 12-days long continuous data sample on Fig. 3. We observe that both deviations start at about 2.3 × 10 -16 , and reach a maximum of 4 × 10 -16 at 10 s integration time. A clear signature of periodic daily perturbation is to be seen with a local maximum of 8 × 10 -18 at 20 000 s integration time. The periodic perturbation is due to the temperature of the room that fluctuates daily with a typical amplitude of 1.6 K during this campaign. The signature seen on the Allan deviation is in agreement with the expectation, knowing the temperature fluctuation and the length of the fiber patch-cord. The residual instability after post-processing can be estimated using the link model derived in [START_REF] Newbury | Coherent transfer of an optical carrier over 251km[END_REF], under the assumption that the fiber noise is reciprocal. It results in a post-processed stability of a few 10 -23 at 1 s integration time and around 10 -30 at long-term, which is well below the noise floor of our system. Thus we expect to reach the noise floor of around 10 -17 at 1 s and 10 -21 at long-term. An upper value of this residual noise could also be evaluated from the residual noise of the international frequency comparisons at LPL or University of Strasbourg after post-processing. Uncompensated fiber paths between the Faraday mirror and the coupler in the short 5-m link are also limiting the stability, since this part is only roughly protected against temperature fluctuations. Considering an uncompensated length of 1 m, and estimating residual daily temperature variations of 0.5 K or below, we estimate this limitation to be 4.3 × 10 -19 at half a day averaging time and around 10 -19 at 1 day.

The end-to-end data and their analysis for the two long haul links can be found elsewhere [START_REF] Chiodo | Cascaded optical fiber link using the internet network for remote clocks comparison[END_REF], [START_REF] Xu | Studying the fundamental limit of optical fiber links to the 10 -21 level[END_REF]. For the four-spans link to Strasbourg, the modified Allan deviation is below 10 -15 at 1 s integration time and the long-term stability is typically 10 -19 or below. For the link to LPL, the modified Allan deviation is 6 × 10 -17 at 1 s integration time and the long-term stability is typically 10 -20 at one day integration time. According to [START_REF] Chiodo | Cascaded optical fiber link using the internet network for remote clocks comparison[END_REF], [START_REF] Xu | Studying the fundamental limit of optical fiber links to the 10 -21 level[END_REF], the frequency bias introduced by the cascaded links and hybrid links are respectively -4.8×10 -20 with an uncertainty of 9×10 -20 and 4.2×10 -21 with an uncertainty of 8×10 -22 . The combined uncertainty is therefore conservatively set to 1 × 10 -19 . Finally the combined uncertainty of the short and long links and the repeater stations is conservatively set to 2 × 10 -19 .

We now discuss the uptime of the multi-branch and composite fiber link, for a duration of 27 continuous days. The various servo loops of the links can indeed unlock or experience some cycle slips due to some short excess noise or during the polarisation optimization in the repeater stations [START_REF] Chiodo | Cascaded optical fiber link using the internet network for remote clocks comparison[END_REF]. The data are summarized in Table I. The uptime of the 30-m active compensation link connecting the two rooms in our laboratory is 100%. The uptime of the post-processing for the 5-m short link connecting the two stations is assessed by applying a data filtering with a bandwidth of 1 Hz, i.e., about 10 times the standard deviation for integration time ranging from 1 s to 10 s. Since it was installed in a quite noisy environment and not enough protected, it could experience some excess environmental noise which unlocks the tracking used to count the beatnote, leading to an uptime of 98.5%. This shows the need for an operational room with minimum human intervention.

For the two long-haul links, we assess the uptime from the end-to-end measurements. The hybrid fiber link is assessed from a simple data selection within a 1 Hz bandwidth and removal of outliers for one of the two-way beatnote signals recorded at SYRTE (the TWB3 observable, see details reported in [START_REF] Xu | Studying the fundamental limit of optical fiber links to the 10 -21 level[END_REF]). In the case of our hybrid setup, this is a robust strategy as the short-term noise is about 4 × 10 -17 , compared to the hop in frequency of 5×10 -15 induced by a single cycle slip. Thus each single cycle slip can be well removed with this procedure and we end up with a conservative estimate of 95%.

Regarding the four-span cascaded link, we developed a more subtle strategy for selecting the data to be kept for a clock comparison. We follow a three-step procedure to select the data. We first apply a very large data filtering with a bandwidth of 10 Hz and produce a boolean quality factor that is set to zero if a point has been removed. Then we process the data and look at the time evolution of three observables, as displayed on Fig. 4. The rolling mean of the relative frequency fluctuation y is performed on a short time of 9 s and allows for detection of outliers. The rolling standard deviation of y is calculated with a longer period of 2750 s in order to detect anomalous noise and reject the period during which the link is not optimally working, and finally the rolling standard deviation of the previously defined quality factor allows for rejection of segment of data where the link is not stable and there is a lot of outliers. This method of data selection was introduced for the purpose of clock comparison and applied on 6 months of data. It seems to be robust as similar performances of the link is found without specific adjustment of the selection procedure parameters. An example of the time evolution of these three observables are displayed on Fig. 4 and shows outliers and period of anomalous noise, which are filtered from the data set.

IV. DISCUSSION

We have shown an easy way to link two repeater laser stations together, with a variant of two-way method, that allows us to match technical requirements of a specific clock comparison campaign. The passive listening of the frequency noise on the short link that connects the two stations reveals asymmetry. These fluctuations are removed from the clock comparison in post processing. It shows that, at the level of exquisite frequency control that optical frequency combs and links are able to achieve, every single detail matters to keep the same level of accuracy. In the specific case of achieving star-liked networks, more complex multi-branches Michelson interferometers were achieved in our groups for the REFIMEVE+ (REseau FIbré MEtrologique à Vocation Européenne+) network, so-called 'super-station', or multibranches eavesdropping setups as reported in [START_REF] Bercy | Ultra-stable optical frequency dissemination on a multi-access fibre network[END_REF]. These are clearly better and more robust solutions. However we observe that, if the short link between the two stations would not been monitored, it would be the dominant source of noise and inaccuracy of our set-up. So we believe that this work shines light on the high level of preparation that is needed at the enduser level for a full benefit of the quality of service of optical frequency dissemination network such as REFIMEVE+. The combined uptime of the two-branch link is 90%. If the uptime may be compared, this is quite interesting to note that the 1410 km link is now extremely stable and robust, thanks to our efforts for automation and remote monitoring abilities. The uptime of the end-to-end is now probably limited by the uptime of servo loops and set-ups in the clocks and combs in our laboratory that is not fully equipped with automated loops and supervision tools. As compared to previous clock comparisons campaign, for which the uptime of the Paris-Strasbourg link was roughly 90% but the uptime of the whole Paris-Strasbourg-Pars link was less than 25%, the uptime of our link was considerably improved, and it is comparable to the uptime of atomic fountains as operated at SYRTE and PTB. So now, the strongest limit on the combined uptime of an optical clock comparison is likely to be the uptime of the neutral atom optical clocks themselves.

Further study on the procedure to validate data need to be carried out. A more precise analysis is especially needed to assess carefully the impact on frequency stability and accuracy of the holes made in the raw data set. With the procedure described in this paper, we report a combined accuracy level for the two long-haul links of 1×10 -19 . This can be compared to the level of accuracy available at remote end points, as in Strasbourg for instance.

Given that the frequency of the beat note can be up to 55 MHz with the currently used K+K FXE counters, and that the long term stability of commercially available solutions as provided by GPS receiver as the Oscilloquartz OSA 5201 are about 3 × 10 -13 after one day integration time as recorded versus UTC(OP), one obtain a similar contribution to the uncertainty budget of 8.3 × 10 -20 . As the performances of the optical clock will improve, that the comparison campaign will last longer, and due to the lack of reliability of some of these commercial GPS receiver, timing solutions provided inband through the fiber will be needed at some crucial point of the network. White-rabbit solution for instance, considering its excellent level of performance at the level of 1 × 10 -14 at 100 s integration time after 500 km fiber link, and considering its excellent performance-to-cost ratio, seems to us a very promising solution [START_REF] Kaur | Time and frequency transfer over a 500 km cascaded White Rabbit network[END_REF].

By contrast the requirement on time synchronisation of the measurements over such an optical frequency dissemination network could relax with the time. Basically the inaccuracy imprint on the measurement scales as the ultra-stable laser drift times the de-synchronisation of the measurement, and can be rapidly a dominating source [START_REF] Lisdat | A clock network for geodesy and fundamental science[END_REF]. It is quite interesting to notice that the drift of the laser improved by several orders of magnitude over the last five years, from typically a few Hz/s down to sub-mHz/s for the state-of-the-art silicon cavity [START_REF] Zhang | Ultrastable Silicon Cavity in a Continuously Operating Closed-Cycle Cryostat at 4 K[END_REF]. With active de-drift system, 1∼2 mHz/s are achievable nowadays, and it makes the need of accuracy for the synchronisation of the measurement up to the 50 ms level for the same level of accuracy of 1 × 10 -19 .

V. CONCLUSION

We have shown a two-branch long-haul fiber link that enables the comparison of multiple optical clocks at their best level of performances to date. We have shown that the two-branch fiber link is not limiting the uptime of the clock comparison. This work represents a decisive step towards the realization of a robust fiber network for clock comparisons, contributing to progress along the roadmap towards the redefinition of the SI second. It also impacts the future use of the REFIMEVE+ network, since the benefit of ultra-stable signal also need good experimental control at the end-user side. We have also presented discussions on the requirements for RF reference signal and time reference signal, and show that the most stringent need are for the RF accuracy. Solutions based on White-Rabbit and/or other in-band techniques, far from being a burden, may open new opportunities for time and RF standard dissemination.

Fig. 1 .

 1 Fig. 1. Sketch of the experiment. Upper right link is a four-span link with the first span labelled as N and the 4th span labelled as N + 3. The lower link on the right is a two-span links with the two spans labelled as M and M + 1. OC: optical couplers. AOM: Acousto-Optic Modulator. FM: Faraday Miror. PLL: Phase Lock Loop. PD: Photodiode. The light red boxes are thermally controlled boxes.

Fig. 2 .

 2 Fig. 2. Histogram of the frequency fluctuations df occurring in the 5-m long fiber link, after 12 days of operation, for 1 s averaging time, in Hz. A slight asymmetry is observed. The red line is the mean m = 7.9 × 10 -19 . The green light is the median n = 1.7 × 10 -17 . Bins width = 45 mHz.

Fig. 3 .

 3 Fig. 3. Frequency instability of the 5 m link connecting the two repeater laser stations. Allan standard deviation and modified Allan deviation of the relative frequency fluctuations recorded with 1 s gate reporting in Λ-mode.

Fig. 4 .

 4 Fig. 4. Time evolution of the three observables used for selecting data in the 4-span cascaded link (a) rolling mean of the frequency fluctuations (b) rolling standard deviation of the frequency fluctuations and (c) rolling standard deviation of the quality factor (see text). Red lines indicates the filtering limits. Anomalous noise and outliers are clearly seen.

TABLE I UPTIME

 I OF ALL THE ELEMENTS OF THE TWO-BRANCH FIBER LINK .

		stabilized link passive	hybrid	cascaded
	length (m)	30	5	43 000 1 410 000
	uptime (%)	100	98.5	95	96.3

ACKNOWLEDGMENT

We thank Rodolphe Le Targat and Hector Alvarez-Martinez for the operation of the ultra-stable laser and the comb at start of our experiment. We thank Michel Abgrall and Baptiste Chupin for allowing us to measure GPS receiver versus UTC(OP). We would like to thank E. Camisard, T. Bono and L. Gydé for their support and for facilitating the access to the RENATER network and facilities, P. Gris and B. Moya of Uni. Strasbourg. We acknowledge funding support from the Agence Nationale de la Recherche (Labex First-TF ANR-10-LABX-48-01, Equipex REFIMEVE ANR-11-EQPX-0039, Idex PSL ANR-10-IDEX-0001-02), the European Metrology Programme for Innovation and Research (EMPIR) under SIB-02 OFTEN. EMPIR is jointly funded by the EMPIR participating countries within EURAMET and the European Union.


M Siddorn

F.-X Coudert

K E Evans

A Marmier

A systematic typology for negative Poisson's ratio materials and the prediction of complete auxeticity in pure silica zeolite JST †

Single crystals can commonly have negative Poisson's ratio in a few directions; however more generalised auxeticity is rarer. We propose a typology to distinguish auxetic materials. We characterise numerous single crystals and demonstrate that partial auxeticity occurs for around 37%. We find average auxeticity to be limited to α-cristobalite and no example of complete auxeticity. We simulate two hundreds pure silica zeolites with empirical potentials and quantum chemistry methods, and for the first time identify complete auxeticity in a zeolite network, JST.

Introduction

The main aims of this study are to develop a convenient typology of auxetic behaviour in materials, to characterise the Poisson's ratio of pure silica zeolites and to identify structures with exceptional values. The Poisson's ratio for anisotropic materials is complex, a function of three variables, two defining a longitudinal direction, one a transverse one. The adjective "auxetic", describing the existence of a negative Poisson's ratio (NPR) 1 , is too limited to describe fundamentally different situations, from single crystals where NPR occurs for very specific directions, to isotropic foams where NPR is present for all directions. Therefore an important objective has been to develop a finergrained typology of auxetic properties to discriminate between the relatively common existence of negative Poisson's ratios in a few narrowly defined combinations of longitudinal and transverse directions and the rarer, more comprehensive cases.

Background

Of the four elastic constants used to describe isotropic materials, Young's modulus (E), bulk modulus (K), shear modulus (G) and Poisson's ratio (ν), it is the Poisson's ratio that has historically been the least explored 2,3 . It can be associated with some interesting and unusual properties, particularly when in a range not normally encountered. Defined as the ratio of transverse to longitudinal strain in a structure or material, it has been accepted theory that the Poisson's ratio can have negative values for over 150 years but it is only since 1989 that it has been actively studied. Due to their unusual nature, auxetic materials could be used in many applications where the benefits of their properties have only recently become apparent. As Evans showed, highly negative ν (-1/2 > ν) can produce large values for indentation resistance and fracture toughness 2 . Of recent particular interest is the manufacture of blast curtains for defence purposes. When an object hits an auxetic fabric the surrounding material closes in on the impacted area rather than away from it, causing an increase in local density and potentially stopping fragments of shrapnel or projectiles 4 . The motivation to study NPR in the family of silicate materials specifically is two-fold. Firstly, zeolites in particular have historically received considerable interest due to their very low density and potential use as catalysts 5 or molecular sieves 6,7 . Early studies of their mechanical properties have been far from systematic and in light of recent results 8,9 it is timely to revisit auxeticity in silicates in general, and zeolites in particular, to compare their mechanical properties with those of other crystalline materials. Secondly, it has been recognised that, singularly among single crystals, α-cristobalite is auxetic in a large directions range 10 ; a logical place to search for materials with equally exceptional properties is amongst other silicates, especially largely unexplored pure silica zeolites. From Brillouin scattering on single crystals, the MFI zeolite was found to be auxetic in the (001) plane 11 . Auxetic zeolites are of interest due to the prospect that 'tuneable' molecular sieves could be created, for which the porosity could be altered by means of an applied stress 12 . Moreover, adsorption of specific guests in their nanopores could provide a way to tune their mechanical properties 13 . The negative Poisson's ratios observed in MFI and other zeolites has been explained using two dimensional mechanisms of rotating rigid bodies, and later more accurately using rotating semi-rigid bodies 14,15 . Pure

Typology of auxetic materials

Previous authors have realised that auxeticity can take different forms 16,17 and introduced the concepts of partial auxeticity and complete auxeticity for cubic crystals. The scheme we propose is more detailed and not restricted to cubic crystals. A material which exhibits auxetic behaviour may experience a negative Poisson's ratio in a very narrow range of transverse directions, for a narrow range of longitudinal directions of deformation. Conversely, it could be isotropic and exhibit negative Poisson's ratios for all transverse directions, for every longitudinal deformation; or it could be anywhere in-between. In order to more accurately discuss auxetic materials and their nature it is therefore desirable to establish a typology in which all categories of auxetic materials can be identified, described and distinguished. The following section describes such a general typology of the auxetic behaviour of a material.

Basic principles

In this context, we use the term 'direction' to refer to the transverse direction for which Poisson's ratio is observed and 'axis of deformation' to refer to the longitudinal axis of loading. Our classification is based on a type, a class and a corresponding numerical indicator; Table 1 describes the 10 possible combinations, from 0 (non auxetic), to 3C (auxetic for all axes and directions). 'Type' describes the axes of deformation, it indicates whether there is NPR around none (type 0), at least one (type 1), an average of (type 2), or every axis of deformation (type 3). 'Class' refines that information by considering NPR in the transverse directions, for at least one (class A), for an average (class B) or for every direction (class C). This is moderately abstract, and to illustrate the criteria which must be satisfied for a material to fall into the respective categories, Fig. 1 shows the elastic transverse deformations around specific axes in four examples (in other words, the axis of deformation is perpendicular to the page). The plots (b), (c), and (d) in Fig. 1 give examples of materials which meet the minimum requirement to be classified as auxetic type 1A, 1B and 1C respectively. While it is enough for a material to respect one of these three criteria along one axis of deformation only to be auxetic type 1, in order to be classified as type 2 the property must be present along an average of all axes of deformation (the way in which we decide if this average condition is met is a matter of choice and is discussed in more detail later in this section). To be classified as type 3 a material must meet the class criteria, shown in Fig. 1, for every axis of deformation. Such extreme auxeticity is rare in anisotropic materials but may be found in isotropic foams. It should be pointed out that a material that falls into a particular category may also fall into one or more others. For example, it is entirely possible for a material to be both 2B and 3A simultaneously. At the simpler level, lower classifications (that is to say less auxetic) are explicitly implied by all higher classifications. Both the type (1, 2 or 3), and the class (A, B or C), carry these implications, provided the alternate value does not increase. Consequently, if a material is 2C it follows it must also be 1C, a classification with a lower type, and additionally 2B and 2A, classifications with lower classes. A 2C material does not necessarily imply 3A, as although the class has been reduced, the type has been increase; a 3A classification could still be satisfied, depending on the material. When establishing the typology of an anisotropic material, the Poisson's ratio must be calculated for several directions around many axes. The Poisson's ratio is readily given by the opposite of the ratio of the two rotated elastic compliances, ܵ ଵଶ ᇱ and ܵ ଵଵ ᇱ .

The rotated compliances can themselves be derived from a set of elastic constants in the principle axes, be they in the form of the compliances matrix ܵ, or more commonly the stiffness matrix ܥ (for more detail on the notations of elasticity,

Physical Chemistry Chemical Physics Accepted Manuscript

This journal is © The Royal Society of Chemistry 2015 Phys. Chem. Chem. Phys., 2015, 00, 1-3 | 3

including relationship between order 4 tensors and order 2 matrices, see for instance [START_REF] Nye | Physical properties of crystals[END_REF] or [START_REF] Marmier | [END_REF] ). The elastic tensor/matrices have either been measured experimentally by methods such as Brillouin scattering or acoustic microscopy, or been calculated from the second derivatives of potential energy models, usually numerically, following an optimisation stage.

Calculating the Poisson's ratio for off-axis combination of longitudinal axis and transverse direction is not especially complex, but often tedious and repetitive. The method has been systematically implemented in the ElAM software package [START_REF] Marmier | [END_REF] . For a typical calculation, the three angles (θ, φ, χ) describing the axis vector and direction vector each varies every degree (see Fig. 2). This generates an excess of 23 million sets of off-axis elastic constants, from which the extrema and averages necessary to establish the auxetic typology can be determined.

Additional considerations on averages

Class A and class C can be determined unambiguously by the minimum or maximum ν around a single axis (if maximum ν is negative then class C, if minimum ν is negative then class A). Class B relies on ν being averaged, and different options for this are possible. The average could be considered as a direct mean (if the mean ν is negative then class B is satisfied), or a simpler median average (if more than half of the ν values are negative then class B is satisfied). When many single crystals are arranged in a random orientation, it is the average of the Poisson's ratio values that will affect the polycrystalline properties 20 rather than the total number of direction exhibiting a particular Poisson's ratio in each single crystal. Therefore, in this paper the mean Poisson's ratio is used to ascertain if a material satisfies class B, as this will give a greater idea as to which are likely to exhibit auxetic behaviour in a polycrystalline structure.

The type section of a classification is derived from the number of axes which exhibit the three different classes. The more generally a class is exhibited, the higher the type of this class will be. Type 1, the simplest case, is satisfied if there is any behaviour of a class exhibited by the material. Therefore, if at least one axis has a particular class then type 1 of this class will be satisfied.

In order for a material to satisfy type 2 classifications, what is considered as the 'average axis' must meet the criteria of a class.

As the average needed is of classes, rather than simple numerical values, the decision process is more subtle than that for determining class B auxeticity and it is valid to consider two distinct criteria. The first of these is whether or not a randomly chosen axis is expected to be of a certain class. This can be expressed as the median of the class when considering all axes, where the classification will be satisfied if more than half of the axes examined exhibit a class. The second criteria considered deals with numerical values which directly relate to the satisfaction of the three classes. Classes A, B, and C are satisfied when the minimum, average, and maximum Poisson's ratios are negative respectively. A material can be considered to be type 2 if the average of these numerical values for all axes is also negative. With this in mind, the classification for type 2 is split into two separate variants: 2i and 2ii. Type 2i is satisfied if more than half of the axes examined show behaviour of a particular class. Type 2ii is satisfied if the mean of relevant values to the class are negative, so that a material is 2iiA if the mean of the minimum ν is negative, 2iiB if the mean of the mean ν is negative, and 2iiC if the mean of the maximum ν is negative. Type 3 is satisfied if a material exhibits a class completely, when any axis of deformation is considered. Therefore, if every axis examined is of a particular class then type 3 of this class will be satisfied for the material.

To allow comparisons between materials and enable us to identify materials which may be of interest, each classification in this typology is matched with a value relevant to the criteria used to establish its satisfaction. Table 2 details the values which accompany each classification of the auxetic behaviour of a system. they should exhibit isotropic auxetic behaviour when part of randomly arranged polycrystalline structures.

Case Study for general materials

To demonstrate the effectiveness of the classification system when used as a comparison criterion between different groups of materials and to establish a baseline of auxeticity, we first study a wide range of different materials. The experimental elastic constants of the same 471 materials used in the work on the correlation between extreme Poisson's ratios and anisotropy by Lethbridge 8 are used to derive auxetic typologies, as they are considered to be representative of the auxetic properties of single crystals in general. Table 3 shows the number of materials for each typology. Crystals exhibiting weak auxetic behaviour of type 1A are surprisingly common, 37.2% of those studied. However, the number of materials satisfying increasing levels of auxeticity decreases rapidly, with only 3.2% 1B, 6.6% 2iA, and 7.9% 2iiA. The 2iB, 2iiB and 2C categories are populated by a single crystal, α-cristobalite, with auxetic coverage of 99% (percentage of axes with negative average PR) -0.14 (average of average) and -0.002 (average of maximum). Three single crystals belong to the 3A category: α-cristobalite again, and two forms of Beryllium Copper alloy with auxetic values (maximum of the minimum) of -0.063, -0.003, and 0.0 respectively. None of the 471 single crystals has 3B or 3C auxeticity. We do not observe any correlation between simpler, non-elastic tensor related properties and the presence or extent of auxeticity. In particular, and somewhat counterintuitively, crystal systems, space groups or even atomic structures are not predictors of auxeticity: for instance, 57% of the BCC metals are auxetic, but 43%, with the same atomic structure, are not 21 .

α-Cristobalite

As seen from this meta-analysis, the auxeticity of α-cristobalite is very peculiar, not because of the low values its Poisson's ratio reaches (-0.51, quite modest), but because of the large range of directions for which it has a negative value. The silicon dioxide polymorph cristobalite has two crystalline phases, a low temperature (α) form and a high temperature (β) form. Brillouin spectroscopy on a single crystal of α-cristobalite showed it to be highly anisotropic and yet still have an aggregate isotropic Poisson's ratio of -0.133, -0.191 and -0.163 for the Reuss, Voigt and Hill averages respectively 10 . This study also showed that the crystal had a shear modulus of roughly 2.4 times that of the bulk modulus. Due to the instability of β-cristobalite, its elastic properties have only been derived through the use of computational modelling. A molecular dynamics study showed that the high temperature phase of cristobalite remains averagely auxetic 22 . Structurally α-cristobalite is analogous to a tetragon formed by four smaller tetrahedra which are able to rotate and dilate.

The unique auxetic nature of α-cristobalite is highlighted by our classification system. Within the auxetic typology it satisfies all classifications except 3B and 3C. Considered in context as a member of the materials previously examined in this study, it is the only one of the 471 to exhibit 2iB, 2iC, 2iiB, and 2iiC auxeticity and one of the three materials to satisfy type 3A. Through use of the classification system proposed α-cristobalite is confirmed as extreme, even among other materials exhibiting a negative Poisson's ratio.

Atomistic simulations of all-silica-zeolites

Methods

In order to extend our sample of materials, we calculated the elastic properties of all-silica zeolites as well as α-quartz and αcristobalite using classical and ab initio simulations. Crystallographic information files (.cif) of the 206 known zeolite frameworks were obtained from the international zeolite association [START_REF] Baerlocher | Database of zeolite structures[END_REF] . The classical calculations were performed with GULP v3.1 [START_REF] Gale | [END_REF] , using the five well established potential models of Catlow, Gale, Matsui, Sastre and Van Beest [25][26][27][28][29] . These models have different characteristics (partial/formal charges, shell models) and have been shown in a recent review to perform well for elastic properties of α-quartz 30 . The first principles calculations were performed with the CRYSTAL09 code 31 within the Density Functional Theory Theory (DFT) approach, using the B3LYP hybrid functional 32 with empirical correction for the dispersive interactions 33 , full symmetry and all electron localized basis sets 34 . The stiffness matrixes for each structure † were used to build a database for the software ElAM [START_REF] Marmier | [END_REF] to examine. ElAM calculates the off axis elastic properties, allowing for the identification of interesting properties such as a negative Poisson's ratio in particular axes as well as establishing the auxetic typology. Not all pure silica frameworks could be relaxed. Depending on the potential models, between 186 and 189 were optimised with classical simulations. DFT simulations led to elastic constants for a subset of 121 frameworks, the remaining being too large for the computer resources available 9 .

This journal is © The Royal Society of Chemistry 2015 Phys. Chem. Chem. Phys., 2015, 00, 1-3 | 5

Validation

To ensure that the results are reasonable, it is critical to validate the models by comparison with experimental data. As little is available for the elastic constants of pure silica zeolites (we are only aware of MFI 11 ), we have to rely on the elastic constants of the "dense" silicas α-quartz and α-cristobalite. Detailed comparisons with experimental lattice and elastic constants for α-quartz, α-cristobalite and MFI are available in tabulated form in table S1(a-c) † . The root-mean-squared errors (RMSE) vary between 4 GPa and 29 GPa (average elastic constant at 40 GPa). DFT seems less accurate with RMSE at 29 GPa (quartz) and 15 GPa (cristobalite). The classical models fare equally well for quartz (4-16 GPa) and cristobalite (6-14 GPa). They are slightly less accurate for MFI (10-29 GPa), probably because the experiments were performed on crystals still containing the organic templates in their pores. The Catlow model is the most consistent with RMSE around 10 GPa for all three crystals. It is notoriously difficult to simulate α-quartz with DFT and elastic constants depend on precise details of the model. Another possible source of error is due to the fact that experimental elastic constants have been obtained at room temperature, while calculated ones are at 0 K.

α-Cristobalite

To further test the potential models used and to gain a better understanding of α-cristobalite, the typology classifications for this crystal from differing potential models are compared. The elastic constants of α-cristobalite are first calculated using the six potential models, resulting in six distinct stiffness matrices for comparison with experimental results. A typology classification is then created from each set of elastic data, and the associated values are shown in Table 4, where the values which satisfy typology classifications are in bold. This table shows that subtle differences in elastic stiffness values strongly influence the classification, as α-cristobalite does not reach '2C' or '3A' for either of the potential models, and that the minimum Poisson's ratio ('1A' value) is less by generally a third. We conclude that used together, not in isolation, the six models provide a reasonable, if somewhat conservative, platform to explore the occurrence of NPR in pure silica zeolites.

Auxetic typology of all silica zeolites

Table 5 shows the percentages of pure silica zeolites that exhibit a certain type of auxeticity for all six models and contrasts them with the base line of materials established previously. pores and low density are no more auxetic than normal materials. On closer inspections however, pure silica zeolites appear to have a significantly higher level of average auxeticity with increased percentages satisfying 2iiA and 1B for most models. This somewhat vindicates the hypothesis that zeolites are good candidates for auxeticity, but barely. Following a trend similar to that of the reference sample, the proportions decrease rapidly above these classifications, with few materials displaying type 3 or class C auxeticity, with one notable exception, discussed later. Comparison of the minimum and maximum calculated Poisson's ratio with elastic anisotropy confirms the work of Lethbridge et al. 8 It was found that the extreme Poisson's ratio all lie on two families of curves approximately symmetrical around a single point of intersection at A* = 1, where A* is the Ledbetter anisotropy measure defined as the square of the maximum shear-sound-wave velocity divided by the square of the minimum shear-sound-wave velocity 35 . Fig. 3 shows the graph of minimum and maximum Poisson's ratio plotted against the elastic anisotropy of the structure, and the pure silica zeolites follow the same patterns as the sample of 471 single crystals. This suggests that overall, zeolite frameworks behave no differently than denser silica crystals.

A very small number of frameworks are 2B auxetic: with Catlow GOO, JST, NPO and VFI, with Gale JST and OSI, with Sastre ABW, IWV, JST and NPO, with BKS AFO, BOF, GON and JST, with Matsui MSO and JST, and with DFT JST only.

Only one framework, JST, is 3C but crucially with all 6 models, including DFT. ; it is interesting that 5 are samarium based, which strongly suggests a chemical bonding origin for their generalised auxeticity, as discussed in 36 .

A completely auxetic zeolite: JST

In contrast, the remarkable agreement between the 6 models suggests that the completely auxetic character of JST derives chiefly from its structure, and not from details of the bonding. Pure silica JST is cubic (space group PA-3, number 205), its primitive cell, depicted in Fig. 4(a), contains 96 oxygen and 48 silicon. Due to this complexity it can be described in many ways. The following description is especially relevant to its auxeticity and helps understand the underlying mechanism. JST has not been synthesised yet, it might be necessary to obtain further confirmation of its framework's extraordinary auxeticity by mechanistic methods such as Finite Element Analysis or fabrication of macroscale models. We have designed a typology for auxetic materials, classifying materials with distinctions in the number of axes where a negative Poisson's ratio can be found, and the degree for which each axis is auxetic. The system has been demonstrated as a tool for comparing groups of materials and their auxetic properties, and the benefits of distinguishing between levels of auxeticity have also been explored, further highlighting the extraordinary properties of α-cristobalite. The elastic properties of a large number of pure silica zeolite frameworks have been calculated and classified with the proposed typology. After comparison of the auxetic classifications with a database of 471 general materials, the pure silica zeolites are shown to be marginally more auxetic, but follow the same trends. We found that the JST frameworks has great potential for complete auxeticity, and have proposed an explanatory mechanism.

Fig. 1 .

 1 Fig. 1. Cross section of materials of various auxetic typologies subjected to a tensile load in the direction normal to the page. The dotted line represents the original cross section, and the full line the final cross section. a) Conventional behaviour with only positive Poisson's ratio. b) The material meets the minimum requirements to be class 1A as this axis has at least some NPR c) The material meets the minimum requirements to be class 1B as the mean of the Poisson's ratios around this axis is negative. d) The material meets the minimum requirements to be class 1C as ν is negative for all transverse directions.

Fig. 2 .

 2 Fig. 2. Angles and vectors describing longitudinal axis (a) and transverse direction (b)

Fig. 3 .

 3 Fig. 3. Minimum and maximum Poisson's ratio vs Ledbetter Anisotropy for pure silica zeolites (crosses) and other materials (squares). Both ν and A* are unitless.

  It is first useful to abstract the bent Si-O-Si bonds by straight Si-Si segments as in Fig.4(b). JST can then be seen as face centred cubic, where the motifs are four three-dimensional 6 sided stars. These stars (12 silicon each), are composed of a quasiplanar central hexagon surrounded by 6 sides alternatingly pointing up and down at angles of ±60.4° (in the Catlow model). The central hexagons lie in the four ሼ111ሽ planes, and connect to each other through connecting triangles; in other words, each silicon belongs to a single star. If one assumes very stiff Si-Si connection, a tensile deformation flattens the stars, and makes them larger, leading to in-plane NPR. This simple mechanistic analysis is confirmed by studying the deformation of a unit cell. Under a 1% (100) strain, the Si-Si bonds remain largely unchanged, most bonds deforming by 0.02% with a maximum of 0.18%, while the out-of-plane side angle varies from 60.4° to 58.9°, a change of 2.4%. At this stage, JST has only been recently produced, and in a complex gallogermanate form with cationic templates in its pores37 . While pure silica This journal is © The Royal Society of Chemistry 2015 Phys. Chem. Chem. Phys., 2015, 00, 1-3 | 7

Fig. 4 .

 4 Fig. 4. (a) Unit cell of pure silica zeolite JST (Si in light grey, O in dark); (b) Representation of JST, seen in the [111] direction. The vertices correspond to Si atoms, and the edges to Si-Si links.

Table 1 . Typology of auxetic materials

 1

	Code	Poisson's ratio is negative
	0	never
	1A	for at least one direction around at least one axis
	1B	averaged over all directions around at least one axis
	1C	for all directions around at least one axis
	2A	for at least one direction around an average of axes
	2B	averaged over all directions around an average of axes
	2C	for all directions around an average of axes
	3A	for at least one direction around all axes
	3B	averaged over all directions around all axes
	3C	for all directions around all axes

Table 2 .

 2 Auxetic typology values

	Code	Value
	0	N/A
	1A	Minimum of the minimum ν
	1B	Minimum of the average ν
	1C	Minimum of the maximum ν
	2iA	Percentage of axes of class A
	2iB	Percentage of axes of class B
	2iC	Percentage of axes of class C
	2iiA	Average of the minimum ν
	2iiB	Average of the average ν
	2iiC	Average of the maximum ν
	3A	Maximum of the minimum ν
	3B	Maximum of the average ν
	3C	Maximum of the maximum ν

Of particular interest are materials which meet the criteria for satisfying 2iiB. Because these materials show a negative Poisson's ratio when averaged over all axes and all directions, 4 | Phys. Chem. Chem. Phys., 2015, 00,

1-3

This journal is © The Royal Society of Chemistry 2015

Table 3 .

 3 Auxetic classifications for case study on general materials

	Code	Number	Percentage
	0	296	62.8
	1A	175	37.2
	1B	15	3.2
	1C	11	2.3
	2iA	31	6.6
	2iB	1	0.2
	2iC	1	0.2
	2iiA	37	7.9
	2iiB	1	0.2
	2iiC	1	0.2
	3A	3	0.6
	3B	0	0.0
	3C	0	0.0

Table 4 -

 4 Associated values for the typology classifications of αcristobalite when using the elastic constants calculated with various potential models and from experimental results. Values corresponding to auxetic behaviour in bold (negative extremum or average above 50%, see Table2).

		Cat.	Gale	BKS	Sas.	Mat.	DFT	Exp.
	1A	-0.32	-0.35	-0.46	-0.34	-0.47	-0.41	-0.51
	1B	-0.09	-0.12	-0.17	-0.04	-0.17	-0.13	-0.26
	1C	0.05	0.04	0.02	0.02	0.01	0.04	-0.07
	2iA	59%	65%	83%	51%	90%	75%	100%
	2iB	39%	48%	61%	8%	66%	49%	99%
	2iC	0%	4%	0%	0%	0%	0%	53%
	2iiA	-0.09	-0.11	-0.19	-0.01	-0.20	-0.15	-0.28
	2iiB	0.02	-0.00	-0.04	0.13	-0.04	0.01	-0.14
	2iiC	0.13	0.10	0.12	0.27	0.13	0.16	0.00
	3A	0.12	0.10	0.03	0.28	0.01	0.07	-0.07
	3B	0.12	0.10	0.13	0.28	0.13	0.21	0.00
	3C	0.18	0.15	0.23	0.41	0.25	0.39	0.10

Table 5

 5

		. Percentage of auxetic zeolites (from between 186 and 189
	zeolites, depending on model) in each auxetic typology, compared with
	equivalent percentages for reference materials (471 single crystals)
		Cat.	Gale	BKS	Sas.	Mat.	DFT	Ref.
	0	67.4	83.4	51.3	78.6	54.0	74.4	62.8
	1A	32.6	16.6	48.7	21.4	46.0	25.6	37.2
	1B	11.8	7.0	21.4	8.0	24.9	2.5	3.2
	1C	1.1	0.5	7.0	0.5	6.9	0.8	2.3
	2iA	10.7	5.9	22.5	8.6	24.9	1.7	6.6
	2iB	1.0	0.5	2.7	0.5	1.1	0.8	0.2
	2iC	0.5	0.5	0.5	0.5	0.5	0.8	0.2
	2iiA	12.8	6.4	25.1	8.6	27.5	0.8	7.9
		1.6	1.0	1.6	1.6	0.5	0.8	0.2
	2iiC	0.5	0.5	0.5	0.5	0.5	0.8	0.2
	3A	0.5	0.5	0.5	0.5	0.5	0.8	0.6
	3B	0.5	0.5	0.5	0.5	0.5	0.8	0.0
	3C	0.5	0.5	0.5	0.5	0.5	0.8	0.0
	Category 1A bears the most relevant comparison to previous
	studies, as it corresponds to a material exhibiting at least some
	negative Poisson's ratio. The proportions vary between 16.6%
	and 48.7 %, which correlate well with the 37.2% base line.
	They are significantly lower than the 69% of 1A for cubic
	metals 21 . One could be tempted to draw the counterintuitive
	conclusion that, as a class of materials, zeolites, despite their

6 | Phys. Chem. Chem. Phys., 2015, 00, 1-3 This journal is © The Royal Society of Chemistry 2015

Table 6

 6 Lattice parameter a (Å), elastic constants Cij (GPa), anisotropy index A* and auxeticity typology values 1A, 2B and 3C for JST

		Cat.	Gale	Sas.	BKS	Mat.	DFT
	a	15.21	15.43 15.12 15.25 15.83	15.38
	C11	36.5	34.1	19.6	15.6	15.4	29.41
	C12	-5.7	-6.28	-7.7	-13.1	-8.8	-8.4
	C44	21.4	20.4	13.9	14.2	13.1	18.2
	A*	1.01	1.01	1.02	1.01	1.08	1.04
	1A	-0.19	-0.23	-0.66	-5.42	-1.33	-0.40
	2B	-0.19	-0.23	-0.65	-5.34	-1.32	-0.39
	3C	-0.19	-0.22	-0.65	-5.26	-1.31	-0.38

Table 6

 6 displays the lattice parameter, elastic coefficients and typology values for JST simulated with the 6 models. It is striking to note that all models are in very good agreement and predict complete auxeticity (3C) and almost complete elastic isotropy (A*~1). On the other hand, the actual values for the Poisson's ratio, here in the form of auxeticity typology values, are less well defined, and vary from -0.2 to -5.4. The most unusual feature of the elastic coefficients is that C 12 is negative. While clearly uncommon, this has previously been observed experimentally by Brillouin scattering for Sm x La 1-x S36 .

Negative C 12 and sufficiently large C 44 are necessary and sufficient conditions for a cubic crystal to be completely auxetic † . Goldstein et al. list 6 possible completely auxetic cubic crystals (table S6 of the supporting material of

17 )

Physical Chemistry Chemical Physics Accepted Manuscript

Acknowledgements

This work was supported by the Engineering and Physical Sciences Research Council [Grant number EP/G064601/1] ], and GENCI-IDRIS HPC resources [project i2014087069]. Notes and references a College of Engineering, Mathematics and Physical Science, University of Exeter, EX4 4QF, UK. b PSL Research University, Chimie ParisTech -CNRS, Institut de Recherche de Chimie Paris, 75005, Paris, France † Electronic Supplementary Information (ESI) available: Detailed exposition of the auxetic typology and elastic constants for all calculated zeolites. See DOI: 10.1039/c000000x/ ‡ Corresponding author: a.s.h.marmier@exeter.ac.uk


Charlotte Brabant
email: charlotte.brabant@yahoo.fr

Alvarez-Vanhard Emilien

Houet Thomas
email: thomas.houet@univ-rennes2.fr

Improving the classification of urban tree diversity from Very High Spatial Resolution hyperspectral images: comparison of multiples techniques

Keywords: Tree species, urban area, reduction dimension, classification, learning methods. I

published or not. The documents may come

Improving the classification of urban tree diversity from Very High Spatial Resolution hyperspectral images: comparison of multiples techniques

Charlotte Brabant, Emilien Alvarez-Vanhard, Thomas Houet

INTRODUCTION

In the context of urban sprawling and increasing and population, urban vegetation has several interests: it provides many services such as reducing the urban heat island (UHI) through shade, evapotranspiration and tree photosynthesis [START_REF] Akbari | Cool surfaces and shade trees to reduce energy use and improve air quality in urban areas[END_REF][START_REF] Bolund | Ecosystem services in urban areas[END_REF]. It improves air quality or even habitat biodiversity [START_REF] Santamouris | Using cool pavements as a mitigation strategy to fight urban heat island-A review of the actual developments[END_REF]. Urban vegetation also has disservices such as allergies due to pollen, constitutes a refuge for harmful bird species for inhabitants (noise, droppings, etc...). The intensity of these services and disservices depends on the plant species, its location and structure [START_REF] Escobedo | Spatial heterogeneity and air pollution removal by an urban forest[END_REF][START_REF] Mccarthy | Drivers of variability in water use of native and non-native urban trees in the greater Los Angeles area[END_REF][START_REF] Mcpherson | Million trees Los Angeles canopy cover and benefit assessment[END_REF]. For example, the effectiveness of UHI mitigation depends on the water needs of the species concerned [START_REF] Lemonsu | Inclusion of vegetation in the Town Energy Balance model for modelling urban green areas[END_REF].

Mapping tree species therefore responds to major urban environmental issues. Remote sensing applied in urban areas makes it possible to overcome the lack of reliable and reproducible information on vegetation to evaluate more precisely ecosystem services [START_REF] Davies | Mapping an urban ecosystem service: quantifying above-ground carbon storage at a city-wide scale[END_REF][START_REF] Weng | Remote sensing of urban environments: Special issue[END_REF]. The existing literature has showed the potential of Very High Spatial Resolution (VHSR) multispectral data to discriminate trees species [START_REF] Immitzer | Tree Species Classification with Random Forest Using Very High Spatial Resolution 8-Band WorldView-2 Satellite Data[END_REF][START_REF] Key | A Comparison of Multispectral and Multitemporal Information in High Spatial Resolution Imagery for Classification of Individual Tree Species in a Temperate Hardwood Forest[END_REF][START_REF] Pu | A comparative analysis of high spatial resolution IKONOS and WorldView-2 imagery for mapping urban tree species[END_REF], as well as VHRS hyperspectral imageries [START_REF] Alonzo | Urban tree species mapping using hyperspectral and lidar data fusion[END_REF][START_REF] Voss | Seasonal Effect on Tree Species Classification in an Urban Environment Using Hyperspectral Data, LiDAR, and an Object-Oriented Approach[END_REF][START_REF] Xiao | Using AVIRIS data and multiplemasking techniques to map urban forest tree species[END_REF][START_REF] Boschetti | Tree species mapping with Airborne hyper-spectral MIVIS data: the Ticino Park study case[END_REF] associated with Lidar or not. However, the high dimensionality of VHRS hyperspectral data can alter the effectiveness of conventional methods to discriminate against urban tree diversity.

The French Space Agency is currently considering developing a new spaceborne VHRS hyperspectral sensor: HYPXIM. The objective of this study is to evaluate (1) the ability of VHSR hyperspectral data to classify urban tree diversity and (2) the most effective methods.

II. DATA

A. Study area

The study site is the city center of Toulouse (43°36'16''N, 1°26'39''E) -France. The urban area is highly artificialized and the vegetation is embedded between buildings and houses, present along roads or in courtyards.

B. Data

Airborne hyperspectral data were acquired using the HySpex sensor (2m spatial resolution, > 480 spectral bands ranging from 0.4 to 2.5µm) in July 2015. In order to assess the potential of the HYPXIM spaceborne sensor (192 bands), HySpex data have been resampled accordingly to HYPXIM specifications at 4m and 8m spatial resolution. A Sentinel 2 imagery at 10m resolution has also been simulated for the sake of the demonstration. Atmospheric corrections were applied using the Cochise method.

A field campaign has been conducted between 31 st may and 2 nd June, 2017. More than 600 trees have been inventoried and geolocated, using a differential GPS (with an absolute geometric accuracy below 10 cm). It is assumed that the tree composition did not changed between the airborne acquisition (2015) and the field survey (2017). In addition, the dataset is very unbalanced: for e.g. Platanus x hispanica have 144 individuals unlike other species like Prunus x avium which have only 7 individuals. Although this unbalanced dataset is representative of the data generally found in urban areas, it can distinguish 12 to 19 families, and from 14 to 27 species respectively, while considering a minimum number of 4 and 10 individuals for each of them.

III. METHODS

Figure 1 summarizes the overall approach to assessing the influence of reduction dimension, learning and classification methods. Their influence is evaluated on the classification performance for two levels of nomenclature (Species and Families), themselves subdivided into 2 sublevels depending on the minimum number of individuals per classes to be considered (respectively 4 and 10).

A. Reduction dimention methods

General methodological approach

Several dimension reduction methods have been applied on the four considered dataset: PCA (Principal Component Analysis), MNF (Minimum Noise Fraction) and an innovative method to select the most discriminant meaningful vegetation indices. 96 spectral vegetation indices (VIs) inventoried by [START_REF] Erudel | Criteria Comparison for Classifying Peatland Vegetation Types Using In Situ Hyperspectral Measurements[END_REF] were computed for the three hyperspectral datasets only. It is assumed that these VIs characterizing the biogeochemical properties of plants can contribute to better discriminate trees. The selection of VI was made using two methods: (1) PCA method is used to eliminate the correlation between VIs. Kaiser's criterion [START_REF] Kaiser | The Application of Electronic Computers to Factor Analysis[END_REF] is used for the selection and components with a eigenvalue greater than 1 are selected; (2) The K-best method uses the F-value of ANOVA (ANalysis Of VAriance) as a ranking criterion. From a correlation matrix (Pearson's r), if the correlation between two VIs is greater than 85% (r > 0,85), the one with the smallest F-value is removed from the selection.

For the sake of the demonstration, corresponding results are provided only for the sub-levels of the 19 and 12 families and compared using the Overall Accuracy (OA) and the Kappa Index of Agreement (KIA).

B. Training data and methods

Each classification is trained and validated using with stratified Kfolds method applied thanks to the in situ dataset. For each individual, vectors were extracted from reduced and unreduced images. It is randomly partitioned into k subsamples with the same proportion of individuals per class. The first method uses Kfolds = 2 (2-folds): the dataset is separated into 50 % training and 50 % validation samples. The second method defines Kfolds = 10 (10-folds), meaning that 10% is used for the validation and 90% for the training, iteratively. This 10-folds methods is applied only for the sublevel with minimum individuals = 10. Each of these two methods is applied 10 times for each dataset.

For the sake of the demonstration, the 2-folds method is used as a training method on all reduced and unreduced image datasets while the 10-folds is applied only on the one exhibiting the best OA with the 2-folds and for the sub-level of the 14 species (i.e. with at least 10 individuals).

C. Influence of classification methods

Two commonly used classifiers were tested and compared: SVM (Support Vector Machine) and RF (Random Forest). SVM classifier is adopted here as it might be sensitive to small and non-homogeneous learning datasets compared to RF for instance [START_REF] Weng | Remote sensing of urban environments: Special issue[END_REF]. SVM parameters (kernel functions, constant C, etc.) have been optimized because they can affect the quality of the training [START_REF] Camps-Valls | Kernel Methods for Remote Sensing Data Analysis[END_REF]. The RF classifier is also known to be a good classifier for hyperspectral data [START_REF] Ham | Investigation of the random forest framework for classification of hyperspectral data[END_REF] because it tends to select or weight characteristics during learning [START_REF] Dalponte | Tree Species Classification in Boreal Forests With Hyperspectral Data[END_REF]. Three parameters of RF have been optimized here to increase the predictive power of the model: the number of tree forests, the number of maximum features, and the criterion.

IV. RESULTS

A. Influence of reduction dimension and classification methods

The number of selected components from MNF and PCA is respectively of 40 and 30 for all hyperspectral images and of 9 and 13 components for Sentinel 2, representing at least 98% of the variance. Using the PCA on VIs provides a selection of 19 vegetation index for HySpex 2m (explaining 86.16% of the variance), 19 for HYPXIM 4m (84.5%) and 17 for HYPXIM 8m (85.89%). When using the ANOVA and the correlation matrix, 45 uncorrelated VIs are selected for HySpex 2m, 39 for HYPXIM 4m and 38 for HYPXIM 8m. Indices with the highest value regardless of the resolution mainly rely to the Chlorophyll content with 24/45 selected indices for HySpex, 19/39 for HYPXIM 4m and 24/38 for HYPXIM 8m.

Reduced datasets using MNF provide best results no matter the nomenclature, the number of individuals, classifiers or datasets used (Table 1). For the "12 Families" sub-level (i.e. individuals > 10), the mean OA is greater than 60%. VI selection have close results with OA = 58.17% (KIA = 0.50) for the HySpex image and RF. PCA provides worst results than any other dimension reduction methods. Reduced datasets provide always better results than nonreduced ones, except for Sentinel 2 image in some cases. Both SVM and RF classifiers offer similar results. Even if RF gives generally the best OA regardless of the nomenclature, the highest OA are obtained with SVM with at least 10 individuals per class, ranging from 59% to 65% with MNF reduced images.

B. Influence of learning methods

The 10-folds learning method applied on MNF image datasets strongly improves OA: +15% for the HySpex 2m, +14% for HYPXIM 4m, +4% for HYPXIM 8m and +9% for Sentinel 2 compared with the 2-folds method (Figure 2). Similar results are obtained with HySpex 2m and HYPXIM 4m with the SVM classifier (OA = 78.38%, KIA = 0.73). Classification accuracies are similar with RF: 75.68% with HySpex 2m and 72.97% with HYPXIM 4m. Those for HYPXIM 8m and Sentinel2 datasets do not exceed 65% (OA = 64.86% and 62.16% respectively with RF).

V. DISCUSSION

This study confirms the literature on the influence of the training data and methods. Interclass and intraclass variability inheriting from an unbalance dataset can have strong effects. Although it is quite common in urban area to observe few dominating species due to urban planning and management strategies and individual choices of inhabitants, a greater number of folds prevent from over-fitting that may inherit from unbalanced urban training datasets. Conversely, depending on the considered objectives, the choice of hyperspectral techniques can be adapted. For instance, we could recommend the use of SVM classifier, applied on MNF reduced image data with a limited number of samples per class if the objective is to map a dominant allergenic tree species. In any case, this study highlights that training samples have to be clearly designed and optimized prior to hyperspectral processing.

Another interesting result comes from the evaluation of the respective influence of hyperspectral methods for improving the classification of tree diversity in cities. Training methods have a greater influence than reduction dimension and classification methods.

This study highlights that VHRS hyperspectral imageries with a spatial resolution of 4 or 8 meters can be very effective for urban tree diversity classification. First of all, it equals or overcomes OA rates obtained in others studies for AVIRIS, AISA or MIVIS hyperspectral sensors exhibiting a finer spatial resolution with a similar / greater number of tree families [START_REF] Alonzo | Urban tree species mapping using hyperspectral and lidar data fusion[END_REF][START_REF] Voss | Seasonal Effect on Tree Species Classification in an Urban Environment Using Hyperspectral Data, LiDAR, and an Object-Oriented Approach[END_REF][START_REF] Xiao | Using AVIRIS data and multiplemasking techniques to map urban forest tree species[END_REF][START_REF] Boschetti | Tree species mapping with Airborne hyper-spectral MIVIS data: the Ticino Park study case[END_REF]. While one can expect an improvement of tree classification with an increasing spatial resolution, results suggest that HYPXIM (4m) can be as appropriate as HySpex (2m) imageries. A finer resolution may involve a greater spectral intra-class variability in tree crowns that can lead to misclassifications [START_REF] Dalponte | Tree Species Classification in Boreal Forests With Hyperspectral Data[END_REF]. Indeed, spectral reflectance may vary within a single tree, depending on biochemistry and water content [START_REF] Cochrane | Using vegetation reflectance variability for species level classification of hyperspectral data[END_REF][START_REF] Ustin | Retrieval of foliar information about plant pigment systems from high resolution spectroscopy[END_REF], solar exposure of leaves, and shaded effects of street orientation or neighbouring building height. Figure 3 shows that at 2 meters resolution, shadows from building or even from London plane alignments are misclassified as Linden species, which is the second dominant tree species in Toulouse.

Finally, obtained results are promising, although there is still room for improvement. Indeed, July may not have been an optimal period for hyperspectral data acquisition and urban tree species differentiation. Spaceborne hyperspectral time series of hyperspectral imageries would undoubtly contribute to a huge step forward into urban tree species classification.

VI. CONCLUSION

The objective of this study was to evaluate the ability of VHSR hyperspectral data to discriminate between urban tree species and/or families and the related most effective methods. It shows that (1) dimensional reduction methods are essential, (2) learning methods contribute significantly to improving the classification of urban tree species and even more (3) the type of classifiers -SVM and RFdoes not have a significant influence. It stressed the importance of the in situ dataset which requires particular attention due to the high diversity of tree species in urban areas and with regard to the expected objectives. A spaceborne hyperspectral VHRS (4 meters) sensor would be particularly promising for urban environmental issues as times series would contribute to a huge step forward. 2) a zoom-in on the "Grand-Rond" urban park exhibiting one of the greatest urban tree diversity.

Fig. 1 .

 1 Fig. 1.

Fig. 2 .

 2 Fig. 2. Results of SVM (top) and RF (bottom) classifiers applied on MNF reduced datasets with 2-Folds and 10-Folds learning methods with the 14 Species nomenclature.

Fig. 3 .

 3 Fig. 3. Results of SVM classifications applied with MNF reduced datasets with 10-Folds method with (a) HySpex 2 meters resolution imageres and (b) HYPXIM 4 meter resolution, on (1) the entire study area and with (2) a zoom-in on the "Grand-Rond" urban park exhibiting one of the greatest urban tree diversity.

TABLE I .

 I CLASSIFICATION RESULTS USING SVM AND RANDOM FOREST ON THE TWO SUB-LEVELS OF FAMILY TREES DIVERSITY, USING VARIOUS UNREDUCED AND REDUCED HYPERSPECTRAL DATASETS AND THE 2-FOLDS LEARNING METHOD.

ACKNOWLEDGMENT

This study has been founded by the National Research Agency of France -ANR HYEP 'Hyperspectral imagerY for Environmental urban Planning" ANR 14-CE22-0016-01.


Julien Champ

Alexis Joly

Pierre Bonnet

Fine-grained Visual Faceted Search

Keywords: H.3.3 [Information Systems]: Information Search and Retrieval search, visual, facet, image, browsing, fine-grained, plant, species, identification

INTRODUCTION

Faceted search has become the de facto search methodology for content rich web applications, in particular online stores and digital libraries. A faceted search system allows each item in the database to be assigned multiple classifications, enabling the classifications to be ordered in multiple ways, rather than in a single, pre-determined, taxonomic order. Facets are typically derived from pre-existing fields in the item's metadata or manually built from the analysis of the items. Using visual facets complementary to textual ones has already been addressed in a few previous works [START_REF] Henry | Multifacet: A faceted interface for browsing large multimedia collections[END_REF][START_REF] Bartolini | Integrating semantic and visual facets for browsing digital photo collections[END_REF][START_REF] Muller | Visualflamenco: Dependable, interactive image browsing based on visual properties[END_REF]. The principle is to automatically compute the facets by image clustering techniques and to display image links rather than textual links in the graphical user interface. As discussed in [START_REF] Henry | Multifacet: A faceted interface for browsing large multimedia collections[END_REF], visual faceted search might overcome several limitations of text oriented facets including incompleteness and inaccuracy of the meta-data. In this paper, we stress another interest of visual faceting in the case of expert data to be searched by any user. We focus in particular on the problem of identifying the species of an observed living organism without any expertise in botany or Permission to make digital or hard copies of part or all of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage, and that copies bear this notice and the full citation on the first page. Copyrights for third-party components of this work must be honored. For all other uses, contact the owner/author(s). Copyright is held by the author/owner(s). MM'14, November 3-7, 2014, Orlando, Florida, USA. ACM 978-1-4503-3063-3/14/11. http://dx.doi.org/10.1145/2647868.2654875.

zoology. This challenge, sometimes known as the taxonomic gap, is recognized as a major issue towards setting up effective and sustainable ecological surveillance tools. Determining the scientific name of a plant or an animal is actually still the primary key to index or access related data. But unfortunately, the usage of expert identification tools such as dichotomous keys and online floras requires a deep knowledge of the morphological description language. Even high-level categories such as the main types of fruit (achene, drupe, capsule, nut, etc.), leaves (compound ,bipenates, etc.) are not understood by most users. In this paper, we introduce two complementary visual faceted search mechanisms in order to widen access to such finegrained attributes. The first one attempts to visually illustrate pre-existing expert facets in the data whereas the other one attempts to discover new facets from the visual contents. Both rely on the offline construction and clustering of a dense visual matching graph computed at the region level (contrary to previous visual faceted search proposals [START_REF] Henry | Multifacet: A faceted interface for browsing large multimedia collections[END_REF][START_REF] Bartolini | Integrating semantic and visual facets for browsing digital photo collections[END_REF][START_REF] Muller | Visualflamenco: Dependable, interactive image browsing based on visual properties[END_REF] that relied on image-centric clustering and global visual features).

AUTOMATIC EXTRACTION OF FINE-GRAINED VISUAL FACETS

Visual facets are automatically extracted from each of the sub-collection of pictures used in our application: leaves, flowers, fruits and trunks (containing themselves categories that might be visually and semantically very similar). This is done through the following steps: Hash-based K-NN search. The approximate K-Nearest Neighbors (K-NN) of each local feature xj belonging to an image I of the dataset are computed efficiently thanks to the hash-based multi-probe search method introduced in [START_REF] Joly | A posteriori multi-probe locality sensitive hashing[END_REF]. Its principle is to train an adaptive search model at indexing time through kernel density estimates computed on exact K-NN samples. This model is used at search time to select the most probable buckets to be visited so as to retrieve on average a fraction α of the real K-NN's (α was set to 0.80 in all our experiments). The advantage of this method compared to other state-of-the-art methods (such as PQ-code [START_REF] Jegou | Product quantization for nearest neighbor search[END_REF]) is that it allows controlling accurately the quality of the retrieved K-NN while being applicable to any quantization function and metrics. In our case, the original scheme is transposed to the use of a more effective hash function (RMMH [START_REF] Joly | Random maximum margin hashing[END_REF]) and to the application of the Hamming Embedding principle [START_REF] Jegou | Hamming embedding and weak geometric consistency for large scale image search[END_REF]. Matching graphs construction. A consistent matching graph is computed efficiently on top of the raw visual NN's using the RANdom Sampling and Search algorithm introduced in [START_REF] Letessier | Consistent visual words mining with adaptive sampling[END_REF]. It allows matching efficiently feature sets with local geometric constraints by iteratively and adaptively selecting features and their neighborhood in the whole collection. As plants are not rigid objects, we however relaxed the affine geometric constraints used in the original method and used instead neighborhood constraints with an adaptive radius computed as the distance of the m-the spatial neighbor of the primary searched local feature. A matching graph is computed in this way for each of the used sub-collection of pictures (i.e. leaves, flowers, fruits and trunks). Matching graphs clustering The matching graphs produced by the previous step are then processed by the MCL algorithm [START_REF] Van Dongen | A cluster algorithm for graphs[END_REF] in order to group image regions in relevant clusters. To derive facets potentially useful at different granularities, MCL is run several times on each matching graph but with different inflation values (typically at 3 granularities).

FACETED SEARCH ENGINE AND GUI

Automatically extracted visual facets are indexed through Apache Solr faceting module as any structured field (the cluster identifier being used as field value). Plant species are considered as the main index entries (i.e. as documents using Solr's terminology). Each species is then associated with a number of multi-valued visual facets, corresponding to the cluster id's of the images of that species in each of the MCL partition (and each organ). Complementary, available text-based facets are also integrated in the faceted search module (in our case <Inflorescence type>, <Infruitescence type>, <Flower colour>, <Leaf type>). Both the automatically extracted visual facets and the text-based facets are displayed to the user using representative images for each facet value (rather than texts like in classical applications). The representative image of each facet value is dynamically computed from the fine-grained matching graph so as to adapt it to the previous user search actions. This is simply done by selecting the most connected node (i.e. image region) among all the nodes belonging to the species that were not filtered before. The images containing the selected nodes are roughly cropped before being displayed on the GUI according to the position and size of the selected image region. Before being presented to the user, all facets (textual and visual) are dynamically reranked according to their filtering efficiency on the remaining species. This is done by computing the following normalized entropy as relevance score of each facet:

S(X) = - n i=1 p(xi).log2(p(xi)) n
where n is the number of values observed in the search results for facet X and xi is the i-th observed value (e.g. value red for facet flower color ). Ranking all facets according to that score allows to highlight the most informative while penalizing the ones with too much possible values. The number n of observed values is displayed to the user on the GUI as well as the number of species matching each facet value.

DATA AND DEMO

Our fine-grained visual faceting scheme is demonstrated on a set of 500 plant species living in France area correspond-ing to the species included in CLEF plant identification task 2014 (within LifeCLEF lab [START_REF] Joly | Lifeclef: Multimedia life species identification[END_REF]). Textual facets were extracted from TelaBotanica1 eFlore sheets (leaf, flower and fruit type, flower color, vernacular name). Full text search is performed on the 500 wikipedia articles of the selected species. The demonstration of the system presented at ACM MM 2014 conference will be organized as a user trial, involving a volunteer tester among the audience for each demo session. The volunteer tester will be provided a picture of a plant or a plant organ and will be in charge of determining what it is. The amount of time to succeed and the logs will be analyzed afterwards and contribute to the evaluation of the system.

http://www.tela-botanica.org/