# Addressing the effect of stacking faults in X-ray diffractograms of graphite through atom-scale simulations

Pascal Puech, Mathilde Jeanningros, David Neumeyer, Marc Monthoux

► To cite this version:

Pascal Puech, Mathilde Jeanningros, David Neumeyer, Marc Monthoux. Addressing the effect of stacking faults in X-ray diffractograms of graphite through atom-scale simulations. Carbon Trends, 2023, 13, pp.100311. 10.1016/j.cartre.2023.100311. hal-04553433

HAL Id: hal-04553433

https://hal.science/hal-04553433

Submitted on 20 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_1.png" /></figure>

# Addressing the effect of stacking faults in X-ray diffractograms of graphite through atom-scale simulations

Pascal Puech *, Mathilde Jeanningros, David Neumeyer, Marc Monthoux *

CEMES, UPR-8011 CNRS, Université de Toulouse, 29 rue Jeanne Marvig, 31055 Toulouse Cedex 4, France

<figure><img src="image_2.png" /></figure>

Check for
updates

## ARTICLE INFO

### Keywords:

X-ray diffraction

Graphite

Structure

Stacking faults, Modeling

## ABSTRACT

Determining the average contribution and actual distribution of rhombohedral stacking sequences within the average graphitic crystallite remained unsolved. To address this issue, we used a bottom-up approach to simulate X-ray diffractograms of graphite crystallites up to one million atoms in various graphene stacking configurations to find out how to match experimental data. The 100–101 2θ range was highly sensitive to the presence of C layers and can be used to accurately characterize the stacking disorder. On the opposite, the 110–112 2θ range was not sensitive to stacking disorder and thus can be taken as a reference to obtain the average in-plane size L<sub>a</sub> of the crystallites, faulted or not. We used the principle of the L'c parameter, introduced in a former work, which corresponds to the height of the coherent sub-domains within the average crystallite obtainable through the shape of the hkl peaks involving both in-plane and out-of-plane directions. Three types of defects within a Bernal sequence were analyzed, which do not affect the peaks the same way: (i) a substitution of A or B by C, (ii) an insertion of C leading to destructive interferences, and (iii) a so-called "block shifting". It is thus possible to go further in the interpretation of X-Ray diffraction. Finally, we discussed how the overall crystallite size L<sub>c</sub> can be obtained from all the 00l peaks.

## 1. Introduction

X-Ray diffraction (XRD) is a powerful way to study graphene-based carbon material structure. Graphene is a monolayer of sp<sup>2</sup> hybridized carbon atoms in honeycomb lattice [1]. The carbon material structure is highly anisotropic and layered. Therefore, how the graphene layers stack is a major issue. Random stacking (turbostratic structure) leads to asymmetric, two-indices *hk* peaks associated with the in-plane crystallite dimensions [2]. This turbostratic structure is typical of most chars, i. e., resulting from the carbonization of organic precursors below ~1600 °C. During carbonization at higher temperatures, AB pairs start to form which can be evidenced by XRD [3]. Beyond 2000 °C, genuine Bernal (hexagonal) structure can form, leading for a while to a quite complex situation where turbostratic stacking, AB pairs, and Bernal stacking may co-exist as so-called Basic Structural Components (BSCs) within the same average crystallite [3]. Afterwards, for graphitizable precursors heat-treated at graphitization temperature (> 2600 °C), the stacking situation in crystallites becomes simpler and tends towards the pure hexagonal structure. We use here the term "graphite" when the graphene layers in a crystallite are stacked according to the

tri-dimensional order with a prevalence of the AB hexagonal sequence even if some rhombohedral stacking sequences are present. This is identical to the IUPAC terminology where graphite corresponds to graphene layers stacked parallel to each other in a three-dimensional crystalline long-range order [4]. In principle, then, the rhomboedral structure (the cell of which is described by an ABC stacking sequence) is not obtained because the hexagonal structure (the cell of which is described by an AB stacking sequence) thermodynamically prevails. However, the equilibrium temperature between both phases is around 450 °C [5]. This can be evidenced on XRD diffractograms [6] and optical spectra [7]. Therefore, for natural graphite or graphitizable carbons graphitized under non-standard conditions, the 100 2θ range remains complex to interpret. In natural graphites indeed, the rhombohedral stacking can coexist with the Bernal one, because graphitization was reached for extended durations (millions of years) at much lower temperatures than in laboratory processes, and/or because geological conditions have induced shearing stresses. The effect of the latter is assessed by grinding graphite which makes the XRD peaks in the 100–101 2θ range become a broad band with complex structuration [8].

Recently, it was proposed to fit experimental X-ray diffractograms

<footnote>* Corresponding authors.
E-mail addresses: pascal.puech@cemes.fr (P. Puech), marc.monthoux@cemes.fr (M. Monthoux).</footnote>
<figure><img src="image_2.png" /></figure>

# Addressing the effect of stacking faults in X-ray diffractograms of graphite through atom-scale simulations

Pascal Puech *, Mathilde Jeanningros, David Neumeyer, Marc Monthoux *

CEMES, UPR-8011 CNRS, Université de Toulouse, 29 rue Jeanne Marvig, 31055 Toulouse Cedex 4, France

<figure><img src="image_3.png" /></figure>

Check for
updates

## ARTICLE INFO

### Keywords:

X-ray diffraction

Graphite

Structure

Stacking faults, Modeling

## ABSTRACT

Determining the average contribution and actual distribution of rhombohedral stacking sequences within the average graphitic crystallite remained unsolved. To address this issue, we used a bottom-up approach to simulate X-ray diffractograms of graphite crystallites up to one million atoms in various graphene stacking configurations to find out how to match experimental data. The 100–101 2θ range was highly sensitive to the presence of C layers and can be used to accurately characterize the stacking disorder. On the opposite, the 110–112 2θ range was not sensitive to stacking disorder and thus can be taken as a reference to obtain the average in-plane size L<sub>a</sub> of the crystallites, faulted or not. We used the principle of the L'c parameter, introduced in a former work, which corresponds to the height of the coherent sub-domains within the average crystallite obtainable through the shape of the hkl peaks involving both in-plane and out-of-plane directions. Three types of defects within a Bernal sequence were analyzed, which do not affect the peaks the same way: (i) a substitution of A or B by C, (ii) an insertion of C leading to destructive interferences, and (iii) a so-called "block shifting". It is thus possible to go further in the interpretation of X-Ray diffraction. Finally, we discussed how the overall crystallite size L<sub>c</sub> can be obtained from all the 00l peaks.

## 1. Introduction

X-Ray diffraction (XRD) is a powerful way to study graphene-based carbon material structure. Graphene is a monolayer of sp<sup>2</sup> hybridized carbon atoms in honeycomb lattice [1]. The carbon material structure is highly anisotropic and layered. Therefore, how the graphene layers stack is a major issue. Random stacking (turbostratic structure) leads to asymmetric, two-indices *hk* peaks associated with the in-plane crystallite dimensions [2]. This turbostratic structure is typical of most chars, i. e., resulting from the carbonization of organic precursors below ~1600 °C. During carbonization at higher temperatures, AB pairs start to form which can be evidenced by XRD [3]. Beyond 2000 °C, genuine Bernal (hexagonal) structure can form, leading for a while to a quite complex situation where turbostratic stacking, AB pairs, and Bernal stacking may co-exist as so-called Basic Structural Components (BSCs) within the same average crystallite [3]. Afterwards, for graphitizable precursors heat-treated at graphitization temperature (> 2600 °C), the stacking situation in crystallites becomes simpler and tends towards the pure hexagonal structure. We use here the term "graphite" when the graphene layers in a crystallite are stacked according to the

tri-dimensional order with a prevalence of the AB hexagonal sequence even if some rhombohedral stacking sequences are present. This is identical to the IUPAC terminology where graphite corresponds to graphene layers stacked parallel to each other in a three-dimensional crystalline long-range order [4]. In principle, then, the rhomboedral structure (the cell of which is described by an ABC stacking sequence) is not obtained because the hexagonal structure (the cell of which is described by an AB stacking sequence) thermodynamically prevails. However, the equilibrium temperature between both phases is around 450 °C [5]. This can be evidenced on XRD diffractograms [6] and optical spectra [7]. Therefore, for natural graphite or graphitizable carbons graphitized under non-standard conditions, the 100 2θ range remains complex to interpret. In natural graphites indeed, the rhombohedral stacking can coexist with the Bernal one, because graphitization was reached for extended durations (millions of years) at much lower temperatures than in laboratory processes, and/or because geological conditions have induced shearing stresses. The effect of the latter is assessed by grinding graphite which makes the XRD peaks in the 100–101 2θ range become a broad band with complex structuration [8].

Recently, it was proposed to fit experimental X-ray diffractograms

<footnote>* Corresponding authors.
E-mail addresses: pascal.puech@cemes.fr (P. Puech), marc.monthoux@cemes.fr (M. Monthoux).</footnote>
following a kind of bottom-up approach. From Fujimoto (who investigated the apparent primitive cell size) [9] to Puech et al. (who revealed the occurrence and signature of AB pairs) [3] and Putman et al. (who investigated the effect of density) [10], simulating diffractograms is becoming a powerful tool to analyze the experimental data. Thanks to the growing availability of Graphics Processing Units computational resources, it is possible to simulate many diffractograms, and to account for crystalline structures involving millions of atoms [3]. This allows going deeper in the understanding of these materials.

In this paper, we investigate the possibility of fully calculating the experimental X-ray diffractogram of a commercial graphite material by using various approaches using the BSC concept introduced earlier [3]. We show how to proceed to obtain the crystallite size, the size distribution of the coherent domains within the average crystallite, and the kind of defects that are required to explain experimental XRD diffractograms obtained from graphite materials.

## 2. Experimental

### 2.1. Computation of the diffractograms from the models

Although examples of the occurrence of extended rhombohedral sequences within an overall crystallite in which hexagonal stacking prevails can be found in the literature [6,11], we considered that this must be a particular case with respect to the standard behavior of graphitizable carbon precursors when subjected to standard graphitization processes. In laboratory conditions, first coherent stacking occur far beyond the equilibrium temperature between the rhombohedral and hexagonal structures. Therefore, the formation of fully developed (i.e., over more than 4 layers) rhombohedral sequences is unfavored. As a consequence, the latter may mostly exist because of the local impossibility, for a layer located between two Bernal sequences being developed within a crystallite, to conform coherently with both so that they could merge as a single Bernal sequence. Therefore, C layers are more likely to be found once in a while as stacking faults within a crystallite than concentrated in full rhombohedral sequences. Correspondingly, the various stacking faults involving layers in C stacking position (as opposed to A and B stacking positions) are the only defects considered in this work, in addition to some random (turbostratic) stacking. The effects of curvature, line defects, or dislocations are not raised but, in graphitized materials, they are believed to be subordinate with respect to the effect of stacking faults. The configurations were obtained by first constructing crystallite models with MATLAB®, and then by calculating the related X-ray diffractograms using the highly parallelized XaNsoNS freeware [12] operated on a mere game computer. For this work, structures typically involving half a million atoms dispatched as graphene disks stacked according to 500 different configurations have been considered. For the calculation, $L_a$ is the diameter of the graphene disk unit in the average crystallite, in which the distance between two in-plane bonded carbon atoms is 0.142 nm. For coherent stacking (involving A, B, and possibly C layers), an intergraphene distance of 0.335 nm is selected. For incoherent stacking (as it happens between stacked entities, see below), a random rotation is introduced (i.e. turbostratic configuration), and the inter-graphene distance is increased to 0.344 nm. As the intensity calculated by XaNsoNS is per carbon atom, the area of peaks corresponding to atom plane families parallel or perpendicular to the z-axis (i.e. *00l* and *hk0*) varies slightly when the dimensions of the crystallite are changed. Those variations are partially due to the volumetric hypotheses on the atom stack considered for the calculation and are quite subordinate (see Fig. S1 in Supplementary Information (SI) for instance).

### 2.2. Experimental X-ray diffraction pattern

XRD diffraction patterns were recorded from a commercial natural graphite (Fisher) and petroleum-pitch-based cokes annealed at 2100 and

2500 °C on a Bruker RX-D8 Advance X-ray diffractometer with a non-monochromated Cu Kα radiation source (average λ = 1.5406 Å). The patterns were corrected from the apparatus function (including the contribution of $K_{α2}$ and $K_{β}$) using the Diffra.EVA software also from Bruker. The cokes were prepared from a coal tar pitch (reference GFEC-1) by Le Carbone Lorraine (now MERSEN) and obtained from the French Carbon Society.

