Introduction

Liquid polymorphism has attracted increasingly scientific interests since the discovery of multiple liquid phases distinguished by density1,2,3. Liquid-liquid phase transition (LLPT) has been suggested for a dozen systems with a locally tetrahedral molecular structure4,5,6,7. For example, Katayama et al.8 adopted an in situ x-ray diffraction technique to observe the LLPT in phosphorus successfully, involving an abrupt, pressure-induced structural change between two distinct liquid forms. Recently, Beye et al.9 have confirmed the existence of the LLPT in silicon through monitoring the electronic structure of silicon changes after strong laser excitation.

Similar to phosphorus and silicon, carbon should also undergo a LLPT owning to its ability to form various covalent bonds10. Unfortunately, it is almost impossible to obtain direct observations for a LLPT under the current experimental condition due to the high melting temperature of graphite. An indirect evidence for the LLPT in carbon was found when Togaya investigated the melting line of graphite by a flash-heating experiment under high pressures11. Discontinuity in the slope of pressure-temperature melting line at about 5.6 GPa suggested the possible phase transition in liquid carbon, based on a fit of the melting line with two straight lines. However, this indirect evidence still had some uncertainties to resolve the issue of a LLPT in liquid carbon12. As direct experiments are difficult, a lot of theoretical work has been conducted to explore such transition in liquid carbon. Glosli and Ree predicted the existence of LLPT between twofold and fourfold coordinated liquid carbon13, using classical molecular dynamics (MD) simulations based on the Brenner's bond-order potential. However, a subsequent first-principles investigation by Wu et al.12 ruled out the occurrence of such a transition in the same temperature and pressure range. Their ab initio results indicated a continuous evolution from a low-density sp/sp2-like liquid to a high-density sp2/sp3 liquid as a function of pressure. On the basis of density-functional based MD simulations, Ghiringhelli and his colleagues proposed that14, with increasing pressure, liquid carbon experience a gradual transformation from a three-fold liquid to a “diamond-like” liquid. Although some work has been done on the structure of liquid carbon, it has long been suspected that carbon might exhibit an LLPT.

In the last decades, low-dimensional carbon materials with a perfect honeycomb structure have attracted considerable attention on the theoretical research and the potential applications at the intersection between chemistry, physics and material sciences owing to their exotic properties. Especially since the experimental preparation of graphene15, it has led to an explosion of interest in two-dimensional (2D) carbon materials due to the excellent electronic, mechanical, thermal and optical performance. However, the research regarding the 2D carbon in liquid form is relatively limited till now. In the previous studies, some other quasi-2D liquids16,17,18 and their phase transition19,20,21,22,23, have been proved to exhibit some special characters different from their bulk counterparts. Then, what about the quasi-2D liquid carbon and its phase transition? Instead of focusing on the phase transition versus pressure in bulk carbon in previous studies12,13,14,24, we report the cooling process of quasi-2D liquid carbon at a constant plane pressure of 5 GPa. Our study indicates a LLPT in liquid carbon, closely related to the followed liquid-solid phase transition (LSPT). This work will further advance the general understanding of structure transitions in liquid carbon.

Results

As shown in Fig. 1(a), the density changes from ρ = 1.35 g/cm3 to ρ = 2.8 g/cm3 during the cooling process. For the high-temperature liquid in stage I (T > 7937 K), the density presents linear change with the temperature, followed by a non-linear dependence in the stage II and subsequently shows a linear relation again in the stage III. Moreover, there is an obvious inflection point after the stage III, which is defined as the starting point for the stage IV. The change of density ρ as a function of temperature T can be fitted as follows:

To further distinguish the segmentation of cooling process, we investigate the coordination fraction as a function of temperature as shown in Fig. 1(b). In this study, the coordination of each atom counts neighbors within a radius of 2 Å. It is reported that the coordination distribution is closely related to the geometry of the local structure, which is extremely useful in the analysis of the various covalent bonding structures13,14. For carbon, the threefold coordinated atoms correspond to the centers of 2D Y-shaped structures, while the twofold coordinated atoms tend to form atomic chains. Our results in Fig. 1(b) predict a marked switching of dominant coordination from two to three in the stage II. In addition, it is clear that the system at 6466 K is fairly diffusive as shown in Fig. 2(a) and hence should be characterized as a liquid. Hence, the stage II is regarded as the stage of LLPT from twofold coordinated liquid to threefold coordinated liquid. In the stage I, the coordination fractions change slowly with the temperature in a linear form and start to change non-linearly from the stage II. The Arrhenius plot of the diffusivities presented in Fig. 2(a) also shows a strongly linear dependence in the stage I. However, with the decrease of temperature, the diffusion coefficient begins to deviate from the linear dependence on the temperature, indicating the changes of the structure. The coordination distribution and the diffusion coefficient both show different change rules between the stages I and II. The break point between the stages I and II obtained by coordination distribution and the diffusion coefficient curves is 7937 K which agrees well with the segmentation based on the change of density. After the stage III, the diffusion coefficient becomes rather small and the structure starts to changes little with the decrease of temperature, which indicates that the liquid carbon has transformed into solid after the stage III. As a result, we define the stage III and IV as LSPT stage and solid stage, respectively. The heat capacity plotted in Fig. 2(b) reaches its peak value in the LSPT stage and then decreases abruptly. Comparing with the sudden dropping after the solidification, the heat capacity before the LSPT stage starts to increase ever since the stage II and take a longer time to reach the peak value. This kind of asymmetry in the change of heat capacity before and after solidification also indicates the existence of the LLPT. The change of surviving structures can also be used to distinguish the four stages, which will be discussed in the following section.

Figure 1
figure 1

Structure transition during the cooling process.

(a) Density as a function of temperature. (b) Coordination fractions versus temperature.

Full size image
Figure 2
figure 2

(a) Arrhenius plot of the diffusivity. The diffusivities are linearly related to temperature dependence in the stage I. The inset shows the mean squared displacement at different temperature. (b) Heat capacity as a function of temperature.

Full size image

The inset in Fig. 3(a) shows the pair-correlation function g(r) for the carbon film at different temperature. An interesting phenomenon is observed that, as the compactness increases with temperature, so does the position of first peak, completely contrary to normal substance where the increase of density will pack the neighboring atoms more close. For carbon film in this study, the position of first peak gradually increases from 1.35 Å to 1.40 Å in the liquid stage as the temperature goes down and then to 1.42 Å after the LSPT. The reason is that the bonds in twofold coordinated structures are shorter than that in the threefold coordinated structures25. As for the second peak in g(r), its position gradually decreases as the temperature drops, due to the reduction of bond angle. At 4000 K, the second peak position decreases to 2.48 Å, nearly equal to the second nearest neighbor (2.46 Å) in graphite. Fig. 3(a) represents the calculated angular correlation functions g(θ) within a radius of 2 Å. It is found that the atomic chains do not show a strictly straight feature, as the bond angles for first neighbors in chain dominant liquid show a flat peak ranging from 110° to 130° rather than around 180°. As the liquid carbon contains more and more Y-shaped structures but less twofold coordinated structures with the decrease of temperature, the peak of g(θ) for first neighbors gets higher and narrower and finally takes shape at 4000 K with a peak position of approximately 115°. To test a longer range order, another g(θ) for the neighbors within 3 Å corresponding to the second minimum of g(r) is provided in Fig. 3(b). It is found that the liquid carbon only forms two distinct peaks at 30° and 60°, while the other larger angles nearly have the similar content. However, the situation is quite different after the LLPT that several new peaks begin to appear every 30 degrees in the range of 90° to 180° and become more and more evident with the decrease of temperature, which indicates the generation of more ordered structures.

Figure 3
figure 3

A series of measuring means to detect the structure changes at different temperatures.

