Juan Angel Sans

Francisco Javier Manjón

André Luis

Jesus Pereira

Javier Ruiz-Fuertes

Catalin Popescu

Alfonso Muñoz

Plácida Rodríguez-Hernández

Julio Pellicer-Porres

Vanesa Paula Cuenca-Gotor

Julia Contreras-García

Jordi Ibañez

Virginia Monteseguro

Unveiling the role of the lone electron pair in sesquioxides at high pressure: compressibility of β-Sb2O3

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The study of the different polymorphs of group-XV sesquioxides (valence 3+) has triggered the interest of the scientific community in the last years [START_REF] Divya | Synthesis, surface and photoluminescence properties of Sm 3+ doped -Bi2O3[END_REF][START_REF] Lakshmi | Carbon nanotube 'wired' octahedral Sb2O3/graphene aerogel as efficient anode material for sodium and lithium ion batteries[END_REF][START_REF] Gerdelidani | Arsenic geochemistry and mineralogy as a function of particle-size in naturally arsenic-enriched soils[END_REF]. Recent studies have predicted the possible structures of their different polymorphs by the analysis on their dynamical stability. Thus, their structures can be understood on the basis of deformations of the original defective-fluorite structure [START_REF] Matsumoto | Structures and energetics of Bi2O3 polymorphs in a defective fluorite family derived by systematic first-principles lattice dynamics calculations[END_REF][START_REF] Matsumoto | Electronic structures of dynamically stable As2O3, Sb2O3, and Bi2O3 crystal polymorphs[END_REF] by applying a vacant array along different crystallographic directions. The different polymorphs of these sesquioxides are likely related to the distortion of the charge density produced by the cationic lone electron pair (LEP), belonging to the outer ns 2 electrons, when there is a strong s-p hybridization between the cation and the anion [6-8]. However, despite the efforts that the scientific community has put to understand the structural stability and formation of the different polymorphs and new phases of these sesquioxides [START_REF] Guo | Pressure-Induced Cubic to Monoclinic Phase Transformation in Erbium Sesquioxide Er2O3[END_REF], the exact role played by the cationic LEP in the different polymorphs remains an open question. It is known that when the distortions are strong enough, the LEP effect can force the structures to form closedmolecular units separated one to each other. This is the case of the most symmetric cubic structures of minerals arsenolite (c-As 2 O 3 ) and senarmontite (-Sb 2 O 3 ), where the large influence of the LEP stabilizes a structure with a large number of voids. The behaviour of arsenolite and senarmontite has been recently studied under compression [START_REF] Guńka | Compressed Arsenolite As4O6 and Its Helium Clathrate As4O6•2He[END_REF][START_REF] Sans | Ordered helium trapping and bonding in compressed arsenolite: Synthesis of As4O6⋅2He[END_REF][START_REF] Cuenca-Gotor | Vibrational and elastic properties of As4O6 and As4O6•2He at high pressures: Study of dynamical and mechanical stability[END_REF][START_REF] Pereira | Structural and vibrational study of cubic Sb2O3 under high pressure[END_REF][START_REF] Zhao | Structural transition and amorphization in compressed α-Sb2O3[END_REF]. While arsenolite shows no pressure-induced transition, a high-pressure (HP) phase of senarmontite has been recently elucidated and labelled alpha-II [space group (s.g.) 114, P4 2 1 c, Z=4], being structurally quite similar to tetragonal -Bi 2 O 3 [START_REF] Zhao | Structural transition and amorphization in compressed α-Sb2O3[END_REF]. The less symmetric structures of group-XV sesquioxides and the behaviour of the cationic LEP in these compounds have also been studied at HP. In particular, the polymorphs belonging to the monoclinic symmetry -Bi 2 O 3 (-As 2 O 3 ), a.k.a. mineral bismite (claudetite I) have revealed a large structural stability at HP up to 20 [START_REF] Lifshitz | Anomalies of Electron Characteristics of a Metal in the High Pressure Region[END_REF] GPa [START_REF] Pereira | Structural study of α-Bi2O3 under pressure[END_REF][START_REF] Soignard | High-pressure behaviour of As2O3: Amorphous-amorphous and crystalline-amorphous transitions[END_REF]. A different behaviour has been observed in another monoclinic phase of As 2 O 3 , a.k.a. claudetite II, that shows a pressure-induced isostructural phase transition (IPT) above 2 GPa [START_REF] Guńka | Cascade of High-Pressure Transitions of Claudetite II and the First Polar Phase of Arsenic(III) Oxide[END_REF]. Finally, sesquioxides with intermediate symmetric structures (tetragonal and orthorhombic), like -Bi 2 O 3 , -Bi 2 O 3 , [START_REF] Cornei | New ε-Bi2O3 Metastable Polymorph[END_REF] and-Sb 2 O 3 are not so well studied at HP despite they can be the key structures to understand the cationic LEP activity and its contribution to HP polymorphs. In these last polymorphs, the LEP effect triggers the formation of empty cavities arranged along a certain crystallographic direction. Among these compounds, only the HP behaviour of -Bi 2 O 3 has been thoroughly studied; giving rise to a 2 nd order ferroelastic IPT above 2 GPa caused by an electronic charge redistribution around the BiO 6 polyhedra [START_REF] Cornei | New ε-Bi2O3 Metastable Polymorph[END_REF][START_REF] Pereira | Isostructural Second-Order Phase Transition of β-Bi2O3 at High Pressures: An Experimental and Theoretical Study[END_REF]. Since the cationic LEP stereoactivity is more intense in compounds with smaller unitcell volumes [START_REF] Norman | Chemistry of arsenic, antimony and bismuth[END_REF], sesquioxides formed with cations with smaller ionic radius than Bi (As or Sb) should exhibit a larger influence of the LEP. Thus, the HP study of orthorhombic -Sb 2 O 3 with similar empty linear channels as -Bi 2 O 3 , can be very interesting to uncover the role played by the LEP in the formation of the different group-XV sesquioxide polymorphs and understand their pressure-induced transformations. In particular, the Sb 2 O 3 compound is used as a retardant in a number of materials [START_REF] Brebu | Effect of flame retardants and Sb2O3 synergist on the thermal decomposition of high-impact polystyrene and on its debromination by ammonia treatment[END_REF][START_REF] Si | Synergistic flame retardant effects and mechanisms of nano-Sb2O3 in combination with aluminum phosphinate in poly (ethylene terephthalate)[END_REF], as oxidation catalyst in several production processes [START_REF] Youk | Preparation and Polymerization of Ethylene 2,6-Naphthalenedicarboxylate Cyclic Oligomers[END_REF][START_REF] Venkateshwarlu | Antimony Trioxide as an Efficient Lewis Acid Catalyst for the Synthesis of 5-Substituted 1HTetrazoles[END_REF], as anode material for Li-and Naion batteries [START_REF] Alifantis | High Energy density lithium batteries: materials, engineering, applications[END_REF][START_REF] Park | Li-alloy based anode materials for Li secondary batteries[END_REF], to increase stability and decrease wear in fluid lubricants [START_REF] Zabinski | Mechanistic study of the synergism between Sb2O3 and MoS2 lubricant systems using Raman spectroscopy[END_REF], to enhance the properties of glasses [START_REF] Ghosh | Transport properties of semiconducting CuO-Sb2O3-P2O5 glasses[END_REF][START_REF] Gómez | Nonresonant third-order nonlinearity of antimony glasses at telecom wavelengths[END_REF], and also as dopant to modify the conductivity of semiconductors, like SnO 2 and ZnO [START_REF] Koivula | The effect of dopant's valence (+III and +V) on the anion/cation uptake properties of antimony-doped tin dioxide[END_REF][START_REF] Samanta | Raman scattering studies of p-type Sb-doped ZnO thin films[END_REF][START_REF] Liu | Donor behaviour of Sb in ZnO[END_REF][START_REF] Kumar | Optical evidence for donor behaviour of Sb in ZnO nanowires[END_REF]. Furthermore, the synthesis, properties, and applications of nanoparticles of antimony oxides have been recently reviewed since nanoparticles possess excellent properties as compared to bulk material, like a higher refractive index, absorbability, abrasive resistance, and proton conductivity [START_REF] Chin | Review on oxides of antimony nanoparticles: synthesis, properties, and applications[END_REF]. In the literature, previous HP Raman scattering (RS) [START_REF] Geng | High-pressure Raman investigation of the semiconductor antimony oxide[END_REF][START_REF] Geng | Experimental Observation of Phase Transition in Sb2O3 under High Pressure[END_REF] and an X-ray diffraction (XRD) studies [START_REF] Zou | Pressure-induced anomalies and structural instability in compressed β-Sb2O3[END_REF], have shown that -Sb 2 O 3 undergoes a structural phase transition above 15 GPa. While those works have studied the overall compression of valentinite up to 40 GPa, showing the presence of two polymorphs up to that pressure, they have not carefully studied the behaviour of the low-pressure (LP) phase of -Sb 2 O 3 prior to the structural phase transition to the HP phase, that was proposed to be of monoclinic symmetry [START_REF] Zou | Pressure-induced anomalies and structural instability in compressed β-Sb2O3[END_REF]. This has motivated us to perform a detailed experimental and theoretical study of the structural, vibrational and electronic properties of -Sb 2 O 3 under compression with a thorough analysis of the theoretical electronic charge topology of the LP phase. We will show that compressed valentinite reveals several internal transitions, prior to its structural phase transition above 15 GPa. Most of these internal transitions will be shown to be related to changes in the cationic LEP under compression and could be considered as IPTs of order higher than 2. The incompleteness of the Landau and Liftshitz's description of solid phase transitions [START_REF] Landau | Statistical physics[END_REF] was already stated by Liftshitz in 1960 with the inclusion of the electronic topological transitions (ETTs) as transitions of order 2 ½ (higher than 2) [39] according to the Ehrenfest's notation [START_REF] Lifshitz | Anomalies of Electron Characteristics of a Metal in the High Pressure Region[END_REF]. The 2 nd order phase transitions requires a soft mode, according to Landau's theory; however, the changes observed in this work are associated to changes in the interaction of the cationic lone electron pairs, i. e. with an electronic origin such as the one observed in ETTs, which are ascribed to transitions of order higher than 2. Recent studies [START_REF] Jones | Topological Catastrophe and Isostructural Phase Transition in Calcium[END_REF] have proved the existence of IPTs of order higher than 1 (continuous), some of them associated with changes in the electronic topology.

