The considerable SiO₂ enrichment that is the predominant feature of magma differentiation is in fact paradoxical. The first crystals to precipitate in a multicomponent system are the most refractory, i.e., they have the highest melting points. In accordance with Lindemann's (1910) early ideas, they are assumed to have the strongest interatomic bonds because melting is usually associated with bond breaking caused by mean square vibrational amplitudes exceeding some critical fraction of bond distances. Now, SiO bonds are commonly supposed to be the strongest in silicate systems because of the small radius (r) and high nominal electrical charge (z) of the Si⁴⁺ ion compared to those of other major cation-building minerals (Fig. 2). On the basis of ionic forces, one would expect SiO₂ to melt at a much higher temperature than CaO and MgO, but the reverse holds true. Whereas the melting temperature of cristobalite, the high-temperature polymorph of SiO₂ is 1999 K, CaO and MgO melt at temperatures higher than 3000 K. As a consequence, one would expect the silica content to decrease progressively upon fractional crystallization, which clearly contradicts the evolutions that take place along the petrogenetic series.

To account for this anomaly, we will first consider the closely related question of the α–β transitions in tectosilicate minerals, which are accompanied by marked changes in macroscopic properties such as volume (Fig. 3). It is now agreed upon that these transitions represent the change from a static to a dynamic structure when temperature increases. For instance, the Raman-active modes of vibration of cristobalite and carnegieite (NaAlSiO₄), which is structurally related to cristobalite, undergo marked modifications at the transition temperatures of 540 and 970 K, respectively (Fig. 4). Then they vary little upon further temperature increases and even upon melting. These modifications are especially clear at low frequencies where the vibrational modes, which involve mainly intertetrahedral motion, merge to form single, broad envelopes similar to those of liquids. Actually, the process at work is free rotation of rigid SiO₄ tetrahedra caused by precession of SiOSi linkages with respect to their average orientation (Dove et al., 1997). Interestingly, the entropy changes associated with the α–β transition of these minerals represent more than half and two thirds of the melting reactions of cristobalite and carnegieite, respectively (Table 1). In other words, setting in partial atomic mobility at the α–β transition through motion of oxygen atoms has an entropy cost comparable to disposing of long-range order at the melting point.

Lattice and molecular dynamics simulations have given a more detailed picture of these transitions for SiO₂ polymorphs (Bourova et al., 2000, 2004). The salient point is that, beyond the scale of SiO tetrahedra, β-cristobalite has a disordered structure at all length scales, which is a dynamical average of three substructures. One of them (α₁) has the same chain arrangement of SiO₄ tetrahedra as α-cristobalite in its tetragonal symmetry (space group P4₁2₁2) with the same value for the intertetrahedral SiOSi angle. The two other substructures α₂ and α₃ have different chains of SiO₄ tetrahedra that are derived from α₁ through a 120° rotation around the <111> axis. Because the three substructures have slightly different free energies, the transition from one to the two others is ensured by the thermal energy of oxygen atoms. The motion of a given oxygen atom takes place mainly along directions perpendicular to other OO bonds, which correspond to the first OO distance. Because of the rigidity of the SiO bond, the vibrational amplitude along these directions is twice larger than in the plane along the Si-Si direction and it keeps increasing with increasing temperature and thermal energy. As for the unit-cell contraction signaled at high-temperatures by the negative thermal expansion coefficient (Fig. 3), it results from the decrease of the first Si-Si distance associated with the large thermal motion of oxygen atoms, which cause the neighboring tetrahedra to come closer to each other.

Although a similar explanation is arrived at for the α–β transition of quartz, an important difference is that the denser arrangement of SiO₄ tetrahedra in quartz compared with cristobalite makes room for only two substructures with similar energies. These are analogous to the Dauphiné twins of α-quartz (α+ and α−), which differ in the orientation of SiO₄ chains around the <001> direction. Turning now to coesite, the densest form of silica with four-fold coordinated Si, we just note that, as experimentally observed (Fig. 3), the monoclinic structure is metastable (with respect to both quartz and SiO₂ liquid) up to the highest temperatures investigated in lattice and molecular dynamics simulations. Thermal expansion is small because a high density limits the extent of rotation of SiO₄ tetrahedra. Above about 1000 K, the structure of coesite nonetheless becomes dynamically disordered and similar to those reported for the β-phases of quartz and cristobalite. Disorder sets smoothly, however, in contrast to its abrupt onset in quartz and cristobalite at their transition temperatures. The radial distribution functions for all bond distances indicate that order then prevails only for the nearest neighbors whereas the angle distributions widen markedly so that the monoclinic form of coesite with a SiOSi angle of 180° is just a time-averaged structure. Because all oxygen displacements are strongly anisotropic in both β-quartz and β-cristobalite, disorder in coesite is only partial and affects primarily intrachannel oxygen.

