Main

Liquid carbon is difficult to produce in the laboratory3,12. It requires temperatures exceeding 4,000 K and pressures of at least several megapascals. In nature, these conditions are present in the interior of large planets such as the ice giants of our solar system, Uranus and Neptune, in which liquid carbon may contribute to the unusual magnetic fields of these planets4,13. Moreover, the equation of state of carbon is of substantial importance to estimate the composition of exoplanets from their observed mass and radius, in particular, for the highly abundant class of sub-Neptunes5. For technology applications, liquid carbon is an important transient state for the synthesis of several advanced carbon materials, such as carbon nanotubes7, nanodiamonds8,14 and Q-carbon15. Liquid carbon may also be key for the synthesis of the BC-8 phase of carbon, which has been predicted for decades at pressures beyond diamond stability16,17 but could not be realized experimentally so far despite extensive efforts18. At the same time, carbon is used in inertial confinement fusion experiments as an ablator material surrounding the deuterium–tritium fuel19. The experimental design that achieved an ignited fusion plasma at the National Ignition Facility20 relies on high-density carbon (diamond) that is subjected to shock compression just above melting in the initial phase of the implosion. This initial compression step is crucial for the subsequent fusion yield21, and a better microscopic understanding of liquid carbon under dynamic compression will help to design more efficient implosions, especially as more amorphous forms of carbon are considered as future ablator materials22.

Modelling the extreme conditions in which liquid carbon prevails is challenging. Planetary interior pressures of around 100 GPa and temperatures approaching 10,000 K result in energy densities that are of the order of those stored in carbon–carbon bonds. Transient chemical bonds persisting from the sp3-bonded diamond lattice are still expected to shape the structure of liquid carbon at these conditions, leading to tetrahedral coordination with four nearest neighbours on average, which is in contrast to simple liquids with icosahedral coordination with up to 12 nearest neighbours3. This complexity inhibits simple approximations and leaves first-principles simulations, usually based on density functional theory with molecular dynamics (DFT-MD), as the sole reliable method to provide predictive abilities. However, even DFT-MD requires assumptions, such as the choice of the exchange-correlation potential, and its computational limitations in system size and simulation times may not capture all relevant physics processes. There have been large discrepancies in the predictions of the melting curve of carbon, including deviations up to a factor of 2 in melting temperatures and fundamental differences in the slope of the melting curve for the diamond phase23. Machine-learning potentials based on DFT try to circumvent the scaling limitations by vastly extending spatial and temporal scales24; however, effects not covered by the training data are not necessarily expected to be captured by the scaled simulations. Although modern DFT-MD methods now seem to converge on high-pressure equilibrium phase diagram calculations with smaller variations25,26,27, these predictions remain to be tested experimentally.

Probing the structural properties of carbon and other low-Z fluids in static high-pressure experiments, for example, using diamond anvil cells, is challenging, and analogues such as glasses are used instead28. The high temperatures required to preserve the liquid state for prolonged periods lead to disintegration of the high-pressure sample confinement. For X-ray probing, the small scattering power of light elements often hampers resolving liquid structures above the background from surrounding material. In turn, most experimental approaches to study liquid carbon have used dynamic techniques such as electrical or optical flash heating and shock compression3. However, as in situ X-ray probing is difficult in these experiments, most results can up to now provide only indirect evidence for the presence of liquid carbon, and a detailed structure measurement has not yet been achieved. Velocimetry measurements in shock compression experiments can provide hints of melting through small changes in the slope of the shock Hugoniot curve29. Using pyrometry on decaying shocks, a temperature plateau was associated with melting, and an anomalously high specific heat in the liquid phase was interpreted to suggest a reconfiguration of atomic packing, from a partially bonded complex fluid to an atomic fluid between 10,000 K and 60,000 K (ref. 30). Measurements of electrical resistivity31 and optical reflectivity32 provided information on the conductivity of liquid carbon before X-ray sources with sufficient flux became available to attempt measurements of structural properties. X-ray absorption spectroscopy at soft X-ray sources using femtosecond flash heating, in which the dynamics can be benchmarked by extreme ultraviolet reflectivity33, provided evidence for π and σ bonds in liquid carbon and some indirect information on structure based on theoretical modelling34,35. However, owing to the short timescales of these studies, the investigated states had not reached thermal equilibrium, and the theoretical methods applied require experimental benchmarking. A direct X-ray scattering measurement of the atomic structure of liquid carbon on the nanosecond timescale realized by laser-driven shock compression could be achieved only at a few distinct points in k-space, which leaves substantial degrees of freedom for models36,37. Synchrotrons and hard X-ray free-electron lasers finally started to allow for the measurement of diffraction patterns from liquids in dynamic compression experiments reaching extreme conditions38,39,40,41,42. New facilities at the European XFEL (EuXFEL) set new standards in this direction43.