(a) and (b) represent the calculated angular correlation functions g(θ) within the radiuses of 2 Å and 3 Å, respectively. Inset in (a) shows the pair-correlation function g(r). Inset in (b) shows a simplified model of carbon island to explain the bond angle distribution. The green and red arrows represent C-C bond within 2 Å and 3 Å, respectively.

Full size image

To explore the microscopic details of the structure evolution during the cooling process, simulation snapshots at different temperature are investigated, as displayed in Fig. 4. It can be found that the structure has been dramatically transformed during the cooling process. In the stage of LLPT, the liquid structure is transformed from the twofold coordinated liquid to the threefold coordinated liquid. In this stage, ring structures are formed by the atomic chains in four different ways, as shown in Fig. 5. The self-kinking mode mostly appears in the high-temperature liquid and rarely appears in a more compact phase with the decrease of temperature. The most common mode of ring formation is induced by threefold coordinated structures. As shown in Fig. 5, there are three different inducing modes distinguished by the number of threefold coordinated structures. It can be found that the Y-shaped threefold coordinated structures provide strong nodes to induce the chains to form rings by coalescence. If new rings are induced by the threefold coordinated structures in old rings, the atomic rings will aggregate to form carbon islands. However, it is worth noting that the rings and islands in this stage are not stable and are easy to break into atomic chains inversely. With the further decrease of temperature, if the size of islands exceeds a critical value, the islands can stable exist as a nucleus, which indicates the beginning of LSPT. The LLPT in liquid carbon can be viewed as the preparation stage of the LSPT. As shown in Fig. 4, the liquid carbon is simultaneously nucleated at individual islands. Subsequently, the islands begin to grow outwards by incorporating the discrete chains and rings. Eventually, the carbon islands will encounter each other and combine to form a complete polycrystalline film. In the stage of LSPT, the thickness of the carbon film is found to suddenly drop, accompanied by an obvious layering phenomenon, as shown in Fig. 6. This should be attributed to that the nucleus can stably exist and grow in the LSPT stage and at the same time more and more disordering chain-like structures transform into the ordering carbon islands parallel to the plane directions.

Figure 4
figure 4

Nucleation and growth in the cooling process.

Brown: threefold coordinated structures; Cyan: twofold coordinated structures; Blue: onefold coordinated structures. (a) 7500 K; (b) 7100 K; (c) 6800 K; (d) 6466 K; (e) 6200 K; (f) 5886 K. Here, only the upper half is displayed.

Full size image
Figure 5
figure 5

Different modes of ring formation.

(a) Self- kinking mode. (b) Inducing mode by one threefold coordinated structure. (c) Inducing mode by two threefold coordinated structures. (d) Inducing mode by three or more threefold coordinated structures. Red: Atoms finally forming a ring. Blue: Initial threefold coordinated structures.

Full size image
Figure 6
figure 6

Thickness of the carbon film as a function of temperature.

(b) Density distribution in the vertical direction.

Full size image

As the threefold coordinated structures are important in the formation of rings, next their surviving will be further discussed, as shown in Fig. 7. The surviving time is defined as the duration of threefold coordinated structures with no change in its composition. In the stage I, there are only a few structures surviving longer than 1 ps. When the temperature decreases to the LLPT stage, the structures surviving for 1 ps start to increase sharply. In this stage, the structures with much longer surviving time also start to increase with the further decrease of temperature. An obvious turning point has been observed in the number of structures surviving for 1 ps between the LLPT and LSPT stages. This is because the surviving time continues to increase in the LSPT stage and more structures can survive for more than 2 ps. In the stage of LSPT, the number of threefold coordinated structures surviving for 2–9 ps reaches peak value sequentially and starts to decrease with the decrease of temperature. Simultaneously, the threefold coordinated structures surviving for longer than 10 ps increase sharply due to the nucleation and growth of islands in the LSPT stage. During the whole cooling process, the surviving time of threefold coordinated structures increase with the decrease of temperature, which provides essential condition to form rings and islands. The longer surviving time means that the liquid carbon in the next moment can inherit more threefold coordinated structures, which provides greater opportunities for the formation of rings.

