Synthesis of paracrystalline diamond

Solids in nature can be generally classified into crystalline and non-crystalline states1–7, depending on whether long-range lattice periodicity is present in the material. The differentiation of the two states, however, could face fundamental challenges if the degree of long-range order in crystals is significantly reduced. Here we report a paracrystalline state of diamond that is distinct from either crystalline or amorphous diamond8–10. The paracrystalline diamond reported in this work, consisting of sub-nanometre-sized paracrystallites that possess a well-defined crystalline medium-range order up to a few atomic shells4,5,11–13, was synthesized in high-pressure high-temperature conditions (for example, 30 GPa and 1,600 K) employing face-centred cubic C60 as a precursor. The structural characteristics of the paracrystalline diamond were identified through a combination of X-ray diffraction, high-resolution transmission microscopy and advanced molecular dynamics simulation. The formation of paracrystalline diamond is a result of densely distributed nucleation sites developed in compressed C60 as well as pronounced second-nearest-neighbour short-range order in amorphous diamond due to strong sp3 bonding. The discovery of paracrystalline diamond adds an unusual diamond form to the enriched carbon family14–16, which exhibits distinguishing physical properties and can be furthered exploited to develop new materials. Furthermore, this work reveals the missing link in the length scale between amorphous and crystalline states across the structural landscape, having profound implications for recognizing complex structures arising from amorphous materials. A study describes the synthesis, structural characterization and formation mechanism of a paracrystalline state of diamond, adding an unusual form of diamond to the family of carbon-based materials.