### 2.3. Methodology

Summing up our previous work and the current ones, a variety of models were considered for calculating diffractograms, as sketched in Fig. 1.

Our previous work [3,13] has shown that a "single-crystallite" approach was powerful to account for X-ray diffractograms of graphitizable cokes annealed from 600 to 2650 °C. More in the detail, this approach consisted of considering each diffractogram to result from a single average crystallite composed of *n* entities, each of them exhibiting one of the 3 BSC types considered in that case (Turbostratic, AB pair, Bernal). Entities are coherent sub-parts – when any - of a crystallite that are rotated randomly to each other. One entity corresponds to a single BSC. To build model *D* in Fig. 1, entities with the Bernal BSC were all given the same single average height $L'_c$. As reminded in the Introduction section, starting from an average crystallite composed of an ensemble of individual graphenes turbostratically stacked, this approach has revealed two clear transitions, (i) around 1600 °C from which AB pairs start to form, and (ii) around 2000 °C from which the Bernal structure starts to develop at the expense of the other BSCs.

In summary, a crystallite is modelled by a cylinder. Its diameter corresponds to $L_a$ and its height to $L_c$. However, several entities of height $L'_c$ can coexist in the crystallite, which can be similar, or not. For example, if a crystallite is only made of individual layers turbostratically stacked, each layer is an entity and the height parameter $L'_c$ is that of one layer; if a crystallite is only made of AB pairs, $L'_c$ is that of 2 layers. For a crystallite involving the Bernal BSC only, i.e., with long-range stacking order, two cases are possible, for which $L_a$ and $L_c$ are however the same, as quantified via the *hk0* and *00l* peaks, respectively. One case is when the whole crystallite is ideally a single entity, hence $L'_c = L_c$ in that case. The other case is when the crystallite is made of several entities, each corresponding to a Bernal stacking, but randomly rotated with respect to each other. For instance, if a random rotation is introduced in the middle of the ideal crystallite, the entity height $L'_c = L_c/2$, whereas $L_a$ and $L_c$ remain unchanged. $L'_c$ only affects peaks involving both in-plane and out-of-plane directions, such as 101 or 112. Therefore, we see that it is necessary to introduce the $L'_c$ parameter to discriminate between both cases, the existence of which explains that the linewidth of peaks such as 101 or 112 is never consistent with the $L_a$ and $L_c$ crystallite sizes obtained from *hk0* and *00l* peaks.

A similar approach is adopted here to analyze experimental diffraction patterns, taking the example of a commercial graphite material. We first tracked, in the diffraction patterns, the regions where specific features can be exploited, i.e. diffractogram features that do not match that of a pure Bernal (hexagonal) BSC, and then we assumed that those features account for the possible contributions of other BSCs, typically the rhombohedral structure, or for the existence of layers in defective stacking (i.e. in C position, by opposition to A and B positions) for the least (it is reminded that a layer B is shifted by a length of a carbon-carbon bond in one of the three directions of such a bond with respect to a layer A, and that, likewise, a layer C is shifted the same way with respect to a layer B. As a result, A and C do not superimpose). Based on these principles, a large variety of stacking configurations were tested (Fig. 1). And then, when necessary, we introduced stacking faults and analyzed their major consequences on the diffractograms to determine which parameters are pertinent to describe X-ray diffractograms of graphitic carbons. In this process, for systems containing a very large number of atoms, we show that carrying-out the calculations by
following a kind of bottom-up approach. From Fujimoto (who investigated the apparent primitive cell size) [9] to Puech et al. (who revealed the occurrence and signature of AB pairs) [3] and Putman et al. (who investigated the effect of density) [10], simulating diffractograms is becoming a powerful tool to analyze the experimental data. Thanks to the growing availability of Graphics Processing Units computational resources, it is possible to simulate many diffractograms, and to account for crystalline structures involving millions of atoms [3]. This allows going deeper in the understanding of these materials.

In this paper, we investigate the possibility of fully calculating the experimental X-ray diffractogram of a commercial graphite material by using various approaches using the BSC concept introduced earlier [3]. We show how to proceed to obtain the crystallite size, the size distribution of the coherent domains within the average crystallite, and the kind of defects that are required to explain experimental XRD diffractograms obtained from graphite materials.

## 2. Experimental

### 2.1. Computation of the diffractograms from the models

Although examples of the occurrence of extended rhombohedral sequences within an overall crystallite in which hexagonal stacking prevails can be found in the literature [6,11], we considered that this must be a particular case with respect to the standard behavior of graphitizable carbon precursors when subjected to standard graphitization processes. In laboratory conditions, first coherent stacking occur far beyond the equilibrium temperature between the rhombohedral and hexagonal structures. Therefore, the formation of fully developed (i.e., over more than 4 layers) rhombohedral sequences is unfavored. As a consequence, the latter may mostly exist because of the local impossibility, for a layer located between two Bernal sequences being developed within a crystallite, to conform coherently with both so that they could merge as a single Bernal sequence. Therefore, C layers are more likely to be found once in a while as stacking faults within a crystallite than concentrated in full rhombohedral sequences. Correspondingly, the various stacking faults involving layers in C stacking position (as opposed to A and B stacking positions) are the only defects considered in this work, in addition to some random (turbostratic) stacking. The effects of curvature, line defects, or dislocations are not raised but, in graphitized materials, they are believed to be subordinate with respect to the effect of stacking faults. The configurations were obtained by first constructing crystallite models with MATLAB®, and then by calculating the related X-ray diffractograms using the highly parallelized XaNsoNS freeware [12] operated on a mere game computer. For this work, structures typically involving half a million atoms dispatched as graphene disks stacked according to 500 different configurations have been considered. For the calculation, $L_a$ is the diameter of the graphene disk unit in the average crystallite, in which the distance between two in-plane bonded carbon atoms is 0.142 nm. For coherent stacking (involving A, B, and possibly C layers), an intergraphene distance of 0.335 nm is selected. For incoherent stacking (as it happens between stacked entities, see below), a random rotation is introduced (i.e. turbostratic configuration), and the inter-graphene distance is increased to 0.344 nm. As the intensity calculated by XaNsoNS is per carbon atom, the area of peaks corresponding to atom plane families parallel or perpendicular to the z-axis (i.e. *00l* and *hk0*) varies slightly when the dimensions of the crystallite are changed. Those variations are partially due to the volumetric hypotheses on the atom stack considered for the calculation and are quite subordinate (see Fig. S1 in Supplementary Information (SI) for instance).

### 2.2. Experimental X-ray diffraction pattern

XRD diffraction patterns were recorded from a commercial natural graphite (Fisher) and petroleum-pitch-based cokes annealed at 2100 and

2500 °C on a Bruker RX-D8 Advance X-ray diffractometer with a non-monochromated Cu Kα radiation source (average λ = 1.5406 Å). The patterns were corrected from the apparatus function (including the contribution of $K_{\alpha 2}$ and $K_{\beta}$) using the Diffra.EVA software also from Bruker. The cokes were prepared from a coal tar pitch (reference GFEC-1) by Le Carbone Lorraine (now MERSEN) and obtained from the French Carbon Society.

### 2.3. Methodology

Summing up our previous work and the current ones, a variety of models were considered for calculating diffractograms, as sketched in Fig. 1.

Our previous work [3,13] has shown that a "single-crystallite" approach was powerful to account for X-ray diffractograms of graphitizable cokes annealed from 600 to 2650 °C. More in the detail, this approach consisted of considering each diffractogram to result from a single average crystallite composed of *n* entities, each of them exhibiting one of the 3 BSC types considered in that case (Turbostratic, AB pair, Bernal). Entities are coherent sub-parts – when any - of a crystallite that are rotated randomly to each other. One entity corresponds to a single BSC. To build model *D* in Fig. 1, entities with the Bernal BSC were all given the same single average height $L'_c$. As reminded in the Introduction section, starting from an average crystallite composed of an ensemble of individual graphenes turbostratically stacked, this approach has revealed two clear transitions, (i) around 1600 °C from which AB pairs start to form, and (ii) around 2000 °C from which the Bernal structure starts to develop at the expense of the other BSCs.

In summary, a crystallite is modelled by a cylinder. Its diameter corresponds to $L_a$ and its height to $L_c$. However, several entities of height $L'_c$ can coexist in the crystallite, which can be similar, or not. For example, if a crystallite is only made of individual layers turbostratically stacked, each layer is an entity and the height parameter $L'_c$ is that of one layer; if a crystallite is only made of AB pairs, $L'_c$ is that of 2 layers. For a crystallite involving the Bernal BSC only, i.e., with long-range stacking order, two cases are possible, for which $L_a$ and $L_c$ are however the same, as quantified via the *hk0* and *00l* peaks, respectively. One case is when the whole crystallite is ideally a single entity, hence $L'_c = L_c$ in that case. The other case is when the crystallite is made of several entities, each corresponding to a Bernal stacking, but randomly rotated with respect to each other. For instance, if a random rotation is introduced in the middle of the ideal crystallite, the entity height $L'_c = L_c/2$, whereas $L_a$ and $L_c$ remain unchanged. $L'_c$ only affects peaks involving both in-plane and out-of-plane directions, such as 101 or 112. Therefore, we see that it is necessary to introduce the $L'_c$ parameter to discriminate between both cases, the existence of which explains that the linewidth of peaks such as 101 or 112 is never consistent with the $L_a$ and $L_c$ crystallite sizes obtained from *hk0* and *00l* peaks.

A similar approach is adopted here to analyze experimental diffraction patterns, taking the example of a commercial graphite material. We first tracked, in the diffraction patterns, the regions where specific features can be exploited, i.e. diffractogram features that do not match that of a pure Bernal (hexagonal) BSC, and then we assumed that those features account for the possible contributions of other BSCs, typically the rhombohedral structure, or for the existence of layers in defective stacking (i.e. in C position, by opposition to A and B positions) for the least (it is reminded that a layer B is shifted by a length of a carbon-carbon bond in one of the three directions of such a bond with respect to a layer A, and that, likewise, a layer C is shifted the same way with respect to a layer B. As a result, A and C do not superimpose). Based on these principles, a large variety of stacking configurations were tested (Fig. 1). And then, when necessary, we introduced stacking faults and analyzed their major consequences on the diffractograms to determine which parameters are pertinent to describe X-ray diffractograms of graphitic carbons. In this process, for systems containing a very large number of atoms, we show that carrying-out the calculations by
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Illustration of the various configurations used in this paper (all models but D) and in our previous work (models A, B and D) [3,13] to model the crystallite(s) from which calculated diffractograms were obtained. Layers with black color stand for Bernal BSC; other colors but red stand for other BSCs (Turbostratic and AB pairs); red layers stand for stacking faults (C layers). Entities with different line drawings indicate different rotational orientations, hence separated by the turbostratic distance (0.344 nm).</figcaption></figure>

considering another approach, the "multiple-crystallite" approach, could
be a way to overcome the difficulty inherent to considering single
crystallites (models G to I in Fig. 1).

2.4. Fitting

For fitting each *hkl* peak of the experimental diffractogram, a pseudo-Voigt function is used, corresponding to a weighted sum (determined by fitting) of a Lorentzian and a Gaussian having the same Full Width at
<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Illustration of the various configurations used in this paper (all models but D) and in our previous work (models A, B and D) [3,13] to model the crystallite(s) from which calculated diffractograms were obtained. Layers with black color stand for Bernal BSC; other colors but red stand for other BSCs (Turbostratic and AB pairs); red layers stand for stacking faults (C layers). Entities with different line drawings indicate different rotational orientations, hence separated by the turbostratic distance (0.344 nm).</figcaption></figure>

considering another approach, the "multiple-crystallite" approach, could be a way to overcome the difficulty inherent to considering single crystallites (models G to I in Fig. 1).

## 2.4. Fitting

For fitting each *hkl* peak of the experimental diffractogram, a pseudo-Voigt function is used, corresponding to a weighted sum (determined by fitting) of a Lorentzian and a Gaussian having the same Full Width at
Half Maximum (FWHM). FWHM<sub>hkl</sub> allows calculating the extension of the (hkl) atomic planes in the average crystallite using the Scherrer law [14]:

$$L_{hkl} = \frac{K\lambda}{FWHM_{hkl} \cos(\theta)}$$

