Andrew Martin Cassidy

Stine Pedersen

Hendrik Bluhm

Valentin Calisti

Thierry Angot

Eric Salomon

Régis Bisson

Liv Hornekaer

Patterned formation of enolate functional groups on the graphene basal plane

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

Patterned formation of enolate functional groups on the graphene basal plane

Andrew Cassidy, a * Stine Pedersen, a Hendrik Bluhm, b Valentin Calisti, c Thierry Angot, c Eric Salomon, c Régis Bisson, c Liv Hornekaer. a Chemical functionalization of graphene is one method pursued to engineer new properties into a graphene sheet. Graphene oxide is the most commonly used chemical derivative of graphene. Here we present experimental evidence for the formation of enolate moieties when oxygen atoms are added to the graphene basal plane. The exotic functional groups are stabilized by simultaneous bond formation between the graphene sheet and the underlying Ir(111) substrate. Scanning tunneling microscopy images demonstrate the patterned nature of C-O bond formation and x-ray photoelectron spectroscopy and high resolution electron energy loss spectroscopy are used to characterize the enolate moiety. The results present a new mechanism for the formation of patterned graphene oxide and provide evidence of a functional group rarely considered for graphene oxide materials.

Graphene is a 2D semi-metal in which all carbon atoms are sp 2 hybridized. [START_REF] Novoselov | Electric field effect in atomically thin carbon films[END_REF] Chemical functionalization of graphene is driven by the desire to engineer the physical and chemical properties of pristine graphene sheets, and to introduce new means of incorporating graphene interfaces with other materials. [START_REF] Bardhan | Enhanced Cell Capture on Functionalized Graphene Oxide Nanosheets through Oxygen Clustering[END_REF][START_REF] Johns | Atomic covalent functionalization of graphene[END_REF][START_REF] Johns | Metal oxide nanoparticle growth on graphene via chemical activation with atomic oxygen[END_REF] The addition of simple atomic species, such as hydrogen, 5 oxygen 6,7 or fluorine, [START_REF] Chen | Layer-dependent fluorination and doping of graphene via plasma treatment[END_REF] has been shown to change the electronic band structure of graphene from that of a semi-metal to that of a semi-conductor or insulator. Graphene oxide is a term used to describe graphene materials that have been subjected to oxidation reactions. Supported graphene oxide materials have been studied as sensor materials, for their potential magnetic properties and as clusters to form graphene oxide quantum dots. [START_REF] Wen | Magnetism-tuning strategies for graphene oxide based on magnetic oligoacene oxide patches model[END_REF][START_REF] Xie | Reactions of SO 2 and NH 3 with epoxy groups on the surface of graphite oxide powder[END_REF][START_REF] Zhang | Reduction of graphene oxide quantum dots to enhance the yield of reactive oxygen species for photodynamic therapy[END_REF][START_REF] Bardhan | Enhanced Cell Capture on Functionalized Graphene Oxide Nanosheets through Oxygen Clustering[END_REF] The growth of metallic nanoparticles on supported graphene sheets is also facilitated by the presence of C-O bonds at the graphene basal plane, [START_REF] Johns | Metal oxide nanoparticle growth on graphene via chemical activation with atomic oxygen[END_REF] with attachment or cleavage of particle growth initiators dependent on the nature of the C-O bond. [START_REF] Tainter | Reactive Force Field Modeling of Zinc Oxide Nanoparticle Formation[END_REF] This may be a viable method for atomic layer deposition on graphene films, allowing for integration of graphene sheets in more complex electronic devices. [START_REF] Vervuurt | Atomic Layer Deposition for Graphene Device Integration[END_REF] Consequently, the nature of the C-O bond formed during synthesis of graphene oxide is of general interest and importance. It is widely assumed that an oxygen atom bonding to the basal plane of a graphene sheet is initially covalently bound through the formation of an epoxy group. [START_REF] Larciprete | Dual path mechanism in the thermal reduction of graphene oxide[END_REF][START_REF] Vinogradov | Impact of atomic oxygen on the structure of graphene formed on Ir(111) and Pt(111)[END_REF][START_REF] Hossain | Chemically homogeneous and thermally reversible oxidation of epitaxial graphene[END_REF][START_REF] Katano | Localized electronic structures of graphene oxide studied using scanning tunneling microscopy and spectroscopy[END_REF] The stability of different oxygen moieties on metal-supported graphene sheets was, however, recently explored for graphene islands on Ru(0001). [START_REF] Novotny | Formation of Supported Graphene Oxide: Evidence for Enolate Species[END_REF] Those authors, building on theoretical results from Jung et al., [START_REF] Jung | Functionalization of Graphene Grown on Metal Substrate with Atomic Oxygen: Enolate vs Epoxide[END_REF] provided the first experimental evidence for the existence of the enolate moiety when oxygen atoms bind to a graphene sheet supported by a metal substrate. Their calculations indicate that interaction with an underlying metallic substrate can guide the formation of enolate groups when an oxygen atom binds to one carbon atom, and the adjacent carbon atom is suitably positioned to bind down to an underlying metal atom. A 370 meV gain in energy was reported for the formation of an enolate group, versus an epoxy group, in certain areas of the moiré pattern formed between the graphene lattice and atoms from the Ru(0001) surface. [START_REF] Novotny | Formation of Supported Graphene Oxide: Evidence for Enolate Species[END_REF] Graphene on Ru(0001), however, is already a strongly coupled system with good evidence for the preexistence of a carbon-metal bond even before oxygen atoms are introduced. [START_REF] Martoccia | A 25 × 25 supercell[END_REF] Hence, it is timely to consider the nature of the C-O bond for other, less strongly coupled systems. Graphene on Ir(111) (Gr/Ir(111)) demonstrates weak van der Waals interaction between the graphene sheet and the Ir(111) surface atoms. [START_REF] Brako | Graphene on the Ir(111) surface: From van der Waals to strong bonding[END_REF] Moreover, single crystal graphene sheets with a low concentration of defects and highly oriented domains are routinely prepared on Ir(111) substrates following chemical vapour deposition. [START_REF] N'diaye | Structure of epitaxial graphene on Ir(111)[END_REF] The attachment of oxygen atoms to Gr/Ir(111) has been reported previously but the possibility of enolate formation was not considered. [START_REF] Larciprete | Dual path mechanism in the thermal reduction of graphene oxide[END_REF][START_REF] Vinogradov | Impact of atomic oxygen on the structure of graphene formed on Ir(111) and Pt(111)[END_REF] Here, we revisit this system adding new highresolution electron energy loss spectroscopy (HREELS) data to demonstrate that at a low flux of oxygen atoms, enolate formation dominates and occurs only at select sites on the Gr/Ir(111) surface, generating a long-range pattern of graphene oxide nano-dots. The slight lattice mis-match between the carbon atoms in the graphene sheet and the Ir(111) surface atoms leads to a long-range repeating moiré pattern. [START_REF] N'diaye | Structure of epitaxial graphene on Ir(111)[END_REF] The moiré unit cell contains: TOP regions where Ir surface atoms are positioned directly beneath the center of a hexagonal ring in the graphene sheet; and HCP and FCC regions, where every second C atom is positioned directly above an Ir atom. HCP and FCC regions differ in the arrangement of Ir atoms in the second and third Ir layer. The remaining carbon atoms, in between is correspondingly reduced in intensity and has shifted to 284.09 eV. The CC binding energy shift following oxidation is likely induced by doping of the graphene sheet. There is also a contribution from a small third peak, CA at 285.84 eV. The data are consistent with reports on graphene oxide by others. [START_REF] Larciprete | Dual path mechanism in the thermal reduction of graphene oxide[END_REF][START_REF] Vinogradov | Impact of atomic oxygen on the structure of graphene formed on Ir(111) and Pt(111)[END_REF] Here, however, we show a strong preference for O addition as one unique type of C-O moiety, as represented by the CO peak. Photoelectrons from the Ir4f core level, from Gr/Ir(111), show two peaks at 60.80 eV and 60.29 eV, for Ir bulk and surface atoms respectively. [START_REF] Bianchi | Surface core level shifts of clean and oxygen covered Ir(111)[END_REF] Thus the graphene sheet is not chemically interacting with the Ir surface. [START_REF] Larciprete | Dual path mechanism in the thermal reduction of graphene oxide[END_REF] After oxidation, a new peak, at 60.60 eV, is introduced to represent Ir surface atoms that are involved in bond formation with the overlying graphene sheet. The pristine Ir surface peak is correspondingly reduced in intensity. There is no evidence in the Ir spectrum for Ir oxide deposits originating from the thermal cracker. After oxidation the O1s signal is fit, predominantly, with 1 Gaussian peak, OA, with a binding energy of 531.07 eV and a second smaller peak, OB, at 528.83 eV. O1s data for O-Ir(111) show intensity at 530 eV and the absence of such a peak in the O-Gr/Ir(111) data illustrates that no Ir-O bond is formed. [START_REF] Bianchi | Surface core level shifts of clean and oxygen covered Ir(111)[END_REF] Fig. 2c shows the evolution of the relative areas for the C1s peaks, with increasing fluence of O atoms. After 30 min dosing the coverage approximates saturation with an upper limit contribution from the CO peak at 17 % and an upper limit 3% contribution from CA. We propose that the CO peak corresponds both to carbon atoms bound to O-atoms, via formation of an enolate, and to neighbouring carbon atoms simultaneously bound to the underlying Ir substrate. [START_REF] Jørgensen | Symmetry-Driven Band Gap Engineering in Hydrogen Functionalized Graphene[END_REF] These percentages concur with analysis of STM images for surface coverage (S1). The size of such a cluster corresponds to the size of the FCC site in the moiré unit cell. This confirms that the most reactive sites for O atom addition are confined to the FCC region of the moiré. Experiments and DFT calculations illustrate that the position of an Ir atom directly beneath a carbon atom, and C-Ir bond formation, significantly increases the energy gain for C-H bond formation in FCC and HCP regions compared to similar bond formation in TOP regions, with a slight energetic preference for FCC regions. [START_REF] Jørgensen | Symmetry-Driven Band Gap Engineering in Hydrogen Functionalized Graphene[END_REF] The values for O-Gr/Cu(111) are included here. b The C=C) stretching vibration does not originate from phonons as in graphene but appears with similar energy losses. c The region between 120-140 meV is complex for oxidized nano-diamond films and peaks assigned to (C-O-O-C), CH3),C=CH), aromatic ring-H) and (C-O-H) in the literature overlap with the clear peak at 130 meV identified in Fig. 3.