Article size (several tens of micrometres) inevitably hampers our understanding of the structure, properties and synthesis mechanism. As such, synthesizing large and high-quality a-D remains both highly demanding and technically challenging. In this study, we attempted to synthesize millimetre-sized non-crystalline diamond using a large-volume multi-anvil press (MAP). Conventional MAP technology could not generate such ultrahigh pressures in a millimetre-sized sample. However, our recently developed ultrahigh-pressure technique using MAP is capable of generating pressures to 50 GPa in millimetre-sized samples 21 . Using this technique, we succeeded in synthesizing p-D with a 1.0 mm diameter at a pressure of 30 GPa. We present details regarding the synthesis and structural analyses of this sample and articulate the synthesis mechanism below.
We selected zero-dimensional fullerene (C 60 ) in a face-centred cubic (fcc) crystal to monitor its phase evolution at 30 GPa and a temperature range of 1,200-1,800 K. Previous investigations have demonstrated that C 60 can transform into polymerized C 60 crystals and disordered sp 2 -sp 3 carbon under HP-HT (13-25 GPa) [22][23][24][25][26] . Here, at an elevated pressure condition of 30 GPa, the X-ray diffraction (XRD) patterns of the recovered samples synthesized at various temperatures are shown in Fig. 1a and Extended Data Fig. 1a. The diffraction signals are distinctly different from the initial fcc C 60 (Fig. 1b), signalling the formation of new phases. Broadened diffraction peaks around ~2.9, ~5.4 and ~8.4 Å −1 are reminiscent of amorphous materials, bearing considerable resemblance to those of a-D (ref. 9 ) and DLC carbon with a high fraction of sp 3 bonds 24 . Moreover, the recovered samples heated above 1,400 K are highly transparent (for example, Fig. 1a inset), suggesting a dense diamond-like (sp 3 -dominated) structure. To quantitatively examine the composition of sp 2 versus sp 3 in the samples, we utilized Raman and electron energy loss spectroscopy (EELS) 27,28 . Raman and EELS characterizations (Fig. 1c, d and Extended Data Fig. 1) indicate that the sample recovered from 1,400 K contains a residual amount (~5.2%) of sp 2 bonds, whereas samples from above 1,400 K are completely sp 3 -bonded, demonstrating successful synthesis of non-crystalline diamond at a pressure much lower than previously reported 8,9 . The densities of the samples at 1,500 and 1,600 K are estimated to be 3.20 and 3.25 g cm −3 (Extended Data Fig. 1), respectively, close to that of a-D (3.30 g cm −3 ) 8,9 . Further increasing the temperature to 1,800 K leads to temperature-induced crystallization (Extended Data Fig. 1a).
High-resolution transmission electron microscopy (HRTEM) images of the recovered samples are shown in Fig. 2 and Extended Data Fig. 2. Lattice fringes associated with LRO are not found in the HRTEM images, indicating the disordered nature of the samples beyond a certain length (~1 nm). We randomly selected small regions (7.0 × 7.0 nm 2 ) and conducted fast Fourier transformation (FFT) to obtain the scattering signals in the momentum space. The diffuse haloes of the FFT patterns corresponding to the selected areas confirm the overall amorphous feature, consistent with the XRD results. On a finer scale, however, visual inspections indicate that the images become more ordered for samples synthesized at elevated temperatures, which is also evidenced by the narrowing FFT rings with increasing temperature. The ordered HRTEM images resemble those of amorphous Si films and bulk metallic glasses with MRO structures 17,29,30 , but they differ from the previously reported a-D described with the CRN model 9 . For the samples annealed at 1,400-1,600 K, a considerable amount of ordered clusters within 0.5-1.0 nm are visible in the HRTEM images. Within the spatial distance of ~1.0 nm, two types of lattice fringes are identifiable, which are found to be close to the atomic arrangements along the [110] and [010] zone axes of CD and HD, respectively, as shown in Fig. 2f. This finding is further validated by the FFT patterns of the selection areas (2.0 × 2.0 nm 2 ), where the intensified 'diffraction' spots matching CD or HD crystalline order are recognizable. Note that both cubic and hexagonal MRO clusters are highly distorted because the inclusive angles (72.1 ± 1.5° for cubic and 89.5 ± 2.5° for hexagonal) of the lattice fringes deviate from those of pristine CD (70.3°) and HD (90°) crystals 31 . The observed morphology of the as-synthesized non-crystalline diamond is clearly distinguished from a-D, hinting at the formation of crystalline MRO structures dominating the samples.
To uncover the hidden structural ordering in the obtained non-crystalline diamond, we relied on large-scale classical molecular dynamics (MD) simulation 32 employing a recently developed realistic angular-dependent potential for carbon 33 and adiabatic-bias MD 34 to accelerate dynamics. Starting from fcc C 60 and mimicking the experimental HP-HT conditions (for example, 1,600 K and 30 GPa), we revealed how the fullerene evolves into the non-crystalline diamond deterministically. On the formation of nearly complete sp 3 bonding at 1,600 K, small but prolific clusters with CD-and HD-like atomic packing gradually take form. These MRO clusters typically encompass and 292 eV are due to transitions of a 1s electron to π* and σ*, corresponding to sp 2 and sp 3 bonding, respectively. The 1,400 K sample has a residual 5.2% of sp 2 bonds, analysed from the ratio of π * and σ * features using single-crystal graphite for calibration 28 .
4-5 atomic shells, henceforth referred to as paracrystallites. With prolonged simulation time, the amount of paracrystallites increases and eventually saturates at a volume fraction of ϕ = 0.7, which is rationalized to be the maximum paracrystal content attainable in the given simulation condition due to thermodynamic constraints. Above a certain threshold (for example, ϕ = 0.3), the paracrystallites start to interconnect and percolate, which is a distinguishing feature of the material. Structural differences of non-crystalline materials are identified by the structure factor S(Q). Subtle but definite differences in S(Q) (for example, a-D versus p-D) can be experimentally discriminated for closely related structures. In our simulation, we accurately determined S(Q) from the atomistic models based on the Ornstein-Zernike equation 35 (see Methods). The comparison of the simulated S(Q) of a-D and p-D (ϕ = 0.7) is provided in Fig. 3a. While both S(Q) profiles are characteristic of non-crystalline materials, the differences are conspicuous, evidenced by peak intensity variations. For a-D, the intensity of the first main diffraction peak at ~2.9 Å −1 is lower than that of the second peak at ~5.4 Å −1 . In comparison to a-D, as the volume fraction of paracrystallites increases, the intensity of the first peak monotonically rises, accompanied by a decrease in the intensity of the second peak (Extended Data Fig. 3a). The intensity ratio of the first peak to the second peak I 1 /I 2 is positively correlated to the degree of MRO arrested in the sample. This relationship is thus used to estimate the content of paracrystallites in the experimental samples. Indeed, the extracted S(Q) of the synthesized samples from the XRD patterns 36 in Fig. 1 exhibits exactly the same predictive trend (Fig. 3b). Plotting I 1 /I 2 against the paracrystal volume fraction ϕ, we find that high proportions of paracrystallites exist in the samples above 1,400 K (Fig. 3c). For instance, the amount of paracrystallites reaches as high as 47% and 52% for samples at 1,500 and 1,600 K, similar to those (46% and 48%) obtained from HRTEM images using autocorrelation function analysis 17 . A satisfactory match between the S(Q) of a paracrystalline model (ϕ = 0.5) and that of the 1,500 K sample is achieved. Figure 3d presents the structural model of a homogeneous p-D with 70% CD-and HD-like paracrystallites (~0.8-1.0 nm in size) and the remaining 30% a-D. The simulated HRTEM image of this atomistic model (Fig. 3e) is in striking agreement with the experimental ones, validating the paracrystalline nature of the samples at 1,500 and 1,600 K.
To further distinguish the structural difference between a-D and p-D and enable identification of CD-or HD-like atomic packing in p-D, we employed the common-neighbour analysis (CNA) 37 and orientational order analysis methods 38 . Figure 4a illustrates the atomic structures of p-D and a-D visualized on the basis of their respective degree of MRO. Through CNA, high fractions of CD-and HD-like clusters are found to be abundant in p-D, but they are starkly absent in a-D. Moreover, the corresponding atomic structure of p-D visualized by the orientational order parameter q 6 (ref. 38 ) reveals the orientational order of the local clusters, supporting the presence of MRO clusters in p-D. To shed light on the size of the paracrystallites, we adopted an orientational correlation function κ(r) (defined in the Methods). Figure 4b shows the radial distribution functions g(r) and the orientational correlation functions of p-D and a-D, respectively. Both a-D and p-D have nearly the same g(r) profile for the first peak, in agreement with their respective tetrahedral atomic environments. However, the intensity of the first peak in the κ(r) of a-D is significantly lower than that of p-D where just two peaks are visible in the profile, suggesting that a-D has SRO within two coordination shells. In contrast, the first five sharp peaks in p-D indicate the strong MRO, and the subsequent three peaks are weakened and broadened because of the absence of LRO. Previously, nano-polycrystalline diamond (NPD) with ultrafine grains has been synthesized from various sp 2 carbons including C 60 under HP-HT conditions 39-43 . Our simulations of NPD with different grain sizes (1.20-2.40 nm) reveal a completely different morphology in NPD due to the presence of LRO (Extended Data Fig. 3). The paracrystallites in p-D are found to be severely distorted in comparison to perfect diamond lattices of which the level of distortion is on par with the grain-boundary atoms of ultrafine NPD (Extended Data Figs. 3c and 4d). This accounts for the absence of Bragg diffraction peaks from its XRD pattern.
We attribute the formation of p-D to two main factors (see Methods). First, we discovered that the CRN-type a-D exhibits strong diamond-like SRO in the first two atomic shells (16 atoms), in comparison to that of amorphous Si (Extended Data Figs. 4 and 5). This is consistent with the argument that sp 3 carbon has the largest tetrahedral parameter among all group IV elements 44 . The pronounced intrinsic diamond-like SRO of a-D greatly facilitates the development of MRO, and consequently

