Digital

Discovery
This article is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported Licence.

View Article Online

PAPER View Journal | View Issue

A scientific machine learning framework to

Cite this: Digital Discovery, 2023, 2,
understand flash graphene synthesis†
1209
Open Access Article. Published on 18 July 2023. Downloaded on 7/2/2025 9:29:31 AM.

Kianoosh Sattari, a Lucas Eddy,bc Jacob L. Beckham,b Kevin M. Wyss,b

Richard Byfield,a Long Qian,b James M. Tour *bde and Jian Lin *a

Flash Joule heating (FJH) is a far-from-equilibrium (FFE) processing method for converting low-value
carbon-based materials to flash graphene (FG). Despite its promises in scalability and performance,
attempts to explore the reaction mechanism have been limited due to the complexities involved in the
FFE process. Data-driven machine learning (ML) models effectively account for the complexities, but the
model training requires a considerable amount of experimental data. To tackle this challenge, we
constructed a scientific ML (SML) framework trained by using both direct processing variables and
indirect, physics-informed variables to predict the FG yield. The indirect variables include current-derived
features (final current, maximum current, and charge density) predicted from the proxy ML models and
reaction temperatures simulated from multi-physics modeling. With the combined indirect features, the
final ML model achieves an average R2 score of 0.81 ± 0.05 and an average RMSE of 12.1% ± 2.0% in
predicting the FG yield, which is significantly higher than the model trained without them (R2 of 0.73 ±
Received 1st April 2023
Accepted 18th July 2023
0.05 and an RMSE of 14.3% ± 2.0%). Feature importance analysis validates the key roles of these indirect
features in determining the reaction outcome. These results illustrate the promise of this SML to
DOI: 10.1039/d3dd00055a
elucidate FFE material synthesis outcomes, thus paving a new avenue to processing other datasets from
rsc.li/digitaldiscovery the materials systems involving the same or different FFE processes.

thermodynamically stable sp2-hybridized graphene sheets.2

1 Introduction Because the whole process is nished in a sub-second scale, the
Despite the vast applications of graphene, scalable synthesis of generated graphene sheets form a metastable state, namely
graphene remains a tremendous challenge. Among the reported turbostratic graphene, which was termed as ash graphene
various types of processing methods,1,2 ash Joule heating (FJH) (FG).2 Such FG remains highly anisotropic in interlayer
was introduced in 2020 to synthesize gram-scale graphene from arrangements.3 This feature makes it highly dispersible in
different carbon feedstocks,3 such as carbon black (CB), solvents and a superior additive for high-performance
metallurgical coke (MC), and waste plastics.4,5 FJH is an elec- composites.3,6
trothermal process in which Joule heating, driven by capacitors The scalability of the FJH makes it a promising method for
with very high discharge rates, affords gross morphological synthesizing the FG, but many unknowns remain in this far-
changes.2 The generated high temperature (>3000 K) breaks the from-equilibrium (FFE) process,7 making it difficult to estab-
chemical bonds and reorganizes the carbon atoms into lish a processing–property relationship.8,9 Recently emerged
data-driven modeling may provide an alternative solution. In
a
the past several years, some models have been demonstrated to
Department of Mechanical and Aerospace Engineering, University of Missouri,
be powerful for tackling a variety of challenges including
Columbia, Missouri 65211, USA. E-mail: [email protected]; [email protected]
b
Department of Chemistry, Rice University, 6100 Main Street, Houston, Texas 77005,
guiding materials synthesis.10–14 Furthermore, we recently con-
USA structed pure data-driven models to discover the parameters
c
Applied Physics Program and Smalley-Curl Institute, Rice University, 6100 Main that controlled the FG yield.15 However, despite reaching an
Street, Houston, Texas 77005, USA impressive accuracy in predicting the FG yield, the model
d
Department of Materials Science and NanoEngineering, Rice University, 6100 Main performance depended on the current parameters measured
Street, Houston, Texas 77005, USA
e
from the reactions. These intermediate parameters were
Department of Computer Science and Engineering, NanoCarbon Center, Welch
therefore unavailable as input parameters for prediction if the
Institute for Advanced Materials, Rice University, 6100 Main Street, Houston, Texas
77005, USA experiments had not yet been performed. As a result, one
† Electronic supplementary information (ESI) available: Features distribution, cannot apply such models to accurately predict the reaction
simulated temperature of reactions inside the quartz tube, ML models' outcome from a new set of direct input parameters such as
hyperparameters, starting materials' physical properties (PDF). See DOI: voltage, pulse duration, and capacitance prior to
https://doi.org/10.1039/d3dd00055a

© 2023 The Author(s). Published by the Royal Society of Chemistry Digital Discovery, 2023, 2, 1209–1218 | 1209
 View Article Online

Digital Discovery Paper

experimentation, which makes them impractical for real important indirect features such as the current parameters of
applications. Thus, developing an ML framework that only uses nal current, maximum current, and charge density were pre-
the direct, controllable experimental parameters to accurately dicted from the proxy ML models. We hypothesize that these
predict reaction outcomes of FJH remains a challenge. current parameters are correlated with the direct experimental
This article is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported Licence.

Normally, a data-driven ML model is a “black-box”, lacking parameters and physical properties of the starting materials. To
interpretability in mapping the relationship of the input and validate this hypothesis, three proxy ML models were trained on
output. Moreover, model training requires considerable these direct parameters to predict those intermediate parameters
amount of data, a crucial aspect that has been a bottleneck for for a new experiment. In this way, the nal ML model does not
many materials processing methods such as FJH for FG rely on any intermediate information to predict the reaction
synthesis.16,17 In contrast, physics-based models can learn the outcome if given a new set of direct experimental parameters.
relationships of the input and output space. Although these Thus, the resulting SML framework is generalizable and needs
Open Access Article. Published on 18 July 2023. Downloaded on 7/2/2025 9:29:31 AM.

