Modeling the polychromism of oxide minerals: The case of alexandrite and cordierite

In this work, we investigate the spectroscopic properties of photochromic alexandrite and cordierite by TD‐DFT. The objective is to assess the TD‐DFT for the simulation of pleochroism (change of color depending on the crystallographic direction of the observation) and the change of color as a function of the light source. For these simulations, we compared an embedding where dangling bonds are saturated by hydrogen atoms and an electrostatic embedding. The electrostatic embedding provided numerically more stable results and allowed a good reproduction of the pleochroism of cordierite, based on a Fe2+‐Fe3+ intervalence charge transfer transition. However, the pleochroism of alexandrite is not as well reproduced, suggesting that TD‐DFT has some difficulties to reproduce the anisotropy of the transition dipole moment, an aspect that is not deeply documented in the literature.


| INTRODUCTION
The quantum chemical modeling of spectroscopic properties in condensed phase is a challenging task since it requires the association of different computational approaches.Until recently, very few quantum chemical methods able to describe excited states were implemented with periodic boundary conditions (PBC).For a long time, the only accurate method to perform excited state calculations under PBC has been GW post-DFT calculations in combination with the Bethe-Salpeter equation.However, due to its large computational cost, the use of GW/BSE calculations to routinely investigate molecular crystals or doped oxides having hundreds of atoms in the unit cell is difficult.The computational cost is exacerbated by the fact that most PBC implementations of GW/BSE use plane-wave basis sets, where the virtual space is generally much larger than when atom-centered basis sets are used.It is important to mention the recent implementation of TD-DFT under PBC in CP2K. 1,2These implementations are still relatively unknown and deserve intensive testing to explore its predictive power.So far, TD-DFT has only been benchmarked on molecular systems. 3 efficient approach to overcome the limited simulation methods for excited states in PBC is to use an embedded cluster.This approach is particularly interesting for localized excited states, such as on a point defect 4 or on few molecules in a molecular crystal. 5,6The idea is to perform a PBC geometry optimization to have an accurate geometry of the system, followed by a cluster extraction.This nonperiodic cluster is then used in a molecular quantum chemical code where a large variety of excited state calculation based on TD-DFT or post-Hartree-Fock methods are available.The difficulties arise from (i) the choice of the cluster size and from (ii) the embedding of this cluster.The first point is very system-dependent, and a dedicated investigation of the influence of the cluster size must be performed.
For the second point, several types of embedding strategies have been proposed in the literature.The simplest embedding consists in not considering any embedding, meaning that the cluster is simulated in vacuum.This approach completely neglects the long-range influence of the solid and frequently generates meaningless orbitals (and electronic transitions) due to dangling bonds created during the cluster extraction.It is also possible to consider an embedding based on a polarizable continuum model, which is mostly done for solvated molecules. 7But this approach does not solve the problems related to atoms with dangling bonds and does not describe the heterogeneity of the environment.In order to go one step further, the outermost atoms of the cluster can be saturated with hydrogen atoms.This embedding procedure has the advantage to be straightforward and to partly remove the unphysical transitions.However, the long-range effect of the solid is still missing.
If the cluster is large, we can assume that the long-range interaction of the environment on the electronic transitions is mainly electrostatic.Thus, an array of point charges can be used to model the electrostatic embedding.The procedure developed by M. Klintenberg et al. and implemented in the code Ewald generates an array of point charges fitted such that they reproduce the electrostatic potential of the infinite periodic solid. 8,9This point charges model has been applied by several groups such as the one of Carlo Adamo and Ilaria Ciofini with the objective to reproduce spectroscopic properties of molecular solids. 5,6This approach has also been used by some of us to model the photochromism of aluminosilicate 10 and by other groups to model several properties such as NMR, 11 aggregation-induced emission, 12 or magnetic interactions. 13The results obtained from this electrostatic embedding are generally in good agreement with the experiments both qualitatively and quantitatively.This electrostatic embedding is now implemented in the FROMAGE package developed by the Crespo-Otero group. 14 thus appears relevant to apply this electrostatic embedding to investigate spectroscopic properties of solid-state systems known to be particularly challenging to model.Here, we set out to simulate the polychromism of two oxide minerals, alexandrite (BeAl 2 O 4 :Cr) and cordierite (Mg 2 Al 4 Si 5 O 18 :Fe).The alexandrite has the fascinating property to feature two different colors as function of the nature of the incoming light (sunlight and artificial light, Figure 1). 15,16The material also shows different colors if observed from different crystallographic directions.The latter effect is referred to as pleochroism in geology.In alexandrite, the colors are due to d-d transitions of the Cr 3+ dopant in the material, generating an absorption band in the visible part of the spectrum.8][19] In this case, the origin of this effect is still a source of debate in the literature, but the most probable origin of the color is an Fe 2+ -Fe 3+ intervalence transition. 20Thus, these two minerals are challenging systems to test embedding methods since they involve d-d localized transitions and d-d through space charge transfer transitions, respectively.Another interesting aspect of the pleochroism is that it requires to investigate the anisotropy of the dipole transition moment, while, in general, only the average total dipole transition is considered through the oscillator strength.
In this work, we present the effect of the choice of the embedding on the computed spectroscopic properties of alexandrite and cordierite.
We used the H-saturated embedding (noted as OH embedding) and the point charges generated by the Ewald package (noted PC).Calculations are done at the TD-DFT level, and thus, the results are sensitive to the choice of the functional.We assessed several functionals, selected to represent the most frequently used functionals for TD-DFT.
Finally, the transitions computed by the most accurate embedding scheme and functional were used to simulate the anisotropy of absorption spectrum and the colors of the two materials.

