Reply to the referee’s comments:

The manuscript from Ulian et al. describes a Density Functional Theory (DFT) based computational study on the high-pressure structural evolution and the elastic properties of 1M and 2M₁ polytypes of KMg₃Si₃AlO₁₀(OH)₂ phlogopite. More specifically, the calculations report the static equation of state parameters (K₀, K’ and V₀) and the full elastic tensor for both polytypes, from which polycrystalline average properties (K_voigt, K_reuss, µ_voigt, µ_reuss, E, υ) and seismic anisotropy are derived. The novel aspect of this work relies on the adoption of a posteriori correction of the hybrid B3LYP functional to include the dispersion effects on the elastic properties of both polytypes.

The data obtained from the theoretical calculations are overall of good quality and valuable to the mineralogical community. Although the effort and success in predicting elastic properties (11 independent elastic constants) of monoclinic compounds are remarkable, the general presentation of the manuscript needs to be improved and parts of the text rephrased to carefully explain and discuss the methodology and new findings, to make it easier to understand to non-expert readers.

In addition, the manuscript in this current form lacks a comprehensive framework describing the relevance of phlogopite in geological sciences and in particular why it is important to have a firm knowledge of both polytypes elastic properties. This is now barely described in the introduction - where the authors only state “it is also important to know the elastic properties of this mineral to understand and explain geophysical observations in subduction zones”. Proper citations necessary to provide context on how the results obtained could help understand geophysical observations (not detailed in the text) are also lacking. Some possible technological applications of phlogopite are mentioned in the introduction section, but besides a few citations, no more details are provided. This impacts directly the conclusions of the manuscript, which is only a short sum up of the main results obtained, without adding any sort of implications concerning Earth Sciences nor Material Sciences. Both the introduction and the conclusion of the manuscript therefore need to be considerably implemented by adding further context and implications.

While this manuscript represents the first attempt to calculate the elastic properties of 2M₁ polytype, the full elastic tensor of phlogopite 1M was already provided by Chheda et al. (2014) (Chheda, T. D., Mookherjee, M., Mainprice, D., dos Santos, A. M., Molaison, J. J., Chantel, J., Manthilake, G., and Bassett, W. A.: Structure and elasticity of phlogopite under compression: Geophysical implications, Physics of the Earth and Planetary Interiors, 233, 1-12, 10.1016/j.pepi.2014.05.004, 2014). Importantly, despite adopting a more sophisticated computational method with respect to Chheda et al. (2014), from the current manuscript it is not clear what the improvement and the advancement with respect to the previous study are. In particular, the study by Chheda et al. (2014) reports the high-pressure evolution of the elastic constants and seismic anisotropy, describing the geophysical implications relevant to subducting zone settings, which are not presented in this manuscript, making it hard to judge the novelty of this work.

In conclusion, at the current stage, I believe the manuscript is unsuitable for publication in a broad-impact journal such as Solid Earth.

Answer: We thank the referee for the general positive opinion on our work. As reported detailed in the following, by point-by-point answers, we amended the manuscript accordingly to better elucidate the scope of the present research and to ameliorate the presentation and discussion of the results to meet the criteria of Solid Earth.

Specific comments:

1.Introduction

Line 44: when writing that the knowledge of phlogopite elastic properties could “explain geophysical observations in subduction zones” some context and citations should be provided. Which seismological/geophysical observations are attributed to phlogopite?

Answer: we thank the referee for pointing out this issue. We rewrote the section, including more information and citations related to both the geological/geophysical and technological implications of a more detailed knowledge of the elastic properties of phlogopite. The new paragraph now states:

“Micas, and then phlogopite, could be used as a geothermobarometer to explain petrogenetic processes occurring in high-pressure and high-temperature (non-ambient) conditions (see Guidotti and Sassi, 1976). It is also a possible carrier of water and potassium in the Earth’s mantle (Sudo and Tatsumi, 1990). Moreover, as suggested by Melzer and Wunder (Melzer and Wunder, 2001), since phlogopite may be formed in the mantle wedge that overlies a subducting slab, it could play a relevant role in the element ratios of large ions (K, Rb, Cs) in the subduction zone.

From the technological perspective, phyllosilicates are considered low-cost insulating materials that could be used in two-dimensional optoelectronic applications, e.g., 2D transistors and heterostructures (Ulian and Valdrè, 2025), thanks to their easy exfoliation, which is related to the perfect cleavage on their (001) planes.

Thus, a detailed knowledge of the elastic properties of phlogopite is of utmost importance for both geological studies, e.g., to explain geophysical observations involving phlogopite-bearing rocks (Geng et al., 2024;Van Reenen et al., 2023;Wang et al., 2023), and to devise new technological uses of this mineral for advanced applications (Cadore et al., 2022).”

Line 52: here you should provide more details on why it is important to include long-range interactions in the physical treatment of phlogopite. This would particularly benefit the non-expert readers.

Answer: according to the reviewer’s suggestion, we included the following sentence to better elucidate this important point:

“Albeit the general consensus is that the TOT layers in micas are held together by the interlayer cations (K+), hence the interaction is mostly of Coulombic (electrostatic) nature, previous works on muscovite and phlogopite showed that the contribution from van der Waals interactions is non negligible in determining the crystal structure, elastic and thermodynamic properties of these minerals (Ulian and Valdrè, 2015c;Ulian and Valdrè, 2023a).”

Line 55: even if a complete knowledge of the elastic behaviour of both 1M and 2M₁ polytypes is given, how can this improve the interpretation of geophysical observations in subduction zones? The stability of polytypes cannot be predicted by thermodynamics, since the presence of either 1M or 2M₁ is not a matter of free energy differences, but rather kinetics. So even if some geophysical observations are attributed to phlogopite, it is extremely difficult to point out which polytype is the main responsible or, if both polytypes are present simultaneously, to try to understand the relative amounts.

Answer: we agree with the referee that, in general, the variations of the elasticity of phlogopite attributed to the two polytypes could be very small. As also reported in next comments, the two polytypes show almost identical bulk moduli (within the fitting uncertainty), however, the properties (Young’s modulus, compressibility, shear modulus, …) derived from the elastic tensor showed small differences according to the direction of observation. This behaviour was indeed the one we expected, and the hypothesis behind our study. We included the following sentence at the end of the paragraph pointed out by the referee to better explain the scope of our work:

“The differences in the elastic behaviour of the two phlogopite polytypes could be indeed subtle, especially in terms of the second-order elastic moduli tensor, where small variations in the directional properties between the 1M and 2M₁ are hypothesised.”

2.Computational methods

Line 67-68: the paper of Grimme (2006), which is extremely important as it contains the theoretical details of the treatment of dispersion contributions, is not present in the References section.

Answer: The reference section now includes the fundamental work of Grimme (2006).

Line 68: the parameters that define the dispersion contribution to the total energy are not only functional-dependent, but also compound-dependent, so care should be taken when using the parameters adopted by other authors on different compounds with respect to the one investigated in this study. Since Civalleri et al. (2008) worked on “NH3, acetylene, CO2, urea, urotropine, propane, benzene, naphthalene, formamide, formic acid, 1,4-dichloro-benzene, 1,4-dicianobenzene, succinic anhydride and boric acid”, can the parameters adopted for such compounds be directly employed for phlogopite as well? Civalleri et al. (2008) used the parameters reported in Table 1 of Grimme (2006), therefore one straightforward way to remove any doubt about the parameters employed for the calculations on phlogopite would be to add a table in the supplementary material with the C₆ parameter of K, Mg, Si, Al, O and H, the scaling factor s₆ adopted for the B3LYP functional, and the van der Waals radii of the various atoms. Moreover, the B3LYP-D* correction employed by Civalleri et al. (2008) consists of an empirical rescaling of the scaling factor of the B3LYP functional and of “the van der Waals radii of heavy atoms and hydrogen”. As reported by Civalleri et al. (2008): “Proposed scaling factors were determined from a manual procedure by progressively increasing the atomic van der Waals radii and trying to find the best agreement between computed and experimental data”. This means that the parametrization they performed is compound-specific and calibrated over experimental results on mostly organic compounds that do not contain K, Mg, Si or Al. Again, I suggest you provide more details on the parameters that were employed for your correction either in the supplementary material or directly in the computational methods section.

Answer: According to the referee’s comment, we included the following sentence immediately after introducing the B3LYP-D* approach:

“According to the proposed methodology, the original B3LYP-D2 parameters were modified as s₆ = 1, the van der Waals radius of the hydrogen atom, R_vdw(H), was set to 1.30, and the R_vdW of the heavier atoms were scaled by a factor of 1.05, while the C₆ values were the same proposed in the DFT-D2 scheme (see Table S1 in the Supplementary Materials).”,

and an additional table (Table S1) was included in the Supplementary Materials.

Line 72: For O basis set, I think it is better to cite this work: “L. Valenzano, F.J. Torres, K. Doll, F. Pascale, C.M. Zicovich-Wilson, R. Dovesi, "Ab Initio study of the vibrational spectrum and related properties of crystalline compounds; the case of CaCO3 calcite", Z. Phys. Chem. 220, 893-912 (2006). DOI 10.1524/zpch.2006.220.7.893”, as it is the first study in which the different contraction exponents and coefficients for O were tested and that led to the 8-411d11G set.

Answer: According to the referee’s suggestion, we replaced the reference of Valenzano et al. (2007) with the one from 2006.

3.Results and discussion

Lines 87-99: this section should go in the methods paragraph.

Answer: According to the referee’s suggestion, we moved this paragraph in the Computational Materials section.

Lines 120-121: it is necessary to specify how many, and what volume states were considered, as this not only can affect the E-V fitting but also prevents the reproducibility of the calculations.

Answer: in the revised manuscript, we modified the sentence to include the number of unit-cell volumes and the relative variation concerning the equilibrium cell. The text now states:

“The compressional behaviour of phlogopite was modelled considering 10 different unit-cell volumes, both smaller (compressed, 7 volumes) and larger (expanded, 3 volumes), within 88% and 108% of the equilibrium geometry volume (V_eq).”

Line 130 equation 1: the term (Ƞ²-1)³ should be a factor, not an exponent. Also, it would be better to specify the natural variable of energy: E(V).

Answer: we amended the term in Eq.(1) and we included the natural variable of energy, as per the suggestion. Now Eq.(1) is:

Line 132: the term “dilaton” can be avoided as it’s rarely used and can be confusing.

Answer: According to the suggestion, we removed the term from the text, stating that it is only a dimensionless parameter.

Lines 133-134: “The pressure values at each unit cell volume (reported in Table 2)…” Table 2 does not report any pressure value, nor unit cell volumes other than the equilibrium one. Are you referring to Table S2?

Answer: We thank you for pointing this out. Indeed, we referred to Table S1 for the 1M polytype and Table S2 for the 2M₁ one. We amended the sentence as:

“The pressure values at each unit cell volume (reported in Tables S1 and S2 for the 1M and 2M1 polytypes, respectively) were calculated ...”

Line 136 equation 3: As already commented for equation 1, the natural variable of pressure should be specified: P(V).

Answer: we amended equation (3) including the natural variable of pressure, as also done for Eq.(1). Now the formula is written as follows:

Line 137: “Table 3” could be a typo, probably you were referring to Table 2 instead?

Answer: Indeed, it was a typo. We amended the reference to the table, now correctly pointing to Table 2.

Lines 139-140, with reference to Figure 2: “In general, there is a fine agreement between the relative variation of the cited structural properties obtained from the present theoretical simulations and the experimental ones from X-ray diffraction”. Please be more specific and define the differences between experiments and calculations, and between these calculations and the previous calculations from Chheda et al. (2014). You need to explain the deviation of computed volumes with respect to experiments at high pressures in Figure 2 (panels A, B and D)? Also, add Chheda et al. (2014) results to Figure 2 and provide the fit of P-V data.

Answer: According to the referee’s suggestion, we added the results of Chheda and collaborators (2014) in Figure 2, alongside the experimental ones of Comodi et al. (2004) and Pavese et al. (2003). In addition, we included a supplementary figure (Figure S1) that reports the absolute values of the unit cell volume and lattice parameters as a function of pressure for the phlogopite-1M. This new figure better highlights the differences between our theoretical approach and that proposed by Chheda and co-workers. Finally, the section related to the experimental/theoretical comparison was expanded, according to the referee’s comment.

Line 141: (comment on Table 2) I think that a valuable comparison that should have been provided to support the importance of dispersion correction in evaluating EoS parameters would have been a calculation of static EoS at B3LYP level with no dispersion correction, to see how different the K₀ (DFT) and K₀ (DFT/D*) actually are and how impactful dispersion effects are on compressibility. This comparison has been done in a previous paper by the authors (Ulian et al., 2021).

Answer: According to the referee’s suggestion, we included in Table 2 and Table 3 the EoS fit results obtained for both polytypes, using the DFT/B3LYP approach without the modified DFT-D2 correction. It can be noted the significant drop in the bulk and axial moduli due to the absence of long-range interactions in the theoretical framework, which is in line with our previous results related to other phases (Ulian et al., 2021).

Line 141: (comment on Table 2) The theoretical bulk moduli of 1M and 2M₁ polytypes provided by this work are almost identical, so how different is their compressional behaviour?

Answer: we agree with the referee, the two moduli are almost the same between the 1M and 2M₁ polytypes. Hence, on a general basis, the two polytypes behave mostly the same under hydrostatic compression. However, the interesting result resides on the axial moduli, related to the a ,b and c behaviours upon compression. Indeed, while the M₀(a) and M₀(b) values are the same between the 1M and 2M₁ phlogopite models, the M₀(c) modulus is significantly higher in the latter polytype, with an increase stiffness of about 10%. We included the following sentences to better explain the different elastic behaviour:

“In general, it can be noted that the K₀ values are almost the same between the two polytypes. However, the 1M and 2M₁ phlogopite models showed a slightly different behaviour in terms of axial compression. In fact, while the M₀(a) and M₀(b) values are the same between the two phlogopite models, the M₀(c) modulus is significantly higher in the 2M₁ polytype, with a stiffness increase of about 10% with respect to the 1M one.”

Line 151: Please rephrase “Similar figures were calculated..”

Answer: according to the referee’s suggestion, we rephrased the sentence as:

“Similar results were obtained by Chheda et al. (2014) from DFT simulations and by Pavese et al. (2003) from XRD experiments. Also, the recalculated values of Comodi et al. (2004) with the linearized formulation lead to a bulk modulus within 0.57% the value obtained from the PV EoS fitting.”

Line 193: I believe it is not necessary to specify here which keyword allows to perform calculation of the elastic tensor in CRYSTAL.

Answer: according to the suggestion, we removed from the text the reference to the ELASTCON keyword used in CRYSTAL to perform the calculation of the elastic tensor.

Line 198: “whose values were reported”. Maybe “are reported” is better?

Answer: we changed the verb to “are reported” as suggested by the reviewer.

Line 200: “the present simulations were in good agreement with the experimental results” by looking at Table 4, the presented results are rather different from the references provided… Again, please discuss the differences. Also, it would be useful to have a figure displaying the C_ij evolution with pressure. See for example Chheda et al. (2014)

Answer: We agree with the referee that the theoretical C_αβvalues are slightly different from the experimental ones, and we extended the discussion to better show the observed variations. However, the evolution of the elastic tensor components with pressure was not considered in the present work, which was more focused on the crystal-structure properties of phlogopite under pressure. This would be the topic of a subsequent work.

Lines 210-211: in the text the authors mention that there is a systematic overestimation due to 1) absence of thermal effects in the calculations and 2) presence of Pulay forces…. I have some concerns about this: thermal effects may explain the overestimation with respect to experimental values, however the C_ij reported in this study differ also from those of Chheda et al. (2014). Even if the overestimation with respect to plane waves calculations is to be attributed to the use of GTOs, PAW results are still in slightly better agreement with experiments. A routine is currently implemented in CRYSTAL that allows to remove eventual BSSE via a geometrical counterpoise method. Could that improve your results and mitigate the effect of the basis set? As reported in a previous comment in the equation of state section, it would have been nice to see a comparison between a B3LYP-D* corrected and B3LYP non-corrected simulation. Also, if an overestimation is “systematic” an estimate of such overestimation should be reported to quantify the expected mismatch.

Answer: We thank the referee for pointing this out, as this is one of the key concepts that must be clear to readers. Let’s focus on the theoretical aspects. It is possible to see from the results reported in Table 4 that the calculation of stresses/forces is rather sensitive to the employed approach. LDA is known to be an overbinding functional, thus, the elastic tensor components tend to be overestimated. Conversely, GGA functionals are underbinding and generally result in the opposite behaviour as shown by Chheda and collaborators (2014). Then, LDA and GGA represent the extremes of a range in which the elastic moduli should be found in athermal conditions. Hybrid functionals like B3LYP, being two rungs above GGAs in the Jacob’s ladder conceived by Perdew, should provide elastic moduli that, ideally, fall between the LDA and GGA results. Hence, the results obtained from hybrid functionals should agree better with the experiments.

