Thermal equation of state of silicon carbide

The significance of SiC is manifested through its broad applications at extreme conditions (high pressure, temperature, radiation, voltage, etc.). For instance, in addition to being widely used for polishing, drilling, cutting, and grinding, SiC is a promising candidate material to encapsulate the fissile fuel in nuclear reactors because of its high melting point, high thermal conductivity, and small neutron capture cross section.¹ Furthermore, SiC is a robust material for fabricating composite materials. A strong composite consisting of aluminum oxide and SiC whiskers was first developed more than three decades ago.² Nowadays, more effort has been put into discovering its advantages in aviation and aerospace industries,³ because its unique properties, such as high strength, excellent thermal shock resistance, and low density comparable to aluminum, meet numerous special technical requirements of aerospace industry. SiC is also an energy related material. As a reinforcing phase, it has been a long-term focus of investigation for developing novel superhard composite materials that are suitable to drill wells for the exploration and extraction of natural gas and crude oil.^4–10 

Thermoelasticity is a fundamental property of condensed matter. However, except a few theoretical computations,^11,12 the experimental investigation of thermal elastic properties is scarce for SiC. In the present work, we study SiC at simultaneously high P-T conditions by using a large volume press integrated with in-situ synchrotron X-ray radiation. By fitting the P-V-T data to a modified high-T Birch-Murnaghan equation of state (EOS), we derived a set of thermal EOS parameters for SiC, such as temperature derivative of the bulk modulus and pressure derivative of the thermal expansion. Furthermore, a thermal pressure approach was used to calculate the temperature derivative of the bulk modulus at a constant volume, (∂K_T/∂T)_V, which is experimentally difficult to be measured.

Additionally, computational simulations were conducted with density functional theory (DFT). The effect of temperature was included through the harmonic lattice vibration contribution on the Helmholtz free energy. The experimentally obtained values were compared with those determined from the simulations, and a strong validation has been established between those two.

At ambient conditions, the crystal structure of SiC is FCC (face-centered cubic) with a space group of Fm3m (225). The starting materials of SiC powders with gain sizes of 10 μm or smaller were purchased from Sigma-Aldrich. The in-situ high P-T synchrotron XRD experiment was carried out via using a DIA-type cubic anvil apparatus at the beam line X17B2 of National Synchrotron Light Source, Brookhaven National Laboratory. The sample under varied P-T conditions was tested by a white radiation, and the diffracted x-rays were recorded by a 13-element detector at a fixed Bragg angle of 2θ = 6.4900°. Details of the experimental setup as well as the preparation of the cell assembly have been described elsewhere.¹³ Here, only a brief description is given. The pre-compressed boron epoxy (4:1 in weight ratio) with a cubic shape was employed as a pressure-transmitting medium to maintain a quasi-hydrostatic environment surrounding the sample. The SiC sample sandwiched by salt (NaCl) powders was compacted into a sleeve made from boron nitride (BN). Then, this BN sleeve was placed into an amorphous carbon tube that was used as a furnace to heat the sample. Finally, this sample assembly was inserted into the central hole in the boron epoxy cube.

In the present study, we used salt as an internal pressure marker to determine the sample pressure. At each P-T condition, XRD pattern and Decker's EOS of salt were used to calculate the sample pressure.¹⁴ The uncertainty of the sample pressure induced from the statistical variation of the X-ray peak positions of salt is typically smaller than 0.2 GPa in the present P-T range. A W/Re25%-W/Re3% thermocouple in direct contact with the sample was used to measure the sample temperature. Uncertainty of the sample temperature was approximated to be ∼10 K. The temperature gradients over the whole sample length and along the sample radial direction are of the order of 20 K, and less than 5 K, respectively, at the highest reached temperature 1100 K.¹³ Furthermore, we chose the positions near the thermocouple junction to collect the X-ray patterns of the sample and salt in order to minimize the uncertainty of the measurements. All of the data were collected from a single high P-T experimental run, which can avoid the systematic errors that are typically present if the data are obtained from several separate experimental runs.

