Pressure induced elastic softening in framework aluminosilicate- albite (NaAlSi3O8)

Albite (NaAlSi3O8) is an aluminosilicate mineral. Its crystal structure consists of 3-D framework of Al and Si tetrahedral units. We have used Density Functional Theory to investigate the high-pressure behavior of the crystal structure and how it affects the elasticity of albite. Our results indicate elastic softening between 6–8 GPa. This is observed in all the individual elastic stiffness components. Our analysis indicates that the softening is due to the response of the three-dimensional tetrahedral framework, in particular by the pressure dependent changes in the tetrahedral tilts. At pressure <6 GPa, the PAW-GGA can be described by a Birch-Murnaghan equation of state with  = 687.4 Å3,  = 51.7 GPa, and  = 4.7. The shear modulus and its pressure derivative are  = 33.7 GPa, and  = 2.9. At 1 bar, the azimuthal compressional and shear wave anisotropy  = 42.8%, and  = 50.1%. We also investigate the densification of albite to a mixture of jadeite and quartz. The transformation is likely to cause a discontinuity in density, compressional, and shear wave velocity across the crust and mantle. This could partially account for the Mohorovicic discontinuity in thickened continental crustal regions.

Plagioclase feldspar is one of the most important mineral solid-solution series between the end-members-albite (NaAlSi 3 O 8 ) and anorthite (CaAl 2 Si 2 O 8 ). Plagioclase is ubiquitous in crustal rocks in the Earth and terrestrial planetary bodies, including Moon 1 . Sodic-plagiocase or albite is a major constituent of several rock types includinggranite, granodiorite, diorite, tonalite, and basalt. The modal abundances of plagioclase in these rocks vary between 16 and 60% 2 . In hydrothermal 'Kluftalbite' veins the modal abundance of albite rich plagioclase could be nearly 100% 3 . In subducted oceanic crust, albite is widely distributed in metamorphosed mid-ocean ridge basalt (MORB) i.e., zeolite to garnet granulite facies rocks with a modal abundance ranging between ~18-22% 4 . Albitization is a common metasomatic process in granitic and mafic lithologies transforming Ca-bearing plagioclase to almost pure albite 5 . Albite is also abundant in unconsolidated sediments, either as detrital grains of various compositions or grows as authigenic crystals of pure albite 6 , and has been reported in oil shales associated with organic matter 7 .
Albite also occurs in wide variety of planetary materials including chondrites ~10% and Martian meteorites 8,9 . Sodic-plagioclase has also been identified in the exposed rocks and soils at several the landing sites of the Mars Exploration Rovers 10 and as a widespread component of andesitic lava flows on the surface 11 . In addition, the Thermal Emission Spectrometer on the Global Surveyor spacecraft orbiting Mars have also identified sodic-plagicoclase as dust particles in the atmosphere 12 . There is also speculation that albite may be present on Mercury 13,14 .
To understand the dynamics of the planetary crust it is important to gain constraints on the physical properties of planetary crusts and its constituent rocks and minerals. Physical properties of plagioclase feldspar vary strongly as a function of chemistry between the albite and anorthite end-members. For instance, the elastic parameter such as bulk modulus increases by 60% from albite to anorthite 15 whereas the transport properties such as viscosity vary over three orders of magnitude at constant temperature 16 . It is fundamental to address how the crystal structure and chemistry influences the physical property across the plagioclase solid solution. Albite also forms a solid solution with orthoclase (KAlSi 3 O 8 ) i.e., the alkali feldspar solid solution series [(Na,K)AlSi 3 O 8 ]. Albite is therefore the key model crystal structure for the aluminosilicate feldspar group of mineral. The crystal structure of albite consists of linked corner-sharing TO 4 tetrahedral units (where T = Al, Si) forming a three-dimensional aluminosilicate framework, often referred to as crankshaft structure.
Despite being a major component in continental and oceanic crust and occurring in wide variety of Earth and planetary settings, our knowledge of fundamental physical properties, e.g. elasticity, of albite is rather limited. This is primarily owing to its complex crystal structure and triclinic space group symmetry. In a pioneering experimental study, ultrasonic velocities were measured for plagioclase single-crystals with 9-60 mol % of anorthite. Ultrasonic measurements at 1 bar and 298 K, in adequate directions of the single-crystals allowed for the determination of 13 independent elastic constants i.e., the elastic stiffness tensor with monoclinic rather than the true triclinic symmetry were determined 17 . The reported values were later revised, but still used monoclinic symmetry 18 . More recently, the full elastic tensor for albite 19 , plagioclase feldspars 20 , and alkali feldspars 21 have been reported at ambient pressure and room temperature. In addition, room pressure elastic constants for plagioclase feldspar have also been predicted using ab initio methods 22 . With increasing depth in the crust, i.e., with increasing temperature and pressure, the crystal structure of albite becomes thermodynamically unstable, while denser mineral phases such as jadeite and quartz replace albite 23 . However, the transition to from albite to a mixture of jadeite and quartz also requires elevated temperatures. It is therefore likely that albite may persist as a metastable phase to greater depths 24 . Thus, it is important to know the pressure-dependence of this crystal structure and its physical properties, including elasticity. Several high-pressure studies have therefore been devoted to elucidating the crystal structure and compressibility 24-26 of albite. These studies have indicated a pressure-induced softening of bulk moduli at around 6 GPa. We still do not have any understanding of how the full elastic stiffness tensor and the bulk and shear moduli change upon compression. In order to understand the elasticity and anisotropy of planetary crusts, we have performed first principles calculations of albite with triclinic symmetry at high-pressure. The equation of state and elastic parameters of albite feldspar at high-pressures will be of crucial importance for the thermodynamic database to predict phase relations and velocity discontinuities across crust and mantle [27][28][29][30][31] .