The experiments reported here were performed using the DiPOLE 100-X high-energy laser at the High Energy Density-Helmholtz International Beamline for Extreme Fields (HED-HIBEF) instrument of EuXFEL43,44. Figure 1 shows the experimental setup. The DiPOLE 100-X laser was used to drive shock waves into glassy carbon samples, which generate high-pressure states with simultaneous heating due to the shock-induced entropy increase. Being a diode-pumped laser system, DiPOLE 100-X features an energy and temporal pulse shape stability at a sub-per-cent level, which allows for highly reproducible shot-to-shot drive conditions. The bright X-ray pulses delivered by EuXFEL with a photon energy of 18 keV were used for in situ X-ray diffraction (XRD) from the shock-compressed sample to monitor the microscopic structure. A velocity interferometer system for any reflector (VISAR) was used to capture the shock dynamics and determine the pressure, together with the diffraction data (Methods). Figure 2 shows integrated lineouts of single-shot diffraction patterns acquired before shock breakout time for increasing pressures. Ambient glassy carbon shows amorphous features associated with sp2 bonds. At (76 ± 8) GPa, the sp2 signature vanishes, and we observe a partial transformation of the glassy carbon to crystalline diamond. Remnants of the amorphous structure are still present at these conditions, possibly because of temperatures not considerably exceeding the glass transition temperature45. This changes at (83 ± 9) GPa, at which the diamond peaks substantially intensify and sharpen in comparison with the lower pressure conditions. This is consistent with the formation of larger crystallites and the probed sample volume being nearly fully composed of diamond. At (106 ± 11) GPa, the crystalline features start to diminish, together with the appearance of broader liquid correlation peaks. We interpret this as a coexistence state between diamond and liquid carbon; these features are also present at (126 ± 12) GPa but with lower diamond content. At (160 ± 14) GPa, we observe a purely liquid state.

Fig. 1: Schematic of the experimental setup.
figure 1

Glassy carbon samples were subjected to shock compression with the DiPOLE 100-X drive laser. The microscopic structure is probed by a bright X-ray pulse of EuXFEL, and two area detectors collect the XRD patterns. The shock dynamics are captured by a VISAR.

Full size image
Fig. 2: XRD patterns of shock-compressed glassy carbon.
figure 2

The drive pressure increases from top to bottom, and the curves are shown with a constant offset of 10 a.u. The X-ray probing times of the lineouts are within 1 ns before shock release to probe mostly homogeneous conditions and avoid pressure release states. At ambient conditions, broad amorphous structures are present. Diamond formation is observed above about 76 GPa, the coexistence of diamond and liquid carbon from about 100 GPa and complete melting at about 160 GPa. The quoted pressure uncertainties are dominated by the error estimations of the shock velocity measurement (Methods). a.u., arbitrary units.

Source Data

Full size image

From the diffraction patterns, we can extract the static structure factor profiles S(k) for the covered scattering wavenumbers k, which allows for a direct comparison with ab initio simulations based on DFT-MD (Methods). In agreement with the simulations, we find a complex liquid with broad liquid correlation peaks forming around the positions of crystalline diamond (Fig. 3), which is compatible with transient bonds resulting in approximately fourfold coordination on average27. More simplistic models with higher coordination numbers, such as Lennard–Jones, result in a first correlation peak between the two observed, and such a structure37 is inconsistent with our observation. Calculating the Fourier transform of S(k) allows us to determine the radial distribution function g(r) and with that the first and second coordination numbers of the liquid carbon state. For the case of complete melting, Fig. 4a shows a range of experimental reconstructions with varying kcutoff (Methods). The span of results is in good agreement with the DFT-MD simulation that fits best to the corresponding XRD pattern. Again, a simple Lennard–Jones liquid does not match. Although the height and width of the first correlation peaks vary between the chosen cutoffs, the integrated area underneath the peaks and thus the extracted first and second coordination numbers remain rather constant. This is reasonable because the structural information is encoded in the XRD pattern, in which simple liquids with high coordination numbers are incompatible. For complete melting, we find a first coordination number of 3.78 ± 0.15 and a second coordination number of 17 ± 2, which is in agreement with several DFT-MD predictions of the bonded liquid27,46,47,48 and our DFT-MD simulations (first coordination number of 3.66 ± 0.05).

