Discrete breathers for understanding reconstructive mineral processes at low temperatures

arXiv:nlin/0404030v4 [nlin.PS] 27 Sep 2006 Discrete breathers for understanding reconstructive mineral processes at low temperatures JFR Archilla∗ , J Cuevas Grupo de F´ısica No Lineal. Universidad de Sevilla. Departamento de F´ısica Aplicada I. ETSI Inform´atica Avda. Reina Mercedes, s/n. 41012-Sevilla, Spain MD Alba, M Naranjo and JM Trillo Departamento de Qu´ımica Inorg´anica. Universidad de Sevilla. Instituto de Ciencia de Materiales de Sevilla. Consejo Superior de Investigaciones Cient´ ificas. P.O. Box 874, 41080-Sevilla (Spain) September 20, 2006 Abstract Reconstructive transformations in layered silicates need a high temperature in order to be observed. However, very recently, some systems have been found where transformation can be studied at temperatures 600◦C below the lowest experimental results previously reported, including sol-gel methods. We explore the possible relation with the existence of intrinsic localized modes, known as discrete breathers. We construct a model for nonlinear vibrations within the cation layer, obtain their parameters and calculate them numerically, obtaining their energies. Their statistics shows that although there are far less breathers than phonons, there are much more above the activation energy, being therefore a good candidate to explain the reconstructive transformations at low temperature. Keywords: discrete breathers; reconstructive transformations; intrinsic localized modes PACS: 63.20.Pw, 63.20.Ry, 63.50.+x, 66.90.+r, 82.20.-w ∗Corresponding author. Email: [email protected] 1 1 Introduction During the last decade, some of the present authors have achieved the synthesis of crystalline high-temperature polymorphs of rare earth (RE) disilicates (RE2Si2O7) at non-expected temperatures, significantly lower than those previously reported, through a reconstructive process (LTRT) from clay minerals as the silicon source.1–4 This finding is of general importance in the development of advanced structural ceramics5or the storage of radioactive wastes3and should allow completion of the available Si O2-RE2O3phase diagram.6Although no precise explanation has been found up to now, some of the present authors had previously suggested a chimie douce mechanism based on the diffusion of RE ions into the interlayer space of the expandable clay minerals.3,4 MacKay and Aubry7have suggested that a possible effect of the existence of localized nonlinear vibrations, named discrete breathers (DBs) could be an apparent violation of Arrhenius’ law, i.e., the phenomenon of chemical reactions taking place at much lower temperatures than expected. Although this hypothesis is adventurous, it is worth exploring. Moreover, experimental evidence of DBs has already been found in several systems such as antiferromagnets,8waveguide arrays,9molecular crystals10 or Josephsonjunctions.11 Moving DBs have been proposed as an explanation of dark tracks in muscovite12 and there existence in muscovite has just been proven through a sputtering experiment.13 With this aim, we have made calculations and shown that the contribution of DBs can provide an interpretation for LTRT in clay minerals. And in order to give experimental support to the hypothesis of DBs, it will be shown that the LTRT phenomenon is not exclusive of expandable clay minerals, as expected by the previously suggested ”chimie douce” mechanism,3,4 but also extensible to non-expandable layered silicates, such as mica muscovite. The layout of this work consists of: Section 2: Some structural considerations on the reconstructive nature of transformation of layered silicates into disilicates; Section 3: A report on a new experiment on LTRT performed by the authors on mica muscovite; Section 4: Argumentation about the difficulty to explain by the conventional chemical kinetics model the latter experiment; Section 5: Description of an alternative model based on DBs with numerical calculations; Section 6: After a summary of breather statistics theory, the description of our numerical simulations and the consequences on the reaction rate, we conclude with the possibility of explaining, for the first time, the LTRT phenomenon in the synthesis of high-temperature polymorphs of silicates by the contribution of DBs. The article in itself is ended 2 with a summary. Two appendices give some detail on phonon and breather statistics respectively. 2 From layered silicates to disilicate crystal structures The synthesis of disilicates from layered silicates, expandable as clay minerals or non-expandable as mica, actually means a reconstructive transformation as shown below. Layered silicates are made up from two basic building blocks: a sheet of edge-sharing [SiO4] units, the tetrahedral sheet, and another one of edgesharing [MO6], the octahedral sheet. There are three main groups of layered silicate minerals, according to the combinations of tetrahedral and octahedral sheets: 1:1, 2:1 and 2:1:1. In 2:1, one octahedral sheet is sandwiched between the apices of two tetrahedral ones. In these, also called T-O-T silicates, layers are either held together by weak van der Waals forces if they are neutral, or may have cations between them for charge balance if substitutions in either tetrahedral or octahedral sheet result in a net layer charge. The 2:1 layered silicates are classified as trioctahedral or dioctahedral, after the full occupation of the octahedral sheet by Mg(II) or two thirds by Al(III). Talc (trioctahedral) and pyrophillite (dioctahedral) are minerals with non-charged layers. In the case of low charge, it results that the clay minerals have the capacity to expand by taking up H2O molecules in the interlayer space. For high charge, there is mica: phlogopite (trioctahedral) and muscovite (dioctahedral). Muscovite is a mica which layer charge comes from the isomorphic substitution of silicon by aluminium in the tetrahedral sheet.14 The potassium located in the interlayer space, for charge balance, cannot be hydrated; thus, muscovite does not expand. Its structure is depicted at the left of Fig. 1. At its right, the interlayer space of muscovite is illustrated and it can be observed that both surfaces of the upper and lower tetrahedral sheet are formed by the basal oxygen atoms from the [SiO4] tetrahedra, which form a rough hexagonal honeycomb structure. The interlayer balancing sheet is therefore sandwiched in between, potassium occupying the dimples left at the centre of each pair of hexagonal cells. In real crystals, although the onsite potential created by the silicate layer tends to preserve the symmetry, there are always distortions like tetrahedral rotation (the actual situation in muscovite is shown in this figure). Each potassium ion is surrounded by six other potassium ions in the interlayer sheet with a ditrigonal symmetry. 3 Figure 1: Crystal unit cell of muscovite ICSD 34406. The circles represent the potassium ions forming the interlayer sheet. (a = 5.19 ˚ A; b = 9.02 ˚ A; c = 20.0˚ A; β=95.7◦) The distance between the potassium and the basal oxygens layers is 1.45 ˚ A. Essentially, in both clay minerals and mica [SiO4] tetrahedra are linked into infinite two-dimensional networks by sharing three oxygens. However, disilicates, or pyrosilicates, are the simplest of the condensed forms, where only two tetrahedra share one edge and constitutes the anion Si2O76− (Fig 2). The transformation of any layered-silicate into disilicate involves the rupture of two silicon-oxygen bonds by each [SiO4], whatever the reaction mechanism might be, the transformation being reconstructive. It is well-known that natural pyrosilicates show a wide range of Si-O-Si angle, from 131◦to 180◦,15 and that the activation energy is reduced if the surface and strain energy terms are diminished by good lattice matching across the interface between the new and parent phase. However, this is not the present case. It rather seems that the disruption of the tetrahedral sheet could be the consequence of localized nonlinear vibration modes as commented in the introduction and explained in the following sections. 3 Experimental RE-disilicates synthesis The method used by us to synthesize RE-disilicates consists of a hydrothermal reconstructive process at low temperatures in an isolated reaction vessel 4 Figure 2: Structure of lutetium disilicate. The circles represent the lutetium ions. constructed in our laboratory. A layered silicate and an aqueous solution are the silicon and RE(III) sources, respectively. Up to now, a set of expandable clay minerals had been studied, rendering conclusions on the relationship between mineralogical compositions and reactivity.3The reaction temperatures were always below the critical one of water; thus both vapour and liquid phases coexist throughout the whole reaction. Reconstructive structural changes occurring in the layered silicate are always analyzed studying the long-range order by X-ray powder diffraction (XRD), the chemical environment of the main constituent elements of the lattice by magic-angle spinning nuclear magnetic resonance spectroscopy (MAS-NMR) and the microstructural and microchemical composition by electron microscopy (SEM) and by energy-dispersive X-ray (EDX). As it has been already mentioned, it has been suggested7that DBs might bring about an apparent violation of Arrhenius law, leading to chemical reactions being observable at much lower temperatures than expected.3It has also been suggested that DBs in the interlayer potassium sheet may be responsible for the dark lines observed in crystals of mica muscovite12,16 and the existence of moving breathers in mica has been recently proven through a sputtering experiment.13 It implies a LTRT process occurring in a nonexpendable layered silicate, which contradicts the mechanism published by some of the present authors to explain the synthesis of RE-disilicates as associated to the capacity for expanding.3In order to confirm experimentally the hypothesis that the LTRT phenomenon is not an exclusive feature of clay minerals, but can also be attributed to non-expandable layered silicates, we 5 have performed the hydrothermal synthesis of Lu-disilicate from muscovite, under the same experimental conditions as those used for clay minerals. When muscovite is hydrothermally treated in stainless steel reactors at 300 ◦C for 72 hours with lutetium nitrate 0.05 M solution, it gives rise to Lu2Si2O7. Figure 3a shows the XRD diagram of the untreated muscovite. It reveals numerous hkl basal reflections compatible with the 2M1polytype and a perfect ordering of the layers. After the hydrothermal treatment, the XRD pattern shown in figure 3b exhibits a number of specific reflections which are consistent with the development of a new crystalline phase Lu2Si2O7 (JCPDS file number 76-1871). The SEM photograph of the mica shows big flakes whereas the sample submitted to hydrothermal treatment, in addition reveals irregular, rough particles, which correspond to Lu2Si2O7. The 29Si MAS NMR spectra of both expandable and non–expandable layered silicates show an evolution from Q3silicon environment to Q1environment.17 This result supports that the LTRT phenomenon is common to expandable and non-expandable layered silicates, both containing an interlayer sheet of cations to balance the layer charge. Note that although the hydrothermal treatment of layered silicates, muscovite and others, is, effectively, at the origin of the development of a new crystalline phase, as a necessary experimental condition, it is not sufficient. In the case of a silicon source different from a layered silicate, which has been used by the authors for the first time, the process does not occur. An example is illustrated in the Fig. 4. This figure shows the 29Si MAS NMR of SiO2submitted to a hydrothermal treatment at 300 ◦C (b) and the typical spectrum of the lutetium disilicate (a). The signal of the SiO2after hydrothermal treatment remains in the chemical shift range typical of Q4 environment (typical of SiO2), which is clearly different for Q1environment (typical for the new phase Lu2Si2O7). Thus, it demonstrates that under hydrothermal condition, SiO2does not transform in a new phase. 4 The conventional chemical kinetics approach Transformation processes of minerals in which there is a major reorganization with bonds and, even, change in the chemical composition are classified as reconstructive. These obey mechanisms which involve very high activation energies when the rupture of strong bonds are involved. Synthesis of RE-disilicates from layered silicates requires the rupture of silicon-oxygen bonds, which are considerably stronger than the bonds between any other element and oxygen. Silicate minerals make up the vast majority of rocks 6 Figure 3: XRD pattern and SEM micrography of untreated (a) and hydrothermally treated (b) muscovite. m=muscovite, *=Lu2Si2O7. The composition of the rough particle inside the circle is compatible with Lu2Si2O7 and their reconstructive transformation processes show activation energies as high as 200 kJ/mol or even higher.18 Therefore, these transformations can be observed in silicate-based minerals at temperatures higher than approximately 1000 ◦C and are apparently impossible at lower ones. It is well known that the over-all effect on temperature on the reaction rate constant kis expressed by Arrhenius law: k=Aexp(−Ea/RT) (1) where Aand Eaare the frequency factor and the activation energy, respectively. According to it, and without any further insight into the mechanism, the rate of a reaction is about 109times faster at 1000 ◦C than it is at 300 ◦C for an activation energy of 200 kJ/mol. It can be shown that the second parameter in the Arrhenius law, namely, the frequency factor A, does not support LTRT synthesis of RE-disilicates as well. Usually, reaction rates are described within the frame of the quasi– equilibrium activated state model.19 According to the latter, the necessary 7 Figure 4: 29Si MAS NMR spectra: a) Lu2Si2O7b) SiO2treated with 50 ml of Lu(NO3)3at 300 ◦C condition for a transformation to take place at a measurable rate is that a sufficient number of atoms have enough energy to achieve the transition state. This energy is supplied by thermal fluctuations. The rate of reaction is then simply the number of activated complexes passing per second over the potential barrier. Applying the conventional transition state (activated complex) theory20–25 to our case, the simplest formulation of the mechanism can be cast in the following form: S (layered silicate) +RE(III)(aq) ↔[S-RE(III)]∗ [S-RE(III)]∗→RE2Si207(2) The rate constant kfor the reaction can be derived by assuming that the transition state (or activated complex) is in equilibrium with the reactants. If C∗represents the concentration of the transition state then the equilibrium constant is: K∗=C∗ [RE(III)] (3) The rate constant kand the equilibrium constant K∗are related by the Eyring equation: k=ν K∗,(4) with ν= kBT/h, h being the Planck constant. To be precise, the rhs of the above expression should be multiplied by a factor η, the transmission coefficient, which is the probability that the complex will dissociate into products instead of back into reactants. For most 8 reactions ηis between 0.5 and 1.0.19 Through a thermodynamic formulation of K∗, it results: k=ν η exp(−∆G0∗/RT) = ν η exp(∆S0∗/R) exp(−∆H0∗/RT) = Aexp(−Ea/RT),(5) where the superindex 0stands for normal conditions. In liquid and solid systems, the p∆V0∗term is negligible and ∆H0= ∆E0∗=Ea. At 300◦C we have calculated the factor Aby substitution of all the parameters in the above expression and it takes the usual value 1013 −1014 s−1for a first order reaction. It does not explain the observation of LTRT synthesis of RE-disilicates as previously concluded from Ea. It is well known that, in parallel to the reorganization of the clay, there must be nucleation of RE2Si2O7crystals and that the activation energy is reduced if the surface and strain energy terms are diminished by good lattice matching across the interface between the new and parent phase.4 However, this is not the present case. It rather seems that the disruption of the tetrahedral sheet is the consequence of localized nonlinear vibration modes, as suggested in Ref. 7. If the vibration modes were delocalized, the relationship shown in Ref. 18 between bond angle, Si-O distance and free energy would be incompatible with an appreciable parent structure-directing character. Ytrium disilicate (Y2Si2O7) has four polymorphs, namely y,β,γand δ. Our LTRT synthesis from the two layered silicates saponite and laponite have shown a structure-directing character of the parent clay by only giving y-Y2Si2O7and δ-Y2Si2O7, the lower and higher temperature polymorphs.4 The relative position of the two tetrahedra in the disilicate unit of the yand γpolymorphs are similar to their position in the tetrahedral sheet of saponite and laponite, which is ∼141◦. The Si-O-Si bond angles in Y2Si2O7polymorphs are 134◦(y), 180◦(β), 170◦(γ) and 158◦(δ). The percentages of yand δphases in the cases of saponite and laponite are explained in Ref. 4 as related with the presence of Al(III) in the precursor framework. Moreover, the importance of maintaining the local Si-O-Si bond angle of precursor structure is shown by the fact that the δ–polymorph has been synthesized at more than 365◦C below the stability range shown in the phase diagram: y-Y2Si2O7→β-Y2Si2O71050(50) ◦C β-Y2Si2O7→γ-Y2Si2O71350(50) ◦C β-Y2Si2O7→δ-Y2Si2O71500(50) ◦C 9 single, exact breathers obtained numerically in the previous section, which are continuation from a single excited oscillator at the anticontinuous limit. Therefore, the numerical curve Pb(E) can be seen as a numerical spectrum for the different breather energies and forms of vibration. The ideal objective, as with other types of spectra, would be to know each type of breather, with its dispersion curve E(νb), and the relative probability of its appearance and to be able to reproduce exactly the numerical spectrum. In principle, all the breather types could be obtained exactly with different conditions at the anticontinuous limit with different frequencies and by path continuation by changing the frequency and studying the different branches at the possible bifurcations. This is a long and difficult task that we do not pursue here. Instead, we try to fit approximately the numerical Pb(E), with a small number of breather types, each one characterized by its minimum energy ∆, maximum energy EM, which can be ∞, and parameter z(see B.1 and B.2), each breather type with a different probability to occur. In this way, we know that we cannot fit exactly the numerical spectrum because we are most probably substituting a number of breather types with an average one. In any case, our numerical Pb(E) is also an approximation, as we would need a extremely large number of simulations to obtain the actual curve. Note that the breather spectrum will not appear in an experimental one for three reasons: a) The number of breathers is about 10−3the number of phonons, and the spectrum is basically dominated by the one-phonon transitions; b) Breathers are localized and, therefore, they cannot be excited by infrared, raman of neutron spectra.36 Figures 9 and 10 show the numerical and theoretical, probability densities and cumulative probabilities, respectively. The parameters of the breathers are: ∆ (kJ/mol) 23.9 36.6 41.4 62.2 67.3 82.9 z1.50 1.17 3.00 0.52 2.07 1.80 EM(kJ/mol)) – 46.9 – – – 94.4 Probability 0.103 0.026 0.281 0.097 0.202 0.290 6.3 Effect on the reaction rate According to Ref. 18 the lowest estimates of the activation energy for a reconstructive transformation as the one described above are Ea= 100 − 200 kJ/mol. Let nph(Ea)≃exp(−βEa) and nb(Ea) be the mean number of phonons and breathers, respectively, per site, with energies E≥Ea. Then, nb(Ea) = hnbiCb(Ea), with hnbi ∼ 0.92 ·10−3, the mean number of 16 0 20 40 60 80 100 120 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 E (kJ/mol) P(E) (kJ/mol)−1 Figure 9: Breathers spectra, i.e., breather probability densities Pb(E) obtained numerically (line with dots) and theoretically (continuous line). The latter is obtained by considering six different types of breathers. See text. breathers per site obtained numerically, and Cb(Ea), the cumulative probability described in the previous subsection, obtained with six different types of breathers in order to fit the numerical probability density. The ratio of the number of breathers to phonons becomes nb(Ea)/nph(Ea)∼104−105. The reaction rate constant with breathers would be kb=Abnb(Ea), with Ab, the frequency factor for breathers Ab, which should be different from A, the frequency factor for phonons. We will assume that Ab=Afor the purpose of comparison. The ratio of reaction rates for breathers and phonons would be kb/k =nb(Ea)/nph(Ea)∼104−105for Ea= 100 −200 kJ/mol. In other words, as the three days experimental time leads to about 30% of the transformation performed, the time without breathers to obtain the same result, would be 104−105times larger and, thus, completely unobservable. Two factors are likely to increase further the reaction rate with breathers. First, since a discrete breather is strongly localized, it seems much more capable of delivering the energy for breaking a Si-O bond, which implies that Abshould in effect be much larger than A. Second, for larger systems than the one used in our simulations, the fluctuations may have larger energies and, therefore, excite other types of breathers with higher energies, as, for 17 0 20 40 60 80 100 120 0 0.2 0.4 0.6 0.8 1 1.2 1.4 E (kJ/mol) P(E) (kJ/mol)−1 Figure 10: Breather cumulative probabilities Cb(E) obtained numerically (line with dots) and theoretically (continuous line). The latter considering six different types of breathers. See text. example, breathers with frequencies above the phonon band, which have energies between 240 and 500 kJ/mol. These breathers would increase the fraction of breathers above the activation energy and therefore the reaction rate. 7 Summary and conclusions Low temperature reconstructive transformations (LTRT) have been achieved in layered silicated by some of the authors at temperatures about 600◦C lower than previously reported. This is a phenomenon for which there is presently no plausible explanation since the bonds involved are the same as in other transformations. New experiments performed by some of the authors on mica muscovite, a non-expandable silicate have discarded their previous hypothesis of LTRT been caused by the expansion of the intersheet layer. We have constructed a model for breathers in the cation later, for which we have obtained reasonable parameters, and with a mixture of numerics and theory we have estimated their effect in the reaction rate. The results are 18 that they would increase enormously the reaction rate and, thus, explain the observed LTRT. This can be easily explained by the fact that, although there are much less breathers than phonons, there are many more with energies above the expected activation energy. Certainly, the statistical theory of breathers is only an approximation, and the numerics cannot be precise at larger energies for which there are so few breathers, except if an enormous number of simulations could be performed. However, an established fact in breather theory is that large breathers have longer life time than small ones and thus tend to overpopulate the regions of high energies if compared to Maxwell-Boltzmann statistics for phonons. Moreover, as they are localized, it seems that they can deliver more easily the required energy to break bonds. The sum of this facts, i.e., localization, much higher number of breathers above a given activation energy and apparent diminution of the activation energy suggests that DBs are good candidates to explain LTRTs. Acknowledgments JFRA and JC acknowledge sponsorship by the Ministerio de Educaci´on y Ciencia, Spain, project FIS2004-01183. MDA, MN and JMT acknowledge sponsorship from the same body, projects MAT2002-03504 and CTQ200405113. JFRA acknowledges the hospitality and the spectra performed at CNRS-LADIR. All the authors acknowledge Prof. R. Livi, from Florence University for useful discussions. Appendices A Phonon statistics A.1 Quantum statistics of one oscillator A quantum harmonic oscillator with frequency ω, equal to its classical one, has energies En= (n+1 2)~ω. If in contact with a thermal bath at temperature T, the probability that it has energy Enis given by P(En) = Aexp(−β En) = Aexp(−β(n+1 2)~ω), with β= 1/kBT, kBbeing the Boltzmann constant. A can be obtained by the normalization condition P∞ n=0 Pn=A Z = 1, where Z=P∞ n=0 exp[−β(n+1 2)~ω] is known as the partition function for the oscillator. Zis a geometric series which can be easily summed leading 19 to: Z=exp(−β~ω 2) 1−exp(−β~ω).(A-1) Therefore A= 1/Z and Pn≡P(En) = exp[−β(n+1 2)~ω] Z= exp(−βn~ω)[1 −exp(−β~ω)] (A-2) Note that nis the excitation number of the oscillator, however, in solid state physics, where ωis the frequency of a normal mode of the solid, it is customary to speak of nas the number of phonons with frequency ω. Once known Pna number of quantities can be readily calculated. The mean energy hEi=Z−1P∞ n=0(n+1 2)~ωexp[−β(n+1 2)~ω] = −Z−1∂Z/∂β = −∂log(Z)/∂β, which leads to: hEi=1 2+1 exp(β~ω)−1~ω . (A-3) Consequently, the mean excitation number or mean number of phonons is: hni=1 exp(β~ω)−1.(A-4) For high temperatures ~ω/kBT << 1, the mean energy becomes hEi ≃ kBTwhich is the classical mean energy of a harmonic oscillator and hni ≃ kBT/~ω. Of particular importance for the present problem is the cumulative probability C(Ea), i.e., the probability that the oscillator has energy Eaor higher above the ground state ~ω/2, i.e., that it can deliver the energy Ea. Let mbe the minimum integer so as m~ω≥Ea, i.e., m=⌈Ea/~ω⌉, with ⌈x⌉, the ceiling function that rounds xtowards plus infinity. Then C(Ea) = Z−1P∞ n=mexp[−β(1 2+n)~ω] = Z−1P∞ n=0 exp[−β(m+1 2+n)~ω] =Z−1exp(−βm~ω)×P∞ n=0 exp[−β(1 2+n)~ω], therefore: C(Ea) = exp(−βm~ω) = exp(−β⌈Ea ~ω⌉~ω).(A-5) If mis of the order of a few tens, the expression above approaches to the classical expression Cclass(Ea) = exp(−βEa). Note, however, that the quantum probability is somewhat smaller. The ratio between the quantum probability and the classical one is between exp(−β~ω) and 1. For a frequency as the one given here for the phonons in the K+plane ω0= 3.16 ·1013s−1and T= 573K this ratio is exp(−β~ω0)≈0.65 and 0.81 for T= 1173K. 20 A.2 Normal modes and phonons Let us consider a solid in the linear approximation, with Nfdegrees of freedom, with Nf= 3 ×Nafor a three-dimensional solid with Naatoms. There are Nfnormal modes with frequencies ωiand wave numbers ki, each one equivalent to a independent harmonic oscillator, and therefore, the previous subsection can be applied to it. The properties of the solid are simply the sum of the properties of the isolated oscillators, just adding the subindex i to the formulae in the preceding section and summing up, i.e., the excitation number hniiof the mode i(or the number of phonons) and the energy of the solid are given by: hnii=1 exp(β~ωi) + 1 ;E= Nf X i=1 (hnii+1 2)~ωi(A-6) The number of modes Nph(Ea) with energy larger or equal to Eaabove their ground state is Nph(E≥Ea) = Nf X i=1 exp(−β⌈Ea ~ωi⌉~ωi).(A-7) If Eais about a few tens larger than any ~ωi, then we simply have: Nph(E≥Ea)≃Nfexp(−βEa),(A-8) which is the classical expression. Again, the ratio between the quantum and classical expressions of Nph(E≥Ea) is somewhat smaller than the unity. Therefore, the cumulative probability C(Ea) i.e., the fraction of modes with energies greater of equal to Ea, becomes: C(Ea)≃exp(−βEa).(A-9) B Breather statistics B.1 Breathers with hard on-site potential The breather statistics theory developed in Ref. 34 for 2D breather in a system with hard on–site potential is based is some simple hypotheses, which, in principle, can be fairly general: 1. An established fact is that breathers in two and three dimensions have a minimum energy ∆.33 21 2. The rate of creation of breathers with energy E,B(E), is proportional to exp(−βE), since breathers form from fluctuations through an activation process. 3. The probability per unit time that a breather with energy Eis destroyed, D(E), is postulated to be inversely proportional to (E−∆)z, with za constant (which means as other constants hereafter that it does not change with the energy E) that depends on the system. This law is the simplest mathematical expression that takes into account that large breathers have longer lives than smaller ones, with the lifetime of breathers with minimum energy ∆ equal to zero. It has to be considered as an approximation as it is not derived from first principles. Let Pb(E) dEbe the probability of existence (or the mean fraction) of breathers with energy between Eand E+ dE. The rate of destruction of breathers with energy Eis proportional to D(E) and Pb(E), therefore, exp(−βE) = A Pb(E)(E−∆)−z,Abeing a constant, or, A Pb(E)=(E− ∆)zexp(−βE). Since R∞ ∆Pb(E)dE= 1, Acan be obtained using the change of variable y=β(E−∆): A=Z∞ ∆ (E−∆)zexp(−βE) = exp(−β∆) βz+1 Z∞ 0 yzexp(−y)dy= exp(−β∆)β−(z+1)Γ(z+ 1) ,(B-1) with Γ(z+ 1) = R∞ 0yzexp(−y)dy, the Gamma function. Thus, the probability of breathers with energy between Eand E+ dEis given by: Pb(E) = 1 A(E−∆)zexp(−βE) = (E−∆)zexp(−βE) exp(−β∆)β−(z+1)Γ(z+ 1) = βz+1 Γ(z+ 1)(E−∆)zexp[−β(E−∆)] .(B-2) The mean energy is given by: hEi=Z∞ ∆ E Pb(E)dE= ∆ + Z∞ ∆ (E−∆)Pb(E)dE= ∆ + Z∞ ∆ βz+1 Γ(z+ 1)(E−∆)z+1 exp[−β(E−∆)]dE= ∆ + 1 βΓ(z+ 1) Z∞ 0 yz+1 exp(−y)dE= ∆ + Γ(z+ 2) βΓ(z+ 1) = ∆ + (z+ 1) kBT . (B-3) 22 The cumulative probability Cb(E), i.e., the probability that a breather has energy higher than E, is given by: Cb(E) = Z∞ E Pb(E)dE=Z∞ E βz+1 Γ(z+ 1)(E−∆)zexp[−β(E−∆)] dE= 1 Γ(z+ 1) Z∞ β(E−∆) yzexp(−y)dy=Γ(z+ 1, β(E−∆)) Γ(z+ 1) ,(B-4) where Γ(z+ 1, x) = R∞ xyzexp(−y)dyis the first incomplete Gamma function.35 The energy for which the probability is maximum is given by E(Pmax) = ∆+zkBT, which shows that breathers tend to populate higher energies than phonons. As an example, Fig. 11 shows Pb(E) and Cb(E) for breathers with ∆ = 20 kJ/mol and z= 2 and the equivalent magnitudes for phonons. The large energies the larger values of Cb(E) soon compensate for the much smaller number of breathers than phonons (around 10−3). Although the extrapolation of Cb(E) to large energies has to be done with caution and a more elaborate theory has yet to be developed, the basic fact that breathers tend to populate higher energies than phonons can be accepted. Piazza et al34 have succeeded in fitting the cumulative probability in Eq. B-4 with the observed one in numerical experiments, which proves that the hypotheses above are reasonable. B.2 Breathers with maximum energy Hereafter we modify slightly their theory developed above. The on-site potential for the 2D system in Ref. 34 is hard, i.e., the frequency of the isolated oscillators increases with the frequency, with the consequences that the breather frequency lies above the phonon band, its energy increases with its frequency and can be considered as unbounded. For soft on-site potentials or, as in the present paper, potentials with both soft and hard parts, breather energies may have an upper limit because the breather frequency enters the phonon band or because a bifurcation where the breather disappears or transforms into a different one, as a multibreather or a breather with different symmetries. The changes are obtained by introducing an upper limit EMin the integrals with respect to the energy. Then, Eq. (B-1) becomes (with y=β(E−∆)): A=ZEM ∆ (E−∆)zexp(−βE) = exp(β∆) βz+1 Zβ(EM−∆) 0 yzexp(−y)dy= exp(β∆)β−(z+1)γ(z+ 1, β (EM−∆)) ,(B-5) 23 0 20 40 60 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 E (kJ/mol) P(E) (kJ/mol)−1 phonons breathers 0 50 100 10−10 10−8 10−6 10−4 10−2 100 E (kJ/mol) C(E) phonons breathers Figure 11: Comparison of the phonon and breather probabilities densities (left) and cumulative probabilities (right). The breather values have been obtained for z= 2 and ∆ = 20 kJ/mol. The temperature is T= 600 K and kBT≈5 kJ/mol. where γ(z+ 1, x) = Rx 0yzexp(−y) dyis the second incomplete gamma function.35 Therefore, the probability density becomes: Pb(E) = 1 A(E−∆)zexp(−βE) = βz+1(E−∆)zexp[−β(E−∆)] γ(z+ 1, β (EM−∆)) .(B-6) The cumulative probability becomes: Cb(E) = ZEM E Pb(E) dE=ZEM E βz+1(E−∆)zexp[−β(E−∆)] γ(z+ 1, β (EM−∆)) dE= 1 γ(z+ 1, β(EM−∆)) Zβ(EM−∆) β(E−∆) yzexp(−y)dy= 1 −γ(z+ 1, β(E−∆)) γ(z+ 1, β(EM−∆)) .(B-7) For EM>> ∆, the expressions above for Pb(E) and Cb(E) transform into the expressions in Eqs. (B-3,B-4) 24 In a system like ours there are different types of breathers and the probabilities or cumulative probabilities calculated above correspond to each type with different minimum and maximum energies ∆ and EM(or without maximum energy), and parameter z. The total number of breathers and the relative probability of each type are unsolved questions. The latter probably depends on the temperature, the breather energies, the phase space occupied by each breather and its equivalent ones through symmetries, and the breather profile, which might be excited more or less easily by the phonons. 25