| COMPUTATIONAL DETAILS
All geometry optimizations in periodic boundary condition were performed with the CRYSTAL17 code 21 along with the global hybrid functional PBE0. 22A 12 Â 12 Â 12 k-point mesh was used for the full geometry optimization of the bulk system (including cell parameters and atomic positions in the cell).For defect containing systems, only the atomic positions in the cell were relaxed and optimized cell parameters of the pristine materials were kept.For alexandrite, the geometry of the structure containing the Cr defect was based on 1 Â 2 Â 2 supercell and the k-point mesh was reduced to 2 Â 2 Â 2. For cordierite, the geometry of the structure containing the Fe defects was based on 1 Â 2 Â 2 supercell and the k-point mesh was reduced to the Γ-point.The following basis sets were used for Fe ([6s4p3d1f]/ (17s11p6d1f)), Si ([4s3p1d]/(20s12p1d)), 23 Al ([4s3p1d]/(17s9p1d)), 24 O ([3s2p1d]/(10s4p1d)), 25 Mg ([4s4p1d]/(13s8p1d)), and Be ([4s3p1d]/(10s4p1d)) atoms.The convergence criterion for the SCF cycle was fixed at 10 À8 Ha per unit cell.All TD-DFT calculations were performed on a cluster extracted from the PBC calculations were done with the Gaussian16 code. 26Two sizes of cluster are considered (see results and discussion section), and the functionals PBE0, 22 M06-2X, 27 HSE06, 28,29 CAM-B3LYP, 30 and ωB97XD 31 have been tested.
Photos of an alexandrite exposed to two different light spectra.(B) Photos of a cordierite observed in two different directions. 16he clusters are always O-terminated.The OH embedding is obtained by saturating the outermost O atoms with H atoms.The PC embedding is obtained by adding a first layer of cations simulated by pseudo-potentials without associated basis sets, followed by a cloud of point charges (generated by Ewald package 8,9 ) that reproduces the electrostatic potential of the crystal.This approach was already tested for the simulation of the photochromism in natural aluminosilicates. 10,32 used a previously reported approach to simulate the perceived color. 33,34First, the absorption spectra were turned into transmission spectra.The resulting RGB coordinates were then calculated based on the model of the 1931's International Commission on Illumination (CIE).Data for a standard observer and the spectral density of sunlight at 6504 K (D65 illuminant) or for indoor light spectrum (A illuminant) were used.

