Theoretical OH stretching vibrations in dravite

Density functional theory is used to investigate the vibrational stretching properties of OH groups in ideal and chemically disordered dravite models. Different schemes of cationic occupancy are considered, including the occurrence of vacancies at the X site and Mg–Al inversion between the Y and Z sites. The harmonic coupling between different OH groups is found to be smaller than 1 cm⁻¹, indicating that the OH stretching dynamic in dravite can be described by considering a collection of nearly independent single OH oscillators. Their harmonic stretching frequency is linearly correlated with the corresponding OH bond length and most of the bands observed in the experimental vibrational spectra can be interpreted as consequences of the cationic occupancy of the sites coordinated to the OH group. The ^V(OH) and ^W(OH) stretching frequencies are affected by the Mg–Al inversion and by the presence of vacancies at the X site. In this last case, the frequencies depend on the isolated or more concentrated character of the X vacancy distribution along the c axis. Based on theoretical stretching frequencies, new interpretations are proposed for some of the bands experimentally observed in synthetic samples of dravite and magnesio-foitite.

1 Introduction

Tourmalines form a super-group of complex hydrous borosilicates with the following general formula: XY₃Z₆(T₆O₁₈)(BO₃)₃V₃W (Henry et al., 2011). According to the chemical systematics reported by Henry et al. (2011), the most common ions (or vacancy) at each site are X = Na⁺, Ca²⁺, K⁺, and vacancy; Y = Fe²⁺, Mg²⁺, Mn²⁺, Al³⁺, Li⁺, Fe³⁺, and Cr³⁺; Z = Al³⁺, Fe³⁺, Mg²⁺, and Cr³⁺; T = Si⁴⁺, Al³⁺, and B³⁺; B = B³⁺; V = OH⁻ and O²⁻; and W = OH⁻, F⁻, and O²⁻. Tourmaline samples thus correspond to solid solutions of a multitude of end-members (Henry et al., 2011). Even in samples close to end-member compositions, disordered cation distributions occur, as revealed by structure refinements based on X-ray diffraction data (e.g., Hawthorne et al., 1993; Bosi et al., 2016). Due to its crystal–chemical versatility, tourmaline is an important indicator of petrogenetic conditions (Van Hinsberg et al., 2011) and has promising environmental applications (Liang et al., 2021).

Information about the local structure of tourmaline-group minerals can be obtained from the spectroscopic study of the OH stretching modes which lead to characteristic bands in the Fourier transform infrared (FTIR) or Raman spectra (e.g., Gonzalez-Carrenño et al., 1988; Fantini et al., 2014; Bosi et al., 2016; Berryman et al., 2016; Watenphul et al., 2016). Compared with lattice modes, the high frequency and weak dispersion of OH stretching modes result in a localized character typical of hard modes (e.g., Salje and Bismayer, 1997), which facilitates their use as probes of short-range order in minerals. Accordingly, the OH stretching spectra of chemically complex natural tourmaline samples could provide relevant information to identify specific cationic configurations and to relate them to the corresponding end-members (e.g., Watenphul et al., 2016). Two categories of structural OH groups can be distinguished in tourmaline (e.g., Gatta et al., 2014). The ^W(OH) shares its oxygen atom (O1) with the three octahedral Y sites and points to the X site along the [001] direction (Fig. 1). The 3 times more abundant ^V(OH) groups are coordinated to two Z and one Y octahedral sites with a weak H bonding interaction with the distant O5 oxygen bridging two tetrahedral sites (Gatta et al., 2014). Due to the sensitivity of OH stretching frequencies to their atomic-scale environment, a diversity of vibrational spectra is observed in the typical range of OH stretching frequencies (3000–3800 cm⁻¹) reflecting the variability in tourmaline chemical composition (e.g., Watenphul et al., 2016). For example, dravite, with nominal chemical composition ^XNa^YMg₃^ZAl₆(^TSi₆O₁₈)(BO₃)${}_{\mathrm{3}}^{\mathrm{V}}$(OH)₃^W(OH), displays a series of weak narrow bands between ∼ 3600 and 3800 cm⁻¹ and a broad and asymmetric stronger band spread between 3400 and ∼ 3600 cm⁻¹ (e.g., Fantini et al., 2014; Bosi et al., 2016; Berryman et al., 2016; Watenphul et al., 2016). Whereas the former bands have been related to the ^W(OH) site and the latter to ^V(OH) site, their assignment to specific atomic-scale environments is still debated, and two different models have been proposed to interpret the ^V(OH) stretching spectrum of tourmaline.

Figure 1 Dravite cluster displaying the main sites exerting a control on the OH stretching frequencies: lateral view (a) and projection along the [001] direction (b). The three symmetrically equivalent ^V(OH) groups are coordinated to one Y and two Z sites, whereas the ^W(OH) group is coordinated to three Y sites. The purple grey triangles at the base of the Z and Y sites correspond to BO₃ units. O and H atoms belonging to OH groups are represented as red and white spheres, respectively (radii are only illustrative).