Therefore, this work shows the importance of the cationic LEP activity in relation to the rest of electrons in the molecular unit in order to understand the properties of this quasi-molecular material. Moreover, thanks to the combination of single-crystal and powder XRD measurements, we will propose an orthorhombic HP phase [s.g. Pcc2] above 15 GPa, that is different to the monoclinic HP phase previously proposed [START_REF] Zou | Pressure-induced anomalies and structural instability in compressed β-Sb2O3[END_REF]. Consequently, this study allows completing the understanding of the polymorphism of group XV sesquioxides and becomes interesting for the future analysis of other LEP-based molecular and quasi-molecular solids under extreme conditions.

Results

Structural properties of -Sb 2 O 3 under compression

The structure of -Sb 2 O 3 [s.g. 56, Pccn, Z=4] at room pressure has been traditionally interpreted as formed by SbO 3 E units (see the SbO 3 polyhedral unit at the bottom part of Fig. 1), where the E refers to the cationic LEP. The link of different SbO 3 E units along the c-axis (two-fold rotation axis) form rods composed by infinite (Sb 2 O 3 ) ∞ chain molecules in such a way that neighbour LEPs form empty linear channels parallel to the rods. In this context, -Sb 2 O 3 can be considered as a quasi-molecular solid (with infinite molecules along the c-axis) or an acicular solid since it grows in form of (Sb 4 O 6 ) ∞ needles extended along the c-axis (see Fig. 1). Please do not adjust margins Please do not adjust margins lines in the top right part of Fig. 1). On the other hand, each of the external O atoms of each Sb 2 O 3 unit along the chain, hereafter named O2, is bonded to a Sb atom of the same unit (bond distance of 1.98 Å) and to one Sb atom of the neighbour Sb 2 O 3 unit (bond distance of 2.02 Å). Additionally, each O2 atom is weakly bonded (bond distance of 2.5 Å) to a Sb atom in a neighbour chain (see almost vertical dashed lines in the top part of Fig. 1). This weak bond provides the link between the different chains or rods and leading to the formation of the 3D structure of valentinite. In summary, by taking these weak extra bonds into account, the structure of valentinite could be interpreted as formed by rods or needles with SbO5E units, where the coordination of Sb is 3 (with 2 extra-long bonds) and the coordination of O1 and O2 atoms is 2 (with 2 extralong bonds) and 2 (with 1 extra-long bond), respectively (see the polyhedral SbO 5 E unit at the top part of Fig. 1).

Figure 2.

Pressure dependence of the experimental (symbols) and theoretical (solid lines) unit-cell volume (left) and lattice parameters (right) in -Sb2O3. Filled squares and empty circles correspond to data of experiment 1 and experiment 2, respectively. The 3rd order BM EoS fit of experimental volume data (dashed line) is also shown. Experimental data fit of lattice parameters to the modified Murnaghan's EoS (dashed lines) are also shown.

Two angle-dispersive powder HP-XRD measurements of -Sb 2 O 3 at room temperature were conducted with a different pressuretransmitting medium (PTM). In experiment 1 (2), a 16:3:1 methanolethanol-water (MEW) mixture (helium) was used. New reflections, not corresponding to the LP phase of valentinite, were only observed above 15 GPa in both experiments (see Fig. S2 and Fig. S3). This result suggests that quasi-hydrostatic conditions were maintained up to 15 GPa in both experiments and points to a phase transition above this pressure, in good agreement with the literature [35-37]. The phase transition is more evident in the case of enhanced quasi-hydrostatic conditions (experiment 2) than in the case of a quasi-hydrostatic environment (experiment 1). A pure HP phase was recorded above 16 GPa in experiment 2, while experiment 1 revealed the coexistence of LP and HP phases from 16 GPa up to 25 GPa. The analysis of the HP phase will be discussed below. Le Bail analysis (partial Rietveld refinement) of powder XRD patterns of the LP phase of valentinite in experiment 1 (experiment 2) were conducted up to 14.6 GPa. For the refinement, the atomic coordinates of oxygen were fixed to those obtained by ab initio simulations in order to get more accuracy in the determination of the coordinates of Sb atoms. An example of the 2D diffraction pattern of experiment 2 and details of partial Rietveld refinement at selected pressures are shown in Fig. S5 and Table S1. The analysis of the LP phase of -Sb 2 O 3 shows a smooth behaviour in the pressure dependence of the unit-cell volume up to 14.6 GPa (Fig. 2(left)), without any jump in volume (associated to a 1st order phase transition) or any apparent change in compressibility (associated to a 2 nd order phase transition). In fact, the data were easily fitted to a 3rd order Birch-Murnaghan (BM) equation of state (EoS) using EoSFIT software [START_REF] Angel | EosFit7c and a Fortran module (library) for equation of state calculations[END_REF]43]. Interestingly, the unit-cell volumes obtained at similar pressures in both experiments using MEW and helium are almost identical. This result contrasts with the large volume difference observed in arsenolite (c-As 2 O 3 ) at similar pressures with these two PTMs due to a pressure-induced helium trapping in this molecular solid [START_REF] Guńka | Compressed Arsenolite As4O6 and Its Helium Clathrate As4O6•2He[END_REF]. The bulk modulus obtained for -Sb 2 O 3 in experiment 1(2) is B 0 =24(2) GPa (B 0 =27.4(8) GPa). These values are in good agreement with our theoretical simulations (B 0 =30.9(5) GPa). The difference between both experiments is likely due to the lack of experimental data below 3 GPa in experiment 2 (due to the helium loading procedure). The agreement between experiments and calculations supports the validity of our theoretical calculations performed under the GGA-PBEsol approximation and will allow us to extend our analysis to different intrinsic properties of this compound. Such a point is crucial since we will show below that there are several IPTs of order higher than 2. The HP dependence of the lattice parameters of -Sb 2 O 3 up to 14.6 GPa (Fig. 2(right)) exhibits a monotonic sublinear behaviour. Both experimental and theoretical data can be fitted to a modified 2 nd order BM EoS (see axial compressibilities and bulk moduli at room pressure in Table S3). The similarity of data from both experiments suggest that quasi-hydrostatic conditions are well maintained below 14.6 GPa in both XRD measurements. Consequently, the slight differences between the lattice parameters obtained at LPs for both experiments are likely due to an overpressure of the Cu signal in experiment 2 (range between 4 and 6 GPa) since the Cu powder was too close to the gasket. Note that the ratio between the different lattice parameters (Fig. S6) exhibits a monotonic variation without any indication of a possible 2 nd order IPT up to 14.6 GPa, unlike its counterpart -Bi 2 O 3 . [START_REF] Cornei | New ε-Bi2O3 Metastable Polymorph[END_REF][START_REF] Pereira | β-Bi2O3 under compression: Op cal and elas c proper es and electron density topology analysis[END_REF] The c-axis shows the lowest axial compressibility at room pressure because its pressure behaviour is dominated by the longitudinal compressibility of the rods, mainly related to short Sb-O2 bond distances inside the rods that link different Sb 2 O 3 units along the caxis. We can also relate the compressibility of the c-axis to the lateral compressibility of the SbO 3 E or SbO 5 E polyhedral units. The a-axis exhibits the highest compressibility because its change is dominated by the compression of the voids between rods; i.e. the compression of the cationic LEPs mainly oriented along the a-axis. In any case, a small contribution of the transversal compression of the relatively rigid rods along the a-axis cannot be neglected due to the bending of the long Sb-O2 bonds that link the rods mainly directed along the baxis (see almost vertical dashed lines in Fig. 1). Finally, the b-axis shows an intermediate compressibility mainly given by the compression of the long Sb-O2 bonds and in a smaller extent to the transversal compression of the rods related to the longitudinal compressibility of SbO 3 E (or SbO 5 E) units along the b-axis. Note that the b-axis compressibility is not affected by the compression of the short Sb-O1 bonds (Fig. S9). Please do not adjust margins Please do not adjust margins Fig. S7 shows the comparison between our theoretical and experimental (obtained from partial Rietveld refinement of HP-XRD patterns) Sb atomic coordinates. These results further confirm that we can rely on the HP evolution of the interatomic distances, polyhedral volume, and compressibilities of the different units obtained from our theoretical calculations. Interestingly, the theoretical pressure dependence of the atomic coordinates of Sb, O1 and O2 atoms exhibit indications of the anomalous compressibility of valentinite. While some coordinates (Sb x , Sb y and O1 z ) show a monotonic trend with compression, the rest of free atomic coordinates, especially those corresponding to external O2 atoms, show a rather different behaviour depending on the pressure range.