Article
the formation of diamond paracrystallites. Second, the successful synthesis of p-D depends highly on the structure of the C 60 precursor. We argue that the formation of p-D is a result of a nucleation-controlled process, where the formation of orientationally correlated sp 3 bonds between neighbouring C 60 units offers fertile sites for nucleation. Our simulated results (see Methods and Extended Data Figs. 6-10 for details) demonstrate that the formation of p-D goes through three main stages: inter-buckyball bond connection (15% sp 3 bonds) during compression to 30 GPa at room temperature; kinetically controlled polymerization and collapse of buckyballs during heating to 1,600 K at 30 GPa; and the formation of p-D from a-D during isothermal annealing at 1,600 K and 30 GPa. In particular, the sp 3 bonds formed in the early stage of polymerization are orientationally ordered in the local contact regions (Extended Data Fig. 6), which persists into the a-D state at high temperatures and contributes to the formation of p-D. Abundant and homogeneous sp 3 bonding of C 60 by polymerization under high pressure and low temperature prevents the otherwise heterogeneous nucleation and growth of crystals, allowing the formation of p-D in a divide-and-conquer manner. To further understand the critical role of C 60 in the formation of p-D, we examined other sp 2 -dominated carbons in experiments (carbon onion and type-1 glassy carbon) as a precursor in the same HP-HT conditions, which failed to synthesize p-D in our experiments (Extended Data Fig. 1a). Additionally, we surmise that the bundled single-walled carbon nanotubes with well-aligned carbon bonds may provide abundant sp 3 bonding sites through intermolecular polymerization under high pressure 45 , like in C 60 , allowing for the synthesis of p-D.
The unique MRO structures of p-D elicit new physical properties. The mechanical properties of p-D were assessed by the measurements of Vickers and nanoindentation hardness. p-D exhibits highly isotropous Vickers hardness H v of 116.1 ± 3.6 GPa and nanoindentation hardness H N of 102.1 ± 1.4 GPa (Extended Data Fig. 11), comparable with natural diamond 42 . The obtained Young's modulus E (937.2 ± 10.1 GPa), shear modulus G (437.9 ± 4.7 GPa) and bulk modulus K (363.3 ± 3.9 GPa) of p-D are close to those of diamond (1140 GPa, 535 GPa and 443 GPa) 46 , consistent with our ab initio predictions (Extended Data Fig. 11e), making it an appealing candidate for ultrahard materials. Moreover, a measurement of thermostability of p-D in the air (Extended Data Fig. 12) gives an onset oxidation temperature of 950 K, much higher than those (673-873 K) of DLC films 47 , nano diamond 48 , chemical-vapour-deposited diamond 49 and NPD. These outstanding combinations of mechanical properties and oxidation resistance endow p-D with great potential for niche technological applications. Last, we point out that the paracrystalline state diamond may exist in other forms (for example, with different types and compositions of paracrystallites contingent on the staring material and synthesis conditions), waiting to be explored to obtain novel carbon materials with unusual properties.

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41586-021-04122-w.