Fig. 3: DFT-MD fits of the experimentally obtained liquid structure.
figure 3

The different datasets are shown with a constant offset of 1 between the curves. For the data with solid–liquid coexistence, the crystalline peaks have been omitted for fitting, and the liquid structure includes thermal diffuse scattering from diamond. The grey-dashed curve shows a hypernetted chain calculation of the liquid structure assuming a Lennard–Jones potential at similar conditions, which results in a first coordination number of around 11 and shows the fundamental difference to the observed structure, which has approximately fourfold coordination.

Source Data

Full size image
Fig. 4: Radial distribution function and phase diagram.
figure 4

a, Experimental reconstructions of the radial distribution function for complete melting at (160 ± 14) GPa and comparison with the simulation of the best fit to the corresponding XRD pattern and a simple Lennard–Jones liquid at 3.8 g cm−3 and 7,000 K. b, Our results for the fluid–solid coexistence and the pure fluid shown in context with the predictions of the carbon melting line and coordination numbers for liquid carbon. The temperature uncertainties shown arise from the error estimations of the DFT-MD fits, and the pressure uncertainties correspond to the values quoted in Fig. 2.

Source Data

Full size image

By fitting with DFT-MD, we also infer estimates for the temperature and the density of the state reached within the probed volume. For the cases with solid–liquid coexistence, we fit a combination of liquid structure and thermal diffuse scattering of diamond, which can also be obtained from DFT (Extended Data Fig. 4). Theoretical predictions expect the correlation peaks of S(k) to move to higher k with increasing density and broaden for higher temperatures27 (Extended Data Fig. 3). Hence, we can provide experimental benchmarks for the melting temperature, the volume change from solid to liquid and the associated latent heat through the liquid structure at melting. The pressures extracted from the DFT-MD fits in Fig. 3 match with those obtained experimentally from VISAR and XRD reasonably well within the measurement uncertainty. Only for the highest pressure case, there is a small discrepancy, but still within the margins, because the uncertainties in temperature (±200 K) and density (±0.05 g cm−3) of the DFT-fit result in a pressure error of  around 8 GPa from the simulations. In the following, we use the pressures inferred from VISAR and the density from XRD. At (106 ± 11) GPa, the positions of the crystalline diamond peaks result in a density of 3.91 g cm−3, whereas the density of the liquid is best matched by 3.62 g cm−3 with a temperature of 6,557 K. The volume change of about (7 ± 1)% between the two phases is in reasonable agreement with DFT-MD predictions along the melting curve in this pressure regime (8% at 104 GPa and 6% at 183 GPa; ref. 46). Our data are also compatible with the pressure–temperature slope in ref. 25 of about 11.2 K GPa−1. Using this value together with the experimentally determined volume change in the Clausius–Clapeyron relation, the entropy of melting and the latent heat are estimated to be about 20 J mol−1 K−1 and about 130 kJ mol−1, respectively. The purely fluid state at (160 ± 14) GPa is best matched by a density of 3.79 g cm−3 and a temperature of 7,314 K.

Figure 4b shows these findings in the phase diagram of carbon with different DFT-MD predictions of the melting curve25,46 and first coordination numbers27,46. The temperature values of our results are inferred from the comparison of the experimental and simulated structure curves. Thus, we assume that the structural representation of liquid carbon in our simulations is correct, which is corroborated by the excellent agreement with our diffraction data (for example, in comparison with the simple Lennard–Jones liquid). Although our deduced temperatures are not free from assumptions, we determine the liquid structure at melting with high precision, which will be a valuable benchmark for all future simulations of the melting transition of carbon in this regime. In general, we find high consistency with the more recent DFT-MD equation-of-state calculations in ref. 25, whereas other simulations46 would require higher temperatures for melting, which is inconsistent with our results. Several melting curves based on more approximate models49,50 do not match our data. Our measurements are expected to achieve thermodynamic equilibrium at the highest temperatures, given the robust diamond formation before melting, the exceeding of the glass transition temperature for glassy carbon and the observed consistency with equilibrium melting models, with that in best agreement25, also having similar simulation settings to those used in our DFT-MD calculations. Moreover, we show the predicted coordination numbers of distinct DFT-MD simulations in refs. 27,46. Again, we find very good agreement with our measurements. The higher values predicted in ref. 46 match slightly better, but ref. 27 is also compatible within the experimental and numerical uncertainties.