In particular, the coordinates of O2 reveal three different pressure ranges: i) from 0 to 4 GPa, ii) from 4 to 10 GPa and iii) from 10 to 15 GPa. In addition, the coordinate Sb z remains almost constant up to 2 GPa and then increases up to 10 GPa, where another slight change of tendency is observed. In summary, the atomic coordinates of compressed valentinite indicate structural changes at 2, 4 and 10 GPa that are not detected in the evolution of the unit-cell volume and lattice parameters.

To conclude this section, we will now discuss the nature of the HP phase observed in -Sb 2 O 3 above 14.6 GPa. At this pressure, new reflections appear at 10.8, 12.2, 13.2, and 16 degrees that reflect displacive shifts of some crystallographic planes of the LP phase. We remark that the most intense reflections of the LP phase remain after the phase transition to the HP phase is finished. This result indicates that the HP phase may maintain a direct relation with the original LP phase. A tentative assignment of the HP phase of valentinite has been performed by comparison with probable structures. The best concordance was achieved with an orthorhombic phase (s.g. Pcc2, see Fig. S10). This structure was predicted to be dynamically stable in As 2 O 3 [START_REF] Matsumoto | Electronic structures of dynamically stable As2O3, Sb2O3, and Bi2O3 crystal polymorphs[END_REF] and is a translationengleiche subgroup of index 2 to s.g. Pccn of -Sb 2 O 3 . In order to corroborate this assignment, we have performed single-crystal HP-XRD measurements of mineral valentinite and ab initio theoretical simulations of the proposed HP phase.

The structure of -Sb 2 O 3 at ambient conditions, with voids extending along the (001) direction provides the compound with a laminar character and a tendency to exfoliate along the (100) plane. This tendency to exfoliate would be avoided if we could get a perfectly oriented (100) single crystal; however, the mineral nature of this sample makes very difficult to obtain a pure (100) -Sb 2 O 3 single crystal, increasing the probability to form deformations in its crystalline structure. Under compression, the reflections at high 2θ angles quickly fade probably as a result of the slide of planes in our single crystalline sample. This fact produces a huge increase of the internal reflections R int value, which is of ~30% at 16 GPa, thus preventing from a complete structural determination of the HP phase of -Sb 2 O 3 . Despite the difficulty to solve the structure at 16 GPa, we observe that the new structure also presents an orthorhombic symmetry, a = 4.590 [START_REF] Zhao | Structural transition and amorphization in compressed α-Sb2O3[END_REF] Å, b = 11.113 [START_REF] Watson | Origin of the Lone Pair of α-PbO from Density Functional Theory Calculations[END_REF] Å and c = 4.9730( 16) Å, very similar to that of the LP phase. The intensities and new reflections that emerge above 16 GPa in the powder XRD experiment cannot be explained with s.g. Pccn. However, they are readily indexed in the non-centrosymmetric s.g. Pcc2. Both phases can be distinguished by systematic extinctions related to their symmetry elements. In particular, certain conditions for systematic extinctions (h00, l = 2n; 0k0, k = 2n; hk0, h + k = 2n) do not apply in s.g. Pcc2. Unfortunately, due to the orientation of our crystal the h00 reflections are not accessible. Furthermore, according to the powder diffraction pattern at 19.3 GPa (Fig. S10), the 0k0 and the hk0 reflections are extremely weak hindering their localization in our single-crystal patterns (Fig. S11). Despite these limitations, the integration of reflections from single-crystal XRD measurements using the CrysAlis program [START_REF] Pereira | β-Bi2O3 under compression: Op cal and elas c proper es and electron density topology analysis[END_REF] indicates that data are compatible with s.g. Pcc2.

The strongest support to the structural determination of the HP phase comes from single-crystal data. Even though the internal R int value is above 30%, SHELXT [46] directly displays two sites for the Sb atom conforming a cationic structure almost identical to that of the LP structure and in perfect agreement with the structural determination from powder XRD. After several refinement iterations of the cationic structure with program SHELXL [47], the R1 value goes down to a minimum of 35%, but it is not reduced when we introduce the O positions (Table S2). Such a drawback, which is the result of the poor R int value prevents any total structural determination but confirms the structural determination made by powder HP-XRD measurements. The two atoms (Sb and O2) that are in 8e crystallographic sites in the LP phase split into two crystallographic sites each with 4e symmetry in the HP phase and the atom (O1) that is in the 4d Wyckoff site in the LP phase splits into a 2b and a 2c crystallographic site in the HP phase. The determination of four out of the six positions of the HP phase with single-crystal XRD is further validated by our ab initio theoretical simulations ( A classification of the vibrational spectrum of valentinite is provided by the theoretical phonon density of states (Fig. S12), which allows dividing it into three main ranges: i) low-frequency (below 100 cm -1 ), where the vibrations are dominated by the motion of the Sb cations; ii) middle-frequency (100-300 cm -1 ), where the vibrations are associated to the motion between polyhedral SbO 3 units as a group (both Sb and O atoms move with predominance of Sb movement below 230 cm -1 and of O motion above 230 cm -1 ); and iii) highfrequency (above 300 cm -1 ), where there are different kinds of mixed bending and stretching vibrations of SbO 3 units which mainly involve motion of oxygen atoms. Please do not adjust margins Please do not adjust margins New Raman-active modes are observed above 14 GPa, thus confirming the transition to the HP phase already found by XRD measurements. Unfortunately, the low quality of the RS data above 14 GPa, even from single-crystal data, does not help to get an unambiguous assignment of the HP phase of β-Sb 2 O 3 .

The pressure dependence of the experimental and theoretical Raman-active modes of the LP phase in valentinite (Fig. 3) shows a rather good agreement between all data, especially above 150 cm -1 where Raman-active modes are more separated. The comparison of theoretical and experimental frequencies and pressure coefficients of the Raman-active modes in valentinite (see Tables S3 andS4) has allowed us to assign in a tentative way the symmetry of many observed Raman peaks. Below 150 cm -1 , the symmetry assignment of the Raman-active modes is doubtful for most modes due to the overlapping of many theoretical modes. The pressure dependence of the theoretical IR-active and silent mode frequencies is shown in Fig. S18 and their frequencies and pressure coefficients are summarized in Table S6. The most striking feature of the pressure dependence of the vibrational modes of valentinite is that many modes show a strong non-linear pressure dependence showing changes around 2, 4 and 10 GPa. These changes agree with those observed in the free atomic coordinates of the Wyckoff positions, as already commented.

Discussion