Figure 7
figure 7

Surviving threefold coordinated structures at different temperature.

The surviving time ts is defined as the duration of threefold coordinated structures.

Full size image

It is should be noted that not all the carbon rings are hexatomic and there is a fair amount of other ring types. As shown in Fig. 8, the number of rings shows a strong dependence on the temperature. In the stage of LLPT, the sharp increase of threefold coordinated structures induce a large number of rings. Quantity variances start to emerge for different ring types in this stage and the 6-membered rings start to have a larger number than the other types of ring. During the LSPT, the ring number continues to increase sharply, dominated by the increase of 6-membered ring. In the following stage IV, the solid phase still seems unstable and the ring number continues to change, due to the different crystal orientation of every nucleus in the LSPT stage. When the growing islands encountering each other to form a complete polycrystalline film, the system must undergoes an optimization process to release the stress on the boundaries of islands. In this stage, the 6-membered rings increase but the other ring types decrease with the temperature except for some transient reverse transformations from 6-membered rings to 5/7-membered rings. These reverse transformations will effectively reduce the stress concentration in local uncomfortable topology pattern, which usually appear under great external tension26,27. After the reverse transformations, the optimization process will result in more 6-membered rings and more optimal structures. During the whole cooling process, it is found that the number of 5-membered and 7-membered rings changes synchronously.

Figure 8
figure 8

The number of rings as a function of temperature.

Full size image

Discussion

Different from previous studies on the LLPT versus pressure in bulk carbon12,13,14,24, we have used MD method to study the cooling process of quasi-2D liquid carbon in this study. Looking closely at the presented data, the cooling process can be divided into four stages: (I) high-temperature liquid; (II) LLPT (twofold to threefold); (III) LSPT; (IV) optimization stage. The high-temperature liquid is dominated by twofold coordinated structures with very little threefold coordinated structures, while the low-temperature liquid contains mostly threefold coordinated structures with very little twofold structures. In the LLPT stage, the system is transformed from twofold coordinated liquid to threefold coordinated liquid. In this stage, the atomic chains form ring structures mostly induced by threefold coordinated structures. The Y-shaped threefold coordinated structures provide strong nodes to induce the chains to form rings by coalescence. Carbon islands are formed by the combination of rings. The LLPT in liquid carbon can be treated as the preparation stage of the LSPT. Once the size of islands exceeds a critical value, the islands can stable exist as a nucleus and the LSPT starts. The liquid carbon is simultaneously nucleated at individual islands. Eventually, the carbon islands will encounter each other and combine to form a complete polycrystalline film. The surviving time of the threefold coordinated structures is also investigated. Our results show that the surviving time increases with the decrease of temperature, which provides greater opportunities for the formation of rings. The inheritance of the threefold coordinated structures provides essential condition to form rings and islands. These findings provide physical and dynamic insights into the freezing transition in liquid carbon.

Methods

In this work, classical MD simulations are used to study the cooling process of liquid carbon using the MD package LAMMPS28. The semiempirical long-range carbon bond-order potentials (LCBOP models), which have been widely used to describe all carbon phase14,29,30,31,32,33,34,35, are employed to describe the carbon-carbon interaction. The carbon potential used in this study can closely match the ab initio MD results for the liquid structure of carbon and account well for many of the known thermodynamic properties of liquid carbon at high pressures and temperatures29,30,31. The Velocity-Verlet algorithm is used with an integration time step of 1 fs. We build a quasi-2D configuration consisting of 3200 carbon atoms initially arranged into 20 × 20 × 1 diamond lattice with the periodic boundary conditions in the plane directions. According to our modeling scheme, the initial system is firstly heated to 9000 K followed by a relaxation at that temperature for 106 steps to obtain a well-equilibrated liquid carbon. Then the liquid carbon is cooled from 9000 K to 4000 K at a quenching rate of 1 K/ps. All the MD simulations are carried out at a constant plane pressure of 5 GPa controlled via the Nóse-Hoover barostat.