concluded that the same mechanism leads to the stabilization of enolate formation upon formation of C-O bonds for Gr/Cu, Gr/Ni and Gr/Ru. [START_REF] Novotny | Formation of Supported Graphene Oxide: Evidence for Enolate Species[END_REF][START_REF] Jung | Functionalization of Graphene Grown on Metal Substrate with Atomic Oxygen: Enolate vs Epoxide[END_REF] Novotny et al. demonstrated that enolate formation was favoured in FCC regions of the graphene Ru(0001) moiré. We suggest that the CA peak is caused by either a small amount of epoxy group formation [START_REF] Larciprete | Dual path mechanism in the thermal reduction of graphene oxide[END_REF][START_REF] Vinogradov | Impact of atomic oxygen on the structure of graphene formed on Ir(111) and Pt(111)[END_REF] or carbonyl formation, which we justify using HREELS data. 1 for peak identifications.

HREELS data, are shown in Fig. 3. The spectrum for Gr/Ir(111) shows two distinct energy loss peaks; at 104 meV for losses incurred by scattering from the ZO phonons in graphene and a second peak at 199 meV, corresponding to scattering from the LO/TO phonons. [START_REF] Endlich | Screening of electronphonon coupling in graphene on Ir(111)[END_REF] After oxidation a new set of peaks is introduced, Fig. 3b. We interpret this spectrum in Table 1, in light of HREELS for H-Gr/Ir(111) reported by Kyhl et al. [START_REF] Kyhl | Hydrogenated graphene on Ir(111): A high-resolution electron energy loss spectroscopy study of the vibrational spectrum[END_REF] and the calculated vibrations of the enolate group reported by Jung et al. [START_REF] Jung | Functionalization of Graphene Grown on Metal Substrate with Atomic Oxygen: Enolate vs Epoxide[END_REF] The LO/TO and ZO phonons modes for graphene persist and in addition signals arise from scattering loses associated with: the vibration of the newly formed C-Ir bonds, reported in relation to Fig. 2, at 61 meV and 78 meV; LO/TO and ZO phonons associated with the newly formed graphane-like clusters reported in relation to Fig. 1 at 152 meV and 181 meV respectively. Three peaks, at 120 meV, 130 meV and 219 meV, can be associated with vibration of the newly formed C-O bonds and the intensity of these peaks varies with the angle at which scattered electrons are detected. In particular the peak at 130 meV is prominent in the specular direction but reduced in out-of-specular directions, indicating formation of a very polar bond.

Where available, experimental reports on the characterization of oxidized carbon substrates by HREELS are listed in Table 1. [START_REF] Ghodbane | HREELS investigation of the surfaces of nanocrystalline diamond films oxidized by different processes[END_REF][START_REF] Shpilman | Oxidation of diamond films by atomic oxygen: High resolution electron energy loss spectroscopy studies[END_REF][START_REF] Acik | The Role of Oxygen during Thermal Reduction of Graphene Oxide Studied by Infrared Absorption Spectroscopy[END_REF] There is consensus that the energy loss peak at 219 meV can be associated with vibration of a C=O bond indicating that we have formation of carbonyl species. CA at 285.84 eV in XPS data discussed in Fig. 2. Energy losses at 120-135 meV have consistently been associated with vibrations from ether and epoxy groups but the strong polar nature of the peak at 130 meV observed here, and the simultaneous formation of the C-Ir bond, points to the identification of an enolate group. The vibrational energies for the 6 vibrational modes of the enolate moiety calculated by Jung et al. are consistent with this assignment. [START_REF] Jung | Functionalization of Graphene Grown on Metal Substrate with Atomic Oxygen: Enolate vs Epoxide[END_REF] Jung et al. report a substrate dependent vibration between 25-29 meV that is likely hidden by the elastic peak in Fig. 3b. The vibration at 49-58 meV for collective motions of the entire CO group normal to the graphene plane overlaps with C-Ir, reported by Kyhl et al. [START_REF] Kyhl | Hydrogenated graphene on Ir(111): A high-resolution electron energy loss spectroscopy study of the vibrational spectrum[END_REF] We have demonstrated that O atoms covalently bind to the graphene basal plane selectively at sites dictated by the underlying Ir(111) substrate. The spectroscopic evidence supports the view that these new bonds take the form of enolates and microscopy demonstrates the patterned nature of bond formation. Future work will focus on exploiting the potential of using this new functional group to perform further chemistry at the graphene basal plane, and to engineer interfaces with other materials.

Fig. 1 .

 1 Fig. 1. Scanning tunnelling microscopy images of graphene on Ir(111) exposed to a low flux of O atoms at room temperature. The schematic and overlying lattice in the top left identify the repeating dark depressions as "TOP" regions in the moiré pattern. Newly formed C-O bonds appear throughout as circular ring features and are exclusively found at FCC and HCP regions; V=450 mV, I=0.20 nA. The inset shows an atomically resolved area where new C-O bonds have been formed; a moiré lattice has been superimposed; V=-16 mV, I=0.45 nA.

Fig. 2 .

 2 Fig. 2. X-ray photoelectron spectroscopy data for (a) pristine graphene on Ir(111) and (b) oxidized graphene on Ir(111). Individual peak details are given in Table S1 (ESI). cc) Time evolution of the relative area of the C1s peaks with increasing fluence of O atoms.

Figure 3 .

 3 Figure 3. High-resolution electron energy loss spectroscopy data recorded (a) on pristine graphene on Ir(111) and (b) on graphene on Ir(111) exposed to a low flux of O atoms. In both cases curves are shown for data collected in the specular and out of specular directions following irradiation with a 20 eV electron beam. The

Table 1 .

 1 Vibrational features in HREELS data

	Vibration mode		O-Gr/Ir(111) (meV) O-Gr/Cu(111) (meV) 20	H-Gr/Ir(111) (meV) 25	Oxidized diamond films
	Phonon Ir			41		40
	C-Ir			61		52
	C-Ir			78		72
	ZO graphene		104		100	100(C=C)) 26 b
	C-O-C) epoxy		120		57 26 ,125 27 ,135 27 ,140 26 c
	Enolate			130	28,29,54,125,145 a
	LO/TO	phonons	in	152		158
	graphane-like clusters
	ZO phonons in graphane-	181		178
	like clusters
	TO/LO graphene		199		195	200-205(C=C)) 27 b
	(C=O) carbonyl		219		220 26 ,220 27

Jung et al. and Novotny et al., COMMUNICATION https://doi.org/10.1039/C8CP05730C Phys. Chem. Chem. Phys., 2018, 20, 28370-28374 | 28372 a Jung et al. calculated 6 vibrational modes for enolate groups at two distinct sites for free standing graphene, HOPG, Gr/Cu(111) and Gr/Ni(111) substrates.

Acknowledgments

We acknowledge financial support from the Danish Council for Independent Research (grant no. 5137-00049B), Innovation Fund Denmark (NIAGRA) (grant no. 5184-00018B), The European Research Council (CoG GRANN, grant no. 648551), and VILLUM FONDEN via the Centre of Excellence for Dirac Materials (grant no. 11744). Affiliation with the Center for Integrated Materials Research (iMAT) at Aarhus University is gratefully acknowledged. We thank The Advanced Light Source, which is supported by the Director, Office of Science, Office of Basic Energy Sciences of the U.S. Department of Energy at LBNL under contract no. DE-AC02-05CH11231.

Conflicts of interest

There are no conflicts to declare.


Vincent Calcagno
email: vincent.calcagno@inra.fr

Frédéric Hamelin

Ludovic Mailleret

Frédéric Grognard

Andrew Higginson

Frederick Adler

How optimal foragers should respond to habitat changes? On the consequences of habitat conversion

Keywords:

The Marginal Value Theorem (MVT) provides a framework to predict how habitat modifications related to the distribution of resources over patches should impact the realized fitness of individuals and their optimal rate of movement (or patch residence times) across the habitat. Most MVT theory has focused on the consequences of changing the shape of the gain functions in some patches, describing for instance patch enrichment. However an alternative form of habitat modification is habitat conversion, whereby patches are converted from one existing type to another (e.g. closed habitat to open habitat). In such a case the set of gain functions existing in the habitat does not change, only their relative frequencies does. This has received comparatively very little attention in the context of the MVT. Here we analyze mathematically the consequences of habitat conversion under the MVT. We study how realized fitness and the average rate of movement should respond to changes in the frequency distribution of patch-types, and how they should covary. We further compare the response of optimal and non-plastic foragers. We find that the initial pattern of patch-exploitation in a habitat, characterized by the regression slope of patch yields over residence times, can help predict the qualitative responses of fitness and movement rate following habitat conversion. We also find that for some habitat conversion patterns, optimal and non-plastic foragers exhibit qualitatively different responses, and that adaptive foragers can have opposite responses in the early and late phases following habitat conversion. We suggest taking into account behavioral responses may help better understand the ecological consequences of habitat conversion.

Introduction

In most analyses of the Marginal Value Theorem [START_REF] Charnov | Optimal foraging the marginal value theorem[END_REF][START_REF] Stephens | Foraging theory[END_REF], including recent re-analyses [START_REF] Calcagno | The marginal value theorem in a nutshell[END_REF][START_REF] Calcagno | The functional response predicts the effect of resource distribution on the optimal movement rate of consumers[END_REF][START_REF] Calcagno | How optimal foragers should respond to habitat changes: a reanalysis of the marginal value theorem[END_REF], emphasis is on understanding how changes in the shape of gain functions in patches impacted the optimal residence times and movement rate. Yet an alternative form of habitat alteration is to change the relative frequency of different patch categories, without modifying the categories themselves. For instance one might increase the quality of habitat not by making individual patches richer, but rather by making the richest patches more frequent (and the poorest patches rarer). In the context of the MVT, this corresponds to a situation where the travel time and the gain 1/14 functions are unchanged, but the frequency distribution of the different patch categories does vary. Whereas changes in the individual gain functions can readily describe scenarios of habitat enrichment [START_REF] Calcagno | The functional response predicts the effect of resource distribution on the optimal movement rate of consumers[END_REF], modifications of the relative patch frequencies more aptly capture scenarios of habitat conversion (transformation of patches from one category to another). This form of habitat change is illustrated in Fig. 1a. For instance, one might think of converting patches of closed habitat (e.g. forest) into open habitat (e.g. clearings or crop) in a landscape mosaic, or changing the relative frequency of disturbed versus pristine feeding sites [START_REF] Laidlaw | Effects of habitat disturbance and protected areas on mammals of peninsular malaysia[END_REF][START_REF] Presley | Effects of habitat conversion on temporal activity patterns of phyllostomid bats in lowland amazonian rain forest[END_REF][START_REF] Stephens | Predicting risk of habitat conversion in native temperate grasslands[END_REF]. Alternatively, this can describe the extinction of predators from some patches, or their experimental extirpation, turning hazardous places into risk-free areas [START_REF] Kotler | Risk of predation and the structure of desert rodent communities[END_REF]. Habitat conversion is actually one pervasive aspect of the current global biodiversity crisis, impacting many different types of ecosystems [START_REF] Faust | Pathogen spillover during land conversion[END_REF][START_REF] Hoekstra | Confronting a biome crisis: global disparities of habitat loss and protection[END_REF][START_REF] Nowakowski | Thermal biology mediates responses of amphibians and reptiles to habitat modification[END_REF]. Nonetheless, scenarios of habitat conversion have not received much attention in the MVT literature (see e.g. [START_REF] Nonacs | State dependent behavior and the marginal value theorem[END_REF]), probably in part because it requires considering entire habitats, and prevents focusing on individual patches. It further renders classical MVT graphical arguments largely inefficient. Indeed, visualizing optimal residence times as points where the marginal rate of gain equals the long-term average rate of gain (E * n ; [START_REF] Stephens | Foraging theory[END_REF]) is of limited help to predict what happens when patch frequencies are modified, as Fig. 1b illustrates.