A thorough analysis of the experimental and theoretical volumes of the SbO 5 E polyhedral units and of the empty cavities of the LP phase of -Sb 2 O 3 under compression (Fig. 4) shows that the bulk compressibility of valentinite is strongly dominated by the behaviour of its structural voids. The compression of the theoretical empty cavity volume shows a monotonic variation up to 14.6 GPa that fits perfectly to a single 3 rd order BM EoS (B 0 =16.9(2) GPa). On the contrary, a fit of the theoretical SbO 5 E polyhedron volume as a function of pressure to a single 3 rd order BM EoS in the same pressure range is not possible, thus revealing the existence of two different compressibility regions below and above 2 GPa. Fit to a 3 rd order BM EoS below (above) 2 GPa yields B 0 =48(5) GPa (B 0 =102(10) GPa). Given the existence of 2 nd order IPTs in other sesquioxides, like β-Bi 2 O 3 , we wonder whether the compressibility changes in the SbO 5 E polyhedron at 2 GPa could be considered as hidden or internal 2 nd order IPT. These adjectives indicate that this transition is not reflected in the bulk behaviour, already analysed, because it remains hidden by the strong compressibility of the linear empty cavity that dominates the unit-cell volume compressibility of valentinite.

The F-f plot, being F the reduced pressure and f the Eulerian strain, is a more sensitive tool than the unit-cell volume vs. pressure plot to observe slight changes in the compressibility [START_REF] Voit | Vibrational spectroscopic and quantum chemical study of antimony(III) oxide[END_REF][START_REF] Bader | Properties of atoms in molecules: Atoms under pressure[END_REF]. Therefore, we have analysed the F-f plot for the theoretically simulated bulk structure, polyhedral SbO 5 E unit, and the empty cavity under compression (Fig. S8) in order to get a better illustration of the different compressibility ranges of valentinite. The F-f plots have unveiled more than a critical pressure value along the stability pressure range of valentinite. The F-f plot of the bulk structure reveals a compressibility change around 10 GPa that is observed also Please do not adjust margins Please do not adjust margins in the polyhedral unit and the empty cavity. On the other hand, the F-f plot of the polyhedral unit confirms a compressibility change at 2 GPa. Finally, it must be noted that the F-f plot of the empty cavity indicates an additional change around 4 GPa. These pressure values agree with those already found in our study of the lattice dynamics and of the theoretical dependence of the free atomic coordinates. Complementary information about the different compressibility regions of valentinite is provided by the HP behaviour of the Sb-O interatomic distances inside the SbO 5 E polyhedral units (Fig. S8). A strong decrease of the theoretical long Sb-O2 distance (theoretically around 2.42 Å and experimentally around 2.5 Å at 0 GPa) is observed in the range from 0 to 4 GPa. Additionally, there is a change in the sign of the slope of the short Sb-O1 bond in the same pressure range. Note that this is the bond from the central O1 atom to the two Sb atoms of the same Sb 2 O 3 molecule. The changes in Sb-O distances indicate that the main effect of pressure up to 2-4 GPa is to bring neighbour rods closer; while above this pressure range, pressure mainly tends to increase the packing inside the rods by compressing the atoms inside the SbO 3 molecular units. Additional evidence of the extraordinary changes observed at 2, 4, and 10 GPa in valentinite are provided by a detailed comparison of the RS measurements (in single crystals and powders) and theoretical lattice dynamics calculations (Figs. 3, S15 and S16). At LPs, most experimental Raman-active frequencies show a blueshift with increasing pressure. This is also the case of the most intense Raman peak, located near 140 cm -1 at room pressure and tentatively attributed to the B 1g 2 mode (Fig. S19). However, this mode shows a negative pressure coefficient above 2 GPa (Fig. S18), which agrees with previous HP-RS measurements [START_REF] Geng | High-pressure Raman investigation of the semiconductor antimony oxide[END_REF][START_REF] Geng | Experimental Observation of Phase Transition in Sb2O3 under High Pressure[END_REF]. In fact, changes in the pressure coefficients of vibrational frequencies around 2 GPa are also observed in other theoretical Raman-active modes (Fig. 3) as well as in IR-active and silent modes (Fig. S18). In particular, the A g 2 and B 3g 8 modes show negative (positive) pressure coefficient in β-Sb 2 O 3 below (above) 2 GPa. Changes in slope of vibrational frequencies around 4 GPa are experimentally evidenced by the experimental mode near 225 cm -1 , assigned to the B 2g 4 mode. Note that this mode has a smaller pressure coefficient than the adjacent A g 4 mode so both modes overlap above 5 GPa and become separated again above 12 GPa (see detail in Fig. S26). Similarly, the theoretical B 3g 1 (A g 2 ) mode shows positive (negative) pressure coefficient below (above) 4 GPa, and the converse above this pressure. Finally, it must be stressed that many Raman-active modes (A g 1 , B 2g 4 , B 3g 6 , B 3g 7 , B 2g 8 , B 1g 7 ) and IR-active and silent modes show a change in slope at 10 GPa. In particular, most of these modes show a softening of their frequencies above this pressure, likely related to the increase of the short Sb-O2 bond distance (Fig. S9). This behaviour suggests the existence of a general instability of the LP phase of β-Sb 2 O 3 above 10 GPa (see also the strange kink in theoretical frequencies near 15 GPa). In this regard, we must recall that a complete softening (frequency going down to zero) of a vibrational mode is a signature of a "real" 2 nd order phase transition according to Landau theory of phase transitions, as it occurs in β-Sb 2 O 3 and β-Bi 2 O 3 [START_REF] Pereira | Isostructural Second-Order Phase Transition of β-Bi2O3 at High Pressures: An Experimental and Theoretical Study[END_REF][START_REF] Pereira | β-Bi2O3 under compression: Op cal and elas c proper es and electron density topology analysis[END_REF]. In this context, the calculated phonon dispersion curves of β-Sb 2 O 3 (Fig. S28) indicate that there is no complete softening of any vibrational mode along the Brillouin zone in the whole range up to 15 GPa. This result implies that we cannot consider the changes observed in the compressibility of valentinite at 2, 4 and 10 GPa as belonging to a 2 nd order IPT, so they must be ascribed to IPTs of order higher than 2. To conclude the discussion regarding the vibrational properties of valentinite under compression, we must mention that the changes in compressibility at different pressures have been confirmed not only by changes in slopes of the frequencies of vibrational modes, but also by the accurate analysis of the pressure dependence of the full width half-maximum (FWHM) of several high-intensity Raman-active modes, as obtained from HP-RS measurements in single-crystal samples (Fig. S27). The HP evolution of the FWHM in two experimental Raman-active modes (one in the low-frequency region and another in the high-frequency region) show that changes in their FWHM occur at roughly the same critical pressures where changes in the F-f plots are observed. In summary, all HP-RS results give experimental support to the conclusions obtained by the analysis of the theoretical simulations regarding the existence of hidden or internal structural changes in valentinite around 2, 4, and 10 GPa. In order to find additional evidence for the different ranges of compressibility of valentinite, we have studied the evolution of the topology of the electronic charge distribution by analyzing the atomic Bader volumes and charges in β-Sb 2 O 3 (Fig. 5). According to Bader's theory, the change in the topological distribution of the charge around an atom may reveal a variation of the molecular geometry since their electronic distribution should be invariant to pressure Please do not adjust margins Please do not adjust margins the compressibility changes observed in valentinite. Consequently, these results show that the anomalous compressibility of the polyhedral SbO 5 E unit around 2-4 GPa may be related to a charge transfer between Sb and O2 atoms. This is a clear indication that the different rods start to interact in a more intense way above this pressure range. Thus, we can consider that above this pressure range the coordination of the Sb atom reaches a real 3+1 coordination with a loss of the isolated rod character at LP; i.e. a quasi-1D to a 3D transition occurs between 2 and 4 GPa in valentinite. Therefore, changes between 2 and 4 GPa could be considered to be related to internal or hidden 2 nd order IPTs not reflected in the bulk compressibility of β-Sb 2 O 3 . Noteworthy, a 2 nd order IPT has also been found in β-Sb where the 2 nd order IPT was also revealed by the analysis of the evolution of the Bader charge under compression.