models are highly interpretable, they are oen difficult to be only limited training samples.21
constructed from complex systems due to a lack of information This SML framework has three advantages over our previ-
about the behavior of the system. Thus, the approximations are ously reported ML model.15 First, the models are able to make
needed to construct physics-based modeling while they can predictions about the reaction outcome without using any
result in inherent model bias. Therefore, hybrid models that intermediate parameters. This facilitates the use of our
combine data-driven and physics-based modeling can be prediction model in a model-based optimization algorithm to
benecial in successful model training with limited experi- optimize the FG yield in just a few iterations. Second, the
mental data while offering high explainability.18–20 These physics-informed descriptors bring additional information to
models can be constructed by modifying the cost functions the model, making the black-box ML models more generaliz-
within data-driven ML models. This modication can adjust the able and accurate in addition to improving the model inter-
model to obey the outputs of the physics-based models. Daw pretability. Third, a general methodology of using separate ML
et al. designed a physics-guided neural networks (PGNNs) models to predict unknown, intermediate reaction parameters
framework that leveraged the output of the physics-based from known direct ones is proposed to solve the challenge of
model and observational features by modifying the loss func- lacking enough input features, particularly related to experi-
tion of the neural network.21 Raissi et al. introduced physics- ments. Thus, such an approach can be readily applied to other
informed neural networks (PINNs) that obeyed physics laws materials processed by the same or different methods. For
described by partial differential equations.22 The additional instance, our developed framework could be employed to
information gained from the physical laws can train the predict the yield of the precious metals recovered from elec-
networks with much less data than needed in pure ML models, tronic wastes as well as the removing efficiency of hazardous
thus broadening the applications where data generation is heavy metals from them by using the same FJH process.28
costly.17 Rabczuk's group constructed deep neural networks Furthermore, the framework holds a potential for guiding
(DNNs) to approximate partial differential equations for structure- and phase-controlled synthesis of inorganic nano-
predictions of materials properties.23,24 They included the crystals through ultrafast FJH.29,30 It can be done by some
system energy in the loss function of the models. However, in exploratory studies to identify the specic important features
the FJH process, they are not practical since there are no dened and collect more experimental data from the new reactions.
physical rules that can well describe the FFE reactions. Another
method of including physics laws into the ML models is to
extract physics-informed features from the experiments or 2 Results and discussion
theory, which are used as the model input to boost the predic- This work used a dataset consisting of 173 separate FJH reac-
tion accuracy.20,25 Sun et al. synergized the indirect physics- tions reported in our previous work.15 The starting materials
informed descriptors with other direct variables in the ML were carbon black (CB), metallurgical coke (MC), plastic waste-
framework to develop materials with superior properties.26 To derived pyrolysis ash (PA), and waste tire-based carbon black
develop thermo-responsive materials, Huang et al. developed (TCB). The structures of the nal products were assessed by
a framework where ML models were informed with physico- wide-area Raman mapping. We applied custom-written scripts
chemical descriptors derived from quantum chemistry calcu- to analyze the collected Raman spectra, which were used to
lation.27 Such physics-based descriptors can serve as the estimate the FG yield. The high-throughput mapping assay
indirect input features to introduce partial physical information guarantees the comprehensive representation of the overall
to the ML framework. property of each sample. In the following sections, we rst
To better understand the FJH process for FG synthesis, herein, analyze the dataset and explain how to quantify the FG yield. We
we demonstrate a scientic machine learning (SML) framework then elaborate the SML framework. Lastly, we present the
that is trained with both direct experimental parameters and model performance in predicting the FG yield.
indirect physics-informed ones. The goal is to predict the yield of
FG. To estimate the reaction temperature from the direct exper-
imental parameters (such as pulse time, voltage, capacitance, and 2.1 Analysis of input and output data
physical information of the input materials), we performed an Raman spectroscopy has been considered a powerful technique
electrical–thermal multi-physics simulation by COMSOL. Other for characterization of carbon structures.31,32 Fig. 1a shows

1210 | Digital Discovery, 2023, 2, 1209–1218 © 2023 The Author(s). Published by the Royal Society of Chemistry
 View Article Online

Paper Digital Discovery

Raman spectra of amorphous carbon and synthesized FG. The of all the 173 reactions. Specically, Fig. 1b shows the distri-
spectrum of amorphous carbon shows two main peaks: D-band bution of average I2D/IG with a mean of 0.66 and a standard
at ∼1350 cm−1 and G-band at ∼1600 cm−1. The Raman spec- deviation of 0.17. Fig. 1c represents a histogram of average ID/IG
trum of FG has a G-peak at ∼1580 cm−1 and a 2D band at with a mean of 0.54 and a standard deviation of 0.14. Fig. 1d
This article is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported Licence.

∼2700 cm−1. The existence of this 2D band suggests formation represents the average FWHM of the 2D band with a mean of
of a graphitic lattice.32 This resonance-enhanced single- 43.88 cm−1 and a standard deviation of 11.55 cm−1. Finally,
Lorentzian 2D band has a narrow full-width at half-maximum Fig. 1e shows a histogram of the FG yield with a mean of 54%
(FWHM) of ∼16 cm−1. The I2D/IG peak intensity ratio reaches and a standard deviation of 27%. Fig. S1† shows the yield
up to 17. Both of them suggest good FG crystallinity.33 From distribution of the FG synthesized from the four starting
each sample, we collected 100 Raman spectra, which was then materials.
averaged to mitigate the variance in the collected individual Fig. 1f and g represent high correlation of I2D/IG with the FG
Open Access Article. Published on 18 July 2023. Downloaded on 7/2/2025 9:29:31 AM.

spectrum. Then, the FG yield can be calculated from these yield, showing a Pearson's r value of 0.73. Fig. 1f shows little
averaged spectra.15 Fig. 1b–e represent the histograms and dependence of the FG yield on ID/IG, while the value of FWHM
statistics distribution of the collected samples for each reaction can well distinguish the samples with a high FG yield (Fig. 1g).