Sample synthesis
The starting material, C 60 , was purchased from Alfa Aesar. HP-HT experiments were performed using an ultrahigh-pressure MAP 50 at the Bayerisches Geoinstitut of the University of Bayreuth in Germany and the Center for High Pressure Science & Technology Advanced Research (HPSTAR) in China. The synthetic pressure was set to 30 GPa, implemented using carbide anvils with 3-mm truncation and 7-mm MgO + 5 wt% Cr 2 O 3 octahedron. The temperature range was from 1,200 to 1,800 K, generated using a LaCrO 3 heater. The pressure was calibrated at 2,300 K by evaluating the Al 2 O 3 content in bridgmanite 51 , and the temperature was monitored using a type-D thermocouple. The starting material was pre-compressed into Pt capsules with diameters of 0.8 and 1.2 mm and a height of 1.6 mm at a pressure of 200 MPa. The samples were compressed to 30 GPa in 7 h at room temperature and heated to the designated temperatures with a heating rate of 100 K s −1 .
In situ Raman spectroscopy (532 nm) of C 60 under high pressure was conducted without pressure-transmitting medium by using a symmetrical DAC with a culet size of 300 μm. The samples and ruby spheres for pressure calibration were loaded into the DAC sample chamber of 120 μm in diameter with a pre-indented rhenium gasket.

Structure characterization
To collect high-quality XRD data, the powder mode of a Bruker D8 Venture diffractometer equipped with Mo Kα radiation was employed. The samples were fixed on an organic XRD holder and rotated 720° in total within an acquisition time of 10 min. A background file was obtained with an empty organic holder with the same collection condition. To extract the structure factors from the XRD data, we followed the optimization method proposed in ref. 36 . The microscopic structures of the samples were characterized by the FEI Talos-F200s transmission electron microscope operated with an accelerating voltage of 200 keV. The bonding analysis of the samples was performed by using Raman spectroscopy and EELS. The Raman spectra were acquired with two different wavelengths (325 nm and 532 nm) of laser excitation by using a Renishaw-type micro-Raman spectroscopy system. The EELS data were obtained using the ARM-200F ( JEOL) microscope equipped with a spherical aberration corrector (CEOS).