In conclusion, our pioneering experiments substantiate the view of liquid carbon as a complex liquid with approximately fourfold coordination at pressures of about 100 GPa. Overall, the structures predicted by modern DFT-MD simulations are in agreement with the experimental data, which underlines the predictive power of this method at pressures around 100 GPa and elevated temperatures. It should be noted that all curves shown here were collected as single-shot events with a repetition rate in the range of minutes. However, both the drive laser and the X-ray probe can run at 10 Hz. Thus, future experiments can obtain higher precision by the accumulation of data and determine the liquid structure of a plethora of compounds made out of light elements at extreme pressure and temperature conditions. This could lead to substantial progress in models for the interior and the evolution of icy giant planets, the classification of exoplanets made of similar constituents, defining new processes of materials synthesis through extreme conditions, and designing the best ablator materials for inertial confinement fusion.

Methods

Laser-driven shock compression and experimental pressure determination

The samples were subjected to shock compression using frequency-doubled pulses of the DiPOLE 100-X laser (515 nm wavelength) and a phase plate that produces a smoothed laser spot of approximately 250 μm in diameter. XRD was recorded by two Varex 4343CT flat panel detectors. For more details on the experimental configuration and geometry, see ref. 43. Quasi-flattop pulses of 10 ns with 35 J maximum energy and 5 ns with 28 J maximum energy resulted in steady shock compression waves. The thickness of the samples was either 60 μm or 92 μm. For each drive condition, we obtained a time series with intervals between different X-ray probe timings of 0.5–1 ns, and the data used for XRD fitting were within the last nanosecond before the shock release for the 6–10 ns of total transit time. We chose the timing before release to avoid the large density–pressure gradients afterwards. Signal from any remaining cold material can be accounted for by subtracting the ambient pattern that is recorded for each sample before the laser shot, with a scaling of the probe timing relative to the recorded shock release. When diamond peaks are present, the density of the crystallites can be determined with high precision and shows constant density within the measurement uncertainty of 0.01 g cm−3 as long as the shock propagates inside the sample (Extended Data Fig. 1). Thus, the assumption of a planar steady shock is reasonably justified (Extended Data Fig. 2). Once the shock releases on the sample rear side, the density of the diamond crystallites approaches ambient density after a few nanoseconds and reaches even lower values afterwards because of the residual high temperatures. The different time series do not show substantial effects of X-ray preheating, as the XRD features associated with the weak sp2 bonds diminish proportionally to the shock distance travelled and are not markedly affected right after the impact of the drive laser. As glassy carbon is not transparent to optical light, the VISAR system can determine only shock transit times, and the assumption of a steady shock propagation provides the shock velocity. In situ XRD allows for determining density as long as diamond and/or liquid carbon are present. With the obtained density and shock velocity, the shock pressure PS can be determined by the Rankine–Hugoniot relations51:

$${P}_{{\rm{S}}}={\rho }_{0}\left(1-\frac{{\rho }_{0}}{{\rho }_{{\rm{S}}}}\right){v}_{{\rm{S}}}^{2}.$$
(1)

The uncertainties of the resulting pressures are dominated by the uncertainty of the shock velocity due to the quadratic scaling. For the lower pressure conditions, there is very good overlap with gas gun Hugoniot measurements on glassy (vitreous) carbon up to 85 GPa (ref. 52; Extended Data Table 1 and Extended Data Fig. 1). For higher pressures, there are no existing Hugoniot data in the literature.

Liquid diffraction analysis

The total scattered X-ray intensity I(k) in our XRD patterns, which are corrected for transmission through aluminium filters in front of the detectors, is given by

$$I(k)=\alpha ({f}^{2}(k)S(k)+{I}_{{\rm{inc}}}(k)),$$
(2)