Fig. 1 (a) Raman spectra of flash graphene (FG) synthesized from carbon black and amorphous carbon. (b–e) Statistical distribution of I2D/IG (b),
ID/IG (c), FWHM of the 2D band (d), and FG yield (e). Distribution of I2D/IG versus FG yield in correlation with (f) ID/IG and (g) FWHM. (h) Distribution
of FG yield synthesized from four starting materials. (i) Statistical comparison on the mean FG yield from four starting materials. (***), (**), and (*)
show significant differences at 0.001, 0.01, and 0.05 levels, respectively. (j) t-SNE plots of features in correlation with the four starting materials.
The features include resistance drop, voltage drop, maximum current, charge density, I2D/IG, ID/IG, FWHM, and reaction yield.

© 2023 The Author(s). Published by the Royal Society of Chemistry Digital Discovery, 2023, 2, 1209–1218 | 1211
 View Article Online

Digital Discovery Paper

Most samples have average FWHM values of >40 cm−1 and I2D/ temperature evolution is simulated from the direct parameters
IG > 0.75. We also analyzed the FG yield from different starting by multi-physics simulation to obtain the maximum tempera-
materials. As illustrated in Fig. 1h, the highest FG yield of 72% ture (TSim.). Thus, compared to our previous model that predicts
and the lowest yield of 37% were obtained for CB and MC, the FG yield,15 more physics-informed input features are used to
This article is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported Licence.

respectively. Fig. 1i shows the statistical comparison of the FG improve the prediction accuracy and generalizability of the nal
yield obtained from the four starting materials. Except for MC model. In the following sections, we will elaborate the proxy
versus TCB, all other two-way comparisons show signicant models, the multi-physics simulation, and the overall archi-
differences at a set 0.05 signicance level. tecture of the nal prediction model.
We hypothesized that the measured parameters including 2.2.1 Proxy models for predicting current parameters. The
resistant drop, voltage drop, nal current, maximum current, current–time curves are measured from the FJH process. Three
charge density, I2D/IG, ID/IG, FWHM, and reaction yield would parameters of IMax, IF/IMax, and CDIT can be extracted from these
Open Access Article. Published on 18 July 2023. Downloaded on 7/2/2025 9:29:31 AM.

depend on the starting material. To test the hypothesis, we curves (Fig. 3a). The distributions of these current parameters
applied t-distributed stochastic neighbor embedding (t-SNE),34 depending on the starting materials were analyzed (Fig. S2†).
a non-linear dimension reduction method, to project all of Signicantly higher IMax values could be realized in the reaction
them in 2D space (Fig. 1j). This analysis shows that those ob- outcomes using MC as the staring material than those in the
tained from MC and CB are clustered and separated from the reactions using other starting materials (Fig. S2a†). But the
others, which indicates that there do exist combination of the higher IMax values do not simply lead to a higher FG yield for the
parameters for achieving the highest FG yield in CB (Fig. 1h). MC samples, as shown in Fig. 1h. Fig. S3a† shows plots of the
The signicant difference in the FG yield from different staring FG yield vs. IMax grouped by the starting materials. Corre-
materials indicates that besides the one-hot encoded material spondingly, Pearson's r values between IMax and the FG yield for
type, inclusion of physical information about the starting CB, PA, and TCB are 0.41, 0.62, and 0.66, respectively, indicating
materials like particle size (MPS), resistance (MR), surface area that they have high correlations, while the correlation of IMax
(MSA), and percentage of sp2 carbon (Msp2) in the input features and the FG yield is not signicant for MC (Fig. S3†). The positive
would greatly increase model accuracy. All these physical correlations between IMax and FG yield for CB, PA, and TCB
properties of the starting materials are tabulated in Table S1.† show that the IMax should pass a threshold value of 1000 (A g−1)
for these samples to reach a higher FG yield.
To train the proxy models that predict these three current
2.2 Model construction and performance parameters, the direct reaction parameters, including the
In our previous study,15 we conducted an exploratory analysis properties of starting materials and FJH parameters, serve as
aiming to rank importance of the input features in determining the inputs of the models which were trained by a ve-fold cross-
the FG yield. However, the trained model in that study relies on validation approach. To test the models, 20% of the total
intermediate experimental variables, thus cannot be applied to samples were used as the never-seen samples. The optimized
a new set of hypothesized experimental parameters if the hyperparameters for these three proxy XGBoost models are lis-
experiment has not been implemented. That is because the ted in Table S2.† It is worth mentioning that the inputs to the
intermediate variables must be rst extracted from the experi- proxy models can be hypothesized for predicting reaction
mentally obtained current–time curves. To address this limita- outcome of a new experiment without performing it. As a result,
tion, herein, we have utilized only the direct experimental the trained models can be used to predict the three current
variables combined with physics-informed features—which parameters for a new reaction. Fig. 3b–d shows comparison of
were derived from the direct experimental variables—as the the predicted three current parameters from the proxy models
input features to train a new SML framework. This framework, versus their true values, from which their Pearson's r values can
shown in Fig. 2, no longer relies on the intermediate features as be calculated to evaluate performance of the proxy models.
the inputs. Instead, it uses only direct experimental variables as Pearson's r values of 0.80, 0.78, and 0.77 were obtained for IMax,
input to predict the FG yield for a new set of hypothesized IF/IMax, and CDIT, respectively. The high correlations between
experimental parameters without conducting the prior experi- the predicted and the true values show that the proxy models
ment. They include direct reaction parameters such as the can predict the output IMax, IF/IMax, and CDIT from the direct
properties of starting materials including particle size (MPS), parameters so that no prior-measurement on the current–time
resistance (MR), and the percent of sp2 content (Msp2) and FJH curves for a hypothesized FJH experiment would be needed.
controllable parameters including charge density released from In our recent work,15 we conducted an exploratory analysis to
capacitance (CD0), heat (H), pulse time (t), atmosphere type examine the inuence of various input features on the predic-
(Atm), and pretreatment voltage (VPre). Using these direct tion of FG yield. Through the importance analysis, we deter-
parameters, three proxy models based on XGBoost were trained mined that the normalized maximum current, the ratio of nal
to predict three intermediate parameters of maximum current to maximum current, and the charge density exhibited were
normalized by mass (IMax), ratio of nal current to maximum most important. Additionally, we experimented additional
current (IF/IMax), and charge density (CDIT, total charge inte- features obtained from current–time (I–T) curves, such as the
grated from the current–time curve and then normalized to mean current in the I–T curve, the number of local maxima/
mass). In this way, measurement of the current–time curves minima in the I–T curve, and the current-derived properties
from a hypothesized experiment is no longer needed. Third, the before normalization. However, these attempted features did