| Structural information
The alexandrite and the cordierite crystallize in the Pmna and Cccm orthorhombic space groups, respectively.Alexandrite is a network of edge-shared octahedron around Al 3+ ions and corner-shared tetrahedron around Be 2+ ions (Figure 2A).Cordierite is an aluminosilicate mineral based on corner-shared tetrahedron around Al 3+ and Si 4+ ions and edge-shared octahedron around Mg 2+ ions.The structure of this last material is highly anisotropic with channels along the c crystallographic direction (Figure 2B).Atomic positions and cell parameters are provided in Supporting Information.
In terms of computed cell parameters, the discrepancy between theory and experiment is below 0.3% for alexandrite and below 0.4% for cordierite, confirming the reliability of the level of theory chosen to determine the geometry in PBC.
The Cr-doped alexandrite is obtained by substituting one Al 3+ ion by a Cr 3+ ion in a 1 Â 2 Â 2 supercell, leading to a doping content of 0.89 at.% in the overall system.This doping is higher than in the sample used in the work of R. Powell et al. (0.09 at.%), which provides the experimental reference in this article. 35Nevertheless, even at 0.89 at.%, the interaction between Cr 3+ ions can be considered negligible.In octahedron geometry, the Cr 3+ spin (noted S) is equal to 3/2.The average Cr-O bond length (1.95 Å) is computed to be slightly larger than the average Al-O bond length (1.89 Å) in agreement with the larger ionic radius of Cr 3+ compared to Al 3+ .Two clusters centered on Cr 3+ are then extracted from this structure.A small cluster including only the six O 2À ions of the first coordination sphere of Cr 3+ and a another larger cluster including also the first polyhedrons of Al 3+ and Be 2+ around the Cr 3+ that corresponds to the second coordination sphere of Cr 3+ (Figure 3A,B).
The Fe-doped cordierite is built by substituting one Mg 2+ by one Fe 2+ in an octahedron and one Al 3+ by one Fe 3+ in a tetrahedron having an edge shared with the Fe 2+ octahedron.We used a 1 Â 2 Â 2 supercell leading to a total Fe content of 0.4 at.% on the overall system.The system is computed in a high spin state (S(Fe 2+ ) = 2 and S(Fe 3+ ) = 5/2 thus S(total) = 9/2) based on experimental spectroscopic characterizations. 36The computed Fe 2+ -O bond length   3C) and a larger cluster (Figure 3D) including also the first Al 3+ , Si 4+ , and Mg 2+ polyhedrons around the Fe 2+ and Fe 3+ defects.

| Spectroscopic properties
Figure 4A presents the computed and experimental transition energies of alexandrite.The absorption spectrum of this mineral is characterized by two spin-allowed transitions, the 4 A 2 ! 4 T 1 ($420 nm experimentally) and the 4 A 2 ! 4 T 2 ($580 nm experimentally). 37Due to the octahedron distortion of the alexandrite structure, there is a degeneracy removal of the 4 T 1 and 4 T 2 states.In practice, it means that several transitions are obtained.These transitions are reported in Figure 4 as small horizontal lines defining the vertical spread of the computed transitions.The results for the large OHsaturated clusters suffered from significant spin contamination and were, therefore, judged unreliable.Hence, we do not report them.No noticeable spin contamination was obtained for the large cluster saturated by point charges.
The nature of the embedding (OH vs PC) was found to have only a small effect, with the large cluster slightly reducing the energy of the excited states.The largest effect comes from the functionals.It is Computed transition energies as a function of the cluster size (small and large) and of the nature of the embedding (H-saturated noted OH and point charges noted PC).Each dark line corresponds to a computed transition energy.The experimental values are centered on the maximum absorption wavelength, and the size of the box includes the full width at half maximum.(A) For alexandrite and (B) for cordierite.known that d-d transition energies are strongly sensitive to the exchange-correlation functional due to the large electronic correlation in the d-shell of the 4th row of the periodic table. 38A large fraction of exact exchange seems necessary to reproduce the electronic transition for instance as provided by M06-2X (54% of exact exchange).
The computed energies and dipole transition moment obtained with this functional for the large PC-embedded cluster are used for further analysis in the next subsection.
The computational and experimental transition energies of the cordierite are represented in Figure 4B.Experimentally, the absorption is characterized by a large band centered at 575 nm.Computationally, several electronic transitions correspond to an intervalence charge transfer.The energy range presented in Figure 4B corresponds to the range of computed energies for each system and each functional.As matter of illustration, the electron variation density of the most intense intervalence transition computed with CAM-B3LYP on the large cluster with PC embedding is given in Figure 5.The associated d CT charge transfer index computed is 1.1 Å, notably lower than the 2.8 Å Fe-Fe distance, but the large q CT value (1.7 e) still confirms the charge transfer character of this transition. 39,40e behavior of the transition energies computed on cordierite as function of the cluster size and of the embedding is notably different than for the alexandrite.For the OH embedding, the computed transition energies are very different from the small and the large clusters.
However, whatever the functional considered, none of the OH-saturated clusters provides transition energies in agreement with experiment.On the contrary, the PC embedding generally provides transition energies much less sensitive to the cluster size.The most accurate functionals seem to be the range-separated ones that are designed for simulating charge transfer transitions.Interestingly, the M06-2X functional that was the best for simulating the localized d-d transition of Cr 3+ in alexandrite provides here the worst results for the d-d charge transfer transition.We here find that the "charge transfer" character of the transition is more important than its "d-d" character since the embedding has a larger influence on the simulation than the choice of the functional.Thus, a range-separated functional has to be used, as they are the only ones adapted for charge transfer transitions The computed energies and dipole transition moment obtained with the CAM-B3LYP functional for the large cluster embedded in point charges are used in the next subsection.

