Mechanisms and conditions of Ti and Cr incorporation in mantle phlogopite: the results of atomistic simulation

Modeling of eight mechanisms for the incorporation of Ti4+ and Cr3+ impurity components into phlogopite was carried out by a semi-empirical method using the GULP (General Utility Lattice Program) software. The calculation of thermodynamic mixing properties in the range of 1–7 GPa and 373–1573 K and the analysis of the structure geometry for the simulated solid solutions provided the following energy-preferred schemes of isomorphic substitution: VI(Mg2+) + 2IV(Si4+) = VI(Ti4+) + 2IV(Al3+) and VI(Mg2+) + 2IV(Al3+) = VI(□) + 2IV(Ti4+), VI(Mg2+) + IV(Si4+) = VI(Cr3+) + IV(Al3+), and 3VI(Mg2+) = VI(Al3+) + VI(Cr3+) + VI(□). It is shown the scheme 2VI(Mg2+) = VI(Ti4+) + VI(□) illustrating entrance of Ti with the formation of a vacancy is realized in the case of microconcentrations of Ti only. Accumulation of high-Ti contents associates with the formation of a vacancy in the octahedral site. This provides incorporation of Ti via the schemes VI(Mg2+) + 2IV(Al3+) = VI(□) + 2IV(Ti4+) and (Mg, Fe2+) + 2OH− = Ti4+  + 2O2− only. It is shown that incorporation of high-Cr concentrations (> 5.5 wt % Cr2O3) is accompanied by an increase in the number of vacancies in the octahedral site with an increase in the proportion of the dioctahedral component K(Al, Cr, □)2AlSi3O10(OH)2.