Another evidence of the occurrence of subtle compressibility changes in valentinite under compression comes from the analysis of the pressure dependence of the orientation and distortion of the cationic LEPs, which has been studied by the topology of the electronic localization function (ELF). The pressure dependence of the LEP orientation, given by the angle formed by the interacting LEPs and its common bond critical point (BCP), reveals a drastic change at 2 GPa (Fig. 6). The change in the orientation of the LEPs around 2 GPa agrees nicely with the features observed in the analysis of the atomic Bader volumes and charges and with those observed in the structural and vibrational properties of valentinite. Therefore, we can conclude that the anomalous behaviour of valentinite at 2 GPa is caused by a sudden reorientation of the cationic LEPs that form the structural voids or channels in valentinite. At room pressure, LEPs are oriented towards the centre of the empty cavity so they interact with their neighbouring LEPs belonging to different a-c crystallographic planes. The strong compression of the empty cavity along the a-axis (shown by the strong decrease of the inter-rod bond distances) shortens the distance between neighbour LEPs belonging to the same a-c crystallographic plane in such a way that these LEPs start to interact. This leads to an abrupt change in the LEP-BCP-LEP angle near 2 GPa. On the other hand, the LEP-LEP asymmetry, measuring the distortion of the two LEP-BCP distances (Fig. 6), indicates a negligible distance asymmetry between the two interacting LEPs below 2 GPa. This asymmetry increases above 2 GPa up to around 4 GPa when it becomes constant. This behaviour can be explained by the strong decrease of the voids volume and the strong reduction of the long interatomic Sb-O2 distances at LPs. Up to 2 GPa, the cation LEPs suffer a gradual reorientation remaining barely undistorted; the relative asymmetry with respect to the room pressure value is 1 (no asymmetry). Above 2 GPa, the strong compression between neighbour rods causes a sudden reorientation of the cationic LEPs and the effect of pressure begins to distort the cationic LEPs themselves. This distortion is reflected in the increase of the relative asymmetry of the LEP-BCP distances (deviation with respect to 1 that corresponds to the same LEP-BCP distance in the two interacting LEPs). This means that the distortion of the LEP cloud is larger for one LEP than for its neighbour one. The difference in the LEP distortions reaches a maximum near 4 GPa; i.e., the LEP cloud considered as an ellipsoid becomes more and more stretched up to 4 GPa. Above this pressure, the LEP cloud cannot be further compressed and the relative LEP asymmetry remains constant. This means that the effect of pressure above 4 GPa is reflected on the decrease of the long interatomic Sb-O distances and specially in the decrease the second shortest Sb-O2 distance (Fig. S8). The compression of the internal Sb-O1 bonds above 4 GPa indicates that a real 3D behaviour is found above this pressure. The anomalous hidden or internal behaviour at 2 and 4 GPa cannot be interpreted as pure 2 nd order IPTs due to the lack of discontinuities in thermodynamic parameters, which avoids their assignment in the Another point is that by definition a phase transition implies changes in atomic bonds (coordination, symmetry or arrangement); however, to our knowledge, the changes of interaction between unpaired electrons, like LEPs, in inorganic solids under compression, has never been studied up to our knowledge. In valentinite, the different interactions between neighbour LEPs may cause subtle changes in the properties of valentinite. However, these changes cannot be assigned to a 2 nd order IPT in the Ehrenfest's classification and consequently should be attributed to a transition of order higher than 2. Finally, we would like to remark that all F-f plots (Fig. S7) show a tendency change at 10 GPa, in perfect correlation with a general softening of many vibrational modes (Figs. 3 andS17). A response to compression not observed by the analysis of the electronic topology.

As already commented, the general softening of Raman-active and IR-active modes observed above 10 GPa may be the result of an ongoing instability, probably the onset of the phase transition experimentally observed above 15 GPa. However, our ab initio calculations do not show any complete softening of the theoretical vibrational modes up to 14 GPa. Note that the B 2u 1 mode near 25 cm - 1 at room pressure shows an almost complete freezing near 15 GPa. This indicates that the structure becomes fully unstable above 14 GPa when the transition occurs. For this reason, the existence of an electronic-driven IPT around 10 GPa has also been considered. Our theoretical calculations of the electronic band structure at several pressures (Fig. S29) show a progressive decrease of the bandgap of valentinite. However, neither a change from semiconducting-tometallic nor a change in the topology of the Fermi level around 10 GPa has been observed. Consequently, the origin of the instability occurring above 10 GPa in valentinite is still unknown and will deserve further studies. In any case, the theoretical interatomic distances in Fig. S9 allows an interesting reading. If we make the hypothesis that the 3+1 coordination of Sb starts at 2 GPa (when the large Sb-O2 distance decreases below 2.36 Å) the increase in coordination to 3+1+1 of Sb would start above 10 GPa (when the much larger Sb-O1 distance decreases below 2.38 Å). Moreover, Please do not adjust margins Please do not adjust margins there is a change in the theoretical slope of the two shortest Sb-O2 bond distances above 10 GPa (Fig. S9). The shortest one exhibits a reduction of the slope while the largest one exhibits an increase of the slope; i.e. it enlarges at a higher rate what suggests a structural instability. Consequently, the increase of coordination of Sb above 10 GPa, thanks to the interaction of the central parts of neighbour Sb 2 O 3 units in the same rod (involving large Sb-O1 bonds), lead to unusual changes of the Sb-O2 distance and could trigger the instability of the structure of valentinite. In the HP phase, the polyhedral unit reaches a fivefold coordination confirming the trend observed with compression. In summary, our analysis shows that the anomalous compressibility of valentinite at 2 and 4 GPa is related to changes of the cationic LEPs that affect the connectivity of external O2 atoms, while that occurring at 10 GPa is not related to cationic LEPs but likely related to changes in the connectivity of internal O1 atoms that in turn affect in a strange way the compressibility of Sb-O2 bonds. S7 and Table S1). This procedure was performed because lighter elements are more sensitive to errors in their determination by XRD measurements. Small variations of the free Wyckoff sites may lead to strong changes in interatomic distances, which can avoid an adequate precision for the experimental determination of bond lengths. The background of each XRD powder pattern has been fitted to a shifted Chebyshev function of 1 st kind and the peak shape has been analysed using the profile coefficients for Simpson's rule integration of pseudovoigt function implemented in GSAS package. For the single crystal XRD experiment a (100)-oriented single crystal was loaded into a Merrill-Bassett diamond-anvil cell, equipped with 500 μm culet diamonds. The sample with a surface of 15040 μm 2 was loaded in the pressure chamber made in the 200 μm hole of a tungsten gasket that had been preindented to 80 μm in thickness. The pressure transmitting medium was He and a ruby chip was included for pressure determination. The XRD data were collected at the Beamline I19, at Diamond Light Source, using the Newport kappa-geometry 4-circle diffractometer fitted with a Dectris Pilatus 300 K pixel-array photon-counting detector. The dataset was collected using a wavelength of λ = 0.4895 Å, with a step size and exposure time of 0.2° and 0.4 seconds respectively. The data were collected using a series of ω scans to give the optimal completeness when using a diamond anvil cell. prescription. To ensure a high convergence of 1-2 meV per formulaunit in the total energy as well as an accurate computation of the forces on the atoms, the integrations over the Brillouin zone (BZ) were performed using dense meshes of special k-points. At selected volumes, the structure was fully relaxed to the optimized configuration through the calculation of the forces on atoms and the stress tensor. In the optimized configurations, the forces on the atoms were less than 0.002 eV/Å and the deviations of the stress tensor from a diagonal hydrostatic form were smaller than 1 kbar (0.1 GPa). The set of energies E, volumes (V) and pressures (P) obtained was fitted with an equation of state (EoS) to obtain the bulk modulus and its pressure derivative. The interatomic distances were extracted with VESTA software [START_REF] Momma | VESTA 3 for three-dimensional visualization of crystal, volumetric and morphology data[END_REF]. This journal is © The Royal Society of Chemistry 20xx

Experimental

Please do not adjust margins Please do not adjust margins

Lattice-dynamic calculations of the phonon modes were carried out at the zone centre ( point) of the BZ with the direct force-constant approach (or supercell method) [76]. These calculations provide not only the frequency of the normal modes, but also their symmetry and their polarization vectors. This allowed us to identify the irreducible representations and the character of the phonon modes at the point.

Conclusions

In this experimental and theoretical approach to the behaviour under compression of -Sb 2 O 3 compound, we have confirmed the 1st order phase transition above 15 GPa of previous works and have proposed that the HP phase of valentinite has an orthorhombic symmetry (s.g. Pcc2) that is a low-symmetry modification of the valentinite structure. The proposed HP phase is supported by ab initio simulations, with good quality factors in the analysis of the powder XRD and with its cationic structure confirmed by single crystal XRD.

Thanks to the analysis of the pressure dependence of Ramanactive mode frequencies in valentinite and its comparison with ab initio calculations, this work provides a solution to the misunderstanding on the irreducible representation of the vibrational modes of -Sb 2 O 3 found in the literature. Moreover, we have presented a tentative assignment of the symmetry of the observed first-order Raman-active modes of -Sb 2 O 3 .