Property measurements
The Vickers hardness H v and fracture toughness K IC were measured using a micro-Vickers hardness tester (Q10A + , Qness) with a Vickers diamond indenter at 4.9 and 9.8 N for a holding time of 15 s. H v and K IC were determined from the following equations 42 : where F (N) is the applied load, L (μm) is the arithmetic mean of the two diagonals of the Vickers indentation, C (μm) is the average length of the radial cracks and E is Young's modulus of p-D. The nanoindentation hardness (H N ) was measured with a three-sided pyramidal Berkovich diamond indenter (Keysight Nano Indenter G200), and the samples were loaded to the maximum load of 50 mN at a loading rate of 2 mN s −1 with a holding time of 10 s. The Young's modulus (E) was derived from the loading/unloading-displacement curves according to the Oliver-Pharr method 52 . The values of Vickers hardness, nanoindentation hardness and Young's modulus (E) are given as an average value taken from more than five tests. The nanoindentation data of p-D were calibrated with natural diamond (100) surface.
The thermostability was evaluated using a thermogravimeter and differential scanning calorimeter, NETZSCHSTA 449 C, from room temperature to 1,773 K in the air with a heating rate of 5 K min −1 .

Structural model of a-D.
In this work, a-D was obtained by conducting extensive ab initio MD simulations, following the protocol introduced in our previous work 9 . To yield better statistics, our simulation system was enlarged to 500 carbon atoms. Ab initio MD was performed with the density-functional-theory-based Vienna Ab-initio Simulation Package 53 . The projector-augmented wave potential 54 with a valence configuration of 2s2p and the generalized gradient approximation were used in the simulation. The kinetic energy cutoff was set to 400 eV, and the simulation was conducted on the gamma point only. In brief, liquid carbon was slowly quenched from 5,000 to 300 K in an NPT (constant number of atoms, pressure and temperature) ensemble under a hydrostatic pressure condition of 50 GPa. The cooling rate was 2.5 × 10 13 K s −1 . A phase transition during the quench process leads to the formation of a-D. The as-obtained high-pressure a-D structure was gradually relaxed to ambient pressure using a conjugate-gradient geometric optimization method. To cross-check the structure resulting from ab initio modelling, a different route was employed to generate a-D. In this approach, we applied the bond-switching algorithm 55 , also known as the WWW method, to generate ideal CRN configurations with different sizes employing the Keating interatomic potential 56 . The 500-atom CRN structure was further subjected to simulated annealing via ab initio MD to minimize local stresses and to achieve a more realistic structural model.