Results
The compression behavior of the albite could be described with a Birch-Murnaghan equation of state 32 . The unit-cell volume and cell axes predicted with LDA and GGA bracket the previous X-ray diffraction results (Fig. 1). LDA typically overbinds, whereas the GGA under-binds, i.e., unit-cell volume data are slightly larger than the experimental results (Fig. 1). The inter-axial angles from our structural optimizations of albite using GGA are in excellent agreement with the previous X-ray diffraction studies 26 (Fig. 1). The unit-cell volume predicted by GGA at 0 GPa shows V GGA 0 > V exp 0 by ~3.4%, whereas unit-cell volume predicted by LDA shows V LDA 0 < V exp 0 by ~4.0% (Table S1, Fig. 1). The bulk moduli predicted using GGA at 0 GPa shows K GGA 0 < K exp 0 by ~1.1% whereas the bulk moduli predicted using LDA shows K LDA 0 > K exp 0 by ~14.5% (Table S1, Fig. 1). Albite (triclinic symmetry) has 21 non-zero independent elastic constants, c ij -three diagonal and compressional elastic stiffness components with, i = j, and i = 1-3; three diagonal and shear elastic stiffness components with i = j, and i = 4-6, and 15 off-diagonal elastic stiffness components with i ≠ j and i = 1-6 33 . The agreement between the predicted elastic constants at room pressures from this study and previous experimental results are  quite good 19 (Fig. 2). For most of the components of full elastic moduli tensor, the experimental result at ambient condition is bracketed by the predictions based on LDA and GGA. The variation of elastic stiffness components with pressure can be explained with a finite-strain formulation 29,34 . Typically, the elastic moduli stiffen (i.e., > 0 dc dP ij ) upon compression, however, in the case of albite, most of the predicted elastic moduli based on LDA and GGA exhibit anomalous behavior, (i.e., < 0 dc dP ij ) between 4-6 and 6-8 GPa respectively (Fig. 2). The predicted anomalous behavior in the elastic parameters is in good agreement with the single-crystal X-ray diffraction study at high-pressure 26