A deep analysis of the compressibility of the theoretically simulated SbO 5 E polyhedral unit, empty linear cavities, and the unit cell volume (F-f plots), as well as the pressure dependence of Sb and O free atomic coordinates reveals changes around 2, 4 and 10 GPa that correlate with changes observed in the vibrational modes of valentinite under compression. These changes were mainly observed in the HP behaviour of the Bader charge and volume around Sb and O2 and they have been interpreted as changes in the orientation and distortion of cation LEPs, respectively. In particular, the changes observed at 2 and 4 GPa can be considered as internal IPTs in valentinite, which are hidden (not evidenced in the bulk compressibility) due to the dominance of the compressibility of the voids of the structure in the bulk compressibility. The transitions at 2 and 4 GPa have been identified with an order higher than 2 because they involve a change of interaction between the neighbouring LEPs (unpaired electrons). These transitions trigger the formation of new bonds, thus increasing the Sb coordination changes from 3 at room pressure to 3+1 above 2 GPa. Hence, they constitute a sequence of transitions from an acicular quasimolecular solid to a 3D solid. Finally, the F-f plots exhibit a change of the compressibility of the whole rod-like structure above 10 GPa that is consistent with the general softening of vibrational modes. We have interpreted this instability caused by the strong interaction between the central O1 atoms and Sb atoms of neighbour Sb 2 O 3 units when the large Sb-O1 distance decreases below 2.38 Å; this leads to an increase in coordination of Sb atoms to 3+1+1 above 10 GPa. Therefore, the coordination of Sb clearly increases to 5 in the high-pressure phase above 15 GPa.

We hope that this work will promote further studies in valentinite for an unambiguous identification of the cause of the changes observed above 10 GPa and the determination of its high pressure phase. We hope that this work, uncovering the interaction between unpaired LEPs in valentinite at 2 GPa, will stimulate further work in the study of subtle pressure-induced phase transitions in molecular and quasi-molecular solids, especially in compounds with cations featuring stereochemically active LEPs, like in group-XV sesquioxides, where a detailed search could lead to the discovery of new and hitherto unexpected internal transitions.

Figure 1 .

 1 Figure 1. (left) 3D structure layout of the low-pressure (LP) phase of -Sb 2 O 3 at ambient pressure (labelling Sb-cyan, O1-red and O2-orange). (right) Projection of the structure in the b-c plane. Lower parts of the figure show the SbO 3 polyhedral units and the two separated chain-like rods within each unit cell extended along the c axis can be observed (LP configuration below 2 GPa). Upper parts of the figure show the SbO 5 polyhedral units so that each rod gets linked with its four neighbor rods (enhanced bonds in the HP configuration above 2 GPa).

  This journal is © The Royal Society of Chemistry 20xx

Figure 4 . 6 |

 46 Figure 4. Pressure dependence of the experimental (solid symbols) and theoretical (empty symbols) volume of the SbO 5 E polyhedral unit (blue) and the empty cavity (red) in -Sb 2 O 3 . The 3 rd order BM EoS fits (solid lines) of theoretical data are also shown.

Figure 5 .

 5 Figure 5. Evolution of the theoretical Bader volumes (top) and charges (bottom) of the non-equivalent atoms in -Sb 2 O 3 . Dashed straight lines are guides to the eyes in order to help visualize the different tendencies observed.

Figure 6 .

 6 Figure 6. Evolution of the LEP-LEP asymmetry (red) and the angle formed by neighboring LEPs and the BCP (blue) under compression. The location of the BCP in the structure is displayed in the inset of the figure for each pressure range proposed.

[ 53 ]

 53 . The Sb volume and charge reveal a change of tendency above 2 and 10 GPa that agrees roughly with the pressure where the compressibility of the SbO 5 E polyhedral volume changes. The O2 volume and charge show changes in tendency between 2 and 4 GPa, which are complementary to that observed in Sb. Finally, the O1 volume and charge follow almost linear trends with a small change of trend near 2 GPa. These results are consistent with the O1 central position in the rod and suggest that O1 atoms are barely involved in This journal is © The Royal Society of Chemistry 20xx

2 O 3 [ 19 , 44 ]

 231944 and also in Sb 2 Te 3 [54] and BiMn 2 O 5 [55],

Table 2

 2

) that confirm the orthorhombic HP phase of β-Sb 2 O 3 as a distortion of its LP phase instead of the monoclinic structure proposed by Zou et al.

[37]

Vibrational properties of -Sb

2 O 3 under compression

  Pressure dependence of the theoretical (lines) and experimental Raman-active mode frequencies of -Sb 2 O 3 obtained from powder (circles) and single crystals (squares) below 400 cm -1 (bottom) and above 400 cm -1 (top). Each colour corresponds to an assigned symmetry: green (A g ), red (B 1g ), black (B 2g ) and blue (B 3g ), whose labels are situated near the theoretical curves. The experimental Raman frequencies of the HP phase observed above 15 GPa (empty triangles) are also shown. vibrational modes of pyramidal SbO 3 (or SbO 3 E) molecules by many authors. On this basis, Cody et al. [49] suggested that internal modes should be similar in the different crystalline structures of Sb 2 O 3 , senarmontite and valentinite; while external modes should depend on the particular space group of the crystalline structure due to the different arrangement of molecules inside the unit cell of both compounds [49]. Later, Miller and Cody used the molecule O 2 Sb-O-SbO 2 to describe the internal modes of vitreous Sb 2 O 3 and assigned internal vibrational modes of this molecule [47]. Finally, a better approach to interpret the RS and IR spectra of valentinite has been performed by Voit et al., who used Sb 10 O 15 units to model the vibrational spectra [49]. Our calculations do not support the separation of vibrational modes of valentinite into external and internal modes of the SbO 3 molecules, as suggested by the previous works, since all molecules are linked together in a complex double chain structure which forms infinite (Sb 2 O 3 ) ∞ chain molecules. Details of the vibrational modes of valentinite and its associated atomic motions are provided in Figs. S11-S13 and S18-S22.

		Raman and infrared (IR) scattering spectroscopies are very sensitive
		tools, much local than XRD measurements, that can show changes in
		the configuration of interatomic bonds. Consequently, they might
		reflect the anomalous compressibility of the LP phase of valentinite
		already commented. According to group theory, the
		centrosymmetric orthorhombic phase of β-Sb 2 O 3 , with D 2h (mmm)
		point group, should have 60 zone-centre phonon modes belonging
		to the following irreducible representation:  = 7A g + 7A u +7B 1g + 7B 1u
		+ 8B 2g + 8B 2u + 8B 3g + 8B 3u .[48] Three of these modes (B 1u , B 2u , and
		B 3u ) correspond to acoustic phonons, while the seven Au modes are
		silent. This results in 30 Raman-active (gerade) optical modes ( R =
		7A g + 7B 1g + 8B 2g + 8B 3g ) and 20 IR-active (ungerade) optical modes
	The RS spectrum of β-Sb 2 O 3 at room pressure is well known	( IR = 6B 1u + 7B 2u + 7B 3u ).
	[35,36,49-52] and has been traditionally divided into two different
	regions separated by a phonon gap: i) the low-frequency region
	(from 50 to 400 cm

This journal is © The Royal Society of Chemistry 20xx J . N a m e ., 2013, 00, 1-3 | 5

Please do not adjust margins Please do not adjust margins Figure 3.

-1 ) and ii) the high-frequency region (from 400 to 800 cm -1 ). The most intense peaks are located in the low-frequency region that is dominated by a peak near 140 cm -1 . Vibrational modes of β-Sb 2 O 3 above (below) 400 cm -1 were assigned to internal (external)

  Ehrenfest's classification[START_REF] Jaegger | The Ehrenfest Classification of Phase Transitions: Introduction and Evolution[END_REF]. New and extended definitions of phase transitions given by Ma and Wang [57] introduces type-I transitions. Such a type of transitions only requires two stable states without variable discontinuity requirement and are called continuous transitions. Unfortunately, the definition of a different initial and final stable state is quite difficult to identify in our case.