The differences between the three silica polymorphs are sketched in terms of potential energy wells on Fig. 5 with the usual assumption that pressure mostly affects the attractive part of interatomic potentials. Energy minima are thus shifted toward the shorter distances at which strained bonds would match the externally imposed forces. At room temperature the potential of cristobalite has three wells but a single minimum. For quartz, the shift of the energy minima toward higher energies at higher pressure results in the existence of only two wells, which have similar energy minima. At the still higher pressure of formation of coesite, there finally remains a single energy minimum. This sketch has another interesting feature. If one uses Landau's criteria δ²ϕ/δr² = 0 and δ³ϕ/δr³ > 0 for equilibrium phase transitions (taking the distance r as the relevant order parameter), one sees on Fig. 5 that the temperatures of these transitions decrease from cristobalite to coesite in the order of the reported melting points. For coesite, the melting point would in particular be lower than the temperature of the onset of dynamical disorder and will remain much lower than would be guessed from ionic bonding.

In summary, the anomalously low melting temperature of SiO₂ polymorphs originates in the fact that significant atomic mobility can occur in open three-dimensional networks without requiring SiO bond breaking. Another manner to view this effect is to invoke the profound contrast between the high rigidity of Si–O bonds and the high flexibility of SiOSi angles, which was well illustrated by early quantum-mechanical calculations made for H₆Si₂O₇ clusters (Gibbs et al., 1981). As shown on Fig. 6, similar energy variations are caused by 0.02 Å changes of SiO bond lengths and 20° variations of SiOSi angles. As a matter of fact, this result represents a simple justification for the rigid unit model (e.g., Swainson and Dove, 1993) as applied to silicates whereby SiO₄ tetrahedra are taken as rigid units between which tilting results from precession of SiO bonds around their average orientation. With increasing vibrational amplitudes, the almost constant length of the very strong SiO bonds then imposes changes in either SiOSi angles or Si-Si distances.

Of course the directionality of SiO bonds becomes less important with decreasing SiO₂ content and, thus, with depolymerization of the silicate network so as to vanish for the isolated tetrahedra of the orthosilicate stoichiometry. This strong bond directionality has nonetheless fundamental implications for atomistic simulations of the structure and properties of crystalline or amorphous silicates. For ionic substances, a common starting point in theoretical simulations is for instance the Born–Mayer potential:

Silicon and carbon belong to the same IV-A group of the periodic table. Both have four electrons in their p shell so that their common tetrahedral coordination has long been explained by the formation of four hybrid orbitals sp³ with either other S or C atoms or with different atoms such as oxygen for Si and hydrogen for C. In silicates, the p shell of silicon is thus filled by outer electrons of the four oxygen donors, which are each directed toward one of the apexes of a regular [SiO₄]⁴⁻ tetrahedron. The differing principal quantum number n notwithstanding, the analogy between the bond strength-bond energy relationships in the s-p combination for carbon and silicon is almost complete. For both elements, the maximum bond strength (i.e., the angular dependence of the magnitude of a bond orbital) is observed not far from an s:p = 1:3 overall proportion. Although the sp³ terminology implies that the weights of the respective s and p components of a hybrid sp³ orbital are indeed in a 1:3 ratio, this condition for bond formation is neither strict nor unique. In addition to the appropriate directionality of the atomic orbitals (AOs), two other conditions have to be satisfied:

• • overlapping of discrete AOs on different atoms bonded in the species along the appropriate direction;
• • delocalization of electrons between adjacent basins.

Despite the limited basis set of atomic orbitals used in early calculations, correct energies were obtained for molecular orbitals by Collins et al. (1972), who found a considerable participation of the 3d orbital of silicon in the valence orbitals of oxygen to form π-bonding e orbitals and π-bonding t₂ orbitals in the valence region. Interestingly, Collins et al. also noted that the [SiO₄]⁴⁻ ion cannot exist in the gaseous state because of the presence of 24 electrons in unbound molecular orbitals. Since then the great many ab initio calculations performed for SiO bonding have not changed much the general picture of molecular orbitals in SiO₂. New Density Functional Theory (DFT) calculations for instance confirm former conclusions about the inner orbitals (orbital coefficients and ensuing character assignation, and energy; see Table 2). Rather than reviewing these features, we will just stress two significant points.

