- Email: [email protected]
Isotope effect on diamond's elastic-stiffness coefficients: An ab initio study
Computational Condensed Matter xxx (2017) 1e9
Contents lists available at ScienceDirect
Computational Condensed Matter journal homepage: http://ees.elsevier.com/cocom/default.asp
Isotope effect on diamond's elastic-stiffness coefficients: An ab initio study Daryl G. Clerc a, *, Hassel Ledbetter b a b
ArTek Product Development, Inc., Red Bud, IL 62278, United States Mechanical Engineering, University of Colorado, Boulder, CO 80309, United States
a r t i c l e i n f o
a b s t r a c t
Article history: Received 19 October 2016 Received in revised form 6 March 2017 Accepted 8 March 2017 Available online xxx
The effects of isotopic substitution on diamond's elastic-stiffness coefficients are studied theoretically by analyzing the zero-point motion and anharmonicity associated with lattice vibrations. Coefficients c11 ; c12 ; c44 ; and bulk modulus B are reported as purely theoretical functions of x, where x denotes the atomic fraction of 13C in 12C13 (1x)Cx. Second-order and third-order force constants are computed at the ab initio level and used as input to these expressions. As x increases, the predicted values of c11 , c12 , c44, and B undergo essentially linear increases: c11 ðxÞ ¼ c11 ð0Þð1 þ 0:000049xÞ; c12 ðxÞ ¼ c12 ð0Þð1 þ 0:00025xÞ; c44 ðxÞ ¼ c44 ð0Þð1 þ 0:000021xÞ and BðxÞ ¼ Bð0Þð1 þ 0:000088xÞ: Thus, compared to the values at x ¼ 0, the values of c11 , c12 , c44 ; and B are predicted to change by only 0.0049%, 0.025%, 0.0021%, and 0.0088%, respectively, at x ¼ 1. Our calculations also resolve a large discrepancy between two reported measurements of c12 , and provide a general method that can be used for arbitrary crystals having diamond's space group. © 2017 Elsevier B.V. All rights reserved.
Keywords: Diamond Elastic-stiffness coefficients Isotope effect Ab initio calculations
1. Introduction A material's elastic and anelastic properties enter quintessentially into a material's equation of state, which encompasses an enormous range of thermophysical properties. Diamond's elastic constants assume special importance because they relate directly to diamond's extreme properties: high hardness, high thermal conductivity, high resistance to extension and shear, low thermal expansivity. Less directly, elastic constants relate to diamond's other physical properties: low friction, wide optical transparency, low dielectric constant, and others. The recent discovery of superconductivity in boron-doped diamond [1] emphasizes further the importance of elastic constants. Specifically, elastic constants determine accurately the Debye temperature, which figures prominently in BCS-theory superconductivity and relates to an enormous variety of mechanicaldphysicaldthermal properties. This study has three principal purposes. First, extend our previous studies of second-order and third-order elastic properties of diamond [2]. Second, build on research of Vogelgesang and coworkers in which a quantum-mechanical analysis of zero-point
* Corresponding author. E-mail address: [email protected] (D.G. Clerc).
motion and lattice anharmonicity yielded the isotopic dependence of the bulk modulus [3], and extend their results by deriving expressions for c11 ðxÞ and c12 ðxÞ: Third, clarify discrepancies between theory and measurement. Vogelgesang and coworkers predicted a very small isotope effect in the bulk modulus: only a 0.1% increase at x ¼ 1 compared with the value at x ¼ 0. That result is quite inconsistent with an ultrasonic measurement of c12(x) by Hurley and colleagues in which c12 at x ¼ 0.99 was found to increase by 87% compared to the value at x ¼ 0, inferring a 17% increase in the bulk modulus at x ¼ 1 compared to x ¼ 0 [4]. This discrepancy was mentioned by Plekhanov in an extensive review of isotope effects [5]. 2. Computational method Calculations used the ab initio implementation of densityfunctional theory embodied in GAPSS (Gaussian Approach to Polymers, Surfaces, and Solids). An overview of the GAPSS program, along with descriptions of the basis sets, computational parameters and methods for establishing convergence of computed physical properties, among other details, were reported in our previous study of diamond [2]. In this study, calculations use two distortions types labeled I and III, where the labels are consistent with the previous study [2]. The
http://dx.doi.org/10.1016/j.cocom.2017.03.003 2352-2143/© 2017 Elsevier B.V. All rights reserved.
Please cite this article in press as: D.G. Clerc, H. Ledbetter, Isotope effect on diamond's elastic-stiffness coefficients: An ab initio study, Computational Condensed Matter (2017), http://dx.doi.org/10.1016/j.cocom.2017.03.003
D.G. Clerc, H. Ledbetter / Computational Condensed Matter xxx (2017) 1e9
and
1 2; 1 2; 0
=
[7]. The lattice parameters a, b, and c measured at
equilibrium (zero strain) are related according to ao ¼ bo ¼ co and a ¼ b ¼ g ¼ 60o . The cell contains two carbon atoms C1 and C2 located at ð0; 0; 0Þ and 1 4; 1 4; 1 4 [7]. The observed bonding =
=
=
corresponding Eulerian tensor elements are ½ε; ε; ε; 0; 0; 0 and ½ε; ε; ε; 0; 0; 0, using the six-component Voigt notation: ½ε1 ; ε2 ; ε3 ; ε4 ; ε5 ; ε6 . Elastic coefficients and force constants are converted from an energy basis to a pressure basis by dividing by the volume of the rhombohedral unit cell (described below in the Theory section). All elastic properties are calculated in the Eulerian framework (infinitesimal strain). For comparison, the elastic properties associated with distortion I are also calculated in the Lagrangian framework (finite strain). The connection between the two is described in Appendix A. Physical properties and statistical error were calculated as follows. Equilibrium total energy, lattice parameter, and second-order force constant were determined by fitting the computed total energies to Hooke's law over the harmonic region. The extent of the harmonic region was assigned at the strains beyond which the total energy deviated from Hooke's law. Force constants were then determined by fitting all energies to a third-order expression while holding equilibrium energy and lattice parameter fixed at their harmonic values, allowing only the second-order and third-order force constants to vary in a two-parameter fit. Statistical errors in the reported force constants were evaluated using five sets of computed data points. One set consisted of the force constant calculated using N points, where N ¼ 32 (distortion I) and N ¼ 27 (distortion III). The four remaining sets consist of N-8 data points. In each set, eight unique combinations of points are removed from the anharmonic region: four from the compression side (negative strain) and four from the expansion side (positive strain). Reported properties are arithmetic averages of the five sets; reported statistical error is the root-mean-square deviation from the mean for each distortion type. Effects of internal relaxation were found to be negligible. In some cases, internal relaxation introduces additional degrees of freedom because symmetry departs from tetrahedral. For distortion III, the four bond lengths remain equal (particular to the value of strain), and the bond angles depart from the tetrahedral value (109.5 ). Therefore, to quantify effects of internal relaxation, the total energy was minimized for distortion III, with respect to the three internal coordinates of the interior carbon atom C2, for zero strain and for strain near the harmonic limit ( ±0.02). Fractional coordinates at the minimum were within 0.2% of the original values, and the total energy at the minimum decreased by less than 5 105 au, only slightly above the level of numerical precision (2 105 au). Consistent with these findings, changes in the nearest-neighbor bond length and changes in the bond angles, compared to the original values, were very small, within 0.09% and 0.8%, respectively. These results indicate that the level of internal relaxation is small, consistent with the high degree of point symmetry retained by C2 during distortion III. This finding agrees with the prediction that diamond's internal-strain parameter is especially small, much smaller than other group-IV elements: Si, Ge, and a-Sn [6].
=
2
symmetry between C2 and its four nearest neighbors (the four vertices of the asymmetric subunit) is tetrahedral with bond length pffiffiffi 3ao =4. The equilibrium volume of the rhombohedral cell is pffiffiffi Vo ¼ 16d3o =3 3; which is one fourth the volume of the f.c.c. super cell. The labeling convention used to define strain in the rhombohedral cell is depicted by the two smaller diagrams in Fig. 1 and described in Appendix B. The change in a material's electronic energy relative to the equilibrium (unstrained) configuration is denoted by DU, defined by the equation
DU ¼
1X 1X cij εi εj þ cijk εi εj εk þ …; 2 6 i;j
(1)
i;j;k
in which terms are shown through third order. Here, DU is the change in total electronic energy (kinetic, Coulomb, and exchangecorrelation) per unit volume at zero-temperature. Parameters εi ; εj ; and εk denote strain in fractional units, cij and cijk denote the second- and third-order elastic coefficients, and indices fi; j; kg take the values one through six. In diamond, symmetry dictates three independent cij ðc11 ; c12 ; c44 Þ and six independent cijk ðc111 ; c112 ; c123 ; c144 ; c166 ; c456 Þ. For distortions that retain the orthogonality of the f.c.c. unit cell, that is ε4 ¼ ε5 ¼ ε6 ¼ 0, Eq. (1) takes the following form:
DU ¼
i 1h 2 c11 ε1 þ ε22 þ ε23 þ 2c12 ðε1 ε2 þ ε1 ε3 þ ε2 ε3 Þ 2 1h þ c111 ε31 þ ε32 þ ε33 þ 3c112 ε21 ε2 þ ε21 ε3 þ ε1 ε22 6 i þ ε3 ε22 þ ε23 ε1 þ ε23 ε2 þ 6c123 ε1 ε2 ε3 ;
(2)
where fourth- and higher-order terms are omitted. For distortion I (uniform dilation), Eq. (2) simplifies as follows:
DU
ðI Þ
¼
1 1 ð3c11 þ 6c12 Þε21 þ ð3c111 þ 18c112 þ 6c123 Þε31 ; 2 6 (3a)
kI ¼ 3c11 þ 6c12 ;
(3b)
3. Theory
gI ¼ 3c111 þ 18c112 þ 6c123 ;
(3c)
To quantify the effect of isotopes on the elastic constants of a material, it is necessary to go beyond an ordinary classical treatment to a quantum-mechanical approach that accounts for zeropoint motion of the nuclei. In this approach, the vibrational energy can be related to the isotopic mass and to the bond force constants that, in turn, can be related to the elastic coefficients. Fig. 1 shows the unit cell used in our calculations. That unit cell is the asymmetric rhombohedral subunit of the standard diamond f.c.c. unit cell (space group Fd3m, No. 227) whose vertices are at Cartesian fractional coordinates ð0; 0; 0Þ; 1 2; 0; 1 2 ; 0; 1 2; 1 2 ,
thus yielding definitions of the corresponding second-order and third-order force constants kI and gI , respectively. For distortion III, the analysis in Appendix C shows that DU can be expressed solely in terms of ε1 , akin to Eq. (3a). Specifically, the four nearest-neighbor bond lengths are equal at each value of strain, and the contribution of changes in bond angles is negligible compared to the contribution of changes in bond length. As a result, the angular dependence in DU can be omitted, thus yielding the following expressions for change in total electronic energy and associated force constants of distortion III.
=
=
=
=
Please cite this article in press as: D.G. Clerc, H. Ledbetter, Isotope effect on diamond's elastic-stiffness coefficients: An ab initio study, Computational Condensed Matter (2017), http://dx.doi.org/10.1016/j.cocom.2017.03.003
D.G. Clerc, H. Ledbetter / Computational Condensed Matter xxx (2017) 1e9
DU
ðIIIÞ
¼
1 1 ð3c11 2c12 Þε21 þ ðc111 þ 6c112 6c123 Þε31 ; 2 6
(4a)
kIII ¼ 3c11 2c12 ;
(4b)
gIII ¼ c111 þ 6c112 6c123 :
(4c)
Eqs. (3a) and (4a) may be written in terms of changes in total electronic energy and changes in nearest-neighbor distance, as follows. Distortions that preserve the orthogonality of the FCC unit cell yield εi ≡Da=ao ¼ Dd=do ≡u=do . Here, do denotes the equilibrium nearest-neighbor distance and the notation u≡Dd is adopted for consistency with the Brillouin and Raman study of isotopically substituted diamond reported by Vogelgesang et al. [3] Using DU ¼ Vo DU, the following expressions result for distortions I and III:
DU
ðI Þ
DU ðIIIÞ ¼
can be approximated by an expression analogous to Eq. (8) because of negligible angular contributions, as discussed in Appendix C, and set equal to the second-order coefficient in Eq. (5b), that is 8d poffiffiffi ð3c11 2c12 Þ. The result is the following pair of equations: 3 3
3c11 þ 6c12 ¼
pffiffiffi ZgI2 3 3kI ; 1 qffiffiffiffiffiffiffiffiffiffiffiffi 4do 8 2k5I M
(10a)
3c11 2c12 ¼
pffiffiffi 2 ZgIII 3 3kIII ffi 1 qffiffiffiffiffiffiffiffiffiffiffiffiffi 4do 8 2k5III M
(10b)
Solving Eqs. (10a) and (10b) for c11, in an isotopically pure initial state ð12 C; x ¼ 0Þ and in an isotopically substituted state (12 C and 13 C,
0 x 1), yields the following two expressions:
pffiffiffi " Z 3 c11 ðxÞ ¼ ðkI þ 3kIII Þ pffiffiffi 16do 8 2
8do 8 ¼ pffiffiffi ð3c11 þ 6c12 Þu2 þ pffiffiffi ð3c111 þ 18c112 3 3 9 3 þ 6c123 Þu3 ;
(5a)
8do 8 pffiffiffi ð3c11 2c12 Þu2 þ pffiffiffi ðc111 þ 6c112 6c123 Þu3 : 3 3 9 3 (5b)
The effects of isotopic substitution may be incorporated using the analysis of Vogelgesang et al. [3] After reversing the sign of the third-order coefficient to account for the opposite sign convention of this study, their Hamiltonian for the asymmetric unit is as follows:
b2 b ðIÞ ¼ P þ 2kI u2 þ 2gI u3 : Н 3 M
(6)
the vibrating pair of atoms at the origin and at (1/4,1/4,1/4). In that analysis, the expectation value of the Hamiltonian was variationally minimized with respect to a parameter in their wavefunction to yield the following expression:
(8)
Here the first term denotes the zero-point energy, d a displacement that is dimensionally equivalent to u; and K a renormalized stiffness constant defined by:
!
ZgI2 KI ¼ kI 1 qffiffiffiffiffiffiffiffiffiffiffiffi : 8 2k5I M
c11 ð0Þ ¼
gI2 3=2
kI
pffiffiffi " Z 3 ðkI þ 3kIII Þ pffiffiffi 16do 8 2
þ
2 3gIII
!
3=2
kIII
gI2
2 3gIII þ 3=2 3=2 kI kIII
!
# 1 pffiffiffiffiffiffiffi ; Mx
1 pffiffiffiffiffiffiffiffiffi M12
(11a)
#
(11b) Here the isotopic mass is given by Mx ¼ xM13 þ ð1 xÞM12, where M12 and M13 denote the atomic masses of 12 C and 13 C. Subtracting Eq. (11b) from Eq. (11a) yields the following expression for the x-dependence of c11 :
pffiffiffi 3Z pffiffiffi c11 ðxÞ ¼ c11 ð0Þ þ 128 2do
gI2 3=2
kI
þ
2 3gIII
!
3=2
kIII
sffiffiffiffiffiffiffiffiffi # " 1 M12 pffiffiffiffiffiffiffiffiffi 1 Mx M12 (11c)
Here u denotes displacement associated with uniform expansion. M denotes atomic mass, d the change in the nearest-neighbor distance, k the force constant between the central atom and each of its four nearest neighbors, g the lowest-order anharmonic contrib the momentum associated with the kinetic energy of bution, and P
rffiffiffiffiffiffiffi 2kI 2 EðdÞ ¼ Z þ 2KI ðd do Þ þ … : M
3
Similarly solving Eqs. (10a) and (10b) for c12 yields the following two expressions:
c12 ðxÞ ¼
pffiffiffi " 3 3 Z ðkI kIII Þ pffiffiffi 32do 8 2
pffiffiffi " 3 3 Z ðkI kIII Þ pffiffiffi c12 ð0Þ ¼ 32do 8 2
gI2
2 gIII 3=2 3=2 kI kIII
gI2 3=2
kI
3 3
2 gIII 3=2
kIII
!
# 1 pffiffiffiffiffiffiffi ; Mx
1 pffiffiffiffiffiffiffiffiffi M12
(12a)
# (12b)
Subtracting Eq. (12b) from Eq. (12a) yields the following expression for c12 ðxÞ:
pffiffiffi 3 3Z pffiffiffi c12 ðxÞ ¼ c12 ð0Þ þ 256 2do
gI2 3=2
kI
2 gIII 3=2
kIII
!
sffiffiffiffiffiffiffiffiffi # " 1 M12 pffiffiffiffiffiffiffiffiffi 1 Mx M12 (12c)
(9)
Eq. (9) reduces to KI ¼ kI in the absence of third-order effects ðgI ¼ 0Þ. With third-order contributions, M1=2 dependence follows. Expressions that connect the elastic coefficients in Eqs. (5a) and (5b) to isotopic mass may be derived, as follows. Given that Eqs. (5a) and (8) each pertain to distortion I that preserves the orthogonality of the unit cell, the second-order coefficient 8d poffiffiffi ð3c11 þ 6c12 Þ in Eq. (5a) is equivalent to the second-order co-
!
An expression for the x-dependence of the bulk modulus follows from B ¼ ðc11 þ 2c12 Þ=3:
# pffiffiffi " ZgI2 1 3 pffiffiffiffiffiffiffi ; BðxÞ ¼ k pffiffiffi 12do I 8 2 k3=2 Mx I
(13a)
# pffiffiffi " ZgI2 1 3 pffiffiffiffiffiffiffiffiffi ; k pffiffiffi Bð0Þ ¼ 12do I 8 2 k3=2 M12 I
(13b)
efficient 2KI in Eq. (8). Similarly, a force constant for distortion III Please cite this article in press as: D.G. Clerc, H. Ledbetter, Isotope effect on diamond's elastic-stiffness coefficients: An ab initio study, Computational Condensed Matter (2017), http://dx.doi.org/10.1016/j.cocom.2017.03.003
4
D.G. Clerc, H. Ledbetter / Computational Condensed Matter xxx (2017) 1e9
Fig. 1. Lattice parameters of diamond's rhombohedral unit cell are depicted within the outlined f.c.c. unit cell. The asymmetric unit consists of atoms at the origin and the central (shaded) atom at (1/4, 1/4, 1/4). For simplicity, only the central atom and its four nearest neighbors are shown. The equilibrium geometry is a ¼ b ¼ c≡ao and a ¼ b ¼ g ¼ 60o . The Eulerian coordinate system is also indicated (axes 1, 2, and 3). The two smaller diagrams on the right-hand side depict the strains associated with distortions I and III.
pffiffiffi 3Z BðxÞ ¼ Bð0Þ þ pffiffiffi 96 2do
gI2 3=2
kI
!
sffiffiffiffiffiffiffiffiffi # " 1 M12 pffiffiffiffiffiffiffiffiffi 1 Mx M12
(13c)
The x-dependence of the lattice parameter from the Vogelgesang et al. study [3] is included here as an ancillary note:
ZgI ; aðxÞ ¼ a∞ qffiffiffiffiffiffiffiffiffiffiffiffiffiffi 6k3I Mx
(14a)
ZgI ffi; að0Þ ¼ a∞ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 6k3I M12
(14b)
sffiffiffiffiffiffiffiffiffi # " ZgI 1 M12 aðxÞ ¼ að0Þ þ pffiffiffi 3=2 pffiffiffiffiffiffiffiffiffi 1 Mx M12 6kI
(14c)
Here a∞ denotes the lattice parameter corresponding to infinite atomic mass (negligible zero-point motion). 4. Results and discussion Table 1 shows our principal results: lattice parameter, bulk modulus, second-order force constants kI and kIII , and third-order force constants gI and gIII . Computed values of the pair kI and gI , and the pair kIII and gIII , were obtained by fitting the computed data points to the following equations:
DU
ðI Þ
DU
ðIIIÞ
¼
¼
1 2 1 3 kε þ gε ; 2 I 1 6 I 1 1 1 k ε2 þ g ε3 ; 2 III 1 6 III 1
(15a)
(15b)
obtained respectively by substituting Eqs. (3b) and (3c) into Eq. (3a), and by substituting Eqs. (4b) and (4c) into Eq. (4a). Total energy curves are shown as functions of strain for distortions I and III in Figs. 2(A) and 3(A), respectively. Linear fits of total energy to the cube of strain are shown in Figs. 2(B) and 3(B). Our principal results for lattice parameter and second-order properties are consistent with measured values. The calculated lattice parameter, 0.35545 nm, is quite accurate, only 0.3% below the measured value of 0.35665 nm [8] and only 0.07% different
from our previously calculated value of 0.35569 nm [2]. The bulk modulus, B ¼ kI =9 ¼ 447 GPa differs from the measured value 443 GPa [9] by only 0.9%. The computed force constants kI ¼ 4026 GPa and kIII ¼ 2783 GPa are consistent with values obtained by substituting measured values of c11 and c12 [9] into Eqs. (3b) and (4b) yielding 3981 GPa and 2989 GPa. Comparison shows only a small difference for kI (1.1%), but more substantial difference for kIII (-6.9%). Our previously reported value of gI , 55,020 (3520) GPa [2], differs markedly from our current value, 36,195 (1692) GPa, because of the extent of anharmonic strain used in the computations. The previous value was computed using only eight data points whose maximum strain is 2.9%. In contrast, the current value was computed using 32 data points across which the maximum strain doubles to 5.8%. Fitting our current data to Eq. (3a) using the previous narrower range yielded larger values of gI , similar to the previous result. However, when additional data points farther into the anharmonic region are included, the magnitude of gI decreases as does the statistical-fitting error. The high quality of the fit across the anharmonic region, shown in Fig. 2(A) and 2(B), indicates that our current result for gI is more accurate. Table 1 also compares gI and gIII to measured values. These measured values are calculated from third-order elastic coefficients (TOECs) reported in studies based on Raman spectroscopy [10,11] and shock wave compression [12]. Compared to our result, values of gI from Raman spectroscopy are far larger in magnitude, by ~60%. In contrast, the value of gI from shock wave compression is ~33% larger. With respect to gIII , the Raman values of 20,492 GPa and 26,423 GPa are 47% to 90% larger than our result, 13,925 GPa. However, the shock compression value of 14,170 GPa is very similar to our gIII , differing by only 1.8%. The above analysis also yielded two interesting observations. First, fitting the total energy to Lagrangian strain as defined in Eq. (A.5a), instead of Eulerian strain, yields gI ¼ 48,087 (1884) GPa, a result quite consistent with shock compression: gI ¼ 47,910 GPa, differing by only 0.4%. Second, using the narrow 2.9% range of strain produces larger magnitude values of gI (~55,000 GPa) consistent with the Raman results, whereas using the 5.8% range yields values more consistent with shock compression. The effect of isotopic substitution may be inferred qualitatively from Eqs. (11c), (12c), (13c), and (14c). The value of c11 predicted by Eq. (11c) depends on the force constants through a positive quan3=2
tity, gI2 =kI
3=2
2 =k þ3gIII III , thereby indicating that c11 increases with
Please cite this article in press as: D.G. Clerc, H. Ledbetter, Isotope effect on diamond's elastic-stiffness coefficients: An ab initio study, Computational Condensed Matter (2017), http://dx.doi.org/10.1016/j.cocom.2017.03.003
D.G. Clerc, H. Ledbetter / Computational Condensed Matter xxx (2017) 1e9
5
Fig. 2. The computed total energy versus strain is shown for distortion type I. (A) The solid line fitting computed data points (diamonds) was generated using Eq. (3a) with computed values of the lattice parameter kI , and gI reported in Table 1. The dashed line is the Hooke-law result generated by fitting only the second-order term. Pairs of dashed vertical lines indicate the computed harmonic limit. (B) Third-order contributions to the total energy in Eq. (3a) are plotted versus the cube of strain. Open diamonds represent expansion (positive strain) and filled diamonds represent contraction (negative strain). The line represents the best fit to the data, using the value of gI reported in Table 1.
increasing x. In contrast, the value of c12 from Eq. (13c) depends 3=2
3=2
2 =k upon gI2 =kI gIII III , whose sign depends upon the relative values of the force constants. For diamond, the sign is positive, and thus c12 increases with increasing x. The value of B and a given by 3=2
3=2
Eqs. (13c) and (14c) depend respectively upon gI2 =kI and gI =kI , which are positive and negative quantities, respectively, thereby indicating that B increases with increasing x and a decreases with increasing x. Concerning the third independent elastic-stiffness coefficient c44, one can deduce its value by assuming that the Zener elastic anisotropy, 2c44 =ðc11 c12 Þ; is invariant for isotopic substitution. Using Eqs. (11c) and (12c) we obtain the following expression for the cubic-face-plane shear modulus:
pffiffiffi c44 ð0Þ 3Z pffiffiffi 256 2do c11 ð0Þ c12 ð0Þ sffiffiffiffiffiffiffiffiffi # ! " 2 9gIII gI2 1 M12 3=2 pffiffiffiffiffiffiffiffiffi 1 3=2 Mx M kIII kI 12
c44 ðxÞ ¼ c44 ð0Þ þ
(16)
Isotopic substitution was found to induce very small changes in the elastic coefficients, the bulk modulus and the lattice parameter. Using Eqs. (11c), (12c), (13c), (14c), and (16), with computed values of c11 ð0Þ; c12 ð0Þ; Bð0Þ; að0Þ; and c44 ð0Þ; from Ref. [2], which are 1043 GPa, 128 GPa, 433 GPa, 0.35569 nm, and 534 GPa, respectively, fractional changes associated with full substitution 13C relative to 12 C have the following values:
½c11 ð1Þ c11 ð0Þ =c11 ð0Þ ¼ 4:9 105 ;
(17a)
½c12 ð1Þ c12 ð0Þ =c12 ð0Þ ¼ 2:5 104 ;
(17b)
½c44 ð1Þ c44 ð0Þ =c44 ð0Þ ¼ 2:1 105 ;
(17c)
½Bð1Þ Bð0Þ =Bð0Þ ¼ 8:8 105 ;
(17d)
½að1Þ að0Þ =að0Þ ¼ 1:5 105 :
(17e)
The magnitudes of the above changes are very small: varying from only 0.0015% to 0.025%. Essentially the same result occurs
Please cite this article in press as: D.G. Clerc, H. Ledbetter, Isotope effect on diamond's elastic-stiffness coefficients: An ab initio study, Computational Condensed Matter (2017), http://dx.doi.org/10.1016/j.cocom.2017.03.003
6
D.G. Clerc, H. Ledbetter / Computational Condensed Matter xxx (2017) 1e9
Fig. 3. The computed total energy versus strain is shown for distortion type III. (A) The solid line fitting computed data points (diamonds) was generated using Eq. (4a) with computed values of the lattice parameter kIII , and gIII reported in Table 1. The dashed line is the Hooke-law result generated by fitting only the second-order term. Pairs of dashed vertical lines indicate the computed harmonic limit. (B) Third-order contributions to the total energy in Eq. (4a) are plotted versus the cube of strain. Open diamonds represent expansion (positive strain) and filled diamonds represent contraction (negative strain). The line represents the best fit to the data, using the value of gIII reported in Table 1.
1=2
after accounting for the small nonlinearity (Mx dependence). For example, values of ½c11 ðxÞ c11 ð0Þ=c11 ð0Þ at x ¼ 0, x ¼ 1 2, and x ¼ 1 are respectively 0%, 0.0025% and 0.0049%. These three values average 0.0025%, which is comparable to the 0.0049% reported for Eq. (17a). We find that the Cauchy discrepancy c12/c44 stays the same within calculation error. Often taken as a clue to interatomicbonding nature, departure from central forces, diamond’s negative departure from unity is the largest for any cubic element and often called a “negative Cauchy pressure” [13]. Our predicted 0.0088% increase in B is comparable to the 0.18% increase reported by Vogelgesang et al. [3], but it is inconsistent with a measured increase of 17% reported by Hurley et al. [4]. The difference originates in the value of c12 that was used to calculate B. Specifically, the measured 87% increase in c12 at x ¼ 1 is enormously larger than our predicted 0.025% increase. Our result for c12, a very small change at x ¼ 1 of 0.025%, is consistent with a simple theory based on Heisenberg's theorem. Gilman, in describing diamond's extreme hardness, derived a basic bulk-modulus relationship B ¼ 3e2 =40pro4 , where e denotes the valence charge of the atom and ro the equilibrium radius of the valence electron cloud around a single atom [14]. Thus, the physical
property that governs the bulk modulus is atomic volume, not atomic mass, and efforts to increase the bulk modulus, which reflects mechanical hardness, may focus on reducing ro . The additional neutron that exists in 13C, being uncharged, does not significantly affect the electronic charge distribution and hence does not significantly alter ro. As a result, we expect the corresponding elastic constants and bulk modulus to be affected little by 13 C12C substitution. 5. Conclusions The effects of isotopic substitution of 13C on the elastic constants of diamond were studied theoretically using predicted properties taken from our ab initio studies, and using an analysis of zero-point motion and the associated anharmonicity of lattice vibrations. We predict that changes in c11 , c12 , c44 , and B with respect to x are essentially linear: c11 ðxÞ ¼ c11 ð0Þð1 þ 0:000049xÞ; c12 ðxÞ ¼ c12 ð0Þð1 þ 0:00025xÞ; c44 ðxÞ ¼ c44 ð0Þð1 þ 0:000021xÞ; and BðxÞ ¼ Bð0Þð1 þ 0:000088xÞ Compared with values at x ¼ 0, c11 ; c12 ; c44 ; and B are predicted to change by only 0.0049%, 0.025%, 0.0021%, and 0.0088%, respectively, at x ¼ 1. Our predicted 0.0088% increase in B is comparable to the 0.18% increase reported by
=
Please cite this article in press as: D.G. Clerc, H. Ledbetter, Isotope effect on diamond's elastic-stiffness coefficients: An ab initio study, Computational Condensed Matter (2017), http://dx.doi.org/10.1016/j.cocom.2017.03.003
D.G. Clerc, H. Ledbetter / Computational Condensed Matter xxx (2017) 1e9
7
Table 1 Previously reported and present results. Quantity Lattice parameter Bulk modulus Second-order force constants
a (nm) B (GPa) kI (GPa) kIII (GPa) gI (GPa) gIII (GPa)
Third-order force constants a b c d e f g h i
Prior studies
This study
0.35665a, 0.35569b 443d 3981e 2989e 58,788g, 54,309h, 47,910i 20,492g, 26,423h, 14,170i
0.35545c 447 (3) 4026 (25)f 2783 (11)f 36,195 (1692) 13,925 (556)
Reference [8]. Reference [2]. The computed f.c.c. lattice parameter is the average for distortions I (0.35545 nm) and III (0.35563 nm). Reference [9]. Calculated using Eqs. (3b) and (4b), and the measured c11 and c12 from Ref. [9]. pffiffiffi Terms of Eq. (6), equivalent values of kI and kIII , are 1844 J/m2 and 1365 J/m2, respectively - obtained by multiplying the GPa values by 16do =3 3. Calculated using Eqs. (3c) and (4c), and three Raman-measured TOECs in Ref. [10]. Calculated using Eqs. (3c) and (4c), and three Raman-measured TOECs in Ref. [11]. Calculated using Eqs. (3c) and (4c), and the average of two sets of three shock-wave-measured TOECs in Ref. [12].
Vogelgesang et al. [3]. The difference between our predicted increase in B and the measured 17% increase in B reported by Hurley et al. [4] originates in the measured 87% increase in c12 at x ¼ 1, which is far larger than our predicted 0.025% increase. The method used in this study may also be used for other crystals having diamond's space group, and as such could lead to further studies within a similar theoretical framework. Acknowledgments This research received partial support from Los Alamos National Laboratory (Contract No. 96628-001-04 30). The authors appreciate the encouragement of this study by Dr. Albert Migliori (Los Alamos National Laboratory).
1 2
h1 ¼ h2 ¼ ε1 þ ε21 ;
III :
1 2
h3 ¼ ε1 þ ε21 ¼ h1 þ ε21 z h1 ;
orthogonality of the solid, the derivatives in Eq. (A.3) that possess cross terms vanish, yielding the following relationship:
hpp
1 ¼ 2
vup vxp
2
vup þ2 vxp
! (A.4)
The remaining terms can be related to Eulerian strain by using Eq. (A.1) with i ¼ j to yield εii ¼ vui =vxi and hence vup =vxp ¼ εpp : That result can be substituted into Eq. (A.4) to yield the connection between Lagrangian and Eulerian strain for distortions I and III, as follows:
I:
1 2
h1 ¼ h2 ¼ h3 ¼ ε1 þ ε21 ;
h4 ¼ h5 ¼ h6 ¼ 0;
h4 ¼ h5 ¼ h6 ¼ 0:
(A.5a)
(A.5b)
Appendix A. Eulerian and Lagrangian frameworks Strain εi relates directly to the strain denoted by εij using the familiar Voigt notation. In turn, εij is connected to the Eulerian or classical definition of strain by the following relationship:
εij ¼
1 vui vuj þ 2 vxj vxi
! (A.1)
Here u denotes infinitesimal displacement in the strained solid. Lagrangian strain hpq accounts for finite deformation by including quadratic strain effects:
hpq ¼
1 vXi vXi dpq 2 vxp vxq
(A.2)
Here, Xi ði ¼ 1; 2; 3Þ denotes position in the strained solid with implicit summation over index i, and xp and xq denote position in the unstrained solid, where p; q2f1; 2; 3g. By using displacement coordinates, ui ¼ Xi xi , Eq. (A.2) can be written as:
hpq ¼
1 vui vui vui vu þ diq þ i dip þ dip dip dpq 2 vxp vxq vxp vxq
(A.3)
Appendix B. Conventions used to define strain The smallest periodic unit in diamond is the rhombohedron, depicted in Fig. 1. The vertices of the rhombohedral cell are defined by four atoms: one at the origin plus three at each face center of the standard FCC supercell. The interior of the rhombohedral cell contains one additional atom at (1/4, 1/4, 1/4). The lattice vectors of the rhombohedral cell are defined as follows in a Cartesian basis:
! a ¼ 〈a cosðp=4Þ; 0; a sinðp=4Þ〉;
(B.1a)
! b ¼ 〈b cosðp=4Þ; b sinðp=4Þ; 0〉;
(B.1b)
! c ¼ 〈0; c sinðp=4Þ; c sinðp=4Þ〉:
(B.1c)
If we define b ¼ ao ð1 þ ε1 Þ and a ¼ ao ð1 þ 4ε1 Þ, where 4 is a constant, then Eqs. (B.1a)e(B.1c) simplify as follows:
ao ! a ¼ pffiffiffi h1 þ ε1 ; 0; 1 þ 4ε1 i; 2
(B.2a)
Because the two types of strain used in this study preserve the Please cite this article in press as: D.G. Clerc, H. Ledbetter, Isotope effect on diamond's elastic-stiffness coefficients: An ab initio study, Computational Condensed Matter (2017), http://dx.doi.org/10.1016/j.cocom.2017.03.003
8
D.G. Clerc, H. Ledbetter / Computational Condensed Matter xxx (2017) 1e9
! ao b ¼ pffiffiffi h1 þ ε1 ; 1 þ ε1 ; 0i; 2
(B.2b)
ao ! p c ¼ ffiffiffi h0; 1 þ ε1 ; 1 þ 4ε1 i 2
(B.2c)
Euclidian inner products and vector magnitudes have the following forms:
! ! ! ! ao a $ b ¼ b $ c ¼ ð1 þ ε1 Þ2 ; 2 2
(B.3a)
! ! ! d a¼ C 2 a ao ¼ pffiffiffi hð1þε1 Þðf1 þf2 1Þ;ð1þε1 Þðf2 þf3 Þ;ð1þ4ε1 Þðf1 þf3 1Þi; 2 (C.1b)
! ! ! d b¼ C 2 b ao ¼ pffiffiffi hð1þε1 Þðf1 þf2 1Þ;ð1þε1 Þðf2 þf3 1Þ;ð1þ4ε1 Þðf1 þf3 Þi; 2 (C.1c)
2
! ! ao a $ c ¼ ð1 þ 4ε1 Þ2 ; 2 ao ! ! j a j ¼ j c j ¼ pffiffiffi 2
(B.3b)
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð1 þ ε1 Þ2 þ ð1 þ 4ε1 Þ2 ;
! b ¼ ao ð1 þ ε1 Þ
(B.3c)
(B.3d)
the relations cosðaÞ ¼ b b$b c , cosðbÞ ¼ b a $b c , and b where the carat denotes unit magnitude, we arrive at cosðgÞ ¼ b a $ b, the following equations: Using
ð1 þ ε1 Þ ; cosðaÞ ¼ cosðgÞ ¼ pffiffiffiqffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ð1 þ ε1 Þ2 þ ð1 þ 4ε1 Þ2
cosðbÞ ¼
ð1 þ 4ε1 Þ
(B.4a)
2
ð1 þ ε1 Þ2 þ ð1 þ 4ε1 Þ2
:
! ! ! d c¼ C 2 c ao ¼ pffiffiffi hð1þε1 Þðf1 þf2 Þ;ð1þε1 Þðf2 þf3 1Þ;ð1þ4ε1 Þðf1 þf3 1Þi 2 (C.1d) Using f1 ¼ f2 ¼ f3 ¼ 1/4, Eqs. (C.1aed) simplify as follows:
! ao d O ¼ pffiffiffi h1 þ ε1 ; 1 þ ε1 ; 1 þ 4ε1 i; 2 2
(C.2a)
! ao d a ¼ pffiffiffi h 1 ε1 ; 1 þ ε1 ; 1 4ε1 i; 2 2
(C.2b)
! ao d b ¼ pffiffiffi h 1 ε1 ; 1 ε1 ; 1 þ 4ε1 i; 2 2
(C.2c)
! ao d c ¼ pffiffiffi h1 þ ε1 ; 1 ε1 ; 1 4ε1 i: 2 2
(C.2d)
(B.4b)
The three sets of preceding equations show that the geometry of diamond's rhombohedral unit cell, for a constrained distortion such as I or III, is defined completely by the single Eulerian strain coordinate ε1 . For distortion I ð4 ¼ 1Þ, Eqs. (B.3a) and (B.3b) reduce to ! ! ! j a j ¼ b ¼ j c j and Eqs. (B.4a) and (B.4b) yield a ¼ b ¼ g ¼ 60o . Relations for distortion III may be obtained through the substitution 4 ¼ 1.
Appendix C. Effects of internal relaxation
These four bonds have the same length, ! ! ! ! j d O j ¼ j d a j ¼ j d b j ¼ j d c j in distortion I ð4 ¼ 1Þ and in distortion III ð4 ¼ 1Þ. For either type of distortion, bond-stretching contributions to the potential energy can be written in terms of one variable ε1 , analogous to the use of the single variable in the potential energy shown in Eq. (3a). The angles between the central atom and the four nearest neighbors can be found analogous to qOa ¼ ; ! ! ! ! ! ! ð d O ; d a Þ ¼ arccosð d O $ d a = d O d a Þ; the following equations emerge:
This appendix shows that distortions I and III lead to equations for the internal coordinate which are consistent with the form of the potential energy used in Eq. (3a). This internal coordinate is the position of the shaded central atom C2 in Fig. 1, defined in term of ! ! ! ! the three primitive lattice vectors as C 2 ¼ f1 a þ f2 b þ f3 c ; where ff1 ; f2 ; f3 g are constants whose equilibrium values are f1 ¼ f2 ¼ f3 ¼ 1=4. Combining this definition with the definitions of the lattice vectors in Eqs. (B.2a), (B.2b), and (B.2c) yields the following vectors for the bonds connecting C2 to its four nearestneighbors:
! ! dO ¼ C2 ao ¼ pffiffiffi hð1 þ ε1 Þðf1 þ f2 Þ; ð1 þ ε1 Þðf2 þ f3 Þ; ð1 þ 4ε1 Þðf1 þ f3 Þi; 2 (C.1a)
cosðqOa Þ ¼ cosðqOc Þ ¼ cosðqab Þ ¼ cosðqbc Þ ð1 þ 4ε1 Þ2 ; ¼ ð1 þ 4ε1 Þ2 þ 2ð1 þ ε1 Þ2 cosðqOb Þ ¼ cosðqac Þ ¼
ð1 þ 4ε1 Þ2 2ð1 þ ε1 Þ2 ð1 þ 4ε1 Þ2 þ 2ð1 þ ε1 Þ2
(C.3a)
:
(C.3b)
For distortion I ð4 ¼ 1Þ, the six cosines are equal to 1/3 irrespective of the level of strain, corresponding to the tetrahedral angle of 109.47. For distortion III ð4 ¼ 1Þ, the six bond angles split into two sets given by Eqs. (C.3a) and (C.3b), with the angles being equal within each set, the particular value depending upon the level of strain. Thus, the change in potential energy associated with distortion III is a function of bond angle as well as bond length. Angular contributions to the potential energy were found to be negligible, as reported in the Computational Method section.
Please cite this article in press as: D.G. Clerc, H. Ledbetter, Isotope effect on diamond's elastic-stiffness coefficients: An ab initio study, Computational Condensed Matter (2017), http://dx.doi.org/10.1016/j.cocom.2017.03.003
D.G. Clerc, H. Ledbetter / Computational Condensed Matter xxx (2017) 1e9
References [1] E. Ekimov, V. Sidorov, E. Bauer, et al., Nature 428 (2004) 542e545. [2] D.G. Clerc, H. Ledbetter, J. Phys. Chem. Sol. 66 (2005) 1589e1597. [3] R. Vogelgesang, A.K. Ramdas, S. Rodriquez, M. Grimsditch, T.R. Anthony, Phys. Rev. B 54 (6) (1996) 3989e3999. Eq. (6) of this study, with opposite sign in the third-order coefficient, corresponds to Eq. (17). Similarly, Eqs. (14a) and (14b) of this study both correspond to Eq. (22), and Eq. (14c) corresponds to Eq. (29). The isotope effect on the bulk modulus (0.18%) is reported on page 3997 in Eq. (38). [4] D.C. Hurley, R.S. Gilmore, W.F. Banholzer, J. Appl. Phys. 76 (12) (1994) 7726e7730. See page 7730, paragraph 3. [5] V. Plekhanov, Mater. Sci. Eng. R. 35 (2001) 139e237. See page 153, paragraph 1. [6] W. Weber, Phys. Rev. B 15 (10) (1977) 4789e4803. See Table II on page 4793. [7] G. Burns, A. Glazer, Space Groups for Solid State Scientists, Academic Press,
[8] [9] [10] [11] [12] [13] [14]
9
Inc., San Diego, second ed., 1990. See page 152, Fig. 6e5(g) and page 157, Fig. 6e8. T. Sato, K. Ohashi, T. Sudoh, et al., Phys. Rev. B 65 (2002) 092102. See page 3, Fig. 2 (0 K values). H. McSkimin, P. Andreatch, P. Glynn, J. Appl. Phys. 43 (3) (1972) 985e987. See page 987, Table II. M.H. Grimsditch, E. Anastassakis, M. Cardona, Phys. Rev. B 18 (2) (1978) 901e904. See page 904, Table II. E. Anastassakis, A. Cantarero, M. Cardona, Phys. Rev. B 41 (11) (1990) 7529e7535. See page 7533, Table IV. J.M. Winey, A. Hmiel, Y.M. Gupta, J. Phys. Chem. Sol. 93 (2016) 118e120. See page 119, Table 2: sets A and B. H. Ledbetter, S. Kim (Eds.), Handbook of Elastic Properties of Solids, Liquids, and Gases, Academic, San Diego, 2001. J.J. Gilman, Mater. Sci. Eng 7 (1971) 357e361. See page 359, Eq. (20).
Please cite this article in press as: D.G. Clerc, H. Ledbetter, Isotope effect on diamond's elastic-stiffness coefficients: An ab initio study, Computational Condensed Matter (2017), http://dx.doi.org/10.1016/j.cocom.2017.03.003
Contents lists available at ScienceDirect
Computational Condensed Matter journal homepage: http://ees.elsevier.com/cocom/default.asp
Isotope effect on diamond's elastic-stiffness coefficients: An ab initio study Daryl G. Clerc a, *, Hassel Ledbetter b a b
ArTek Product Development, Inc., Red Bud, IL 62278, United States Mechanical Engineering, University of Colorado, Boulder, CO 80309, United States
a r t i c l e i n f o
a b s t r a c t
Article history: Received 19 October 2016 Received in revised form 6 March 2017 Accepted 8 March 2017 Available online xxx
The effects of isotopic substitution on diamond's elastic-stiffness coefficients are studied theoretically by analyzing the zero-point motion and anharmonicity associated with lattice vibrations. Coefficients c11 ; c12 ; c44 ; and bulk modulus B are reported as purely theoretical functions of x, where x denotes the atomic fraction of 13C in 12C13 (1x)Cx. Second-order and third-order force constants are computed at the ab initio level and used as input to these expressions. As x increases, the predicted values of c11 , c12 , c44, and B undergo essentially linear increases: c11 ðxÞ ¼ c11 ð0Þð1 þ 0:000049xÞ; c12 ðxÞ ¼ c12 ð0Þð1 þ 0:00025xÞ; c44 ðxÞ ¼ c44 ð0Þð1 þ 0:000021xÞ and BðxÞ ¼ Bð0Þð1 þ 0:000088xÞ: Thus, compared to the values at x ¼ 0, the values of c11 , c12 , c44 ; and B are predicted to change by only 0.0049%, 0.025%, 0.0021%, and 0.0088%, respectively, at x ¼ 1. Our calculations also resolve a large discrepancy between two reported measurements of c12 , and provide a general method that can be used for arbitrary crystals having diamond's space group. © 2017 Elsevier B.V. All rights reserved.
Keywords: Diamond Elastic-stiffness coefficients Isotope effect Ab initio calculations
1. Introduction A material's elastic and anelastic properties enter quintessentially into a material's equation of state, which encompasses an enormous range of thermophysical properties. Diamond's elastic constants assume special importance because they relate directly to diamond's extreme properties: high hardness, high thermal conductivity, high resistance to extension and shear, low thermal expansivity. Less directly, elastic constants relate to diamond's other physical properties: low friction, wide optical transparency, low dielectric constant, and others. The recent discovery of superconductivity in boron-doped diamond [1] emphasizes further the importance of elastic constants. Specifically, elastic constants determine accurately the Debye temperature, which figures prominently in BCS-theory superconductivity and relates to an enormous variety of mechanicaldphysicaldthermal properties. This study has three principal purposes. First, extend our previous studies of second-order and third-order elastic properties of diamond [2]. Second, build on research of Vogelgesang and coworkers in which a quantum-mechanical analysis of zero-point
* Corresponding author. E-mail address: [email protected] (D.G. Clerc).
motion and lattice anharmonicity yielded the isotopic dependence of the bulk modulus [3], and extend their results by deriving expressions for c11 ðxÞ and c12 ðxÞ: Third, clarify discrepancies between theory and measurement. Vogelgesang and coworkers predicted a very small isotope effect in the bulk modulus: only a 0.1% increase at x ¼ 1 compared with the value at x ¼ 0. That result is quite inconsistent with an ultrasonic measurement of c12(x) by Hurley and colleagues in which c12 at x ¼ 0.99 was found to increase by 87% compared to the value at x ¼ 0, inferring a 17% increase in the bulk modulus at x ¼ 1 compared to x ¼ 0 [4]. This discrepancy was mentioned by Plekhanov in an extensive review of isotope effects [5]. 2. Computational method Calculations used the ab initio implementation of densityfunctional theory embodied in GAPSS (Gaussian Approach to Polymers, Surfaces, and Solids). An overview of the GAPSS program, along with descriptions of the basis sets, computational parameters and methods for establishing convergence of computed physical properties, among other details, were reported in our previous study of diamond [2]. In this study, calculations use two distortions types labeled I and III, where the labels are consistent with the previous study [2]. The
http://dx.doi.org/10.1016/j.cocom.2017.03.003 2352-2143/© 2017 Elsevier B.V. All rights reserved.
Please cite this article in press as: D.G. Clerc, H. Ledbetter, Isotope effect on diamond's elastic-stiffness coefficients: An ab initio study, Computational Condensed Matter (2017), http://dx.doi.org/10.1016/j.cocom.2017.03.003
D.G. Clerc, H. Ledbetter / Computational Condensed Matter xxx (2017) 1e9
and
1 2; 1 2; 0
=
[7]. The lattice parameters a, b, and c measured at
equilibrium (zero strain) are related according to ao ¼ bo ¼ co and a ¼ b ¼ g ¼ 60o . The cell contains two carbon atoms C1 and C2 located at ð0; 0; 0Þ and 1 4; 1 4; 1 4 [7]. The observed bonding =
=
=
corresponding Eulerian tensor elements are ½ε; ε; ε; 0; 0; 0 and ½ε; ε; ε; 0; 0; 0, using the six-component Voigt notation: ½ε1 ; ε2 ; ε3 ; ε4 ; ε5 ; ε6 . Elastic coefficients and force constants are converted from an energy basis to a pressure basis by dividing by the volume of the rhombohedral unit cell (described below in the Theory section). All elastic properties are calculated in the Eulerian framework (infinitesimal strain). For comparison, the elastic properties associated with distortion I are also calculated in the Lagrangian framework (finite strain). The connection between the two is described in Appendix A. Physical properties and statistical error were calculated as follows. Equilibrium total energy, lattice parameter, and second-order force constant were determined by fitting the computed total energies to Hooke's law over the harmonic region. The extent of the harmonic region was assigned at the strains beyond which the total energy deviated from Hooke's law. Force constants were then determined by fitting all energies to a third-order expression while holding equilibrium energy and lattice parameter fixed at their harmonic values, allowing only the second-order and third-order force constants to vary in a two-parameter fit. Statistical errors in the reported force constants were evaluated using five sets of computed data points. One set consisted of the force constant calculated using N points, where N ¼ 32 (distortion I) and N ¼ 27 (distortion III). The four remaining sets consist of N-8 data points. In each set, eight unique combinations of points are removed from the anharmonic region: four from the compression side (negative strain) and four from the expansion side (positive strain). Reported properties are arithmetic averages of the five sets; reported statistical error is the root-mean-square deviation from the mean for each distortion type. Effects of internal relaxation were found to be negligible. In some cases, internal relaxation introduces additional degrees of freedom because symmetry departs from tetrahedral. For distortion III, the four bond lengths remain equal (particular to the value of strain), and the bond angles depart from the tetrahedral value (109.5 ). Therefore, to quantify effects of internal relaxation, the total energy was minimized for distortion III, with respect to the three internal coordinates of the interior carbon atom C2, for zero strain and for strain near the harmonic limit ( ±0.02). Fractional coordinates at the minimum were within 0.2% of the original values, and the total energy at the minimum decreased by less than 5 105 au, only slightly above the level of numerical precision (2 105 au). Consistent with these findings, changes in the nearest-neighbor bond length and changes in the bond angles, compared to the original values, were very small, within 0.09% and 0.8%, respectively. These results indicate that the level of internal relaxation is small, consistent with the high degree of point symmetry retained by C2 during distortion III. This finding agrees with the prediction that diamond's internal-strain parameter is especially small, much smaller than other group-IV elements: Si, Ge, and a-Sn [6].
=
2
symmetry between C2 and its four nearest neighbors (the four vertices of the asymmetric subunit) is tetrahedral with bond length pffiffiffi 3ao =4. The equilibrium volume of the rhombohedral cell is pffiffiffi Vo ¼ 16d3o =3 3; which is one fourth the volume of the f.c.c. super cell. The labeling convention used to define strain in the rhombohedral cell is depicted by the two smaller diagrams in Fig. 1 and described in Appendix B. The change in a material's electronic energy relative to the equilibrium (unstrained) configuration is denoted by DU, defined by the equation
DU ¼
1X 1X cij εi εj þ cijk εi εj εk þ …; 2 6 i;j
(1)
i;j;k
in which terms are shown through third order. Here, DU is the change in total electronic energy (kinetic, Coulomb, and exchangecorrelation) per unit volume at zero-temperature. Parameters εi ; εj ; and εk denote strain in fractional units, cij and cijk denote the second- and third-order elastic coefficients, and indices fi; j; kg take the values one through six. In diamond, symmetry dictates three independent cij ðc11 ; c12 ; c44 Þ and six independent cijk ðc111 ; c112 ; c123 ; c144 ; c166 ; c456 Þ. For distortions that retain the orthogonality of the f.c.c. unit cell, that is ε4 ¼ ε5 ¼ ε6 ¼ 0, Eq. (1) takes the following form:
DU ¼
i 1h 2 c11 ε1 þ ε22 þ ε23 þ 2c12 ðε1 ε2 þ ε1 ε3 þ ε2 ε3 Þ 2 1h þ c111 ε31 þ ε32 þ ε33 þ 3c112 ε21 ε2 þ ε21 ε3 þ ε1 ε22 6 i þ ε3 ε22 þ ε23 ε1 þ ε23 ε2 þ 6c123 ε1 ε2 ε3 ;
(2)
where fourth- and higher-order terms are omitted. For distortion I (uniform dilation), Eq. (2) simplifies as follows:
DU
ðI Þ
¼
1 1 ð3c11 þ 6c12 Þε21 þ ð3c111 þ 18c112 þ 6c123 Þε31 ; 2 6 (3a)
kI ¼ 3c11 þ 6c12 ;
(3b)
3. Theory
gI ¼ 3c111 þ 18c112 þ 6c123 ;
(3c)
To quantify the effect of isotopes on the elastic constants of a material, it is necessary to go beyond an ordinary classical treatment to a quantum-mechanical approach that accounts for zeropoint motion of the nuclei. In this approach, the vibrational energy can be related to the isotopic mass and to the bond force constants that, in turn, can be related to the elastic coefficients. Fig. 1 shows the unit cell used in our calculations. That unit cell is the asymmetric rhombohedral subunit of the standard diamond f.c.c. unit cell (space group Fd3m, No. 227) whose vertices are at Cartesian fractional coordinates ð0; 0; 0Þ; 1 2; 0; 1 2 ; 0; 1 2; 1 2 ,
thus yielding definitions of the corresponding second-order and third-order force constants kI and gI , respectively. For distortion III, the analysis in Appendix C shows that DU can be expressed solely in terms of ε1 , akin to Eq. (3a). Specifically, the four nearest-neighbor bond lengths are equal at each value of strain, and the contribution of changes in bond angles is negligible compared to the contribution of changes in bond length. As a result, the angular dependence in DU can be omitted, thus yielding the following expressions for change in total electronic energy and associated force constants of distortion III.
=
=
=
=
Please cite this article in press as: D.G. Clerc, H. Ledbetter, Isotope effect on diamond's elastic-stiffness coefficients: An ab initio study, Computational Condensed Matter (2017), http://dx.doi.org/10.1016/j.cocom.2017.03.003
D.G. Clerc, H. Ledbetter / Computational Condensed Matter xxx (2017) 1e9
DU
ðIIIÞ
¼
1 1 ð3c11 2c12 Þε21 þ ðc111 þ 6c112 6c123 Þε31 ; 2 6
(4a)
kIII ¼ 3c11 2c12 ;
(4b)
gIII ¼ c111 þ 6c112 6c123 :
(4c)
Eqs. (3a) and (4a) may be written in terms of changes in total electronic energy and changes in nearest-neighbor distance, as follows. Distortions that preserve the orthogonality of the FCC unit cell yield εi ≡Da=ao ¼ Dd=do ≡u=do . Here, do denotes the equilibrium nearest-neighbor distance and the notation u≡Dd is adopted for consistency with the Brillouin and Raman study of isotopically substituted diamond reported by Vogelgesang et al. [3] Using DU ¼ Vo DU, the following expressions result for distortions I and III:
DU
ðI Þ
DU ðIIIÞ ¼
can be approximated by an expression analogous to Eq. (8) because of negligible angular contributions, as discussed in Appendix C, and set equal to the second-order coefficient in Eq. (5b), that is 8d poffiffiffi ð3c11 2c12 Þ. The result is the following pair of equations: 3 3
3c11 þ 6c12 ¼
pffiffiffi ZgI2 3 3kI ; 1 qffiffiffiffiffiffiffiffiffiffiffiffi 4do 8 2k5I M
(10a)
3c11 2c12 ¼
pffiffiffi 2 ZgIII 3 3kIII ffi 1 qffiffiffiffiffiffiffiffiffiffiffiffiffi 4do 8 2k5III M
(10b)
Solving Eqs. (10a) and (10b) for c11, in an isotopically pure initial state ð12 C; x ¼ 0Þ and in an isotopically substituted state (12 C and 13 C,
0 x 1), yields the following two expressions:
pffiffiffi " Z 3 c11 ðxÞ ¼ ðkI þ 3kIII Þ pffiffiffi 16do 8 2
8do 8 ¼ pffiffiffi ð3c11 þ 6c12 Þu2 þ pffiffiffi ð3c111 þ 18c112 3 3 9 3 þ 6c123 Þu3 ;
(5a)
8do 8 pffiffiffi ð3c11 2c12 Þu2 þ pffiffiffi ðc111 þ 6c112 6c123 Þu3 : 3 3 9 3 (5b)
The effects of isotopic substitution may be incorporated using the analysis of Vogelgesang et al. [3] After reversing the sign of the third-order coefficient to account for the opposite sign convention of this study, their Hamiltonian for the asymmetric unit is as follows:
b2 b ðIÞ ¼ P þ 2kI u2 þ 2gI u3 : Н 3 M
(6)
the vibrating pair of atoms at the origin and at (1/4,1/4,1/4). In that analysis, the expectation value of the Hamiltonian was variationally minimized with respect to a parameter in their wavefunction to yield the following expression:
(8)
Here the first term denotes the zero-point energy, d a displacement that is dimensionally equivalent to u; and K a renormalized stiffness constant defined by:
!
ZgI2 KI ¼ kI 1 qffiffiffiffiffiffiffiffiffiffiffiffi : 8 2k5I M
c11 ð0Þ ¼
gI2 3=2
kI
pffiffiffi " Z 3 ðkI þ 3kIII Þ pffiffiffi 16do 8 2
þ
2 3gIII
!
3=2
kIII
gI2
2 3gIII þ 3=2 3=2 kI kIII
!
# 1 pffiffiffiffiffiffiffi ; Mx
1 pffiffiffiffiffiffiffiffiffi M12
(11a)
#
(11b) Here the isotopic mass is given by Mx ¼ xM13 þ ð1 xÞM12, where M12 and M13 denote the atomic masses of 12 C and 13 C. Subtracting Eq. (11b) from Eq. (11a) yields the following expression for the x-dependence of c11 :
pffiffiffi 3Z pffiffiffi c11 ðxÞ ¼ c11 ð0Þ þ 128 2do
gI2 3=2
kI
þ
2 3gIII
!
3=2
kIII
sffiffiffiffiffiffiffiffiffi # " 1 M12 pffiffiffiffiffiffiffiffiffi 1 Mx M12 (11c)
Here u denotes displacement associated with uniform expansion. M denotes atomic mass, d the change in the nearest-neighbor distance, k the force constant between the central atom and each of its four nearest neighbors, g the lowest-order anharmonic contrib the momentum associated with the kinetic energy of bution, and P
rffiffiffiffiffiffiffi 2kI 2 EðdÞ ¼ Z þ 2KI ðd do Þ þ … : M
3
Similarly solving Eqs. (10a) and (10b) for c12 yields the following two expressions:
c12 ðxÞ ¼
pffiffiffi " 3 3 Z ðkI kIII Þ pffiffiffi 32do 8 2
pffiffiffi " 3 3 Z ðkI kIII Þ pffiffiffi c12 ð0Þ ¼ 32do 8 2
gI2
2 gIII 3=2 3=2 kI kIII
gI2 3=2
kI
3 3
2 gIII 3=2
kIII
!
# 1 pffiffiffiffiffiffiffi ; Mx
1 pffiffiffiffiffiffiffiffiffi M12
(12a)
# (12b)
Subtracting Eq. (12b) from Eq. (12a) yields the following expression for c12 ðxÞ:
pffiffiffi 3 3Z pffiffiffi c12 ðxÞ ¼ c12 ð0Þ þ 256 2do
gI2 3=2
kI
2 gIII 3=2
kIII
!
sffiffiffiffiffiffiffiffiffi # " 1 M12 pffiffiffiffiffiffiffiffiffi 1 Mx M12 (12c)
(9)
Eq. (9) reduces to KI ¼ kI in the absence of third-order effects ðgI ¼ 0Þ. With third-order contributions, M1=2 dependence follows. Expressions that connect the elastic coefficients in Eqs. (5a) and (5b) to isotopic mass may be derived, as follows. Given that Eqs. (5a) and (8) each pertain to distortion I that preserves the orthogonality of the unit cell, the second-order coefficient 8d poffiffiffi ð3c11 þ 6c12 Þ in Eq. (5a) is equivalent to the second-order co-
!
An expression for the x-dependence of the bulk modulus follows from B ¼ ðc11 þ 2c12 Þ=3:
# pffiffiffi " ZgI2 1 3 pffiffiffiffiffiffiffi ; BðxÞ ¼ k pffiffiffi 12do I 8 2 k3=2 Mx I
(13a)
# pffiffiffi " ZgI2 1 3 pffiffiffiffiffiffiffiffiffi ; k pffiffiffi Bð0Þ ¼ 12do I 8 2 k3=2 M12 I
(13b)
efficient 2KI in Eq. (8). Similarly, a force constant for distortion III Please cite this article in press as: D.G. Clerc, H. Ledbetter, Isotope effect on diamond's elastic-stiffness coefficients: An ab initio study, Computational Condensed Matter (2017), http://dx.doi.org/10.1016/j.cocom.2017.03.003
4
D.G. Clerc, H. Ledbetter / Computational Condensed Matter xxx (2017) 1e9
Fig. 1. Lattice parameters of diamond's rhombohedral unit cell are depicted within the outlined f.c.c. unit cell. The asymmetric unit consists of atoms at the origin and the central (shaded) atom at (1/4, 1/4, 1/4). For simplicity, only the central atom and its four nearest neighbors are shown. The equilibrium geometry is a ¼ b ¼ c≡ao and a ¼ b ¼ g ¼ 60o . The Eulerian coordinate system is also indicated (axes 1, 2, and 3). The two smaller diagrams on the right-hand side depict the strains associated with distortions I and III.
pffiffiffi 3Z BðxÞ ¼ Bð0Þ þ pffiffiffi 96 2do
gI2 3=2
kI
!
sffiffiffiffiffiffiffiffiffi # " 1 M12 pffiffiffiffiffiffiffiffiffi 1 Mx M12
(13c)
The x-dependence of the lattice parameter from the Vogelgesang et al. study [3] is included here as an ancillary note:
ZgI ; aðxÞ ¼ a∞ qffiffiffiffiffiffiffiffiffiffiffiffiffiffi 6k3I Mx
(14a)
ZgI ffi; að0Þ ¼ a∞ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 6k3I M12
(14b)
sffiffiffiffiffiffiffiffiffi # " ZgI 1 M12 aðxÞ ¼ að0Þ þ pffiffiffi 3=2 pffiffiffiffiffiffiffiffiffi 1 Mx M12 6kI
(14c)
Here a∞ denotes the lattice parameter corresponding to infinite atomic mass (negligible zero-point motion). 4. Results and discussion Table 1 shows our principal results: lattice parameter, bulk modulus, second-order force constants kI and kIII , and third-order force constants gI and gIII . Computed values of the pair kI and gI , and the pair kIII and gIII , were obtained by fitting the computed data points to the following equations:
DU
ðI Þ
DU
ðIIIÞ
¼
¼
1 2 1 3 kε þ gε ; 2 I 1 6 I 1 1 1 k ε2 þ g ε3 ; 2 III 1 6 III 1
(15a)
(15b)
obtained respectively by substituting Eqs. (3b) and (3c) into Eq. (3a), and by substituting Eqs. (4b) and (4c) into Eq. (4a). Total energy curves are shown as functions of strain for distortions I and III in Figs. 2(A) and 3(A), respectively. Linear fits of total energy to the cube of strain are shown in Figs. 2(B) and 3(B). Our principal results for lattice parameter and second-order properties are consistent with measured values. The calculated lattice parameter, 0.35545 nm, is quite accurate, only 0.3% below the measured value of 0.35665 nm [8] and only 0.07% different
from our previously calculated value of 0.35569 nm [2]. The bulk modulus, B ¼ kI =9 ¼ 447 GPa differs from the measured value 443 GPa [9] by only 0.9%. The computed force constants kI ¼ 4026 GPa and kIII ¼ 2783 GPa are consistent with values obtained by substituting measured values of c11 and c12 [9] into Eqs. (3b) and (4b) yielding 3981 GPa and 2989 GPa. Comparison shows only a small difference for kI (1.1%), but more substantial difference for kIII (-6.9%). Our previously reported value of gI , 55,020 (3520) GPa [2], differs markedly from our current value, 36,195 (1692) GPa, because of the extent of anharmonic strain used in the computations. The previous value was computed using only eight data points whose maximum strain is 2.9%. In contrast, the current value was computed using 32 data points across which the maximum strain doubles to 5.8%. Fitting our current data to Eq. (3a) using the previous narrower range yielded larger values of gI , similar to the previous result. However, when additional data points farther into the anharmonic region are included, the magnitude of gI decreases as does the statistical-fitting error. The high quality of the fit across the anharmonic region, shown in Fig. 2(A) and 2(B), indicates that our current result for gI is more accurate. Table 1 also compares gI and gIII to measured values. These measured values are calculated from third-order elastic coefficients (TOECs) reported in studies based on Raman spectroscopy [10,11] and shock wave compression [12]. Compared to our result, values of gI from Raman spectroscopy are far larger in magnitude, by ~60%. In contrast, the value of gI from shock wave compression is ~33% larger. With respect to gIII , the Raman values of 20,492 GPa and 26,423 GPa are 47% to 90% larger than our result, 13,925 GPa. However, the shock compression value of 14,170 GPa is very similar to our gIII , differing by only 1.8%. The above analysis also yielded two interesting observations. First, fitting the total energy to Lagrangian strain as defined in Eq. (A.5a), instead of Eulerian strain, yields gI ¼ 48,087 (1884) GPa, a result quite consistent with shock compression: gI ¼ 47,910 GPa, differing by only 0.4%. Second, using the narrow 2.9% range of strain produces larger magnitude values of gI (~55,000 GPa) consistent with the Raman results, whereas using the 5.8% range yields values more consistent with shock compression. The effect of isotopic substitution may be inferred qualitatively from Eqs. (11c), (12c), (13c), and (14c). The value of c11 predicted by Eq. (11c) depends on the force constants through a positive quan3=2
tity, gI2 =kI
3=2
2 =k þ3gIII III , thereby indicating that c11 increases with
Please cite this article in press as: D.G. Clerc, H. Ledbetter, Isotope effect on diamond's elastic-stiffness coefficients: An ab initio study, Computational Condensed Matter (2017), http://dx.doi.org/10.1016/j.cocom.2017.03.003
D.G. Clerc, H. Ledbetter / Computational Condensed Matter xxx (2017) 1e9
5
Fig. 2. The computed total energy versus strain is shown for distortion type I. (A) The solid line fitting computed data points (diamonds) was generated using Eq. (3a) with computed values of the lattice parameter kI , and gI reported in Table 1. The dashed line is the Hooke-law result generated by fitting only the second-order term. Pairs of dashed vertical lines indicate the computed harmonic limit. (B) Third-order contributions to the total energy in Eq. (3a) are plotted versus the cube of strain. Open diamonds represent expansion (positive strain) and filled diamonds represent contraction (negative strain). The line represents the best fit to the data, using the value of gI reported in Table 1.
increasing x. In contrast, the value of c12 from Eq. (13c) depends 3=2
3=2
2 =k upon gI2 =kI gIII III , whose sign depends upon the relative values of the force constants. For diamond, the sign is positive, and thus c12 increases with increasing x. The value of B and a given by 3=2
3=2
Eqs. (13c) and (14c) depend respectively upon gI2 =kI and gI =kI , which are positive and negative quantities, respectively, thereby indicating that B increases with increasing x and a decreases with increasing x. Concerning the third independent elastic-stiffness coefficient c44, one can deduce its value by assuming that the Zener elastic anisotropy, 2c44 =ðc11 c12 Þ; is invariant for isotopic substitution. Using Eqs. (11c) and (12c) we obtain the following expression for the cubic-face-plane shear modulus:
pffiffiffi c44 ð0Þ 3Z pffiffiffi 256 2do c11 ð0Þ c12 ð0Þ sffiffiffiffiffiffiffiffiffi # ! " 2 9gIII gI2 1 M12 3=2 pffiffiffiffiffiffiffiffiffi 1 3=2 Mx M kIII kI 12
c44 ðxÞ ¼ c44 ð0Þ þ
(16)
Isotopic substitution was found to induce very small changes in the elastic coefficients, the bulk modulus and the lattice parameter. Using Eqs. (11c), (12c), (13c), (14c), and (16), with computed values of c11 ð0Þ; c12 ð0Þ; Bð0Þ; að0Þ; and c44 ð0Þ; from Ref. [2], which are 1043 GPa, 128 GPa, 433 GPa, 0.35569 nm, and 534 GPa, respectively, fractional changes associated with full substitution 13C relative to 12 C have the following values:
½c11 ð1Þ c11 ð0Þ =c11 ð0Þ ¼ 4:9 105 ;
(17a)
½c12 ð1Þ c12 ð0Þ =c12 ð0Þ ¼ 2:5 104 ;
(17b)
½c44 ð1Þ c44 ð0Þ =c44 ð0Þ ¼ 2:1 105 ;
(17c)
½Bð1Þ Bð0Þ =Bð0Þ ¼ 8:8 105 ;
(17d)
½að1Þ að0Þ =að0Þ ¼ 1:5 105 :
(17e)
The magnitudes of the above changes are very small: varying from only 0.0015% to 0.025%. Essentially the same result occurs
Please cite this article in press as: D.G. Clerc, H. Ledbetter, Isotope effect on diamond's elastic-stiffness coefficients: An ab initio study, Computational Condensed Matter (2017), http://dx.doi.org/10.1016/j.cocom.2017.03.003
6
D.G. Clerc, H. Ledbetter / Computational Condensed Matter xxx (2017) 1e9
Fig. 3. The computed total energy versus strain is shown for distortion type III. (A) The solid line fitting computed data points (diamonds) was generated using Eq. (4a) with computed values of the lattice parameter kIII , and gIII reported in Table 1. The dashed line is the Hooke-law result generated by fitting only the second-order term. Pairs of dashed vertical lines indicate the computed harmonic limit. (B) Third-order contributions to the total energy in Eq. (4a) are plotted versus the cube of strain. Open diamonds represent expansion (positive strain) and filled diamonds represent contraction (negative strain). The line represents the best fit to the data, using the value of gIII reported in Table 1.
1=2
after accounting for the small nonlinearity (Mx dependence). For example, values of ½c11 ðxÞ c11 ð0Þ=c11 ð0Þ at x ¼ 0, x ¼ 1 2, and x ¼ 1 are respectively 0%, 0.0025% and 0.0049%. These three values average 0.0025%, which is comparable to the 0.0049% reported for Eq. (17a). We find that the Cauchy discrepancy c12/c44 stays the same within calculation error. Often taken as a clue to interatomicbonding nature, departure from central forces, diamond’s negative departure from unity is the largest for any cubic element and often called a “negative Cauchy pressure” [13]. Our predicted 0.0088% increase in B is comparable to the 0.18% increase reported by Vogelgesang et al. [3], but it is inconsistent with a measured increase of 17% reported by Hurley et al. [4]. The difference originates in the value of c12 that was used to calculate B. Specifically, the measured 87% increase in c12 at x ¼ 1 is enormously larger than our predicted 0.025% increase. Our result for c12, a very small change at x ¼ 1 of 0.025%, is consistent with a simple theory based on Heisenberg's theorem. Gilman, in describing diamond's extreme hardness, derived a basic bulk-modulus relationship B ¼ 3e2 =40pro4 , where e denotes the valence charge of the atom and ro the equilibrium radius of the valence electron cloud around a single atom [14]. Thus, the physical
property that governs the bulk modulus is atomic volume, not atomic mass, and efforts to increase the bulk modulus, which reflects mechanical hardness, may focus on reducing ro . The additional neutron that exists in 13C, being uncharged, does not significantly affect the electronic charge distribution and hence does not significantly alter ro. As a result, we expect the corresponding elastic constants and bulk modulus to be affected little by 13 C12C substitution. 5. Conclusions The effects of isotopic substitution of 13C on the elastic constants of diamond were studied theoretically using predicted properties taken from our ab initio studies, and using an analysis of zero-point motion and the associated anharmonicity of lattice vibrations. We predict that changes in c11 , c12 , c44 , and B with respect to x are essentially linear: c11 ðxÞ ¼ c11 ð0Þð1 þ 0:000049xÞ; c12 ðxÞ ¼ c12 ð0Þð1 þ 0:00025xÞ; c44 ðxÞ ¼ c44 ð0Þð1 þ 0:000021xÞ; and BðxÞ ¼ Bð0Þð1 þ 0:000088xÞ Compared with values at x ¼ 0, c11 ; c12 ; c44 ; and B are predicted to change by only 0.0049%, 0.025%, 0.0021%, and 0.0088%, respectively, at x ¼ 1. Our predicted 0.0088% increase in B is comparable to the 0.18% increase reported by
=
Please cite this article in press as: D.G. Clerc, H. Ledbetter, Isotope effect on diamond's elastic-stiffness coefficients: An ab initio study, Computational Condensed Matter (2017), http://dx.doi.org/10.1016/j.cocom.2017.03.003
D.G. Clerc, H. Ledbetter / Computational Condensed Matter xxx (2017) 1e9
7
Table 1 Previously reported and present results. Quantity Lattice parameter Bulk modulus Second-order force constants
a (nm) B (GPa) kI (GPa) kIII (GPa) gI (GPa) gIII (GPa)
Third-order force constants a b c d e f g h i
Prior studies
This study
0.35665a, 0.35569b 443d 3981e 2989e 58,788g, 54,309h, 47,910i 20,492g, 26,423h, 14,170i
0.35545c 447 (3) 4026 (25)f 2783 (11)f 36,195 (1692) 13,925 (556)
Reference [8]. Reference [2]. The computed f.c.c. lattice parameter is the average for distortions I (0.35545 nm) and III (0.35563 nm). Reference [9]. Calculated using Eqs. (3b) and (4b), and the measured c11 and c12 from Ref. [9]. pffiffiffi Terms of Eq. (6), equivalent values of kI and kIII , are 1844 J/m2 and 1365 J/m2, respectively - obtained by multiplying the GPa values by 16do =3 3. Calculated using Eqs. (3c) and (4c), and three Raman-measured TOECs in Ref. [10]. Calculated using Eqs. (3c) and (4c), and three Raman-measured TOECs in Ref. [11]. Calculated using Eqs. (3c) and (4c), and the average of two sets of three shock-wave-measured TOECs in Ref. [12].
Vogelgesang et al. [3]. The difference between our predicted increase in B and the measured 17% increase in B reported by Hurley et al. [4] originates in the measured 87% increase in c12 at x ¼ 1, which is far larger than our predicted 0.025% increase. The method used in this study may also be used for other crystals having diamond's space group, and as such could lead to further studies within a similar theoretical framework. Acknowledgments This research received partial support from Los Alamos National Laboratory (Contract No. 96628-001-04 30). The authors appreciate the encouragement of this study by Dr. Albert Migliori (Los Alamos National Laboratory).
1 2
h1 ¼ h2 ¼ ε1 þ ε21 ;
III :
1 2
h3 ¼ ε1 þ ε21 ¼ h1 þ ε21 z h1 ;
orthogonality of the solid, the derivatives in Eq. (A.3) that possess cross terms vanish, yielding the following relationship:
hpp
1 ¼ 2
vup vxp
2
vup þ2 vxp
! (A.4)
The remaining terms can be related to Eulerian strain by using Eq. (A.1) with i ¼ j to yield εii ¼ vui =vxi and hence vup =vxp ¼ εpp : That result can be substituted into Eq. (A.4) to yield the connection between Lagrangian and Eulerian strain for distortions I and III, as follows:
I:
1 2
h1 ¼ h2 ¼ h3 ¼ ε1 þ ε21 ;
h4 ¼ h5 ¼ h6 ¼ 0;
h4 ¼ h5 ¼ h6 ¼ 0:
(A.5a)
(A.5b)
Appendix A. Eulerian and Lagrangian frameworks Strain εi relates directly to the strain denoted by εij using the familiar Voigt notation. In turn, εij is connected to the Eulerian or classical definition of strain by the following relationship:
εij ¼
1 vui vuj þ 2 vxj vxi
! (A.1)
Here u denotes infinitesimal displacement in the strained solid. Lagrangian strain hpq accounts for finite deformation by including quadratic strain effects:
hpq ¼
1 vXi vXi dpq 2 vxp vxq
(A.2)
Here, Xi ði ¼ 1; 2; 3Þ denotes position in the strained solid with implicit summation over index i, and xp and xq denote position in the unstrained solid, where p; q2f1; 2; 3g. By using displacement coordinates, ui ¼ Xi xi , Eq. (A.2) can be written as:
hpq ¼
1 vui vui vui vu þ diq þ i dip þ dip dip dpq 2 vxp vxq vxp vxq
(A.3)
Appendix B. Conventions used to define strain The smallest periodic unit in diamond is the rhombohedron, depicted in Fig. 1. The vertices of the rhombohedral cell are defined by four atoms: one at the origin plus three at each face center of the standard FCC supercell. The interior of the rhombohedral cell contains one additional atom at (1/4, 1/4, 1/4). The lattice vectors of the rhombohedral cell are defined as follows in a Cartesian basis:
! a ¼ 〈a cosðp=4Þ; 0; a sinðp=4Þ〉;
(B.1a)
! b ¼ 〈b cosðp=4Þ; b sinðp=4Þ; 0〉;
(B.1b)
! c ¼ 〈0; c sinðp=4Þ; c sinðp=4Þ〉:
(B.1c)
If we define b ¼ ao ð1 þ ε1 Þ and a ¼ ao ð1 þ 4ε1 Þ, where 4 is a constant, then Eqs. (B.1a)e(B.1c) simplify as follows:
ao ! a ¼ pffiffiffi h1 þ ε1 ; 0; 1 þ 4ε1 i; 2
(B.2a)
Because the two types of strain used in this study preserve the Please cite this article in press as: D.G. Clerc, H. Ledbetter, Isotope effect on diamond's elastic-stiffness coefficients: An ab initio study, Computational Condensed Matter (2017), http://dx.doi.org/10.1016/j.cocom.2017.03.003
8
D.G. Clerc, H. Ledbetter / Computational Condensed Matter xxx (2017) 1e9
! ao b ¼ pffiffiffi h1 þ ε1 ; 1 þ ε1 ; 0i; 2
(B.2b)
ao ! p c ¼ ffiffiffi h0; 1 þ ε1 ; 1 þ 4ε1 i 2
(B.2c)
Euclidian inner products and vector magnitudes have the following forms:
! ! ! ! ao a $ b ¼ b $ c ¼ ð1 þ ε1 Þ2 ; 2 2
(B.3a)
! ! ! d a¼ C 2 a ao ¼ pffiffiffi hð1þε1 Þðf1 þf2 1Þ;ð1þε1 Þðf2 þf3 Þ;ð1þ4ε1 Þðf1 þf3 1Þi; 2 (C.1b)
! ! ! d b¼ C 2 b ao ¼ pffiffiffi hð1þε1 Þðf1 þf2 1Þ;ð1þε1 Þðf2 þf3 1Þ;ð1þ4ε1 Þðf1 þf3 Þi; 2 (C.1c)
2
! ! ao a $ c ¼ ð1 þ 4ε1 Þ2 ; 2 ao ! ! j a j ¼ j c j ¼ pffiffiffi 2
(B.3b)
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð1 þ ε1 Þ2 þ ð1 þ 4ε1 Þ2 ;
! b ¼ ao ð1 þ ε1 Þ
(B.3c)
(B.3d)
the relations cosðaÞ ¼ b b$b c , cosðbÞ ¼ b a $b c , and b where the carat denotes unit magnitude, we arrive at cosðgÞ ¼ b a $ b, the following equations: Using
ð1 þ ε1 Þ ; cosðaÞ ¼ cosðgÞ ¼ pffiffiffiqffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ð1 þ ε1 Þ2 þ ð1 þ 4ε1 Þ2
cosðbÞ ¼
ð1 þ 4ε1 Þ
(B.4a)
2
ð1 þ ε1 Þ2 þ ð1 þ 4ε1 Þ2
:
! ! ! d c¼ C 2 c ao ¼ pffiffiffi hð1þε1 Þðf1 þf2 Þ;ð1þε1 Þðf2 þf3 1Þ;ð1þ4ε1 Þðf1 þf3 1Þi 2 (C.1d) Using f1 ¼ f2 ¼ f3 ¼ 1/4, Eqs. (C.1aed) simplify as follows:
! ao d O ¼ pffiffiffi h1 þ ε1 ; 1 þ ε1 ; 1 þ 4ε1 i; 2 2
(C.2a)
! ao d a ¼ pffiffiffi h 1 ε1 ; 1 þ ε1 ; 1 4ε1 i; 2 2
(C.2b)
! ao d b ¼ pffiffiffi h 1 ε1 ; 1 ε1 ; 1 þ 4ε1 i; 2 2
(C.2c)
! ao d c ¼ pffiffiffi h1 þ ε1 ; 1 ε1 ; 1 4ε1 i: 2 2
(C.2d)
(B.4b)
The three sets of preceding equations show that the geometry of diamond's rhombohedral unit cell, for a constrained distortion such as I or III, is defined completely by the single Eulerian strain coordinate ε1 . For distortion I ð4 ¼ 1Þ, Eqs. (B.3a) and (B.3b) reduce to ! ! ! j a j ¼ b ¼ j c j and Eqs. (B.4a) and (B.4b) yield a ¼ b ¼ g ¼ 60o . Relations for distortion III may be obtained through the substitution 4 ¼ 1.
Appendix C. Effects of internal relaxation
These four bonds have the same length, ! ! ! ! j d O j ¼ j d a j ¼ j d b j ¼ j d c j in distortion I ð4 ¼ 1Þ and in distortion III ð4 ¼ 1Þ. For either type of distortion, bond-stretching contributions to the potential energy can be written in terms of one variable ε1 , analogous to the use of the single variable in the potential energy shown in Eq. (3a). The angles between the central atom and the four nearest neighbors can be found analogous to qOa ¼ ; ! ! ! ! ! ! ð d O ; d a Þ ¼ arccosð d O $ d a = d O d a Þ; the following equations emerge:
This appendix shows that distortions I and III lead to equations for the internal coordinate which are consistent with the form of the potential energy used in Eq. (3a). This internal coordinate is the position of the shaded central atom C2 in Fig. 1, defined in term of ! ! ! ! the three primitive lattice vectors as C 2 ¼ f1 a þ f2 b þ f3 c ; where ff1 ; f2 ; f3 g are constants whose equilibrium values are f1 ¼ f2 ¼ f3 ¼ 1=4. Combining this definition with the definitions of the lattice vectors in Eqs. (B.2a), (B.2b), and (B.2c) yields the following vectors for the bonds connecting C2 to its four nearestneighbors:
! ! dO ¼ C2 ao ¼ pffiffiffi hð1 þ ε1 Þðf1 þ f2 Þ; ð1 þ ε1 Þðf2 þ f3 Þ; ð1 þ 4ε1 Þðf1 þ f3 Þi; 2 (C.1a)
cosðqOa Þ ¼ cosðqOc Þ ¼ cosðqab Þ ¼ cosðqbc Þ ð1 þ 4ε1 Þ2 ; ¼ ð1 þ 4ε1 Þ2 þ 2ð1 þ ε1 Þ2 cosðqOb Þ ¼ cosðqac Þ ¼
ð1 þ 4ε1 Þ2 2ð1 þ ε1 Þ2 ð1 þ 4ε1 Þ2 þ 2ð1 þ ε1 Þ2
(C.3a)
:
(C.3b)
For distortion I ð4 ¼ 1Þ, the six cosines are equal to 1/3 irrespective of the level of strain, corresponding to the tetrahedral angle of 109.47. For distortion III ð4 ¼ 1Þ, the six bond angles split into two sets given by Eqs. (C.3a) and (C.3b), with the angles being equal within each set, the particular value depending upon the level of strain. Thus, the change in potential energy associated with distortion III is a function of bond angle as well as bond length. Angular contributions to the potential energy were found to be negligible, as reported in the Computational Method section.
Please cite this article in press as: D.G. Clerc, H. Ledbetter, Isotope effect on diamond's elastic-stiffness coefficients: An ab initio study, Computational Condensed Matter (2017), http://dx.doi.org/10.1016/j.cocom.2017.03.003
D.G. Clerc, H. Ledbetter / Computational Condensed Matter xxx (2017) 1e9
References [1] E. Ekimov, V. Sidorov, E. Bauer, et al., Nature 428 (2004) 542e545. [2] D.G. Clerc, H. Ledbetter, J. Phys. Chem. Sol. 66 (2005) 1589e1597. [3] R. Vogelgesang, A.K. Ramdas, S. Rodriquez, M. Grimsditch, T.R. Anthony, Phys. Rev. B 54 (6) (1996) 3989e3999. Eq. (6) of this study, with opposite sign in the third-order coefficient, corresponds to Eq. (17). Similarly, Eqs. (14a) and (14b) of this study both correspond to Eq. (22), and Eq. (14c) corresponds to Eq. (29). The isotope effect on the bulk modulus (0.18%) is reported on page 3997 in Eq. (38). [4] D.C. Hurley, R.S. Gilmore, W.F. Banholzer, J. Appl. Phys. 76 (12) (1994) 7726e7730. See page 7730, paragraph 3. [5] V. Plekhanov, Mater. Sci. Eng. R. 35 (2001) 139e237. See page 153, paragraph 1. [6] W. Weber, Phys. Rev. B 15 (10) (1977) 4789e4803. See Table II on page 4793. [7] G. Burns, A. Glazer, Space Groups for Solid State Scientists, Academic Press,
[8] [9] [10] [11] [12] [13] [14]
9
Inc., San Diego, second ed., 1990. See page 152, Fig. 6e5(g) and page 157, Fig. 6e8. T. Sato, K. Ohashi, T. Sudoh, et al., Phys. Rev. B 65 (2002) 092102. See page 3, Fig. 2 (0 K values). H. McSkimin, P. Andreatch, P. Glynn, J. Appl. Phys. 43 (3) (1972) 985e987. See page 987, Table II. M.H. Grimsditch, E. Anastassakis, M. Cardona, Phys. Rev. B 18 (2) (1978) 901e904. See page 904, Table II. E. Anastassakis, A. Cantarero, M. Cardona, Phys. Rev. B 41 (11) (1990) 7529e7535. See page 7533, Table IV. J.M. Winey, A. Hmiel, Y.M. Gupta, J. Phys. Chem. Sol. 93 (2016) 118e120. See page 119, Table 2: sets A and B. H. Ledbetter, S. Kim (Eds.), Handbook of Elastic Properties of Solids, Liquids, and Gases, Academic, San Diego, 2001. J.J. Gilman, Mater. Sci. Eng 7 (1971) 357e361. See page 359, Eq. (20).
Please cite this article in press as: D.G. Clerc, H. Ledbetter, Isotope effect on diamond's elastic-stiffness coefficients: An ab initio study, Computational Condensed Matter (2017), http://dx.doi.org/10.1016/j.cocom.2017.03.003