Structural model of p-D and phase transition from C 60 to p-D under HP-HT.
Both classical MD 32 and ab initio MD simulations were carried out to study the structure evolution of C 60 under HP-HT conditions and simulate the formation of p-D. For classical MD, a newly developed angular-dependent interatomic potential (C-ADP) 33,57 was used to describe carbon. The development of the realistic C-ADP potential was based on mapping the complex energy landscape obtained with density-functional-theory-based calculations of more than 5,000 atomic configurations of C (>10 6 atoms in the fitting database) 58,59 .
The empirical potential has been tested to be reliable in describing carbon phase behaviours over a wide temperature and pressure range (0-100 GPa) without compromising the computational speed, enabling us to simulate large samples across an extended temporal scale (in the sub-microsecond regime). The detailed MD simulation steps are as follows. 1) A 36,000-atom fcc-C 60 crystal was created with C 60 buckyballs randomly oriented. The fcc-C 60 crystal was equilibrated at 300 K and zero pressure for 1 ns. 2) The well-equilibrated C 60 crystal was compressed to 30 GPa at room temperature at a constant compression rate of 0.5 GPa ns −1 . 3) The system was brought up to 1,600 K at 30 GPa at a constant heating rate of 1.0 × 10 6 K s −1 . 4) The system was subjected to long-time annealing at 1,600 K and 30 GPa, up to 15 ns. 5) The brute-force MD faces great challenges in extending the timescale to minutes (necessary for the formation of p-D). To circumvent this discrepancy, enhanced sampling approaches must be taken. Here, to bridge the temporal gap between experiment and simulation, we enabled an adiabatic-bias MD (ABMD) technique 34 to accelerate the dynamics, which was tested to be a more effective approach than other barrier-crossing techniques such as metadynamics 60 and the extended Lagrangian method 61 .
In the ABMD simulations stated above, a biasing potential was added to act on a collective variable, driving the system away from the amorphous state and gearing towards the formation of paracrystallites. In this work, the mean value of the order parameter q 6 (that is, , where N refers to the total number of atoms) was used as a collective variable. The parameter q 6 for a particular atom is a coarse-grained quantity that measures the extent to which the orientation of the atoms in the coordination shells matches the orientation of the central atom, revealing the crystalline MRO. The definition of q i 6 is given in the following section. The biasing potential V in the ABMD simulation was constructed as: . η t ( ) is a white noise added to the minimum position of the potential. The parameters used in the bias potential were: K = 5 × 10 5 eV and the target order parameter q = 0.6 t 6 . The ABMD method was implemented in the PLUMED package 62 .
To corroborate the classical MD results of C 60 collapse and polymerization during the early stage of compression, we also conducted ab initio MD simulation to study the phase behaviour of fcc C 60 (720 atoms). To this end, the fcc-C 60 crystal was first slowly compressed to 30 GPa by geometrical optimization. Then, the phase evolution of the compressed C 60 at 30 GPa from 300 to 1,500 K was monitored by conducting ab initio NPT simulation at a heating rate of 2 × 10 14 K s −1 .
On compression to 30 GPa at 300 K, both classical and ab initio MD simulations demonstrate that the C 60 molecules essentially maintain the cage structure but undergo significant deformations (see Extended Data Fig. 6 for details). The compressed C 60 molecules are inter-connected by sp 3 bonding (the fraction of sp 3 bonds is up to ~15% at 30 GPa; Extended Data Fig. 6a). This finding is consistent with our experimental results, where the in situ Raman spectra of C 60 indicate the polymerization and deformation of C 60 molecules under compression (up to ~32 GPa; Extended Data Fig. 6b). With increasing temperature, temperature-induced sp 2 -sp 3 transition (25% to 70%) mainly occurs between 500 and 1,200 K (Extended Data Figs. 7a and 8). There is a percolation transition of sp 3 bonds at 500-800 K, which corresponds to the polymerization of C 60 molecules. Hence, from a materials synthesis point of view, a hybrid sp 2 -sp 3 type of amorphous carbon with a controllable fraction of mixed sp 2 /sp 3 bonding can be obtained by heating to different temperatures (most preferably, above 800 K after the percolation transition takes place). On heating to 1,600 K, the proportion of sp 3 bonds exceeds 93%. Note that the sp 3 fraction was estimated on the basis of the number of four-fold-bonded carbon atoms (that is, CN = 4) within a cutoff distance of 1.85 Å. This value may be slightly underestimated in comparison to that analysed with the maximally localized Wannier function 9 . The structural transition of fcc C 60 under 30 GPa is confirmed by ab initio MD simulation (Extended Data Fig. 7b).
Isothermal annealing at high temperatures (for example, 1,600 K) further promotes the sp 2 -sp 3 transition (Extended Data Fig. 7a). The existence of parallel sp 3 bonds, due to inter-buckyball connections in the early stage of polymerization (identifiable in Extended Data Fig. 6), persists into the a-D annealed at high temperatures. To illustrate this, we calculate the bond-orientation function of the a-D obtained at 1,600 K, defined as , where r is a unit vector with a random orientation, r ij refers to an atomic bond between atom i and its nearest neighbour j, and the angular bracket refers to the ensemble average. If all atomic bonds in the system are randomly distributed, one would expect a uniform distribution of r Φ(ˆ) at all directions. Nonetheless, in compressed C 60 , despite the seemingly amorphous feature, the sp 3 bonds in the a-D have preferred orientations, manifested by the three-dimensional bond-orientation distribution plot and the two-dimensional pole figure (Extended Data Fig. 9). The sp 3 bonds with orientational order in a-D provide fertile sites for nucleation and contribute to the nucleation of p-D during isothermal annealing.
Extending the isothermal annealing time up to 15 ns at 1,600 K and 30 GPa, our classical MD simulation reveals the gradual transition from a-D to p-D via the nucleation and growth mechanism. As shown in Extended Data Fig. 10 and Supplementary Video 1, the fraction of CD-and HD-like paracrystallites gradually increases with prolonged simulation time. Extended Data Fig. 10p and Supplementary Video 1 show the three-dimensional structural model of a p-D with a moderate fraction of paracrystallites (20%) achieved through classical MD. The trajectories of the paracrystallites are illustrated in Extended Data Fig. 10. To further increase the volume fraction of p-D, we resorted to ABMD and the final results are presented in Fig. 3. The appearance of the paracrystallites is dynamically changing, from which one can estimate the critical nucleus size of the paracrystallites (see also Supplementary Video 1).
The critical nuclei are found to contain approximately 20-30 atoms (~3.0 Å), suggesting that the paracrystallites belong to supercritical nuclei. Thermodynamically, these supercritical nuclei have certain probabilities of transitioning back to the amorphous state. Owing to the dynamic fluctuations, the maximum amount of paracrystallites is theoretically limited within a certain fraction (that is, ϕ = 0.7).
Structural model of NPD. Classical MD was used to generate atomic structures of NPD. To construct the initial structure, the Voronoi tessellation method was used to generate nearly equal-sized nanograins with desirable grain sizes (1.2-2.4 nm). The grains were filled with carbon atoms in diamond cubic arrangements. The as-constructed samples, typically consisting of ~100,000 atoms, were thermally annealed at 1,600 K and 30 GPa for 0.2 ns, followed by relaxations to ambient conditions at 300 K and 0 GPa. The exact grain size of the NPD was computed post-mortem, based on: where N is the total number of atoms, ϕ is the fraction of atoms associated with nanocrystals and n g is the number of grains.
Structural analysis. For SRO, we proposed a local order parameter s to gauge the similarity between the first two atomic shells of each atom of a-D and that of CD or HD crystal. Inclusion of the second shell in the consideration of the local atomic environment is justified if one is to perform a Voronoi tessellation 63 to determine atomic coordination. The similarity order parameter s between cluster A and a reference cluster B is defined as: where N A and N B denote the number of atoms in clusters A and B, respectively; r i and r j are the Cartesian coordinates of the ith and jth atoms in the corresponding cluster; σ is a Gaussian smearing parameter, which is set to 0.15 in this work; T is a transformation matrix that accounts for affine scaling transformation and cluster rotation, where l is a scaling parameter and α, β and γ are Euler angles for rotation. Having obtained s CD and s HD , we define s = max(s CD , s HD ). On the basis of the above definition, s lies between 0 and 1 (for example, two identical clusters would yield s = 1 regardless of their respective orientations). A large s value implies a high similarity between the local cluster and the reference crystal (that is, enhanced SRO). The similarity order parameter also provides a facile tool to quantify lattice distortions of the paracrystallites. For MRO, a local Steinhardt order parameter q 6 was used to evaluate the orientational order associated with each atom, which contains ordering information including a few atomic shells. For atom i, the order parameter q i 6 was computed by: where N c is the total number of atoms included in the coordination shells (that is, MRO clusters, typically up to four atomic shells in this work), and q 6m (i) is the Steinhardt vector 64 of atom i, calculated from: where k loops over all atoms in the nearest neighbours of atom i (in this work, CN = 4), and Y 6m is one of the sixth-order spherical harmonics.
While the order parameter q 6 yields information about the degree of structural ordering in the system, it is incapable of revealing the exact type of atomic packing (for example, CD versus HD). To mitigate this insufficiency, the CNA method was adopted 65 . The second-nearest neighbours of each atom in CD have fcc packing with a CNA index of 421, whereas for HD, the 12 second-nearest neighbours are arranged in hexagonal close-packing that can be differentiated with a CNA index of 422 (see Extended Data Fig. 5). By performing CNA, each atom can thus be exclusively labelled as the CD-like, HD-like or disordered type.
Having determined the packing type of the atoms, we used the following algorithm to compute the volume fraction of the paracrystallites in the sample ϕ. First, all atoms identified to have CD-or HD-like packing via CNA analysis, together with their 16 nearest neighbours, are considered as paracrystalline atoms. Second, other four-fold-coordinated atoms that are directly bonded to the aforementioned paracrystalline atoms are included in the calculation of ϕ, to account for the 'boundary' portion of the paracrystallites. Following this algorithm, we computed the volume fraction of the crystalline component of NPD (d = 1.56 nm) to be 94% for benchmarking.
To mathematically quantify the spatial extension of the paracrystallites, we invoked a two-point orientational correlation function, κ(r), given by: where the bracket indicates that the calculation is averaged over all atomic pairs. The κ(r) essentially measures the orientation relationship of two coordination polyhedral at a distance r. In two limiting cases, κ(r) = 1 for perfect crystals, and κ(r) = 0 for totally uncorrelated structures. Empirically, we consider κ(r) = 0.3 as the threshold to estimate the correlation length of the paracrystallites.
Structure factor calculation. Our ab initio modelling provides a reliable and predictive structural model for a-D. Nonetheless, the simulated sample has a finite-size effect. If S(Q) is calculated on the basis of the Fourier transformation of the pair distribution function g(r), truncation ripples will be inevitably incurred in the calculated S(Q). To circumvent this issue, we accurately calculated the structure factors based on the Ornstein-Zernike theory 35 , which relies on solving for the short-ranged direction correlation function c(r) from the Ornstein-Zernike equation: where h(r) = g(r) − 1.
In practice, we followed Baxter's factorization 66 method together with Dixon-Hutchinson's numerical algorithm 67 to obtain c(r) from g(r). Having derived c(r), S(Q) was computed from its Fourier transform C(Q) following: In this way, we were able to accurately analyse the differences between the structure factors of a-D and p-D without being encumbered by artefacts due to the finite-size effect.