The first is the basic role of the medium permittivity. If the [SiO₄]⁴⁻ group is placed in a reaction field with a permittivity ɛ of 3.85, which is close to that of molten silica, important changes are observed, at the same theory level, in the orbital populations and energies of both valence and inner regions (Table 2) such that the group becomes “stable” to all extent. If we then perform a Mulliken population analysis on our unstable gaseous ionic cluster we obtain a charge of 0.252 e localized on the central Si and of −1.063 e on each of the four bonded oxygen. But when the ions are placed in the dielectric medium the charge localized on Si increases to 1.917 e and that of the four bonded oxygen to −1.479 e. As defined by Stewart et al. (1980) in terms of formal and residual charges, the fractional ionic character of the bond is reduced to 0.063 for the isolated cluster (in simple terms, an essentially covalent ionic cluster), whereas f_i = 0.479 in the dielectric medium. The permittivity of silica thus has a strong stabilizing effect on [SiO₄]⁴⁻ groups.

The second point concerns bonding between adjacent [SiO₄]⁴⁻ groups and the ensuing anomalously low melting temperature of the SiO₂ polymorphs discussed in the previous section. The particular strength of Si–O bonds has the consequence that electron delocalization associated with polymerization through creation of SiOSi linkages can be expected neither to be large nor to impose strict geometrical constraints on the arrangements of the [SiO₄]⁴⁻ building blocks. In other words, the form of MOs eventually encompassing adjacent atomic centers is not per se indicative of a bonding condition. The Extended Hückel Molecular Orbital calculations aimed at establishing the extent of bond overlap population n(SiO) and its effect on the cluster geometry carried out by Louisnathan and Gibbs (1972) on both the silicate ion [SiO₄]⁴⁻ and the orthosilicic acid group [H₄SiO₄]° showed a negligible effect of bond overlap on angular distortion. More advanced calculations in terms of electron localization/delocalization through Electron Localization Function or the Laplacian of the electron density (see Gatti, 2005; for an extensive treatment) have been more recently carried out by Gibbs et al. (2003) on geometry optimized H₆Si₂O₇ molecule (i.e., for the structure yielded by the minimum of the energy hypersurface) and in coesite to show the extent of bonding and non-bonding regions when one varies the SiOSi angle between 110 and 180°.

Of course, silicon is not the only network-forming cation in silicate crystals and melts. The most important other element to be considered in this respect is aluminum. But for Al³⁺ to substitute for Si⁴⁺ the deficit of the nominal 3+ charge of Al ions must be compensated, as done in feldspars, by association with either an alkali or half an alkaline earth. When different alkali and alkaline earth cations are present, however, the question arises to determine which ones are the most effective to ensure charge compensation of nominally trivalent aluminum in tetrahedral coordination. To give a firmer theoretical basis to this question ab initio calculations have been made on the association with various clusters of SiO₄ and AlO₄ tetrahedra of Na, K, Ca and Mg, the ad-atoms (i.e., atoms initially placed at a certain arbitrary distance from the cluster). In a first study (Gatti et al., 2012), we considered 4- to 6-membered rings in which either one or two Al atoms were substituted for Si. Specifically, the energy transfer between the free ion M was investigated in terms of a homogeneous gas-phase reaction of the type:

The structure of each complex can be determined from the true minimum of the potential energy surface and its energy computed with reference to that of its Al-free counterpart. In this respect, thermal contributions to the thermodynamic functions are negligible compared to the total electronic energy. Barring some stabilizing effects of vibrational motion, the calculated energy changes thus simply reflect electronic energy transfers. For 4-membered rings, the geometry remains largely unaffected by the presence of the ad-atom, the position of the latter within the ring nonetheless depending on its nature (Fig. 7a). More complicated effects are observed in larger rings for which the easily bent SiOSi angles allow greater flexibility and, thus, more deformation (Fig. 7b). From a geochemical standpoint, energy calculations have the important consequence that they would make calculations of equilibrium constants possible for the exchange reactions between the various Al-M complexes. Instead of assuming that only one or two M cations act as charge compensators, the actual contribution of each possible cation could then be predicted along with its effects on relevant physical properties.

For this purpose, however, it is necessary to perform accurate electron localization/delocalization analyses to determine how the valence electrons redistribute upon alkali or alkaline earth complexation, i.e., how the electrons of the various molecular orbitals redistribute among the adjacent basins of the ring network (Gatti et al., 2012). As found earlier for crystal-like aluminosilicate complexes (Vayssilov et al., 2000), one observes that the stabilization of O1s orbitals in the SiAl cluster exceeds by far the amount of electronic energy gained in complexation, which indicates that a substantial redistribution is operating throughout the complex. The Al-free 4-membered ring may be interpreted as stabilizing a Pauling double bond because it has a low-lying molecular orbital essentially constituted of the s orbitals of oxygen with some contribution from s orbitals of silicon. In a 3-d representation, the electron cloud of this σ-bonding molecular orbital appears delocalized on the four oxygen donor centers on the xy plane of the ring when arranged in a standard configuration.