In the “traditional” model (e.g., Gonzalez-Carrenño et al., 1988; Fantini et al. 2014; Berryman et al., 2016; Kutzschbach et al., 2017; Pieczka et al., 2020), the ^V(OH) stretching frequency is dominantly determined by the nature of cations coordinated to the OH group. This model is based on the observation that for a constant H bonding environment, the OH stretching frequency is primarily affected by the characteristics of the chemical bond that the oxygen atom shares with neighboring cations, through bond valence requirements, and to a lesser extent by the perturbations arising from substitutions at more distant sites (e.g., Hawthorne, 2016; Bronzova et al., 2019). A different interpretation of the ^V(OH) bands has however been proposed by Watenphul et al. (2016) and further used by Bosi et al. (2016). In this interpretation, it is suggested that the ^V(OH) stretching mode frequency is determined by the chemistry of the cluster formed by the YZZ sites which are symmetrically related by the 3-fold rotation axis in the rhombohedral unit cell (Fig. 1) and not by the chemistry of a sole YZZ species. This assumption implies a collective response of the three ^V(OH) groups of the cluster to chemical disorder (Hawthorne, 2016). To illustrate the difference between the two models (Bosi et al., 2016), a ^YMg^ZAl^ZAl-^YMg^ZAl^ZAl-^YAl^ZAl^ZAl cluster is expected to lead to two different ^V(OH) stretching bands in the traditional model (one for the ^YMg^ZAl^ZAl environment and one for the ^YAl^ZAl^ZAl environment), whereas the “YZZ-triplet” model predicts a single ^V(OH) band for this cluster. In this latter case, experimental observations suggest that the frequency related to a 2^YMg^ZAl^ZAl-^YAl^ZAl^ZAl cluster in dravite would be smaller than that related to a 3^YMg^ZAl^ZAl cluster by ∼ 40 cm⁻¹ (Watenphul et al., 2016).

According to Watenphul et al. (2016), the theoretical basis proposed for the YZZ-triplet model relates to the vibrational properties of mixed crystals, as reviewed by Chang and Mitra (1971). A property shared by the mixed crystals investigated by Chang and Mitra (1971) is their lattice-type dynamic, which in its simplest expression is analogous to that of a linear chain of point masses connected to each other by springs. Depending on the interplay of atomic masses and spring constants, these systems can present a continuous shift of long-wavelength phonon frequencies as a function of their chemical composition (one-mode behavior; e.g., Na_xK_1−xCl) or two sets of phonon frequencies close to those of the end-members (two-mode behavior; e.g., CdS_xSe_1−x). However, the influence of chemical disorder on the vibrational properties of crystals with molecular units, such as tourmaline-group minerals, could differ from that observed on the crystals with a lattice-type dynamic investigated by Chang and Mitra (1971). When the frequency of the internal vibrations is significantly higher than that of the lattice phonon branches, intermolecular interactions act as a perturbation of the single molecule properties (e.g., Farmer, 1974; Califano et al., 1981). Collective properties are then related to the harmonic coupling between single oscillators, as well illustrated by the OH stretching phonons in, for example, brucite (Mitev et al., 2010). In this picture, the local or collective character of the response to chemical disorder depends on the respective magnitude of the intermolecular coupling and chemically induced local perturbation. In the weak coupling limit, chemical disorder can differently influence the vibration of single units, whereas a more collective response is expected when the vibrational coupling is stronger than the perturbation.

In the present study, we theoretically investigate at a quantum chemical level the influence of chemical disorder on the OH stretching modes of dravite. The modeling is based on density functional theory (DFT) which is a level of theory accurate enough to assign the OH stretching bands of hydrous minerals to specific vibrational modes involving the coupled motions of OH groups (e.g., Balan et al., 2001, 2002, 2005). Among the various end-members of tourmaline, we selected dravite because of the well-tested performances of DFT in the modeling of Al and Mg hydrous minerals (e.g., Balan et al., 2001, 2002, 2005; Mitev et al., 2010; Blanchard et al., 2018). We first confirm that the OH stretching modes in ideal dravite are almost fully decoupled from the other vibrational modes, and we determine the magnitude of the harmonic vibrational coupling between single OH oscillators. We then determine the theoretical stretching frequencies of OH groups in models with cationic occupancies that vary from ideal dravite, showing that their response to a chemical perturbation is consistent with the “traditional” interpretation proposed for the ^V(OH) stretching bands. Based on the theoretical frequencies, revised interpretations are proposed for some of the ^W(OH) and ^V(OH) bands experimentally observed in synthetic samples.

2 Methods