where α is a scaling factor from atomic units to measured counts on the detector. The static structure factor S(k) was then determined by using the atomic form factors f(k) and incoherent scattering functions Iinc(k) for carbon as tabulated in ref. 53. The normalization of the experimental S(k) curve was included into the fitting procedure to the DFT-MD simulations by minimizing the function

$${\left(\frac{I(k)}{\alpha }-{p}_{1}{I}_{{\rm{inc}}}(k)-{p}_{2}{f}^{2}(k)S{(k)}_{\rho ,T}^{{\rm{DFT}}}\right)}^{2},$$
(3)

where the ratio of the scaling parameters for incoherent and coherent scattering, p1 and p2, respectively, resulted in a constant value for all analysed datasets as the incoherent scattering is not affected by structural changes. For fitting the DFT-MD simulations, we cut the measured XRD pattern at kcutoff = 7 Å−1, because at higher angles, the applied filter and geometry corrections for the detector sensitivity become more severe43 and would increase the uncertainty of our analysis. The radial distribution function

$$g(r)=1+\frac{1}{2{{\rm{\pi }}}^{2}nr}{\int }_{0}^{{k}_{{\rm{cut-off}}}}[S(k)-1]\,k\sin (kr){\rm{d}}k$$
(4)

with atomic number density n was obtained from S(k) by linear extrapolation of the experimentally obtained curve for (160 ± 14) GPa to k = 0. The coordination numbers

$${N}_{{\rm{C}}}=4{\rm{\pi }}n{\int }_{{r}_{1}}^{{r}_{2}}{r}^{2}g(r){\rm{d}}r$$
(5)

were then obtained by integrating between the corresponding minima of g(r) at r1 and r2. For the first coordination number, r1 = 0 and r2 is given by the first minimum of g(r). For the second coordination number, r1 is the first minimum of g(r) and r2 is the second. The provided error estimations of the coordination numbers were determined by reasonable variations of kcutoff between 6.8 Å−1 and 7.5 Å−1 for calculating g(r) and the density uncertainty.

DFT-MD simulations

All DFT-MD simulations in this work were performed with the Vienna ab initio simulation package (VASP)54,55,56. The electronic and ionic parts were decoupled by the Born–Oppenheimer approximation and, for fixed ion positions, the electronic problem was solved in the finite temperature DFT approach57 using a projector-augmented wave pseudopotential (labelled PAW_PBE C_h)58,59 and the Perdew–Burke–Ernzerhof functional60 for the exchange-correlation contribution. A 2 × 2 × 2 Monkhorst–Pack sampling61 was used for the k-space of a simulation box with 64 atoms, and a plane wave cutoff energy of 1,000 eV was used, in which the resulting structure curves show excellent agreement with calculations using 216 atoms, which have been performed for selected parameters. The molecular dynamics time step was t = 0.2 fs. A Nosé–Hoover thermostat was used with the Nosé mass set to 0.5 am, corresponding to 67 time steps (0.37 × 1015 Hz). The number of the considered MD steps (for the calculation of the structure factors) is generally in the range of 10,000. To perform the least-square fitting to the diffraction data, the static ion–ion structure factor was computed on a density and temperature grid ranging from 3.6 g cm−3 to 3.9 g cm−3 and from 6,000 K to 8,000 K with four and five density and temperature increments, respectively. The simulations provide the three-dimensional particle density distribution n(r, t), and by Fourier transform and averaging, the structure factor S(k) can be obtained37. Close to the melting line, finite size effects that can seed crystallization features were circumvented by training a high-dimensional neural network potential. The forces and energies predicted by DFT were learnt by a Behler–Parrinello high-dimensional neural network potential62 implemented in the n2p2 software package63,64. For more details on this method, see ref. 65. The temperature control in all molecular dynamics simulations was performed by a Nosé–Hoover thermostat66,67. A finer resolution of static structure factors in the density–temperature plane was achieved by computing a neural-network-based three-dimensional representation of S(k) in the k–ρ–T space as suggested in ref. 68. The neural network is a feedforward neural network with linear connections and layers with 3, 64, 1,024, 1024, 1,024, and 1 neurons and ReLU activation between the layers. A benchmark of the interpolation grid result at 3.7 g cm−3 and 8,000 K with a DFT-MD simulation using a box with 216 atoms at these conditions is shown in Extended Data Fig. 5.