Figure 1.

Habitat conversion in the MVT. (a) A patchy habitat with two sorts of patches. Large light disks represent pristine patches, small dark disks disturbed patches. Starting from the situation on the left, with disturbed patches in frequency 1/10, some pristine patches are converted into disturbed patches, resulting in the habitat shown on the right, with 3/10 of disturbed patches. (b) Corresponding changes in residence times under the MVT. It is assumed that pristine patches have more favorable gain functions (light curves) than perturbed patches (dark curves), and the optimal residence times are such that gain functions are tangent to a line of slope E * n (dotted lines). Increasing the frequency of perturbed patches decreases E * n , so that the optimal residence time on both patch-types slightly increase. However, the average residence time decreases, and thus the overall movement rate increases, which is difficult to deduce from graphical intuition.

In a previous mathematical reanalysis of the MVT [START_REF] Calcagno | How optimal foragers should respond to habitat changes: a reanalysis of the marginal value theorem[END_REF], we only briefly mentioned the effect of changing the relative patch frequencies p i (Section 5 in [START_REF] Calcagno | How optimal foragers should respond to habitat changes: a reanalysis of the marginal value theorem[END_REF]). We remarked that following changes in the frequency distribution of patch types, all optimal residence times should vary in the same direction, in opposite direction of the long-term realized rate of gain E * n . However we did not analyze further how the average residence time t * j and the overall rate of movement in the habitat would respond to such changes in the frequency of different patch types. We simply observed that "improving habitat quality by manipulating relative patch frequencies decreases all patch residence times, i.e. increases the movement rate". The last part of the statement is overly simplified. Indeed, the fact that residence time decreases for each individual patch, while true, does not guarantee that the average residence time decreases, since the relative frequencies of 2/14 different patch types was changed. For instance, if residence times are longer on the best patches, then increasing the frequency of the latter mechanically increases the average residence time, and this might counteract the previous change in behavior observed on any particular patch. This is the kind of scenario illustrated in Fig. 1 (transition from left to right).

One would expect that manipulating the resource distribution over patches through changes in the gain functions [START_REF] Calcagno | The marginal value theorem in a nutshell[END_REF][START_REF] Calcagno | The functional response predicts the effect of resource distribution on the optimal movement rate of consumers[END_REF] or through habitat conversion would yield consistent predictions, but this remains to be established. Furthermore, the adaptive response of optimal foragers following habitat change may counteract the direct effects of habitat conversion, and therefore adaptive foragers might in some conditions present qualitatively different responses from static (non-plastic) foragers. To address these points, we here extend our previous reanalyses with a treatment of the effects of changing the relative abundances of different patch-types.

The MVT and notations

We build on [START_REF] Calcagno | The functional response predicts the effect of resource distribution on the optimal movement rate of consumers[END_REF][START_REF] Calcagno | How optimal foragers should respond to habitat changes: a reanalysis of the marginal value theorem[END_REF] and consider the general (heterogeneous) marginal value theorem according to which the optimal residence time t * i on patches of type i is defined implicitly by

   dFi(ti) dti ti=t * i = Fj (t * j ) Tj +t * j i ∈ Ω t * i = 0 i / ∈ Ω (1)
where F i is the gain function in patches of type i, T i is the average time to reach a patch of type i (travel time, usually regarded as the same for all patches) and Ω refers to the set of patches that are effectively exploited in the habitat. Brackets are used to denote spatial averages over the entire habitat (Ω ∪ Ω C ), i.e.

y j = s j=1
p j y j with p j the frequency of patches of type j , for a total of s different patch-types

(j ∈ (1, s)).
Quantity F j (t * j ) / T j + t * j is the long term average rate of gain realized in the habitat (called E * n for short). A well-known consequence of eq. ( 1) is that all exploited patches should be left at the same quitting rate (instantaneous rate of gain at the time the individual leaves) equal to the long term average rate of gain in the habitat E * n .

Following [START_REF] Calcagno | The functional response predicts the effect of resource distribution on the optimal movement rate of consumers[END_REF], we remark that unexploited patches can be gotten rid of by restricting the system to the s patches that are effectively exploited. This implies rescaling the p j to new p j accordingly and increasing the travel times to some effective values Tj , given by

p j = p j k∈Ω p k (2)
and

Tj = T j + k / ∈Ω p k T k k / ∈Ω p k k / ∈Ω p k k∈Ω p k = T j + k / ∈Ω p k T k k∈Ω p k (3)
In these conditions, instead of (1), we can study the slightly simpler equation

dF i (t i ) dt i ti=t * i = F j (t * j ) Tj + t * j (4) 3/14
where there are only s = card(Ω) patch-types to consider, all effectively exploited (i.e. t * i > 0 for all i).

For simplicity we will henceforth restrict our attention to the set of patches that are effectively exploited. We will drop the tildes and ignore Ω, but one should remember that T j , s and p j are intended as their modified values introduced above.

An important quantity to characterize the optimal MVT strategy is the patch-exploitation pattern, quantified by the regression slope (within a given habitat) of patch yields (F j (t * j )) over patch residence times (t * j ). We will call it ρ IN T RA (see [START_REF] Calcagno | The marginal value theorem in a nutshell[END_REF]).

Depending on patch characteristics (shape of the gain functions) and travel time, ρ IN T RA can be positive or negative [START_REF] Calcagno | The marginal value theorem in a nutshell[END_REF][START_REF] Calcagno | How optimal foragers should respond to habitat changes: a reanalysis of the marginal value theorem[END_REF], and it will prove important in predicting the consequences of habitat conversion.

3 Habitat conversion: manipulating the patch frequency distribution

3.1 Which are the best patches?

Differentiating (4) with respect to p i we get

∂E * n ∂p i = 1 T j + t * j 2 T j + t * j F i (t * i ) -F j (t * j ) (T i + t * i )
or, isolating E * n :

∂E * n ∂p i = E * n F i (t * i ) F j (t * j ) - T i + t * i T j + t * j
Of course the total variation in E * n will further depend on how the other patch-types change in frequency; see next Section). Still, requiring the partial derivative to be positive, we get a definition of a "good" (or "better than average") type of patch as one that, if made more frequent in the habitat, would tend (marginally) to increase the realized fitness E * n . This yields the condition

F i (t * i ) T i + t * i > F j (t * j ) T j + t * j = E * n (5)
Note that this criterion to be a good patch is close to the criterion to be effectively exploited (i.e. in Ω). Indeed, only patches such that F i (t * i )/t * i > E * n should be exploited and the others should be entirely disregarded [START_REF] Charnov | Optimal foraging: some theoretical explorations[END_REF]. This condition to be an exploited patch is, as expected, less stringent that the condition to be a good patch (remark that the value of T i does not appear at the denominator, in contrast to eq. ( 5)).

From [START_REF] Charnov | Optimal foraging the marginal value theorem[END_REF], the good patches are characterized by a large F i (t * i ) and/or a small T i + t * i . Two scenarios can thus be distinguished regarding the identity of those patches. First, they can have both shorter exploitation time and greater absolute gains, in which case the two effects described above work jointly, and the best patches are unambiguously the best. Alternatively, the best patches can have greater absolute yields, but longer exploitation time, or smaller yield and shorter exploitation time. In this case the two effects work against one another, and the best patches are not the best in all aspects. In this situation their identity will depend quantitatively on the relationship between absolute yield and exploitation time. Examples of the two situations are shown in Fig. 2. [START_REF] Charnov | Dimensionless invariants from foraging theory's marginal value theorem[END_REF] in green. Three types of patchexploitation patterns are useful to distinguish: (a) Best patches are the best in all aspects: they have both greater absolute yield and shorter residence time. This is detectable as a negative association of patch yield and residence time within the habitat (ρ IN T RA < 0). In (b) and (c) the best patches are not the best in all aspects, corresponding to a positive association of patch yield and residence time within the habitat (ρ IN T RA > 0). In (a) they have greater absolute yield but longer residence time, while in (c) they have shorter residence time but lower yield. In (b) the association is strong, i.e. the regression slope of yield over residence time is steeper than the long-term average rate of gain E * n . The situation in Fig. 1 also fell into this category. In (c) the association is weak, i.e. ρ IN T RA is shallower than E * n . The three situations were obtained with the standard gain function n(1 -exp(-ht)). In (a) patches differ in resource accessibility (h); in (b) they differ in resource content (n); and in (c) they differ in both.

Response of realized fitness to habitat conversion

Any particular pattern of habitat conversion will induce changes in at least two of the p j , while keeping the following constraint satisfied:

s j=1 p j = 1
In order to describe the action of habitat conversion, we introduce a variable x that encodes the specific pattern of conversion (i.e. which patch-types were converted into which). Changing the value of x may alter every patch relative frequency (holding the above constraint satisfied), in a way that describes the type of habitat conversion considered (for example, an increase in the frequency of poor versus good patches). The definition of x can be understood as a direction in the space of the p j along which habitat change operates. We consequently let all p j be functions of x, and the above constraint implies

d dx s j=1 p j = s j=1 dp j dx = 0
The total variation of E * n with respect to x is thus obtained as the sum of the partial effects of each frequency change, i.e.:

dE * n dx = s j=1 dp j dx ∂E * n ∂p j
Considering we are at an optimum with respect to all residence times, we can omit any variation in the t * j from the derivative of E * n (see [START_REF] Calcagno | How optimal foragers should respond to habitat changes: a reanalysis of the marginal value theorem[END_REF]), and it follows

dE * n dx = E * n   1 F j (t * j ) s j=1 dp j dx F j (t * j ) - 1 T j + t * j s j=1 dp j dx (T j + t * j )  
which we may want to rewrite as, in terms of relative variations:

d ln E * n dx = ∂ ln F j (t * j ) ∂x - ∂ ln T j + t * j ∂x
Here, and throughout the rest of the paper, partial derivatives indicate direct variation through the p j , omitting the indirect variation occurring through changes in the optimal residence times, i.e. ∂/∂x = dp j /dx ∂/∂p j .

5/14