Discussions
The compressibility and elastic stiffness tensor of albite at high-pressure can be explained by the response of the alumino-silicate framework of tetrahedral units upon compression. In the albite crystal structure there are four distinct tetrahedral sites 26,35 , T 1o , T 1m , T 2o and T 2m (Fig. 3). The aluminum atoms occupy the T 1o site based on neutron diffraction 35 and subsequent empirical and DFT studies indicated that the energetics of Al ordered in T 1o sites are of the order of 30 meV lower than T 2o sites 36 . In a recent study, it has been noted that variation in Al/Si order and exact chemistry of alkali site (Na or K) has only minor effect on the aggregate elasticity 21 . However, elasticity varies significantly across Na-Ca plagioclase series 22 Fig. 3; Supplementary Movies-1-4). In feldspar, four distinct tilts are recognized, i.e., φ 1 to φ 4 (Supplementary information). In triclinic albite, owing to lower symmetry, there are two distinct components for φ 1 and φ 2 tilts, i.e., φ 1 is decomposed into φ 1o and φ 1m . Similarly the tilt φ 2 is decomposed into φ 2o and φ 2m . In the present study all the tilt component shows anomalous behavior between 4-6 and 6-8 GPa for LDA and GGA respectively. This is a clear demonstration of the "cause and effect" relationship between the pressure dependent changes in the internal structure and the elasticity (Fig. 4). The flexibility of the alumino-silicate framework of tetrahedral unit results in significant anisotropy in elasticity. The a* direction in albite is significantly softer than the b* and c* direction. The pressure dependence of the P-wave velocity predicted from DFT is in very good agreement with previous measurements to pressures of 2.5 GPa 37,38 (Fig. 5). Recent elasticity measurements predicts stiffer velocities for the b* and c* directions 19 (Fig. 5). Although, the earlier work on feldspar megacrysts 37,38 are "softer" relative to recent work 19 in directions perpendicular to cleavage, i.e., b* and c* directions, the comparison is in better accord with the "true" single crystal behavior in the non-cleavage direction (100). It is likely that in the earlier experiments to determine the elastic constants 38,39 , the closing of the cracks at low-pressure lead to significantly greater pressure derivatives (Fig. 5). Agreement between all the previous experimental studies and present results are excellent for the softer a* direction. Substantial drop in compressional wave velocity, V P occur between 4 and 8 GPa for the a* and c* crystallographic directions. These correspond to the discontinuous behavior of the individual elastic constants, c 11 , and c 33 (Fig. 2). The pressure dependence for shear elastic constants c 44 and c 66 also low values between 4 and 8 GPa, which would correspond to velocity drops for V S for propagation along a* and c* directions. Pressure-induced softening of the shear elastic constants have also been reported for the tetrahedraly coordinated silica (SiO 2 ) polymorphs that shares corners and forms 3-D framework, as in quartz [40][41][42][43][44] and coesite 45,46 .
The full elastic anisotropy compares quite well with recent experimental study 19 and the DFT study 22 (Fig. 5). The P-wave azimuthal anisotropy, AV P decreases from 43% at 1 bar to 34% at 8 GPa to, whereas the S-wave anisotropy AV S increases slightly from 46% at 1 bar to 49% at 8 GPa. At 1 bar, the anisotropy of P-wave to S-wave velocity ratios i.e., V P /V S1 and V P /V S2 are 68 and 59% respectively. Upon compression of albite to 8 GPa, the anisotropy in V P /V S1 and V P /V S2 reduces to 58 and 36% respectively.
The slight discrepancy elasticity and anisotropy between the previous DFT study and the present study is likely related to the differences in the computational parameters used. For instance, a higher energy cut-off of 800 eV is used in this study compared to 582 eV. In addition, the present study uses projector augmented wave method (Method section) as opposed to the pseudopotential method used in earlier DFT study 22 . The Mohorovicic discontinuity (Moho) marks the first major discontinuity separating the earth's crust from the underlying denser mantle. The discontinuity is likely to be related to the changes in composition across crust and underlying mantle 48,49 . Across Moho the rock type is likely to change from gabbro consisting of pyroxene and plagioclase minerals to eclogite consisting of denser omphacite pyroxene and garnet. Several mineralogical transformations play important role in eclogitization and the proportion of mineral phase changes gradually across the transition rendering a continuous change in density and physical properties 49 . In a recent study, it has been suggested that for the thickened continental crust, the Moho could, in part, be explained by the transformation of albite (ab) to a mixture of jadeite (jd) and quartz (qtz) 50 . In this study, we used elasticity results of albite (ab) from this study and recent studies 19,26 and combined them with the recent elasticity results of jadeite (jd) 51,52 and quartz (qtz) 47 to evaluate the discontinuity in compressional velocity across the univariant transformation: NaAlSi 3 O 8 (ab) = NaAlSi 2 O 6 (jd) + SiO 2 (qtz). We used a thermodynamic code, Perple X 28 and the thermodynamic database 27 to predict the phase boundary and the changes in the physical properties such as density and P-, and S-wave velocity across the boundary (Fig. 6). In our analysis, the pressure derivatives of the velocities for albite i.e., dV dP P and dV dP S are from this study. The pressure derivatives of the velocities for jadeite 51,52 and quartz 60 are estimated from previous results. The temperature derivatives of the velocities i.e., dV dT P and dV dT S are derived from Perple X 28 . We find that at a depth of around ~40 km, the average P-and S-wave velocity of mineral assemblage consisting of jadeite and quartz is ~1.0 km/s greater than that of albite (Fig. 6). This is also very similar to the known contrast between crust and the mantle. The weighted average crustal P-wave velocity across several continental crustal settings is ~6.45 (± 0.2) km/s where as the underlying average mantle velocity is ~8.09 (± 0.2) km/s 61 has a similar velocity contrast across the crust-mantle discontinuity.
Alkali oxide-Na 2 O is a minor component in the deep Earth mineralogical model 49 , metasomatism in subduction zones by Na-rich fluids could stabilize jadeite rich rocks-jadeitites and albite rich rocks-albitite which occur together with serpentinite melange [62][63][64] . In the deep Earth, beyond the thermodynamic stability field of albite and jadeite, alkalis including sodium could be incorporated into denser mineral phases such as Na-hollandite  53 . Mineral abbreviationsqtz, coe, and st refers to quartz, coesite, and stishovite with SiO 2 stoichiometry i.e., silica polymorphs; ab and na-holl refers to albite and sodium hollandite with same stoichiometry NaAlSi 3 O 8 ; ne and na-cf refers to nepheline and sodium calcium ferrite structured phase with NaAlSiO 4 stoichiometry and jd refers to jadeite with NaAlSi 2 O 6 stoichiometry. Density and bulk modulus data are derived from-qtz 47 ; coe 54 ; st 55 ; ab 20,26 , this study; na-holl 56 ; ne 57 , na-cf 58,59 . and calcium ferrite (cf) structured phase with NaAlSiO 4 stoichiometry 53 . Upon compression, a series of mineral transformation and densification is likely to modify albite to jadeite + quartz/coesite/stishovite and eventually to Na-holandite 53 . Similarly, jadeite could also be transformed to a mixture of NaAlSiO 4 (cf phase) and stishovite 57 . Both hollandite and cf phases have crystal structure consisting of SiO 6 and AlO 6 units that forms tunnels where alkali atom such as Na and K reside 49,[65][66][67] . The density and elasticity of sodium bearing phases including albite exhibit a positive correlation. Recent experimental studies have indicated that sodium could also be incorporated in ringwoodite 68 . How sodium is partitioned between major mantle minerals i.e., magnesium silicates and minor aluminosilicate at high pressure remains unknown and will be important to constrain the fate of sodium in the deep Earth.