61,62]. Partial

  Rietveld refinement was performed by fitting the atomic coordinates of Sb and fixing the atomic coordinates of O1 and O2 to those obtained by ab-initio simulations corresponding to each pressure, in order to obtain a better fit of atomic coordinates of Sb atoms (see Fig.

	Experimental Details
	Synthetic 99.9% pure white powder of -Sb 2 O 3 samples used in
	the present experiments were acquired from Sigma-Aldrich Inc.,
	while fully transparent single crystals of valentinite were
	obtained from Lac Nicolet Mine in Quebec (Canada). Angle-
	dispersive powder HP-XRD measurements at room temperature
	up to 25 GPa were conducted in a membrane-type diamond
	anvil cell (DAC) at the BL04-MSPD beamline of ALBA
	synchrotron in two different experiments. The sample was
	loaded inside DAC using either a MEW mixture (experiment 1)
	or (experiment 2) as a PTM. The pressure inside the
	cavity was calibrated through the EoS of copper [58]. Incident
	monochromatic beam with wavelength of 0.4246 Å was focused
	to 20 x 20 m 2 and cleaned using a pinhole of 50 m [59].
	Images covering a 2θ range up to 20° were collected using a
	SX165 CCD located at 240 mm from sample. Integrated XRD
	patterns profiles as a function of 2θ were obtained with the
	Fit2D software [60]. Le Bail analysis (partial Rietveld
	refinement) of experimental diffraction patterns corresponding
	to experiment 1 (experiment 2) were carried out with GSAS
	package software [

66,67].

  The data were integrated with the program CrysAlisPro[45] which incorporates routines that omit regions of the detector shaded by the diamond-anvil cell from integration. Room-temperature unpolarized HP-RS measurements on powder (up to 22 GPa) and single crystal -Sb 2 O 3 (up to 15 GPa) samples, using MEW as a PTM, were carried out with a Horiba Jobin Yvon LabRAM HR spectrometer equipped with a thermoelectrically cooled multichannel CCD detector that enables a spectral resolution better than 2 cm -1 [63]. RS measurements were excited using the 6328 Å line of a He:Ne laser. Special attention was paid during the experimental measurements to avoid the undesired laser heating effects already observed in the cubic polymorph of the antimony oxide[START_REF] Sui | Laser effects on phase transition for cubic Sb 2 O 3 microcrystals under high pressure[END_REF]. The ruby fluorescence scale[START_REF] Mao | Calibration of the ruby pressure gauge to 800 kbar under quasi-hydrostatic conditions[END_REF] was used to determine the pressure. It must be noted that MEW and helium provide nearly hydrostatic conditions up to 10 GPa and 20 GPa, respectively, and quasi-hydrostatic conditions up to the maximum pressure attained in our experiments [

	Theoretical Details
	Ab-initio total-energy calculations were performed for -Sb 2 O 3
	within the density functional theory (DFT) [68] using the plane-wave
	method and the pseudopotential theory with the Vienna Ab-initio
	Simulation Package (VASP) [69]. The projector-augmented wave
	scheme (PAW) [70-73] was employed and the basis set of plane
	waves was extended up to an energy cutoff of 520 eV in order to
	achieve highly converged results and accurate description of the
	electronic properties. The exchange-correlation energy was taken in
	the generalized gradient approximation (GGA) with the PBEsol [74]

This journal is © The Royal Society of Chemistry 20xx J . N a m e ., 2013, 00, 1-3 | 3

This journal is © The Royal Society of Chemistry 20xx J . N a m e ., 2013, 00, 1-3 | 9

Acknowledgements

This research was partially supported by MINECO (FIS2017-83295-P, MAT2016-75586-C4-1/2/3-P, PGC2018-097520-A-100, RED and PID2019-106383GB-41/42/43), by the MALTA-CONSOLIDER TEAM network (RED2018-102612-T) and by the EU-FEDER. JAS acknowledges Ramon y Cajal fellowship (RYC-2015-17482) for financial support. VM thanks the Juan de la Cierva Program for financial support (FJCI-2016-27921). We also thank ALBA (2016091901 experiment) and DIAMOND light source (CY22223 experiment) for providing beamtime, and Mark Warren, Sarah Barnett and Dominik Daisenberger for his help in the single-crystal X-ray diffraction experiment under pressure.

Author Contributions

Conceptualization has been performed by JAS. Wirting has been performed by JAS and FJM. Investigation on experimental methods was conducting by JAS, FJM, JR, CP, JP, VPCG, JI and VM. Investigation on theoretical calculations was done by AM, PRH, JCG and VM. Methodology was discussed by JAS and FJM.

Conflicts of interest

There are no conflicts to declare.


Nicolas ; Gross
email: nicolas.gross@inrae.fr

Ecodiv Inrae

Yoann Pinguet

Yoann Le Bagousse-Pinguet

Pierre Liancourt

Hugo Saiz

François Munoz

Unveiling ecological assembly rules from trait distributions 5 6

Keywords: dispersal limitation, community assembly rules, environmental filtering, 59 functional diversity, niche differentiation, skewness-kurtosis relationship, stochasticity, trait 60 distributions

Deciphering the effect of neutral and deterministic processes on community assembly is critical 45 to understand and predict diversity patterns. The information held in community trait 46 distributions is commonly assumed as a signature of these processes, but empirical and 47 modelling attempts have most often failed to untangle their confounding, sometimes opposing, 48 impacts. Here, we simulated the assembly of trait distributions through stochastic (dispersal 49 limitation) and/or deterministic scenarios (environmental filtering, niche differentiation). We 50 characterized the shape of trait distributions through the skewness-kurtosis relationship. We 51 identified commonalities in the co-variation between the skewness and the kurtosis of trait 52 distributions with a unique signature for each simulated assembly scenario. Our findings were 53 robust to variation in the composition of regional species pools, dispersal limitation, and 54 environmental conditions. While ecological communities can exhibit a high degree of 55 idiosyncrasy, identification of commonalities across multiple communities can help to unveil 56 ecological assembly rules in real-world ecosystems.

(ranging from -3 to +3) to avoid edge effects (Denelle et al. 2019). We set "m" to either low or 188 high dispersal limitation as explained in the first simulation experiment.

189

The total number of individuals in the external species pool was set to 25,000. This pool 190 delineates the 'regional species pool' in our simulations. To evaluate how the richness of the 191 species pool influenced community assembly, the pool of immigrants varied from 10 to 500 192 species. Richness at the regional level was fixed for each simulation run. The total number of 193 individuals per community was set to 250 (Table S1), and we simulated 100 communities in 194 each run. We chose these parameter values to approximate "realistic" vegetation sampling in 195 the context of dryland (Gross et al. 2017, Table S1). Within communities, the number of 196 individuals per species could vary according to the simulated assembly processes. The trait 197 values of the regional species pool followed a uniform distribution, such that all trait values had 198 equal probability of being selected. We randomly assigned trait values to each species. During 199 this procedure, the model allowed different species to exhibit similar trait values, therefore 4, see also SKR plots in Fig. S1). The neutral and 255 stabilizing scenarios both generated weaker SKRs than the disruptive and directional scenarios 256 (significantly lower R², Fig. 4b). These first two scenarios produced distributions that clustered 257 in the skewness-kurtosis space, and tended to overall converge toward a single skewness-258 kurtosis coordinate as species richness in the regional pool increased (from dark to bright blue 259 dots, Fig. 4a). In contrast, both the directional and the disruptive scenarios yielded strong SKRs, 260 i.e. characterized by significantly higher R² than neutral and stabilizing scenarios (Fig. 4b).

F o r R e v i e w O n l y F o r R e v i e w O n l y 13 287
Furthermore, increasing dispersal limitation marginally influenced the Y-intercept of the SKRs.

288

The mean overlap among scenarios remained extremely low under high dispersal limitation (Y-289 intercept: Mo = 0.03, Fig. 4b). We observed significant differences between the stabilizing, 290 neutral and disruptive scenarios, although the directional scenario overlapped with the neutral 291 one (overlap = 0.18). These results suggested that the SKR approach is therefore able to 292 differentiate a neutral assembly scenario from the effect of environmental filtering (stabilizing 293 scenario) and niche differentiation (disruptive scenario) even under a high stochasticity. Changing environmental conditions across communities affected local mean trait value for all 297 the studied scenarios except the neutral one (Fig. S3). Furthermore, the variance was constant 298 and low for the stabilizing scenario while all moments varied under the disruptive scenario (Fig.

S3

). For the disruptive scenario, increasing the environmental constraint (by changing the Env 300 parameter from 0 to -2 or from 0 to +2) reduced the variance and increased the skewness and 301 the kurtosis (Fig. S3a). Nevertheless, the SKR parameters calculated across communities along

328

The SKR framework shifts our perspective to focus on commonalities in trait 329 distributions -rather than differences -to unveil assembly rules. For instance, the disruptive

398

Second, the distributions of trait values at large spatial scale (e.g. at the scale of the 399 regional species pool) are important drivers determining local assembly processes and trait 400 diversity (Carmona et al. 2016;Le Provost et al. 2017;Spasojevic et al. 2018). To evaluate 401 how regional trait pools influence SKR patterns, we conducted a sensitivity analysis using 402 different trait distributions at the regional scale (see Fig. S5 for details). Our results show that 403 the regional trait pool influenced the SKR parameters for each scenario. It supports the need to 404 account for the regional trait pool to provide reliable predictions on trait diversity within and 405 across communities (Carmona et al. 2016;Spasojevic et al. 2018). Nevertheless, we also 406 observed that the relative differences among contrasting assembly scenarios remained 407 consistent despite variation in the composition of the regional trait pool, highlighting the 408 robustness of the SKR approach.

409

Third, downscaling from a global / regional trait pool to the local community scale may 410 require the development of a spatially-explicit framework, as contrasting assembly processes 586 Zhang, D., Peng, Y., Li, F., Yang, G., Wang, J., Yu, J., et al. (2019) We examined whether the relationship between the skewness and the kurtosis of trait 600 distributions (Fig. 1a) can help deciphering the signatures of contrasting assembly processes.

F o r R e v i e w O n l y F o r R e v i e w O n l y F o r R e v i e w O n l y F o r R e v i e w O n l y F o r R e v i e w O n l y F o r R e v i e w O n l y F o r R e v i e w O n l y F o r R e v i

601

The approach is inspired by optimization procedures increasingly used in physics, climatology 602 and economy (e.g. Cristelli et al. 2012). We apply the approach to diagnose assembly rules

200

  allowing functional redundancy to occur among species(Munoz et al. 2018). For simplicity,201we did not consider intraspecific trait variability in our simulations. Conspecific individuals 202 thus displayed the same trait value. ,000 communities: 500 runs where species richness varied from 10 to 500 206 species in the species pool × 100 communities per run × 4 scenarios × 2 levels of dispersal 207 limitation (low vs. high) × 2 environmental contexts (a fixed environment vs. an environmental 208 gradient). For each run, we recorded the total number of species in a given community (species 209 richness). We quantified species relative abundance as the relative frequency of each species in 210 each community (i.e., the number of individual of a species divided by the total number of 211 individual in the community).

  four moments of the trait distribution associated to each simulated 213 community: the mean, variance, skewness and kurtosis. Then, we calculated the parameters of 214 the SKRs across the 100 simulated distributions obtained from each run (R², slope, Y-intercept 215 and the distance to the lower boundary; see Box 1). For each community, we also computed 216 classical indices of trait diversity using the dbFD function in the "FD" package in R [functional 217 dispersion (FDis), the Rao index, functional evenness (FEve); Laliberté and Legendre 2010]. 218 The four distribution moments and all trait diversity metrics were abundance-weighted. While 219 the SKR parameters were estimated across the 100 communities generated by each simulation 220 run, the four distribution moments, and taxonomic and trait diversity indices were calculated at 221 the community level, and thus generated 100 values for each run. 222 We evaluated whether the parameters of the SKRs, the four distribution moments and 223 all diversity metrics can be used to discriminate assembly processes from trait distributions. For 224 instance, lower variance, FDis and Rao values compared to the neutral scenario are expected 225 for the stabilizing and directional filtering, i.e., trait convergence. Conversely, higher variance, 226 FDis and Rao values compared to the neutral scenario are expected for the disruptive filtering, 227 i.e., trait divergence. To test for significant differences among scenarios, we used the "overlap" 228 package in R (Ridout & Linkie 2009). For each diversity metric, we calculated the overlap 229 between each pair of scenarios and the mean overlap considering all scenarios together. Overlap 230 values ranged from 0 to 1. An overlap < 0.05 indicated that the metric under consideration 231 significantly discriminate assembly scenarios.

235

  Fig.2b). Overall, our results showed limited ability of these metrics of trait diversity to identify

