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 [Formula: see text]) 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.

Introduction

From the mineralogical perspective, calcite (CaCO3, space group , rhombohedral-I class) is an anisotropic crystalline phase formed by layers of anions, with covalent C–O bonds, alternately stacked to layers of Ca2+ ions held together by ionic interactions[1]. In addition, this mineral can be considered a model of heterodesmic structure, because it presents different types of bonds between the atoms and/or atomic groups. A graphical representation of the crystal structure of calcite is presented in Fig. 1.

Figure 1

Rhombohedral cell of calcite, viewed along two crystallographic directions, [100] and [001] in the upper and lower panels, respectively. In the upper panel, the stacking of and Ca2+ layers are clearly visible, whereas in the lower panel, the polyhedrons indicate the octahedral coordination of the Ca2+ ions in between the oxygen belonging to the carbonate groups.

Due to its structural, elastic, and optical properties, calcite is an important mineral in the scientific and technological fields, with wide applications in mineralogy, geology, chemistry, materials science, optics, and instrumentation engineering[2-4]. Calcite is one of the most birefringent materials (δ = 0.172, with δ being birefringence), where its transparent variety is known as Iceland spar. In general, birefringence derives from the anisotropic refractive indexes of the mineral, which in turn are related to its dielectric (second-rank) tensor. When a crystal is deformed by a stress (described with a second-rank tensor), its optical properties change according to the photoelastic (or elasto-optics, p) and piezo-optics (π) fourth-rank tensors, which correlate the variations of the inverse of the dielectric tensor with the applied deformation. The relationship between the second-rank stress tensor σ to the second-rank pure strain tensor η is given by the generalized Hooke’s law,andwhere C and S are the components of the fourth-rank stiffness and compliance tensors C and S, respectively, with S = C–1. The indexes i, j, k and l represent Cartesian directions x, y and z. Both stress and strain tensors are symmetric, hence their independent components are 6 and the 81 elastic components reduces to 36. By using the Voigt’s notation, pairs of Cartesian indexes can be mapped by a single index, so that just two indexes v, u = 1, …, 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 and . In terms of crystal lattice, the components of the stiffness matrix are defined aswith Ω the unit cell volume and E the energy of the system. This relationship shows that the terms C and C 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 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 aswhere 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 . Given the stress–strain relationship previously introduced, the piezo-optic tensor π, i.e., the fourth-rank tensor whose π components correlate the variations with the stress σ, can be obtained withand

It is worth mentioning that, although the elastic tensor of a material is symmetric, the photoelastic one is generally not, thus p ≠ p and π ≠ π 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-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-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 p41 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 piezo-optic 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%, ΔΩ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 Ca2+ ions alternately stacked (along the z direction) with layers of CO32– anions (see Fig. 1). In this regard, the proposed approach including the dispersive forces correctly describes the Ca2+–O2– attractive 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

Calculated structure (unit cell parameters a and c, volume Ω, C–O bond distance dC–O and Ca–O interaction distance dCa–O), dielectric tensor components (static 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).

	Experimental	PBEc	PBE-D2d	PBE-D2e	B3LYPc	B3LYP-D* d	B3LYP-D*e
a (Å)	4.991a	5.039	5.020	5.032	5.037	5.028	5.041
c (Å)	17.062a	17.402	16.955	17.118	17.330	16.968	17.116
Ω (Å3)	368.1a	382.7	370.1	375.4	380.8	371.5	376.7
ΔΩ (%)	0a	+ 4.0	+ 0.5	+ 2.0	+ 3.5	+ 0.9	+ 2.3
dC–O (Å)	1.2840a	1.2978	1.2952	1.2963	1.2878	1.2860	1.2870
dCa–O (Å)	2.3590a	2.3900	2.3607	2.3733	2.3907	2.3695	2.3818
ϵxx0	8.5b	–	11.26	10.53	7.8	8.54	8.40
ϵzz0	8.0b	–	9.15	8.11	6.4	7.50	7.00
ϵxx∞	2.7b	–	2.80	2.78	2.6	2.59	2.57
ϵzz∞	2.2b	–	2.20	2.19	2.1	2.08	2.06
nω	1.640–1.660b	–	1.674	1.668	–	1.609	1.603
nε	1.486b	–	1.483	1.479	–	1.441	1.437
Δ	0.1540–0.1740b	–	0.191	0.189	–	0.168	0.166