Depending on which patch-types are converted into which (i.e. on the habitat conversion pattern x), realized fitness may increase or decrease. One might of course want to define "habitat quality" a priori, irrespective of its actual impact on the realized fitness of an optimal forager. Here we instead take the general (and internally consistent) approach of requiring that habitat quality (at least as perceived by the individual) increases E * n ( [START_REF] Calcagno | How optimal foragers should respond to habitat changes: a reanalysis of the marginal value theorem[END_REF]). Therefore a specific habitat conversion pattern will be said to increase overall habitat quality if and only if it causes an increase in realized fitness.

From the above equation E * n increases with x if and only if

∂ ln F j (t * j ) ∂x > ∂ ln T j + t * j ∂x (6)
Quite intuitively, this shows that habitat conversion can increase the quality of the habitat in two different ways. First, by making more frequent the patch-types that provide the greatest absolute gains (positive left-hand side in ( 6)). Second, by making more frequent the patches that do not take long to exploit (those that are easy to reach and/or with short residence times; negative right-hand side). The net effect of habitat conversion on realized fitness rests on the balance of these two effects.

Here, it is important to stress that the variation in realized fitness of an optimal forager is not influenced by the plastic response of residence times (see also [START_REF] Calcagno | How optimal foragers should respond to habitat changes: a reanalysis of the marginal value theorem[END_REF]). It follows that our definitions of patch and habitat quality, and thus the response of fitness to habitat conversion, are equally valid for optimal (perfectly plastic) or static (non-plastic) foragers, i.e. those that would keep their patch-exploitation pattern unchanged following habitat conversion.

Comparing eq. ( 6) with the condition for E * n to increase when manipulating patch characteristics (equation ( 6) in [START_REF] Calcagno | How optimal foragers should respond to habitat changes: a reanalysis of the marginal value theorem[END_REF]), we see that it is almost identical: the first term is the same, and the second is similar, only lacking the penalty term T j / T j + t * j . This is because when varying T j as a patch attribute, we only operate on a fraction of the entire time it takes to exploit one patch (while the optimal residence time itself was out of external control). Here, in contrast, when converting patches from one type to another, the entire duration of the exploitation cycle (travel and residence time) is impacted, and the penalty factor therefore vanishes to unity.

Remark that since our sensitivity analysis approach considers infinitesimal changes in the p j , conclusions are unaffected by the implicit treatment of unexploited patches (see [START_REF] Calcagno | The marginal value theorem in a nutshell[END_REF] and (3) above). Indeed, infinitesimal changes generically do not change Ω, and the latter can be regarded as a constant. Of course, to predict the consequences of sustained changes (i.e. to integrate over x), one should remember to update Ω appropriately, when a patch-type leaves or enters the set of exploited patches.

Response of average movement rate to habitat conversion

We now consider the effect of habitat conversion on the optimal average movement rate. In the context of the MVT, the latter is defined as the rate of patch switching 1/ T j + t * j , i.e. the inverse of the average time needed to exploit one patch, following [START_REF] Calcagno | The functional response predicts the effect of resource distribution on the optimal movement rate of consumers[END_REF]. For a systematic (directional) forager, this rate of patch-switching controls the linear speed through the habitat [START_REF] Arditi | Optimal foraging on arbitrary food distributions and the definition of habitat patches[END_REF], and for a random-walker it sets the timescale of the exploration and is as such expected to be proportional to the diffusion coefficient [START_REF] Turchin | Quantitative analysis of movement[END_REF]. In what follows we will thus study the variation of T j + t * j . This quantity is closely related to another quantity of interest in the field of behavioral ecology, the average residence time t * j .

Obviously, the total variation in T j + t * j can be split into two components. First, there is the direct effect of changing the relative frequency of the different patches, 6/14 which mechanically alters the average. This effect is

∂ T j + t * j ∂x = s j=1 dp j dx T j + t * j
Second, there is an indirect effect whereby the resulting change in E * n (computed in the previous section) impacts all individual optimal residence times t * j (but not travel times, obviously). It follows that

d T j + t * j dx = s j=1 dp j dx T j + t * j + s j=1 p j dt * j dx
Replacing the total variations of optimal residence times with their expression (see [START_REF] Calcagno | How optimal foragers should respond to habitat changes: a reanalysis of the marginal value theorem[END_REF]), some calculations (Appendix) yield the following criterion for the average movement rate to increase with x:

∂ ln F j (t * j ) ∂x > ∂ ln T j + t * j ∂x 1 + |H| F j (t * j ) (7)
where, as in [START_REF] Calcagno | How optimal foragers should respond to habitat changes: a reanalysis of the marginal value theorem[END_REF], H is the harmonic mean of the second time derivatives of gain functions (H < 0). This is very similar to eq [START_REF] Charnov | Dimensionless invariants from foraging theory's marginal value theorem[END_REF], except that the relative change in patch exploitation time (right-hand side) is multiplied by a constant greater than one (parenthesis). The latter depends on the shape of the gain functions (height and second time-derivative).

Therefore the impact of the change in patch exploitation time is amplified compared to eq. ( 6).

The previous result holds true for an optimal forager that has quickly enough adjusted its movement strategy following habitat conversion. If, on the extreme opposite, we consider a non-plastic forager (or, similarly, a plastic forager before it had time to update its strategy), the observed indirect variation of optimal residence times (dt * j /dx) would be nil, and the criterion for average movement rate to increase is then simply:

0 > ∂ ln T j + t * j ∂x (8)
Figure 3. Predicting the impacts of habitat conversion: a summary. The graph summarizes mathematical predictions on the effect of habitat conversion on realized fitness and on average movement rate, for an optimal forager (top) or a non-plastic forager (bottom). The plane shown represents all possible scenarios of habitat conversion, as characterized by ∂ ln T j + t * j /∂x (x-axis) and ∂ ln F j (t * j ) /∂x (y-axis). Three areas can be identified that yield qualitatively different predictions (legend on the left). Note that fitness and movement rate can covary negatively only in the gray dotted area, that is much smaller for an optimal forager. The pattern of patch-exploitation prior to habitat conversion (see Figure 2) can be used to determine in which part of the plane we are in (see legend). Circled points locate a specific scenario of habitat conversion that will be simulated numerically in Fig. 4.

Discussion

Consequences of habitat conversion for fitness and movement rate

Combining ( 6) and ( 7), we obtain a simple graphical representation of what impact any pattern of habitat conversion would have on fitness and movement rate for an optimal forager, and from ( 6) and ( 8) the corresponding for a non-plastic forager (Figure 3). All conversion patterns can be sufficiently characterized with a pair of quantities: the relative change in average patch-exploitation time (∂ ln T j + t * j /∂x) and the relative change in patch yield (∂ ln F j (t * j ) /∂x). This reveals some constraints on the possible patterns of co-variation of fitness and movement rate following habitat conversion.

First, we can see that for all types of foragers, plastic or non-plastic, habitat conversion patterns that decrease habitat quality (realized fitness) will necessarily decrease average movement rate as well (gray area). Conversely, if habitat conversion results in higher average movement rate, then it should increase overall fitness as well 8/14 (white dotted area). In other words, if fitness decreases or average movement rate increases, fitness and movement necessarily co-vary positively. However, there is the possibility of a negative covariation for some habitat conversion patterns (gray dotted area).

Importantly, the likelihood of a negative co-variation of fitness and movement is much more important for a non-plastic forager compared to an optimal forager (the white dotted area in larger in Fig. 3). For optimal foragers, a negative covariation of fitness and movement rate only occurs for a narrow set of habitat conversion patterns (gray dotted area in Fig. 3), feasible only if

|H| F j (t * j )
is large enough. In graphical terms, this means that the gain functions should be sufficiently curved (concave), and consistently so (H is a harmonic mean, very sensitive to low values), relative to their height.

As a consequence, there is an entire set of habitat conversion patterns, comprised between the lines of eq. ( 7) and eq. ( 8) in Figure 3, for which optimal foragers and non-plastic foragers will exhibit qualitatively different responses. This difference between optimal and non-plastic foragers would also manifest itself as a difference between the initial and final responses of an optimal forager. Indeed, following habitat conversion, an adaptive foraging that gradually (rather than instantaneously) adjusts its patch-exploitation pattern [START_REF] Mcnamara | Optimal foraging and learning[END_REF], would initially exhibit a negative covariation of fitness and movement rate, but ultimately reverse to a positive covariation, once it has fully updated its strategy. Numerical examples of these situations (see Appendix for simulation methods) are provided in Figure 4.

Predicting the consequences of habitat conversion from the initial pattern of patch exploitation

We can obtain further predictions under the reasonable and common assumption that all patches have on average the same travel time T , or equivalently that travel times show no consistent response to habitat conversion (see Appendix). This occurs for instance if travel times are controlled by the spatial location of patches. In these conditions optimal movement rate is entirely governed by average optimal residence time, and the two are just inversely related [START_REF] Calcagno | The functional response predicts the effect of resource distribution on the optimal movement rate of consumers[END_REF].

In this case, the partial derivative on the right-hand side of eq. ( 6), ( 7) and ( 8) is entirely controlled by the residence times on converted patches. It is therefore possible to discriminate different portions of Figure 3 (i.e. different possible habitat conversion patterns) from the observed pattern of patch exploitation in the (unperturbed) habitat (Fig. 2). In practice, one must compute the regression slope of patch absolute yields

(F i (t * i )) over residence times (t * i ) in the initial habitat, that we called ρ IN T RA ( [2]).
This regression slope provides the generic value of

E * n ∂ ln F j (t * j ) /∂ ln T j + t * j
following habitat conversion (an exact value with two patch-types, and a value representative of most possible habitat conversion patterns in the general case). It follows that the three contrasted patch-exploitation patterns shown in Figure 2, discriminable from the regression slope ρ IN T RA , generically produce habitat conversion patterns that map onto different parts of Figure 3.

On Fig. 3, one can see that patch-exploitation pattern (a), i.e. a negative regression of patch yield over residence time, ensures a positive covariation of fitness and movement following habitat conversion, for optimal or non-plastic foragers alike. The same is true for patch-exploitation pattern (c), except that changes will have the opposite sign. However, for patch-exploitation pattern (b) (strong positive regression 9/14 slope), whilst a negative co-variation is ensured for non-plastic foragers, no firm prediction can be achieved for an optimal forager (Fig. 3). Indeed, prediction further requires knowledge of which patches were converted, and of the curvature of fitness functions (parameter H in eq. 7), to tell apart the gray and white dotted areas. Furthermore, it is for patch-exploitation pattern (b) that one would expect qualitative differences between optimal and non-plastic foragers, or non-monotonous dynamical responses from adaptive foragers (see Figs. 3 &4). Figure 4. A simulated example of habitat conversion. In a two patch-type habitat (pristine versus perturbed patches; same gain functions as in Fig. 1 & 2), between time 0 and time 1 (shaded parts) the frequency of pristine patches is gradually increased from 1/5 to 3/4. The variation of fitness (left) and average movement rate (right) are shown through time for two foragers: an optimal forager that gradually adjusts its patch-exploitation pattern to the new habitat conditions, and a non-plastic forager, that maintains the same patch exploitation pattern. Note that the two foragers initially present the same qualitative responses (negative covariation of fitness and movement), but ultimately qualitatively different responses (the optimal forager adopting greater movement rate). In our simulation habitat conversion caused a 61% increase in F j (t * j ) and a 3% increase in T j + t * j . As an indication these values were located as a circle on Fig. 3; simulation results conform to mathematical predictions.