Relevant properties were obtained using the density functional theory (DFT) and the generalized gradient approximation (GGA) to the exchange–correlation functional (XC functional) as proposed by Perdew, Burke and Ernzerhof (PBE) (Perdew et al., 1996). The modeling was done using a plane-wave and pseudo-potential scheme, as implemented in the PWscf and PHonon codes from the Quantum ESPRESSO package (Giannozzi et al., 2009; http://www.quantum-espresso.org, last access: 3 April 2022). The ionic cores were described by optimized norm-conserving pseudo-potentials from the SG15 ONCV library (Hamman, 2013; Schlipf and Gigy, 2015). The wave functions and the charge density were expanded in plane waves with 80 and 480 Ry cutoffs, respectively. The discrete sampling of the Brillouin zone for electronic integration was performed using a $\mathrm{3}\times \mathrm{3}\times \mathrm{3}$ k-point grid for the primitive rhombohedral cell. It was restricted to a single k point with coordinates (0.0, 0.0, −1.1063832) in 2$\mathit{\pi }/a$ units for the hexagonal cells and to the Γ point for the $\mathrm{1}\times \mathrm{1}\times \mathrm{2}$ hexagonal supercells. Relaxation of atomic positions was performed until the residual forces on atoms were less than 10⁻⁴ Ry a.u.⁻¹.

The relaxed R-centered hexagonal cell parameters (a=16.11 Å, c= 7.28 Å) of ideal dravite (space group R3m) overestimate by ∼ 1 % the corresponding experimental values (e.g., a= 15.9297 Å, c= 7.2064 Å; Bosi et al., 2016). This overestimation is consistent with previous studies of other minerals performed at the same theoretical level (e.g., Balan et al., 2019). Models corresponding to chemically disordered dravite have been built by varying the site occupancies of a single 3YZZ-X cluster in the hexagonal cell of dravite. The structure of the chemically disordered models (Supplement) was then obtained by relaxing all the atomic coordinates. The cell parameters were kept equal to those of ideal dravite to avoid unwanted distortions of the periodic lattice due to the locally induced symmetry lowering. Due to the moderate value of the c cell parameter of the hexagonal cell (7.28 Å), some interaction can however occur between the periodic images of the modified sites along the c axis. Accordingly, the most pertinent configurations of disordered dravite have been also investigated using $\mathrm{1}\times \mathrm{1}\times \mathrm{2}$ hexagonal supercells, which contain 323 or 324 atoms depending on the model.

For the models with heterovalent substitutions, the macroscopic neutrality of the periodic system was ensured by a compensating homogeneous electrostatic background spread over the cell. This electrostatic correction avoids the divergence of the total energy but is not expected to affect the microscopic properties of ionic systems (e.g., Leslie and Gillan, 1985).

Harmonic vibrational properties (atomic displacements and frequencies of the transverse optical modes at the Brillouin zone center) were calculated from the second-order derivatives of the total energy with respect to atomic displacements using the linear response theory (Baroni et al., 2001) as implemented in the PHonon code (Giannozzi et al., 2009). The eigendisplacements of OH stretching modes are reported in the Supplement. Note that coordinates of eigendisplacement vectors correspond to the coordinates of the dynamical matrix eigenvectors divided by the square root of the relevant atomic mass.

3 Results

3.1 Theoretical properties of OH stretching modes in ideal dravite

To ascertain that the OH stretching modes in dravite are decoupled from other vibrational modes occurring at lower frequency, the full dynamical matrix of the primitive rhombohedral cell was diagonalized, leading to four OH stretching modes among a total of 159 vibrational modes. These four modes correspond to the ^W(OH) stretching mode at 3856 cm⁻¹ and to one and two ^V(OH) modes belonging to the A₁ and E irreducible representations, respectively, as determined by the group theory analysis of Watenphul et al. (2016). The A₁ mode at 3653.4 cm⁻¹ corresponds to the in-phase displacement of the three symmetrically equivalent ^V(OH) groups and the doubly degenerated E modes at 3653.5 cm⁻¹ to their out-of-phase displacements. The infrared or Raman intensities are carried by the A₁ mode because the three ^V(OH) groups are almost parallel to each other. Then, the diagonalization of a partial dynamical matrix was performed by only considering the displacements of the oxygen and hydrogen atoms belonging to the OH groups, setting to zero all other coefficients. The OH stretching frequencies obtained using the partial matrix differ by less than 0.15 cm⁻¹ from those resulting from the full computation (Supplement Table S1), which attests to the quality of the approximation. Accordingly, the OH stretching modes form a subset of vibrational modes decoupled from all the other modes. This also implies that the mass of the other atoms does not play any significant role in the OH stretching dynamic. This property is specific to the stretching modes, which are well separated from the other modes occurring at lower frequency (e.g., Balan et al., 2005, 2017). The other OH vibrations, corresponding to two libration and three hindered translation displacements per OH group, cannot be determined in such a straightforward manner because they can couple to other atomic motions due to their lower frequency.

The OH stretching vibrations, being decoupled from the other degrees of freedom, can be described by solely considering the OH stretch coordinates (e.g., Califano et al., 1981; Mitev et al., 2010). The zone-center transverse optical stretching modes then relate to a harmonic cell Hamiltonian with diagonal terms depending on the sum of single bond contributions and out-of-diagonal terms accounting for the bilinear coupling between OH groups. In such local mode description, the splitting of the A₁ and E stretching modes of three symmetrically equivalent ^V(OH) groups is determined by a single parameter λ measuring their harmonic coupling: λ= (ω(A₁)-$\mathit{\omega }\left(E\right))/\mathrm{3}$ (see Supplement). If the dominant interactions occur within the 3(YZZ) cluster, the local mode description is formally equivalent to that of the potential coupling of three equivalent X–H groups in a C_3V molecule, where X is a heavier atom (e.g., Halonen and Child, 1983). Accordingly, the very small splitting (∼ 0.1 cm⁻¹) observed between the A₁ and E ^V(OH) stretching modes denotes the weakness of the coupling interaction between the ^V(OH) oscillators in dravite. Inspection of the eigendisplacement vectors also indicates that the ^W(OH) stretching is not significantly coupled to the ^V(OH) stretching. It is thus arguable that the OH stretching dynamic in dravite is better described by a collection of nearly independent oscillators than by collective atomic motions typical of lattice modes.

The hexagonal cell offers the possibility to assess the vibrational interaction between OH groups belonging to different ^W(OH)-^Y(Mg)₃–^V(OH)₃–^Z(Al)₆ clusters. For that purpose, the full dynamical matrix of the hexagonal cell was diagonalized, leading to 12 OH stretching modes among a total of 483 vibrational modes. These 12 modes correspond to 3 ^W(OH) stretching modes and to 3 A₁ and 6 E ^V(OH) modes. Vibrational interactions lead to a spreading over ∼ 0.04 and ∼ 1.4 cm⁻¹ of the three ^W(OH) and nine ^V(OH) frequencies of the full system, respectively (Table S2). Using a partial dynamical matrix reduced to the four OH groups belonging to a single cluster, the ^W(OH) stretching mode is calculated at 3856.2 cm⁻¹, whereas the A₁ and E ^V(OH) modes are calculated at 3653.2 and 3652.3 cm⁻¹ (Table S2). These frequencies differ by less than 0.2 cm⁻¹ from the mean frequency of the modes with corresponding symmetry in the full dynamical system. They are also within 1.5 cm⁻¹ from those obtained on the rhombohedral cell using a comparatively denser k-point grid.

3.2 Influence of cationic ordering and X site vacancies on OH groups in dravite

The considered models of non-ideal dravite (Table 1) include vacant X site, Mg for Al substitution at Z site, Al for Mg substitution at Y site and Al–Mg inversion involving Y and Z sites. Models combining a vacant X site with one or two Al substitutions in nearby Y sites or Al–Mg inversion have also been considered. In these cases, two configurations occur in the $\mathrm{1}\times \mathrm{1}\times \mathrm{2}$ hexagonal supercell depending on the relative position of the vacant X site and modified Y and Z sites along the c axis. All the models have a local charge disequilibrium smaller than or equal to one (Table 1). Note that the local charge imbalance induced by heterovalent substitutions or X vacancies (Table 1) is compensated by a countercharge homogeneously spread over the system. Although not explicitly accounting for coupled substitutions, this approach makes it possible to disentangle the influence on the OH stretching modes of individual contributions to the chemical disorder.

Table 1 Local charge imbalance, theoretical OH stretching frequencies (ω_OH) and theoretical O–H distances (d(OH)) in dravite models obtained using the hexagonal cell or a larger $\mathrm{1}\times \mathrm{1}\times \mathrm{2}$ hexagonal supercell.

n: number of symmetrically equivalent OH groups, leading to n modes with frequencies within 1.3 cm⁻¹.

As supported by the weakness of the vibrational coupling between single OH oscillators, the OH stretching frequencies in non-ideal dravite models are obtained by restricting the calculation to the dynamical matrix coefficients involving the O and H atoms belonging to the four OH groups of the chemically disordered cluster. In the vacancy-bearing models using a $\mathrm{1}\times \mathrm{1}\times \mathrm{2}$ supercell, only the stretching modes of ^W(OH) or ^V(OH) groups of the perturbed cluster and occurring as nearest neighbors of the vacant site are reported in Table 1. Based on the analysis of ideal dravite properties, it can be assumed that the OH stretching frequencies are within ∼ 1 cm⁻¹ from those obtained using the full dynamical matrix. This approach is similar to the partial Hessian techniques which are commonly used to analyze the vibrational properties and contribution to zero point energy of selected parts of extended systems (e.g., Li and Jensen, 2002; Svelle et al., 2008; Ghysels et al., 2010; De Moor et al., 2011; Tao et al., 2018; Pascale et al., 2018; Kraka et al., 2020). Their accuracy has been demonstrated for highly localized vibrations, such as OH stretching modes (e.g., Li and Jensen, 2002; Ghysels et al., 2010).

Influence of cationic substitutions at Y and Z sites

Cationic substitutions affecting the Y or Z sites and vacancies at the X site (Fig. 1) modify the geometry of the cluster, and the specific contribution of an OH group to a given mode can be identified by inspecting the related eigendisplacement vectors (Supplement). When the 3-fold symmetry is preserved (Fig. 1), the ^V(OH) stretching modes correspond to three modes at similar frequencies belonging to the A₁ or E representations as in ideal dravite (Table 1). When the symmetry is reduced to a mirror symmetry, one mode corresponds to the OH group located on the mirror plane. The two other OH groups lead to two modes at similar frequencies corresponding to their in-phase or out-of-phase motion (Table 1). Among these two modes, the spectroscopic activity is concentrated in the mode with the in-phase motion of OH groups. For a lower symmetry, the three OH groups are no more equivalent, and three modes are observed at different frequencies (Table 1). It is also noteworthy that the presence of chemical disorder has an effect on the OH stretching frequencies whose magnitude can significantly exceed the weak coupling between single OH oscillators (Table 1). Therefore, it is not possible to establish a symmetry-based filiation between the vibrational patterns in most of the disordered models and those of ideal dravite.

Considering the results obtained for the large $\mathrm{1}\times \mathrm{1}\times \mathrm{2}$ hexagonal supercell (Table 1), the theoretical stretching frequency of the ^W(OH) group decreases to 3808 cm⁻¹ when Al is substituted for Mg at the Y site. The stretching of the ^V(OH) coordinated to three Al cations is calculated at 3573 cm⁻¹. The two other ^V(OH) groups coordinated to one ^YMg and two ^ZAl cations (as in ideal dravite) are symmetrically equivalent by a mirror symmetry. The two nearly degenerated modes, corresponding to their in-phase or out-of-phase motions, are observed at ∼ 3650 cm⁻¹. This frequency is close to that determined for the similar environment in ideal dravite (3653 cm⁻¹). In the case of the substitution of Mg for Al at the Z site, the cluster looses its 3-fold axial and mirror symmetries. The theoretical stretching frequency of the ^W(OH) group is calculated at 3860 cm⁻¹, close to the frequency of ideal dravite (3856 cm⁻¹) in which the ^W(OH) group has a similar local environment. This highlights the fact that the ^W(OH) site is not significantly influenced by the Z site occupancy. The stretching frequency of the ^V(OH) group coordinated to ^YMg^ZMg^ZAl cations increases to 3692 cm⁻¹. The two remaining ^V(OH) groups associated with ^YMg^ZAl^ZAl configurations are no longer symmetrically equivalent (Fig. 1) and lead to modes at 3655 and 3657 cm⁻¹, which is still close to the ^V(OH) stretching frequency in ideal dravite (3653 cm⁻¹). These results attest to a dominant influence of the occupancy of the nearest-neighbor cationic sites and to a smaller influence (∼ 3 cm⁻¹) of the occupancy of the next-nearest-neighbor sites on the OH stretching frequency.

Influence of Al–Mg inversion

A series of models with different configurations corresponding to Al–Mg disorder over the Y and Z sites (Fig. 1) has been investigated in restricting the calculation to the simple hexagonal cell. Among these configurations (Supplement Table S2), only the one with a Mg–Al inversion at neighboring Y and Z sites has been found to be energetically competitive with the ideal dravite configuration, with an increase in the total energy limited to 13.4 kJ mol⁻¹ (Table S2). The other configurations are significantly less stable (by more than 40 kJ mol⁻¹; Table S2) and therefore less likely to occur at thermodynamic equilibrium, even at temperatures as high as 800 ^∘C. This temperature is beyond the low-pressure stability domain of tourmaline (Van Hinsberg et al., 2011) and corresponds to that of the dravite thermal treatment of Bosi et al. (2016). It is higher than the synthesis temperature (600 ^∘C) of Berryman et al. (2016).

The properties of the most stable Mg–Al inversion in which ^YAl and ^ZMg occur at neighboring sites have then been computed using the larger $\mathrm{1}\times \mathrm{1}\times \mathrm{2}$ hexagonal supercell. The ^W(OH) frequency (3811 cm⁻¹; Table 1) is very close to that determined for the Al substitution at Y site (3808 cm⁻¹; Table 1). The two ^V(OH) groups coordinated to two ^ZAl and one ^YMg lead to frequencies of 3652 and 3654 cm⁻¹, still within 2 cm⁻¹ from those in ideal dravite (3653 cm⁻¹). In contrast, the frequency of the ^V(OH) group in a ^YAl^ZMg^ZAl environment is calculated at 3611 cm⁻¹, significantly lower than that determined for the ^YMg^ZAl₂ environment (3653 cm⁻¹). The lower frequency in the inversion model is correlated with a longer OH bond (0.975 vs. 0.972 Å in ideal dravite; Fig. 2), which likely arises from shorter Al–OH distances (1.95 and 1.97 Å) in the inverted configuration than in ideal dravite (1.99 Å). Accordingly, iso-chemical configurations corresponding to various distributions of Al and Mg cations over the Y and Z sites cannot be assumed to be energetically degenerated.

Influence of X site vacancies

In the dravite model with a vacant X site (^X□), the vibrational pattern is similar to that of ideal dravite, but OH stretching frequencies are shifted to lower values (Table 1). The ^W(OH) and ^V(OH) stretching modes are computed at 3754 and 3628 cm⁻¹, respectively, corresponding to a significant decrease in the stretching frequency of the OH groups in the nearest environment of the vacant site. The presence of a vacancy combined with Al for Mg substitutions at Y sites or Al–Mg inversion also reduces the stretching frequencies with respect to the similar environments in Na-bearing systems (Table 1). For the ^W(OH) groups, the frequency decrease ranges from approximately −80 to −100 cm⁻¹.

Finally, the interactions between chemically disordered clusters can be assessed by comparing the stretching frequencies obtained for the $\mathrm{1}\times \mathrm{1}\times \mathrm{2}$ hexagonal supercell with those obtained for the same configurations but using the smaller simple hexagonal cell (Fig. 3). The differences are within ± 15 cm⁻¹ for the majority of the OH groups. However, larger variations in the stretching frequency, up to 27 cm⁻¹, are observed for both the ^W(OH) and ^V(OH) groups in the models combining Al substitutions at Y sites and vacancies at X sites, as well as for the ^V(OH) group in ^YAl^ZAl₂-^XNa and ^YAl^ZAl^ZMg-^X□ configurations. These more significant frequency variations are related to a change in the OH bond length (reaching −0.0018 Å), which indicates that they originate from static interactions between modified clusters in adjacent cells. These observations are particularly significant for the vacancy-bearing models because they suggest that the ordering scheme of X vacancies along the c axis could have a significant effect on the stretching frequencies of the ^W(OH) and ^V(OH) groups.

4 Discussion

4.1 General properties of OH stretching modes in dravite

The marginal coupling between single OH stretching vibrations in dravite implies that their frequency is essentially determined by the force constant of the OH bond. Consistently, the OH stretching frequencies in dravite are linearly correlated with the corresponding bond lengths (Fig. 2). The slope of the correlation (−15447 cm⁻¹ Å⁻¹; Fig. 2) is consistent with that theoretically determined for an OH⁻ ion submitted to an external electric field (−16000 cm⁻¹ Å⁻¹) by Hermansson (1993), as well as with those theoretically determined, for example, for OH defects in quartz (−17270 cm⁻¹ Å⁻¹; Jollands et al., 2020) and corundum (−18878 cm⁻¹ Å⁻¹; Balan, 2020). Through the modifications of the OH bond length, the stretching frequencies depend on the chemical nature of atoms in their local environment, with a lesser influence of the occupancy of sites beyond the first coordination sphere. However, it is noteworthy that the expected relation between the OH stretching frequency and the sum of nearby cationic charges is not strictly followed. This is particularly the case for the Mg–Al inversion, in which the stretching frequency of the ^V(OH) group in a ^YAl^ZAl^ZMg-^XNa environment (3611 cm⁻¹) is significantly lower than that corresponding to a ^YMg^ZAl^ZAl-^XNa environment (3653 cm⁻¹). Similar frequency variations are observed for the equivalent configurations with X vacancies (Table 1).

Figure 2 Relation between theoretical OH stretching frequency (ω_OH) and OH bond length (d(OH)) in dravite models ($\mathrm{1}\times \mathrm{1}\times \mathrm{2}$ supercell and hexagonal cell; Table 1).

Figure 3 Difference between the OH stretching frequencies computed using a simple hexagonal cell and a larger $\mathrm{1}\times \mathrm{1}\times \mathrm{2}$ supercell ($\mathrm{\Delta }{\mathit{\omega }}_{\mathrm{OH}}={\mathit{\omega }}_{\mathrm{OH}}\left(\text{simple}\right)-{\mathit{\omega }}_{\mathrm{OH}}\left(\text{supercell}\right)$) reported as a function of the frequencies computed for the $\mathrm{1}\times \mathrm{1}\times \mathrm{2}$ supercell (Table 1). Squares: ^V(OH) groups; circles: ^W(OH) groups. Full symbols: models with a vacant X site; open symbols: models with full Na occupation of X sites.

Modifications of Y or Z site occupancies, which shift the stretching frequency of the neighboring OH group by typically more than 20 cm⁻¹, are significantly stronger than the vibrational coupling between single OH groups (<1 cm⁻¹). Thus, these modifications cannot be considered as a weak perturbation of the vibrational properties of the ideal system. This shift however still corresponds to a weak perturbation of the internal stretching of single OH groups (>3400 cm⁻¹), supporting the use of OH stretching modes as a local probe of chemical disorder. This behavior differs from that of lattice excitations in mixed crystals, which present a collective response to local mass or force constant variations (Chang and Mitra, 1971). In particular, it does not appear necessary to consider a collective response of all three ^V(OH) groups belonging to a given 3(YZZ) cluster in the interpretation of the ^V(OH) stretching bands. The present theoretical results support a traditional interpretation of observed ^V(OH) frequencies based on the cationic occupancy of the sites coordinated to the OH group (e.g., Berryman et al., 2016; Kutzschbach et al., 2017; Piecza et al., 2020).

Although the OH groups of a specific cluster are not vibrationally coupled to each other or to those of other clusters, the local modification of cationic occupancies can affect the properties of more distant OH groups through electrostatic interactions or propagation of structural distortions. The influence of chemical disorder on OH groups in the medium range of the modified sites can be assessed by examining the bond length of the ^W(OH) and ^V(OH) groups associated with the ideal ^YMg₃-^XNa and ^YMg^ZAl₂-^XNa environments, respectively, in the chemically disordered models. In this case, departure from the bond lengths determined in ideal dravite can be ascribed to the chemical modification of sites beyond their first coordination sphere. Using the correlation between stretching frequency and bond length, it then possible to assess the related distribution of stretching frequencies (Fig. 4). For the $\mathrm{1}\times \mathrm{1}\times \mathrm{2}$ supercell models (Table 1), most of the stretching frequencies of the ^W(OH) in ^YMg₃-^XNa environment occur within a ∼ 15 cm⁻¹ range. A similar range is observed for the ^V(OH) groups in ^YMg^ZAl₂-^XNa environment. Excluding the specific interactions involving the ordering of X site vacancies, these ranges are consistent with the magnitude of interactions inferred for the modified clusters by comparing the stretching frequencies in the $\mathrm{1}\times \mathrm{1}\times \mathrm{2}$ supercell and the smaller hexagonal cell (Fig. 3).

Figure 4 Distribution of the theoretical OH stretching frequencies determined from the length of OH bonds in ideal local configurations in chemically disordered models ($\mathrm{1}\times \mathrm{1}\times \mathrm{2}$ supercells). The theoretical frequencies of ideal dravite are reported. The observed distributions attest to the static interaction between the sites with modified occupancies and their medium-range environment.

The distribution of frequencies related to the occupancy of sites beyond the nearest-neighbor environment of OH groups can contribute to the broadening of the OH stretching signals for ∼ 15 cm⁻¹. This value has been assessed from a limited number of cationic configurations, and it is likely that the higher number of configurations in real samples further increases the broadening. In addition, extended chemical variations as observed in solid solutions can lead to variations in the cell geometry (e.g., Bosi and Lucchesi, 2004). Local strain contributions not accounted for in the present modeling (performed at constant volume) can also indirectly influence the ^V(OH) and ^W(OH) stretching frequencies, contributing to their distribution in chemically complex samples. In comparison, the homogeneous broadening of the OH stretching bands is expected to be of smaller magnitude. This broadening is due to the finite lifetime of the vibrational excitations. It arises from anharmonic phonon–phonon interactions related to coupling terms involving third- and higher-order derivatives of the total energy with respect to atomic displacements. In general, the high frequency of OH stretching modes limits energy decay processes and favors dephasing processes related to phonon–phonon scattering (e.g., Bonn et al., 1996). These anharmonic effects usually induce a moderate broadening of the OH stretching bands as a function of temperature (e.g., Zhang et al., 2007; Balan et al., 2010). The weak temperature dependence of the OH stretching line shape in the Raman and infrared spectra of tourmaline reported by Zhao et al. (2012) and Watenphul et al. (2017) confirms that homogeneous broadening mechanisms should not play a dominant role in the observed linewidth. Accordingly, inhomogeneous contributions related to disorder have a most likely dominant control on the width of the FTIR or Raman OH stretching bands in tourmaline.

4.2 Interpretation of the dravite vibrational spectra in the OH stretching range

Turning to the interpretation of experimentally observed bands, we mainly focus on synthetic samples because their controlled chemical composition reduces the number of possible configurations, facilitating the interpretation. However, the proposed interpretations should also hold for natural samples close to end-members such as those investigated by Bosi et al. (2016) and Watenphul et al. (2016). Note that reported experimental frequencies can also vary within a few wavenumbers depending on the studies and investigated samples.

The ^W(OH) bands reported by Berryman et al. (2016) at 3778 and 3739 cm⁻¹ in synthetic dravite (sample EB28 with composition Na_0.84K_0.02□ _0.13Mg_2.52Al_6.73Si_5.72O₁₈(BO${}_{\mathrm{3}}{)}_{\mathrm{3}}^{\mathrm{V}}$(OH)${}_{\mathrm{3}}^{\mathrm{W}}$(OH)) can be ascribed to ^YMg₃-^XNa and ^YAl^YMg₂-^XNa environments, respectively. A consensus exists on these assignments (Berryman et al., 2016; Watenphul et al., 2016; Bosi et al., 2016), which correspond to an upshift of theoretical frequencies with respect to their experimental counterparts amounting to ∼ 75 cm⁻¹ (Fig. 5). Similar systematic shifts have been observed in previous theoretical studies of OH defects in quartz, corundum and diopside with values ranging from ∼ 10 to 125 cm⁻¹, depending on the system (Jollands et al., 2020; Balan et al., 2020; Balan, 2020). Such a systematic difference between the theoretical and observed OH stretching frequencies results from the approximate cancelation of errors introduced by the use of the PBE XC functional and the neglect of OH anharmonicity (Balan et al., 2007).

Figure 5 Relation between theoretical (ω_OH, Table 1) and experimental (ω_exp) OH stretching frequencies in dravite. The theoretical frequencies are those of the larger $\mathrm{1}\times \mathrm{1}\times \mathrm{2}$ supercell models. Experimental frequencies are indicated (sample EB28 of Berryman et al., 2016, except the 3677 cm⁻¹ frequency from Gourdant and Robert, 1997). The dotted line corresponds to the equation ${\mathit{\omega }}_{\mathrm{OH}}={\mathit{\omega }}_{\exp }+\mathrm{75}$ cm⁻¹.

Assuming that the shift is systematic and holds for the other types of OH environments, the present results support the following assignments. The band at 3677 cm⁻¹ reported in synthetic dravite samples by Gourdant and Robert (1997) can be ascribed to ^W(OH) in a ^YMg₃-^X□ environment. Similarly, the band at 3640 cm⁻¹ reported by Berryman et al. (2016) can be related to ^W(OH) in a ^YAl^YMg₂-^X□ environment. These two interpretations are consistent with those previously proposed from the analysis of experimental data. Based on the present results, the weaker band at 3620 cm⁻¹ (Berryman et al., 2016) could however better correspond to ^V(OH) groups in a ^YMg^ZAl^ZMg-^XNa environment than to ^W(OH) groups in a ^YAl₂^YMg-^X□ environment. This last ^YAl₂^YMg-^X□ environment in dravite is expected to lead to a band at lower frequency, overlapping with the strong band related to ^V(OH) groups (Fig. 5).

Still assuming that the systematic difference between theoretical and experimental frequencies is valid for the ^V(OH) groups, the three components identified in the broad ^V(OH) band by Berryman et al. (2016) can be assigned to specific environments of OH groups. The highest component at 3577–3578 cm⁻¹ is ascribed to a ^YMg^ZAl^ZAl-^XNa environment, and the broad component at 3549 cm⁻¹ arises from contributions of ^YAl^ZAl^ZMg-^XNa (i.e., Mg–Al inversion) and ^YMg^ZAl₂-^X□ environments. The component at 3506–3513 cm⁻¹ is associated with both ^YAl^ZAl₂-^XNa and ^YAl^ZAl^ZMg-^X□ environments. Finally, the weaker component at 3446 cm⁻¹ could be assigned to a ^YAl^ZAl₂-^X□ environment.

The present theoretical results also suggest that the ^W(OH) and ^V(OH) stretching frequencies are sensitive to the ordering of X vacancies along the c axis. Accordingly, the frequencies determined for a more dilute configuration corresponding to the $\mathrm{1}\times \mathrm{1}\times \mathrm{2}$ hexagonal supercell were compared to observations made on dravite samples. In contrast, those obtained using the smaller hexagonal cell should be more relevant to samples with a high concentration of X vacancies. This is the case of magnesio-foitite, which is an Al-rich Mg–Al tourmaline with vacant X site. The Raman spectrum of the synthetic MF2 sample of Berryman et al. (2016) with a Na_0.01□_0.98Mg_1.97Al_7.24Si_5.83O₁₈(BO${}_{\mathrm{3}}{)}_{\mathrm{3}}^{\mathrm{V}+\mathrm{W}}$(OH)₄ composition displays two well-defined bands at 3619 and 3657 cm⁻¹ and a strong band at 3551 cm⁻¹. A broader asymmetric band has been fitted using two components at 3459 and 3511 cm⁻¹, the latter having a full-width at half maximum (FWHM) of ∼ 100 cm⁻¹. A comparison of theoretical and experimental frequencies can thus be done using the frequencies of vacancy-bearing models determined using the simple hexagonal cell (Fig. 6). In this comparison, the potential effect of the small unit-cell contraction experimentally observed between dravite and magnesio-foitite (−0.7 %; Berryman et al., 2016) is neglected. On the basis of the correlation between theoretical and experimental frequencies established for dravite, the two bands at 3619 and 3657 cm⁻¹ could be related to ^W(OH) groups in ^YAl₂^YMg-^X□ and ^YAl^YMg₂-^X□ environments, respectively. In this case, the upshift of the ^W(OH) frequency due the X site vacancy ordering explains why the dravite band at 3640 cm⁻¹ is not observed in synthetic Mg-foitite, which has X sites almost vacant (Veličkov, 2002; Berryman et al., 2016). Compared with dravite, the lack of ^YAl^ZMg₂-^XNa environments and the upshift of the ^YAl₂^YMg-^X□ ^W(OH) frequency lead to the well-defined band at 3619 cm⁻¹. This is a further indication that the 3620–3622 cm⁻¹ component in dravite does not correspond to ^W(OH) groups in an ^YAl₂^YMg-^X□ environment. On the same level, the band at 3551 cm⁻¹ can be ascribed to ^V(OH) in a ^YMg^ZAl₂-^X□ environment, and the component at 3459 cm⁻¹ corresponds to ^V(OH) in a ^YAl^ZAl₂-^X□ environment. Finally, the broad component at 3511 cm⁻¹ is more difficult to assign to a single environment. Contributions from the ^YAl^ZAl^ZMg-^X□ environment can be expected (Fig. 6), but its width indicates that multiple configurations, for example, related to ordered X vacancies interacting with disordered Al–Mg occupations over the Y and Z sites along the c-axis, likely contribute to this component.

Figure 6 Theoretical OH stretching frequencies of models with vacant X sites obtained using the simple hexagonal cell (Table 1) (full symbols) reported as a function of the experimental (ω_exp) OH stretching frequencies of the Mg-foitite MF2 sample of Berryman et al. (2016). For comparison with Fig. 5, the theoretical frequencies obtained for the larger $\mathrm{1}\times \mathrm{1}\times \mathrm{2}$ supercell (corresponding to a more dilute vacancy concentration) are reported as a function of the experimental frequencies of dravite (open symbols). The dotted line corresponds to the equation ${\mathit{\omega }}_{\mathrm{OH}}={\mathit{\omega }}_{\exp }+\mathrm{75}$ cm⁻¹.

Although the present assignments (Table 2) partly differ from those previously proposed by Berryman et al. (2016) and Bosi et al. (2016) using the “traditional” approach, they could help resolve the apparent inconsistency between the relative area of the spectroscopic components identified in the ^V(OH) band and a realistic inter-site distribution of cations (Bosi et al., 2016). The fit reported by Berryman et al. (2016) leads to a relative area reaching 37 % for the high-frequency component at 3577–3578 cm⁻¹, and its counterpart reported by Bosi et al. (2016) at 3563 cm⁻¹ in a natural dravite sample has a relative area of 63 %, far from the occurrence probability expected for ^YMg^ZAl^ZMg-^XNa environments (17 %) in this sample. The present assignment of the high-frequency components of the ^V(OH) band to the local ^YMg^ZAl₂-^XNa environment is thus more consistent with a higher occurrence probability of the ideal environment in dravite samples close to end-member compositions (estimated at 39 % in Bosi et al., 2016). In the same line, the assignment of the bands at 3619 and 3657 cm⁻¹ to ^W(OH) groups in ^YAl₂^YMg-^X□ and ^YAl^YMg₂-^X□ environments, respectively, is consistent with the Al-rich character of magnesio-foitite.

Table 2 Proposed interpretations of experimental bands in synthetic dravite and magnesio-foitite. Experimental frequencies (cm⁻¹) are those of samples EB28 and MF2 of Berryman et al. (2016), except the 3677 cm⁻¹ frequency from Gourdant and Robert (1997). The broad component at 3511 cm⁻¹ in Mg-foitite likely corresponds to multiple configurations.

In natural samples, the F for OH substitution and the partial deprotonation of OH groups related to Ca²⁺ incorporation at the X site also contribute to challenges in constraining compositions from OH stretching band areas. In addition, the fit of broad and asymmetric bands such as the ^V(OH) band depends on the assumptions made on the number and line shape of individual components and on the baseline subtraction (Meier, 2005). The corresponding uncertainty in the relative area of individual components thus complicates the use of vibrational spectroscopy to infer reliable inter-site distribution of cations in natural or synthetic samples when the corresponding signals overlap and form broad spectroscopic bands. In this context, the theoretical shifts of ^V(OH) frequencies predicted for specific YZZ or X occupancies could be used as independent constraints to fit the experimental ^V(OH) signal.

5 Conclusions

We have theoretically analyzed the OH stretching properties of a series of ideal and chemically disordered dravite models using DFT. The results show that the OH stretching dynamic in this system can be understood by considering an ensemble of weakly coupled single oscillators.

Localization is a general feature of OH stretching modes related to the low mass of the H atom, which ensures their significant frequency separation and related minimal coupling with the other vibrational modes. This frequency separation is commonly observed in hydrous minerals and hydrous defects in oxides and nominally anhydrous minerals (e.g., Farmer, 1974; Wöhlecke and Kovács, 2001; Libowitzky and Beran, 2006). It supports the use of partial Hessian techniques to analyze the OH stretching modes in extended systems while keeping the computational cost at a reasonable level (e.g., Ghysels et al., 2010).

The weak harmonic coupling observed between single ^V(OH) oscillators in dravite is consistent with their relatively distant positions within a single cluster (∼ 6.5 Å) or between adjacent clusters (∼ 5.5 Å). For comparison, OH groups located in the same oxygen plane in talc are ∼ 5.3 Å apart, and it is well established that their stretching frequency is locally determined by the occupancy of neighboring sites in chemically disordered samples (e.g., Farmer, 1974; Petit et al., 2004; Blanchard et al., 2018).

The present results thus support a traditional interpretative model which assumes that the local environment of a single OH group has a determining influence on its stretching frequency (e.g., Berryman et al., 2016; Kutzschbach et al., 2017; Piecza et al., 2020). Vacancies at the X site were however found to affect the ^W(OH) and ^V(OH) frequencies on a longer scale, contributing to the differences observed between the dravite and Mg-foitite spectra (Berryman et al., 2016). The same theoretical approach could be used in conjunction with a combination of other experimental techniques sensitive to short- or long-range crystal order to unravel the complex cationic ordering in other tourmaline groups or other complex hydrous silicates such as amphiboles (Hawthorne, 2016).