{'fraction_non_alphanumeric': 0.029326481923288676, 'fraction_numerical': 0.01714075977243531, 'mean_word_length': 4.452371422853712, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Structural stability and lattice dynamics of SiO2 cristobalite

Among the phases of SiO2 are alpha and beta cristobalites, which have a long and somewhat controversial history of proposed structural assignments and phase-transition mechanisms. Recently, Zhang and Scott found new indications that the higher-temperature beta phase has space group I-42d and, by assuming a group-subgroup relationship between phases, they argued that the lower-temperature alpha phase should have lower symmetry than that of the widely-accepted P4_12_12 space group. With this motivation, we use first-principles calculations to investigate the energy, structure, and local stability of P4_12_12 and I-42d structures. We also compute the frequencies of the zone-center phonon modes in both structures, as well as certain zone-boundary modes in the I-42d structure, and compare with experiment. We then argue that the various P4_12_12 and I-42d enantiomorphs can be grouped into three clusters, each of which is identified with a three-dimensional manifold of structures of P2_12_12_1 symmetry in which the P4_12_12 and I-42d appear as higher-symmetry special cases. We find that there are relatively high energy barriers between manifolds, but low barriers within a manifold. Exploring the energy landscape within one of these manifolds, we find a minimal-energy path connecting P4_12_12 and I-42d structures with a surprisingly low barrier of ~5 meV per formula unit. Possible implications for the phase-transition mechanism are discussed.


I. INTRODUCTION
The fact that SiO 2 can exist in numerous crystalline and amorphous forms, and its status as one of the most prevalent compounds on earth, has stimulated a long history of experimental and theoretical investigation. Here we focus on the α ("low") and β ("high") cristobalite phases, which are stable near the melting temperature and are metastable at room temperature.
The structure of the higher-temperature β phase has a history of controversy. Early indications of a cubic structure with 180 • bond angles (space group F d3m) 1 were challenged by others 2,3 who hypothesized that the true β-phase structure has lower symmetry and that the apparent cubic structure arises from averaging over spatial domains or dynamical fluctuations. In particular, Wright and Leadbetter 3 argued for a tetragonal structure belonging to space group I42d (D 12 2d ). While some subsequent works have provided support for this identification, 4,5,6 other authors maintain that it is better to describe the β phase as a dynamically disordered one having overall F d3m symmetry but with a large population of rigid-unit-mode (RUM) fluctuations. 7,8 To some degree, the argument may be semantic; if the fluctuations have strong short-range correlations in space and time, it is difficult to distinguish this picture from one of dynamic domains of a lower-symmetry structure. 9 Thus, for example, either picture may be able to explain the fact that there are two more first-order lines in the Raman and infrared spectra than would be expected from F d3m symmetry, 6 and the question of which description is "correct" might depend on the time and length scales of the experimental probes in question.
In contrast, the assignment of the α-cristobalite phase to the tetragonal P 4 1 2 1 2 (D 4 4 ) space group 10 has un-til recently been noncontroversial. However, based on a reexamination of Raman and infrared vibrational spectroscopies, Zhang and Scott 6 have recently raised new questions about the identity of the α phase. By using Raman spectroscopy to study small single crystals of βcristobalite, these authors argued that the β structure must be D 2d , not cubic, and assuming a group-subgroup relationship for the β-to-α transition, concluded that the α phase should have some lower symmetry such as D 2 instead of D 4 . The apparent D 4 symmetry of α-cristobalite could result from spatial or dynamic averaging over D 2 domains, in analogy to what had been proposed for the β phase. To support their assumption that a groupsubgroup relationship should hold, Zhang and Scott also pointed to the temperature dependence of the optical phonon frequencies near the phase transition as being inconsistent with a reconstructive phase transition 11 and as suggesting a nearly second-order behavior, although arguing in the opposite direction are the facts that the latent heat and volume change at the transition are quite substantial. 12 In their paper, Zhang and Scott 6 reexamined earlier Raman and infrared spectroscopic measurements not only on the α-and β-cristobalite SiO 2 , 13 but also on α and β AlPO 4 (Ref. 14) and α BPO 4 (Ref. 15) cristobalites. Note that the replacement of Si atoms by Al and P (or B and P) atoms immediately reduces the symmetry according to P 4 1 2 1 2 (D 4 4 ) → C222 1 (D 5 2 ) for the α phase 16 and I42d (D 12 2d ) → I4 (S 2 4 ) for the β phase. Also of possible relevance is the pressure-induced phase transition from α-SiO 2 to a high-pressure monoclinic cristobalite phase. 17 The relationship of these other cristobalites to the α and β phases of SiO 2 is an interesting avenue for future exploration, but falls outside the scope of the present work.
In the present work, we have carried out first-principles calculations for candidate α and β cristobalite structures in the framework of density-functional theory (DFT) in order to check the stability of both phases and to explore the energy landscape connecting them. We have also calculated phonon frequencies and infrared activities for both α and β phases, and explored how the phonon modes in the different phases are related to each other and to those of the high-symmetry cubic phase. Our calculations are effectively zero-temperature ones, and thus cannot properly treat the temperature-induced α-β cristobalite phase transition. Nevertheless we hope that the information obtained from these calculations can eventually be built into a realistic statistical-mechanical theory that correctly takes the RUM fluctuations into account in its description of the α and β phases at experimentally relevant temperatures.
The paper is organized as follows. In Sec. II we give a brief review of α and β cristobalite structures and describe the methods used in the calculations. Then, in Sec. III, we present the results of our calculations of structural and lattice vibrational properties of the two phases and of the energy landscape connecting them. We discuss those results in Sec. IV. Finally, we summarize the work in Section V.

A. Cristobalite structures
In order to describe the structures of the SiO 2 α and β cristobalite phases, it is easiest to start by considering the "ideal cristobalite" structure, which is constructed by placing Si atoms in a diamond structure with oxygen atoms located midway between each pair of nearestneighbor Si atoms. This structure has the space group F d3m (O 7 h ) and has two formula units per primitive unit cell. Each Si atom with its four surrounding O atoms forms a tetrahedron, so the whole structure can be visualized as a network of tetrahedra connected at their apices.
The generally accepted structure of α-cristobalite is arrived at by starting from the ideal structure and making nearly rigid rotations of the tetrahedra about [100] and [010] axes (in the original diamond cubic frame), leading to a tetragonal structure with its axis alongẑ. This is illustrated in Fig. 1(a), but in the conventional tetragonal frame, related to the original cubic frame by a 45 • rotation aboutẑ. The tetrahedral rotations are also accompanied by small strains and tetrahedral translations needed to keep the apices coincident, as would be FIG. 1: (Color online) Projection on x-y plane of the (a) α and (b)β structures, proposed as candidates for α and β cristobalite phases respectively. Darker shading is used to represent more distant tetrahedra so that the spiral structure of the connected tetrahedrons becomes evident; double vertical lines indicate that the adjoining tetrahedra are actually disconnected because they are separated in the z-direction. expected from enforcement of the rigid-unit constraints. The space group of the structure is P 4 1 2 1 2 (D 4 4 ), and since the four rotations shown in Fig. 1(a) are all different, the number of formula units per primitive unit cell is now increased to four.
As mentioned earlier, diffraction experiments on the β-cristobalite phase tend to give inconclusive results because of spatial and dynamical averaging. Nevertheless, based on the comparison of structure factors predicted by various disorder models and the ones obtained in their xray diffraction experiments, Wright and Leadbetter 3 concluded that the β-cristobalite has local I42d (D 12 2d ) spacegroup symmetry. Their proposed structure can also be constructed from the ideal structure, but this time by rotating all the tetrahedra around theẑ axis, yielding the structure shown in 1(b). The number of formula units per primitive unit cell remains at two as in the ideal structure (although it can alternatively be described, as in Fig. 1(b), by a doubled conventional cell containing four formula units). Again, the structure is highly consistent with the rigid-unit constraints.
Because we do not want to presuppose an identification of a particular experimentally observed phase with a particular crystal structure, we henceforth adopt a notation in which the phases are identified by labels "α" and "β" without tildes, whereas the putative crystal structures shown in 1(a) and (b) will be referred to as "α" and "β" structures, respectively. Our working hypothesis is that the α and β phases have microscopic crystal structures of typeα andβ respectively, but we adhere to a distinction in the notation in order to discriminate clearly between the specified structures used in our calculations and the hypothetical identification of these with experimental phases.

B. Computational methods
The calculations were carried out using the ABINIT implementation 32 of density-functional theory with Perdew-Burke-Ernzerhof 33 version of the generalized gradient approximation (GGA) for electron exchange and correlation. Since it is the smallest unit cell that contains bothα andβ structures, all calculations were performed on the four-formula-unit computational cell shown in Fig. 1, even though the primitive cell is smaller in theβ structure. The Brillouin zone was sampled by a 4 × 4 × 4 Monkhorst-Pack grid. 34 Structural properties were computed using projector augmented-wave 35 potentials converted from ultrasoft pseudopotentials 36 with a plane-wave cutoff of 22 Ha unless otherwise specified, while phonon frequencies, eigenvectors, and Born charges were computed 37 using norm-conserving Trouiller-Martins pseudopotentials 38 at an energy cutoff of 50 Ha (after repeating the structural relaxation using these potentials). The acoustic sum rule was imposed on the force constants, and charge neutrality was imposed on the Born charges. Throughout the paper, the symmetry analysis associated with crystal space groups has been carried out using the Bilbao package. 39,40 III. RESULTS

A. Structural properties ofα andβ structures
We started our calculations by considering the ideal cubic structure and relaxing its volume, obtaining a c =7.444Å for the lattice constant of its eight-formulaunit cubic cell. Then, working in the frame of the fourformula-unit tetragonal cell, we found the phonon frequencies at the Γ point of its Brillouin zone, corresponding to phonons at the Γ point and one X point [namely (2π/a c )(001) or equivalently (2π/a c )(110) in the cubic frame] of the primitive two-formula-unit fcc cell. For the "ideal structure" of space group F d3m, the symmetry decomposition of these phonons into irreducible representations is (The translational T 1u mode has been omitted.) The E u mode and all X modes are doubly degenerate; the T 1u , T 2u and T 2g modes are triply degenerate; and A 2u is non-degenerate. We found that the triply-degenerate T 2u mode at Γ is unstable with an imaginary frequency of i83 cm −1 . All other optical phonons have real frequencies, the lowest being at 250 cm −1 . Furthermore, one of the doublydegenerate (X 4 ) modes is unstable with a frequency of i53 cm −1 . We thus conclude that the ideal cristobalite structure is unstable with respect to these distortions. 35  Next we imposed distortions corresponding to these unstable modes and did a full relaxation of the structure subject to the symmetry constraints of the resulting space group. The unstable (i53 cm −1 ) mode at X leads to the space group P 4 1 2 1 2 (or P 4 3 2 1 2) which corresponds toα-cristobalite, while the (i83 cm −1 ) mode at Γ takes us to the space group I42d ofβ-cristobalite. The energy of the relaxed ground state as a function of volume per formula unit is shown for both cases in Fig. 2, with the energy of the cubic phase also shown for reference. The corresponding structural parameters at the energy minimum are given in Table I.
From Fig. 2 it is clear that theα andβ structures indeed have lower energies than the ideal cristobalite when the volume becomes smaller than some critical volume V 0 ∼ 55Å 3 . (Above this volume, the imposed distortions disappear during relaxation and the structure returns to the ideal one.) We find that both theα andβ structures have a quite similar dependence of energy on volume. According to our calculation, the relaxedβ structure has a slightly lower energy than that of theα structure (12 meV per formula unit). This appears to be in conflict with the experimental situation, since the α phase is experimentally more stable at lower temperatures. However, when we repeated our calculations using a local-density approximation (LDA) exchange-correlation functional, 41 theβ structure was found to be lower by 1 meV per formula unit. We thus conclude that the small energy difference between the two nearly-degenerate structures is a quantity that is too delicate to be reliably obtained by our DFT calculations. A similar discrepancy between the results from LDA and GGA functionals was found in Ref. [24].
We also analyzed what happens to the bond lengths and angles in theα andβ structures as a function of volume. The details of this analysis are deferred to App. A. Briefly, for V < V 0 , the O-Si-O bond angles and Si-O bond lengths inside the tetrahedra are found to remain almost constant, while the Si-O-Si bond angles change by ∼35 • . The details are shown to be very close to the predictions of a picture of tilting of perfectly rigid tetrahedra. The fact that the three phases are indistinguishable for V > V 0 is also easily explained, since the tilts of rigid tetrahedra can only decrease the volume of the ideal structure. Thus, for V > V 0 the tetrahedra cannot remain rigid and the Si-O bond length must increase, and only when V becomes smaller than some volume V 0 will one of the RUMs condense in the structure in order to maintain the preferred Si-O bond lengths.

Phonons at Γ inα cristobalite
We next repeated the calculation of the phonon frequencies for the fully relaxedα andβ cristobalite structures. We did this in order to compare with experimental measurements and also to check the stability of the structures and to investigate, at least at harmonic order, the nature of the energy landscape around these structures. An analysis extending beyond the harmonic approximation will be presented in Sec. III C 2.
The decomposition of the optical Γ phonons into irreducible representations for theα structure in space group (The translational A 2 and E zero modes have been omitted.) Only the E modes are doubly degenerate; all others are non-degenerate. Tables II and III present the phonon frequencies at the Γ point for the fully relaxedα-cristobalite structure. All phonon frequencies are positive, although some appear to be rather low in frequency. For the infrared (IR) active modes shown in Table II, the transverse mode frequencies were computed initially, and their mode dynamical charges were also computed using where ξ i,λβ is an eigenvector of the dynamical matrix, Z * i,αβ is the Born atomic charge tensor, and M i is the mass of the i-th atom in amu. The norms of the modecharge vectorsZ * λ = [ α (Z * λ,α ) 2 ] 1/2 are also given in the Table. The longitudinal dynamical matrix was then constructed and diagonalized using standard methods, 37 and the resulting LO mode frequencies are presented in the last column of Table II. It can be seen that there are wide variations in the mode dynamical charges, and consequently, large variations in the LO-TO splittings. Irrep

Phonons at Γ and M inβ-cristobalite
Similar calculations of phonon frequencies were also carried out for theβ structure proposed by Wright and Leadbetter 3 for β-cristobalite. Since the primitive cell of theβ structure contains only two formula units while theα structure contains four, it should be kept in mind that the Γ point of theα structure maps not only into the Γ point of theβ structure, but also into a second point that would be denoted as X = (2π/a c )(110) in the original fcc frame, or (2π/a)(100) (where a ≃ a c / √ 2) in the rotated frame of Fig. 1; we shall refer to this as the M point in accordance with the conventional labeling of the body-centered-tetragonal (bct) primitive cell in the latter frame. The decompositions of the Γ and M phonons into irreducible representations for theβ structure in space group I42d are (Translational B 2 and E zero modes have been omitted.) All M-point modes and Γ-point E modes are doubly de- generate, while other modes are non-degenerate.
In Tables IV and V we present the our results for the Γ-point and M-point phonon modes, respectively, of β-cristobalite. The frequencies given for the IR-active modes at Γ are the transverse ones only. The tables also show the correspondences between the Γ modes in theα-cristobalite structure and the Γ and M modes in theβ-cristobalite structure, as determined by comparing phonon eigenvectors.

Relation to unstable phonons in the cubic phase
The triply degenerate Γ-point mode of the cubic structure having imaginary frequency i83 cm −1 , which condensed to form theβ structure, now has positive frequencies of 289 cm −1 for the non-degenerate A 1 mode and 126 cm −1 for the E doublet in theβ structure. This same triplet corresponds, in theα structure, to the lowestfrequency phonon of frequency 29 cm −1 , which has symmetry B 1 , and to an E doublet at 127 cm −1 having almost the same frequency as in theβ structure. The doubly-degenerate unstable mode of the cubic structure at i53 cm −1 , which condensed to form theα structure, now appears in theα structure at frequencies 197 cm −1 and 103 cm −1 with symmetries A 1 and B 1 , respectively. In theβ structure, on the other hand, the same doublet appears as the lowest-frequency phonon in that structure, namely the doublet at 35 cm −1 with symmetry M 3 M 4 . These relations between the unstable phonons in the "ideal structure" and the phonons inα andβ structures are shown in Table VI.

Comparison with experiment for α-cristobalite
In view of the recent questions that have been posed about the identity of the α-cristobalite phase, 6 we have carried out a more detailed analysis of the phonons in theα structure. In particular, we have calculated the LO frequencies of the Γ-point phonons inα-cristobalite as a function of the angle at which the limitq → 0 is taken. It turns out that the labels A 2 and E are not well-defined at arbitraryq because of mixing between modes of these symmetries. Moreover, it can happen that if one starts VII: Left: Computed mode frequencies and irreps for α-cristobalite, with direction of dynamical polarization in parentheses for IR-active modes. The modes that are adiabatically connected asq is rotated fromx toẑ appear on the same line. All modes other than A2 modes are Raman-active. Right: Tentative assignments to measured mode frequencies in powder samples.

Present theory
Experimental datâ a Not pure LO at allq. b LO-TO splitting is negligible. c Part of structured peak. d Inactive in β phase. e Inβ structure the A2 component also becomes Raman active. f Corresponds to RUM mode in ideal cristobalite. g Observation from [13]. h Observation of 50 cm −1 mode is from [42], all others from [43]. i Empirical assignments from [13].
with an E mode atq ẑ and follows the branch asq is rotated, one arrives at an A 2 mode whenq lies in the xy plane, or vice versa. Experiments have typically been done on powder samples, so that one should in principle average the phonon spectrum over all possible directions for q → 0. Moreover, some phonon modes with E symmetry have a very small LO-TO splitting, so they would most likely appear in experiment as a single line.
For all these reasons, a direct comparison of experimental data with our results as presented in Tables II and  III is problematic. Nevertheless, we attempt such a com-parison in Table VII. Despite the difficulties, the agreement with experimental data is rather good, with a few exceptions that will be discussed shortly. We generally underestimate the experimental frequencies by ∼20 cm −1 for lower frequency phonons and by ∼35 cm −1 for higher frequencies, but otherwise our results reproduce the experimental pattern of frequencies, and the irrep assignments are also consistent with those obtained from empirical models. 13 Furthermore, the identification of the modes that are not expected to be active in the β phase (fifth and ninth columns of Table VII) because they cor-respond to zone boundary modes in the β phase (see Table V) or because they are inactive Γ-point modes (see Table IV) agrees well with the results reported in Ref. [13].
The first anomaly is related to the experimentally observed IR mode with a frequency of 798 cm −1 in the α phase that does not disappear upon transition to the β phase. Finnie et al. [44] explained this by suggesting that a two-phonon process in the β phase replaces the fundamental mode of the α phase. Our calculations identify two almost-degenerate IR modes that are close to this frequency, an A 2 mode at 751 cm −1 and an E mode 752 cm −1 . In theα structure both of these modes are IR active, but the Born charge of the E mode is only 0.02 while that of A 2 is 0.52, which means that the E mode in theα structure is almost invisible. In theβ structure, the A 2 mode disappears since it is no longer at Γ. On the other hand, the E mode remains at the Γ point and its Born charge is increased to 0.51. These results suggest a possible explanation for the "anomaly," namely that there are two IR modes in the α phase; one of them is much weaker than the other, but upon the transition to the β phase, the stronger one disappears by symmetry while the weaker one greatly increases its IR activity. The reason why the 752 cm −1 E phonon in theα structure acquires a larger Born charge upon converting to theβ structure is that it gets some admixture of the 440 cm −1 E mode, which has a much larger Born charge (0.74) than that of the 752 cm −1 mode (0.02).
The second anomaly is related to the 1076 cm −1 mode that is still present upon the transition to the β phase in the form of a fairly broad feature (see Fig. 1 in Ref. [43]), whereas it would be expected to vanish by symmetry. Swainson et al. [13] attributed this mode to a possible higher-order process. We think that it could also be related to the fact that the 1050 cm −1 A 2 mode that is Raman inactive in theα-phase becomes Raman active iñ β-phase.
We also predict several phonon modes that are not seen experimentally, such as the Raman-active modes at 127, 358, 378, 584, 752 and 1048 cm −1 . Since we have not computed Raman matrix elements, it is possible that the Raman intensities are small for these modes. We also find one weak IR-active mode at 752 cm −1 that is not seen in the experiments.
A very low-frequency phonon at 50 cm −1 has been reported in α cristobalite. 42,45,46 We believe this most likely corresponds to the B 1 phonon that we have calculated to appear at 29 cm −1 , corresponding closely to the RUM mode that takes the ideal cubic cristobalite structure into theβ structure. The same conclusion regarding the lowest-frequency B 1 phonon was reached in Ref. [13]. The minimal-energy path between α and β phases that is related to this low-frequency phonon is discussed in Sec. III C 2.
C. SiO2 cristobalite stability analysis As shown in Sec. III B, all calculated optical phonons inα andβ cristobalites have ω 2 > 0, indicating that the relaxed structure is stable with respect to those modes. In view of the suggestion in Ref. [6] that the α phase might locally have D 2 rather than D 4 point-group symmetry, we checked carefully for instabilities leading from theα structure to D 2 structures, but found none. The possible subgroups of P 4 1 2 1 2 (D 4 4 ) having D 2 symmetry (without reduced translational symmetry) are C222 1 (D 5 2 ) and P 2 1 2 1 2 1 (D 4 2 ), and the phonon distortions leading to these symmetry-lowered structures are the ones of B 2 and B 1 symmetry respectively. The lowest-frequency mode of B 2 symmetry is at 275 cm −1 , so there is certainly no sign of an instability there. On the other hand, the lowest-frequency B 1 phonon is nearly soft at 29 cm −1 , suggesting that theα structure is nearly unstable to a spontaneous transformation into the P 2 1 2 1 2 1 structure.
To check this possibility more carefully, we started from the relaxedα structure and followed the distortion corresponding to the 29 cm −1 B 1 phonon, and confirmed that the energy increases monotonically (no double-well structure). Moreover, starting from one of these structures having a small amount of the 29 cm −1 mode frozen in, a subsequent relaxation inside the resulting space group P 2 1 2 1 2 1 lead to a recovery of the starting space group P 4 1 2 1 2. We thus conclude, at least within our zero-temperature first-principles calculations, that theα structure is locally stable, i.e., does not spontaneously lower its symmetry to D 2 .
Nevertheless, the presence of several modes of quite low frequency in theα-cristobalite structure may be suggestive of low-energy pathways leading from theα to thẽ β structure or between domains of theα structure. For example, we have shown above that the lowest-frequency 29 cm −1 mode in theα-cristobalite structure corresponds to a phonon of the "ideal structure" that leads to theβ structure, and vice versa. This might suggest that there is a relatively low energy barrier in the configuration space that connects one structure to the other. Other phonons from the unstable triplet and doublet in the "ideal structure" have frequencies that are higher, but still low enough to suggest that there is a low energy barrier for creation of the domains. In order to clarify these issues, we shall explore the energy landscape around thẽ α andβ structures in more detail in Sec. III C 2. First, however, we begin with a general discussion of RUMs in the cristobalite phases in the next subsection.

Rigid unit mode analysis
Here we analyze the RUMs present in the highsymmetry cubic structure, 47 but constrained to maintain the periodicity of the four-formula-unit (Z=4) cell of thẽ α-cristobalite structure (i.e., containing two unit cells of theβ-cristobalite structure). When these constraints are imposed, we find that there are five linearly independent RUM distortions as shown in Fig. 3. The first three of these distortions, which we denote asβ 1 ,β 2 , andβ 3 , consist of tetrahedral rotations of alternating signs about a single Cartesian axis, and carry the system into the I42d symmetry of theβ phase. The last two, which we denote asα 1 andα ′ 1 , consist of a pattern of rotations around the [110] and [110] axes in the frame of Fig. 1, together with small translations of the tetrahedra needed to keep them connected at their apices, and carry the system into the P 4 1 2 1 2 (or, forα ′ 1 , into the enantiomorphic P 4 3 2 1 2) symmetry of theα-cristobalite structure. Not shown in Fig. 1 are RUM rotationsα 2 andα ′ 2 associated with a second X point, andα 3 andα ′ 3 associated with a third X point.
Within the context of an ideal rigid-unit geometry (in which no additional relaxations are allowed), one can make the following mathematical analysis. The freezing in of the idealβ 1 RUM leads to aβ structure oriented as in Sec. III A, whereas the freezing in of the idealα 1 RUM leads to anα structure as in that section. In the idealβ structure, all five of the modes shown in Fig. 3 remain as true RUMs -i.e., the tetrahedra can remain undistorted to first order in the mode amplitudes. Thus, all five modes are expected to have low frequencies in a more realistic description. However, modesα 2 ,α ′ 2 ,α 3 , andα ′ 3 are no longer RUMs when a finiteβ 1 RUM is present.
In theα 1 structure, only theα 1 ,α ′ 1 , andβ 1 distortions remain as true RUMs. However, theβ 2 andβ 3 modes at least share the same translational symmetry, and so may be expected to have somewhat low frequencies. The remainingα 2 ,α ′ 2 ,α 3 , andα ′ 3 modes are incompatible both in the RUM sense and in their translational periodicity. 48 Not surprisingly, when we impose the translational periodicity consistent with the five modes shown in Fig. 1, we confirm that these five distortions correspond quite closely to the five unstable phonon modes that we found in the ideal structure. The unstable Γ modes correspond toβ 1 ,β 2 andβ 3 , while the unstable X modes correspond toα 1 andα ′ 1 . They also correspond closely to the lowfrequency phonons in theα andβ structures as discussed in Sec. III B 3.
Extending the mathematical analysis of the compati-bility of RUMs, it can be shown that there is an entire three-dimensional subspace of rigid-unit structures (i.e., with the tetrahedral rigidity condition satisfied exactly) in which finite rotations of type (α 1 ,α ′ 1 ,β 1 ) are simultaneously present, and having the space group P 2 1 2 1 2 1 that is induced if any two of them are present. In a similar way, there are two additional three-dimensional (3D) manifolds (α 2 ,α ′ 2 ,β 2 ) and (α 3 ,α ′ 3 ,β 3 ) corresponding to different choices of the X point and thus having different Z=4 supercells. The three subspaces meet only at a single point (the cubic phase with all angles vanishing), and RUMs selected from different 3D manifolds are always incompatible with each other in the sense that the perfect tetrahedral rigidity cannot be preserved when imposing both. This picture has important consequences for our understanding of the possible paths connecting domains of theα andβ structures, as discussed below.

Energy landscape inside 3D manifolds
After we have explained the origin of the low-energy phonons in theα andβ structures by relating them to RUM modes, we would now like to explore the energy landscape around these structures. To do so, we begin by finding a configuration space containing both structures. Since there is no group-subgroup relation between theα andβ structures, we seek a maximal common subgroup of both structures. In the present case, this leads to the space group P 2 1 2 1 2 1 (D 4 2 ). In the P 2 1 2 1 2 1 configuration space, theα andβ structures represent two special points, and we know that the energy has local minima at these points because all computed phonon frequencies were found to be positive there. But then we also expect that there must be at least one saddle point connecting these points. To search for this saddle point, we started from the midpoint between thẽ α andβ structures in the 12-dimensional P 2 1 2 1 2 1 configuration space (described by nine internal coordinates and three cell parameters), and identified the unit vector e pointing between the two structures. We then applied a simple saddle-point search strategy in which component of the force vector parallel toê was reversed in sign before executing the steepest-descent update. This algo- rithm can be expected to succeed if the saddle point is not too far from the midpoint and if the principal axis of the negative Hessian eigenvalue at the saddle point is roughly parallel toê. In the present case, it led to a rapid convergence on the desired saddle point. Surprisingly, we find that the saddle point has a very low energy, only 5 meV per formula unit above that of theα structure, or 17 meV above that of theβ structure. The three points representing theα andβ structures and the saddle point determine a plane in the 12dimensional configuration space. To confirm that the path running through the saddle point encounters only a single barrier, we have plotted the structural energy in this plane (without relaxation of other coordinates) in Fig. 4. We have somewhat arbitrarily carried out a linear transformation on the coordinates in such a way that theα andβ structures lie at (0, 0) and (1, 0) respectively, while the saddle point lies at (0.5, 1), in Fig. 4. The results confirm the picture of a simple barrier of 5 meV encountered when going from theα to theβ structure.
Note that a transformation that would lead from thẽ α to theβ structure along a straight line in configuration space would have an enormously higher barrier of 195 meV per formula unit. This is because the straightline path is a poor approximation to a RUM. If instead we follow a curved minimum-energy path fromα through the saddle toβ and compute the relaxed Si-O bond lengths and O-Si-O bond angles along this path, we find that these remain almost constant. This strongly suggests that this minimum-energy path may be well approximated by some RUM-like distortion.
In the Sec. III C 1, we pointed out that within the framework of ideal rigid-unit rotations, there is an entire three-dimensional subspace of structures for which the tetrahedral rigidity conditions are satisfied exactly, in which finite rotations of all three types are present. We label an arbitrary configuration in this 3D manifold by (α 1 ,α ′ 1 ,β 1 ), where by convention the order of opera- tions isα 1 followed byα ′ 1 and thenβ 1 . The space group at a generic point in this 3D manifold is P 2 1 2 1 2 1 , the same one we have just been discussing. It thus seems likely that the minimum-energy path in Fig. 4 may correspond approximately to a path from the point (α 1 ,0,0) to the point (0,0,β 1 ) and lying, at least approximately, in the two-dimensional (2D) subspace (α 1 ,0,β 1 ).
To test this conjecture, we first created an ideal rigidunit structure for each pair of angles (α 1 , β 1 ) on a twodimensional mesh. We then used our first-principles calculations to relax each structure subject to the constraint that these two angle variables should not change. Technically, we did this by carrying out the minimization of the energy in the ten-dimensional subspace orthogonal to the two-dimensional surface for each starting point (α 1 , β 1 ). We typically found that these relaxations were small, confirming the approximate validity of the RUM picture.
The energy surface determined in this way is plotted as a function of rotation angles α 1 and β 1 in Fig. 5. The minima corresponding to theα structure are immediately visible near the left and right sides of the figure, while those corresponding to the slightly lowerenergyβ structure appear near the top and bottom. The minimum-energy path appears to be roughly circular on this plot, and four equivalent saddle points are apparent at α 1 ≃ ±19 • and β 1 ≃ ±13 • . These saddle points are equivalent to the one identified in Fig. 4, with a bar-rier height of 5 meV per formula unit relative to theα structure. We thus confirm the presence of a very lowenergy barrier between these structures, and identify it as approximating a certain path in the space of rigidunit rotations. A video animation showing the evolution of the structure along this path is provided in the supplementary material. 49 It is important to note that, according to the simplified model of Eq. (1) of Ref. [7], the energy would remain exactly zero on the entire (α 1 , β 1 ) surface of Fig. 5 since the ideal rigid-unit structures satisfy the rigidity conditions analytically. The RUM framework envisages extensions to make the model more realistic; one way to do this is by adding an energy term that depends on the relative tilts of neighboring tetrahedra. 50 We tried this by introducing a simple double-well potential model that penalizes the departure of the Si-O-Si bond angles from a preferred bending angle. In this model the change of total energy per formula unit is where the sum runs over all N Si-O-Si bond angles φ i in the unit cell. We found that we could obtain an optimal fit 51 to the results of our first-principles calculations using parameters E 0 = 83 meV per formula unit and φ 0 = 145 • . The energy landscape of the fitted model looks very similar to the results plotted in Fig. 5. In particular, the overall circular aspects of the energy landscape and minimal-energy path in Fig. 5 are reproduced. However, the model unfortunately assigns identical energies to theα 1 andβ 1 structures, and moreover predicts a path connecting them on which the energy remains completely flat. This happens because, for any given pair of angles (α 1 ,β 1 ) on or near this path, one can find a small α ′ 1 such that the rigid-unit structure (α 1 ,α ′ 1 ,β 1 ) has all its Si-O-Si bond angles exactly equal to φ 0 . Therefore our simplified model of Eq. (7), or any other model that depends solely on the Si-O-Si angles, predicts a zero-barrier path betweenα 1 andβ 1 structures. This behavior is reminiscent of an early model of Nieuwenkamp 52 for β cristobalite, in which the Si-O-Si bond was assumed to rotate freely on an annulus lying in the plane that is equidistant between Si atoms.
To further test the model of Eq. (7), we performed firstprinciples calculations on a mesh of (unrelaxed) structures in the 2D space (α 1 ,α ′ 1 ,β 1 =0). The resulting energy landscape is shown in Fig. 6. The barrier betweenα 1 and α ′ 1 structures is now about 35 meV, substantially higher than for the path connectingα andβ structures. When we use the same fitting parameters obtained earlier, we again get very good overall agreement; the energy landscape obtained from our model has the same diamondlike appearance as in Fig. 6, and saddle points appear in very similar locations. Moreover, the barrier of 41 meV predicted by the model is in quite good agreement with the first-principles value of 35 meV. However, in this case Energy as a function of rotation angles φα 1 and φ α ′ 1 , corresponding to rotationsα1 andα ′ 1 shown in Fig. 3. The origin corresponds to "ideal" cristobalite. The four minima (filled circles) correspond to theα ′ 1 structure (top and bottom) and to theα1 structure (left and right). The energy difference separating contours is 3 meV per SiO2 formula unit. the picture presented by Fig. 6 is somewhat misleading, because it turns out that the entire minimum-energy path lying in the β 1 =0 plane is unstable, and falls to lower energy asβ 1 is turned on. Thus, the apparent saddle points in Fig. 6 are actually stationary points with two negative eigenvalues in the 3D (α 1 ,α ′ 1 ,β 1 ) space. Within the model of Eq. (7), in fact, the lowest-energy path connecting theα 1 andα ′ 1 structures is actually completely flat, being composed of a segment connectingα 1 toβ 1 and then another connectingβ 1 toα ′ 1 . This observation agrees with our first-principles calculations, since if we start from the purported saddle-point configuration and do a structural relaxation subject to the constraint that α 1 = α ′ 1 , the structure is found to converge to theβ 1 structure as expected.
implies that applied pressure would tend to favor theα phase and raise the α-to-β transition temperature. At the saddle point, we find that the volume per formula unit is 46.8Å 3 , which is just slightly larger than for either of the parent-phase structures. This finding may be of interest for future studies of the pressure-dependence of the phase-transition mechanism.

Domain walls
The barriers discussed in Sec. III C 2 refer to transformation pathways in which the crystal remains periodic and transforms homogeneously, and the energy barriers are given per unit cell. It would also be of interest to consider the energies of domain walls between variousα and β structures. This is beyond the scope of the present investigation, but the results for homogeneous transformations may give some hints as to what could be expected. For example, we speculate that domain walls connecting α andβ structures belonging to the same 3D rigid-unit manifold will probably have a rather low energy per unit area, while those connecting structures belonging to different 3D manifolds would be expected to have much higher energies.

IV. DISCUSSION
In this section we give a brief overview of several previously proposed models of α and β cristobalite phases, and discuss how the results of our calculations relate to those models.
The RUM model of Ref. [7] describes the β phase as an average cubic structure that has strong dynamical fluctuations occurring simultaneously into RUMs in all allowed regions of the Brillouin zone. A simplified version of this picture would be one in which the tetrahedra are assumed to be completely free to pivot around their apices, as in Eq. (1) of Ref. [7]. In general, the simultaneous excitation of more than one RUM will have an associated energy cost because the tetrahedra typically cannot remain perfectly rigid while undergoing both kinds of distortion simultaneously. However, as an exception, we have identified 3D rigid-unit manifolds within which the geometrical constraints can simultaneously be satisfied. Within the model of Eq. (1) of Ref. [7], or the split-atom model of Ref. [54], the energy landscape within this special 3D manifold would be completely flat, and one would expect that freezing in of one RUM of typeα 1 ,α ′ 1 orβ 1 would have no consequence on the energy profile of one of these other RUM distortions.
However, once one goes beyond the simplest versions of the model and includes terms that depend on the Si-O-Si bond angles at the apices, our calculations indicate that the RUM distortions of typeα 1 ,α ′ 1 andβ 1 become coupled and have a rich energy landscape. As a step in this direction, the more sophisticated split-atom model hav-ing an additional energy term depending on Si-O-Si bond angles 50 should provide an improved description. However, we note that even this model, or any model based solely on Si-O-Si bond angles, still has a nonphysical behavior in that it would necessarily predict zero-energy barriers between theα 1 andβ 1 structures, as discussed at the end of Sec. III C 2. Nevertheless, we believe that the split-atom and similar models can provide important complementary information to ours, since they are not restricted to periodic supercell structures as ours are.
Among the models of cristobalite phase transitions is also the model of Hatch and Ghose. 55 They argue that the β phase is dynamically and spatially fluctuating between the 12 different possibleα domains having P 4 1 2 1 2 space-group symmetry. The counting arises because there are three different X points; each exhibits a doublet of degenerate modes leading to enantiomorphic α 1 andα ′ 1 structures (see Fig. 3); and the tetrahedra can rotate by ±φ. The model is based on symmetry arguments and assumes that all of the barriers separating these 12α structures are small. However, our work suggests that the barriers separating different types of α domains have very different barriers. Furthermore, their model does not take into account the fact that thẽ β structure is easily accessible with a very low barrier, suggesting that fluctuations into theβ structure may be more important than some of the otherα structures.
Finally, O'Keeffe and Hyde 4 do discuss a path connectingα andβ structures, but it is of a different type than those discussed above since it connectsα andβ structures belonging to different 3D rigid-unit manifolds. In our notation, their path would connectα 1 orα ′ 1 toβ 2 orβ 3 , etc. Such a path would involve the simultaneous application of RUM rotations that are incompatible with each other, and as such would be expected to have a high energy barrier.
To clarify our view of the cristobalite phase transitions, we start by emphasizing once again the existence of three distinct 3D rigid-unit manifolds, as described above at the end of Sec. III C 1. To review, one of these is described by rotation angles (α 1 ,α ′ 1 ,β 1 ) giving rise to structures of space group P 2 1 2 1 2 1 whose translational periodicity is that corresponding to the X point (2π/a c )(001) or equivalently (2π/a c )(110). This manifold contains theα 1 ,α ′ 1 , andβ 1 structures, and their partners with reversed sense of rotation, as shown schematically in Fig. 7. The second and third 3D subspaces are described by rotations (α 2 ,α ′ 2 ,β 2 ) and (α 3 ,α ′ 3 ,β 3 ), with periodicities set by X points (2π/a c )(010) = (2π/a c )(101) and (2π/a c )(100) = (2π/a c )(011), respectively. We have found that these three subspaces are essentially incompatible, in the sense that it is not possible to combine rotations taken from any two of them into a combination that preserves the rigid-unit constraints. This occurs in part because these three 3D rigid-unit subspaces have incompatible translational symmetries, but also because of incompatibilities in the patterns of rotations.
The structure of the space sketched in Fig. 7 is intended to reflect a three-level hierarchy of energies and energy barriers as suggested by our analysis. In the model of Eq. (7), the energy is degenerate for all six of the structures shown in Fig. 7, as well as on the solid curves connecting them. 49 According to our first-principles results, this picture is modified so that theα andβ structures are local minima, with low-energy saddle points (∼5 meV) between them (see Fig. 5). The low curvature of the energy surface along these curves is reflected in the presence in Table VI of a very soft 29 cm −1 B 1 mode starting from theα structure, and a 35 cm −1 M 3 M 4 doublet starting from theβ structure. 56 While our calculations have thẽ α structures at a slightly higher energy than theβ ones, this is presumably reversed in the true physical system. The next energy scale in the hierarchy is that associated with the direct paths betweenα structures in the same 3D manifold, indicated by the dashed lines in Fig. 7. As shown in Fig. 6, this energy is on the order of ∼35 meV, so that the true minimum-energy path between neighboringα structures goes instead through (or perhaps nearly through) theβ structures.
Finally, the highest energies are associated with the barriers separating any of the structures in Fig. 7 from any of the structures in the other two 3D subspaces. These barriers are on the order of 80 meV, the energy needed to pass through the undistorted cubic phase. While not enormously larger than the 35 meV mentioned above, this is high enough that we do not expect these barriers to be especially relevant for the phase transitions in this system.
We can now speculate on the nature of the phase transition between α and β cristobalites. We propose that in the lower-temperature α phase, the system is locally frozen onto one of the minima of typeα in one of the 3D manifolds, but with substantial fluctuations along the low-energy paths leading to the two neighboringβ structures in the same manifold. Then, in the highertemperature β phase, we speculate that the system instead shifts over and condenses locally onto one of thesẽ β structures, but with substantial fluctuations along the low-energy paths leading to the four neighboringα structures, all in the same 3D manifold. The fact that there are four low-energy paths to fluctuate along, instead of two, is consistent with the fact that the β phase (being the higher-temperature phase) has higher entropy. If the system were truly to freeze onto a singleβ structure, it would be globally tetragonal, with space group I42d. However, it is also possible that the system forms on some larger scale into spatiotemporal domains composed ofβ structures from all three of the 3D manifolds, giving an overall average F d3m structure in accord with the picture espoused in Refs. [7,8].
Let us return for a moment to the recent work of Zhang and Scott, 6 who argued that their Raman studies of single crystals of β-cristobalite were inconsistent with O h symmetry. Assuming D 2d symmetry instead for the β phase, these authors then noted that D 4 is not a subgroup of D 2d , and thus that the existence of a groupsubgroup relation for the phase transition would rule out the assignment of the α phase to the D 4α structure. On this basis, they suggested that a lower symmetry, such as D 2 , should be considered for α-cristobalite. Our view, instead, is that a group-subgroup relation does not have to hold for the transition, since the transition is known to be of first order, and thus assignments of D 2d and D 4 for the α and β phases respectively are not inconsistent. As pointed out in the introduction, while certain spectroscopic signatures of the transition are indicative of a weakly first-order transition, the volume change and latent heat at the transition are substantial. The transition may perhaps be described as a reconstructive transition in the sense of Tolédano and Dmitriev, 11 although in the present case the rearrangements of atoms can occur very gently, because of the existence of very low-barrier paths of D 2 symmetry connecting the D 4 (α) and D 2d (β) structures. The situation may be somewhat analogous to the tetragonal-to-orthorhombic and orthorhombic-torhombohedral transitions in ferroelectric perovskites such as BaTiO 3 and KNbO 3 , where the presence of low-barrier paths of monoclinic symmetry is associated with the weakly first-order nature of the transitions. 57 Unfortunately our calculations are carried out at 0 K with crystal periodicity imposed. It is therefore difficult to draw any firm conclusions about the nature of the phase transitions between cristobalite phases, especially if fluctuations are as important as we think they are, and much of what we have said above must remain speculative. Nevertheless we hope that the results of our calculations will be of use in developing improved models that may allow for realistic finite-temperature modeling of the phase transitions in this system, ultimately leading to a resolution of the controversies that have surrounded this system over the years.

V. SUMMARY
Based on first-principles calculations, we have performed a detailed analysis of theα (P 4 1 2 1 2) andβ (I42d) structures of cristobalite SiO 2 . In particular, we have confirmed that both structures are locally stable against all possible distortions associated with Γ-point modes of the four-formula-unit conventional cell. We have calculated phonon frequencies for theα andβ structures, compared these to the experimental values, and discussed how the phonons in these two structures are related to each other. We have also tried to resolve some experimental anomalies that were found in spectroscopic studies of the cristobalite phases. Finally, we have explored the energy landscape connecting theα andβ structures. We have emphasized the existence of three distinct 3D manifolds of structures, each of which contains bothα andβ structures that can be connected to each other within the manifold by paths with a surprisingly small barrier of 5 meV per formula unit, while paths connecting different manifolds have a much higher barrier. While our calculations do not properly treat fluctuations, we nevertheless have speculated on the possible consequences of our findings for the understanding of the α-β phase transition in cristobalite SiO 2 .
In a picture in which the rigid-unit constraints are perfectly enforced, it turns out that the structures ofα and β symmetries are completely determined by a single parameter, which can be taken to be the volume V per formula unit relative to the corresponding value V 0 in the ideal cubic structure. (That is, V 0 is the volume below which rigid distortions start to appear, as explained in Sec. III A.) In this Appendix, we check to see how closely our structures, as optimized from the first-principles calculations, match with this elementary model.
The solid curves in Figs. 8 and 9 show the mathematical predictions of this elementary model, obtained by applying rotations of typesα 1 andβ 1 (see Fig. 3) in such a way as to keep the tetrahedra perfectly rigid. (For V > V 0 , the elementary model cannot be satisfied, and the ideal cubic parameters are plotted instead.) The symbols shown in Figs. 8 and 9 denote the results of our first-principles calculations where, for each specified value of V , the volume was treated as a constraint while all other structural parameters were relaxed. The fit was optimized by choosing a common V 0 = 55.1Å 3 for both α andβ structures. For reference, the first-principles equilibrium volumes are 45.7 and 46.7Å 3 for theα and β structures, respectively. We find that the agreement is extraordinarily good for all of the internal parameters, but that there are some significant discrepancies in the c/a ratios. At first sight this may seem contradictory: why are the c/a ratios off by many percent, while the Si-O bond lengths agree within ∼0.05%? The answer is connected with the presence of volume-preserving tetragonal distortions of low energy cost. In such a distortion, each tetrahedron is stretched slightly along c and compressed in a (or vice versa), and it happens that the tetrahedral angle of arccos(1/ √ 3) is precisely the one at which Si-O bond lengths are preserved to first order in the distortion amplitude. While the O-Si-O bond angles do change at first order, this may entail a smaller energy cost than for bond-length changes. As expected from this analysis, we find that our first-principles O-Si-O bond angles differ from the ideal ones by ∼4%. In short, it appears that it is energetically more important to preserve bond lengths than bond angles, and that for geometrical reasons this translates into an enhanced freedom for the c/a ratio to deviate from the ideal rigid-unit geometry.