The above considerations will hold if we exactly use the same computational parameters in the simulations with different functionals: basis set quality, k-point sampling, convergence criteria, and so on. Thus, a different basis set type, such as localised Gaussian-type orbitals (our work) and plane waves (Chheda et al., 2014), is expected to lead to small or large variations. Compared to PW, GTOs are affected by BSSE, which in turn provide some bias in the calculation of forces/stresses due to Pulay stress, as explained with more details in the cited references (Ulian and Valdrè, 2018; Ulian et al., 2021). Generally, the more complete the basis set is, the lower the BSSE. However, considering the size of the phlogopite models, especially the 2M₁ polytype, very large basis sets like the triple-ζ one proposed by Peintinger, Oliveira and Bredow (POB-TZVP) are computationally very expensive. This is why we chose a double-ζ basis set that was already tested in previous works on phyllosilicates, which provided results with adequate accuracy compared to experimental findings. Unfortunately, we could not employ the geometrical counterpoise (gCP) approach proposed by the research group of Prof. Grimme because of it is limited to a selection of general basis sets, such as the cited POB-TZVP. Moreover, gCP has been tested only for molecular crystals such as urea, and only recently we begun analysing its suitability for inorganic solid systems, such as realgar (see https://doi.org/10.1107/S1600576724000025). We are currently continuing these tests with other systems, however, a thorough analysis of the accuracy of the gCP method for minerals as complex as phlogopite was not within the purpose of the present work.

As suggested by the referee, we included the elastic tensor of both phlogopite polytypes calculated with the uncorrected B3LYP functional. It is possible to note that they are more in line with the experimental findings, even without the gCP correction. However, we must remember that our results are referred to athermal conditions, i.e., no zero-point energy and thermal contributions were included in the simulations. If we consider both BSSE and thermal effects, we expect that the elastic moduli will significantly drop below the experimental ones. From this perspective, albeit being larger, we suggest that the B3LYP-D* results are more in line with the elastic wave propagation measurements.

Regarding the term “systematic” we used when discussing the discrepancy between the previous theoretical results of Chheda et al. (2014) and the experimental measurements, the quantification of the overestimation is not trivial, as it depends on both BSSE and thermal effects. A crude estimation could be given from the absolute differences between the theoretical/experimental results, which are now reported in the text.

All this discussion was reported in the revised manuscript.

Line 245: where are equations (14) from?

Answer: The equations are from the work of Ranganathan and Ostoja-Starzewski (2008). This reference was included in the amended text before the mentioned equations.

Lines 277 – 280: A table with a comparison between the numerical values of V_P and V_S predicted in this work and those obtained by Chheda et al. (2014) (and maybe also Alexandrov et al., 1974) would make it easier and more straightforward to compare the presented results and the literature data.

Answer: following the referee’s suggestion, we included a table (Table 6) in the revised manuscript to ease the comparison between the present results and those reported by Chheda et al. (2014) and Alexandrov et al. (1974).

Figures/Tables captions

Table 1: report error bars in the experimental data

Answer: according to the reviewer’s suggestion, we included the errors associated with the different values, when available from the cited references.

Figure 2: “Evolution of (a) unit cell volume V/V₀…” should be “Evolution of normalized (a) volume (V/V₀), and lattice cell parameters (b) (a/a₀)…”

Answer: as suggested by the referee, we included the term “normalized” in the caption of Figure 2.

Table 4: there is an inconsistency between how the components of the elastic tensor are labelled: C_αβ in equations 4 and 5 whereas C_ij in Table 4. Be consistent with the terminology.

Answer: we thank the referee for pointing this out. We amended this inconsistency by labelling the elastic tensor components using the “C_αβ” notation throughout the manuscript.

Supplementary materials

Table S1: if these are the volume states used for the E – V fitting, why is the equilibrium volume not included in the table? Usually when you perform EoS calculations with CRYSTAL the equilibrium volume is always included by default regardless of how many and which E – V points you consider. Same thing for Table S2 on the 2M₁ polytype. In this second case, I guess that the data at P = 0 GPa are those reported in Table 1 of the manuscript, but for the sake of clarity and completeness I would leave them also in the tables provided in the supplementary section, and report in the computational methods section of the manuscript at least the P – V conditions you considered for your static EoS.

Answer: According to the referee’s suggestion, we included the crystallographic properties of both 1M (Table S2 in the revised manuscript) and 2M₁ polytypes (Table S3) obtained in equilibrium conditions (0 GPa).

REVISED MANUSCRIPT

High-pressure behaviour and elastic constants of 1M and 2M1 polytypes of phlogopite KMg₃Si₃AlO₁₀(OH)₂

Gianfranco Ulian¹, Francesca Ranellucci¹, Giovanni Valdrè¹

¹Dipartimento di Scienze Biologiche, Geologiche e Ambientali, Università di Bologna “Alma Mater Studiorum”, Plesso di Mineralogia, Piazza di Porta San Donato 1, 40126 Bologna, Italy

Correspondence to: Giovanni Valdrè (giovanni.valdre@unibo.it)

Abstract. In the present work, the elastic properties of both 1M and 2M1 phlogopite polytypes, KMg₃Si₃AlO₁₀(OH)₂ (monoclinic crystal system) were investigated from PV equation of state fitting and by analysis of the fourth-rank elastic tensor. The analysis was performed within the Density Functional Theory framework, using all-electron Gaussian-type orbitals basis sets and the B3LYP functional corrected a posteriori to include long-range interactions (B3LYP-D*). In general, the elastic properties of the two polytypes were strongly anisotropic, with the axial moduli ratio M(a) : M(b) : M(c) being close to 4 :  4 : 1. The volume-integrated third-order Birch-Murnaghan equation of state fitting parameters at 0 K were K₀ = 57.9(2) GPa, K’ = 8.29(7) and V₀ = 489.82(3) ų for phlogopite-1M, which were very close to those of the 2M₁ polytype, i.e., K₀ = 58.3(1) GPa, K’ = 8.71(8) and V₀ = 978.96(9) ų. The monoclinic elastic tensors obtained for the two polytypes of phlogopite, which have never been experimentally reported for both minerals so far, were in line with the PV behaviour of the mineral, providing further data related to the directional dependence of the elastic properties and seismic wave propagation. The elastic properties from both PV hydrostatic compression and from the elastic moduli tensor were discussed against the available experimental and theoretical data in the scientific literature, extending the knowledge on this important trioctahedral phyllosilicate.

1 Introduction

Phlogopite [commonly abbreviated with Phl, K(Mg,Fe)₃Si₃AlO₁₀(OH)₂, with Mg/Fe ratio greater than 2], is a ubiquitous mica mineral that can be found in different geological environments and pressure-temperature conditions, i.e. igneous, sedimentary, and metamorphic rocks (Icenhower and London, 1995;Trønnes, 2002), and one of the water and potassium reservoirs in the Earth’s mantle (Virgo and Popp, 2000;Tutti et al., 2000).

From the crystallographic point of view, like other phyllosilicates, phlogopite is made of a 2:1 layer formed by two tetrahedral sheets (labelled as T) that sandwich an octahedral one (O). These T-O-T (or TOT) layers present strong (covalent/ionic) in-plane bonds and a negative charge because of Al^III/Si^IV substitutions in the T sheets, which is balanced by potassium cations in the interlayer. Thus, the structure is held together along the c crystallographic axis by a mix of Coulombic (electrostatic) and weak (van der Waals) interactions (Ventruti et al., 2009). Phlogopite presents different polytypes due to different arrangement of the TOT layers along the [001] direction, with the 1M and the 2M₁ being the most common in nature both belonging to the monoclinic crystal system (Lacalamita et al., 2012), while a third polytype, 3T (trigonal), is known to be very rare (Gatta et al., 2011). A graphical representation of the most common polytypes is shown in Fig.1.

Phlogopite is also a widely adopted mineral in the technological field, for instance in the manufacture of ceramics and glasses (King et al., 2000;Ariane et al., 2023) and, more recently, gained attention as a potential dielectric (insulating) material for optical and electronic applications because it presents a band gap of ca. 5–6 eV that can be modulated by microchemical variations, e.g., by Fe^II/Mg^II substitutions (Ulian and Valdrè, 2023b).

Figure 1. (a) Structural model of the phlogopite-2M₁ model viewed along the [100] direction. In the bottom right corner, a graphical representation of the tetrahedral rotation angle α_r is shown. (b) Schematic view of the stacking of the TOT layers in the 1M (left) and 2M₁ (right) polytypes. The black arrows link two points related by symmetry.

Micas, and then phlogopite, could be used as a geothermobarometer to explain petrogenetic processes occurring in high-pressure and high-temperature (non-ambient) conditions ( see Guidotti and Sassi, 1976). It is also a possible carrier of water and potassium in the Earth’s mantle (Sudo and Tatsumi, 1990). As suggested by Melzer and Wunder (2001), since phlogopite may be formed in the mantle wedge that overlies a subducting slab, it could play a relevant role in the element ratios of large ions (K, Rb, Cs) in the subduction zone.

Also, the bulk chemical composition of altered rocks such as mica-amphibole-rutile-ilmenite-diopside (MARID) is very similar to the chemistry of phlogopite in an olivine – quartz – alkali ternary [(Mg,Fe)₂SiO₄–SiO₂(quartz)–K₂O+Na₂O] (alkali) ternary (Sweeney et al., 1993). From the technological perspective, phyllosilicates are considered low-cost insulating materials that could be used in two-dimensional optoelectronic applications, e.g., 2D transistors and heterostructures (Ulian and Valdrè, 2025), thanks to their easy exfoliation, which is related to the perfect cleavage on their (001) planes.

Thus, detailed knowledge of the elasticity of phlogopite is of utmost importance for both geological studies, e.g., to explain geophysical observations involving phlogopite-bearing rocks (Geng et al., 2024;Van Reenen et al., 2023;Wang et al., 2023), and to devise new technological uses of this mineral for advanced electronic and optical applications (Cadore et al., 2022). Many experimental studies were focused on the equation of state of phlogopite (Pavese et al., 2003;Comodi et al., 2004;Gatta et al., 2011), but mainly on the 1M polytype. As also observed by Chheda et al. (2014), no measurement of the full monoclinic elastic constant tensor has ever been reported in the scientific literature. The only experimental evidence of this elastic property is related to elastic wave propagation (EWP) experiments that provided just the 5 independent constants, considering the mineral as pseudo-hexagonal (Aleksandrov et al., 1974;Alexandrov and Ryzhova, 1961). The full elastic tensor of phlogopite-1M was only recently provided from ab initio simulations within Density Functional Theory (DFT) at the local density approximation (LDA) and generalized gradient approximation (GGA) levels (Chheda et al., 2014). However, the authors did not employ any correction to include long-range interactions in the physical treatment of phlogopite. Albeit the general consensus is that the TOT layers in micas are held together by the interlayer cations (K⁺), hence the interaction is mostly of Coulombic (electrostatic) nature, previous works on muscovite and phlogopite showed that the contribution from van der Waals interactions is non negligible in determining the crystal structure, elastic and thermodynamic properties of these minerals (Ulian and Valdrè, 2015a;Ulian and Valdrè, 2023b). Furthermore, no studies on the elastic properties of phlogopite-2M₁ were ever carried out and reported. The differences in the elastic behaviour of the two phlogopite polytypes could be indeed subtle, especially in terms of the second-order elastic moduli tensor, where small variations in the directional properties between the 1M and 2M₁ are hypothesised.

For all these reasons, the present work is focused on the elastic properties of phlogopite, considering both 1M and 2M₁ polytypes. To this aim, ab initio DFT simulations were performed (i) to understand the compressional behaviour of the mineral through an equation of state fitting, (ii) to obtain the full monoclinic elastic moduli and (iii) to assess the role of polytypism on the elastic behaviour of phlogopite. All the data were compared and discussed with both experimental and theoretical results found in the scientific literature to extend the knowledge on this important phyllosilicate.

2 Computational Methods

As previously introduced, the simulations reported in the present work were carried out within the Density Functional Theory framework. In detail, the hybrid B3LYP functional (Becke, 1993;Lee et al., 1988), which includes 20% of exact Hartree-Fock exchange and some non-local contribution to the exchange-correlation terms, was employed throughout the present study. B3LYP is considered a suitable choice when dealing with the simulation of the structural and elastic properties of minerals due to its high accuracy, and it was already adopted with success for the description of several other minerals (Prencipe et al., 2009;Ottonello et al., 2010;Belmonte, 2017;Ulian and Valdrè, 2023a;Ulian and Valdrè, 2019) and phlogopite (Ulian and Valdrè, 2023b). However, since most generalized-gradient approximation functionals, including their hybrid counterparts, do not adequately include van der Waals (long-range) interactions, the DFT-D2 correction (Grimme, 2006) was adopted, employing the parameters of Civalleri and co-workers (2008) for the B3LYP functional (B3LYP-D* method). According to the proposed methodology, the original B3LYP-D2 parameters were modified as s₆ = 1, the van der Waals radius of the hydrogen atom, R_vdw(H), was set to 1.30, and the R_vdW of the heavier atoms were scaled by a factor of 1.05, while the C₆ values were the same proposed in the DFT-D2 scheme (see Table S1 in the Supplementary Materials). The all-electron Gaussian-type orbitals (GTO) used to construct the Kohn-Sham orbitals within the linear combination of atomic orbitals (LCAO) approach were the same used in a recent work on the 1M polytype of phlogopite (Ulian and Valdrè, 2023b). Si, Al, Mg, K, O and H atoms were described with 88-31G*  (Nada et al., 1996), 85-11G* (Catti et al., 1994), 8-511d1G (Valenzano et al., 2007), 86-511G (Dovesi et al., 1991),  8-411d11G (Valenzano et al., 2006), and 3-1p1G (Gatti et al., 1994), respectively.

The accuracy of the Coulomb and exchange series was controlled by five thresholds set to 10^–8 (ITOL1 to ITOL4) and 10^–16 (ITOL5) for the structural relaxation procedures, whereas the convergence on the total energy was set to 10^–8 Hartree. The reciprocal space was sampled with a 5×5×5 Monkhorst–Pack mesh (Monkhorst and Pack, 1976), corresponding to 39 independent k points.

The lattice and internal geometry were optimised within the same run using a numerical gradient for the unit cell parameters and an analytical gradient for the atomic positions. The Hessian matrix was upgraded with the BFGS algorithm (Broyden, 1970b, a;Fletcher, 1970;Goldfarb, 1970;Shanno, 1970). The tolerances for the maximum allowed gradient and the maximum atomic displacement were set to 10^–5 Ha bohr^–1 and 4 × 10^–5 bohr, respectively.

All the simulations were performed with the CRYSTAL17 code (Dovesi et al., 2018), whereas the QUANTAS software (Ulian and Valdrè, 2022, 2024a) was employed to post-process the elastic data. Graphical representations have been carried out with the molecular graphics program VESTA (Momma and Izumi, 2008).

The initial model of the crystal structure of 2M₁- phlogopite was created from the single-crystal X-ray diffraction data of Lacalamita et al. (2012), which was subsequently optimised within the B3LYP-D* level of theory. The mineral belongs to the monoclinic system (space group C2/c), presenting two independent tetrahedral sites (labelled as T) and two octahedral sites (M), with partial occupancies of the cationic sites with several possible atoms. For example, the TO₄ tetrahedral sites in the unit cell are generally occupied by about 70% by Si and 30% by Al (Laurora et al., 2007;Lacalamita et al., 2012), which can be translated in 1/4 sites being randomly occupied by aluminium. However, as also explained in a previous study on the 1M polytype (Ulian and Valdrè, 2023b), theoretical modelling requires each site to be deterministically occupied by a single kind of atom; in the example cited above, a T site could present either Al or Si but not both at the same time. Thus, to preserve the monoclinic lattice of phlogopite when the Al^III/Si^IV substitutions are included in the structure, the 2M₁ polytype with formula KMg₃(AlSi₃)O₁₀(OH)₂ was described with the P2/c space group (23 inequivalent atoms, 88 atoms in the unit cell, Z = 4 unit formulas), which is a sub-group of the C2/c symmetry. Within the P2/c space group, phlogopite-2M₁ presented four inequivalent T sites (Si^IV in T1–T3, Al^III in T4) and three non-equivalent O sites (M1–M3). Hence, each T site was related to four SiO₄ or AlO₄ tetrahedra, whereas all octahedral sites are related to two MgO₆ octahedra.

3. Results and discussion 3.1 2M₁-Phlogopite crystal structure

Table 1 reports the lattice parameters, tetrahedral (T_thick) and octahedral (M_thick) sheet thicknesses, interlayer distance (I_thick), and selected polyhedral properties of the 2M₁ phlogopite polytype, obtained from DFT simulations with (B3LYP-D*) and without (B3LYP) the correction to include the effects of long-range interactions. Previous XRD refinements are also shown in Table 1 for a direct comparison.

Table 1. Equilibrium geometry of 2M₁-phlogopite as obtained from DFT/B3LYP and B3LYP-D*, compared to previous experimental refinements.

Notes: a – X-ray diffraction data of Laurora et al. (2007); b – single-crystal XRD refinement of Lacalamita et al. (2012) performed on sample BU1_14. T_thick, M_thick, and I_thick are tetrahedral sheet thickness (calculated from the z coordinates of the basal and apical oxygen atoms), the octahedral sheet thickness (calculated from the z coordinates of the apical and hydroxyl O atoms) and interlayer thickness (calculated from the z coordinates of the basal oxygen atoms), respectively. TQE and OQE are the tetrahedral and octahedral quadratic elongations, respectively (Robinson et al., 1971). T and M represent the generic tetrahedral and octahedral sites, respectively.

The simulations at the DFT/B3LYP-D* level of theory are in good agreement with the experimental crystallographic data reported in literature, showing a small underestimation of the unit cell volume (∆V = –1.0%), which is the result of a slight expansion of the tetrahedral and octahedral sheet thicknesses (∆T_thick = 1.5% and ∆O_thick = 1.8%, respectively) and a contraction of the interlayer distance (∆I_thick = –2.5%). These results for the 2M₁-phlogopite are in agreement with previous simulations of the 1M polytype (Ulian and Valdrè, 2023b) and other phyllosilicates and layered minerals (Ulian et al., 2013;Ulian and Valdrè, 2015c, b;Ulian and Valdrè, 2019), where the inclusion of long-range interactions in the physical treatment is of utmost importance to properly simulate the crystal-chemistry and properties of this kind of structure. When the DFT-D2 correction is not implemented in the simulation framework, a dramatic increase in the unit cell volume can be noted (+6 % with respect to the DFT/B3LYP-D* results, +5 % from the XRD data), which is the result of the increase of the lattice parameters a (+0.9 %), b (+0.9 %), and c (+3.6 %) and the shrinking of the β angle (–2.1 %). Considering the absence of any thermal effects in the Density Functional Theory simulations, the overestimation of the lattice constants and unit cell volume is a further sign of the relevant role played by van der Waals interactions in micas.

3.2 Compressional behaviour of phlogopite 1M and 2M₁ polytypes

The compressional behaviour of phlogopite was modelled considering 10 different unit-cell volumes, both smaller (compressed, 7 volumes) and larger (expanded, 3 volumes), within 88% and 108% of the equilibrium geometry volume (V_eq). Then, the internal coordinates and lattice parameters were optimized keeping the selected volume fixed during the procedure, constraining the space group of 1M-phlogopite (P2) and 2M₁-phlogopite (P2/c), a procedure known as “symmetry preserving, variable cell-shape structure relaxation” that was proposed by Pfrommer and collaborators (1997). The equilibrium geometry of 1M-phlogopite was retrieved from the work of Ulian and Valdrè (2023b) that was obtained with the same computational settings here adopted for the 2M₁ polytype.

The dependence of the total energy E of the mineral at the different volumes V, i.e., the E(V) curve, was described in terms of a volume-integrated 3^rd-order Birch-Murnaghan equation of state, BM3 (Birch, 1947), as reported by Hebbache and Zemzemi (2004):

where η is a dimensionless parameter, K₀ is the bulk modulus at 0 K, K’ its pressure first derivative and V₀ the volume at zero pressure. The pressure values at each unit cell volume (reported in Tables S2 and S3 for the 1M and 2M1 polytypes, respectively) were calculated using the fitting parameters K₀, K’ and V₀ obtained from the previous equation and employed in the well-known P-V formulation of the BM3 (Birch, 1947):

with P₀ the reference pressure (0 GPa). In Table 2, the BM3 equation of state parameters for both 1M and 2M₁ phlogopite are reported. A graphical representation of the normalised volume and axis lengths as a function of calculated pressure is presented in Fig.2., alongside the experimental results from X-ray diffraction (Comodi et al., 2004;Pavese et al., 2003) and the theoretical calculations of Chheda et al. (2014) carried out at the LDA and GGA level of theory.

Figure 2. Evolution of normalised (a) unit cell volume V/V₀, (b) a/a₀, (c), b/b₀ and (d) c/c₀ as a function of pressure for the phlogopite polytypes 1M and 2M₁. The experimental data of Pavese et al. (2003) and Comodi et al. (2004) are reported for a comparison.

A direct comparison between the different experimental and computational approaches is however possible on the phlogopite-1M, which is better depicted in Fig.S1 (Supplementary Materials) considering the absolute values of the cell volume and lattice vectors. In general, compared to the high-pressure XRD refinements, the crystal structure from the DFT/B3LYP-D* approach well agrees with the results of Pavese et al. (2003) at low pressure (P < 5 GPa), however there is a slight overestimation of the unit cell volume at high pressure (see Fig.2a and Fig.S1a). Analysing the lattice vectors, it can be noted that the a and b lengths are very close to the experimental ones at each pressure state, whereas the c-axis shows a stiffer behaviour and an overestimation that increases with pressure. By considering the results of Comodi et al. (2004), the present B3LYP-D* data are slightly underestimated in terms of unit cell volume, which is due to the lower values of the a- and b-axes, whereas the evolution of the c lattice parameters is correctly described by the proposed computational approach (see Fig.S2b-d). It is interesting to note that if the long-range interactions are not included in the theoretical framework, the B3LYP functional leads to a general overestimation of the unit cell volume and lattice parameters. The a and b vectors are closer to the results of Comodi and collaborators (2004) than those of Pavese et al. (2003), however, the c lattice parameter is severely overestimated when compared to both XRD results. All these considerations are reflected in the PV fitting using the BM3 equation of state, where the B3LYP-D* approach resulted in a small overestimation of the bulk modulus K₀ when compared to the experimental data (see Table 2), and the theoretical K’ and V0 values between the ranges obtained from the XRD refinements. Conversely, the B3LYP data shows a large underestimation of the bulk modulus of about 18 GPa, and an overestimation of the cell volume.

Table 2. Equation of state parameters and associated standard deviations (σ) obtained for the 1M and 2M₁ polytypes of phlogopite.

Notes: GTO – Gaussian-type orbitals; PAW – projector-augmented wave; SCXRD – single-crystal X-ray diffraction; SXRPD – synchrotron X-ray powder diffraction.

Compared to previous theoretical results, the crystal structure of the phlogopite 1M polytype and the related PV fitting obtained from the uncorrected B3LYP functional are almost identical to the GGA values of Chheda et al. (2014). Instead, the B3LYP-D* approach, which includes the modified DFT-D2 correction, provided crystallographic results that are between those calculated using the LDA (lower bound in Fig.S1) and the GGA functionals (upper bound). Interestingly, the LDA approach led to a correct description of the c lattice parameter with pressure, whereas the generalised-gradient approximation functional resulted in the same overestimation noted with the B3LYP approach. Thus, GGA and hybrid B3LYP functionals provide an adequate description, i.e., within the experimental information, of the a and b lattice parameters, where strong ionic/covalent bonds are responsible of the TOT structure; though, both approaches fail in properly describing the c-axis behaviour. This is a further confirmation of the hypothesis that the forces acting along the [001] direction of micas are not only electrostatic/ionic, but there is also a contribution from weak van der Waals interactions (Ulian and Valdrè, 2015c, 2023a), as demonstrated using the DFT-D2 correction. The authors are aware that the B3LYP-D* approach, whose DFT-D2 parameters were empirically recalibrated on molecular crystals (Civalleri et al., 2008), is just one of the possible methods to properly treat long-range interactions, and other more ab initio methods could be employed, e.g., the vdW functionals developed of Dion and collaborators (2004). However, the proposed theoretical framework was already suitable to provide new insights into this fundamental topic regarding the bonding nature of mica phyllosilicates.

The axial compressibility, M₀, was obtained using a linearized form of the third-order Birch-Murnaghan equation of state, as described by Angel et al. (2014), which was also employed to re-calculate the experimental compressibility reported by Pavese et al. (2003) and Comodi et al. (2004) including the uncertainties on pressure and lattice parameters. The results are shown in Table 3, alongside the theoretical results of Chheda and collaborators (2014). The results are consistent for all the simulations, as the bulk modulus  calculated from the axial compressibility is within 2.64% (1M) and 1.67% (2M₁) from the K₀ value obtained from the equation of state fit. In general, it can be noted that the K₀ values are almost the same between the two polytypes. However, the 1M and 2M₁ phlogopite models showed a slightly different behaviour in terms of axial compression. In fact, while the M₀(a) and M₀(b) values are the same between the two phlogopite models, the M₀(c) modulus is significantly higher in the 2M₁ polytype, with a stiffness increase of about 10% with respect to the 1M one. Similar results were obtained by Chheda et al. (2014) from DFT simulations and by Pavese et al. (2003) from XRD experiments. Also, the recalculated values of Comodi et al. (2004) with the linearized formulation lead to a bulk modulus within 0.57% the value obtained from the PV EoS fitting.

Table 3.  Axial compressibility parameters obtained from linearized equation of state fit.

Notes: a – present work; b – Pavese et al. (2003); c – Comodi et al. (2004); d – Chheda et al. (2014). GTO – Gaussian-type orbitals; PAW – projector-augmented wave; SCXRD – single-crystal X-ray diffraction; SXRPD – synchrotron X-ray powder diffraction.

Crystal structure data of both 1M and 2M₁ polytypes are reported in Tables S2 and S3 (Supplementary Materials), respectively, providing further insights into the pressure effects on the phlogopite internal geometry.

The TO₄ tetrahedra showed a stiff behaviour, with zero-pressure bulk moduli equal to 282(3) GPa and 285(2) GPa for the 1M and 2M₁ polytypes, respectively, whereas the octahedral MO₆ units were more compressible, presenting a bulk modulus of 143(1) GPa for both phlogopite polytypes (Fig.3a). This means that the TOT layer of the two minerals presented similar elastic properties. As also explained in previous works on phyllosilicates (Comodi et al., 2004;Chheda et al., 2014;Ulian and Valdrè, 2015a), the variation of the tetrahedral rotation angle α_r (Donnay et al., 1964a, b) is the main mechanism that accommodated the compression of the mineral. The α_r value is greater than 0 for both phlogopite polytypes at equilibrium, which means that the TO₄ mesh is not hexagonal (ideal) but forms a di-trigonal ring (see Fig.1), and α_r increases with pressure as reported in Fig.3b. This trend is in line with the simulations performed by Chheda and collaborators (2014), and with the general compressional behaviour of trioctahedral micas reported by several authors (Comodi and Zanazzi, 1995;Gatta et al., 2015;Gatta et al., 2011;Pavese et al., 1999;Pavese et al., 2003;Pawley et al., 2002). In addition, it is interesting to note that the octahedral flattening angle ψ became smaller by increasing pressure in the 1M polytype, whereas a direct proportionality was observed for the 2M₁ one.

Figure 3. (a) Variation as a function of pressure of the octahedral flattening angle ψ (red symbols) and tetrahedral rotation angle α_r (blue symbols) in the 1M and 2M₁ phlogopite models. (b) Evolution of the tetrahedral TO₄ and octahedral MO₆ polyhedral volume with pressure.

3.3 Elastic constants and derived quantities

To provide a complete analysis of the elastic behaviour of phlogopite polytypes, the second-order elastic tensor of the two models was also calculated. The approach relies on stress-strain relationships based on total energy E(V, ε) calculations through a Taylor expansion in terms of the strain components truncated at the second order:

where ε is the strain, σ is the stress, C is the 6 × 6 matrix representation of the elastic tensor (Voigt’s notation), V is the volume of the strained unit cell and V₀ is the volume at equilibrium. In this expression, α, β = 1, 2, 3, … 6. The adiabatic second-order elastic moduli are related to the second derivatives of the total energy E on the strain according to:

The C_αβ values were calculated by imposing a certain amount of strain ε along the crystallographic direction corresponding to the component of the dynamical matrix. This procedure is automated in the CRYSTAL code (Perger et al., 2009), and in the present simulations the default values that control the calculation of the elastic moduli were employed. To calculate the elastic tensor, the crystal was oriented with the b and c crystallographic axes parallel to the y and z Cartesian axes, respectively. Within this crystal orientation and considering the monoclinic lattice, the 6 × 6 matrix C of the elastic moduli in Voigt’s notation is given by:

whose values are reported in Table 4, considering both B3LYP and B3LYP-D* approaches. As expected, the inclusion of long-range interactions leads to a stiffer behaviour of phlogopite, especially for the diagonal terms C_αα.

Table 4. Elastic moduli C_αβ (in GPa) of phlogopite polytypes obtained in the present DFT simulations, compared to previous theoretical and experimental data.

Notes: GTO – Gaussian-type orbitals; PAW – projector-augmented wave; EWP – Elastic wave propagation. a – present work; b – Chheda et al. (2014); c – Aleksandrov et al. (1974).

Very few experimental data were reported in the scientific literature because suitable single crystals of phlogopite are rare, and many of these specimens are required for a complete characterisation of the elastic tensor of low-symmetry phases. Nevertheless, the present simulations are in good agreement with the experimental results of Alexandrov and Ryzhova (1961) and Aleksandrov et al. (1974), where the authors measured the elastic moduli from two crystals, one from the Slyudyanka and one from the Aldan region, near Lake Baikal (Russia). It is worth noting that only 5 independent elastic tensor components related to the hexagonal symmetry were reported instead of the 13 moduli required for the actual monoclinic symmetry of phlogopite. Also, two pairs of elastic moduli were equal, i.e., C₁₁ = C₂₂ and C₄₄ = C₅₅, and the off-diagonal components C₁₅, C₂₅, C₃₅ and C₄₆ were zero due to symmetry constraints. The present theoretical results at the B3LYP-D* level of theory indeed showed that C₁₁ ≈ C₂₂ and C₄₄ ≈ C₅₅, suggesting the near-isotropic elastic behaviour of the basal plane. Vaughan and Guggenheim (1986) measured the the full monoclinic elastic tensor for the dioctahedral equivalent of phlogopite, i.e., muscovite KAl₂Si₃AlO₁₀(OH)₂, finding a similar elastic behaviour with a small, but not insignificant, deviation from the hexagonal symmetry. Within this computational framework, the calculated C_αα moduli of phlogopite are overestimated by about 31 GPa (α = 1, 2, 3) and 13 GPa (α = 4, 5, 6) on average, but the trends of the diagonal components C₁₁ ≈ C₂₂ > C₃₃ and C₄₄ ≈ C₅₅ < C₆₆ are the same as those experimentally obtained from elastic wave propagation experiments (Aleksandrov et al., 1974;Alexandrov and Ryzhova, 1961). The uncorrected B3LYP functional resulted smaller absolute deviations from the experimental data (ca. 4 GPa and 2 GPa for the uniaxial and biaxial tensor components, respectively) and the same C₁₁ ≈ C₂₂ > C₃₃ relationship, however, the C₄₄ modulus was smaller than the C₅₅ one, and also much smaller than the one calculated at the B3LYP-D* level. Regarding the off-diagonal components, both approaches resulted in C₁₂ > C₁₃ although the uncorrected functional showed C₁₃ > C₂₃, while the inclusion of the van der Waals interactions provided C₁₃ ≈ C₂₃, which is in better agreement with the experimental evidence. This further confirms that the proper treatment of long-range interactions is relevant for the fundamental analysis of the elastic properties of layered silicates, especially when using Gaussian-type orbitals basis sets.

The present results at the B3LYP(-D*) level of theory are also in line with the simulations performed by Chheda and collaborators (2014) using LDA and GGA functionals and plane-waves basis sets, albeit some notable difference is present (see Table 4). Considering the cited work, it is known that the local density approximation generally underestimates lattice vectors and overestimates the forces, resulting in a stiffer behaviour of simulated phases. Conversely, GGA functionals are generally underbinding, hence the calculated elastic properties can be slightly underestimated. Then, hybrid functionals like B3LYP, which include a small fraction of the exact Hartree-Fock exchange to the total energy of the system, are expected to provide results that are between the LDA/GGA extremes if, and only if, other computational parameters (basis set type and quality, k-point sampling, converge criteria, etc.) are the same. Another source of discrepancy between the present elastic results and those of Chheda et al. (2014) is indeed the basis set type. Localised Gaussian-type orbitals generally provide slightly larger C_αβ values than those obtained with plane-waves basis sets in density functional theory simulation because of the Pulay forces that arise from the Hellmann-Feynman theorem and the use of atom-centred GTOs. As previously explained (Ulian et al., 2021;Ulian and Valdrè, 2018), the Pulay forces are associated with the basis set superposition error (BSSE), which is due to the incompleteness of GTOs and that artificially increases the calculated stress/force within a solid phase. The authors are aware that an approach like the geometrical counterpoise (gCP) method (Brandenburg et al., 2013;Kruse and Grimme, 2012) devised for molecule and molecular crystals could at least partially remove this BSSE-related error. However, the current implementation is limited to some general-purpose GTO basis sets that do not include the one employed in the present work, and a thorough analysis of its suitability for inorganic crystals and minerals is still ongoing (Ulian and Valdrè, 2024b). Besides the Pulay stress, the general overestimation of the theoretical results, both the present one and those of Chheda et al. (2014), could be also due to the athermal conditions, i.e., the DFT simulations were carried out at 0 K without zero-point energy and any thermal contribution, while the elastic wave propagation measurements are typically carried out at room temperature (about 300 K). Regarding the present work, the use of gCP and including temperature effect, e.g., with the quasi-harmonic approximation, can provide results that are in more agreement with the experimental findings. 

For practical applications, it is commonly preferred the use of the polycrystalline averages, namely the bulk (K) and shear (μ) moduli, calculated with the Voigt–Reuss–Hill approach (Hill, 1952;Nye, 1957) according to the formulas:

where S_ij = C_ij^–1 are the component of the compliance tensor, i.e., the inverse of elastic moduli tensor C, and the subscripts V and R indicate the Voigt (upper) and Reuss (lower) bound of the isotropic elastic properties. The Young's modulus E and the Poisson’s ratio υ were calculated from the following relations:

and

The polycrystalline mechanical properties were reported in Table 5. The K_R values obtained for both phlogopite 1M and 2M₁ polytypes were 58.8 GPa and 59.4 GPa, respectively, which are consistent with those calculated through the BM3 fitting procedure (57.9 GPa and 58.3 GPa, respectively).

Table 5. Isotropic elastic properties (bulk K, shear μ, and Young’s E moduli, and Poisson’s ratio υ) of phlogopite, alongside the universal anisotropy index A^U and the percentage of anisotropy of the bulk (A_K) and shear (A_μ) moduli.

a – present work; b – Chheda et al. (2014); c – Aleksandrov et al. (1974). Subscripts V and R indicate the Voigt (upper) and Reuss (lower) bounds, respectively. ^†Values calculated in the present work from the second-order elastic moduli reported by Aleksander et al. (1974).

The elastic anisotropy of phlogopite was described using the universal anisotropic index, A^U, as proposed by Ranganathan and Ostoja-Starzewski (2008):

For all but isotropic systems, for which A^U = 0, the ratio between the Voigt and Reuss moduli is not equal to unity, thus A^U > 0. The percentage anisotropy in compression (A_K) and shear (A_μ) was also calculated according to the following formulas (Ranganathan and Ostoja-Starzewski, 2008):

The anisotropies obtained via these expressions were reported in Table 5, showing that the A^U values of the two phlogopite polytypes are very similar (about 3.8 at the B3LYP-D* case), and between the values calculated from the theoretical results of Chheda and collaborators (2014) using the LDA (2.9) and GGA (4.4) functionals. There is a slight discrepancy between the different DFT approaches and the experimental results (A^U = 11.8), which is related to the high anisotropy of the shear moduli. It is suggested that the high μ_V/μ_R ratio is probably due to the very small values of the C₄₄ = C₅₅ elastic tensor components that were experimentally measured by Aleksandrov et al. (1974) with elastic wave propagation methods. Similar figures were found for the bulk and shear anisotropies, with A_K ≈ 14% and A_μ ≈ 26% calculated at the B3LYP-D* level of theory, intermediate values between those obtained using LDA and GGA functionals by Chheda and co-workers (2014). Conversely, the uncorrected B3LYP functional led to an anomalously high anisotropy index (35.2) due to the very large difference between the Voigt and Reuss shear moduli (A_μ ≈ 77%).

To further assess the validity of the linear model used to calculate the axial compressibilities, the axial moduli M(a), M(b) and M(c) were calculated from the elastic compliance tensor components according to the expressions reported by Mookherjee et al. (2016):

For instance, the calculated moduli were M(a) = 344.3 GPa, M(b) = 365.2 GPa and M(c) = 88.0 GPa for phlogopite-1M, and M(a) = 360.5 GPa, M(b) = 353.9 GPa and M(c) = 89.1 GPa for the 2M₁ polytype, values consistent with those reported in Table 3 obtained from the equation of state fit and that provide another measure of the anisotropy of the mineral, with a ratio M(a) : M(b) : M(c) ≈ 4 : 4 : 1 for both polytypes.

Directional single-crystal elastic properties (Young’s modulus E, linear compressibility β = K^–1, shear modulus μ and Poisson’s ratio υ) were calculated in the Cartesian space using the QUANTAS code (Gaillac et al., 2016;Ulian and Valdrè, 2022). The interested reader can find in the cited works all the formulations employed to calculate the spatial dependence of the cited properties. Differently from the other elastic properties reported above, where no significant variation between the 1M and 2M₁ phlogopite polytypes was evinced, the directional properties showed some effect due to the different TOT stacking (see Fig.S1 in the Supplementary Materials). For example, the shape of the Young’s modulus on the (xy) plane (Fig.S1a) is slightly more compressed along NE-SW and NO-SE directions in the polytype-1M, while it appeared more isotropic in phlogopite-2M₁. The shear moduli on the (xz) plane of the latter polytype (see Fig.S1c) is slightly canted, with the oblate maxima direction almost parallel to the east-west direction.

Finally, the phase (seismic) velocities were calculated by solving Christoffel’s equation (Musgrave, 1970) using the routines implemented in the QUANTAS code (Jaeken and Cottenier, 2016;Ulian and Valdrè, 2022). The results are shown in Fig.4, where the compressional wave velocity (v_P, P-wave, primary), the anisotropy of the shear velocities (S-wave, secondary), and the polarization of the fast shear wave (v_S1) are reported for both phlogopite polytypes. From the stereographic projections, particularly the S-wave anisotropy, it is immediately visible an east-west mirror plane normal to the b-axis that is related to the monoclinic crystal system of the mineral. This figure is in very good agreement with the theoretical analysis of Chheda et al. (2014), with some degree of deviation of the absolute values of the seismic velocities and anisotropy due to the stiffer elastic moduli obtained with the combination of the Gaussian-type orbitals basis sets and hybrid B3LYP-D* functional.

Figure 4. Lamber equal area (stereographic) projection of the upper hemisphere of the seismic anisotropy calculated for (a) 1M and (b) 2M₁ polytypes of phlogopite. In each panel, the left picture is related to the primary velocity, the central one is the S-wave anisotropy (related to the fast and slow shear waves) and the rightmost stereoplot is the polarization of the v_S1 (fast secondary) shear velocity. All the quantities were calculated at the DFT/B3LYP-D* level of theory.

Table 6 reports the maximum and minimum v_P values, the percentage of v_P anisotropy (Av_P) and the range of S-wave anisotropy (Av_S), as calculated from the present GTO/B3LYP(-D*) simulations for both phlogopite polytypes, together with previous experimental (Aleksandrov et al., 1974) and theoretical (Chheda et al., 2014) data related to the phlogopite 1M. The results are in general good agreement, with a slight overestimation of the maximum P-wave velocity and a smaller minimum primary wave velocity obtained from the uncorrected B3LYP functional. The shear wave anisotropy at the B3LYP-D* level of theory is within the values obtained by Chheda and collaborators (2014), with the maximum Av_S value being underestimated. This latter datum is instead more than double when calculated from the B3LYP results, and very overestimated compared to the experimental value reported by Aleksandrov et al. (1974).

Table 6. Maximum and minimum compressional wave velocity (v_P, km s^–1), P-wave anisotropy (Av_P, %), and maximum and minimum S-wave anisotropy (Av_S, %) of phlogopite.

a – present work; b – Chheda et al. (2014); c – Aleksandrov et al. (1974).

As mentioned in the introduction, the knowledge of the elastic properties could help explain and interpret some geophysical observations reported in previous analyses. For instance, a large delay time (δt) of about 1–2 s was measured in certain subduction zones, like the Ryukyu and Calabrian arcs, in addition to the trench-parallel shear wave splitting (Long and Silver, 2008). Chheda and collaborators (2014) suggested that peridotitic rocks with olivine as the major component would require a very large thickness ranging from 100 to 150 km to justify this large δt value. However, phyllosilicates like talc (Mainprice et al., 2008;Ulian et al., 2014), muscovite (Ulian and Valdrè, 2015a), antigorite (Mookherjee and Capitani, 2011) and phlogopite present significantly higher anisotropic behaviours than olivine. Thus, their presence in the subduction layer would require smaller thickness (~10–15 km) to explain both the trench parallel shear wave polarization anisotropy and the large delay time (Chheda et al., 2014). Indeed, According to previous investigations (Trønnes, 2002), phlogopite has a wide PT stability, thus it can be found in different settings, e.g., in the upper mantle and subduction zone, the latter being also characterised by hydrating conditions that could further stabilise this trioctahedral mica and other hydrous phases.

1. Conclusions

In the present work, using ab initio Density Functional Theory simulations, the elastic properties of two phlogopite polytypes, namely the 1M and the 2M₁, were investigated considering hydrostatic compression up to about 15 GPa and the monoclinic elastic tensor. The following conclusions can be drawn:

(1) despite the different stacking of the TOT layers, the elastic properties of phlogopite were not severely affected by polytypism. Considering the numerical errors of the fitting procedure, the bulk modulus K₀ obtained from the third-order Birch-Murnaghan equation of state was in the range 55 – 60 GPa for both polytypes.

(2) the compressional behaviour of phlogopite is highly anisotropic, with the a and b crystallographic axes being about four times less compressible than the c one, as noticed from the linear moduli. This agrees with the experimental observations made on dioctahedral and trioctahedral micas, and it is due to several factors, i.e., the high compressibility of the interlayer region, the low compressibility of the tetrahedral and octahedral polyhedral, and the accommodation of the stress within the TOT layers via variation of the tetrahedral rotation angle α_r.

(3) the elastic tensors of phlogopite-1M and phlogopite-2M₁ have comparable C_ij components (in Voigt’s notation), with a difference in sign and absolute value in the C₂₅ one. This is probably the reason of the small differences observed in the directionally dependent elastic properties.

(4) the DFT simulations were able to correctly describe the behaviour of the known 1M polytype of phlogopite, extending the knowledge to the 2M1 one. In particular, the use of the DFT-D2 correction provided a robust description of the physical forces acting along [001], namely the stacking direction of the TOT + i layers.

It is also known that the elastic properties of micas can be deeply affected by the presence of cationic/anionic substitutions in the T and O sheets, e.g., Fe^II substituting Mg^II. These are currently under investigation and will be the subject of a future work.

From the mineralogical, petrological, and geophysical perspectives, the knowledge of the physical properties of phlogopite is of utmost importance for deepening the understanding of the thermodynamic and kinetic stability of altered rock assemblages, such as the mica–amphibole–rutile–ilmenite–diopside system. Indeed, the phlogopite chemical composition is similar to that of the ternary system olivine–quartz–alkali, which is also close to the bulk composition of the cited altered rocks. Also, the elasticity of the trioctahedral mica is relevant both for bulk and 2D applications of phlogopite in materials science, owing to the dielectric properties of the mineral that can be tuned by the exfoliation in the basal plane and cationic substitutions, particularly in the octahedral sheet. First principles simulations can help in shedding light on these topics, providing useful resources for the interpretation of geophysical phenomena and the development of more, sustainable materials for advanced applications.

Acknowledgements

The authors wish to thank the University of Bologna for supporting the present research. All the simulations were performed with the computational resources (HPC cluster) and software license of the Interdisciplinary Research Centre of Biomineralogy, Crystallography and Biomaterials, Department of Biological, Geological and Environmental Sciences, University of Bologna.

Competing Interests

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Author Contributions

Conceptualization, G.U. and G.V.; methodology, G.U. and F.R.; validation, G.U. and G.V.; formal analysis, G.U.; investigation, G.U. and G.V.; data curation, G.U. and F.R.; writing—review and editing, G.U., F.R. and G.V.; visualization, G.U. and F.R.; supervision, G.V. All authors have read and agreed to the published version of the manuscript.

Data Availability

The results of the present work are reported in the manuscript and in a dedicated dataset published at the following link: https:/doi.org/10.17632/r4wmxz7kcc.1.

References

Aleksandrov, K. S., Alchikov, U. V., Belikov, B. P., Zaslavskii, B. I., and Krupnyi, A. I.: Velocities of elastic waves in minerals at atmospheric pressure and increasing precision of elastic constants by means of EVM [in Russian]. In: Izv. Acad. Sci. USSR, Geol. Ser., 1974.

Alexandrov, K. S., and Ryzhova, T. V.: Elastic properties of rock-forming minerals. II. Layered silicates, Bulletin of the Academy of Sciences of the U.S.S.R, Geophysics Series, 9, 1165-1168, 1961.

Angel, R. J., Gonzalez-Platas, J., and Alvaro, M.: EosFit7c and a Fortran module (library) for equation of state calculations, Zeitschrift Fur Kristallographie, 229, 405-419, 10.1515/zkri-2013-1711, 2014.

Ariane, K., Tamayo, A., Chorfa, A., Rubio, F., and Rubio, J.: Optimization of the nucleating agent content for the obtaining of transparent fluormica glass-ceramics, Ceramics International, 49, 9826-9838, 10.1016/j.ceramint.2022.11.156, 2023.

Becke, A. D.: Density-Functional Thermochemistry .3. The Role of Exact Exchange, Journal of Chemical Physics, 98, 5648-5652, 10.1063/1.464913, 1993.

Belmonte, D.: First Principles Thermodynamics of Minerals at HP-HT Conditions: MgO as a Prototypical Material, Minerals, 7, 183, 10.3390/Min7100183, 2017.

Birch, F.: Finite elastic strain of cubic crystal, Physical Review, 71, 809-824, 1947.

Brandenburg, J. G., Alessio, M., Civalleri, B., Peintinger, M. F., Bredow, T., and Grimme, S.: Geometrical correction for the inter- and intramolecular basis set superposition error in periodic density functional theory calculations, Journal of Physical Chemistry A, 117, 9282-9292, 10.1021/jp406658y, 2013.

Broyden, C. G.: The convergence of a class of double-rank minimization algorithms: 1. General considerations, IMA Journal of Applied Mathematics, 6, 76-90, 10.1093/imamat/6.1.76, 1970a.

Broyden, C. G.: The convergence of a class of double-rank minimization algorithms: 2. The new algorithm, IMA Journal of Applied Mathematics, 6, 222-231, 10.1093/imamat/6.3.222, 1970b.

Cadore, A. R., De Oliveira, R., Longuinhos, R., Teixeira, V. D. C., Nagaoka, D. A., Alvarenga, V. T., Ribeiro-Soares, J., Watanabe, K., Taniguchi, T., Paniago, R. M., Malachias, A., Krambrock, K., Barcelos, I. D., and De Matos, C. J. S.: Exploring the structural and optoelectronic properties of natural insulating phlogopite in van der Waals heterostructures, 2D Materials, 9, 035007, 10.1088/2053-1583/ac6cf4, 2022.

Catti, M., Ferraris, G., Hull, S., and Pavese, A.: Powder Neutron-Diffraction Study of 2M₁ Muscovite at Room Pressure and at 2 Gpa, European Journal of Mineralogy, 6, 171-178, 1994.

Chheda, T. D., Mookherjee, M., Mainprice, D., dos Santos, A. M., Molaison, J. J., Chantel, J., Manthilake, G., and Bassett, W. A.: Structure and elasticity of phlogopite under compression: Geophysical implications, Physics of the Earth and Planetary Interiors, 233, 1-12, 10.1016/j.pepi.2014.05.004, 2014.

Civalleri, B., Zicovich-Wilson, C. M., Valenzano, L., and Ugliengo, P.: B3LYP augmented with an empirical dispersion term (B3LYP-D*) as applied to molecular crystals, CrystEngComm, 10, 405-410, 10.1039/b715018k, 2008.

Comodi, P., and Zanazzi, P. F.: High-Pressure Structural Study of Muscovite, Physics and Chemistry of Minerals, 22, 170-177, 1995.

Comodi, P., Fumagalli, P., Montagnoli, M., and Zanazzi, P. F.: A single-crystal study on the pressure behavior of phlogopite and petrological implications, American Mineralogist, 89, 647-653, 10.2138/am-2004-0420, 2004.

Dion, M., Rydberg, H., Schröder, E., Langreth, D. C., and Lundqvist, B. I.: Van der Waals Density Functional for General Geometries, Physical Review Letters, 92, 246401, 10.1103/PhysRevLett.92.246401, 2004.

Donnay, G., Morimoto, N., and Takeda, H.: Trioctahedral one-layer micas. I. Crystal structure of a synthetic iron mica, 17, 1369-1373, https://doi.org/10.1107/S0365110X64003450, 1964a.

Donnay, G., Morimoto, N., and Takeda, H.: Trioctahedral one-layer micas. II. Prediction of the structure from composition and cell dimensions, 17, 1374-1381, https://doi.org/10.1107/S0365110X64003462, 1964b.

Dovesi, R., Roetti, C., Freyria Fava, C., Prencipe, M., and Saunders, V. R.: On the elastic properties of lithium, sodium an potassium oxide. An ab initio study., Chemical Physics, 156, 11-19, 1991.

Dovesi, R., Erba, A., Orlando, R., Zicovich-Wilson, C. M., Civalleri, B., Maschio, L., Rerat, M., Casassa, S., Baima, J., Salustro, S., and Kirtman, B.: Quantum-mechanical condensed matter simulations with CRYSTAL, Wiley Interdisciplinary Reviews-Computational Molecular Science, 8, E1360, 10.1002/Wcms.1360, 2018.

Fletcher, R.: A new approach to variable metric algorithms, The Computer Journal, 13, 317-322, 10.1093/comjnl/13.3.317, 1970.

Gaillac, R., Pullumbi, P., and Coudert, F. X.: ELATE: an open-source online application for analysis and visualization of elastic tensors, Journal of Physics-Condensed Matter, 28, 275201, 10.1088/0953-8984/28/27/275201, 2016.

Gatta, G. D., Merlini, M., Rotiroti, N., Curetti, N., and Pavese, A.: On the crystal chemistry and elastic behavior of a phlogopite 3T, Physics and Chemistry of Minerals, 38, 655-664, 10.1007/s00269-011-0438-z, 2011.

Gatta, G. D., Lotti, P., Merlini, M., Liermann, H.-P., Lausi, A., Valdrè, G., and Pavese, A.: Elastic behaviour and phase stability of pyrophyllite and talc at high pressure and temperature, Physics and Chemistry of Minerals, 42, 309-318, 10.1007/s00269-014-0721-x, 2015.

Gatti, C., Saunders, V. R., and Roetti, C.: Crystal-field effects on the topological properties of the electron-density in molecular-crystals - the case of urea, Journal of Chemical Physics, 101, 10686-10696, 10.1063/1.467882, 1994.

Geng, X., Tian, S., Xu, W., Chen, L., Lu, N., Liang, Z., Hu, W., and Liu, J.: A Two-Stage Geodynamic Model for Post-Collisional Potassic-Ultrapotassic Magmatism in Southeast Tibet, Journal of Geophysical Research: Solid Earth, 129, 10.1029/2024JB028887, 2024.

Goldfarb, D.: A family of variable-metric methods derived by variational means, Mathematics of Computation, 24, 23-26, 10.1090/S0025-5718-1970-0258249-6 1970.

Grimme, S.: Semiempirical GGA-type density functional constructed with a long-range dispersion correction, Journal of Computational Chemistry, 27, 1787-1799, 10.1002/jcc.20495, 2006.

Guidotti, C. V., and Sassi, F. P.: Muscovite as a petrogenetic indicator mineral in pelitic schists, Neues Jahrbuch fuer Mineralogie, Abhandlungen, 123, 97-142, 1976.

Hebbache, M., and Zemzemi, M.: Ab initio study of high-pressure behavior of a low compressibility metal and a hard material: Osmium and diamond, Physical Review B, 70, 10.1103/Physrevb.70.224107, 2004.

Hill, R.: The elastic behaviour of a crystalline aggregate, Proceedings of the Physical Society, London, Section A, 65, 349-354, 1952.

Icenhower, J., and London, D.: An experimental study of element partitioning among biotite, muscovite, and coexisting peraluminous silicic melt at 200 MPa (H₂O), American Mineralogist, 80, 1229-1251, 10.2138/am-1995-11-1213, 1995.

Jaeken, J. W., and Cottenier, S.: Solving the Christoffel equation: Phase and group velocities, Computer Physics Communications, 207, 445-451, 10.1016/j.cpc.2016.06.014, 2016.

King, T. T., Grayeski, W., and Cooper, R. F.: Thermochemical reactions and equilibria between fluoromicas and silicate matrices: A petromimetic perspective on structural ceramic composites, Journal of the American Ceramic Society, 83, 2287-2296, 10.1111/j.1151-2916.2000.tb01549.x, 2000.

Kruse, H., and Grimme, S.: A geometrical correction for the inter- and intra-molecular basis set superposition error in Hartree-Fock and density functional theory calculations for large systems, Journal of Chemical Physics, 136, 10.1063/1.3700154, 2012.

Lacalamita, M., Mesto, E., Scordari, F., and Schingaro, E.: Chemical and structural study of 1M- and 2M (1)-phlogopites coexisting in the same Kasenyi kamafugitic rock (SW Uganda), Physics and Chemistry of Minerals, 39, 601-611, 10.1007/s00269-012-0515-y, 2012.

Laurora, A., Brigatti, M. F., Mottana, A., Malferrari, D., and Caprilli, E.: Crystal chemistry of trioctahedral micas alkaline and subalkaline volcanic rocks: A case study from Mt. Sassetto (Tolfa district, central Italy), American Mineralogist, 92, 468-480, 10.2138/am.2007.2339, 2007.

Lee, C. T., Yang, W. T., and Parr, R. G.: Development of the Colle-Salvetti Correlation-Energy Formula into a Functional of the Electron-Density, Physical Review B, 37, 785-789, 10.1103/PhysRevB.37.785, 1988.

Long, M. D., and Silver, P. G.: The subduction zone flow field from seismic anisotropy: a global view, Science, 319, 315-318, 10.1126/science.1150809, 2008.

Mainprice, D., Le Page, Y., Rodgers, J., and Jouanna, P.: Ab initio elastic properties of talc from 0 to 12 GPa: Interpretation of seismic velocities at mantle pressures and prediction of auxetic behaviour at low pressure, Earth and Planetary Science Letters, 274, 327-338, 10.1016/j.epsl.2008.07.047, 2008.

Melzer, S., and Wunder, B.: K-Rb-Cs partitioning between phlogopite and fluid: Experiments and consequences for the LILE signatures of island arc basalts, Lithos, 59, 69-90, 10.1016/S0024-4937(01)00061-5, 2001.

Mookherjee, M., and Capitani, G. C.: Trench parallel anisotropy and large delay times: Elasticity and anisotropy of antigorite at high pressures, Geophysical Research Letters, 38, L09315, 10.1029/2011gl047160, 2011.

Momma, K., and Izumi, F.: VESTA: a three-dimensional visualization system for electronic and structural analysis, J. Appl. Crystallogr., 41, 653-658, 2008.

Mookherjee, M., Tsuchiya, J., and Hariharan, A.: Crystal structure, equation of state, and elasticity of hydrous aluminosilicate phase, topaz-OH (Al2SiO4(OH)(2)) at high pressures, Physics of the Earth and Planetary Interiors, 251, 24-35, 10.1016/j.pepi.2015.11.006, 2016.

Musgrave, M. J. P.: Crystal Acoustics: introduction to the study of elastic waves and vibrations in crystals, Holden-Day, San Francisco, CA, USA, 1970.

Nada, R., Nicholas, J. B., McCarthy, M. I., and Hess, A. C.: Basis sets for ab initio periodic Hartree-Fock studies of zeolite/adsorbate interactions: He, Ne, and Ar in silica sodalite, International Journal of Quantum Chemistry, 60, 809-820, 10.1002/(SICI)1097-461X(1996)60:4<809::AID-QUA3>3.0.CO;2-0, 1996.

Nye, J. F.: Physical properties of crystals, Oxford University Press, Oxford, 1957.

Ottonello, G., Civalleri, B., Ganguly, J., Perger, W. F., Belmonte, D., and Zuccolini, M. V.: Thermo-chemical and thermo-physical properties of the high-pressure phase anhydrous B (Mg₁₄Si₅O₂₄): An ab-initio all-electron investigation, American Mineralogist, 95, 563-573, 10.2138/am.2010.3368, 2010.

Pavese, A., Ferraris, G., Pischedda, V., and Mezouar, M.: Synchrotron powder diffraction study of phengite 3T from the Dora-Maira massif: P-V-T equation of state and petrological consequences, Physics and Chemistry of Minerals, 26, 460-467, DOI 10.1007/s002690050208, 1999.

Pavese, A., Levy, D., Curetti, N., Diella, V., Fumagalli, P., and Sani, A.: Equation of state and compressibility of phlogopite by in-situ high-pressure X-ray powder diffraction, European Journal of Mineralogy, 15, 455-463, 10.1127/0935-1221/2003/0015-0455, 2003.

Pawley, A. R., Clark, S. M., and Chinnery, N. J.: Equation of state measurements of chlorite, pyrophyllite, and talc, American Mineralogist, 87, 1172-1182, 2002.

Perger, W. F., Criswell, J., Civalleri, B., and Dovesi, R.: Ab-initio calculation of elastic constants of crystalline systems with the CRYSTAL code, Computer Physics Communications, 180, 1753-1759, DOI 10.1016/j.cpc.2009.04.022, 2009.

Pfrommer, B. G., Côté, M., Louie, S. G., and Cohen, M. L.: Relaxation of Crystals with the Quasi-Newton Method, Journal of Computational Physics, 131, 233-240, 10.1006/jcph.1996.5612, 1997.

Prencipe, M., Noel, Y., Bruno, M., and Dovesi, R.: The vibrational spectrum of lizardite-1T [Mg(3)Si(2)O(5)(OH)(4)] at the Gamma point: A contribution from an ab initio periodic B3LYP calculation, American Mineralogist, 94, 986-994, 2009.

Ranganathan, S. I., and Ostoja-Starzewski, M.: Universal elastic anisotropy index, Physical Review Letters, 101, 055504, 10.1103/PhysRevLett.101.055504, 2008.

Shanno, D. F.: Conditioning of quasi-Newton methods for function minimization, Mathematics of Computation, 24, 647-656, 10.1090/S0025-5718-1970-0274029-X 1970.

Sudo, A., and Tatsumi, Y.: Phlogopite and K‐amphibole in the upper mantle: Implication for magma genesis in subduction zones, Geophysical Research Letters, 17, 29-32, 10.1029/GL017i001p00029, 1990.

Sweeney, R. J., Thompson, A. B., and Ulmer, P.: Phase relations of a natural MARID composition and implications for MARID genesis, lithospheric melting and mantle metasomatism, Contributions to Mineralogy and Petrology, 115, 225-241, 10.1007/BF00321222, 1993.

Trønnes, R. G.: Stability range and decomposition of potassic richterite and phlogopite end members at 5-15 GPa, Mineralogy and Petrology, 74, 129-148, 10.1007/s007100200001, 2002.

Tutti, F., Dubrovinsky, L. S., and Nygren, M.: High-temperature study and thermal expansion of phlogopite, Physics and Chemistry of Minerals, 27, 599-603, 10.1007/s002690000098, 2000.

Ulian, G., Tosoni, S., and Valdrè, G.: Comparison between Gaussian-type orbitals and plane wave ab initio density functional theory modeling of layer silicates: Talc Mg₃Si₄O₁₀(OH)₂ as model system, Journal of Chemical Physics, 139, 204101, 10.1063/1.4830405, 2013.

Ulian, G., Tosoni, S., and Valdrè, G.: The compressional behaviour and the mechanical properties of talc [Mg₃Si₄O₁₀(OH)₂]: a density functional theory investigation, Physics and Chemistry of Minerals, 41, 639-650, 10.1007/s00269-014-0677-x, 2014.

Ulian, G., and Valdrè, G.: Density functional investigation of the thermo-physical and thermo-chemical properties of 2M(1) muscovite, American Mineralogist, 100, 935-944, 10.2138/am-2015-5086, 2015a.

Ulian, G., and Valdrè, G.: Structural, vibrational and thermophysical properties of pyrophyllite by semi-empirical density functional modelling, Physics and Chemistry of Minerals, 42, 609-627, 10.1007/s00269-015-0748-7, 2015b.

Ulian, G., and Valdrè, G.: Density functional investigation of the thermophysical and thermochemical properties of talc Mg3Si4O10(OH)(2), Physics and Chemistry of Minerals, 42, 151-162, 10.1007/s00269-014-0708-7, 2015c.

Ulian, G., and Valdrè, G.: Second-order elastic constants of hexagonal hydroxylapatite (P6₃) from ab initio quantum mechanics: comparison between DFT functionals and basis sets, International Journal of Quantum Chemistry, 118, e25500, 10.1002/qua.25500, 2018.

Ulian, G., and Valdrè, G.: Equation of state and second-order elastic constants of portlandite Ca(OH)₂ and brucite Mg(OH)₂, Physics and Chemistry of Minerals, 46, 101-117, 10.1007/s00269-018-0989-3, 2019.

Ulian, G., Moro, D., and Valdrè, G.: Elastic properties of heterodesmic composite structures: The case of calcite CaCO₃ (space group R-3c), Composites Part C: Open Access, 6, https://doi.org/10.1016/j.jcomc.2021.100184, 2021.

Ulian, G., and Valdrè, G.: QUANTAS, a Python software for the analysis of solids from ab initio quantum mechanical simulations and experimental data, J. Appl. Crystallogr., 55, 386-396, 10.1107/S1600576722000085, 2022.

Ulian, G., and Valdrè, G.: The effect of long-range interactions on the infrared and Raman spectra of aragonite (CaCO₃, Pmcn) up to 25 GPa, Scientific Reports, 13, 2725, 10.1038/s41598-023-29783-7, 2023a.

Ulian, G., and Valdrè, G.: Crystal-chemical, vibrational and electronic properties of 1M-phlogopite K(Mg,Fe)₃Si₃AlO₁₀(OH)₂ from Density Functional Theory simulations, Applied Clay Science, 246, 107166, 10.1016/j.clay.2023.107166, 2023b.

Ulian, G., and Valdrè, G.: SEISMIC, a Python-based code of the QUANTAS package to calculate the phase and group acoustic velocities in crystals, Computers and Geosciences, 188, 10.1016/j.cageo.2024.105615, 2024a.

Ulian, G., and Valdrè, G.: Crystal structure and elastic and phonon properties of realgar versus pressure, J. Appl. Crystallogr., 57, 220-231, 10.1107/S1600576724000025, 2024b.

Ulian, G., and Valdrè, G.: Crystallographic, electronic and vibrational properties of 2D silicate monolayers, J. Appl. Crystallogr., 58, 349-362, 10.1107/S1600576725000731, 2025.

Valenzano, L., Torres, F. J., Doll, K., Pascale, F., Zicovich-Wilson, C. M., and Dovesi, R.: Ab initio study of the vibrational spectrum and related properties of crystalline compounds; the case of CaCO₃ calcite, Zeitschrift Fur Physikalische Chemie-International Journal of Research in Physical Chemistry & Chemical Physics, 220, 893-912, 10.1524/zpch.2006.220.7.893, 2006.

Valenzano, L., Noel, Y., Orlando, R., Zicovich-Wilson, C. M., Ferrero, M., and Dovesi, R.: Ab initio vibrational spectra and dielectric properties of carbonates: magnesite, calcite and dolomite, Theoretical Chemistry Accounts, 117, 991-1000, 10.1007/s00214-006-0213-2, 2007.

Van Reenen, D. D., Smit, C. A., Huizenga, J. M., Tsunogae, T., and Safonov, O.: Thermo-tectonic evolution of the Neoarchaean Southern Marginal Zone of the Limpopo granulite Complex (South Africa), South African Journal of Geology, 126, 373-406, 10.25131/sajg.126.0027, 2023.

Vaughan, M. T., and Guggenheim, S.: Elasticity of muscovite and its relationship to crystal structure, Journal of Geophysical Research: Solid Earth, 91, 4657-4664, https://doi.org/10.1029/JB091iB05p04657, 1986.

Ventruti, G., Levy, D., Pavese, A., Scordari, F., and Suard, E.: High-temperature treatment, hydrogen behaviour and cation partitioning of a Fe-Ti bearing volcanic phlogopite by in situ neutron powder diffraction and FTIR spectroscopy, European Journal of Mineralogy, 21, 385-396, 10.1127/0935-1221/2009/0021-1903, 2009.

Virgo, D., and Popp, R. K.: Hydrogen deficiency in mantle-derived phlogopites, American Mineralogist, 85, 753-759, 10.2138/am-2000-5-614, 2000.

Wang, J., Wang, Q., Ma, L., Hu, W. L., Wang, J., Belousova, E., and Tang, G. J.: Rapid Recycling of Subducted Sediments in the Subcontinental Lithospheric Mantle, Journal of Petrology, 64, 10.1093/petrology/egad056, 2023.

We thank the referee for the positive comments on our work. Please, find in the following the point-by-point answers to each concern.

1) The mismatch between the calculated volumes and experimental results becomes increasingly large with the increasing pressure. Both structures are stronger in calculations than in experiments. What could be the reason?

Answer: we thank the referee for pointing this out. The main reason of this behaviour is due to the computational framework. We performed DFT simulations in the so-called “static” conditions, meaning that zero-point and thermal effects are neglected. Thus, the minerals exhibited this generally stiffer behaviour that increases with pressure.

2) The authors use the B3LYP-D* hybrid functional with Grimme D2 corrections, which is justified and in line with previous studies. However, it would be beneficial to include a short discussion on the expected accuracy of the B3LYP-D* method for high-pressure phases, referencing recent benchmarks on similar layered silicates.

Answer: according to the referee’s suggestion, we included in Section 2, after the description of the B3LYP-D* approach, some references related to the previous works on talc, pyrophyllite, and muscovite. In the cited references, the accuracy was very high, with results comparable to the experimental data.

3) The analysis of the seismic anisotropy (Fig. 4) is clear and well executed. The authors might consider a more explicit discussion of the potential implications for seismic interpretations in phlogopite-bearing mantle regions or subduction zones, where polytype coexistence could affect anisotropic wave propagation.

Answer: following the referee’s suggestion, we included a small paragraph discussing the possible implications associated with the co-presence of different polytypes, with slightly different off-diagonal elastic moduli, in mantle regions and subduction zones.