302

  Fig.5b). The R² of the SKRs also differed among scenarios with the stabilizing scenario

  distributions across multiple communities unveil assembly rules 310 Our reference neutral scenario generated a wide range of trait diversity patterns, following the 311 view that "simple stochastic models can reproduce natural diversity patterns" (Hubbell 2005).

  assembly scenarios generated a wide range of trait distributions 313(Enquist et al. 2015; Le Bagousse-Pinguet et al. 2017). The ambiguity among scenarios further 314 increased in our simulations as stochasticity increased or when the environmental conditions 315 changed (e.g., sampling performed along environmental gradients). Therefore, our study 316 therefore demonstrates that distinguishing deterministic from purely neutral patterns based on 317 standard trait-based metrics, or by the four moments of trait distributions taken separately, is 318 particularly challenging. 319 The SKR framework provided far less ambiguous insights into underlying assembly 320 processes, as contrasting assembly scenarios exhibited distinct SKR patterns. The SKR 321 approach did not only well discriminate the neutral from trait-based filtering scenarios, but also 322 clearly differentiated each of the filtering scenarios from one another, irrespective of the degree 323 of dispersal limitation, environmental condition or the size of the species pool. The SKR 324 approach appears therefore as a promising tool to unveil assembly processes from trait 325 distributions, even under high stochasticity, a yet important limitation of existing approaches 326 aiming at interpreting trait diversity patterns of real-world communities (Weiher et al. 2011; 327 Münkemüller et al. 2020).

  330 filtering generated a set of symmetric and asymmetric distributions. Under this scenario, 331 increasing environmental constraints simultaneously reduced the variance and increased both 332 the skewness and the kurtosis of trait distributions (Enquist et al. 2015; Le Bagousse-Pinguet 333 et al. 2017), consistently with the environmental filtering hypothesis (Keddy 1992). However, 334 the trait distributions all aligned along a specific SKR. The observed SKR pattern thus described 335 a family of trait distributions, in which each single distribution represents an instance of a more 336 general trait distribution operating at larger spatial scale (e.g., at the biome scale in Gross et al.

  of skewness. Accounting for the SKR when assessing the impact of global change 363 on ecological communities should provide additional insights on how the ecological reassembly 364 under environmental change, and may provide a more robust validation of theoretical 365 predictions (e.g. the "Trait Driver Theory": Enquist et al. 2015, 2017). Focusing on the distance 366 to the lower boundary of the SKR (see Box 1) (Liu et al. 2020) or the deviation of individual 367 communities from empirical to random SKR (Gross et al. 2017) -rather than raw variations in 368 kurtosis -may help diagnosing shifts in community evenness and community disassembly 369 under ongoing environmental change.

  to empirical patterns of trait diversity 372 Our simulations suggest that niche differentiation and dispersal limitation jointly shape the trait 373 diversity in real-world communities. Combining the disruptive scenario with high dispersal 374 limitation was the only case that generated a strong family of trait distributions (i.e., a SKR 375 with a R² > 0.90). This family of trait distributions exhibited lower kurtosis than under the 376 neutral scenario, a pattern similar to the empirical SKRs observed in drylands worldwide (Gross 377 et al. 2017). When simulating along an environmental gradient, the disruptive scenario also 378 reproduced the impact of environmental filtering by reducing variance and increasing the 379 kurtosis within communities (Fig. S3), and the shift in dominance of contrasting functional 380 groups (e.g. the shrub-to-grass shift in dominance commonly observed in drylands; Fig. S4, 381 Bestelmeyer et al. 2018). Altogether, our results corroborate the view that dryland communities 382 are shaped by a combination of these multiple drivers. 383 The modelling platform ecolottery provides a simple and flexible environment to 384 simulate the effect of multiple stochastic and deterministic processes on trait distributions 385 (Munoz et al. 2018). Further implementations are needed to provide reliable and quantitative assembly processes have been shown to act differently on traits 389 describing contrasted axes of functional specialization (e.g. size vs. resource use related-traits, 390 Kraft et al. 2008; Cornwell & Ackerly 2009; Maire et al. 2012). On the one hand, the disruptive 391 scenario may well apply to resource-use traits as contrasting strategies for resources acquisition 392 and utilization may favor species coexistence (Maire et al. 2012). On the other hand, the 393 directional scenario could constitute a valuable hypothesis for size-related traits in light limited 394 environments where asymmetric competition predominates (Schamp et al. 2007; Gross et al. 395 2009). Considering the effect of multiple assembly processes acting on independent functional 396 dimensions would allow a more realistic representation of community assembly and species 397 coexistence (Maire et al. 2012).

  .J.,Boyle, B., Buzzard, V., Duran, S.M., Henderson, A.N., Hulshof, C.M., et 584 al. (2019). Climate shapes and shifts functional biodiversity in forests worldwide.585Proceedings of the National Academy ofSciences, 116,[587][588][589][590][591][592]

603Fig

  Fig.1b).

Figure 1 660 Figure 2

 16602 Figure 1 Characterizing complex trait distributions using the Skewness-Kurtosis Relationship

669 k 670 Figure 3

 6703 Figure 3 Species richness (Sp. Rich.) and commonly-used trait diversity indices (FDis, Rao,

Figure 4

 4 Figure 4 Effect of the different assembly scenarios on the Skewness-Kurtosis Relationship

Figure 5

 5 Figure 5 Effect of the different assembly scenarios on the Skewness-Kurtosis Relationship

The Skewness-Kurtosis Relationship (SKR) 599

  . Trait identity and

		BOX 1:
	587	functional diversity co-drive response of ecosystem productivity to nitrogen
	588	enrichment. Journal of Ecology, 107, 2402-2414.
	589
		F o
		r
		R
		e
		v i
		e
		w
		O n
		l
		y

Acknowledgments

437

We thank the editor Brian Enquist and the three anonymous reviewers whose 438 comments/suggestions helped improve and clarify this manuscript; Miguel Berdugo, Lucas