Introduction
Phlogopite is a trioctahedral mica with a structure composed of three-layer packages of two tetrahedral layers [(AlSi 3 O 10 ) 5− ] interconnected by cations (Fe, Mg) 2+ in the octahedral sites. The three-layer packages are bound by singly charged cations (K + , Na + ) with a coordination number of 12.
Phlogopite is a common mineral in alkaline basalt, ultrapotassic rocks, kimberlite, lamproite, lamprophyre, granulite, carbonatite (Brod et al. 2001), and rocks of ophiolite complexes (Peng et al. 1995). In addition, phlogopite is mica of metasomatized mantle peridotite. At the same time, phlogopite inclusions in diamonds are rarely considered as an indicator of assemblages, since such finds are extremely rare being often interpreted as epigenetic (Meyer 1987;Harris 1992). Nevertheless, classification of diamonds into eclogitic and ultramafic (dunite-harzburgite, wehrlite, and lherzolite) associations should account for minor elements rather than major components in the composition of minerals only. As was shown in (Erlank et al. 1987;Sobolev et al. 1997Sobolev et al. , 1998Sobolev et al. , 2009) the content of such elements as Ti and Cr in phlogopite syngenetic with diamond may indicate that the mineral belongs to the eclogitic or peridotitic assemblage.
The incorporation of Ti 4+ and Cr 3+ in phlogopite and some other trioctahedral micas proceeds via various substitution mechanisms involving both tetrahedral and octahedral sites ( Table 1). The most common substitution mechanisms (1-4) (Dymek 1983;Henry et al. 2005) imply the presence of Al in the tetrahedral site of phlogopite, which is characteristic of igneous rocks including high-pressure ones (Arima and Edgar 1981). The schemes (5), (6), and (7), which are a combination of schemes (1-4) and two major mechanisms that provide incorporation of Al into the octahedral site, VI (R 2+ ) + IV (Si 4+ ) = VI (Al 3+ ) + IV (Al 3+ ) and 3 VI (R 2+ ) = 2 VI (Al 3+ ) + VI (□), are discussed as well and considered as basic ones for Mg-Fe micas from hightemperature metapelites (Thu et al. 2016). The scheme (8) 1 3 8 Page 2 of 15 of subordinate deprotonation substitution (Ti-oxy substitution) is considered to be one of the main mechanisms for the incorporation of Ti into phlogopite at high temperatures of the granulite facies (Dymek 1983;Henry et al. 2005), which has recently received special attention (e.g., (Ventruti et al. 2020) and references therein).
In this paper, we discuss the schemes (1-4) only to illustrate incorporation of Ti 4+ into phlogopite (Table 1), since we focus on high-pressure Ti-bearing mantle phlogopite, which is not characterized by the high concentrations of Al in the octahedral site (Arima and Edgar 1981).
The mechanisms of Cr 3+ incorporation into phlogopite are poorly discussed in literature. Basically, only the schemes (9) and (10) are considered (e.g., (Ferenc et al. 2016) for high-chromium trioctahedral micas (up to 6.54 wt % Cr 2 O 3 ) from Muranska Zdychava listvenite). Cr-bearing phlogopites are reported as inclusions in diamonds as well (up to 3.5 wt % Cr 2 O 3 ) (Sobolev et al. 2009) and in other rocks of different geodynamic environments, in which the concentration of Cr 2 O 3 impurity reaches 2.3 wt % Cr 2 O 3 (Delaney et al. 1980;Shee 1985;Lorand and Cottin 1987;Dawson 2002;Bosi et al. 2012;Guiliani et al. 2016). Thus, in addition to the schemes (9) and (10), here we study the schemes (11) and (12) and simulate two new end-members, using the same approach to Cr 3+ ion incorporation as for Ti end-members. Note that the scheme (11) suggests accumulation of Al 3+ in the octahedral site in association with Cr 3+ entrance. These schemes are acceptable from the point of the local charge balance and correlate with the Ti incorporation schemes; however, additional factors, such as the preference energy for the octahedral coordination (McClure 1957, may impose extra restrictions on such mechanisms.
It is important that the choice of one or another scheme dominating for phlogopite is related to the coexisting mineral assemblage that controls the composition of the studied phase. Due to the wide variations in natural and experimental mineral associations, and, consequently, in the composition of phlogopite in terms of the content of both major and minor components, the data on preferable mechanisms of Ti and Cr incorporation into phlogopite remain ambiguous.
The inconsistency of experimental information, which makes possible to give preference to one or another mechanism for Ti and Cr incorporation into phlogopite, has led to necessity of theoretical models for determination of the most preferable defect formation scheme. The aims of this study are (1) to identify the most realistic schemes for the incorporation of impurity components into phlogopite and (2) to estimate the maximum concentrations of Ti and Cr in the mineral by atomistic computer simulation under various thermodynamic conditions of the Earth's upper mantle.

Simulation methodology
The mechanisms of minor element incorporation into the phlogopite were simulated using the GULP (General Utility Lattice Program) software (Gale and Rohl 2003), which is based on the method of minimizing the structural energy of a crystal using semi-empirical interatomic potentials.
Modeling of structural defects in phlogopite was carried out in two different ways, which allowed us to control the reproducibility and correctness of the results. The first approach was to use the standard procedure for calculation of point defects by the Mott-Littleton two region strategy (Mott and Littleton 1938) implemented in the GULP software. Such approach suggests that a point defect in the crystal structure and the immediate region around the defect defined by a sphere of radius r 1 , participate in the procedure of minimizing the energy of interatomic interaction with the symmetry restrictions completely removed. A screening layer defined by the radius r 2 is introduced between the defect area and the periodic crystal. For a correct calculation of a point defect, it is usually sufficient to use the r 1 value of ~ 6-7 Å (150-300 atoms) and the thickness of the screening layer (r 2 -r 1 ) of 10 Å (1500-2000 atoms). However, the need to calculate the associates of point defects given in Table 1 required an increase in the dimension of these parameters to 10 and 20 Å, respectively. As is shown by our test simulations, the use of spheres of a larger radius turned out to be inexpedient, since in this case the calculation time only increased significantly, with an extremely slight change in the simulation accuracy.
The second way to study of the defect energy in the crystal structure was to simulate the associate in the central regions of supercells of different dimensions of the model structure of phlogopite with the non-translational symmetry removed (the calculation was carried out within the space group P1) directly. Two types of supercells were used: 4 × 2 × 2 (704 atoms) and 6 × 3 × 3 (2376 atoms) (see supplementary information. Fig. S1).
The calculations of defect complexes upon simulation using supercells were carried out for several topologically different configurations of defects separated from each other over the space of the supercell by a certain distance. The choice of topological configurations was based on the selection of atoms with such coordinates that the distance between them was minimal (from 3.2 to 5.5 Å).

Results and discussion
The reliability of model potential The results of the simulation of the crystal structure of phlogopite in comparison with the experimental values from (Redhammer and Roth 2002) are shown in supplementary information (Table S1). As is evident from the table, the unit-cell volume is reproduced quite accurately, and the deviations of the cell parameters are characterized by an error of ≤ 2.7%. The parameters of compressibility, compressional (V P ) and shear (V S ) sound velocities, and structural patterns (unit-cell parameters and volume) of phlogopite were studied in the pressure range from 0 to 12 GPa and compared with the published data (Pavese et al. 2003;Comodi et al. 2004;Chheda et al. 2014). These results are shown in figures ( Fig. 1 and Supplementary information Fig. S2) and demonstrate a good reproducibility. Change in the unit-cell volume is within the limit of estimates available from literature; the bulk and shear modulus retain the slope reported by Chheda et al. (2014), and there is complete agreement with the results (Chheda et al. 2014) for sound velocities (Fig. 1). The results of calculations of change in the cell parameters a and b (Å) are slightly overestimated; however, an error does not exceed 0.06 Å (Fig. S2), which is very close to those obtained by Chheda et al. (2014) using the LDA (local density approximation) method. Thus, the selected set of potentials reproduces the structural patterns and compressibility parameters of phlogopite correctly under the high-pressure conditions (up to 12 GPa) in the range of confidence intervals of quantum mechanical calculations. This allowed us to use the selected set of potentials for a comparative analysis of the energies of defect formation in the crystal structure according to the isomorphic schemes described above.

Mixing properties of simulated Tiand Cr-phlogopite end-members
The results obtained allowed us to estimate thermodynamic mixing properties of solid solution and isomorphic capacity of phlogopite for Ti 4+ and Cr 3+ . ∆H mix and ∆S mix required for calculation of the isomorphism limits as a function of temperature and pressure may be determined using the methods of computer simulation. The mixing entropy ∆S mix was calculated in the pressure range of 1-7 GPa and temperatures of 373-1573 K for each series of solid solutions via the formula: where S conf is the configuration and ∆S vib is the vibrational contributions.
S conf calculated via the formula: where k is the Boltzmann's constant, N is the Avogadro's number, x 1 and x 2 are the mole fractions of solid solution components assuming ideal mixing behavior. The ∆S vib was calculated via the formula: where S ss is the vibrational entropy of the given solid solution, S Phl , S Imp -vibrational entropies of its pure components: Phl-phlogopite; Imp-impurity (Ti or Cr) end-member. The values of S ss , S Phl and S Imp at given P-T conditions were calculated using the GULP software. The results of calculated configuration, vibrational and mixing entropies are available in Supplementary Electronic Table (Table E1). The simulations showed that S conf often provides the main contribution to the mixing entropy. The solid solutions formed via the mechanisms 2 VI (Mg 2+ ) = VI (Ti 4+ ) + VI (□) and VI (Mg 2+ ) + IV (Si 4+ ) = VI (Cr 3+ ) + IV (Al 3+ ), for which, based on the calculations, the vibrational contribution exceeds the configuration one, are exceptions.
The results of calculations of the mixing enthalpy ∆H mix for the considered isomorphic substitution mechanisms are shown in Fig. 2. Assuming ideal mixing behavior, the mixing enthalpy ∆H mix is a parabolic compositional function: where x 1 and x 2 are the mole fractions of solid solution components, Q-interaction parameter. Mixing enthalpy and entropy for pure components are equal zero. Calculations within the limits of infinite dilution were carried out in the approximation of a linear dependence of the interaction parameter Q, which may be represented by the additive compositional function (see supplementary information in Table S2): Compressibility parameters: a unit-cell volume vs. pressure, b bulk and shear modulus (Voight and Reuss bounds) vs. pressure, c compressional V P and shear V S sound velocity vs. pressure, of phlogopite modeled in this study compared with data from X-Ray dif-fraction studies (Pavese et al. 2003;Comodi et al. 2004) and with the results of the first principles simulation (Chheda et al. 2014) via two differet methods: LDA (local density approximation) and GGA (generalized gradient approximation) where Q 1 and Q 2 are the Margules parameters defined as follows: where E defect is the isolated defect energy of impurity incorporation into phlogopite (Imp in Phl) and incorporation of phlogopite into the impurity end-member; E str is the structural energy of corresponding solid solution end-member (Phl, phlogopite; Imp, impurity (Ti or Cr) end-member) per formula. The structural energies and energies of the formation of point defects were calculated using the GULP software without a pressure and temperature applied to the structure for the considered end-members of solid solutions. The results of simulations via three different methods for each considered isomorphic scheme are reported in Table 3. Structural energy is the energy of the formation of a crystal structure with a given charge state of atoms. Within the Mott-Littleton method, the defect energy (E defect ) is the energy difference between the perfect region 1 and 2 ( E p tot ) and the defective case ( E d tot ), rather than the individual contributions: , where x are the atomic cartesian coordinates and ξ are the cartesian displacements, respectively (Gale and Rohl 2003). The defect energy calculated for each mechanism in the selected supercells is an averaged value in the calculation for nonequivalent atomic configurations. The associate energy was normalized to 1 point defect.
As is evident, the data (Table 3) obtained in the 4 × 2 × 2 supercell are consistent with the results of simulation in the 6 × 3 × 3 cell and calculations via the Mott-Littleton model. The difference between the defect energies for all simulated end-members obtained by different methods does not exceed 0.4 eV. Such consistence allows us to consider the results obtained as quite reliable. Note that the defect energy is negative and minimal for the vacancy schemes of the isomorphic incorporation of minor elements into phlogopite.
The results obtained allowed us to evaluate the change in the Gibbs free energy (∆G mix = ∆H mix -T∆S mix ) depending on the compositions of the considered solid solutions in a given P-T range. The coexisting compositions of solid solutions were obtained at a minimum of ∆G mix by searching for the minimum of the fourth-degree polynomial function for each series, and the solvus curves were plotted for some of them (Fig. 3). The solvus lines of solid solutions, for which complete miscibility is fixed at a certain temperature, are the third-degree polynomial functions of a series of points at a given pressure. The critical temperatures were calculated as maxima of the solvus polynomial function.
Let us consider the conditions for the decomposition of each of the considered solid solutions.

Phlogopite-Ti-phlogopite solid solution
The highest isomorphic capacity in the series of the studied solid solutions of Ti-bearing phlogopite was gained via the mechanism VI (Mg 2+ ) + 2 IV (Si 4+ ) = VI (Ti 4+ ) + 2 IV (Al 3+ ) (Fig. 3a). The solvus of such solid solution is asymmetric (see supplementary information Table S3) and its maximum (critical point) is shifted towards Ti-phlogopite, which is consistent with the polarity rule of isomorphic substitutions (an ion with a high charge is more preferably incorporated Fig. 2 Comparison of the dependences of the mixing enthalpy on the composition of solid solutions for the considered schemes of impurity incorporation into mineral than an ion with a lower charge occupying the same crystallographic site upon heterovalent isomorphism). The vacancy mechanism VI (Mg 2+ ) + 2 IV (Al 3+ ) = VI (□) + 2 IV (Ti 4+ ) suggests limited miscibility (Fig. 3b), and the maximum isomorphic capacity of phlogopite detected at 7 GPa is almost identical to the capacity of its Ti end-member (22 mol % TiPhl). The isomorphic capacity in the KMg 3 AlSi 3 O 10 (OH) 2 -KMg 3 AlTi 3 O 10 (OH) 2 solid solution is significantly lower (Fig. 3c). The entry of Ti 4+ ions into phlogopite is extremely limited (microconcentrations) up to a temperature of 1373 K. The solvus curve is characterized by a slight asymmetry being shifted towards the Ti end-member, which is consistent with the polarity rule for isomorphism. The G-x sections demonstrate complete immiscibility over the entire range of the studied P-T parameters for a solid solution simulated via the vacancy mechanism 2 VI (Mg 2+ ) = VI (Ti 4+ ) + VI (□); therefore, this mechanism may provide minor concentrations of impurities only.

Phlogopite-Cr-phlogopite solid solution
Analysis of the G-x sections for a series of solid solutions with KMg 3 CrSi 3 O 10 (OH) 2 and K(Al,Cr,□) 2 AlSi 3 O 10 (OH) 2 end-members showed the complete miscibility even at room conditions (see supplementary information, Figs. S3 and S4). At the same time, it should be noted that a substitution like IV (Al 3+ ) = IV (Cr 3+ ) is most likely hypothetical and obviously unfavorable from the point of the crystal field theory (Cr 3+ ions are preferentially incorporated into the octahedral rather than tetrahedral sites) (McClure 1957). The crystal field stabilization energy (CFSE) is -1.2∆ o in the octahedral site and -0.8∆ t in the tetrahedral site for Cr 3+ ions; since ∆ t = 4/9∆ o , the CFSE in the tetrahedron is 0.35∆ o . Hence, the energy of octahedral site preference in the crystal field for the Cr 3+ is ∆E oct = 0.85. Thus, energy of the octahedral site preference by Cr ions provides an additional, and rather significant, energy effect, which cannot be taken into account within the applied model of classical potentials. Nevertheless, the values of the mixing enthalpy and entropy obtained for simulated solid solutions of phlogopite with K(Mg 2 ,Cr)Al 2 Si 2 O 10 (OH) 2 and K(Cr 2 ,□)AlSi 3 O 10 (OH) 2 end-members allowed us to define the limited miscibility when Cr 3+ ions enter the octahedron. Isomorphic substitution via the scheme VI (Mg 2+ ) + IV (Si 4+ ) = VI (Cr 3+ ) + IV (Al 3+ ) suggests limited miscibility over the entire tested temperature range (Fig. 3d). An increase in pressure leads to an expansion of the stability fields of both solid solutions, and this effect increases with temperature. The critical temperature is estimated as 1075 ± 17 K for the vacancy mechanism (Fig. 3e) in the pressure range of 3-5 GPa, (see supplementary information Table S3); however, at a pressure of 1 GPa, the critical temperature is much higher, 1370 K. Note that the solvus curve is relatively symmetrical and its maximum is close to the composition of 50 mol % K(Cr 2 ,□) AlSi 3 O 10 (OH) 2 (x(Cr) = 0.46-0.50 in the selected pressure range). The depression isomorphism rule characteristic of the simplest isovalent substitutions most often does not work for the studied solid solutions, which are characterized by limited solubility. This shows that the isomorphic capacity at various pressures in each case is controlled by a whole complex of parameters, such as the compressibility of various polyhedra, volume of vacancy regions, and their mutual arrangement.

Change of the structure geometry as a function of impurity concentration
Analysis of the local structure of the simulated solid solutions was carried out with account for changes in both individual and average interatomic distances and volumes of coordination polyhedra within the framework of the phenomenological theory (Urusov 1992 figure (Fig. 4). Figure 5 shows the change in the structure geometry depending on the composition of the solid solution upon incorporation of the considered impurities via the competing substitution mechanisms. Deviations from the additivity of structural parameters depending on the composition of the considered solid solutions are reported in supplementary information (Table S4).
Note that the dependences of structural parameters and cell volume on the composition of solid solution are close to linear upon Ti incorporation into phlogopite via to the non-vacancy substitution mechanisms IV (Si 4+ ) = IV (Ti 4+ ) and VI (Mg 2+ ) + 2 IV (Si 4+ ) = VI (Ti 4+ ) + 2 IV (Al 3+ ) (Fig. 5a). There are slight deviations of these mechanisms from additivity, which do not exceed 0.09 Å and 3.678 Å 3 for bond length and cell volume, respectively (see supplementary information Table S4). At the same time, change in the structural parameters for vacancy mechanisms has a pronounced non-linear character. The mechanism VI (Mg 2+ ) + 2 IV (Al 3+ ) = VI (□) + 2 IV (Ti 4+ ) is characterized by a negative deviation of a, c, and cell volume V, and a positive deviation of b from additivity. Such regularities are rather rare; however, according to (Urusov 1977), they may be one of the reasons for the observed increase in isomorphic capacity with increasing pressure.
The Vegard's and Retger's additivity rules are systematically violated upon the substitution like 2 VI (Mg 2+ ) = VI (Ti 4+ ) + VI (□) (see supplementary information Table S4). There is a regular increase in the unit-cell volume V and c parameter up to 40 mol % of Ti component. A further increase in the Ti component ceases to have a noticeable effect on the value of these parameters. This pattern may explain the detected immiscibility in the K(Mg,Ti,□)AlSi 3 O 10 (OH) 2 -KMg 3 AlSi 3 O 10 (OH) 2 series described in the previous section.
We also note that the vacancy mechanisms are characterized by a decrease in the density of the structure (Fig. 5). At the same time, the opposite effect is observed for non-vacancy mechanisms. According to the mechanism VI (Mg 2+ ) + 2 IV (Si 4+ ) = VI (Ti 4+ ) + 2 IV (Al 3+ ), the structure density increases with increasing concentration of the Ti component of the solid solution, which is naturally explained by an increase in the weight by the substitution of heavy Ti 4+ ions for light Mg 2+ . However, the density in the KMg 3 AlTi 3 O 10 (OH) 2 -KMg 3 AlSi 3 O 10 (OH) 2 series is constant, despite increase in the unit-cell parameters, which is compensated by the difference in interatomic distances: Ti-O bond length is larger than Si-O bond length.
Cr incorporation into phlogopite shows no obvious deviations from the additivity rules, which do not exceed 0.11 Å and 28.13 Å 3 for bond length and unit-cell volume, respectively (see supplementary information Table S4). A linear change in the density of Cr phlogopites is recorded for all mechanisms with the most significant increase for the vacancy mechanisms (Fig. 5b) which is explained by an increase in the weight of heavy substitutional ions of Cr 3+ for light Mg 2+ and Al 3+ , as well as for Ti substitutions.
Let us consider the patterns of the local structure of a solid solution in terms of changes in interatomic distances, volumes of coordination polyhedra, compliances of cation sites (Dollase 1980), and degree of their relaxation λ (Vegard 1928;Urusov 1992) depending on the impurity content. The values of λ for most polyhedra are close to the bond length alternation model, in which it is equal to 1 (Table 4, supplementary information Fig. S5a, b, e, f, Fig. S6a, b, e, f), which demonstrates the constancy of the lengths of individual bonds in the structures of analyzed solid solutions, regardless of the composition of the solid solution. This is explained by the key influence of the size of the common structural unit of an isomorphic mixture, which follows the isomorphism assistance rule (Urusov 1977). The

Comparison of incorporation of Cr and Ti in phlogopite
Comparison of the patterns of Cr 3+ and Ti 4+ incorporation into phlogopite by the competing substitution mechanisms shows that the formation of Cr-bearing solid solutions requires less energy. Although the radii of impurity ions in the octahedral site have very close values (rTi 4+ = 0.75 Å, rCr 3+ = 0.76 Å), incorporation of Ti 4+ ions distort the structure of the solid solution stronger than Cr 3+ ions compared to pure phlogopite (Fig. 5). This is due to the greater difference in charges and electronegativity of the Ti 4+ -Mg 2+ pair compared to the Cr 3+ -Mg 2+ pair. In addition, such distortion results in a stronger change in the structural energy in relation to the energy of pure phlogopite in the case of Ti 4+ incorporation (Table 3).

Comparison of isomorphic substitution schemes
Analysis of thermodynamic mixing properties, as well as structure geometry of the considered solid solutions, allowed us to identify the most energetically preferable schemes.
The most likely scheme of Ti incorporation into phlogopite is VI (Mg 2+ ) + 2 IV (Si 4+ ) = VI (Ti 4+ ) + 2 IV (Al 3+ ), which stimulates complete miscibility at relatively low temperatures in the entire pressure range of phlogopite stability. The simplest isovalent substitution (Si 4+ ) IV = (Ti 4+ ) IV is energetically inferior to that described above, despite the smaller deviations of structural parameters from additivity, and exhibits highly limited miscibility over the entire considered P-T range. The vacancy mechanism 2 VI (Mg 2+ ) = VI (Ti 4+ ) + VI (□) is the most preferable in the case of very low Ti concentrations only, which is explained by the minimum energies of the defect formation. Moreover, such substitution is unlikely in significant scale, as indicated by the facts that the additivity rules, as well as the Vegard's and Retger's rules, do not work, the Margules parameters have extremely high values, and, as a result, the immiscibility is detected. The formation of a vacancy in the octahedral site upon Ti incorporation into the tetrahedron most likely proceeds via the scheme VI (Mg 2+ ) + 2 IV (Al 3+ ) = VI (□) + 2 IV (Ti 4+ ). Such substitution is energetically more favorable in comparison with the scheme 2 VI (Mg 2+ ) = VI (Ti 4+ ) + VI (□) and allows the accumulation of higher Ti concentration (Table 3).
According to the results of the simulations, the schemes 3 VI (Mg 2+ ) = VI (Al 3+ ) + VI (Cr 3+ ) + VI (□) showing miscibility even at the room conditions and VI (Mg 2+ ) + IV (Si 4+ ) = VI (Cr 3+ ) + IV (Al 3+ ) with limited miscibility turn out to be the most energetically preferable. Anomalies with change in the structure geometry were not detected for these mechanisms and they are characterized by the low values of the Margules parameters and mixing enthalpy. The incorporation of chromium into phlogopite is less preferable via the substitution 3 VI (Mg 2+ ) = 2 VI (Cr 3+ ) + VI (□). Despite the absence of anomalies with change in the structure geometry, following the additivity rules, and detection of complete miscibility, this mechanism is the least favorable energetically (Table S2). As was mentioned above, the scheme IV (Al 3+ ) = IV (Cr 3+ ) is obviously unfavorable in accordance with the crystal field theory, despite the extremely low values of the mixing enthalpy (Fig. 2).

Some reasoning for the schemes of isomorphic substitution under the natural conditions
The quantitative estimates of the isomorphic capacity of phlogopite in terms of the content of Ti and Cr impurities in comparison with their calculated concentrations in the studied end-members of solid solutions, as well as the highest Ti and Cr concentrations recorded in phlogopites from inclusions in natural diamonds, other natural mineral assemblages, and those synthesized in different experimental systems are reported in Table 5. As is evident, the results of atomistic simulation do not contradict the regularities in impurity accumulation detected in natural samples and in run products (Table 5). The maximum Ti concentrations calculated for the most energetically preferable stoichiometric mechanism VI (Mg 2+ ) + 2 IV (Si 4+ ) = VI (Ti 4+ ) + 2 IV (Al 3+ ) exceed those in natural assemblages (Table 5). At the same time, the formation of a vacancy in the octahedral site is often reported in the literature. In the absence of X-ray diffraction data, information on the content of Fe 3+ and OH, as well as unconfirmed Ti-oxy substitution (Mg,Fe 2+ ) + 2OH − = Ti 4+ + 2O 2− (Ventruti et al. 2020;Cruciani and Zanazzi 1994;Schingaro et al. 2005;Mesto et al. 2006), the 2 VI (Mg 2+ ) = VI (Ti 4+ ) + VI (□) substitution is discussed in many studies (e.g., (Arima and Edgar 1981;Tronnes et al. 1985;Abrecht and Hewitt 1988;Thu et al. 2016)); based on the results of our study, we suggest only minor Ti incorporation via this scheme. In addition, the TiO 2 content of > 7.81 wt % (> 50 mol % K(Mg,Ti,□) AlSi 3 O 10 (OH) 2 ) should result in the formation of transitional di-/trioctahedral mica varieties, which were not recorded in the literature (Koval et al. 1988;Sobolev et al. 2009).
At the same time, Ti incorporation via the scheme VI (Mg 2+ ) + 2 IV (Al 3+ ) = VI (□) + 2 IV (Ti 4+ ) may result in accumulation of 15.58 wt % TiO 2 in the composition of the K 2 (Mg 5 ,□)Ti 2 Si 6 O 20 (OH) 4 component (Table 5) without significant structural changes and with retention of three-layer packages of divalent cations. It is important that such Al substitution for Ti leads to the complete absence of Al in the tetrahedron, which is unlikely under the natural conditions. Thus, comparing the substitutions 2 VI (Mg 2+ ) = VI (Ti 4+ ) + VI (□) and VI (Mg 2+ ) + 2 IV (Al 3+ ) = VI (□) + 2 IV (Ti 4+ ), we may suggest that the accumulation of Ti will be more energetically favorable in the tetrahedron, while the high Al concentrations may be gained in the octahedral site via one of the schemes: VI (R 2+ ) + IV (Si 4+ ) = VI (Al 3+ ) + IV (Al 3+ ) or 3 VI (R 2+ ) = 2 VI (Al 3+ ) + VI (□). Such multistage substitution can compete with Ti-oxy substitution, and the choice of a certain mechanism will be determined by the conditions of mineral crystallization.
As was shown above, the classical potential may not cover the essential crystal chemistry to make more quantitative predictions. However, the maximum Cr concentrations in phlogopite obtained from the data of atomistic modeling were higher than those in the natural samples and in the previous experimental studies. We should specially note the high Cr 2 O 3 content (up to 6.54 wt %) reported for trioctahedral