where K = 0.89 is used for both the Bernal and rhombohedral structures [15]. It is worth noting that, in our previous work dealing with non-fully graphitized cokes [3], we generated ad hoc functions to get rid of the use of the Scherrer law for two-indice peaks *hk* typical of the turbostratic structure and the related need to determine *K* value, which has remained an issue for decades.

# 3. Results

## 3.1. The single-entity approach

The single-entity approach was first used in our previous work [3]. It considers that the experimental diffractogram accounts for a single average crystallite composed of a single entity with a given BSC (model A in Fig. 1). Of course, this is a textbook case, with little chance to occur, but useful to evidence the basic differences between the diffractogram features of the various BSCs. Therefore, another single-entity model is tested, using a multiple-crystallite approach which consists of calculating diffractograms corresponding to the sum of *n* crystallites, each containing a single entity which the $L_a$ are distributed according to a Gaussian and with a same $L_c$ value (model G in Fig. 1).

The result for both models is reported and compared in Fig. 2, considering the example of $L_a$ = 10 nm for both the single average-

<figure><img src="image_4.png" /><figcaption>Fig. 2. Calculated diffractograms for the Turbostratic, AB pair, Bernal, and rhombohedral BSCs respectively, considering either a single average crystallite with $L_a$ = $L_c$ = 10 nm (blue color, corresponding to model A in Fig. 1), or the Gaussian weighted sum of diffractograms for $L_a$ ranging from 5 to 15 nm (by step of 1 nm) and $L_c$ = 10 nm (red color, corresponding to model G in Fig. 1). The fact that the *00l* peaks for the AB pair BSC are in intermediate 2θ positions between that of the Bernal BSC and that of the Turbostratic BSC is due to the contribution of the turbostratic configuration between the pairs piled up within the 10 nm-high crystallite.</figcaption></figure>
Half Maximum (FWHM). FWHM<sub>hkl</sub> allows calculating the extension of the (hkl) atomic planes in the average crystallite using the Scherrer law [14]:

$$L_{hkl} = \frac{K\lambda}{FWHM_{hkl} \cos(\theta)}$$

where K = 0.89 is used for both the Bernal and rhombohedral structures [15]. It is worth noting that, in our previous work dealing with non-fully graphitized cokes [3], we generated ad hoc functions to get rid of the use of the Scherrer law for two-indice peaks *hk* typical of the turbostratic structure and the related need to determine *K* value, which has remained an issue for decades.

### 3. Results

#### 3.1. The single-entity approach

The single-entity approach was first used in our previous work [3]. It considers that the experimental diffractogram accounts for a single average crystallite composed of a single entity with a given BSC (model A in Fig. 1). Of course, this is a textbook case, with little chance to occur, but useful to evidence the basic differences between the diffractogram features of the various BSCs. Therefore, another single-entity model is tested, using a multiple-crystallite approach which consists of calculating diffractograms corresponding to the sum of *n* crystallites, each containing a single entity which the $L_a$ are distributed according to a Gaussian and with a same $L_c$ value (model G in Fig. 1).

The result for both models is reported and compared in Fig. 2, considering the example of $L_a$ = 10 nm for both the single average-

<figure><img src="image_5.png" /><figcaption>Fig. 2. Calculated diffractograms for the Turbostratic, AB pair, Bernal, and rhombohedral BSCs respectively, considering either a single average crystallite with $L_a$ = $L_c$ = 10 nm (blue color, corresponding to model A in Fig. 1), or the Gaussian weighted sum of diffractograms for $L_a$ ranging from 5 to 15 nm (by step of 1 nm) and $L_c$ = 10 nm (red color, corresponding to model G in Fig. 1). The fact that the *00l* peaks for the AB pair BSC are in intermediate 2θ positions between that of the Bernal BSC and that of the Turbostratic BSC is due to the contribution of the turbostratic configuration between the pairs piled up within the 10 nm-high crystallite.</figcaption></figure>
crystallite model (blue), and the median value of the Gaussian $L_a$ dis- tribution for the multiple crystallite model (red). For all, $L_c$ = 10 nm. The blue curve (single-crystallite approach) is perfectly masked by the red curve (multiple-crystallite approach), meaning that both approaches are equivalent, at least when the size distribution (either through $L_a$, as in Fig. 2, or through $L_c$) is Gaussian. In the case of the AB pair BSC, one can note that, if only the 10 peak region is used, confusion with the Tur- bostratic BSC can be made, hence considering both the 10 and 11 peak regions is needed to discriminate between both BSCs.

What is true for the Gaussian-shaped peak 100 is also true for 002: a Gaussian size distribution gives a diffractogram comparable to that of the average value. As each diffractogram is normalized by the number of carbon atoms, a Gaussian distribution means more small crystallites than large crystallites.

One remark can be made about the peak indexation, which uses three Laue indices: since the beginning of the 20th century, because of the prevalence of the hexagonal structure, the peak related to the inter- graphene spacing in graphenic carbons has always been indexed 002 in all publications to make simple the notation when dealing with car- bon with mixed structures. However, this is inconsistent with the characteristics of the unit cell respective to each structure. Therefore, with the rhombohedral structure, the interlayer spacing is associated with the 003 peak (because there are three graphene layers in the unit cell) while, with the Bernal structure, it is written 002 (because there are two graphene layers in the unit cell). Likewise, in this work, it is written 001 for the turbostratic structure as the unit cell is made of a single

graphene layer only.

## 3.2. The single crystallite/multi-entities approach applied to a Bernal crystallite

Comparing the Bernal and rhombohedral graphite structures (Fig. 2), the 103 peak seems very interesting as this peak is specific to Bernal stacking, and likewise for the 102 peak in the Bernal diffractogram, although the peak intensity is less. The Bernal 112 peak is also inter- esting as it is located at the same position and with the same intensity as the rhombohedral (in that case, the indexation is 113 but corresponds to the same plane family), hence, it could be taken as a reference if needed.

In Fig. 3 are reported many diffractograms corresponding to a single crystallite model for two examples of $L_a$ values (10 and 20 nm, red and blue curves, respectively) and a single $L_c$ value (~12 nm, corresponding to 36 stacked layers). Therefore, for a given $L_a$, the diffractograms result from the same number of atoms. For each $L_a$ value, the average crys- tallite is made of *n* entities, all with the same BSC (Bernal) and same height $L'_c$, while the entities are turbostratically stacked to each other within the average crystallite (model B in Fig. 1). The various dif- fractograms for a given $L_a$ are obtained by varying $L'_c$ from $L'_c$ = $L_c$ down to $L'_c$ = 1 nm corresponding to a minimum of three stacked layers (two layers would correspond to the AB pair BSC type). In the SI is reported the diffractogram for the ABA 3-layer BSC (Fig. S2).

One can note that the 112 peak is clearly visible even with two layers (hence corresponding to the AB pair BSC). For both peaks of interest

<figure><img src="image_5.png" /><figcaption>**Fig. 3.** The pure Bernal structure case: diffractograms for a crystallite with a height $L_c$ ~12 nm and two examples of $L_a$ dimensions, 10 nm (red color) and 20 nm (blue color), based on model B in Fig. 1. The crystallite is composed of similar entities (same BSC type - Bernal - and same height $L'_c$). For a given $L_a$, $L'_c$ is varied from $L'_c$ = $L_c$ (generating the sharpest 112 peak for example) down to $L'_c$ ~ 0.670 nm (2 layers) (generating the broadest 112 peak for example). Thus, there are 6 different $L'_c$ for each color (each $L_a$). Dark symbols in the inserts correspond to complementary calculations with $L_a$ = 30, 40, and 50 nm. FWHM<sub>hkl</sub> for 112 and 103 peaks were extracted and their plots in function of $L'_c$ were fitted with simple laws and reported in the inserts, and likewise for the areas $A_{hkl}$ of the peaks.</figcaption></figure>
crystallite model (blue), and the median value of the Gaussian $L_a$ dis- tribution for the multiple crystallite model (red). For all, $L_c$ = 10 nm. The blue curve (single-crystallite approach) is perfectly masked by the red curve (multiple-crystallite approach), meaning that both approaches are equivalent, at least when the size distribution (either through $L_a$, as in Fig. 2, or through $L_c$) is Gaussian. In the case of the AB pair BSC, one can note that, if only the 10 peak region is used, confusion with the Tur- bostratic BSC can be made, hence considering both the 10 and 11 peak regions is needed to discriminate between both BSCs.

What is true for the Gaussian-shaped peak 100 is also true for 002: a Gaussian size distribution gives a diffractogram comparable to that of the average value. As each diffractogram is normalized by the number of carbon atoms, a Gaussian distribution means more small crystallites than large crystallites.

One remark can be made about the peak indexation, which uses three Laue indices: since the beginning of the 20th century, because of the prevalence of the hexagonal structure, the peak related to the inter- graphene spacing in graphenic carbons has always been indexed 002 in all publications to make simple the notation when dealing with car- bon with mixed structures. However, this is inconsistent with the characteristics of the unit cell respective to each structure. Therefore, with the rhombohedral structure, the interlayer spacing is associated with the 003 peak (because there are three graphene layers in the unit cell) while, with the Bernal structure, it is written 002 (because there are two graphene layers in the unit cell). Likewise, in this work, it is written 001 for the turbostratic structure as the unit cell is made of a single

graphene layer only.

## 3.2. The single crystallite/multi-entities approach applied to a Bernal crystallite

Comparing the Bernal and rhombohedral graphite structures (Fig. 2), the 103 peak seems very interesting as this peak is specific to Bernal stacking, and likewise for the 102 peak in the Bernal diffractogram, although the peak intensity is less. The Bernal 112 peak is also inter- esting as it is located at the same position and with the same intensity as the rhombohedral (in that case, the indexation is 113 but corresponds to the same plane family), hence, it could be taken as a reference if needed.

In Fig. 3 are reported many diffractograms corresponding to a single crystallite model for two examples of $L_a$ values (10 and 20 nm, red and blue curves, respectively) and a single $L_c$ value (~12 nm, corresponding to 36 stacked layers). Therefore, for a given $L_a$, the diffractograms result from the same number of atoms. For each $L_a$ value, the average crys- tallite is made of *n* entities, all with the same BSC (Bernal) and same height $L'_c$, while the entities are turbostratically stacked to each other within the average crystallite (model B in Fig. 1). The various dif- fractograms for a given $L_a$ are obtained by varying $L'_c$ from $L'_c$ = $L_c$ down to $L'_c$ = 1 nm corresponding to a minimum of three stacked layers (two layers would correspond to the AB pair BSC type). In the SI is reported the diffractogram for the ABA 3-layer BSC (Fig. S2).

One can note that the 112 peak is clearly visible even with two layers (hence corresponding to the AB pair BSC). For both peaks of interest

<figure><img src="image_6.png" /><figcaption>**Fig. 3.** The pure Bernal structure case: diffractograms for a crystallite with a height $L_c$ ~12 nm and two examples of $L_a$ dimensions, 10 nm (red color) and 20 nm (blue color), based on model B in Fig. 1. The crystallite is composed of similar entities (same BSC type - Bernal - and same height $L'_c$). For a given $L_a$, $L'_c$ is varied from $L'_c$ = $L_c$ (generating the sharpest 112 peak for example) down to $L'_c$ ~ 0.670 nm (2 layers) (generating the broadest 112 peak for example). Thus, there are 6 different $L'_c$ for each color (each $L_a$). Dark symbols in the inserts correspond to complementary calculations with $L_a$ = 30, 40, and 50 nm. FWHM<sub>hkl</sub> for 112 and 103 peaks were extracted and their plots in function of $L'_c$ were fitted with simple laws and reported in the inserts, and likewise for the areas $A_{hkl}$ of the peaks.</figcaption></figure>
(112 and 103), the inserts give the FWHMs and the plots in function of L'c are fitted by considering simple laws. It is seen that they both relate to the value of L'_c but, as opposed to FWHM112, FWHM103 is only slightly dependent on L_a. Other inserts also report the area A_hkl of the peaks to prove that they are constant, as expected with the formula used for the calculation (see XaNSoNS documentation). It is also observed that 00l peak positions gradually shift towards higher 2θ angles as L'_c is decreased. This is because this decrease results in an increase in the contribution of the turbostratic stacking distance at 0.344 nm between the entities. It is meaningful that the 00l peaks remain symmetric while shifting in position, instead of generating an asymmetry towards lower 2θ angles. This observation can be considered as the signature of the fact that all the entities involved are actually in partial coherence, for belonging to the same average crystallite.

### 3.3. The single-crystallite/multiple-entities approach applied to a rhombohedral crystallite

In Fig. 4 are reported the diffractograms calculated for the rhombohedral structure following the same principle as for Fig. 3. Therefore, they also correspond to model B in Fig. 1, only the BSC type is different.

As stated from Fig. 2, the 105 peak is specific to the rhombohedral structure as it is located at a 2θ position where no feature exists for the hexagonal structure (Bernal). The 113 peak (corresponding to the 112 peak for Bernal) is also interesting as it relates to L'_c, as shown by the FWHM112 = f(L'_c) plot. In the inserts are given the specific parameters to

consider for fitting the peaks, i.e., a constant area and the functions to obtain the FWHMhkl. The remark regarding the 00l position variation is the same as for the Bernal structure.

### 3.4. First attempt to model the commercial graphite

Fig. 5 shows the experimental diffractogram obtained from the commercial graphite and two petroleum-pitch-based cokes annealed at 2100 and 2500 °C respectively. For the graphite, from the 100 peak and the Scherrer formula, a crystallite size L_a = 50 nm may be deduced. However, if the 110 peak is considered, another L_a value, 100 nm, is obtained. This latter value is more reliable for the reason that will be evidenced later on and will be the crystallite size used in the following.

Visually comparing the experimental diffractogram of the graphite sample in Fig. 5 with what a perfect graphite diffractogram should look like (see the calculated diffractogram in Fig. 2), shows that this commercial graphite cannot be merely described in terms of a single average crystallite comprising a single entity with the Bernal BSC. FWHM100 and FWHM101 are equal to 0.08 and 0.05°, respectively (for FWHM101, it is roughly equal to $\sqrt{\frac{20}{L_c^2} + \frac{17}{L_c^2}}$, as reported in the SI of [3] and in agreement with our calculations), which means that a clear separation down to the baseline between the 100 and 101 peaks was expected (see the red line in Fig. 6e for instance), which is not the case. Fig. 6 reports a first fitting attempt by considering an average crystallite of size L_a = 100 nm composed of several entities also with L_a = 100 nm and randomly

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** The pure rhombohedral structure case: diffractograms for a crystallite with a height L_c ~12 nm and two L_a dimensions, 10 nm (red color) and 20 nm (blue color), based on model B in Fig. 1. The crystallite is composed of n entities with the same BSC type (rhombohedral) and same height L'_c. For a given L_a, L'_c is varied from L'_c = L_c (generating the sharpest 105 peak for example) down to L'_c = ~1 nm (3 layers, generating the broadest 105 peak for example). Thus, there are 6 different L'_c for each color (each L_a). Dark symbols in the inserts correspond to complementary calculations with L_a = 30, 40, and 50 nm. FWHMhkl for 113 and 105 peaks were extracted and their plots in function of L'_c were fitted with simple laws and reported in the inserts. Same for the areas A_hkl of the peaks.</figcaption></figure>
(112 and 103), the inserts give the FWHMs and the plots in function of L'c are fitted by considering simple laws. It is seen that they both relate to the value of L'_c but, as opposed to FWHM112, FWHM103 is only slightly dependent on L_a. Other inserts also report the area A_hkl of the peaks to prove that they are constant, as expected with the formula used for the calculation (see XaNSoNS documentation). It is also observed that 00l peak positions gradually shift towards higher 2θ angles as L'_c is decreased. This is because this decrease results in an increase in the contribution of the turbostratic stacking distance at 0.344 nm between the entities. It is meaningful that the 00l peaks remain symmetric while shifting in position, instead of generating an asymmetry towards lower 2θ angles. This observation can be considered as the signature of the fact that all the entities involved are actually in partial coherence, for belonging to the same average crystallite.

### 3.3. The single-crystallite/multiple-entities approach applied to a rhombohedral crystallite

In Fig. 4 are reported the diffractograms calculated for the rhombohedral structure following the same principle as for Fig. 3. Therefore, they also correspond to model B in Fig. 1, only the BSC type is different.

As stated from Fig. 2, the 105 peak is specific to the rhombohedral structure as it is located at a 2θ position where no feature exists for the hexagonal structure (Bernal). The 113 peak (corresponding to the 112 peak for Bernal) is also interesting as it relates to L'_c, as shown by the FWHM112 = f(L'_c) plot. In the inserts are given the specific parameters to