DFT simulations were conducted with the Vienna Ab initio Simulation Package (VASP).^15–18 Potentials of Si and C were selected via using projector-augmented wave (PAW) method^19,20 with both local density approximation (LDA)²¹ and PW91 generalized gradient approximation (GGA).^22,23 The plane wave cutoff energy was chosen to be 520 eV. A 2 × 2 × 2 enlarged supercell of SiC conventional unit cell was constructed with a 4 × 4 × 4 Monkhorst k-point mesh^24,25 to ensure convergence. During the iterative minimization of the total energy, Gaussian smearing was used with a sigma value as small as 0.1 until the convergence criterion of 10⁻⁵eV was met. To include the thermal contribution at ambient pressure, we used the quasi-harmonic approximation functionality of the code package Phonopy.^26–28 At each volume around the energy minimum, force sets were calculated with finite displacement method, and then were used to obtain force constants, dynamical matrix, and phonon frequencies and eigenvectors on a 41 × 41 × 41 q-point mesh. They were further transformed into the temperature dependent vibrational Helmholtz free energies, which were added to the total energy obtained through a VASP static self-consistent calculation.^27,28 At each temperature from 0 K to 1100 K with a step of 5 K, Helmholtz free energies of SiC structures at 7 different volumes around the energy minimum point were fitted to the Birch-Murnaghan EOS, treating K′ as a variable. Temperature-resolved parameters, K and α, were derived from the fits, and because of the small temperature step, those two parameters versus temperature were depicted with smooth curves, as shown in Figs. 4(b) and 4(c).

The measurements were conducted by following the P-T route shown in Fig. 1. The XRD pattern collected at each P-T condition was fitted by Gaussian function to derive the lattice parameters and then to calculate the unit cell volume V of SiC at that P-T point. The data sets of P-V-T are demonstrated in Fig. 2. In the present study, a modified high-T Birch-Murnaghan EOS (truncated to third order)^29–33 was employed to fit the P-V-T data set and derive the thermoelastic parameters, as illustrated in Fig. 2. This modified EOS could be formulated with the following equation:^29–31 

where

and

The physical meanings of the variables in the above equations are given as follows:

• K_To: isothermal bulk modulus at ambient temperature

• K_T: isothermal bulk modulus at a nonambient temperature T

• (∂K/∂T)_P: temperature derivative of the bulk modulus

• (∂K/∂P)_T: pressure derivative of the bulk modulus.

• V_o: unit-cell volume at ambient conditions

• V_T: unit-cell volume at atmospheric pressure and temperature T

• V_PT: unit-cell volume at a nonambient P-T condition

• α(0, T): volumetric thermal expansion at atmospheric pressure, commonly formulated by α(0, T) = a + bT − c/T² (T in Kelvin).³⁴ 

In Equation (1), the thermal effect on material's response to pressure is considered by replacing K_o and V_o/V_P, typical parameters in isothermal EOS, with K_T and V_T/V_PT, respectively. During the data processing, K' is set to be 4. Meanwhile, in order to simplify the data fitting, the higher order terms, such as c/T² in α(0, T), ∂²K/∂T², and ∂²K/∂P∂T, were not taken into account. A least-squares fit to the P-V-T data with Equation (1) yields K_To = 237(2) GPa, (∂K/∂T)_P = −0.037(4) GPa K⁻¹, and α(0, T) = a + bT with a = 5.77(1) × 10⁻⁶ K⁻¹ and b = 1.36(2) × 10⁻⁸ K⁻². Errors, enclosed in the parentheses, come from the least-squares fitting. The pressure derivative of the volumetric thermal expansion, (∂α/∂P)_T, was calculated from the thermodynamic identity

and was found to be −6.53 ± 0.64 × 10⁻⁷ K⁻¹GPa⁻¹. We approximated the uncertainty of (∂α/∂P)_T from the error propagation of K_To and (∂K/∂T)_P.

Another thermoelastic parameter (∂K_T/∂T)_V, the temperature derivative of the bulk modulus at a constant volume, is important and conveys many thermodynamic insights, but is difficult to be determined directly from experimental measurement. Here, we used a thermal pressure approach^31,35–39 to derive it. In this approach, thermal pressure, P_th, is defined as the difference between two pressures with the same volume, one measured at a given temperature and the other calculated from Equation (1) at room temperature. The thermal pressure against temperature is plotted in Fig. 3. Meanwhile, there is an alternative method,^35,37 which can be used to calculate the thermal pressure with a given volume at any temperature above 300 K. It is represented by the following:

Then, fitting the plot of thermal pressure versus temperature using Equation (3) as shown in Fig. 3 yields the values of αK_T(V₃₀₀, T) and (∂K_T/∂T)_V: 0.0032(6) and −0.028(4)GPa K⁻¹, respectively. Furthermore, from the thermodynamic identity