1212 | Digital Discovery, 2023, 2, 1209–1218 © 2023 The Author(s). Published by the Royal Society of Chemistry
 View Article Online

Paper Digital Discovery

This article is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported Licence.
Open Access Article. Published on 18 July 2023. Downloaded on 7/2/2025 9:29:31 AM.

Fig. 2 Schematic and data flow of the proposed SML framework, where the temperature simulated by the multi-physics simulation, predicted
current parameters, precursor information, and direct FJH parameters are used as the input of the final ML model.

represented as TSim.. In the simulation, the direct input mate-

rials and reaction parameters were used. As shown in Fig. S4a
and b,† the FJH quartz tube was simulated as a cylinder with
a diameter of 8 mm and a length of ∼20 mm. TSim. over the
pulse time for all the 173 reactions are shown in Fig. S4c.† It
shows that the relationship between the temperature and pulse
time is not a linear one. There are reactions realizing a higher
temperature in a smaller pulse time.
In future, it is possible to apply multi-scale simulation
techniques such as the one introduced by Talebi et al.35 We
hypothesize that collection of the experimental temperature
during the reaction for validation, these techniques would
render better understanding on the FJH process.
2.2.3 Performance of the nal model. The predicted
current parameters and TSim. were combined with the direct FJH
parameters and precursor information to serve as inputs of six
different regression models including linear regression (LR),
Fig. 3 (a) A represented current–time plot and the current parameters
multilayer perceptron (MLP), Bayesian regression (BR), decision
derived from it. Distributions of predicted and true (b) IMax; (c) IF/IMax;
and (d) CDIT values. Their corresponding Pearson's r values are shown tree (DT), random forest (RF), and eXtreme Gradient Boosting
in the figures. (XGBoost). By using a 5-fold cross-validation method for
training and testing, the optimized hyperparameters for these
models are listed in Table S3.† Fig. 4a and b show the coefficient
not enhance the prediction accuracy and had marginal impor- of determination (R2) and root mean squared error (RMSE) for
tance (close to 0) in predicting the FG yield. all six tested models in predicting the FG yield. Among them,
2.2.2 Simulation of reaction temperature as a physics- the XGBoost model reached the highest average R2 score of 0.81
informed input feature. In an FJH process, the electrical with a standard deviation of 0.05 and the lowest average RMSE
energy is rapidly discharged from capacitors, leading to a time- of 12.1% with a standard deviation of 2.0% on the testing
dependent, spatially distributed temperature prole. While samples for 5 different train-test splits. By trying the 5-split
temperature is an important parameter that controls the FG method, one can identify if the model's performance is
yield, we hypothesize that using it as an input feature would consistent across different subsets of the data or if it is sensitive
improve the predictive accuracy of the model. Deng et al. re- to specic variations. Increasing the number of training
ported the effects of direct reaction parameters like the mass of samples can decrease the variance of R2 over different splits.
the starting materials, physical properties of starting materials, The result shows that the designed framework obtains a slightly
pulse time, pulse voltage, pre-treatment voltage, and the higher accuracy compared to our previous model.15 The differ-
maximum temperature achieved in the FJH process.28 To test ence is that unlike the previous model relying on calculation of
the hypothesis, the electrical–thermal multi-physics package in intermediate variables the current model uses only the direct
COMSOL was applied to simulate the temperature evolving over experimental variables as the input. Taking a XGBoost model
the pulse duration of each reaction. The maximum temperature trained from one of the 5 different splits for example, compar-
of the reaction was then used as an input descriptor, ison of the predicted FG yields versus the true values was shown

© 2023 The Author(s). Published by the Royal Society of Chemistry Digital Discovery, 2023, 2, 1209–1218 | 1213
 View Article Online

Digital Discovery Paper

To test the signicance of including the physics-informed

features as the input to the model, we trained a separate
XGBoost model without using them as the input. As shown in
Fig. 4e, if the TSim. is excluded, the R2 score is reduced to 0.79
This article is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported Licence.

and RMSE is increased to 13.7% for the same testing dataset. If

both the simulated temperature and the predicted current
parameters are excluded, the R2 score is greatly decreased to
0.74 and RMSE is increased to 15.1% (Fig. 4f). This results
because the current parameters may reect the change of the
starting materials' resistance and the contact resistance
between the starting materials and the electrode over the pulse
Open Access Article. Published on 18 July 2023. Downloaded on 7/2/2025 9:29:31 AM.

time. The temperature is a key parameter that determines the

reaction outcome. Consequently, these physics-informed
descriptors can offer complementary information to the
model with increased the prediction accuracy.

2.3 Model interpretation