consider for fitting the peaks, i.e., a constant area and the functions to obtain the FWHMhkl. The remark regarding the 00l position variation is the same as for the Bernal structure.

### 3.4. First attempt to model the commercial graphite

Fig. 5 shows the experimental diffractogram obtained from the commercial graphite and two petroleum-pitch-based cokes annealed at 2100 and 2500 °C respectively. For the graphite, from the 100 peak and the Scherrer formula, a crystallite size L_a = 50 nm may be deduced. However, if the 110 peak is considered, another L_a value, 100 nm, is obtained. This latter value is more reliable for the reason that will be evidenced later on and will be the crystallite size used in the following.

Visually comparing the experimental diffractogram of the graphite sample in Fig. 5 with what a perfect graphite diffractogram should look like (see the calculated diffractogram in Fig. 2), shows that this commercial graphite cannot be merely described in terms of a single average crystallite comprising a single entity with the Bernal BSC. FWHM100 and FWHM101 are equal to 0.08 and 0.05°, respectively (for FWHM101, it is roughly equal to $\sqrt{\frac{20}{L_c^2} + \frac{1.7}{L_c^2}}$, as reported in the SI of [3] and in agreement with our calculations), which means that a clear separation down to the baseline between the 100 and 101 peaks was expected (see the red line in Fig. 6e for instance), which is not the case. Fig. 6 reports a first fitting attempt by considering an average crystallite of size L_a = 100 nm composed of several entities also with L_a = 100 nm and randomly

<figure><img src="image_7.png" /><figcaption>**Fig. 4.** The pure rhombohedral structure case: diffractograms for a crystallite with a height L_c ~12 nm and two L_a dimensions, 10 nm (red color) and 20 nm (blue color), based on model B in Fig. 1. The crystallite is composed of n entities with the same BSC type (rhombohedral) and same height L'_c. For a given L_a, L'_c is varied from L'_c = L_c (generating the sharpest 105 peak for example) down to L'_c = ~1 nm (3 layers, generating the broadest 105 peak for example). Thus, there are 6 different L'_c for each color (each L_a). Dark symbols in the inserts correspond to complementary calculations with L_a = 30, 40, and 50 nm. FWHMhkl for 113 and 105 peaks were extracted and their plots in function of L'_c were fitted with simple laws and reported in the inserts. Same for the areas A_hkl of the peaks.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Fig. 5. (a) Experimental diffractogram of the commercial graphite (blue line), a pitch-coke annealed at 2100 °C (red line), and the same annealed at 2500 °C (magenta line). (b) Close-up of the 42–47° 2θ range. Small humps assigned to a minor contribution of rhombohedral stacking in the diffractogram of graphite are arrowed.</figcaption></figure>

rotated to each other. In this first attempt, all the entities are of the Bernal structure. Then, two approaches are considered: (i) the "single-size entity" approach, in which all the entities within the crystallite are identical in size, hence they exhibit identical $L'_c$ values (model B in Fig. 1); (ii) the "distributed-size entities" approach, in which the entities within the crystallite may exhibit a whole distribution of $L'_c$ values (model C in Fig. 1).

As the 2θ range around 100 and 101 peaks shows dramatic variations with respect to the type of structure present (see Fig. 2), better focus on the 103 (Bernal)/105 (rhombohedral) peak region. First of all, there is no 105 peak, meaning that no rhombohedral BSC type is present as large crystals in significant amount in the material (at least with $L'_c > 2$ nm, corresponding to 6 layers, which is the limit beyond which the 105 peak would discriminate from the background, see insert $FWHM_{105} = f(L'_c)$ in Fig. 4). If we consider only a single $L'_c$ value for all the crystallite entities, the fitting of the 103 and 112 peaks is not accurate (Fig. 6a and g), meaning that considering a Gaussian distribution of $L'_c$ values would not work either (as stated from Fig. 2 with the example of varying $L_a$). On the other hand, considering multiple entities with their $L'_c$ values distributed according to a Voigt function (Fig. 6d and j) provides a fairly good fitting of the 103 and 112 peaks (Fig. 6b and h), consistently with a large number of adjustable parameters. However, it can be noted that the related optimized $L'_c$ distributions (Fig. 6d and j), yet providing a same mean $L'_c$ value, are not identical, which cannot be the reality, as both pertain to the same material. Moreover, none of those two attempts using a single BSC type (Bernal, here), allows fitting the diffractogram in the 100–101 2θ region (Fig. 6f and l).

To account for the difference between the experimental diffractograms of the commercial graphite (Fig. 5) and that of ideal graphite (Fig. 2, second top), some contribution of AB pairs could be assumed. Such an assumption was actually necessary to fit pitch-based coke diffractograms whatever the annealing temperature [3] (two of them, at 2100 and 2500 °C, are shown in Fig. 5). However, that AB pairs remain in large amount (~50 %) at the highest annealing temperatures, say,

2500 °C and above, is doubtful (see Fig. 6a in [3]). Correspondingly, the presence of AB pairs should broaden the base of the 112 peak (as seen for the 2100 °C-annealed coke in Fig. 5 for instance), but this is not what is seen on the diffractograms of both the commercial graphite or the 2500 °C-annealed coke (Fig. 5). Another observation which suggests that considering that high temperature-annealed cokes are composed of pure Bernal and AB pair BSCs only is unsatisfactory to follow the plot of the $L_a$ values calculated from the 10(0) and 11(0) peaks respectively (Fig. S3 in SI). As soon as $L_a$ are larger than ~30 nm, they are clearly underestimated when calculated from the 10(0) peak. This means that, for the highest annealing temperatures, some structural feature should be present that affects the 10(0) peak but not the 11(0) peak. On the other hand, rhombohedral stacking is known to be often present in graphitized carbons [6,11], but a significant contribution of the rhombohedral BSC is excluded since the 105 peak is not formed. Alternatively, the occurrence of stacking faults as ABC triplets could be envisaged, the same way as the existence of AB pairs was revealed in the coke series [3,13] but they can be excluded here for the same reason (see the calculated diffractogram of ABC triplets in Fig. S2). Consequently, only stacking faults as randomly dispersed mispositioned graphene layers (C layers) should be introduced in the models, which makes full physical sense, as previously discussed in Section 2.1. One way to do so is to reconsider the ideality of the Bernal BSC, initially designed as intrinsically perfect, by introducing layers in C position within an ABAetc stacking sequence.

### 3.5. Introducing defects in Bernal crystallites with small $L'_c$ entities (single-crystallite approach)

To fully understand the effect of stacking faults as dispersed C layers, we first use the one-crystallite approach. The crystallite is made of n entities with same BSC (Bernal) and same $L'_c$. $L_c$ varies slightly so that to include an integer number of $L'_c$, hence corresponding to heights of ~33 to 39 layers), and $L_a = 20$ nm (which minimizes calculation times while
<figure><img src="image_8.png" /><figcaption>**Fig. 5.** (a) Experimental diffractogram of the commercial graphite (blue line), a pitch-coke annealed at 2100 °C (red line), and the same annealed at 2500 °C (magenta line). (b) Close-up of the 42–47° 2θ range. Small humps assigned to a minor contribution of rhombohedral stacking in the diffractogram of graphite are arrowed.</figcaption></figure>

rotated to each other. In this first attempt, all the entities are of the Bernal structure. Then, two approaches are considered: (i) the "single-size entity" approach, in which all the entities within the crystallite are identical in size, hence they exhibit identical $L'_c$ values (model B in Fig. 1); (ii) the "distributed-size entities" approach, in which the entities within the crystallite may exhibit a whole distribution of $L'_c$ values (model C in Fig. 1).

As the 2θ range around 100 and 101 peaks shows dramatic variations with respect to the type of structure present (see Fig. 2), better focus on the 103 (Bernal)/105 (rhombohedral) peak region. First of all, there is no 105 peak, meaning that no rhombohedral BSC type is present as large crystals in significant amount in the material (at least with $L'_c > 2$ nm, corresponding to 6 layers, which is the limit beyond which the 105 peak would discriminate from the background, see insert $FWHM_{105} = f(L'_c)$ in Fig. 4). If we consider only a single $L'_c$ value for all the crystallite entities, the fitting of the 103 and 112 peaks is not accurate (Fig. 6a and g), meaning that considering a Gaussian distribution of $L'_c$ values would not work either (as stated from Fig. 2 with the example of varying $L_a$). On the other hand, considering multiple entities with their $L'_c$ values distributed according to a Voigt function (Fig. 6d and j) provides a fairly good fitting of the 103 and 112 peaks (Fig. 6b and h), consistently with a large number of adjustable parameters. However, it can be noted that the related optimized $L'_c$ distributions (Fig. 6d and j), yet providing a same mean $L'_c$ value, are not identical, which cannot be the reality, as both pertain to the same material. Moreover, none of those two attempts using a single BSC type (Bernal, here), allows fitting the diffractogram in the 100–101 2θ region (Fig. 6f and l).