Conclusions

We envisioned habitat change in the Marginal Value Theorem from the perspective of habitat conversion, i.e. changes in the relative frequencies of the different categories of patches. This differs from the classical approach of varying the shape of the gain functions [START_REF] Calcagno | The functional response predicts the effect of resource distribution on the optimal movement rate of consumers[END_REF][START_REF] Calcagno | How optimal foragers should respond to habitat changes: a reanalysis of the marginal value theorem[END_REF]. The predictions obtained are in good agreement with previous re-analyses of the MVT [START_REF] Calcagno | The marginal value theorem in a nutshell[END_REF][START_REF] Calcagno | How optimal foragers should respond to habitat changes: a reanalysis of the marginal value theorem[END_REF]. For example:

• An increase in habitat quality may have any (negative, null or positive) impact on average movement rate (or average residence time).

• An increase in average movement rate with habitat quality requires less stringent conditions [START_REF] Calcagno | The marginal value theorem in a nutshell[END_REF][START_REF] Calcagno | How optimal foragers should respond to habitat changes: a reanalysis of the marginal value theorem[END_REF]. While a negative correlation of patch yield and residence time ensures that an increase in habitat quality will increase average movement (patch-exploitation pattern (a); leftmost halves of Fig. 3), a positive correlation does not guarantee a decrease in movement rate.

• For an optimal forager, if there is a strong positive correlation of patch yield and residence time (pattern (b)), one cannot predict the response of average movement to habitat quality based solely on simple average quantities (average patch yield, average residence time, average rate of gain E * n ). We further need to know the harmonic mean of the gain function second-time derivatives (|H|; see eq. ( 7)). The greater |H|, i.e. the more curved the gain functions, the more likely it is to observe a decrease in movement with habitat quality. Consistent results were obtained with the other form of habitat change (see Theorem 3 and eq. ( 17) in [START_REF] Calcagno | How optimal foragers should respond to habitat changes: a reanalysis of the marginal value theorem[END_REF]).

10/14

It thus seems one can safely transpose results obtained with one approach to the other. We can expect existing predictions regarding the consequences of changing the resource distribution in a habitat that were derived from changing the gain functions [2-4, 6, 21] to extend to scenarios of habitat conversion as well. This also suggests one might obtain an integrated mathematical treatment of the consequences of habitat changes under the MVT.

In this work we also obtained novel predictions that can guide our prediction and interpretation of actual responses to habitat conversion. For instance, the pattern of patch-exploitation in the pristine habitat, that can be characterized from a simple regression line of of patch yields over patch residence times (ρ IN T RA ; Fig. 2; see also [START_REF] Calcagno | The marginal value theorem in a nutshell[END_REF]), can be used to predict what type of responses to expect following habitat conversion (Fig. 3). Furthermore, only some types of habitat conversion patterns, occurring for patch-exploitation pattern (b), are expected to generate qualitatively different responses for plastic and non-plastic foragers. Last, adaptive foragers could for some habitat conversion patterns exhibit qualitatively different responses in the short and in the long-term (Fig. 4), making it challenging to extrapolate ultimate consequences from initial responses.

While the theoretical consequences of habitat conversion are usually considered at the scale of populations and communities [START_REF] Mouquet | Extending the concept of keystone species to communities and ecosystems[END_REF][START_REF] Tilman | Habitat destruction and the extinction debt[END_REF], we suggest the MVT can offer a framework to apprehend these questions at the individual and behavioral scales. For instance, dynamic changes in the movement strategies of individuals on the short-term may change the relative exploitation of different parts of the habitat, as well as the realized level of migration and connectivity. Depending on their sign, behavioral responses might mitigate, or on the contrary amplify, the population and community responses expected under the assumption of static migration rates and connectivity levels. A better integration of the behavioral and population scales might therefore improve our ability to predict the ecological consequences of habitat change in general, and habitat conversion in particular.

n dx 1 H
where H is the harmonic mean of the second time-derivatives, i.e.