Ranking importance of the input features to the well-trained
model in predicting the FG yield would offer additional infor-
mation about the reaction. The selected features included the
CD0, MPS, MR, Msp2, predicted IMax, predicted IF/IMax, predicted
CDIT, TSim., t, VPre, Atm, and H. A Pearson's correlation map
between these quantitative features is shown in Fig. 5a. Low
Pearson's r values between any two features indicate that they
Fig. 4 Performance of the ML models in predicting the FG yield. (a) R2
are quite independent features for the model to afford accurate
scores and (b) RMSE of the predicted FG yield by the six ML models
when using five different train-test splitting ways. The error bars prediction. For instance, the correlation of the chosen physical
represent the standard deviations from these five testing ways. (c) Plot properties of the starting materials is low, indicating that they
of predicted FG yields by the XGBoost model vs. their true values from offer complementary information of the materials properties
different starting materials. (d) Relative error distribution of the pre- when serving as the input features. In contrast, the surface area
dicted FG yields shown in (c). Plot of the predicted FG yields by the
has a high Pearson's r value of 0.9 with the particle size, thus we
XGBoost model vs. their true values after excluding (e) TSim. and (f) both
TSim. and predicted current parameters from the direct input excluded it from the nal input features. Fig. 5b shows the
parameters. ranking of the features. CD0 and TSim. were ranked the top 2
important features in determining the FG yield, which explains
why they play a critical role in the model accuracy (Fig. 4). Other
features such as the predicted current parameters also have
in Fig. 4c from which an R2 score of 0.84 and RMSE of 11.8% a signicant importance in the nal prediction. In previous
were calculated. As a comparison purposes, we considered works,15,36,37 voltage and CD0 were reported to have effects on the
a base model that predicts the average value of all testing transformation rate. Fig. 5c shows that the FJH reactions with
samples for all the samples. The RMSE for such a naı̈ve model low CD0 values have a lower FG yield. In contrast, the ones
was 29.6% that is signicantly higher than that of XGBoost leading to a high FG yield have high CD0 values. This observa-
predictions. Samples ashed with CB as the starting material tion agrees well the results shown in these works. In addition, it
possessed the highest FG yields, while MC-derived FG had the is found that there is a CD0 threshold value of 100 (C g−1) for
lowest FG yield. Fig. 4d shows the relative error (RE) distribution achieving an FG yield of >50%. This observation agrees well
of the predicted FG yields compared with the true values. It with other FFE processes. For instance, laser-induced synthesis
shows that 71% of the reactions have the predicted yields of of graphene from polymers was only initiated when a laser ux
#10% error of the true values, and only ∼11% of the reactions reaches a threshold value.38 Fig. 5d shows the importance of
show the predicted FG yields with an error of >20%. We further TSim. in predicting the FG yield. It shows that when TSim. exceeds
examined the distribution of the residuals, a difference of the a threshold value as indicated in green yellow, and red colors,
predicted and the true values. The residuals show a biased the FG yield is signicantly higher than those with low TSim.. A
toward negative values for samples with the high FG yields, as decision tree extracted from the XGBoost model supports the
shown in Fig. S5.† This indicates that the model usually predicts hypothesis that high TSim. and CD0 are critical in model accu-
a lower FG yield value for the reactions resulting in a higher FG racy for predicting the FG yield (Fig. S6†). Fig. S7† compares CD0
yield value, while for the training samples with an average FG with C, V0, and m in correlating with the FG yield. It shows that
yield of 54%, the predictions for unsure testing samples are correlation of FG yield with CD0 is higher than that with C, V0,
biased toward the average value. and m, which validates the importance of CD0 in the accurate
prediction of the FG yield.

1214 | Digital Discovery, 2023, 2, 1209–1218 © 2023 The Author(s). Published by the Royal Society of Chemistry
 View Article Online

Paper Digital Discovery

4 Methods and experimental section

4.1 Materials
Four carbon feedstocks were used as the starting materials.
This article is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported Licence.

They are carbon black (Cabot BP2000), metallurgical coke

(SunCoke Energy Inc., 70–100 mesh size, 150–210 mm grain
size), pyrolysis ash (Shangqiu Zhongming EcoFriendly Equip-
ment Co.), and pyrolyzed rubber tire-derived carbon black
(Ergon Asphalt and Emulsion Co.). We ground the materials
using a mortar and pestle before and aer FJH.
Open Access Article. Published on 18 July 2023. Downloaded on 7/2/2025 9:29:31 AM.

4.2 FJH process

A custom FJH apparatus was used for all the 173 experiments.
Precursor powders with a mass between 100 and 400 mg were
sandwiched between two graphite electrodes and compressed
inside a quartz tube with an inner diameter of 8 mm. Then,
a series circuit with eight 6 mF capacitors (Mouser #80-
PEH200YX460BQU2), two 5.6 mF capacitors (80-
ALS70A562QH500), and nine 18 mF capacitors (Mouser #80-
ALS70A183QS400) were used. Arrangement of capacitors was set
to reach the peak capacitance values employed in each ash
reaction. To charge the capacitors, the voltage was supplied by
a DC source consisting of an AC wall outlet fed through an AC–
DC converter. FJH reactions were performed inside a desiccator
lled with argon, air, or light vacuum (10 mm Hg) that was used
as a categorical descriptor for atmosphere type (Atm) among
direct input features. Aer applying the initial voltage, the nal
voltage was recorded aer each reaction. A voltage drop was
then calculated by subtracting the nal voltage from the initial
one. Resistance of the samples were measured before and aer
Fig. 5 Analysis of the input features to the final XGBoost model. (a) each reaction to monitor electrical contact between the elec-
Quantitative correlation map of the input features. (b) Feature trodes and the samples. Pulse time was modulated by insulated
importance of the input features. Predicted FG yields versus the true gate bipolar transistors (IGBTs) using programmable
values when correlated with (c) CD0 and (d) TSim.. In (d) TSim. is in a log millisecond-level delay time. It was connected to a Hall effect
scale.
sensor through an inductor and controlled via custom LabVIEW
scripts. The Hall effect sensor was employed to collect current–
time curves. A custom-written Python script was applied on the
3 Conclusion current–time curves to extract current parameters for the proxy
This study demonstrates an SML framework that bridges a gap model training.
between the input processing parameters with the predicted FG
yield. Herein, a systematic method of using proxy ML models 4.3 Material characterization and analysis
and multi-physics simulation for extracting physics-informed Wide-area Raman spectral mapping was selected as the primary
descriptors, including current-derived properties and simu- method for characterizing the FG property. To mitigate poten-
lated temperature, has been developed. These additional input tial biases caused by spatial variations, mapping was conducted
features prove to play a critical role in improving the prediction on a 1 mm2 region to obtain 64 spectra per mapping. The
accuracy of the nal ML model. Feature importance analysis spectra were acquired using a Renishaw inVia Raman micro-
further validates this conclusion. Besides the TSim. and CD0, the scope equipped with a 50× lens and a 5 mW 532 nm Nd:YAG
selected physical properties of the starting materials are also laser for excitation. Laser focus was maintained across the area
important features. Explainability of the model by the quanti- using Renishaw Wire 5.5 LiveTrack soware. The resulting
tative analysis offers a glimpse on the reaction mechanism spectra were analyzed using the written Python scripts. Before
about the FJH. In summary, development of this SML frame- analysis, each spectrum underwent baseline correction using
work offers a methodology of predicting the outcome of new a RamPy polynomial t and then was smoothened via
experiments, thus saving the cost and time because of per- a Savitzky–Golay lter. A FJ product was categorized as “gra-
forming unnecessary experiments, which would speed up the phene” was based on the following characteristics of its spec-
FG synthesis. Finally, the methodology can be readily applied to trum: I2D/IG > 0.3, 15 cm−1 < FWHM2D < 70 cm−1, and a signal-
other material systems processed by other processing methods. to-noise ratio (SNR2D) of > 8. Others were categorized as