To account for the difference between the experimental diffractograms of the commercial graphite (Fig. 5) and that of ideal graphite (Fig. 2, second top), some contribution of AB pairs could be assumed. Such an assumption was actually necessary to fit pitch-based coke diffractograms whatever the annealing temperature [3] (two of them, at 2100 and 2500 °C, are shown in Fig. 5). However, that AB pairs remain in large amount (~50 %) at the highest annealing temperatures, say,

2500 °C and above, is doubtful (see Fig. 6a in [3]). Correspondingly, the presence of AB pairs should broaden the base of the 112 peak (as seen for the 2100 °C-annealed coke in Fig. 5 for instance), but this is not what is seen on the diffractograms of both the commercial graphite or the 2500 °C-annealed coke (Fig. 5). Another observation which suggests that considering that high temperature-annealed cokes are composed of pure Bernal and AB pair BSCs only is unsatisfactory to follow the plot of the $L_a$ values calculated from the 10(0) and 11(0) peaks respectively (Fig. S3 in SI). As soon as $L_a$ are larger than ~30 nm, they are clearly underestimated when calculated from the 10(0) peak. This means that, for the highest annealing temperatures, some structural feature should be present that affects the 10(0) peak but not the 11(0) peak. On the other hand, rhombohedral stacking is known to be often present in graphitized carbons [6,11], but a significant contribution of the rhombohedral BSC is excluded since the 105 peak is not formed. Alternatively, the occurrence of stacking faults as ABC triplets could be envisaged, the same way as the existence of AB pairs was revealed in the coke series [3,13] but they can be excluded here for the same reason (see the calculated diffractogram of ABC triplets in Fig. S2). Consequently, only stacking faults as randomly dispersed mispositioned graphene layers (C layers) should be introduced in the models, which makes full physical sense, as previously discussed in Section 2.1. One way to do so is to reconsider the ideality of the Bernal BSC, initially designed as intrinsically perfect, by introducing layers in C position within an ABAetc stacking sequence.

### 3.5. Introducing defects in Bernal crystallites with small $L'_c$ entities (single-crystallite approach)

To fully understand the effect of stacking faults as dispersed C layers, we first use the one-crystallite approach. The crystallite is made of n entities with same BSC (Bernal) and same $L'_c$. $L_c$ varies slightly so that to include an integer number of $L'_c$, hence corresponding to heights of ~33 to 39 layers), and $L_a = 20$ nm (which minimizes calculation times while
<figure><img src="image_8.png" /><figcaption>Fig. 6. (a) and (b) Two attempts for fitting (red) the 103 peak of the experimental diffractogram (blue) of the commercial graphite considering L<sub>a</sub> = 100 nm (see text) and the two cases of L'<sub>c</sub> value distribution respectively reported underneath: one (c) according to the "single-size entity" approach (see text, and model B in Fig. 1) which the optimization provides a L'<sub>c</sub> value of 7.23 nm, the other (d) according to the "distributed-size entities" approach (see text, and model C in Fig. 1), the optimization of which provides a distribution of L'<sub>c</sub> with a mean value at 3.70 nm. (e) and (f) show how the experimental diffractogram (blue) is fitted (red) in the 100/101 peak region considering the two approaches above-mentioned (in which the related functions were given in the SI of [3]). From (g) to (l): same as previously, now considering the 112 peak instead of the 103 peak. It appears that considering multiple entities with a distribution of L'<sub>c</sub> (as in (d) and (j)) instead of an identical value for all entities (as in (c) and (i)) allows obtaining a consistent mean L'<sub>c</sub> value (~3.7 nm) from both 103 and 112 peaks, and a better fitting of the peaks 103 (compare (b) with (a)), 112 (compare (h) with (g)), 100, and 101 (compare (f) with (e), and (l) with (k)). In (a), (b), (g), (h), "std" stands for "standard deviation".</figcaption></figure>

generating peaks narrow enough). Then stacking faults are introduced in the Bernal sequences. Various ways to do so are tested as reported in Fig. 7, namely, substitution, insertion, and so-called "block-shifting". Starting from, e.g., an ABABABA sequence, "substitution" is, for instance, to change it into a ABCBABA sequence; "insertion" is, for instance, to change it into an ABCABAB sequence; "block-shifting" is, for instance, to change it into an ABACBCB sequence.

In Fig. 8 is reported the effect, with respect to the diffractogram of the ideal Bernal BSC, of introducing one defect as a single stacking fault in ABABetc Bernal sequences of 7 layers, and building a crystallite of 35 layers with 5 of such defective entities. Crystallites of 36, 33, 39, and 34

layers were also considered so that they can be built with an integer number *n* of defective entities of 7, 9, 11, 13, and 17 layers respectively (model E in Fig. 1). By operating this way, the comparison between the diffractograms are not affected significantly by the difference in *L*<sub>c</sub> from a case to another. As there are *n* possible positions for the stacking fault in each case (as listed in Table S1 in SI), each of the related diffractograms in Fig. 8 results from the average over all the possible diffractograms obtained from all the possible positions for a given defect in entities of a given number of layers (among 7, 9, 11, 13, or 17).

The first observation is that, whatever the number of layers involved or the way the stacking fault is introduced, the 110-112 peak range is

<figure><img src="image_9.png" /><figcaption>Fig. 7. The three various ways considered in this work to introduce stacking faults in a Bernal crystallite.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Fig. 6. (a) and (b) Two attempts for fitting (red) the 103 peak of the experimental diffractogram (blue) of the commercial graphite considering L<sub>a</sub> = 100 nm (see text) and the two cases of L'<sub>c</sub> value distribution respectively reported underneath: one (c) according to the "single-size entity" approach (see text, and model B in Fig. 1) which the optimization provides a L'<sub>c</sub> value of 7.23 nm, the other (d) according to the "distributed-size entities" approach (see text, and model C in Fig. 1), the optimization of which provides a distribution of L'<sub>c</sub> with a mean value at 3.70 nm. (e) and (f) show how the experimental diffractogram (blue) is fitted (red) in the 100/101 peak region considering the two approaches above-mentioned (in which the related functions were given in the SI of [3]). From (g) to (l): same as previously, now considering the 112 peak instead of the 103 peak. It appears that considering multiple entities with a distribution of L'<sub>c</sub> (as in (d) and (j)) instead of an identical value for all entities (as in (c) and (i)) allows obtaining a consistent mean L'<sub>c</sub> value (~3.7 nm) from both 103 and 112 peaks, and a better fitting of the peaks 103 (compare (b) with (a)), 112 (compare (h) with (g)), 100, and 101 (compare (f) with (e), and (l) with (k)). In (a), (b), (g), (h), "std" stands for "standard deviation".</figcaption></figure>

generating peaks narrow enough). Then stacking faults are introduced in the Bernal sequences. Various ways to do so are tested as reported in Fig. 7, namely, substitution, insertion, and so-called "block-shifting". Starting from, e.g., an ABABABA sequence, "substitution" is, for instance, to change it into a ABCBABA sequence; "insertion" is, for instance, to change it into an ABCABAB sequence; "block-shifting" is, for instance, to change it into an ABACBCB sequence.

In Fig. 8 is reported the effect, with respect to the diffractogram of the ideal Bernal BSC, of introducing one defect as a single stacking fault in ABABetc Bernal sequences of 7 layers, and building a crystallite of 35 layers with 5 of such defective entities. Crystallites of 36, 33, 39, and 34

layers were also considered so that they can be built with an integer number *n* of defective entities of 7, 9, 11, 13, and 17 layers respectively (model E in Fig. 1). By operating this way, the comparison between the diffractograms are not affected significantly by the difference in *L*<sub>c</sub> from a case to another. As there are *n* possible positions for the stacking fault in each case (as listed in Table S1 in SI), each of the related diffractograms in Fig. 8 results from the average over all the possible diffractograms obtained from all the possible positions for a given defect in entities of a given number of layers (among 7, 9, 11, 13, or 17).

The first observation is that, whatever the number of layers involved or the way the stacking fault is introduced, the 110-112 peak range is

<figure><img src="image_10.png" /><figcaption>Fig. 7. The three various ways considered in this work to introduce stacking faults in a Bernal crystallite.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Effect of introducing stacking faults as C layers dispersed in Bernal sequences. Each calculated diffractogram is obtained from a crystallite of 33–39 layers of $L_a$ = 20 nm built with *n* entities made of a defective Bernal BSC of either 7, 9, 11, 13, or 17 layers each in which one of the 3 types of stacking fault considered (from top to bottom: substitution, insertion, and block-shifting) is introduced (model E in Fig. 1). Light gray corresponds to the diffractogram of pristine Bernal and in red to blue colors are the diffractograms with stacking faults. Each calculated diffractogram is averaged over all the possible stacking fault positions for each case (as listed in Table S1 in SI).</figcaption></figure>

not affected. This is a highly important point as it means that **110** and **112** peaks are the reliable ones from which $L_a$ and $L'$ can be respectively extracted. On the other hand, although the **103** peak is not affected by the C substitution either, a C insertion results in a splitting of the peak,

and block-shifting results in a decrease in its height. In the **100–101** peak range, two observations can be done:
<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Effect of introducing stacking faults as C layers dispersed in Bernal sequences. Each calculated diffractogram is obtained from a crystallite of 33–39 layers of $L_a$ = 20 nm built with *n* entities made of a defective Bernal BSC of either 7, 9, 11, 13, or 17 layers each in which one of the 3 types of stacking fault considered (from top to bottom: substitution, insertion, and block-shifting) is introduced (model E in Fig. 1). Light gray corresponds to the diffractogram of pristine Bernal and in red to blue colors are the diffractograms with stacking faults. Each calculated diffractogram is averaged over all the possible stacking fault positions for each case (as listed in Table S1 in SI).</figcaption></figure>

not affected. This is a highly important point as it means that **110** and **112** peaks are the reliable ones from which $L_a$ and $L'_c$ can be respectively extracted. On the other hand, although the **103** peak is not affected by the C substitution either, a C insertion results in a splitting of the peak,

and block-shifting results in a decrease in its height. In the **100–101** peak range, two observations can be done:
* in all cases, some intensity is generated in the valley between the 100 and 101 peaks, mostly at the expense of the 101 peak. Again, the perturbance is the most dramatic for the case of C insertion, generating a large base for the 100–101 block. These two features may be a signature of the occurrence of stacking faults within Bernal stacking. This goes well with the experimental diffractogram of the commercial graphite (see Fig. 5), in which (i) the height of 101 peak is significantly lower than its counterpart in perfect graphite (see the ideal Bernal BSC diffractogram in Fig. 2) with respect to the height of 100 peak; (ii) the base of the 101 peak is quite broad; (iii) small peaks (arrowed in Fig. 5) show up revealing the occurrence of rhombohedral sequences or as individual crystallites (however with height < ~6 layers, as the 105 peak is not well formed).

* a C insertion in crystallite entities with small $L'_c$ destroys the 101 peak.

In all cases, though, the perturbance is as less pronounced as the number of layers increases (keeping in mind that only one C layer is considered).

The results reported in Fig. 8 confirm that considering stacking faults as C layers dispersed in Bernal sequences is the right hypothesis, which makes more physical sense than considering the introduction of pure rhombohedral sequences. The effect of stacking faults may now be applied to large crystallites, hence including entities with large $L'_c$, more likely to correspond to the actual situation in graphitized materials.

### 3.6. Introducing defects in Bernal crystallites with large $L'_c$ entities (multiple-crystallite approach)

In Fig. 6, two cases were considered for calculated pure Bernal BSC diffractograms, either the single, average $L'_c$ case (Fig. 6c and i) or the distributed-$L'_c$-value case (Fig. 6d and j). A similar approach is now