| Absorption spectra and polychromism simulations
The absorption spectra of the alexandrite and cordierite were simulated by convoluting TD-DFT computed electronic transitions with Gaussian functions having a full width at half maximum extracted from the experimental absorption spectra.For the pleochroism simulations (i.e., anisotropy of absorption spectrum), a first step consists in the simulation of the absorption spectra for polarized light with an electric field parallel to one of the crystallographic directions (depending to the anisotropy observed experimentally).To simulate the E//a absorption spectrum, for instance, the transition intensities were simulated with an oscillator strength computed only using the coordinate of transition dipole moment parallel to a-axis.The simulated absorption spectra of polarized light are given in Figure 6.
For alexandrite, the relative intensities of the 580 and 420 nm bands are qualitatively well reproduced for the E//a spectrum, although the intensity of the 580 nm band is too low.However, the E//b absorption spectrum is not well reproduced by the simulation.For this direction, the 580 nm band is experimentally found to be more intense than the one at 420 nm.Even though the 580 nm band is computed to be more intense in the E//b spectrum than for E//a, it remains significantly less intense than the 420 nm band.Also with the other functionals the same behavior is obtained (see Figure S1).This highlights the difficulty of TD-DFT to reproduce the dipole transition moment of d-d transitions.Interestingly, transition intensities are pretty well reproduced for cordierite with almost no absorption in the visible part of the spectrum computed for E//a and E//c and a large band centered at 580 nm for E//b, both computationally and experimentally.
To obtain the final color upon an observation along a specific crystallographic direction, let's say a-axis (corresponding to a propagation vector k parallel to a, noted k//a), the absorption spectrum along this direction is obtained by summing the E//b and E//c simulated absorption spectra since only the electric field perpendicular to the propagation vector is observed.The final colors are given in Figure 7A.Even if the experimental intense pleochroism of alexandrite is not well reproduced, we can note a slight change of the simulated color.The discrepancy between simulation and observation comes from the poor simulation of the relative intensity between the 420 and 580 nm bands when changing the direction of observation.This originates from the wrong simulation of the relative intensities between the two absorption bands for an electric field parallel to the b-direction (Figure 6B).For cordierite, the pleochroism from white to F I G U R E 5 Computed variation of the electron density associated to the most intense transition obtained with CAM-B3LYP (large cluster with PC embedding).Green and red areas correspond to electron density decrease and increase, respectively.blue-purple is pretty well reproduced following the good agreement between theory and experiment of the polarized light spectra.
In the specific case of alexandrite, a polychromism is also experimentally observed when changing the nature of the light.To simulate this phenomenon, the full absorption spectrum of this material was modeled considering all the crystallographic direction (taking the standard oscillator strength) and the color was simulated considering either the D65 (sunlight) or the A (indoor light) spectra.The final color is presented in Figure 7B.The experimental color for alexandrite reported by K. Schmetzer is "blue-green" under sunlight and "red-purple" under indoor light which is relatively well reproduced by calculations with a "dark-green" color under sunlight and a "dark-orange" color under indoor light.In that case, the final color is more dominated by the positions of the absorption bands than by their relative intensities, explaining why we obtained a better simulated color than for the pleochroism (in the different crystallographic direction).
F I G U R E 6 (A-C) Simulated and experimental absorption spectra of alexandrite for the three light polarizations.(D-F) Simulated and experimental absorption spectra of cordierite for the three light polarizations. 17,35Simulated colors from experimental and simulated absorption curves are presented in inset.