© 2023 The Author(s). Published by the Royal Society of Chemistry Digital Discovery, 2023, 2, 1209–1218 | 1215
 View Article Online

Digital Discovery Paper

“amorphous carbon”. Spectra lacking a sufficiently strong G V0, while CDIT conveys information about the voltage drop
band (SNRG < 8) were excluded from the analysis as they indi- during the FJH process.
cated poor laser beam focusing. The FG yield was calculated for
each sample by dividing the total number of Raman spectra that 4.7 Evaluation metrics
This article is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported Licence.

were categorized for graphene by the total number of Raman

The coefficient of determination (R2) is used to evaluate the
spectra obtained from that sample. This afforded an approxi-
prediction accuracy of a model as shown in eqn (4). The Pearson
mate numerical measure of the sample's bulk crystallinity.
correlation coefficient (r) dened in eqn (5), on the other hand,
measures how the predicted values catch the trend compared to
4.4 Data inclusion the true values.
At the spectra-level, we included all spectra identied as having P
N
ðyi  y^i Þ2
a G peak with an SNR of >8 (in the range of 1500–1700 cm−1).
Open Access Article. Published on 18 July 2023. Downloaded on 7/2/2025 9:29:31 AM.

2 i¼1
R ¼1 (4)
Spectra not containing a G peak were attributed to poor laser P
N
ðyi  yÞ2
focusing and excluded. i¼1

P
N  
4.5 Training of machine learning models ðyi  yÞ  y^i  y^
i¼1
r ¼ sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
ffi (5)
Six different ML models (LR, MLP-R, BR, DT-R, RF-R, and XGB- P N  2
2P
N
R) were trained to predict FG yield. The Scikit-Learn package ðyi  yÞ y^i  y^
i¼1 i¼1
from Python was used for constructing all the models. We kept
20% of the dataset unseen for testing. Cross-validation was where y is the true values, ŷ is the predicted values, y is the mean
applied to optimize the hyperparameters. To test the accuracy of value, and N is the number of samples in both. In eqn (5), ^y is
the model for different testing samples, we tried 5 different the average of all predicted ŷ.
train/test splits. The results were reported as metrics' mean ± Other evaluation metrics including residuals (R), relative
standard deviation. error (RE) and root mean squared error (RMSE) are dened in
eqn (6)–(8), respectively.
4.6 Feature engineering
R = ŷ − y (6)
Twelve selected features (as tabulated in Table S4†) included the
jy  y^j
charge density (CD0) released from the capacitors, starting RE ¼  100% (7)
y
materials' type (M), particle size (MPS), resistance (MR), surface
area (MSA), and percent of sp2 content (Msp2), predicted vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
u
normalized maximum current (IMax), predicted ratios of nal u1 X N
RMSE ¼ t ðyi  y^i Þ2 (8)
current to the maximum current (IF/IMax), predicted charge N i¼1
density that is dened as area under the current–time curve
normalized by mass (CDIT), simulated temperature (TSim.), where y is the true values, ŷ is the predicted values, and N is the
pulse time (t), pre-treatment voltage (VPre), atmosphere type number of samples.
(Atm), and nominal heat (H) were used as the input features to
the nal ML models. 4.8 FEA simulation on temperature
CD0, CDIT, and H are dened in eqn (1)–(3), respectively. The electrical–thermal multi-physics package in COMSOL
V0  C Multiphysics version 6.0 licensed under University of Missouri
CD0 ¼ (1)
m was applied to simulate the temperature evolving over the pulse
duration of each reaction. In the simulation, we added three
I t physics interfaces nodes—electric currents, heat transfer in
CDIT ¼ (2)
m solids, and the multiphysics node. The starting materials mass
and particle sizes as well as pulse time, voltage, and capacitance
V0 2
H¼ (3) of each reaction were used as the input to the simulation. Also,
MR  t
we considered 140, 130, 120, and 113 (S m−1) for the electrical
where V0 is the voltage, C is the capacitance of the capacitors, m conductivity and 0.4, 1.2, 2.2, and 2.7 (W m−1 K−1) for the
is the mass of the starting materials, MR is the initial resistance thermal conductivity of the starting materials CB, PA, MC, and
of the starting material, and t is the pulse time. M is one-hot TCB, respectively. The applied electrical and thermal conduc-
encoding for the types of the starting materials. It was only tivity values are in the range of reported experimental
used as input to the proxy models and not in the nal model. values.39–41 To set up the electrical boundary conditions, one
CDIT was calculated by trapezoidal integration of the current– side of the simulated cylinder was grounded (0 V) and the other
time curve collected by a Hall effect sensor. Even if CDIT and side was applied with the input voltage values shown in the
CD0 have the same units, they include different information dataset. To set up the heat boundary conditions, we used the
about the reaction. CD0 depends on the initial nominal voltage room temperature as the initial temperature of the system and

1216 | Digital Discovery, 2023, 2, 1209–1218 © 2023 The Author(s). Published by the Royal Society of Chemistry
 View Article Online

Paper Digital Discovery

applied convective heat ux around the cylinder surfaces. Aer plastic waste pyrolysis ash into ash graphene, Carbon,
nding the location with the maximum temperature in each 2021, 174, 430–438.
reaction, we used the nal simulated temperature (end of each 5 P. A. Advincula, V. Granja, K. M. Wyss, W. A. Algozeeb,
pulse time) of the location as the input to the SML as TSim.. W. Chen, J. L. Beckham, D. X. Luong, C. F. Higgs and
This article is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported Licence.