a[19]; b[24]; c[20]; dPresent work (0 K); ePresent work (298 K).

Table 1 reports also the static and high-frequency dielectric constants ( 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 ( = 8.5,  = 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 .

Elastic and photoelastic properties of calcite

We report in Table 2 the elastic (C), photoelastic (p) and piezo-optic (π) 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.

Table 2

Calcite elastic (C, GPa), photoelastic (p, dimensionless) and piezo-optic (π, TPa–1) components of the corresponding fourth-rank tensors, expressed in Voigt’s notation.

	Elastic constants
	C11	C33	C12	C13	C14	C44	C66
PBE-D2 (0 K)	161.47	84.01	66.09	59.75	– 20.80	35.64	47.69
PBE-D2 (298 K)	154.06	80.92	61.56	55.29	– 19.02	33.48	46.25
B3LYP-D* (0 K)	163.94	89.65	63.75	59.43	– 20.60	37.35	50.09
B3LYP-D* (298 K)	156.47	86.03	59.85	55.13	– 18.86	35.24	48.31
Experimental (298 K)a	146.3	85.3	59.7	50.8	– 20.8	34.0	43.3
Experimental (298 K)b	149.4	85.2	57.9	53.5	– 20.0	34.1	45.8

a[26], b[27], c[10], d[9], e[11], calculated using results from Nelson et al.[10].

The elastic constants C11, C33, C44 and C66 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 and piezo-optic π 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 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 p41 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 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 p14 tensor component have a dependence on the wavelength λ.

Figure 2

The independent Pöckels’ tensor components of calcite calculated at finite wavelengths λ at the DFT/PBE-D2 level of theory, (a) p11, (b) p12 and p13, (c) p31 and p13 and (d) p44, p14, p41. The lines connecting the points are meant as a guide for the eye, whereas the dashed vertical line is the λ = 514 nm used in the experiments of Pöckels[9] and Nelson et al.[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 between the calculated room temperature and 0 K photoelastic constants as a function of λ.

In detail, there is a strong variation of the strain-optical components for 300 nm < λ < 500 nm, then the p tensor components asymptotically reach the value calculated at the λ = ∞ limit, with a trend of the form within the 300 nm–1000 nm range. The p41 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 p41 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 ) the choice of the reference system can influence the sign of shifting and rotating piezo-optic components.

Most of the p components show a variation between the values calculated at room and absolute zero temperatures, Δp = p(298 K) – p(0 K) (Fig. 2c,d), which asymptotically reaches zero as a function of λ. The calculated Δp are in absolute terms higher for λ < 500 nm. The exceptions to this trend are given by the p14 and p44 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 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, , 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-42]. It is a self-consistent, perturbative method that describes the relaxation of the crystalline orbitals when affected by an external electric field . 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 and are the unperturbed and perturbed Hamiltonian operators, respectively, the latter being the sum of the static electric field, (frequency ω = 0, wavelength λ = ∞), and the frequency-dependent one, whose amplitude and frequency of oscillation are and ω, respectively. is the electron position vector in real space. However, the scalar potential operator 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 indicating a reciprocal space vector, E(t) and being the b component of the applied electric field and of the gradient vector, respectively. The total energy of the system, Etot, 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 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-46] all the information related to CPHF/KS approach for periodic systems, and the detail of its implementation in the CRYSTAL code[39-42], which allows calculating the electronic contribution to the dielectric tensor for both static (frequency ω = 0, wavelength λ = ∞) and frequency-dependent (ω ≠ 0, λ ≠ ∞) electric fields.

Elastic and photoelastic constants

The elastic moduli C and photoelastic p 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, 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 C66 = (C11 – C12)/2. Similarly, we can define the photoelastic (Pöckels’) tensor p as:where p66 = (p11 – p12)/2. By previous definition, the piezo-optic tensor π has the same structure as that of the photoelastic components:

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 π 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-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.