In 6-membered Al-free rings an orbital of nearly identical energy is observed. This orbital is essentially composed of the s AO of all oxygen donor centers, with minor contributions from the Si atoms. The substitution of one or two Si atoms of the ring with Al has some perturbation effect on the density distribution. Whereas in the case of the silico-aluminate 4-membered ring the apparent delocalization still operates on all donor centers, in the monosubstituted 6-membered ring the delocalization excludes the two donor centers nearest to Al. The Pauling double bond is apparently present in linear chains as well, such as the [Si₂O₇]⁶⁻ dimers. Also in this case all oxygen centers in the cluster are surrounded by the electron cloud, which excludes the central metals (i.e., not only the bridging oxygen but also the singly bonded ones). Further calculations are required to check this point for clusters with different geometries. But the facts that the same Al-compensation mechanism seems to operate over the whole geochemically range of melt composition, and to affect linear chains as well, suggest that the effect is fairly general.

When dealing with volatile solubility, electrostatic, dispersive (van der Waals) and repulsive solute-solvent interactions must be taken into account. All three are quantum-mechanical in origin and may be addressed theoretically with the appropriate choice of functional theory level and basis set. In addition, a fourth factor must be considered, namely, the cavitation energy of the solute, i.e., the energy needed to form cavities in the solvent appropriate for hosting solute species. In physical chemistry these contributions are so well known that solvent properties (molecular radius, dielectric constant, etc.) are conformed to the observed solubility of some noble gas (e.g., He in water or benzene) in the most widely used quantum-mechanical computational packages.

In geochemistry, however, one generally assumes that species such as noble gases do not interact with the aqueous or silicate species among which they are hosted in the available cavities of the solution (e.g., Bouhifd et al., 2008; Chamorro-Perez et al., 1998). This assumption is actually too simple so that it prevents the solubility of these “inert” gases in aqueous and silicate solutions from being understood at a fundamental level. Actually, the bulk Gibbs free energy of solution for a gaseous species in its stable standard state is simply related to the Henry's Law constant K_H by:

Although quantum-mechanical dispersion effects thus have to be addressed, the main factor in eq. (3) is the cavitation energy on which we will thus focus here. The diameter of the solvent is a thermally averaged quantity over various orientations of the solvent with respect to the cavity surface. If we denote it by σ_s, the problem is to evaluate the reversible work W₀ necessary to create in the solvent a cavity of radius r≤σs/2. To solve this problem, we used the formulation based on the scaled particle theory (SPT) of Reiss and coworkers (Reiss et al., 1960) with which Pierotti (1976) obtained:

To make room in the solvent for a solute of diameter σ_i the necessary cavity must have a radius r=σs+σi2. Denoting by Y=πρsσs6 the reduced number density and by ς=σi/σs, Pierotti (1976) rewrote the work as:

With this formalism, one finds that the solvent radius σ_s is satisfactorily reproduced by a simple expression of the form σS=2×1.1383rcoll−0.3836 where r_coll is the collisional radius (i.e., rcoll=V/83; Stearn and Eyring, 1937). The striking feature is that this expression fits nicely the observed solubility trends for noble gases in a variety of solvents that include non-polar (C₆H₆) and polar (H₂O) liquids. The major feature of this rather simple formalism thus is that a volatile species does not dissolve in pre-existing cavities, but in cavities that are formed by interaction with the solute. Assessment of noble gas solubility in various fluids then is a simple but stringent test of the reliability of the parameterization adopted to model the Gibbs free energy of cavitation. Based on the wealth of solubility experiments carried out by geochemists through the decades, it is possible to retrieve not only the thermodynamic parameters of the various solutes, but also the solvent properties needed to determine the solubilities observed at high-temperatures and pressures through inverse methods (Ottonello and Richet, 2014).

Since the physical and chemical properties of macroscopic bodies are determined from those of their atomic constituents, fundamental insights on these properties have to be derived from accurate descriptions of bonding at an atomic level. To illustrate the kind of work to be done to understand the Earth's structure and dynamics at a really fundamental level, this review has dealt with a few important topics but is clearly not comprehensive; volatile dissolution or phase diagram calculations are just two examples among the great many other features that can be investigated. In this respect, it is well known that, over experiments, calculations have the great advantage of being performed as readily under extreme conditions as under room pressure and temperature. The composition ranges studied can also be varied at will to address the very large diversity of natural conditions. Even though the size of the system that can be investigated is still limited in terms of number of atoms and, thus of structural complexity, continuous improvements in computer power will continue to alleviate this problem.