Study of the variation of the optical properties of calcite with applied stress, useful for specific rock and material mechanics

Calcite (CaCO3, trigonal crystal system, space group R3¯c\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R\overline{3}c$$\end{document}) is a ubiquitous carbonate phase commonly found on the Earth’s crust that finds many useful applications in both scientific (mineralogy, petrology, geology) and technological fields (optics, sensors, materials technology) because of its peculiar anisotropic physical properties. Among them, photoelasticity, i.e., the variation of the optical properties of the mineral (including birefringence) with the applied stress, could find usefulness in determining the stress state of a rock sample containing calcite by employing simple optical measurements. However, the photoelastic tensor is not easily available from experiments, and affected by high uncertainties. Here we present a theoretical Density Functional Theory approach to obtain both elastic and photoelastic properties of calcite, considering realistic experimental conditions (298 K, 1 atm). The results were compared with those available in literature, further extending the knowledge of the photoelasticity of calcite, and clarifying an experimental discrepancy in the sign of the p41 photoelastic tensor component measured in past investigations. The methods here described and applied to a well-known crystalline material can be used to obtain the photoelastic properties of other minerals and/or materials at desired pressure and temperature conditions.

www.nature.com/scientificreports/ …, 6 can be used, where 1 = xx, 2 = yy, 3 = zz, 4 = yz, 5 = xz and 6 = xy 5 . With this notation, the fourth-rank elastic tensor C is represented by a symmetric 6 × 6 matrix and the stress-strain relation can be rewritten as σ v = C vu η u and η v = S vu σ u . In terms of crystal lattice, the components of the stiffness matrix are defined as with Ω the unit cell volume and E the energy of the system. This relationship shows that the terms C vu and C uv are equivalent, resulting in the symmetry of the C and S matrices and in the further reduction from 36 to 21 independent elastic components. For all but isotropic (cubic) crystalline materials, the dielectric constants are described by a second-rank tensor ǫ , which determines the optical indicatrix of the system. In the case of calcite, the optical indicatrix is negative, because the extraordinary refractive index, n ε , is lower than the ordinary one, n ω . When the crystal is strained, the variation of the inverse dielectric tensor ǫ −1 can be related to the amount of induced strain by means of the fourth-rank photoelastic (Pöckels') tensor p, whose components are known as elasto-optic or strain-optical coefficients. Hence, the elasto-optic constants can be calculated as where �ǫ −1 ij are the differences between the strained and unstrained inverse dielectric tensor. As in the case of the elastic tensor, it is possible to reduce the fourth-rank Pöckels' tensor to a 6 × 6 matrix, where p vu = ∂�ǫ −1 v ∂η u .
(3) www.nature.com/scientificreports/ Given the stress-strain relationship previously introduced, the piezo-optic tensor π, i.e., the fourth-rank tensor whose π vu components correlate the �ǫ −1 v variations with the stress σ u , can be obtained with and It is worth mentioning that, although the elastic tensor of a material is symmetric, the photoelastic one is generally not, thus p vu ≠ p uv and π vu ≠ π uv and both p and π have 36 independent components. This knowledge finds its usefulness in the determination of the stress state of rocks and materials and/or for the analysis of fractures in rock samples containing calcite, if the photoelastic and/or piezo-optic properties are accurate enough 6,7 .
The photoelasticity of calcite was experimentally determined and widely published in the past literature, where we observed and reviewed quite scattered results, with the same photoelastic tensor component showing high variation up to about 200%, leaving the reader doubtful and confused (see for instance refs [8][9][10][11][12] ). This wide range and scattering of results is associated to the experimental methods employed to characterize photoelasticity and piezo-optics of materials, which are complex and involve the measurements of different quantities, e.g., elastic properties and refractive indexes, on specifically cut and prepared samples [13][14][15] . It must be also emphasized that the just cited properties were determined in several studies carried out by various researchers, who performed the measurements in different experimental conditions and instrumental setups, resulting in scattered data with different accuracies 16 . For example, the calcite photoelasticity was measured by variations of refractive index 8,9 , Brillouin scattering 10 and analysis of Raman spectral intensities 12 , where each technique presented different degrees of accuracy. In fact, as also discussed by Andrushchak and co-workers 17 , some of the employed experimental means (e.g., Brillouin scattering techniques) are accompanied by very high uncertainties, which often lead to ill-defined photoelastic/piezo-optic constants, in particular on their sign. Indeed, this is the case of the photoelastic tensor component p 41 of calcite, whose sign disagrees between the different experiments as reported in the references 10,18 .
Here, we propose a theoretical (Density Functional Theory, DFT) investigation of the photoelastic and piezooptic properties of calcite to further extend the knowledge of the optical properties of this material, and to aid answering the above cited uncertainties and cross-correlating the simulated and experimental data. This work aims also at showing the interested reader a possible way to calculate and model this physical property for other minerals and/or materials of interest for both mineralogical, geological and materials science applications.
Structural and dielectric properties of calcite. We performed first-principle DFT simulations to calculate the photoelastic properties of calcite, both at absolute zero (standard DFT settings) and in typical experimental conditions (298 K, 1 atm), using two well-known approaches, PBE-D2 and B3LYP-D*, which include the effects of long-range interactions (see the Methods section for details).
The necessary starting point to calculate the photoelastic properties of calcite is given by a good description of both the crystal structure and the dielectric properties. The simulation results related to the equilibrium geometry of calcite are reported in Table 1. It can be noted that the inclusion of long-range interactions via a semi-empirical scheme produced a unit cell that is closer to the X-ray diffraction refinement found in literature 19 and reported in Table 1. At 298 K and with the inclusion of van der Waals correction, the unit cell volume is in line with that of the XRD refinement at the same conditions, with only a slight overestimation for both functionals (ΔΩ PBE-D2 = + 2.0%, , dielectric tensor components (static ǫ 0 and high-frequency ǫ ∞ ), refractive index (ordinary n ω and extraordinary n ε ) and birefringence (δ) of calcite, obtained at DFT level with PBE-D2 and B3LYP-D* functionals corrected for long-range interactions in static condition (0 K) and at 298 K (at 1 atm = 0.0001 GPa). a 19 ; b24 ; c20 ; d Present work (0 K); e Present work (298 K). www.nature.com/scientificreports/ ΔΩ B3LYP-D* = + 2.3%). To our knowledge, this is the first time that a proper treatment for long-range interactions is included in the calculation of the total DFT energy for calcite at room temperature. As expected, because of the heterodesmic nature of calcite, the dispersive forces are more effective on the c lattice parameter than on the a = b one, since in the trigonal unit cell of the mineral there are layers of Ca 2+ ions alternately stacked (along the z direction) with layers of CO 3 2anions (see Fig. 1). In this regard, the proposed approach including the dispersive forces correctly describes the Ca 2+ -O 2attractive interactions, resulting in an optimized distance between the layers, and gives a correct description of the unit cell of a crystal, which is mandatory for obtaining reliable structure-dependent properties, such as vibrational 20,21 , elastic 22 and photoelastic ones 23 . Table 1 reports also the static and high-frequency dielectric constants ( ǫ 0 and ǫ ∞ , respectively), the refractive index n and the birefringence δ. At ambient conditions (298 K, 1 atm), the components of the static dielectric tensor along the xx and zz directions are in good agreement with the experimental and theoretical ones tabulated in literature ( ǫ 0 xx = 8.5, ǫ 0 zz = 8.0) 24,25 . In general, PBE-D2 data are slightly higher than those calculated at B3LYP-D* level, with the maximum difference of ca. + 18% on ǫ 0 xx .
Elastic and photoelastic properties of calcite. We report in Table 2 the elastic (C vu ), photoelastic (p vu ) and piezo-optic (π vu ) constants calculated at PBE-D2 and B3LYP-D* level of theory. Let us discuss first the elastic moduli, where it can be noted that at 298 K both functionals corrected for the long-range interactions provide a good description of the stiffness tensor components, with respect to the experimental results of Dandekar 26 and Chen and co-workers 27 , with values slightly increased by about 10%. In fact, it is well-known that, even by including the thermal effects, there is a small overestimation of the elastic moduli because of the Pulay stress, i.e., an effect due to the incompleteness of the basis sets. This effect occurs during the derivation of the basis sets with respect to the position of the atoms and slightly increases the values of the elastic tensor components.
The elastic constants C 11 , C 33 , C 44 and C 66 obtained with the hybrid B3LYP functional are about 2-6% higher than those calculated with PBE, while the other stiffness components show an underestimation of about 1%. This is a common figure of the two approaches because the standard generalized-gradient approximation functional is associated to both underbinding between atoms and underestimated bulk moduli for ionic compounds, as recently shown by Zhang and co-workers 28 . For solids as calcite, the use of hybrid functionals that include a fraction of exact Fock exchange is preferred because they are more accurate in the simultaneous description of both structural (lattice) and energy (e.g., cohesive) properties 22,28,29 .
The photoelastic p vu and piezo-optic π vu components obtained from our DFT simulations (at 0 K and 298 K) for ω = 0 are reported in Table 2, together with those obtained by different experimental measurements 9-12 . In Table 2. Calcite elastic (C vu , GPa), photoelastic (p vu , dimensionless) and piezo-optic (π vu , TPa -1 ) components of the corresponding fourth-rank tensors, expressed in Voigt's notation. a 26 , b27 , c10 , d9 , e11 , calculated using results from Nelson et al. 10 . www.nature.com/scientificreports/ the next paragraphs, we discuss the present simulation results with those experimentally derived by Pöckels 9 and Nelson et al. 10 .
Our theoretical study suggests that the p 41 photoelastic constant is negative in the selected crystal orientation (crystallographic a-axis and c-axis parallel to the -x and z Cartesian directions, respectively). As also suggested by Erba and Dovesi 23 , this discrepancy on the sign of the photoelastic component could arise from the finite wavelength of the laser sources used in the different experimental studies. Indeed, our theoretical photoelastic and piezo-optic constants of calcite reported in Table 2 were obtained at zero electric field frequency (ω = 0, λ = ∞), but it is known that the experiments generally employ monochromatic lasers centred at a specific wavelength, hence the experimentally derived photoelasticity does not correspond to that of the static limit. To check if the finite wavelength could be the same source of discrepancy experimentally observed in calcite, and at the same time to show how it affects the p vu components, we calculated the fourth-rank tensor p at different λ between 300 and 1000 nm (hence, at electric field frequency ω ≠ 0) at the PBE-D2 level of theory, both at room temperature (298 K) and in static conditions (0 K). The results are graphically reported in Fig. 2, showing that all but the p 14 tensor component have a dependence on the wavelength λ.
In detail, there is a strong variation of the strain-optical components for 300 nm < λ < 500 nm, then the p vu tensor components asymptotically reach the value calculated at the λ = ∞ limit, with a trend of the form p vu ( ) = −k + p vu (∞) within the 300 nm-1000 nm range. The p 41 value as a function of λ is always lower than zero (Fig. 2b), and no change of sign in the explored wavelength range was observed. Hence, this analysis showed that for calcite the experimental ambiguity of the sign of the p 41 photoelastic component is not related to the laser source employed in the experiments, but it could be related to some instrumental uncertainties or in the orientation of the crystal during the measurements carried out by Pöckels 9 . Similar discussions were proposed in a relatively recent review by Mytsyk 30 , who explained that for some crystal classes (such as 3m ) the choice of the reference system can influence the sign of shifting and rotating piezo-optic components.
Most of the p vu components show a variation between the values calculated at room and absolute zero temperatures, Δp vu = p vu (298 K) -p vu (0 K) (Fig. 2c,d), which asymptotically reaches zero as a function of λ. The calculated Δp vu are in absolute terms higher for λ < 500 nm. The exceptions to this trend are given by the p 14 and p 44 components, whose temperature differences are almost constants (Fig. 2d).
Our work provides the first ab initio determination of the photoelastic and piezo-optic tensor components at 0 K and at 298 K of calcite, together with their related quantities (crystal structure, optical and elastic properties), showing the relevancy of the use of quantum mechanical methods in this research field. The p and π tensors are very stable with respect to internal DFT parameters, meaning that theoretical simulations performed within this framework are suitable for this kind of characterization. This work further extends the knowledge of important physical properties of anisotropic and heterodesmic crystals, providing a complete description as a function of λ of the photoelastic (strain dependent) and piezo-optical (stress dependent) response of calcite, which can be employed to analyse the state of stress and strain within a geological sample and/or a material and envisage other possible applications. In addition to being useful to clarify potential ambiguities in the sign of photoelastic/piezo-optic constants determined experimentally, the theoretical approach here exploited is indeed extremely convenient and useful to cross-check the data, because experimental acousto-optical techniques are cumbersome and complex. In fact, Andrushchak and collaborators 16 showed that with optical interferometry 57 measurements on 16 specifically cut samples are necessary to obtain the 36 independent photoelastic constants, and, for a trigonal crystal of the 3m class, 11 measurements on 2 samples are still needed. The present approach could be extremely useful also to model the photoelastic and piezo-optic behaviour of other anisotropic crystals, such as those recently developed by Tear and co-workers 31,32 .

Methods
Computational approach. In the present work, the investigation of the photoelasticity of calcite was conducted within the Density Functional Theory (DFT) framework, employing the CRYSTAL17 code 33 . Regarding the Hamiltonian, we used both the standard generalized gradient approximation (GGA) functional developed by Perdew-Burke-Ernzerhof (PBE) 34 , and the hybrid B3LYP 35,36 . We studied the effects of different DFT functionals on the calculated properties of the material and, in this case, both PBE and B3LYP are well-established functionals for solid state physics applications. However, a typical drawback of the selected Hamiltonians is the lack of adequate treatment of long-range interactions, which we here accounted for by means of the DFT-D2 scheme 37 . To this purpose, we employed a specific parametrization of the selected correction for dispersive forces for B3LYP, labelled as B3LYP-D* 38 . The multi-electron wave function was constructed using a linear combination of atomic orbitals, which in turn are described as Gaussian-type functions. We employed for calcite a basis set that was optimized in the previous work of Valenzano and co-workers 20 , more specifically the one labelled as BSD (basis set D) in the cited work. The total energy was calculated by sampling the First Brillouin zone with a 6 × 6 × 6 Monkhorst-Pack mesh (32 irreducible k-points), setting the tolerances for the self-consistent field loop to 10 -8 Ha.
Dielectric tensor. The electronic dielectric tensors at equilibrium, ǫ(0) , and at each deformation step, ǫ(η) , were calculated according to the following expression reported by Erba and Dovesi 23 : with α and Ω being the electronic polarizability tensor and the unit cell volume, respectively. The polarizability is evaluated analytically with a coupled-perturbed Hartree-Fock/Kohn-Sham (CPHF/KS) approach adapted for solid systems as described by Ferrero and co-workers [39][40][41][42] . It is a self-consistent, perturbative method that www.nature.com/scientificreports/  10 , whose respective data at such wavelength are reported together with the theoretical ones. The black symbols are the photoelastic constants calculated in our simulations in static limit (λ = ∞). Panels (e) and (f) report the differences Δp vu between the calculated room temperature and 0 K photoelastic constants as a function of λ. www.nature.com/scientificreports/ describes the relaxation of the crystalline orbitals when affected by an external electric field E . Taking the electronic charge as e = -1 a.u., and as indicated by Gu et al. 43 , the general expression of the Hamiltonian Ĥ from this approach is given by: where Ĥ 0 and Ĥ ′ are the unperturbed and perturbed Hamiltonian operators, respectively, the latter being the sum of the static electric field, E st (frequency ω = 0, wavelength λ = ∞), and the frequency-dependent one, whose amplitude and frequency of oscillation are E ω and ω, respectively. r is the electron position vector in real space. However, the scalar potential operator Ĥ ′ = � E · � r is unbounded and breaks the translational symmetry, hence, for periodic structures, it is replaced by the following formula, as suggested by Maschio et al. 44 : with k indicating a reciprocal space vector, E b (t) and ĥ (E b (t)) being the b component of the applied electric field and of the gradient vector, respectively. The total energy of the system, E tot , under the effect of the electric field perturbation can then be expressed in terms of a perturbative series of the field components, indicated with the indexes i, j and k: where E (0) tot is the unperturbed total energy and the tensors of increasing rank μ, α and β are the permanent electric dipole moment, the polarizability and the hyperpolarizability, respectively, as reported by Ferrero et al. 45 . The interested reader can find in dedicated literature [43][44][45][46] all the information related to CPHF/KS approach for periodic systems, and the detail of its implementation in the CRYSTAL code [39][40][41][42] , which allows calculating the electronic contribution to the dielectric tensor for both static (frequency ω = 0, wavelength λ = ∞) and frequencydependent (ω ≠ 0, λ ≠ ∞) electric fields.
Elastic and photoelastic constants. The elastic moduli C vu and photoelastic p vu constants were calculated using the following three independent strains of the unit cell, i.e., η 1 = η 2 , η 3 and η 4 = η 5 = η 6 , which were sufficient to obtain each independent elastic and photoelastic component for a rhombohedral-I crystal class 23,47 : For each independent strain, we employed 5 deformations with strain amplitude δ ± 0.010 Å (step of 0.005 Å). The internal geometry (atomic positions) was optimized at fixed lattice parameters to include the nuclear relaxation, and then it was calculated the electronic dielectric tensors at both equilibrium, and strained configuration, ǫ(0) and ǫ(η) , respectively 23 .
The crystallographic axes of calcite were oriented with the crystallographic a-axis and c-axis parallel to the -x and z Cartesian directions, respectively. Within this convention, and according to Nye 5 , the independent and non-zero components of the stiffness tensor C for a crystal of rhombohedral-I class are: where the dots in the 6 × 6 matrix are null elements and with C 66 = (C 11 -C 12 )/2. Similarly, we can define the photoelastic (Pöckels') tensor p as: where p 66 = (p 11 -p 12 )/2. By previous definition, the piezo-optic tensor π has the same structure as that of the photoelastic components:  www.nature.com/scientificreports/ As shown in the above formulation, the piezo-optic matrix π can be subdivided in four 3 × 3 sub-matrices A, B, C and D. When the indexes v, u = 1, 2, 3, the π vu components are called the principal coefficients, which describe the relation between the principal refractive indexes and the normal stresses (sub-matrix A). For v = 1, 2, 3 and u = 4, 5, 6 (sub-matrix B), the piezo-optic components are referred as the shifting coefficients, connecting the variations of the principal refractive indexes with shear stresses. Sub-matrices C (v = 4, 5, 6 and u = 1, 2, 3) and D (v, u = 4, 5, 6) contain the so-called rotating and rotating-shifting piezo-optic components, which describe the rotation of the optical indicatrix under the effect of normal and shear stresses, respectively 17 .
The calculation of the structural, dielectric and (photo)elastic properties were conducted both in static conditions, i.e., at 0 K without any thermal contribution, and at room temperature (298 K) at a pressure of 1 atm (= 0.0001 GPa). For the latter, we employed the quasi-harmonic approximation to introduce the temperature effect on the cited properties, as described in previous literature [48][49][50] . Five unit cell volumes, the equilibrium one, two compressed and two expanded, were used for the quasi-harmonic approximation. Finally, we employed the so-called "quasi-static approximation (QSA)" to describe the thermo-elasticity of calcite, which assumes that the stiffness depends only on the thermal expansion of the crystal as suggested by Destefanis et al. 51 , where, as explained, it was shown that QSA provides a qualitatively good description of the thermo-elastic constants.