| CONCLUSION
The color of alexandrite and cordierite changes as function of the direction of observation (pleochroism) or as a function of the spectrum of the incident light.In this work, we take up the challenge to simulate the polychromism of alexandrite and cordierite.From our computational investigation, we reach several conclusions: • First, on the embedding schemes.In agreement with previous investigations, the point charge embedding provides numerically more stable results than caping O 2À with hydrogens. 41In our present study, the ground-state wave functions suffered less of spin contamination, we obtained fewer "artificial" transitions due to dangling bonds, and the overall results are less sensitive toward the size of the cluster.However, the PC embedding cannot be considered as a generally accurate embedding scheme.It lacks polarization of the F I G U R E 7 (A) Simulated and experimental colors of alexandrite and cordierite as a function of the crystallographic observation.(B) Simulated and experimental colors of alexandrite as a function of the nature of the illumination. 15,16nvironment and only accounts for electrostatics.The future of the embedding maybe lays on density (matrix) embedding 42,43 or the use of absolutely localized molecular orbitals (ALMO), 44,45 which go beyond the electrostatics of the crystal structure.
• Second regarding the computed transition energies and intensities.In this work, we faced the well-documented difficulty to treat the d-d electronic transitions by TD-DFT. 38Of course, the best strategy would be to apply post-Hartree-Fock methods, where the main challenge is to deal with large clusters.However, our study shows that small electrostatically embedded clusters give similar results to larger clusters, suggesting that multi-reference computations could be applied.This will be explored in future work.TD-DFT is well documented for its performance for molecules in solution or gas phase where the transition intensities are averaged in all direction of space through the oscillator strength. 46With the example of alexandrite, we herein demonstrate the challenge to reproduce not only the transition energies but also the anisotropy of the dipole transition moment.
• Third, on the final simulated spectra and colors.For alexandrite, the change of color by modifying the light source is well reproduced since it is mainly related to a good reproduction of the transition energies, which is achieved by the M06-2X functional.However, the pleochroism simulation, that is, the change of color as a function of the direction of observation, is less satisfactory since the relative change of transition intensities is not properly captured by M06-2X (and any other functionals).For the cordierite, the agreement between theory and experiment is near quantitative both for the spectra and for the colors, although a range-separated functional (CAM-B3LYP), rather than a global hybrid, was necessary to achieve this.Our results support the hypothesis that the polychromism of cordierite is due to an intervalence charge transfer transitions between Fe 2+ and Fe 3+ .
In summary, the simulation of polychromism of minerals should attract more attention by the theoretical chemistry community as it turns out to be a challenging task that can be used to benchmark excited states methods.It is also a very interesting case to assess the recent implementations of TD-DFT methods in PBC that could involve in the future electron-phonon coupling for a more accurate simulation of transition energies and absorption band shapes.
Representation of the structure of alexandrite (A) and cordierite (B).Oxygen atoms are in red.Pink, light blue, dark blue, and yellow polyhedrons correspond to Be 2+ , Al 3+ , Si 4+ , and Mg 2+ polyhedrons, respectively.For cordierite, four unit cells are presented and the view is along the c-axis.F I G U R E 3 Representation of the small clusters of alexandrite (A) and cordierite (c) and of the large clusters of alexandrite (B) and cordierite (D).Oxygen atoms are in red.Pink, light blue, dark blue, and yellow polyhedrons correspond to Be 2+ , Al 3+ , Si 4+ , and Mg 2+ opaque polyhedrons, respectively.Orange, dark green, and light green transparent polyhedrons correspond to Fe 2+ , Cr 3+ , and Fe 3+ polyhedrons, respectively.

( 2 .
13 Å) is close to the one of Mg-O bond length (2.09 Å) in agreement with their similar ionic radius.The Fe 3+ -O bond length (1.87 Å) is notably longer than that the Al 3+ -O one (1.75Å) again in agreement with their respective ionic radii.The methodology of cluster extraction used for cordierite is similar to the one used for alexandrite: a small cluster having only the first coordination sphere of O 2À ions (Figure