H = 1/ s j=1 p j d 2 F j (t * j ) dt 2 j -1
Replacing the variation of E * n with its expression (given in Section 3.

This, replacing E * n with its expression, yields eq. ( 7). Now, if all patches have on average the same travel time T , or if variation in travel time shows no consistent trend with habitat conversion, we have d T j /dx = 0. Note that under the second scenario (i.e. if all travel times are not equal), achieving a null derivative would in practice requires specific forms of habitat conversion, owing to the constraint that the sum of the p j is constant. The constraint gradually vanishes as the number of patch-types (s) gets large. From d T j /dx = 0 it follows:

∂ T j + t * j ∂x = ∂ t * j ∂x

Numerical simulations

In order to generate the numerical simulations presented in Figure ( 4), we used a simple gradient ascent algorithm. Individuals were assumed to update each residence time (t j , omitting the asterisk as they need not be at optimal value) gradually, in the direction that (locally) increases the long-term average rate of gain (E * n ) and at a rate proportional to the fitness differential, i.e.

13/14 dt j dt = ω ∂E * n ∂t j [START_REF] Kotler | Risk of predation and the structure of desert rodent communities[END_REF] where t denotes ecological time (on which habitat changes take place) and ω is a constant quantifying the speed of behavioral adjustments. Habitat conversion is modeled by specifying how the patch relative frequencies change through time, i.e. by specifying functions p j (t). In the simulations of Fig. 4, the function were taken to be linear, so that relative frequencies changed linearly from their initial to their final values. It was assumed that individuals, prior to the onset of habitat change, had settled at the optimal residence times for the initial habitat. When ω is very large the forager is effectively optimal and immediately adjusts its strategy to match current habitat conditions (in accordance with the MVT). When ω = 0 the individual does not adjust its strategy (non-plastic forager). Intermediate values of ω represent less-than-perfect plastic foragers that gradually adapt to habitat changes.

14/14

Figure 2 .

 2 Figure 2. Patch-exploitation patterns within habitats. Hypothetical habitats are shown, with the best patch type (sensu eq.[START_REF] Charnov | Dimensionless invariants from foraging theory's marginal value theorem[END_REF] in green. Three types of patchexploitation patterns are useful to distinguish: (a) Best patches are the best in all aspects: they have both greater absolute yield and shorter residence time. This is detectable as a negative association of patch yield and residence time within the habitat (ρ IN T RA < 0). In (b) and (c) the best patches are not the best in all aspects, corresponding to a positive association of patch yield and residence time within the habitat (ρ IN T RA > 0). In (a) they have greater absolute yield but longer residence time, while in (c) they have shorter residence time but lower yield. In (b) the association is strong, i.e. the regression slope of yield over residence time is steeper than the long-term average rate of gain E * n . The situation in Fig.1also fell into this category. In (c) the association is weak, i.e. ρ IN T RA is shallower than E * n . The three situations were obtained with the standard gain function n(1 -exp(-ht)). In (a) patches differ in resource accessibility (h); in (b) they differ in resource content (n); and in (c) they differ in both.

  We can directly identify the partial log derivative of F j (t * j ) and T j + t * j in the parenthesis, and upon factoring E * n / |H| out we obtain:Introducing the partial log derivative of T j + t * j , and requiring the parenthesis to be negative (since average movement rate varies in opposite direction of T j + t * j ) yields

														2 in the main text),
	we get
	d T j + t * j dx	=	s j=1	dp j dx	T j + t * j +	E * n H	 	1 F j (t * j )	s j=1	dp j dx	F j (t * j ) -	1 T j + t * j	s j=1	dp j dx	(T j + t * j )	 
	d t * j dx		=	E * n |H| T ∂ ln F j (t *   -∂ ln F j (t * j ) ∂x + |H| E * n s j=1 dp j dx j ) ∂x > ∂ ln T j + t * j ∂x ∂ T j + t * j + ∂x	|H| E * n

j + t * j + ∂ ln T j + t * j ∂x  

Acknowledgments

This work was supported by INRA and Université Côte d'Azur (IDEX JEDI). This preprint has been reviewed and recommended by Peer Community In Ecology (https://dx.doi.org/10.24072/pci.ecology.10005)

Conflict of interest disclosure

The authors of this preprint declare that they have no financial conflict of interest with the content of this article.

Appendix

Derivation of eq. [START_REF] Charnov | Optimal foraging: some theoretical explorations[END_REF] As explained in Section 3.3 in the main text, and following [START_REF] Calcagno | How optimal foragers should respond to habitat changes: a reanalysis of the marginal value theorem[END_REF] to express the indirect variation of t * j , we can formulate the total variation of T j + t * j as :


domain: phys chim

text:
Arthur Marronnier
email: arthur.marronnier@polytechnique.edu

Heejae Lee

Bernard Geffroy

Jacky Even

Yvan Bonnassieux

Guido Roma

Structural Instabilities related to Highly Anharmonic Phonons in Halide Perovskites

Hybrid organic-inorganic perovskites have emerged over the past five years as absorber layers for novel high-efficiency low-cost solar cells that combine the advantages of organic and inorganic semiconductors. Unfortunately, the physics behind the electronic transport in these materials is still poorly understood. Here, employing the linear response approach of Density Functional Theory, we reveal strong anharmonic effects in inorganic CsP bI 3 and a double-well phonon instability at the center of the Brillouin zone for both cubic and orthorhombic phases. Factoring these low-energy phonons into electron-phonon interactions and band gap calculations could help better understand the electrical transport properties in perovskites. Furthermore, previously reported soft phonon modes are here stabilized at the actual lower-symmetry equilibrium structure, which occur in a very flat energy landscape, highlighting the strong competition between the different phases of CsP bI 3 . The perovskite oscillations through the corresponding energy barrier could explain the underlying ferroelectricity and the dynamical Rashba effect predicted in halide perovskites for photovoltaics.

Hybrid Perovskite Solar Cells (PSCs), a new generation of solar cells which combines the advantages of organic and inorganic semiconductors, were born only a few years ago after perovskites were used as absorbing materials in dye-sensitized solar cells (DSCs). [START_REF] Kojima | Organometal Halide Perovskites as Visible-Light Sensitizers for Photovoltaic Cells[END_REF] After a mere 6 years of research, the efficiency of PSCs has jumped from 4% in 2009 to certified efficiencies over 20% 2 in 2015 and a record efficiency of 22.1% 3 in 2016, becoming the first truly low-cost and highly efficient generation of solar cells.

Researchers started phonon calculations of hybrid M AP bI 3 perovskites as early as 2013, when Quarti et al. [START_REF] Quarti | The Raman Spectrum of the CH 3 NH 3 PbI 3 Hybrid Perovskite: Interplay of Theory and Experiment[END_REF] measured its Raman spectrum, assisted by Density Functional Theory (DFT) calculations which showed the importance of the torsional mode of the methylammonium cations as a marker of the orientational disorder of the material. As shown by various groups since then, [START_REF] Leguy | The Dynamics of Methylammonium Ions in Hybrid Organic-Inorganic Perovskite Solar Cells[END_REF][START_REF] Chen | Rotational Dynamics of Organic Cations in the CH 3 NH 3 PbI 3 Perovskite[END_REF][START_REF] Even | Molecular Disorder and Translation/Rotation Coupling in the Plastic Crystal Phase of Hybrid Perovskites[END_REF] this picture is further complicated by the reorientational dynamics of the organic component similar to what is usually observed in plastic crystals, since quasielastic neutron scattering measurements showed that dipolar CH 3 N H + 3 ions reorientate between the faces quite rapidly. Very recently, Whalley et al. [START_REF] Whalley | Phonon Anharmonicity, Lifetimes, and Thermal Transport in CH 3 NH 3 PbI 3 from Many-body Perturbation Theory[END_REF] and Beecher et al. [START_REF] Beecher | Direct Observation of Dynamic Symmetry Breaking above Room Temperature in Methylammonium Lead Iodide Perovskite[END_REF] found double-well instabilities present at the Brillouin Zone (BZ) boundary of M AP bI 3 , much like what we report here in the CsP bI 3 case, as well as short phonon quasiparticle lifetimes and mean free paths. However, contrary to our findings for CsP bI 3 , they did not find any instabilities at the center of the BZ for M AP bI 3 . These structural effects in organohalide perovskites were studied in detail by Amat et al. [START_REF] Amat | Cation-Induced Band-Gap Tuning in Organohalide Perovskites: Interplay of Spin-Orbit Coupling and Octahedra Tilting[END_REF] who showed that the interplay of spin-orbit coupling (SOC)

and octahedra tilting can have a strong impact on their electronic and optical properties, including the band gap.

It is only in the past two years that inorganic (mostly cesium-based) perovskites for photovoltaics have attracted more and more attention from the PSC community after being incorporated into complex but more stable mixed-halide mixed-cation perovskite structures (such as [HC(N H 2 ) 2 ] 0.83 Cs 0.17 P b(I 0.6 Br0.4) 3 ) demonstrating the feasibility of achieving above 25% efficiency four-terminal tandem cells 11 at the beginning of 2016.

In particular, the vibrational properties of inorganic and lead-free perovskite CsSnI 3 have been studied in detail using ab-initio calculations. Lora da Silva et al. employed [START_REF] Da Silva | Phase Stability and Transformations in the Halide Perovskite CsSnI 3[END_REF] the quasiharmonic approximation to study the temperature-dependent lattice dynamics of the four different phases of cesium tin iodide. Patrick et al. reported [START_REF] Patrick | Anharmonic Stabilization and Band Gap Renormalization in the Perovskite CsSnI 3[END_REF] strong anharmonic effects in CsSnI 3 (including soft phonon modes) using self-consistent phonon theory and frozen phonons, and showed that these soft modes were stabilized at experimental conditions through anharmonic phonon-phonon interactions between the Cs ions and their iodine cages.

They further calculated the renormalization of the band gap due to vibrations. These two works demonstrate the important role of temperature in accurately analysing anharmonic perovskite materials like CsSnI 3 .

The role of phonons in the phase transitions of perovskites containing lead and cesium was recognized long ago (e.g. CsP bCl 3 14 ). This suggests that vibrational properties could control the structural transitions also for CsP bI 3 , which is crucially important for most recent PSCs.

The vibrational properties of inorganic CsP bI 3 perovskites have not been extensively studied so far. Stoumpos et al. studied the dependence of CsP bI 3 's band gap on temperature and the reversibility of its phase changes from non-perovskite structures to perovskite structures at high temperature. They showed [START_REF] Stoumpos | The Renaissance of Halide Perovskites and Their Evolution as Emerging Semiconductors[END_REF] that the room temperature stable δ-phase (yellow) converts to the black perovskite α-phase upon heating above 360 • C. On cooling, the perovskite structure converts to the black perovskite β-and γ-phases at 260 and 175 • C, respectively. After a few days, a full recovery of the initial δ-phase is observed. A prediction of possible structural instabilities in CsP bI 3 and RbP bI 3 perovskite was made, on the basis of relatively high values of the Born effective charges, by Brgoch and co-workers [START_REF] Brgoch | Ab Initio Calculations of Band Gaps and Absolute Band Positions of Polymorphs of RbPbI 3 and CsPbI 3 : Implications for Main-Group Halide Perovskite Photovoltaics[END_REF] ; but to the best of our knowledge it has remained a speculation until we present our results here.

In this paper, we study the vibrational properties of both the high-temperature cubic α-phase and the orthorhombic δ-phase of CsP bI 3 , and compare them to the phonons of M AP bI 3 's cubic phase. For both phases, phonon spectra, derived using the linear approach of Density Functional Theory (Density Functional Perturbation Theory or DFPT), showed the existence of instabilities which we investigated in detail using total energy calculations along the corresponding phonon eigenvectors (frozen phonons framework). More information about the DFT and DFPT calculations can be found in the Computational Methods section at the end of the paper. Additional calculations on the influence of spin-orbit coupling on these results are also detailed in the Supporting Information file.

The phonon dispersion relations of inorganic perovskite CsP bI 3 in the cubic high temperature α-phase, with commonly accepted space group P m3m, [START_REF] Trots | High-temperature Structural Evolution of Caesium and Rubidium Triiodoplumbates[END_REF] are reported in figure 1a. The presence of soft modes at different high-symmetry points in the Brillouin zone can be explained by the fact that the α-phase is stable only at high temperature, where both energetic and entropic contributions determine the free energy F, and not at the (zero) temperature of DFT calculations. In principle, these modes are stabilized by taking into account temperature, i.e., the vibrational entropy and its contribution to the free energy.

What is more surprising is the presence of imaginary modes at the Γ-point (one triply degenerated mode at -15 cm -1 ), where one expects to find acoustic modes, linearly going to zero at Γ, as the lowest frequency modes. These instabilities and in particular this anomaly at Γ for the cubic phase of CsP bI 3 have already been reported by Kawai et al. [START_REF] Kawai | The Mechanism of Slow Hot-Hole Cooling in Lead-Iodide Perovskite: First-Principles Calculation on Carrier Lifetime from Electron-Phonon Interaction[END_REF] (see figure 1b), but they offered no explanation or further analysis.

As circled in red on figure 1a, the Acoustic Sum Rule (ASR) was imposed in our calculations through post-processing of the dynamical matrix to ensure translational invariance.

The ASR did bring back to zero the three-fold degenerated acoustic branches at the Γ-point but did not affect the soft mode, confirming its nature as an optical mode. We analyzed the dynamical matrix's eigenvector corresponding to this soft mode, and its anti-symmetrical characteristics (see figure 2) are consistent with the irreducible representations of the P m3m symmetry group as reported by Even et al. [START_REF] Even | Molecular Disorder and Translation/Rotation Coupling in the Plastic Crystal Phase of Hybrid Perovskites[END_REF] The soft modes obtained here are coherent with previous neutron scattering measurements clearly indicating that the energy of the whole acoustic phonon density of states is remarkably low at 80 • C. [START_REF] Even | Molecular Disorder and Translation/Rotation Coupling in the Plastic Crystal Phase of Hybrid Perovskites[END_REF] In particular, it was shown for CsP bCl 3 that the transverse acous- tic (TA) branches always stay under 2 meV (16 cm -1 ) 14 and that the optical phonons close to the R and M points also lie in the same energy range (< 2 meV). Here, our calculations show instabilities of the same order of magnitude (around 20 cm -1 ) coming from optical modes close to the M point (see figure 1a). Similar conclusions were drawn more recently for CH 3 N H 3 P bI 3 9 and CH 3 N H 3 P bBr 3 . [START_REF] Létoublon | Elastic Constants, Optical Phonons, and Molecular Relaxations in the High Temperature Plastic Phase of the CH 3 NH 3 PbBr 3 Hybrid Perovskite[END_REF] We think that these instabilities indicate strong anharmonic effects (also confirmed by the unusually large Debye-Waller factors measured previously for M AP bX 3 compounds [START_REF] Even | Molecular Disorder and Translation/Rotation Coupling in the Plastic Crystal Phase of Hybrid Perovskites[END_REF][START_REF] Baikie | A combined Single Crystal Neutron/Xray Diffraction and Solid-State Nuclear Magnetic Resonance Study of the Hybrid Perovskites CH 3 NH 3 PbX 3 (X = I, Br and Cl)[END_REF] )

and a very flat energy profile around equilibrium in CsP bI 3 . This very flat profile could affect the reported low-to-negative frequency values in two ways: (i) by affecting the evaluation of the frequency values itself (due to numerical noise), and (ii) by preventing the structural relaxation algorithm to find the ground state structure. In order to assess these two hypotheses, we first (i) used smaller and smaller DFPT convergence thresholds and second (ii) performed frozen phonon calculations to investigate anharmonicity up to the fourth order.

The use of very tight thresholds for total energy convergence, forces during relaxations, and self-consistency of the linear response calculations (see Computational Methods section) reduced the number of imaginary phonons at Γ, but one triply degenerated soft mode remained (see figure 1a), compared to 5 as published by Kawai et al. [START_REF] Kawai | The Mechanism of Slow Hot-Hole Cooling in Lead-Iodide Perovskite: First-Principles Calculation on Carrier Lifetime from Electron-Phonon Interaction[END_REF] (figure 1b).

Unfortunately the most widely-used approach to go further than the harmonic approximation and to calculate the free energy F by adding a volume (thus temperature) dependent vibrational free energy, the quasi-harmonic approximation, 21 cannot be straightforwardly applied [START_REF] Da Silva | Phase Stability and Transformations in the Halide Perovskite CsSnI 3[END_REF] due to the presence of these unstable phonon modes. Another technique is to guess the unstable directions using phonon eigenvectors obtained with the linear response approach of DFT. Here, we computed frozen phonon-perturbed total energies for the soft optical mode at Γ. A displacement of the form:

u κ (η) = η • x κν Γ (q Γ ) (1)
was imposed on each ion κ, starting from the relaxed positions, and for various values of the unitless displacement parameter η, where x is the chosen soft mode eigenvector obtained from harmonic DFPT phonon calculations.

Using Laudau theory [START_REF] Cowley | Structural Phase Transitions I. Landau Theory[END_REF] and the adiabatic approximation, which assumes the electronic system instantaneously changes after new ionic positions are reached, the anharmonic energy [START_REF] Fujii | Neutron-scattering Study on Phase Transitions of CsPbCl 3[END_REF]23 as a function of displacement is given by: The results, shown in figure 3, fit well with this fourth-order model, confirming anharmonic effects in cubic CsP bI 3 and the hypothesis (ii) that we were not at the equilibrium structure. Given the very small height of the barrier (7.3 meV), we can assume that the perovskite structure actually oscillates between both equilibrium positions. To get an idea of the order of magnitude of these oscillations, one can write:

U (η) = U 0 + η 2 4 hω ν Γ (q Γ ) + η 4 24 Φ 4,ν Γ + O(η 6 ) (2)
τ = τ 0 e E k B T (3)
where k B is the Boltzmann constant and E the energy barrier. Even if we consider, for τ 0 , the period of the lowest frequency optical phonon (ω ≈ 17 cm -1 ), this gives an oscillation timescale at room temperature of around 0.7 • 10 -12 s. These frozen phonon calculations are performed following the soft phonon mode and thus the new minimum here is not relaxed in volume. After volume relaxation, the energy barrier reads 8.4 meV.

The soft mode (infrared-active Γ - 4 representation, figure 2) that corresponds to the highsymmetry high-energy structure is a polar soft mode induced by the displacement of positively charged cesium ions in one direction, and the displacement of negatively charged iodine ions in the opposite direction. The oscillations of the perovskite along this polar soft mode and between the two low-symmetry equilibrium structures could explain why ferroelectricity has not been observed at a macroscopic scale in these perovskites, considering that a ferroelectric behavior is more likely induced at the scale of these fast fluctuations. Similarly, it was shown recently that for hybrid methylamonium-based perovskites, the coupled inorganic-organic degrees of freedom give rise to a spatially local and dynamical Rashba effect which fluctuates on the subpicosecond time scale typical of the methylammonium cation dynamics. [START_REF] Etienne | Dynamical Origin of the Rashba Effect in Organohalide Lead Perovskites: A Key to Suppressed Carrier Recombination in Perovskite Solar Cells?[END_REF] We think that the oscillations predicted here in the cubic and orthorhombic phases of CsP bI 3 could also be at the origin of a similar dynamical Rashba effect. Additional calculations on the dynamical Rashba effect induced by this symmetry breaking can be found in the Supporting Information file.

Performing new phonon calculations at the new (volume relaxed) equilibrium positions definitively removed the remaining soft mode (see figure 4). In conclusion, both (i) using smaller convergence thresholds and (ii) finding the actual minimum using frozen phonon calculations allowed us to investigate and abolish the soft modes at Γ for the cubic phase of CsP bI 3 . Similar phonon dispersions are obtained when taking into account SOC, as shown in the Supporting Information file, although caution must be observed due to the reduced band gap.

This means that this strongly anharmonic mode at Γ will not condensate at lower temperatures; it does not correspond to the transition to the yellow phase. On the contrary, the does not correspond to a simple distortion of the reference cubic perovskite lattice. [START_REF] Stoumpos | Semiconducting Tin and Lead Iodide Perovskites with Organic Cations: Phase Transitions, High Mobilities, and Near-Infrared Photoluminescent Properties[END_REF] We relaxed the two structures and found a slightly lower energy for the second structure. As the difference in energy was very small, we interpolated a few configurations between them to check for a barrier, and, to our surprise, we found an even lower minimum (see figure 5), with a relatively flat energy landscape. The conclusions are the same when taking into account spin-orbit coupling and lead semicore electrons. This is another signature of a possible anharmonic behavior in CsP bI 3 . A further confirmation of this fact is that, even for the lowest energy structure of figure 5, we again found a few soft modes at the Γ-point, as shown in figure 6. at other high symmetry points. Similar to our work with the cubic phase, we followed the eigenvectors of the soft modes in order to explore the energy landscape. As with the cubic phase, frozen phonon calculations around the lowest energy structure of figure 5 were fitted up to the 4th order and showed (see figure 7) a similar double-well instability with a higher potential barrier (124 meV, meaning an oscillation time of around 0.6 • 10 -10 s at room temperature) and once more a good agreement between frozen calculations and DFPT on the second order coefficient. We recall that these frozen phonon calculations are performed following a soft phonon mode and thus that the new minimum here is not relaxed in volume.

After volume relaxation, the energy barrier reads the much higher value of 2.39 eV. We think that this very significant difference could be the sign of a strong ferroelasticity in the δ-phase of CsP bI 3 .

As in the cubic case, phonons were re-calculated at the (volume-relaxed) new equilibrium positions and all the soft modes at Γ were eliminated (see figure 8). This stable structure has a lower symmetry than the P nma space group and, as such, does not match with the black γ-phase described by Stoumpos, [START_REF] Stoumpos | The Renaissance of Halide Perovskites and Their Evolution as Emerging Semiconductors[END_REF] nor with the structure of the γ-phase of CsSnI 3 . [START_REF] Yamada | Structural phase transitions of the polymorphs of CsSnI 3 by means of rietveld analysis of the X-ray diffraction[END_REF] We then compared the vibrational properties of the inorganic CsP bI 3 to those of M AP bI 3 .

The phonon spectrum of the pseudocubic phase of the hybrid perovskite M AP bI 3 , obtained similarly using DFPT calculations, did not show any soft modes at Γ (figure 9). This spectrum also differs significantly as it contains high-frequency optical modes around 3000 cm -1

which are the signature of the M A + organic molecule vibrations.

This could suggest that M AP bI 3 behaves noticeably more harmonically than CsP bI 3 .

However, as shown very recently, [START_REF] Whalley | Phonon Anharmonicity, Lifetimes, and Thermal Transport in CH 3 NH 3 PbI 3 from Many-body Perturbation Theory[END_REF][START_REF] Beecher | Direct Observation of Dynamic Symmetry Breaking above Room Temperature in Methylammonium Lead Iodide Perovskite[END_REF] this is due to the fact that symmetry breaking is facilitated by the presence of the organic molecule, but a similar behavior, with a double well potential around the most symmetric structure, is also present in M AP bI 3 . This double well instability was found at the edge of the BZ (M and R points) and not at the Γ point.

The barrier was reported to be somewhat higher, between 0.02 and 0.04 eV. Furthermore, we think that our simulation here is not able to capture the dynamical translation-rotation coupling (recently experimentally observed [START_REF] Létoublon | Elastic Constants, Optical Phonons, and Molecular Relaxations in the High Temperature Plastic Phase of the CH 3 NH 3 PbBr 3 Hybrid Perovskite[END_REF] ) which may yield further contribution to the anharmonicity of the phonon modes in M AP bI 3 .

To summarize, through energy landscape and phonon calculations for the inorganic perovskite CsP bI 3 we reveal that the assumed equilibrium structure of the high-temperature α-phase, of space group P m3m, results from a dynamical average between lower symmetry configurations in a very flat energy landscape.

For the orthorhombic low-temperature phase, similarly, we find a rather flat energy landscape between the two experimentally proposed variants of the yellow δ-phase, a lower symmetry equilibrium configuration and double well potentials with energy barriers on the order As for the hybrid M AP bI 3 perovskite, translation-rotation coupling of the organic molecule needs to be considered in order to fully investigate its anharmonicity.

These findings highlight the strong competition between the yellow δ-phase and the black α/β/γ-phases, and suggest that the cation network suffers from dynamical disorder even at relatively low temperature. The unrealized ferroelectricity in these perovskites could be replaced by the polar structural fluctuations we report here, also leading to a dynamical Rashba effect as previously proposed for similar perovskite structures. [START_REF] Etienne | Dynamical Origin of the Rashba Effect in Organohalide Lead Perovskites: A Key to Suppressed Carrier Recombination in Perovskite Solar Cells?[END_REF][START_REF] Even | DFT and k • p modelling of the phase transitions of lead and tin halide perovskites for photovoltaic cells[END_REF][START_REF] Azarhoosh | Research Update: Relativistic origin of slow electron-hole recombination in hybrid halide perovskite solar cells[END_REF] Since SOC almost closes the electronic band gap for the α-phase, the calculation of phonon frequencies close to the metallic transition might well be affected by numerical instabilities (see Supporting Information). The effect of spin orbit coupling on the vibrational properties of halide perovskites, both organic and inorganic, will certainly be of high interest in future studies.

Further work is needed to take into account these low-energy highly occupied phonon states in order to give new theoretical estimates of electronic and optical properties (e.g.

the band gap or the exciton binding energy [START_REF] Bokdam | Role of Polar Phonons in the Photo Excited State of Metal Halide Perovskites[END_REF] ) and to improve electron-phonon coupling calculations. Indeed, electron-phonon coupling was shown to deeply impact the chargecarrier mobilities [START_REF] Wright | Electron phonon coupling in hybrid lead halide perovskites[END_REF][START_REF] Sendner | Optical phonons in methylammonium lead halide perovskites and implications for charge transport[END_REF] in these perovskites. Mixing Cs with other cations, like the organic ones used in most recent versions of hybrid PSCs, affects the coupling in the double-well network, which could help explain the mechanisms through which mixed-cation PSCs have shown interesting improvements in efficiency and stability with respect to M AP bI 3 -based perovskite solar cells.

Computational Methods

Electronic-structure calculations were performed within the Density-Functional Theory (DFT) 33 34 framework, as implemented in the Quantum Espresso code. [START_REF] Giannozzi | Quantum Espresso: a Modular and Open-source Software Project for Quantum Simulations of Materials[END_REF] All total energy and force calculations in this work were performed with the Local Density Approximation (LDA), expanding the wave functions in a plane-wave basis set. Unless otherwise specified, non-relativistic (scalar-relativistic for Pb) and norm-conserving pseudopotentials were used, with the Cs [5s 2 5p 6 6s 1 ], I [5s 2 5p 5 ] and Pb [5d 10 6s 2 6p 2 ] electrons treated as valence states. We also performed additional calculations with fully relativistic US pseudopotentials (for Pb) and PAW (Cs and I) datasets (with the same number of va-lence electrons as in the scalar case) in order to check the influence of spin-orbit coupling (see results in the Supporting Information file). These relativistic pseudopotentials were generated using the PS Library 0.2.2. [START_REF] Corso | Pseudopotentials periodic table: From H to Pu[END_REF] Plane-wave cutoffs of 70 Ry (CsP bI 3 ) and 80 Ry (M AP bI 3 ) were used. The Brillouin zone (BZ) was sampled with Γ-centered Monkhorst-Pack meshes [START_REF] Monkhorst | Special Points for Brillouin-zone Integrations[END_REF] with subdivisions of α-CsP bI 3 : 8x8x8; δ-CsP bI 3 : 10x5x3 and M AP bI 3 : 8x8x8.

Phonon calculations were performed using the linear response approach of DFT, as known as DFPT (Density-Functional Perturbation Theory), [START_REF] Baroni | Phonons and Related Crystal Properties from Density-Functional Perturbation Theory[END_REF] as implemented in the Quantum Espresso code. [START_REF] Giannozzi | Quantum Espresso: a Modular and Open-source Software Project for Quantum Simulations of Materials[END_REF] High-frequency dielectric tensors and Born effective charges were calculated with the linear response and used to calculate long range contributions to the dynamical matrices at Γ. All phonon dispersions were performed using at least 5 q-points/ Å-1 in each direction.

In order to eliminate some of the soft phonon modes, very tight convergence thresholds of 10 -4 Ry/bohr for the force calculations and 10 -14 for the phonon self-consistent algorithm were used.

( a )Figure 1 :

 a1 Figure 1: Phonon modes of cubic CsP bI 3 .

Figure 2 :

 2 Figure 2: Displacement eigenvector corresponding to the soft phonon mode found at Γ for cubic CsP bI 3 . Its irreducible representation is the infrared-active Γ - 4 representation. The blue, purple and green atoms respectively denote Cs, Pb and I.

Figure 3 :

 3 Figure 3: Potential-energy surface from frozen phonon calculations of cubic CsP bI 3 along the eigenvector of the unstable optical phonon at Γ as a function of displacement parameter η. The 3N dimensional displacement needed to reach the new minimum corresponds to around 0.43 Å. ω f rozenphonons and Φ 4 are obtained by fitting equation 2. As expected, the results show good agreement between frozen phonon and DFPT estimations of the second order coefficient (ω).

  remaining phonon instabilities at the M and R points at the edge of the Brillouin zone are related to soft modes which condensate at lower temperatures to yield the β tetragonal and γ orthorhombic black phases of CsP bI 3 .[START_REF] Even | Molecular Disorder and Translation/Rotation Coupling in the Plastic Crystal Phase of Hybrid Perovskites[END_REF][START_REF] Stoumpos | The Renaissance of Halide Perovskites and Their Evolution as Emerging Semiconductors[END_REF]

Figure 4 :

 4 Figure 4: Phonon modes of cubic CsP bI 3 at the new equilibrium position determined in figure 3 after volume relaxation. Using both tight convergence thresholds and the results of our frozen phonon calculations allowed us to remove the soft modes at Γ.

  in 2008 by powder diffraction studies. This structure is referred to as the yellow δ-phase and

Figure 5 :

 5 Figure 5: Equilibrium energy of CsP bI3's orthorhombic δ-phase (T<530K). The 5 intermediate position sets are linearly interpolated between sets 1 ("Moeller") 25 and 7 ("Trots") 17 which are equilibrium positions relaxed from experimental initial structures. The distance is given in 3N dimensions.

Figure 6 :

 6 Figure 6: Phonon modes of orthorhombic δ-CsP bI 3 in the lowest energy structure of figure 5.

Figure 7 :

 7 Figure 7: Potential energy surface from frozen phonon energy calculations of orthorhombic δ-CsP bI 3 along the eigenvector of the lowest soft phonon at Γ as a function of displacement parameter η. The 3N dimensional displacement needed to reach the new minimum corresponds to around 2.2 Å. ω f rozenphonons and Φ 4 are obtained by fitting equation 2. The results show good agreement between frozen calculations and DFPT on the second order coefficient (ω).

Figure 8 :

 8 Figure 8: Phonon modes of orthorhombic δ-CsP bI 3 at the new equilibrium position determined in figure 7 after volume relaxation. Using both tight convergence thresholds and the results of our frozen phonon calculations allowed us to remove the soft modes at Γ.

Figure 9 :

 9 Figure 9: Phonon frequencies of pseudocubic M AP bI 3 .

(

  LICSEN). He graduated from CNAM (Engineer in physical chemistry), Paris 1987. His current research interests are OLEDs, OPVs (design, fabrication and characterization) and organic/inorganic hybrid photovoltaics (SiNWs, ZnO nanostructures, hybrid perovskite). He is the author or co-author of 120 publications in international scientific journals, 1 chapter book and 5 patents. He is an Editor of The European Physical Journal Applied Physics (EPJAP). Yvan Bonnassieux is currently a Professor at Ecole polytechnique (Paris-Saclay University, France) and deputy manager of the LPICM laboratory (laboratory of the interface and thin Layers physics). Former student of the ENS of Cachan, he received is PhD in 1998 and HDR in 2011. His current research themes include compact modelling, design and characterization of components and circuits of organic and large area electronics. He was also an invited professor at Columbia University in 2015. He is the author or joint author of more than 80 publications in scientific magazines, 120 communications in international conferences and 4 patents. Jacky Even has been full Professor at INSA Rennes engineering school since 1999. He received a PhD in Physics from the University of Paris VI in 1992. He was assistant professor at the Physics Department of Rennes University (1992-1999), where he gained broad experience in phase transitions and solid-state chemical reactions in molecular materials, using theoretical and experimental approaches: neutron and X-ray scattering, Raman/FTIR spectroscopy, calorimetry, among others. In 1999, he created FOTON laboratory's simulation team, to address fundamental questions on semiconductors at the atomistic level as well as to perform optoelectronic device simulations. Besides hybrid perovskite materials and colloidal nanoplatelets, his theoretical activity is now dedicated to semiconductor nanostructures, photovoltaic and light-emitting (LED, laser, ...) devices for silicon photonics, and optical telecommunication. Guido Roma is a research-engineer at Service de Recherches de Métallurgie Physique (SRMP), CEA-Saclay (Paris-Saclay University, France). After graduating from the University of Modena and receiving a PhD at the University of Rome Tor Vergata (Italy), he spent four years at CECAM (Centre Européen de Calcul Atomique et Moléculaire, Lyon, France) before joining SRMP in 2000. He holds a habilitation from the University of Paris VI (2012). His main experience is on point defects thermodynamics and kinetics, starting from a first principles description of electronic and vibrational properties, with a focus on insulators and semiconductors. His interests range from irradiation effects on materials for nuclear applications to functional materials for energy, in particular photovoltaics, a field in which he got involved during a stay at the University of Mainz (Germany) between 2011 and 2014.

  for Renewable Energies. Arthur earned a PhD grant from the French Department of Energy (IPEF, state engineer) and started his PhD in the LPICM lab of Ecole Polytechnique in 2015. His current research is focused on anharmonicity, instabilities and oxygen defects in halide perovskites for photovoltaics. Heejae Lee is a Ph.D. candidate at Paris-Saclay University within Ecole Polytechnique Graduate School. He joined LPICM and Prof. Bonnassieux's group after completing a dual master course between Kyung-Hee University (KHU), Korea and Ecole Polytechnique, France. Prof. Jin Jang and Prof. Bonnassieux were the supervisors of his master course.

	His current research is focused on the electrical transport and aging characteristics of hybrid
	perovskite thin films.

Bernard Geffroy is currently a research engineer at CEA Saclay (Paris-Saclay University, France) in the Laboratory of Innovation in Surface Chemistry and Nanosciences

Acknowledgement

Arthur Marronnier's PhD project is funded by the French Department of Energy (MEEM) and by the Graduate School of École des Ponts ParisTech. This work was granted access to the HPC resources of TGCC and CINES under the allocation 2016090642 made by GENCI. Jacky Even's work has been performed within the GOTSOLAR FETOPEN project, which has received funding from the European Union's Horizon 2020 research and innovation program under grant agreement No 687008. The information and views set out in this work are those of the authors and do not necessarily reflect the official opinion of the European Union. Neither the European Union institutions and bodies nor any person acting on their behalf may be held responsible for the use which may be made of the information contained herein.

Author Information

Notes

The authors declare no competing financial interest.

Biographies

Supporting Information Available

The following files are available free of charge.

• SI.pdf: This Supporting Information file is devoted to presenting some checks that we made in order to show that some of the approximations used in this paper do not affect its conclusions. In particular, most of the results presented in the paper have been obtained by neglecting spin-orbit coupling.


Angélique Villière

Fabrice Guillet

Sarah Le Roy

Catherine Fillonneau

Carole Prost

Carole Prost Introduction

la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction of an innovative olfactometry-dedicated software using intuitive aroma wheel interface

Introduction of an olfactometry-dedicated software using intuitive aroma wheel interface

1 ONIRIS -College of Veterinary Medicine and Food Science, Flavor group -GEPEA UMR CNRS 6144, Nantes, France 2 Polytech Nantes, LINA CNRS 6241, Nantes, France  A software with a good usability for judges, overcoming bias of existing systems, providing a complete characterization and including data treatment capabilities 2 .

 A wheel aroma presentation, consistent with those currently used for sensory characterization.

Angélique VILLIÈRE 1 , Fabrice GUILLET 2 , Sarah LE ROY 1 , Catherine FILLONNEAU 1 , Carole PROST 1

Olfactometry is a valuable methodology used to investigate odorant active compounds in complex aroma profiles. Considering the numerous studies using this technique, little is done to improve olfactometric data acquisition. An innovative software dedicated to olfactometric data acquisition has been developed and is herein presented (patent pending 1 ).  Facilitates chemometric approaches to determine contribution of compounds to aromas and apprehend existing interactions.

  Time angelique.villiere@oniris-nantes.fr  Good usability for judges  No oral transmission  No disruption of respiratory rhythm  Continuous evaluation of perceptions Data processing

FERMENTATIVE

Table of results

 of  Direct visualization of individual and panel aromagrams,  Access to detection and intensity data besides odor characterization  Automatic conversion from time to Linear Retention Indices  Ready-to-use concatenated data  Comparison of products, judges potentialities

	N° zone	Apex Retention time (min)	LRI	START time (min)	END time (min)	LRI START	LRI END	Judge 1	Judge 2	Judge 3	Judge 4	Judge 5	Judge 6 Judge 7	Judge 8
	27	8,97	1199	8,93	9,22	1197 1212		amylic	UNKNOWN					red fruits
	28	9,28	1216	9,23	9,43	1213 1225 CHEMICAL	feet	feet	feet	FERMENTATIV E	feet	feet	feet
	29	9,73	1241	9,68	9,85	1239 1248	FRUITY	exotic fruits				FRUITY		exotic fruits
	30	10,00	1256	9,95	10,05	1254 1259		FRUITY
	31	10,36	1277	10,33	10,38	1275 1278			FLORAL

1 Prost C., Bertrand D. 2010. Procédé d'analyse d'une impression sensorielle, patent n°FR 2933788, July 9th 2008 2 Villière A., Le Roy S., et al. 2015. Evaluation of aroma profile differences between sué, sautéed and pan fried onions using an innovative olfactometric approach. Flavour journal 4(24): 2-19.