being used while introducing stacking faults, in order to investigate the consequences on the 100–101 peak region (Fig. 9). From Fig. 6, it was found that the peak fitting was more accurate with considering a $L'_c$ distribution. However, it was pointed out that the best fitting corresponded to distributions that were different for 103 and for 112, which cannot be the reality. The reason may be provided by Fig. 7, which showed that 103 may be more or less affected by the presence of stacking faults, as opposed to 112. Therefore, we will consider here the $L'_c$ values determined from the 112 peak (Figs. 6i–j).

In Fig. 9a are reported the calculated diffractograms in the 100–101 peak region while introducing one of the stacking faults (red, blue, and green lines) within a single crystallite with a single entity (in that case $L_c = L'_c$), compared to ideal Bernal BSC (black line) and ideal rhombohedral BSC (magenta line). Calculations hence use model F in Fig. 1.

To build Fig. 9b, considering a single crystallite with enough layers (hundreds would be needed) so that to allow it to be made of enough entities with a variety of $L'_c$ distributed as shown in Fig. 6j would require unrealistic calculation times. Therefore, we switched to the multiple-crystallite approach (model H in Fig. 1). About one hundred crystallites were considered, made of 3 to 100 layers. Each crystallite has a height $L_c = L'_c$ since it is made of a single, defective entity, i.e. a Bernal BSC including one of the 3 types of stacking faults investigated. For a given stacking fault, each curve in Fig. 9b is then obtained from the sum of those crystallites with different $L_c$ (or $L'_c$), the proportions of which are that of the distribution shown in Fig. 6j (and again reported in insert in Fig. 9b). Because of the large number of crystallites considered for Fig. 9b, averaging the calculations over all the possible stacking fault positions as we did for Fig. 8 is impossible. Therefore, for both the insertion and block-shifting stacking faults, only one configuration was considered so that 25 % of the carbon atoms are mispositioned comparatively to pristine Bernal. For substitution, 10 % of mispositioned carbon atoms were considered, corresponding to changing 10 % of B

<figure><img src="image_11.png" /><figcaption>**Fig. 9.** 100 and 101 peak region: effect, on the calculated diffractograms ($L_a = 20$ nm), of a single position - arbitrarily selected - for the various types of stacking faults (red: substitution, blue: insertion, green: block-shifting) compared to pristine Bernal (black) and pristine rhombohedral (magenta). The crystallites exhibit the $L'_c (=L_c)$ values obtained from Fig. 6: (a) corresponds to the case that considered the best-fitting single average $L'_c$ value, based on model F in Fig. 1; (b) corresponds to the case that considered a Voigt distribution of the $L'_c = L_c$ values, based on model H in Fig. 1, averaged over the about hundred crystallites considered for each line (from 3 to 100 layers).</figcaption></figure>
* in all cases, some intensity is generated in the valley between the 100 and 101 peaks, mostly at the expense of the 101 peak. Again, the perturbance is the most dramatic for the case of C insertion, generating a large base for the 100–101 block. These two features may be a signature of the occurrence of stacking faults within Bernal stacking. This goes well with the experimental diffractogram of the commercial graphite (see Fig. 5), in which (i) the height of 101 peak is significantly lower than its counterpart in perfect graphite (see the ideal Bernal BSC diffractogram in Fig. 2) with respect to the height of 100 peak; (ii) the base of the 101 peak is quite broad; (iii) small peaks (arrowed in Fig. 5) show up revealing the occurrence of rhombohedral sequences or as individual crystallites (however with height < ~6 layers, as the 105 peak is not well formed).

* a C insertion in crystallite entities with small $L'_c$ destroys the 101 peak.

In all cases, though, the perturbance is as less pronounced as the number of layers increases (keeping in mind that only one C layer is considered).

The results reported in Fig. 8 confirm that considering stacking faults as C layers dispersed in Bernal sequences is the right hypothesis, which makes more physical sense than considering the introduction of pure rhombohedral sequences. The effect of stacking faults may now be applied to large crystallites, hence including entities with large $L'_c$, more likely to correspond to the actual situation in graphitized materials.

### 3.6. Introducing defects in Bernal crystallites with large $L'_c$ entities (multiple-crystallite approach)

In Fig. 6, two cases were considered for calculated pure Bernal BSC diffractograms, either the single, average $L'_c$ case (Fig. 6c and i) or the distributed-$L'_c$-value case (Fig. 6d and j). A similar approach is now

being used while introducing stacking faults, in order to investigate the consequences on the 100–101 peak region (Fig. 9). From Fig. 6, it was found that the peak fitting was more accurate with considering a $L'_c$ distribution. However, it was pointed out that the best fitting corresponded to distributions that were different for 103 and for 112, which cannot be the reality. The reason may be provided by Fig. 7, which showed that 103 may be more or less affected by the presence of stacking faults, as opposed to 112. Therefore, we will consider here the $L'_c$ values determined from the 112 peak (Figs. 6i–j).

In Fig. 9a are reported the calculated diffractograms in the 100–101 peak region while introducing one of the stacking faults (red, blue, and green lines) within a single crystallite with a single entity (in that case $L_c = L'_c$), compared to ideal Bernal BSC (black line) and ideal rhombohedral BSC (magenta line). Calculations hence use model F in Fig. 1.

To build Fig. 9b, considering a single crystallite with enough layers (hundreds would be needed) so that to allow it to be made of enough entities with a variety of $L'_c$ distributed as shown in Fig. 6j would require unrealistic calculation times. Therefore, we switched to the multiple-crystallite approach (model H in Fig. 1). About one hundred crystallites were considered, made of 3 to 100 layers. Each crystallite has a height $L_c = L'_c$ since it is made of a single, defective entity, i.e. a Bernal BSC including one of the 3 types of stacking faults investigated. For a given stacking fault, each curve in Fig. 9b is then obtained from the sum of those crystallites with different $L_c$ (or $L'_c$), the proportions of which are that of the distribution shown in Fig. 6j (and again reported in insert in Fig. 9b). Because of the large number of crystallites considered for Fig. 9b, averaging the calculations over all the possible stacking fault positions as we did for Fig. 8 is impossible. Therefore, for both the insertion and block-shifting stacking faults, only one configuration was considered so that 25 % of the carbon atoms are mispositioned comparatively to pristine Bernal. For substitution, 10 % of mispositioned carbon atoms were considered, corresponding to changing 10 % of B

<figure><img src="image_12.png" /><figcaption>**Fig. 9.** 100 and 101 peak region: effect, on the calculated diffractograms ($L_a$ = 20 nm), of a single position - arbitrarily selected - for the various types of stacking faults (red: substitution, blue: insertion, green: block-shifting) compared to pristine Bernal (black) and pristine rhombohedral (magenta). The crystallites exhibit the $L'_c$ ($L_c$) values obtained from Fig. 6: (a) corresponds to the case that considered the best-fitting single average $L'_c$ value, based on model F in Fig. 1; (b) corresponds to the case that considered a Voigt distribution of the $L'_c = L_c$ values, based on model H in Fig. 1, averaged over the about hundred crystallites considered for each line (from 3 to 100 layers).</figcaption></figure>
layers by C layers. For consistency, we did the same for Fig. 9a.

The 100 peak is affected in position and linewidth by the presence of stacking faults. Stacking faults also make the 101 peak systematically lower in height and with a larger base, which reminds the features of the experimental diffractogram of the commercial graphite (Fig. 5). More in the detail, we see that, again, considering a distribution of L'c = Lc (Fig. 9b) is more suitable as it is more likely to generate some intensity in the valley between the 100 and 101 peaks, whatever the defect type, than the "single-crystallite" approach (Fig. 9a). Thus, a combination of stacking faults can explain the shape of the experimental diffractogram of the commercial graphite, in particular in the 100–101 peak region as shown in the left insert in Fig. 10.

In Fig. 10 is reported, as in Fig. 5, the diffractogram of the commercial graphite. L<sub>a</sub>, L<sub>c</sub>, and L'<sub>c</sub> determined from both 100 and 110, 002, and 112 respectively are provided. Using 100 underestimates L<sub>a</sub> by a factor of 2 compared to the value obtained from 110, consistently with Fig. S3 in SI. The reason is obvious from considering the right insert, where both the 100 and 110 peaks are reported using the abscissa 2θ × cos(θ)/Kλ. With such a scale, the FWHMs are directly proportional to 1/L<sub>a</sub> considering the Scherrer law. In addition to the width difference, the shape of 110 is symmetric while that of 100 is not.

In the discussion above, only the cases of stacking faults were able to
be considered, and in addition to this, they were considered indepen-
dently from each other. In practice, that the effect of these various cases
of stacking faults, and others, are combined in experimental diffracto-
grams is about certain. Therefore, in most cases, identifying the defects
present in various graphitized materials by our methodology can only be
approached so far; instead of being precisely and exhaustively deter-
mined, because possible combinations are countless, and considering all
of them is impossible. But it is believed that, in the current state, the
methodology proposed here can be satisfactory, and is much more ac-
curate than any other previous ones.

The *00l* peaks associated with the atomic plane family perpendicular
to the *c* direction of the crystallite are also worth some attention. They

all exhibit an asymmetric shape (Fig. 11). To fit the experimental peaks
and obtain the L<sub>c</sub> of the average crystallite, several approaches are
presented in the literature. The standard approach is to consider the
linewidth of the peaks [16]. The decomposition into symmetric peaks is
another option, where a mixture of crystallites with various inter-
graphene distances (IGDs) arbitrarily selected is considered [17]. Fitting
each *00l* peak with a complex function is always possible but the pa-
rameters used could be physically meaningless. We tried many strategies
which resulted in some variations in ⟨L<sub>c</sub>⟩ and ⟨L'<sub>c</sub>⟩ meaning that these
parameters obtained from the fitting are indicative. But all the strategies
aiming at explaining the peak asymmetry tended to the same conclu-
sions. As we have already shown, for a given number of layers, the shape
of the *00l* peak remains symmetric whatever the Turbostratic/Bernal
proportions in the crystallite, however with a variation of the 2θ position
of the *00l* peak (see Fig. 3 in [3]). Therefore, the reason for the peak
asymmetry may only be found in the distribution of the crystallites. To
stay close to physical reality, here is shown the case in which the ma-
terial is made of a distribution of crystallites with different L<sub>c</sub>, and each
of them possibly made of several entities with the Bernal BSC (model I in
Fig. 1). The *00l* peak characteristics of the crystallites are related to two
parameters, namely the average ⟨IGD⟩ and the crystallite size L<sub>c</sub>. ⟨IGD⟩
for each crystallite may slightly vary because of the variable contribu-
tion of the turbostratic stacking between the entities, or other reasons
responsible to IGD larger than the expected 0.335 nm of ideal graphite,
such as the occurrence of screw dislocations [18]. We consider a dis-
tribution of ⟨IGD⟩ with the associated distribution of crystallite size L<sub>c</sub>. L<sub>c</sub>
is inversely proportional to FWHM<sub>00l</sub>. Both ⟨IGD⟩ and FWHM<sub>00l</sub> pa-
rameters are adjustable positive numbers which vary linearly between
two limit-values without any constraint. The result obtained from an
optimized crystallite distribution (i.e., which is able to fit the three *00l*
peaks at best) succeeded in accounting for each peak asymmetry.
However, a further accurate fitting (Fig. 11) was obtained by consid-
ering the possibility of an experimental mispositioning of 0.33 mm of the
sample with respect to the ideal eucentric position in the X-ray chamber

<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Experimental diffractogram of the commercial graphite. The left insert is a close-up on the 100/101 peak range. The right insert represents the 100 and 110 peaks with abscissa 2θ × cos(θ)/Kλ (derived from the Scherrer law), which makes that the FWHM is directly proportional to 1/L<sub>a</sub>.</figcaption></figure>
layers by C layers. For consistency, we did the same for Fig. 9a.

The 100 peak is affected in position and linewidth by the presence of stacking faults. Stacking faults also make the 101 peak systematically lower in height and with a larger base, which reminds the features of the experimental diffractogram of the commercial graphite (Fig. 5). More in the detail, we see that, again, considering a distribution of L'c = Lc (Fig. 9b) is more suitable as it is more likely to generate some intensity in the valley between the 100 and 101 peaks, whatever the defect type, than the "single-crystallite" approach (Fig. 9a). Thus, a combination of stacking faults can explain the shape of the experimental diffractogram of the commercial graphite, in particular in the 100–101 peak region as shown in the left insert in Fig. 10.

In Fig. 10 is reported, as in Fig. 5, the diffractogram of the commercial graphite. L<sub>a</sub>, L<sub>c</sub>, and L'<sub>c</sub> determined from both 100 and 110, 002, and 112 respectively are provided. Using 100 underestimates L<sub>a</sub> by a factor of 2 compared to the value obtained from 110, consistently with Fig. S3 in SI. The reason is obvious from considering the right insert, where both the 100 and 110 peaks are reported using the abscissa 2θ × cos(θ)/Kλ. With such a scale, the FWHMs are directly proportional to 1/L<sub>a</sub> considering the Scherrer law. In addition to the width difference, the shape of 110 is symmetric while that of 100 is not.

In the discussion above, only the cases of stacking faults were able to
be considered, and in addition to this, they were considered indepen-
dently from each other. In practice, that the effect of these various cases
of stacking faults, and others, are combined in experimental diffracto-
grams is about certain. Therefore, in most cases, identifying the defects
present in various graphitized materials by our methodology can only be
approached so far; instead of being precisely and exhaustively deter-
mined, because possible combinations are countless, and considering all
of them is impossible. But it is believed that, in the current state, the
methodology proposed here can be satisfactory, and is much more ac-
curate than any other previous ones.

The *00l* peaks associated with the atomic plane family perpendicular
to the *c* direction of the crystallite are also worth some attention. They

all exhibit an asymmetric shape (Fig. 11). To fit the experimental peaks
and obtain the L<sub>c</sub> of the average crystallite, several approaches are
presented in the literature. The standard approach is to consider the
linewidth of the peaks [16]. The decomposition into symmetric peaks is
another option, where a mixture of crystallites with various inter-
graphene distances (IGDs) arbitrarily selected is considered [17]. Fitting
each *00l* peak with a complex function is always possible but the pa-
rameters used could be physically meaningless. We tried many strategies
which resulted in some variations in ⟨L<sub>c</sub>⟩ and ⟨L'<sub>c</sub>⟩ meaning that these
parameters obtained from the fitting are indicative. But all the strategies
aiming at explaining the peak asymmetry tended to the same conclu-
sions. As we have already shown, for a given number of layers, the shape
of the *00l* peak remains symmetric whatever the Turbostratic/Bernal
proportions in the crystallite, however with a variation of the 2θ position
of the *00l* peak (see Fig. 3 in [3]). Therefore, the reason for the peak
asymmetry may only be found in the distribution of the crystallites. To
stay close to physical reality, here is shown the case in which the ma-
terial is made of a distribution of crystallites with different L<sub>c</sub>, and each
of them possibly made of several entities with the Bernal BSC (model I in
Fig. 1). The *00l* peak characteristics of the crystallites are related to two
parameters, namely the average ⟨IGD⟩ and the crystallite size L<sub>c</sub>. ⟨IGD⟩
for each crystallite may slightly vary because of the variable contribu-
tion of the turbostratic stacking between the entities, or other reasons
responsible to IGD larger than the expected 0.335 nm of ideal graphite,
such as the occurrence of screw dislocations [18]. We consider a dis-
tribution of ⟨IGD⟩ with the associated distribution of crystallite size L<sub>c</sub>. L<sub>c</sub>
is inversely proportional to FWHM<sub>00l</sub>. Both ⟨IGD⟩ and FWHM<sub>00l</sub> pa-
rameters are adjustable positive numbers which vary linearly between
two limit-values without any constraint. The result obtained from an
optimized crystallite distribution (i.e., which is able to fit the three *00l*
peaks at best) succeeded in accounting for each peak asymmetry.
However, a further accurate fitting (Fig. 11) was obtained by consid-
ering the possibility of an experimental mispositioning of 0.33 mm of the
sample with respect to the ideal eucentric position in the X-ray chamber

<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Experimental diffractogram of the commercial graphite. The left insert is a close-up on the 100/101 peak range. The right insert represents the 100 and 110 peaks with abscissa 2θ × cos(θ)/Kλ (derived from the Scherrer law), which makes that the FWHM is directly proportional to 1/L<sub>a</sub>.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Fig. 11.** From left to right: fits (red lines) of the experimental 002, 004 and 006 peaks from the commercial graphite (blue lines) considering the distribution of crystallites which fits the peaks the best (see the insert in the left figure), combined with a mispositioning of 0.33 mm of the sample with respect to the eucentric position. The value of $\langle L_c \rangle$ used to fit the three peaks is the same, as deduced from the crystallite size distribution, to be compared with the $L_c(\text{FWHM}_{00\text{l}})$ values reported in each figure as obtained for each peak by using the standard Scherrer law considering $K = 0.89$.</figcaption></figure>

(Fig. S4).

Fig. 11 shows one possibility among several which were obtained by varying the starting parameter values, which roughly provided satisfactory results, provided some eucentric mispositioning of the sample (here of 0.33 mm) was considered to achieve a fairly good fitting. However, it is worth noting that, in none of them, 0.335 nm (the ideal IGD in genuine graphite) was the most represented value of the $\langle IGD \rangle$ distribution, indicating that this natural graphite is far from ideal.

From the distribution of $\langle IGD \rangle$ values, it is possible to calculate a $\langle L'_{c} \rangle$ value (see SI) which equals 4.1 nm (otherwise ranging from 3.4 to 4.3 nm in our several tries) while $\langle L_c \rangle$ equals 46.7 nm (otherwise ranging from 38 to 47 nm in our several tries). This $\langle L_c \rangle$ is in the same range and a bit larger than the value $L_c(\text{FWHM}_{002})$ of 32.8 nm obtained using the Scherrer law on the most intense part of the 002 peak, which makes a consistent result. On the other hand, this $\langle L'_{c} \rangle$ value obtained from $\langle IGD \rangle$ is also consistent with the 3.7 nm calculated from the 112 peak. Two statements are worth making about those results: (i) the $L_c$ and $L'_{c}$ values found for the natural graphite are close to that of a 2650 °C-annealed coke (29 and 2.5 nm, respectively [3]), despite their precursor and thermal history are quite different; (ii) it was not necessary to consider stacking faults in the modeling (model I in Fig. 1) to achieve convenient fitting of the 00l peaks, which confirms – logically – that stacking faults do not affect the reciprocal elements involving the c direction in the diffractograms.

Finally, back to the whole graphite diffractogram, it is likely that the remaining discrepancies with respect to perfect fitting could be assigned to the presence of a specific combination of stacking faults, or to other defects not investigated here (e.g., dislocations). Unfortunately, the number of defects and possible configurations is again too huge for being explored thoroughly by mean of calculations using the means currently available. Interestingly, in the coke series, the annealing temperature at which the $L_a$ values obtained by respectively applying the Scherrer law to the 100 and 110 peaks start to differ significantly (Fig. S3 in SI) indicates the temperature threshold at which stacking faults start to occur and contribute significantly to the diffractograms. This threshold appears to be ~2300 °C for the coal-tar-pitch coke series, corresponding to the temperature at which $L_a$ = ~30 nm is reached (Fig. 7a in [3]).

## 4. Conclusions

Bottom-up modeling was used once via a single-crystallite approach to model crystallized graphenic systems with millions of atoms involved in non-defective Bernal stacking [3,13]. In contrast with analyzing X-ray

diffractograms by merely using the Scherrer law, which comes with uncertainties and inconsistencies due to the complexity of graphenic materials at atom scale, such a bottom-up modeling enables to find a consistent way to reliably obtain the averaged well-known $L_a$ and $L_c$ of the crystallites that the material is made of. It also provides $\langle L'_{c} \rangle$, which applies to the average height of the coherent entities stacked randomly (= turbostratically) within the average crystallite.

It is shown here that bottom-up modeling can also be used to model large, defective systems provided a multiple-crystallite approach is used. Based on two examples, it is shown that the X-ray diffractograms of highly graphitized carbon materials, either natural or industrial, exhibit features which make them different from that of genuine graphite.

It is demonstrated that the diffractogram alterations with respect to the genuine graphite (Bernal) structure relate to the occurrence of stacking faults, with respect to the A and B positions corresponding to the hexagonal structure. Such stacking faults are shown to correspond to the presence of graphene layers in C position randomly dispersed within Bernal sequences, preferably to any other possibility such as ABC triplets or full rhombohedral sequences.

Three dispersion configuration for C layers are analyzed, namely insertion, substitution, and block shifting, which alter the graphite diffractogram in different ways. The 110-112 2θ-range is not sensitive to such stacking faults, therefore, meaningful $\langle L_a \rangle$ and $\langle L'_{c} \rangle$ values can be obtained from the 110 and $101_{\text{Bernal}}$ peaks respectively. On the other hand, the 100-101 2θ-range is a fingerprint of the stacking faults in the crystallite as interference effects strongly affect the 100 and $101_{\text{Bernal}}$ peaks and result in significantly increasing the background in the valley between 100 and 101 peaks. The $101_{\text{Bernal}}$ peak is disturbed the most by the insertion of a C layer in the Bernal sequence, as it makes the crystallite correspond to the merge of 2 crystallites with misaligned stacking. As a consequence of the resulting destructive interference effects, the 101 peak is both lowered and broadened. The systematic inconsistency obtained when calculating $L_a$ using the Scherrer law on the 100 and 110 peaks respectively is thereby explained.

In case 00l peaks exhibit some asymmetry towards lower 2θ angles, an explanation is proposed, which is that the material is actually made of a distribution of crystallites of various $L_c$, each of them constituted of Bernal entities with various $L'_{c}$. This enables to obtain a $\langle L_c \rangle$ value which accounts for the three peaks, 002, 004, and 006, whereas merely using the Scherrer law provides inconsistent values for each. Incorporating a small eucentric mispositioning of the sample in the X-ray chamber as an adjustable parameter improves the fitting. Therefore, when performing X-ray diffraction experiments, it is recommended to adjust the specimen eucentric height accurately.
<figure><img src="image_14.png" /><figcaption>**Fig. 11.** From left to right: fits (red lines) of the experimental 002, 004 and 006 peaks from the commercial graphite (blue lines) considering the distribution of crystallites which fits the peaks the best (see the insert in the left figure), combined with a mispositioning of 0.33 mm of the sample with respect to the eucentric position. The value of $\langle L_c \rangle$ used to fit the three peaks is the same, as deduced from the crystallite size distribution, to be compared with the $L_c(\text{FWHM}_{00\text{l}})$ values reported in each figure as obtained for each peak by using the standard Scherrer law considering $K = 0.89$.</figcaption></figure>

(Fig. S4).

Fig. 11 shows one possibility among several which were obtained by varying the starting parameter values, which roughly provided satisfactory results, provided some eucentric mispositioning of the sample (here of 0.33 mm) was considered to achieve a fairly good fitting. However, it is worth noting that, in none of them, 0.335 nm (the ideal IGD in genuine graphite) was the most represented value of the $\langle IGD \rangle$ distribution, indicating that this natural graphite is far from ideal.

From the distribution of $\langle IGD \rangle$ values, it is possible to calculate a $\langle L'_{c} \rangle$ value (see SI) which equals 4.1 nm (otherwise ranging from 3.4 to 4.3 nm in our several tries) while $\langle L_c \rangle$ equals 46.7 nm (otherwise ranging from 38 to 47 nm in our several tries). This $\langle L_c \rangle$ is in the same range and a bit larger than the value $L_c(\text{FWHM}_{002})$ of 32.8 nm obtained using the Scherrer law on the most intense part of the 002 peak, which makes a consistent result. On the other hand, this $\langle L'_{c} \rangle$ value obtained from $\langle IGD \rangle$ is also consistent with the 3.7 nm calculated from the 112 peak. Two statements are worth making about those results: (i) the $L_c$ and $L'_{c}$ values found for the natural graphite are close to that of a 2650 °C-annealed coke (29 and 2.5 nm, respectively [3]), despite their precursor and thermal history are quite different; (ii) it was not necessary to consider stacking faults in the modeling (model I in Fig. 1) to achieve convenient fitting of the 00l peaks, which confirms – logically – that stacking faults do not affect the reciprocal elements involving the c direction in the diffractograms.

Finally, back to the whole graphite diffractogram, it is likely that the remaining discrepancies with respect to perfect fitting could be assigned to the presence of a specific combination of stacking faults, or to other defects not investigated here (e.g., dislocations). Unfortunately, the number of defects and possible configurations is again too huge for being explored thoroughly by mean of calculations using the means currently available. Interestingly, in the coke series, the annealing temperature at which the $L_a$ values obtained by respectively applying the Scherrer law to the 100 and 110 peaks start to differ significantly (Fig. S3 in SI) indicates the temperature threshold at which stacking faults start to occur and contribute significantly to the diffractograms. This threshold appears to be ~2300 °C for the coal-tar-pitch coke series, corresponding to the temperature at which $L_a$ = ~30 nm is reached (Fig. 7a in [3]).

## 4. Conclusions

Bottom-up modeling was used once via a single-crystallite approach to model crystallized graphenic systems with millions of atoms involved in non-defective Bernal stacking [3,13]. In contrast with analyzing X-ray

diffractograms by merely using the Scherrer law, which comes with uncertainties and inconsistencies due to the complexity of graphenic materials at atom scale, such a bottom-up modeling enables to find a consistent way to reliably obtain the averaged well-known $L_a$ and $L_c$ of the crystallites that the material is made of. It also provides $\langle L'_{c} \rangle$, which applies to the average height of the coherent entities stacked randomly (= turbostratically) within the average crystallite.

It is shown here that bottom-up modeling can also be used to model large, defective systems provided a multiple-crystallite approach is used. Based on two examples, it is shown that the X-ray diffractograms of highly graphitized carbon materials, either natural or industrial, exhibit features which make them different from that of genuine graphite.

It is demonstrated that the diffractogram alterations with respect to the genuine graphite (Bernal) structure relate to the occurrence of stacking faults, with respect to the A and B positions corresponding to the hexagonal structure. Such stacking faults are shown to correspond to the presence of graphene layers in C position randomly dispersed within Bernal sequences, preferably to any other possibility such as ABC triplets or full rhombohedral sequences.

Three dispersion configuration for C layers are analyzed, namely insertion, substitution, and block shifting, which alter the graphite diffractogram in different ways. The 110-112 $2\theta$-range is not sensitive to such stacking faults, therefore, meaningful $\langle L_a \rangle$ and $\langle L'_{c} \rangle$ values can be obtained from the 110 and $101_{\text{Bernal}}$ peaks respectively. On the other hand, the 100-101 $2\theta$-range is a fingerprint of the stacking faults in the crystallite as interference effects strongly affect the 100 and $101_{\text{Bernal}}$ peaks and result in significantly increasing the background in the valley between 100 and 101 peaks. The $101_{\text{Bernal}}$ peak is disturbed the most by the insertion of a C layer in the Bernal sequence, as it makes the crystallite correspond to the merge of 2 crystallites with misaligned stacking. As a consequence of the resulting destructive interference effects, the 101 peak is both lowered and broadened. The systematic inconsistency obtained when calculating $L_a$ using the Scherrer law on the 100 and 110 peaks respectively is thereby explained.

In case 00l peaks exhibit some asymmetry towards lower $2\theta$ angles, an explanation is proposed, which is that the material is actually made of a distribution of crystallites of various $L_c$, each of them constituted of Bernal entities with various $L'_{c}$. This enables to obtain a $\langle L_c \rangle$ value which accounts for the three peaks, 002, 004, and 006, whereas merely using the Scherrer law provides inconsistent values for each. Incorporating a small eucentric mispositioning of the sample in the X-ray chamber as an adjustable parameter improves the fitting. Therefore, when performing X-ray diffraction experiments, it is recommended to adjust the specimen eucentric height accurately.
Overall, if using the multiple-crystallite approach in bottom-up modeling allows revealing the existence of stacking faults, identifying which types of stacking faults are actually present in the material is nowadays impossible because making too heavy calculations so far. But we are confident that it will be possible in the future, by involving artificial intelligence in the calculation process.

## Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Data availability

Data will be made available on request.

## Acknowledgements

Agnieszka Dabrowska (University of Warsaw, Poland) is thanked for carrying-out the XRD on the cokes at CEMES.

## Supplementary materials

Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.cartre.2023.100311.

## References

[1] E. Fitzer, K.H. Kochling, K.H. Boehm, H. Marsh, H. Recommended terminology for the description of carbon as a solid (IUPAC Recommendations 1995), Pure Appl. Chem. 67 (1995) 473–506.

[2] B.E. Warren, X-ray diffraction in random layer lattices, Phys. Rev. 59 (1941) 693–698.

[3] P. Puech, A. Dabrowska, N. Ratel-Ramond, G.L. Vignoles, M. Monthoux, New insight on carbonisation and graphitisation mechanisms as obtained from a bottom-up analytical approach of X-ray diffraction patterns, Carbon N Y 147 (2019) 602–611.

[4] A. Bianco, H.M. Cheng, T. Enoki, Y. Gogotsi, R.H. Hurt, N. Koratkar, T. Kyotani, M. Monthoux, C.R. Park, J.M.D. Tascon, J. Zhang, All in the graphene family–A recommended nomenclature for two-dimensional carbon materials, Carbon N Y 65 (2013) 1–6.

[5] J.P. Nery, M. Calandra, F. Mauri, Ab-initio energetics of graphite and multilayer graphene: stability of Bernal versus rhombohedral stacking, 2D Mater. 8 (2021), 035006.

[6] C.E. Precker, P.D. Esquinazi, A. Champi, J. Barzola-Quiquia, M. Zoraghi, S. Muinos-Landin, A. Setzer, W. Böhlmann, D. Spemann, J. Meijer, T. Muenster, O. Baehre, G. Kloess, H. Beth, Identification of a possible superconducting transition above room temperature in natural graphite crystals, New J. Phys. 18 (2016), 113041.

[7] S.L. Ramos, M.A. Pimenta, A. Champi, Multiple-excitation study of the double-resonance Raman bands in rhombohedral graphite, Carbon N Y 179 (2021) 683–691.

[8] A. Milev, M. Wilson, G.K. Kannangara, N. Tran, X-ray diffraction line profile analysis of nanocrystalline graphite, Mater. Chem. Phys. 111 (2008) 346–350.

[9] H. Fujimoto, Theoretical X-ray scattering intensity of carbons with turbostratic stacking and AB stacking structures, Carbon N Y 41 (2003) 1585–1592.

[10] K.J. Putman, M.R. Rowles, N.A. Marks, C. de Tomas, J.W. Martin, I. Suarez-Martinez, Defining graphenic crystallites in disordered carbon: moving beyond the platelet model, Carbon N Y 209 (2023), 117965.

[11] Q.-Y. Lin, T.-Q. Li, Z.-J. Liu, Y. Song, L.-L. He, Z.-J. Hu, Q.-G. Guo, H.-Q. Ye, High-resolution TEM observations of isolated rhombohedral crystallites in graphite blocks, Carbon N Y 50 (2012) 2347–2374.

[12] V.S. Neverov, XaNSoNS: gPU-accelerated simulator of diffraction patterns of nanoparticles, SoftwareX 6 (2017) 63–68.

[13] P. Puech, M. Monthoux, A new insight on the understanding of carbonisation and graphitisation mechanisms, Indian J. Eng. Mater. Sci. 27 (2020) 1095–1099.

[14] P. Scherrer, Bestimmung dergröße und derinneren struktur von kolloidteilchen mittels röntgenstrahlen, Nachrichten von der Königliche Gesellschaft der Wissenschaften zu Göttingen, Mathematisch-Physikalische Klasse (1918) 98–100.

[15] B.D. Cullity, Elements of X-ray Diffraction, 2nd edition, Addison-Wesley, 1978.

[16] H. Li, C. Yang, F. Liu, Novel method for determining stacking disorder degree in hexagonal graphite by X-ray diffraction, Sci. China Ser. B: Chem. 52 (2009) 174–180.

[17] J.B. Aladekomo, R.H. Bragg, Structural transformations induced in graphite by grinding: analysis of 002 X-ray diffraction line profiles, Carbon N Y 28 (1990) 897–906.

[18] J.W. Martin, J.L. Fogg, K.J. Putman, G. Francas, E.P. Turner, N.A. Marks, I. Suarez-Martinez, Graphite forms via annihilation of screw dislocations. arXiv: 2206.09105v (2022).
Overall, if using the multiple-crystallite approach in bottom-up modeling allows revealing the existence of stacking faults, identifying which types of stacking faults are actually present in the material is nowadays impossible because making too heavy calculations so far. But we are confident that it will be possible in the future, by involving artificial intelligence in the calculation process.

## Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Data availability

Data will be made available on request.

## Acknowledgements

Agnieszka Dabrowska (University of Warsaw, Poland) is thanked for carrying-out the XRD on the cokes at CEMES.

## Supplementary materials

Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.cartre.2023.100311.

## References

[1] E. Fitzer, K.H. Kochling, K.H. Boehm, H. Marsh, H. Recommended terminology for the description of carbon as a solid (IUPAC Recommendations 1995), Pure Appl. Chem. 67 (1995) 473–506.

[2] B.E. Warren, X-ray diffraction in random layer lattices, Phys. Rev. 59 (1941) 693–698.

[3] P. Puech, A. Dabrowska, N. Ratel-Ramond, G.L. Vignoles, M. Monthoux, New insight on carbonisation and graphitisation mechanisms as obtained from a bottom-up analytical approach of X-ray diffraction patterns, Carbon N Y 147 (2019) 602–611.

[4] A. Bianco, H.M. Cheng, T. Enoki, Y. Gogotsi, R.H. Hurt, N. Koratkar, T. Kyotani, M. Monthoux, C.R. Park, J.M.D. Tascon, J. Zhang, All in the graphene family–A recommended nomenclature for two-dimensional carbon materials, Carbon N Y 65 (2013) 1–6.

[5] J.P. Nery, M. Calandra, F. Mauri, Ab-initio energetics of graphite and multilayer graphene: stability of Bernal versus rhombohedral stacking, 2D Mater. 8 (2021), 035006.

[6] C.E. Precker, P.D. Esquinazi, A. Champi, J. Barzola-Quiquia, M. Zoraghi, S. Muinos-Landin, A. Setzer, W. Böhlmann, D. Spemann, J. Meijer, T. Muenster, O. Baehre, G. Kloess, H. Beth, Identification of a possible superconducting transition above room temperature in natural graphite crystals, New J. Phys. 18 (2016), 113041.

[7] S.L. Ramos, M.A. Pimenta, A. Champi, Multiple-excitation study of the double-resonance Raman bands in rhombohedral graphite, Carbon N Y 179 (2021) 683–691.

[8] A. Milev, M. Wilson, G.K. Kannangara, N. Tran, X-ray diffraction line profile analysis of nanocrystalline graphite, Mater. Chem. Phys. 111 (2008) 346–350.

[9] H. Fujimoto, Theoretical X-ray scattering intensity of carbons with turbostratic stacking and AB stacking structures, Carbon N Y 41 (2003) 1585–1592.

[10] K.J. Putman, M.R. Rowles, N.A. Marks, C. de Tomas, J.W. Martin, I. Suarez-Martinez, Defining graphenic crystallites in disordered carbon: moving beyond the platelet model, Carbon N Y 209 (2023), 117965.

[11] Q.-Y. Lin, T.-Q. Li, Z.-J. Liu, Y. Song, L.-L. He, Z.-J. Hu, Q.-G. Guo, H.-Q. Ye, High-resolution TEM observations of isolated rhombohedral crystallites in graphite blocks, Carbon N Y 50 (2012) 2347–2374.

[12] V.S. Neverov, XaNSoNS: gPU-accelerated simulator of diffraction patterns of nanoparticles, SoftwareX 6 (2017) 63–68.

[13] P. Puech, M. Monthoux, A new insight on the understanding of carbonisation and graphitisation mechanisms, Indian J. Eng. Mater. Sci. 27 (2020) 1095–1099.

[14] P. Scherrer, Bestimmung dergröße und der inneren struktur von kolloidteilchen mittels röntgenstrahlen, Nachrichten von der Königliche Gesellschaft der Wissenschaften zu Göttingen, Mathematisch-Physikalische Klasse (1918) 98–100.

[15] B.D. Cullity, Elements of X-ray Diffraction, 2nd edition, Addison-Wesley, 1978.

[16] H. Li, C. Yang, F. Liu, Novel method for determining stacking disorder degree in hexagonal graphite by X-ray diffraction, Sci. China Ser. B: Chem. 52 (2009) 174–180.

[17] J.B. Aladekomo, R.H. Bragg, Structural transformations induced in graphite by grinding: analysis of 002 X-ray diffraction line profiles, Carbon N Y 28 (1990) 897–906.

[18] J.W. Martin, J.L. Fogg, K.J. Putman, G. Francas, E.P. Turner, N.A. Marks, I. Suarez-Martinez, Graphite forms via annihilation of screw dislocations. arXiv: 2206.09105v (2022).