J. M. Tour, Waste plastic- and coke-derived ash graphene

Data availability as lubricant additives, Carbon, 2023, 203, 876–885.
6 Y. Wu, P. A. Advincula, O. Giraldo-Londoño, Y. Yu, Y. Xie,
Data and processing scripts for this paper, including the Z. Chen, G. Huang, J. M. Tour and J. Lin, Sustainable 3D
collected experiments are available at GitHub repository at Printing of Recyclable Biocomposite Empowered by Flash
https://github.com/linresearchgroup/SciML_FJH. Graphene, ACS Nano, 2022, 16, 17326–17335.
7 R. Raj, Joule heating during ash-sintering, J. Eur. Ceram.
Open Access Article. Published on 18 July 2023. Downloaded on 7/2/2025 9:29:31 AM.

Author contributions Soc., 2012, 32, 2293–2301.

8 D. Dai, Q. Liu, R. Hu, X. Wei, G. Ding, B. Xu, T. Xu, J. Zhang,
J. L. conceived the idea. K. S. designed the framework, con- Y. Xu and H. Zhang, Method construction of structure-
structed the ML models, performed the temperature simula- property relationships from data by machine learning
tion, and analyzed the data. K. S. wrote the rst manuscript assisted mining for materials design applications, Mater.
which was thoroughly revised by J. L. R. B. assisted K. S. in Des., 2020, 196, 109194.
model development and manuscript writing. J. L. B. provided 9 J. Yu, X. Yong, Z. Tang, B. Yang and S. Lu, Theoretical
discussion on the results. L. E. and K. M. W. performed exper- Understanding of Structure–Property Relationships in
iments for data collection. J. M. T. supervised L. E., J. L. B., and Luminescence of Carbon Dots, J. Phys. Chem. Lett., 2021,
K. M. W., in the experimental design, data collection as well as 12, 7671–7687.
revising the manuscript. All authors discussed and commented 10 K. Sattari, Y. Xie and J. Lin, Data-driven algorithms for
on the manuscript. inverse design of polymers, So Matter, 2021, 17, 7607–7622.
11 Y. Xie, K. Sattari, C. Zhang and J. Lin, Toward autonomous
Conflicts of interest laboratories: Convergence of articial intelligence and
experimental automation, Prog. Mater. Sci., 2023, 132,
Universal Matter Inc. has licensed the FG process from Rice 101043.
University. J. M. T. is a stockholder in that company, but not an 12 P. Raccuglia, K. C. Elbert, P. D. F. Adler, C. Falk,
officer, director, or employee. Conicts of interest are managed M. B. Wenny, A. Mollo, M. Zeller, S. A. Friedler, J. Schrier
through regular disclosure to and compliance with the Rice and A. J. Norquist, Machine-learning-assisted materials
University Office of Sponsored Programs and Research discovery using failed experiments, Nature, 2016, 533, 73–76.
Compliance. 13 Y. Xie, C. Zhang, X. Hu, C. Zhang, S. P. Kelley, J. L. Atwood
and J. Lin, Machine Learning Assisted Synthesis of Metal–
Acknowledgements Organic Nanocapsules, J. Am. Chem. Soc., 2020, 142, 1475–
1481.
J. L. and J. M. T. thank U.S. Army Corps of Engineers, ERDC 14 C. Wen, C. Wang, Y. Zhang, S. Antonov, D. Xue, T. Lookman
(grant number: W912HZ-21-2-0050) for the nancial support. and Y. Su, Modeling solid solution strengthening in high
This work was also partially funded by National Science Foun- entropy alloys using machine learning, Acta Mater., 2021,
dation (award numbers: 1825352 and 2154428), and the Air 212, 116917.
Force Office of Scientic Research (FA9550-22-1-0526). 15 J. L. Beckham, K. M. Wyss, Y. Xie, E. A. McHugh, J. T. Li,
P. A. Advincula, W. Chen, J. Lin and J. M. Tour, Machine
References Learning Guided Synthesis of Flash Graphene, Adv. Mater.,
2022, 34, 2106506.
1 Z. Sun, S. Fang and Y. H. Hu, 3D Graphene Materials: From 16 J. Hoffmann, Y. Bar-Sinai, L. M. Lee, J. Andrejevic, S. Mishra,
Understanding to Design and Synthesis Control, Chem. Rev., S. M. Rubinstein and C. H. Rycro, Machine learning in
2020, 120, 10336–10453. a data-limited regime: Augmenting experiments with
2 K. M. Wyss, D. X. Luong and J. M. Tour, Large-Scale synthetic data uncovers order in crumpled sheets, Sci. Adv.,
Syntheses of 2D Materials: Flash Joule Heating and Other 2019, 5, eaau6792.
Methods, Adv. Mater., 2022, 34, 2106970. 17 G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris,
3 D. X. Luong, K. V. Bets, W. A. Algozeeb, M. G. Stanford, S. Wang and L. Yang, Physics-informed machine learning,
C. Kittrell, W. Chen, R. V. Salvatierra, M. Ren, Nat. Rev. Phys., 2021, 3, 422–440.
E. A. McHugh, P. A. Advincula, Z. Wang, M. Bhatt, H. Guo, 18 B. Kapusuzoglu and S. Mahadevan, Physics-Informed and
V. Mancevski, R. Shahsavari, B. I. Yakobson and J. M. Tour, Hybrid Machine Learning in Additive Manufacturing:
Gram-scale bottom-up ash graphene synthesis, Nature, Application to Fused Filament Fabrication, JOM, 2020, 72,
2020, 577, 647–651. 4695–4705.
4 K. M. Wyss, J. L. Beckham, W. Chen, D. X. Luong, P. Hundi, 19 J. Willard, X. Jia, S. Xu, M. Steinbach and V. Kumar,
S. Raghuraman, R. Shahsavari and J. M. Tour, Converting Integrating scientic knowledge with machine learning for

© 2023 The Author(s). Published by the Royal Society of Chemistry Digital Discovery, 2023, 2, 1209–1218 | 1217
 View Article Online

Digital Discovery Paper