Method
We investigated albite using density functional theory (DFT) [69][70][71] [Hohenberg and Kohn, 1964;Kohn and Sham, 1965]. DFT-based studies of the energetics, equation-of state, and elasticity have been widely used to examine condensed matter including minerals that are stable in the Earth's interior [75][76][77][78][79][80] . We used the local density approximation (LDA) and the semi-local generalized gradient approximation (GGA) [81][82][83] . In addition, we also used the projector augmented wave method (PAW) 84 within the Vienna ab initio simulation package (VASP) [84][85][86][87] . We used an energy cut-off E cut ranging from 400 eV to 1100 eV and k-point sampling ranging from 2 to 46 k-points in the irreducible wedge of the Brillouin zone. We find that the results are fully converged at an energy cut-off E cut of 800 eV and a 3 × 2 × 3 k-point mesh Monkhorst-Pack 88 (Supplementary Figure S2). The total energies and pressures were converged to within 1.06 meV and 0.35 GPa. We used the experimentally determined crystal structure of albite from X-ray and neutron diffraction 35 as the input for the crystal structure and the starting point for a full structural optimization of the crystal structure. We performed all the calculations in the primitive unit cell (Space group: P1, 52 atoms). The static (T = 0 K) calculation leads to primitive symmetry instead of the C1 symmetry observed in experiments. For the determination of the elastic stiffness tensor, we strained the crystal structure and relaxed the internal degrees of freedom while preserving the symmetry of the crystal structure. The elastic constants c ijkl were obtained by relating changes in stress with the applied strain, σ ij = ∑ kl c ijkl ε kl . We applied ± 1% strains to accurately determine the stresses in the limit of small strain (ε kl → 0) 29,34 . The Cartesian frame for the elastic constants has X 2 parallel to crystallographic axis b, X 3 parallel to c*, and X 1 normal to X 2 and X 3 , forming a right-handed system for triclinic symmetry. We used the petrophysical software to determine the elastic anisotropy 89