# Carbon under extreme conditions: Phase boundaries and electronic properties from first-principles theory

See allHide authors and affiliations

Communicated by Russell J. Hemley, Carnegie Institution of Washington, Washington, DC, December 8, 2005 (received for review October 6, 2005)

## Abstract

At high pressure and temperature, the phase diagram of elemental carbon is poorly known. We present predictions of diamond and BC8 melting lines and their phase boundary in the solid phase, as obtained from first-principles calculations. Maxima are found in both melting lines, with a triple point located at ≈850 GPa and ≈7,400 K. Our results show that hot, compressed diamond is a semiconductor that undergoes metalization upon melting. In contrast, in the stability range of BC8, an insulator to metal transition is likely to occur in the solid phase. Close to the diamond/liquid and BC8/liquid boundaries, molten carbon is a low-coordinated metal retaining some covalent character in its bonding up to extreme pressures. Our results provide constraints on the carbon equation of state, which is of critical importance for devising models of Neptune, Uranus, and white dwarf stars, as well as of extrasolar carbon-rich planets.

Elemental carbon has been known since prehistory, and diamond is thought to have been first mined in India >2,000 years ago, although recent archaeological discoveries point at the possible existence of utensils made of diamond in China as early as 4,000 before Christ (1). Therefore, the properties of diamond and its practical and technological applications have been extensively investigated for many centuries. In the last few decades, after the seminal work of Bundy and coworkers (2) in the 1950s and `60s, widespread attention has been devoted to studying diamond under pressure (3). For example, the properties of diamond and, in general, of carbon under extreme pressure and temperature conditions are needed to devise models of outer planet interiors (e.g., Neptune and Uranus) (4–6), white dwarfs (7, 8) and extrasolar carbon planets (9).

Nevertheless, under extreme conditions the phase boundaries and melting properties of elemental carbon are poorly known, and its electronic properties are not well understood. Experimental data are scarce because of difficulties in reaching megabar (1 bar = 100 kPa) pressures and thousands of Kelvin regimes in the laboratory. Theoretically, sophisticated and accurate models of chemical bonding transformations under pressure are needed to describe phase boundaries. In most cases, such models cannot be simply derived from fits to existing experimental data, and one needs to resort to first-principles calculations, which may be very demanding from a computational standpoint.

It has long been known that diamond is the stable phase of carbon at pressures above several gigapascal (2). Total energy calculations (10, 11) based on Density Functional Theory (DFT) predict a transition to another fourfold coordinated phase with the BC8 symmetry¶ at ≈1,100 GPa and 0 K, followed by a transition to a simple cubic phase at pressures >3,000 GPa. These transitions have not yet been investigated experimentally, because the maximum pressure reached so far in diamond anvil cell experiments on carbon is 140 GPa (14). The only experimental information on diamond melting in the approximately 10^{2}- to 10^{3}-GPa range is from recent shock-wave experiments, where melting and a transition to a conducting fluid were observed (15). However, these measurements were insufficient to locate phase boundaries and inconclusive as to whether diamond undergoes an insulator-to-metal transition before or after melting.

A number of experimental studies (16, 17) and *ab initio* simulations (18–20) have focused on the phase boundaries close to the graphite–diamond–liquid triple point (≈4,000 K and 15 GPa). Unlike other solids with the diamond structure (e.g., Si and Ge), at these pressures the melting line of carbon has a positive slope (19). However, at much higher pressure (several hundred GPa), Martin and Grumbach (21) suggested that the slope becomes negative. This result was inferred from the comparison of the specific volumes of liquid and solid carbon computed from *ab initio* molecular dynamics. Empirical potentials (22–25) have also been used to investigate melting of carbon with different degrees of success. In general, empirical potentials lack predictive power over large density variations and, more specifically, important discrepancies between *ab initio* calculations and empirical potentials have been recently found, e.g., for the existence of a liquid/liquid phase transition at low pressure (20, 24, 26).

In this paper, we report on the properties of carbon in the pressure and temperature range of 15–2,000 GPa and 0–10,000 K, as obtained from first-principles calculations based on DFT. We have investigated solid/liquid phase boundaries and analyzed structural and electronic transformations as pressure and temperature are increased. In particular, we have determined diamond and BC8 melting lines and predicted the location of the diamond/BC8/liquid triple point, inferred from the crossing of computed phase boundaries. Structural and electronic properties of diamond and BC8, and the solid/solid transition in this range of temperatures have also been investigated by a combination of different techniques.

## Results

We first computed the equation of state of diamond at zero temperature to test the accuracy of our theoretical and numerical approach. The agreement with the available experimental data in the range 0–140 GPa (14) is excellent. For example, the computed equilibrium volume only differs by 1% from the experimental value, and optical phonon frequencies agree with experiment within 3% (14). This level of accuracy is already obtained with 64-atom supercells and a Γ-point sampling of the Brillouin zone, as evident from convergence tests performed with supercells containing up to 216 atoms and with unit cells with up to 12 × 12 × 12 *k*-point meshes. Consistent with previous studies, we find that at zero temperature the transition pressure from diamond to BC8 is 1,075 GPa. Other candidate structures previously proposed in the literature, such as R8 (27), ST12 (28), and hexagonal diamond (29), were also examined but were found to be unstable relative to BC8 in the pressure range studied here.∥

Calculations of diamond and BC8 melting lines were carried out by direct simulation of solid/liquid coexistence (i.e., by using a so-called two-phase simulation method). This approach has been recently shown to be applicable and efficient within the framework of *ab initio* molecular dynamics simulations of comparable scale (30–32). The two-phase method consists of choosing a given thermodynamic condition in (*P, T*) space, where liquid and solid systems are prepared and equilibrated. The phases are then put in contact and constant pressure–constant temperature simulations are carried out to determine which of the two phases is more stable at that given (*P, T*). As simulation progresses, the more stable phase eventually fills the entire volume available in the simulation. The melting line is located by repeating this process at different (*P*_{n}, *T*_{n}) points. Several complementing criteria, including diffusion coefficients, were used to determine whether melting or crystallization occurs. In addition, correlation functions, coordination numbers, and local order parameters (33) were computed throughout all simulations to monitor the time evolution of the solid/liquid interfaces.

In the simulations of diamond/liquid coexistence and BC8/liquid coexistence, we used 128- and 256-atom cells, respectively, with half of the atoms initially in a solid state and the other half in a liquid state. For selected thermodynamic conditions, we also repeated calculations of diamond/liquid coexistence with 256-atom cells (i.e., doubled in the direction perpendicular to the solid–liquid interface) to assess errors introduced by finite size effects on the calculated melting lines. In addition, we compared two-phase simulations carried out with finite and zero electronic temperature for selected points in phase space; at pressures >1,000 GPa, we assessed the error of the pseudopotential approximation by repeating some of our simulations using a pseudopotential with a smaller cut-off radius (decreased from 1.4 to 1.14 a.u. and, consequently, a higher plane wave cutoff, increased to 70 Ryd). We used *s, p* nonlocal, and *d* local potentials. In contrast to the low-pressure regime, in which the bonding is dominated by *s*–*p* hybridization, at extreme pressures, an accurate description of the *d* channel is important because the coordination of the liquid becomes >4 and configurations other than tetrahedral are present in the molten phase.

Finally, we determined the diamond/BC8 phase line using two independent approaches: (*i*) by combining results from two-phase simulations with equation of state calculations near the diamond/BC8/liquid triple point and (*ii*) from calculations of the free energies of diamond and BC8 based on the quasiharmonic approximation and including quantum ionic motion.

## Discussion

**Phase Boundaries.** The computed melting lines of diamond and BC8 are shown in Fig. 1. The internal consistency of our simulations was verified by computing the slope of melting lines from the Clapeyron equation (*dP*/*dT*_{m} = Δ*H*/*T*_{m}Δ*V*), where the specific volume (Δ*V*) and enthalpy (Δ*H*) are obtained from one-phase simulations, whereas the melting temperature (*T*_{m}) at which these quantities are computed is taken from the two-phase results (Fig. 1).

A detailed error analysis shows that at high-pressure (>500 GPa) errors introduced by finite size effects and by neglecting the finite electronic temperature (*T*_{e}) are of opposite sign. The melting temperature of diamond increases by ≈200 K in going from 128- to 256-atom simulation cells (at *P* = 1,000 GPa), but it decreases by ≈500 K when considering *T*_{e} ≈ 8,000 K in our calculation with 128-atom simulation cells (at 500–750 GPa and 8,000 K). Introducing electronic temperature likely reduces the directional character of the interatomic forces in the liquid, thus increasing the number of effective degrees of freedom of the ions, which, in turn, results in an increase of the ionic entropy and, hence, a reduction of the free energy of the liquid relative to the that of the solid. In contrast, the introduction of an electronic temperature has little or even a null effect on the insulating solid. In addition, we found a decrease of the melting temperature of ≈200 K at 1,000 GPa on decreasing the pseudopotential cutoff radius; this effect diminishes rapidly at lower pressures.

Although we estimate our numerical errors to be approximately ±150/250 K at fixed pressure, the error introduced by DFT is more difficult to quantify. We note that in cases for which the electronic charge density tends to be more homogeneous in the liquid phase than in the solid, the generalized gradient-corrected approximation (GGA) is expected to favor the former phase and, thus, to underestimate the meting temperature, which is the case, for example, for Al, Si, and LiH, and it is likely to be so for diamond and BC8. Although quantitative estimates for carbon are difficult to make in the absence of experimental data and calculations beyond the GGA, it is worth noting that, for Si, the melting temperature at ambient conditions is underestimated by ≈10%, and this error is expected to decrease as pressure is increased (31, 36).

The melting line of diamond obtained in our calculations is consistent with that reported by Grumbach and Martin (21), although the triple point found here is at a lower pressure than one might infer from their results. Our predicted diamond melting temperatures are also lower than those obtained in recent calculations with semiempirical potentials (26), especially in the high-pressure region. The reason for this discrepancy may come from the data set to which the empirical potentials have been fitted, which only contains low pressure, mostly tetrahedral, structures. We note that the gradual change in coordination number found here is the key for understanding the changes in the volume of the liquid relative to the solid as pressure is increased, and it is responsible for the existence of a maximum on the melting line, as we discuss below.

As apparent from Fig. 1, both diamond and BC8 melting lines have maxima slightly >8,000 K, at ≈450 and 1,400 GPa, respectively. The physical origin of these maxima stems from the changes that covalent interactions undergo upon compression. A signature of these changes can already be found at low temperature. Fig. 2 shows the zero temperature phonon dispersion curves in diamond as a function of pressure. As reported in ref. 37, acoustic branches of the phonon dispersion are rather flat near the Brillouin zone boundaries. In addition, the frequencies of the transverse acoustic branches near the L and X symmetry points have a maximum as a function of pressure. Although their softening above ≈250 GPa cannot destabilize the diamond structure at low temperature, it points at a loss of stability of the *sp*^{3} hybridization as a function of pressure. Similar effects appear also in the liquid phase at high temperature, where they are responsible for more dramatic changes in the liquid structure upon compression.

To analyze the structure of the molten phase, we have computed the pair correlation function and the coordination number of the liquid at different pressures (Fig. 3) close to the melting line. Unlike Si and Ge, we find that liquid carbon retains a partial covalent character in its bonding up to extreme pressures. We also find that the coordination of the fluid is lower than that reported in previous *ab initio* works (21). The changes taking place in the liquid, which are responsible for the maxima in the melting lines, are illustrated in Fig. 3, where we report the number of differently coordinated sites as a function of *P*. The fraction of fourfold coordinated atoms has a maximum near 400 GPa, whereas that of fivefold coordinated atoms peaks between 1,000 and 1,500 GPa. These maxima correspond to the maxima on the melting line displayed in Fig. 1.

The intersection of the diamond and BC8 melting lines yields the location of a solid/solid/liquid triple point at 7,445 K and 850 GPa. The phase boundary between the two solid phases close to the triple point can be calculated from the Clapeyron equation with the use of the computed slopes of the melting lines and equilibrium volumes of the three phases: $$mathtex$$$$mathtex$$[1] Here, subscripts D, B and L indicate diamond, BC8, and liquid phases, respectively, and all quantities are evaluated at the triple point. With this approach, difficulties associated with the explicit computation of free energies are avoided and the results do not rely on the harmonic approximation, which may not be valid at high temperature.

We find that the BC8/diamond phase boundary has a negative slope (Fig. 1) near the triple point. The difference between our BC8/diamond slope and that of ref. 21 is not too surprising, given that the phase stability presented in ref. 21 came from estimates based on the Lindeman criterion. The stability of BC8 relative to diamond is apparently enhanced at high temperature, possibly as a result of additional degrees of freedom in bonding angle configurations (as geometrically described in ref. 38).

At low temperature, the diamond–BC8 boundary is evaluated by carrying out free energy calculations in the quasiharmonic approximation. We have computed the full phonon dispersion relations for both the diamond and BC8 phases (39) and obtained the phonon energy and entropy by integration of the specific heat over temperature. The coexistence line between BC8 and diamond is then determined by matching the values of the Gibbs free energies, i.e., *G*_{D}(*T*_{coex}, *P*_{coex}) = *G*_{B}(*T*_{coex}, *P*_{coex}). Our results presented in Fig. 1 show a phase boundary that starts at 1,075 GPa at zero temperature and has a negative pressure versus temperature slope at finite temperature. We note that zero-point energy effects were taken into account and that they reduce the diamond/BC8 transition pressure by 62 GPa at zero temperature relative to a classical total energy calculation.

Finally, we note that by extrapolating the low temperature diamond/BC8 transition line, we obtain an accurate matching with the transition line near the triple point, as determined by Eq. **1** (Fig. 1.) The excellent agreement between the slopes of the phase boundaries obtained with three different methods [(*i*) free-energy matching, (*ii*) two-phase simulation for the location of the diamond/BC8/liquid triple point, and (*iii*) diamond/BC8 slope calculated with Eq. **1**] is evidence for the accuracy and consistency of the computed phase boundaries.

**Changes of Electronic Properties upon Melting.** Finally, we turn to the discussion of electronic properties of carbon under extreme conditions. It is well known that, at zero temperature, the band gap of diamond increases with pressure, and this result is reproduced in our calculations. The GGA approximation is known to underestimate the gap, compared with experiment (40, 41) (5.45 eV) and calculations with many-body corrections (5.6 eV) (42). To study whether diamond metalizes before melting, we have investigated the gap determined by the single-particle energy difference between the highest occupied orbital and the lowest unoccupied orbital of the solid as a function of *P* and *T*. Snapshots from 64-atom molecular dynamics were selected, and electronic structure calculations of the gap and the electrical conductivity [using the Kubo formula (43)] were carried out with dense *k*-point meshes (up to 6 × 6 × 6).

The calculated electronic gap of diamond decreases steadily with temperature at constant pressure [this was shown also to be a trend at ambient pressure and low temperature by experimental work (44)]; however, the zero temperature DFT–GGA gap (ranging from 4.4 eV at P = 0 to 6.9 eV at *P* = 1,100 GPa) is so large that temperature effects are not sufficient to close the gap before melting occurs. For example at *P* = 1,000 GPa, the gap reduces from 6.8 eV at *T* = 0 to 3.6 ± 0.5 eV at the melting temperature (6,750 K).

In our simulations the diamond DFT gap remained finite even at temperatures where the solid is likely to be superheated (i.e., metastable.) Frequency-dependent conductivity calculations confirm that the energy differences between the highest occupied orbitals and the lowest unoccupied orbitals correspond to the optical gaps. The liquid, instead, shows a continuous density of states at the Fermi level and a finite conductivity at zero frequency. Thus, we conclude that diamond remains an insulator in the solid phase, and an insulator-to-metal transition occurs only upon melting. Fig. 4 shows the spread of the maximally localized Wannier functions (45) as a function of pressure at fixed temperatures. The marked difference between the spread found in diamond and the fluid attests to the qualitative difference between their electronic properties.

Our results are consistent with recent shock-wave experiments (15), during which a transition was found to a conducting phase (revealed by a reflectance change) along the Hugoniot upon shock compression at approximately 600–700 GPa. The diamond Hugoniot crosses our calculated melting line at *P* = 720 GPa (Fig. 1), and, as in the case of the experiment, we predict a transition from an insulating solid phase to a conducting liquid phase.

The gap of BC8 has a very different behavior than that of diamond; it is much smaller in the region where BC8 is stable (see the corresponding maximally localized Wannier function spreads in Fig. 4). The DFT–GGA gap of BC8 is <0.45 eV at *T* = 0 at the diamond/BC8 transition point, and it slightly decreases with pressure. This small gap makes it difficult to assert a prediction on the metallic nature of this phase at finite temperature. If disorder reduces the gap at the same rate as in the case of the diamond phase, then BC8 is likely to become metallic in the solid phase at relatively low temperature before melting.

In conclusion, we have reported a first-principles computational study of the phase diagram of carbon, and we have determined solid/liquid and solid/solid phase boundaries up to ≈2,000 GPa and ≈10,000 K. Our results provide a consistent description of elemental carbon in a broad range of temperature and pressures, and a description of its electronic properties is given within the same framework. The diamond/BC8/liquid triple point is found to be at lower pressure than previously thought; in particular, this point is close to recent estimates (46–48) of the Neptune and Uranus core conditions. Although some old estimates locate the core conditions in the diamond stability region determined here (7,000 K, 600 GPa) (46), newer estimates point at conditions were the liquid is stable instead (8,000 K, 800 GPa) (47, 48). Therefore, our data call for a partial revision of current planetary models, especially in the planetary core regions. Overall, this work provides constraints to the carbon equation of states and it may help to interpret future experimental work.

## Methods

We carried out electronic structure calculations and first-principles molecular dynamics by using DFT in the GGA; we used normconserving pseudopotentials to describe the interaction between core and valence electrons, and we expanded single-particle orbitals in plane waves.** Simulations were carried out within the adiabatic Born–Oppenheimer approximation, in which the electronic ground state is computed self-consistently at each ionic step. The adiabatic approach is efficient for handling the small or null electronic band gap of the liquid phase of compressed carbon, and it allowed us to use ionic time steps up to 10 a.u.

## Acknowledgments

We thank Francois Gygi, Tadashi Ogitsu, Eric Schwegler, Xiaofei Wang, and Roberto Car for many useful discussions and Roger Falcone for constant help and support and for a critical reading of the manuscript. This work was performed under the auspices of the U.S. Department of Energy at the University of California/Lawrence Livermore National Laboratory under Contract W-7405-Eng-48. S.A.B. was supported by the Natural Sciences and Engineering Research Council of Canada.

## Footnotes

↵§ To whom correspondence should be sent at the present address: Department of Chemistry, University of California, Davis, CA 95616. E-mail: gagalli{at}ucdavis.edu.

Author contributions: A.A.C. and G.G. designed research; A.A.C. and S.A.B. performed research; A.A.C. and S.A.B. analyzed data; and A.A.C., S.A.B., and G.G. wrote the paper.

Conflict of interest statement: No conflicts declared.

Abbreviations: DFT, Density Functional Theory; GGA, generalized gradient-corrected approximation.

↵¶ BC8 has been experimentally shown to be a metastable phase of silicon (12) and germanium (13) and therefore proposed as a possible high-pressure phase of carbon.

↵∥ The discrepancy with previous reports on the stability of R8 stems from the numerical accuracy of the calculations that have appeared in the literature, which were not fully converged as a function of plane wave cutoff and

*k*-point sampling.↵** Periodic boundary conditions and a plane wave basis with a 50-Ryd energy cutoff and Troullier–Martins pseudopotential (49) were used in this work. Our

*ab initio*molecular dynamics simulations were carried out with the qbox code, written by F. Gygi (Lawrence Livermore National Laboratory). Static calculations of electronic properties and total energy convergence with*k*-points were obtained through the use of the abinit code [a common project of the Université Catholique de Louvain (Louvain-la-Neuve, Belgium), Corning Incorporated, and other contributors].

- Received October 6, 2005.

- Copyright © 2006, The National Academy of Sciences

Freely available online through the PNAS open access option.

## References

- ↵
- ↵
- ↵
- ↵
- Hubbard, W. B. (1981) Science 214.
**,**145–149. - ↵Benedetti, L. R., Nguyen, J. H., Caldwell, W. A., Liu, H., Kruger, M. & Jeanloz, R. (1999) Science 286.
**,**100–102.pmid:10506552 - ↵
- ↵
- ↵Kuchner, M. J & Seager, S. (2005) Astrophys. J., in press..
- ↵Yin, M. T & Cohen, M. L. (1984) Phys. Rev. Lett. 50.
**,**2006–2009. - ↵
- ↵Hu, J. Z., Merkle, L. D., Menoni, C. S. & Spain, I. L. (1994) Phys. Rev. B Condens. Matter 34.
**,**4679–4684. - ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵Galli, G., Martin, R. M., Car, R. & Parrinello, M. (1990) Science 250.
**,**1547–1549. - ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵Kechin, V. V. (2002) Phys. Rev. B Condens. Matter 65, 052102–052105..
- ↵Zemansky, M. W. (1968) Heat and Thermodynamics (McGraw–Hill, New York), 5th Ed..
- ↵
- ↵
- ↵
- ↵
- ↵Schulz, M. & Weiss, H. (1982) Landolt-Börnstein Tables, Numerical Data and Functional Relationships in Science and Technology, ed. Madelung, O. (Springer, Berlin), Vol. New Series 17a..
- ↵
- ↵
- ↵
- ↵Papadopoulos, A. D. & Anastassakis, E. (1990) Phys. Rev. B Condens. Matter 43.
**,**9916–9923. - ↵
- ↵
- ↵Guillot, T. (1999) Science 286.
**,**72–77.pmid:10506563 - ↵
- ↵

## Citation Manager Formats

## Article Classifications

- Physical Sciences
- Physics