engineering and environmental systems, ACM Comput. High-surface-area corundum nanoparticles by resistive
Surv., 2022, 55, 1–37. hotspot-induced phase transformation, Nat. Commun.,
20 M. Arias Chao, C. Kulkarni, K. Goebel and O. Fink, Fusing 2022, 13, 5027.
physics-based and deep learning models for prognostics, 30 B. Deng, Z. Wang, W. Chen, J. T. Li, D. X. Luong, R. A. Carter,
This article is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported Licence.

Reliab. Eng. Syst. Saf., 2022, 217, 107961. G. Gao, B. I. Yakobson, Y. Zhao and J. M. Tour, Phase
21 A. Daw, A. Karpatne, W. Watkins, J. Read and V. Kumar, controlled synthesis of transition metal carbide
arXiv, 2021, preprint, arXiv:1710.11431v3, DOI: 10.48550/ nanocrystals by ultrafast ash Joule heating, Nat.
arXiv.1710.11431. Commun., 2022, 13, 262.
22 M. Raissi, P. Perdikaris and G. E. Karniadakis, Physics- 31 A. C. Ferrari, Raman spectroscopy of graphene and graphite:
informed neural networks: A deep learning framework for Disorder, electron–phonon coupling, doping and
solving forward and inverse problems involving nonlinear nonadiabatic effects, Solid State Commun., 2007, 143, 47–57.
Open Access Article. Published on 18 July 2023. Downloaded on 7/2/2025 9:29:31 AM.

partial differential equations, J. Comput. Phys., 2019, 378, 32 A. C. Ferrari and D. M. Basko, Raman spectroscopy as
686–707. a versatile tool for studying the properties of graphene,
23 C. Anitescu, E. Atroshchenko, N. Alajlan and T. Rabczuk, Nat. Nanotechnol., 2013, 8, 235–246.
Articial Neural Network Methods for the Solution of 33 J. A. Garlow, L. K. Barrett, L. Wu, K. Kisslinger, Y. Zhu and
Second Order Boundary Value Problems, Comput. Mater. J. F. Pulecio, Large-Area Growth of Turbostratic Graphene
Continua, 2019, 59, 345–359. on Ni(111) via Physical Vapor Deposition, Sci. Rep., 2016, 6,
24 E. Samaniego, C. Anitescu, S. Goswami, V. M. Nguyen- 19804.
Thanh, H. Guo, K. Hamdia, X. Zhuang and T. Rabczuk, An 34 L. Van der Maaten and G. Hinton, Visualizing data using t-
energy approach to the solution of partial differential SNE, J. Mach. Learn. Res., 2008, 9, 2579–2605.
equations in computational mechanics via machine 35 H. Talebi, M. Silani, S. P. A. Bordas, P. Kerfriden and
learning: Concepts, implementation and applications, T. Rabczuk, A computational library for multiscale
Comput. Methods Appl. Mech. Eng., 2020, 362, 112790. modeling of material failure, Comput. Mech., 2014, 53,
25 Z. Ren, L. Gao, S. J. Clark, K. Fezzaa, P. Shevchenko, A. Choi, 1047–1071.
W. Everhart, A. D. Rollett, L. Chen and T. Sun, Machine 36 K. S. Ravi Chandran, Transient Joule heating of graphene,
learning–aided real-time detection of keyhole pore nanowires and laments: Analytical model for current-
generation in laser powder bed fusion, Science, 2023, 379, induced temperature evolution including substrate and
89–94. end effects, Int. J. Heat Mass Transfer, 2015, 88, 14–19.
26 S. Sun, A. Tiihonen, F. Oviedo, Z. Liu, J. Thapa, Y. Zhao, 37 J. Y. Huang, S. Chen, Z. F. Ren, G. Chen and
N. T. P. Hartono, A. Goyal, T. Heumueller, C. Batali, M. S. Dresselhaus, Real-Time Observation of Tubule
A. Encinas, J. J. Yoo, R. Li, Z. Ren, I. M. Peters, Formation from Amorphous Carbon Nanowires under
C. J. Brabec, M. G. Bawendi, V. Stevanovic, J. Fisher and High-Bias Joule Heating, Nano Lett., 2006, 6, 1699–1705.
T. Buonassisi, A data fusion approach to optimize 38 J. Lin, Z. Peng, Y. Liu, F. Ruiz-Zepeda, R. Ye, E. L. G. Samuel,
compositional stability of halide perovskites, Matter, 2021, M. J. Yacaman, B. I. Yakobson and J. M. Tour, Laser-induced
4, 1305–1322. porous graphene lms from commercial polymers, Nat.
27 X. Huang, D. Lv, C. Zhang and X. Yao, Machine-learning Commun., 2014, 5, 5714.
reveals the virtual screening strategies of solid hydrogen- 39 S. Khodabakhshi, P. F. Fulvio and E. Andreoli, Carbon black
bonded oligomeric assemblies for thermo-responsive reborn: Structure and chemistry for renewable energy
applications, Chem. Eng. J., 2023, 456, 141073. harnessing, Carbon, 2020, 162, 604–649.
28 B. Deng, D. X. Luong, Z. Wang, C. Kittrell, E. A. McHugh and 40 D. Pantea, H. Darmstadt, S. Kaliaguine, L. Sümmchen and
J. M. Tour, Urban mining by ash Joule heating, Nat. C. Roy, Electrical conductivity of thermal carbon blacks:
Commun., 2021, 12, 5794. Inuence of surface chemistry, Carbon, 2001, 39, 1147–1158.
29 B. Deng, P. A. Advincula, D. X. Luong, J. Zhou, B. Zhang, 41 D. Han, Z. Meng, D. Wu, C. Zhang and H. Zhu, Thermal
Z. Wang, E. A. McHugh, J. Chen, R. A. Carter, C. Kittrell, properties of carbon black aqueous nanouids for solar
J. Lou, Y. Zhao, B. I. Yakobson, Y. Zhao and J. M. Tour, absorption, Nanoscale Res. Lett., 2011, 6, 457.

1218 | Digital Discovery, 2023, 2, 1209–1218 © 2023 The Author(s). Published by the Royal Society of Chemistry