(∂K_T/∂T)_P was calculated to be −0.041(6), which is consistent with the value derived from Equation (1). We recalculated (∂α/∂P)_T by substituting (∂K_T/∂T)_P into Equation (2), resulting in a value of −7.23 ± 0.95 × 10⁻⁷ K⁻¹GPa⁻¹. It is in agreement with the one derived from the Equation (1). All of the results regarding the thermoelasticity of SiC, obtained from current study and previous investigations using different methods, are listed in Table I. The bulk modulus K_o of SiC determined in the present study is slightly smaller than the value obtained using a diamond anvil cell (DAC).^40,41 This difference may be partially due to the inherently different instrumental responses of two experimental setups (DAC and large volume press) to non-hydrostatic stress distribution within samples.

Fig. 4(a) shows the LDA generated curve set of Helmholz free energy against volume at temperatures from 0 K to 1100 K. Table I also contains the values obtained from the simulations. We see that the bulk modulus (K₀) values computed by LDA and PW91 are smaller than the experimental value by 6.7% and 13.9%, respectively. For the absolute values of (∂K/∂T)_P, the two computed results are also smaller than the experimental counterparts. The difference between theoretical and experimental values from Fig. 4 is small and comparable to those in the literature where theoretical values for K_To vary from 200 to 250 GPa.⁴² The variety of functionals and implementation of the first-principles methods lead to this variation. Apart from these variations in computational techniques, the single crystalline setup in the simulation is different from the actual experimental high P-T conditions where microstructures and complexity of stress distribution across the sample exist. Modeling these higher order effects is outside the scope of current state of the art first-principles computations.

As described earlier, we assumed that K′ = 4 in the procedure of experimental data fitting. In Table I, the values of computed K′ are 3.9 (LDA) and 4.0 (GGA). Throughout the temperature range of the simulations, K′ remains between 3.9 and 4.0 for both approximations. Therefore, K′ = 4 was a very reasonable assumption. We also assumed that both K and α linearly depend on T, characterized by (∂K/∂T)_P and a + bT, respectively. Since the simulations give us more sampling points of T, we analyzed the linearity of K and α with the data output of the LDA simulation. In Fig. 4(b), from 300 K to 1100 K, the temperature range explored in the experiment, K does linearly depend on T. On the other hand, a curvature of α in Fig. 4(c) is slightly more noticeable in the same temperature range. This indeed adds imprecision to the estimation of α at a specific T. Table II summarizes this effect. The ambient condition (300 K) values of α determined from experiment or two simulations through a + bT are all greater than 9 × 10⁻⁶ K⁻¹. However, if we use the direct output from Phonopy calculations, described by the curve in Fig. 4(c), the values, 7.19 × 10⁻⁶ K⁻¹ (LDA) and 8.26 × 10⁻⁶ K⁻¹ (PW91), are much closer to the reference value 8.31 × 10⁻⁶ K⁻¹.

In summary, we have performed in situ XRD experiment for SiC under pressures up to 8.1 GPa and temperatures up to 1100 K. By fitting the measured P-V-T data to a modified high-T Birch-Murnaghan EOS, a complete set of thermoelastic parameters has been derived. An alternative method, namely, thermal pressure approach, has been used to calculate αK_T(V₃₀₀, T) and (∂K_T/∂T)_V, complementing the results obtained from the modified Birch-Murnaghan EOS. The bulk modulus K_o determined in the present study is slightly smaller than those measured previously by using a DAC. We also performed DFT calculations with two different approximations, LDA and GGA-PW91. The computed results generally agree with those determined by the experiment, and the small discrepancy of the volumetric thermal expansion is attributed to its deviation from the linear relationship with temperature.

The experimental work was carried out at the beam line X17B2 of National Synchrotron Light Source, Brookhaven National Laboratory, which is supported by the Consortium for Materials Properties Research in Earth Sciences (COMPRES) under NSF cooperative Agreement EAR 01-35554. Theoretical work was supported by NSF CMMI 1234777 and the Ohio Supercomputer Cluster.⁴³ Portion of this research was supported by the Michigan Space Grant Consortium, the Research Faculty Fellowship of Oakland University, and the National Nuclear Security Administration under the Stewardship Science Academic Alliances program through DOE Cooperative Agreement #DE-NA0001982.