HRTEM simulation
HRTEM image simulations were performed on the basis of a multi-slice algorithm as implemented in MULTEM 68 . The sample thickness was typically 80 Å, and a partial coherence illumination mode was used.

Data availability
The data that support the findings of this study are available from the corresponding authors upon request.

Code availability
The software used for data analysis is available from H.S. upon request. in (e, f) indicate cubic-and hexagonal-like MRO clusters, respectively. In subfigure g, the lattice fringes marked with cyan and yellow solid lines indicate that nanosized cubic and hexagonal diamond crystallites are precipitated from the non-crystalline matrix at 1800 K. Fig. 3 | Structure factor S(Q), structural models, and simulated HRTEM images of simulated a-D, p-D with different fractions of paracrystallites, and NPD with different grain sizes. a, S(Q) of simulated a-D and p-D with different fractions of paracrystallites. With increasing ϕ, the intensity of the first peak increases, and the position of the peak shifts to high-Q, which is consistent with the experimental observation that at higher annealing temperatures, the intensity of the first peak increases and the position moves to the right (from 2.89Å −1 at 1200 K to 2.95 Å −1 at 1600 K, Fig. 3b in the main text). b, S(Q) of simulated NPDs with different grain sizes, juxtaposed with that of the 1800 K sample from our experiment. For the simulated NPD samples, the peak intensity decreases and the peak width broadens as the grain size is refined. Bragg peaks are clearly discernable for the samples with a grain size of 1.2 nm. Experimentally, the protruding diffraction peaks in the 1800 K sample corresponds exactly to the characteristic Bragg peaks of NPD, indicative of crystallization. c, In ultrafine NPD, the less-distorted "core" regions of the nanocrystals (blue spheres) are clearly seen. For NPD with d = 1.56 nm, ϕ is estimated to be 94%, meaning 6% atoms are located at the interfaces and cannot be properly assigned to any of the crystalline grains. Fig. 4 | Short-range ordering of a-D, a-Si, p-D, and NPD. a, Distributions of tetrahedrality q t for a-D and a-Si, both of which were obtained from the WWW bond-switching approach and further relaxed with ab initio MD. The tetrahedral order parameter q t is defined as ( ) , where θ ij is the angle formed by two vectors pointing from the center silicon to two of its neighboring silicon atoms. The summation runs over all the combinations of the four nearest neighbors. The narrow distribution of a-D clearly indicates that a-D has a strong tetrahedral SRO than a-Si. b, Bond-angle distribution within the first atomic shell for a-D and a-Si. Compared with a-Si, a-D has a narrower bond-angle distribution. c, Distributions of the second-nearest-neighbour order parameter s of a-D and a-Si. The order parameter s measures the similarity of the local atomic environment (16 atoms) of a central atom with respect to perfect CD or HD lattice arrangements. The higher the s value, the less distortion of the local cluster. Compared with a-Si, the distribution function of s for a-D shifts to high s values with a centroid of 0.42, suggesting much enhanced SRO within the first two atomic shells in a-D. The high similarity of the local atomic environments between a-D and crystalline diamond stems from the strong directional sp 3 bonding of carbon, as evidenced from the narrower bond-angle distribution of a-D in comparison to that of a-Si. d, Distributions s of a-D, p-D, and NPD with a grain size of 1.92 nm. The inset atomic clusters illustrate the cluster configurations (red spheres) corresponding to their degree of similarity with the perfect CD lattice arrangement (green spheres). Paracrystalline diamond has a medium s value (0.51) between that of a-D and NPD, revealing crystalline ordering in p-D. It is worth noting that, for NPD with d = 1.92 nm, the distribution spike at s = 0.96 corresponds to interior atoms in nano-grains that are less distorted. Such atoms are absent in paracrystalline diamond, marking the distinction between p-D and NPD.