The control of the physical, chemical, and electronic properties of laser -induced graphene (LIG) is crucial in the fabrication of flexible electronic devices. However, the optimization of LIG production is time -consuming and costly. Here, we demonstrate state -of-the-art automated parameter tuning techniques using Bayesian optimization to advance rapid si ngle- step laser patterning and structuring capabilities with a view to fabricate graphene -based electronic devices. In particular, a large search space of parameters for LIG explored efficiently. As a result, high -quality LIG patterns exhibiting high Raman G/D ratios at least a factor of four larger than those found in the literature were achieved within 50 optimization iterations in which the laser power, irradiation time, pressure and type of gas were optimized. Human -interpretable conclusions may be deri ved from our machine learning model to aid our understanding of the underlying mechanism for substrate -dependent LIG growth, e.g. high - quality graphene patterns are obtained at low and high gas pressures for quartz and polyimide, respectively. Our Bayesian optimization search method allows for an efficient experimental design that is independent of the experience and skills of individual researchers, while reducing experimental time and cost and accelerating materials research. 1. Introduction Laser materials processing techniques are powerful tools in materials manufacturing [1-3]. Such approaches have also been applied to form graphene [4, 5] , a two -dimensional material with extraordinary physical and chemical properties [6] . Ground- breaking work focused on 2 direct laser writing of graphene oxide (GO) , a graphene derivative that can be mass -produced via wet chemical pathways [4, 7, 8] . Laser reduction of GO removes insulating oxygen functional groups with precision, which enables micro -patterning and structuring capabilities that are essential for the production of graphene -based electronics [9] . Recently, laser -induced graphene (LIG) patterned on polyimides (PI) substrates has shown potential in fabricating flexible ele ctronic devices [10] . Following the discovery of LIG , precise control over its physical, chemical and electronic properties would allow to broaden the scope of its applications to new areas such as micro -supercapacitors [11, 12] and fuel cell technology [13]. Therefore, there is considerable interest in advancing strategies for property engineering of LIG ; in partic ular, optimizing the lasing parameters and the local environment for different substrates allows to control its composition and morphology. An optimization study is typically concerned with optimizing an objective function that yields a response depending on the values of a given set of parameter s. A straightforward strategy for such optimizations is grid search, whereby the parameter space is sub -divided into a grid and configurations on this grid are sampled for evaluation. Although grid search is success ful for experimental design [3, 14, 15] , it is often an inefficient st rategy with real - valued and dependent parameters in larger dimensions . For example, the Taguchi approach constructs the objective landscape with orthogonal arrays, in which parameter inputs are divided a priori into a grid of possible values [15] , while the Fisher method assumes that parameters take discrete values [14] , such that optimal points that lie between grid points are not taken into account and thus crucial features of the objective landscape can be overlooked. Consequently, such traditional strategies require a vast number of experimentations or computations to sample on a grid that is sufficiently fine , and thus are not suitable for optimiz ation studies in which evaluations of the objective function are expensive . The laser reduction parameters are typically optimized using manual trial and error as an exhaustive grid search is infeasible, and the efficiency of the optimization and end resul t vary depending on the expertise of the individual performing the optimization. In contrast, state -of- the-art techniques, such as Bayesian optimization (BO) and related machine learning (ML) techniques , have recently been applied as an alternative approach in material informatics for high-throughput experiment al designs [16- 18]. The advances in materials informatics has been primarily focused on utilizing large databases of computational work to accelerate the discovery of ne w materials [19- 22]. The optimization of the properties of a new material often 3 requires extensive parameter studies . Some groups have reported high- throughput experimental work, where ML models were sequentially updated with new measurements [23-25], following a general Bayesian optimization approach. We streamline optimization of LIG patterning conditions within a minimum number of iterations by adopting this iterative, adaptive approach . Self-adaptive learning systems leverage experience from past applications of machine learning to achieve better results and improve over time [26] . In particular, we investigate the application of techniques for automated parameter tuning through Bayesian model -based optimization [27]. This technique falls under the umbrella of automated parameter tuning, also called hyperparameter configuration, which aims to find the best parameter configuration of an algorithm for a particular application [28] . Techniques for automated parameter tuning usually proceed in an iterative fashion – they predict the configuration to evaluate, and the result of this evaluation informs the predictions for the configuration to evaluate next. At the heart of these techniques are so -called surrogate models, which is cheap to evaluate and allows for a targeted exploration of the parameter space, identifying promising configurations that available resources for evaluations of the underlying process should be directed towards. This family of techniques is often referred to as model -based optimization (MBO), as the optimization process is based on the predictions of the surrogate models that serve as a replacement for the underlying process. Surrogate models are induced using ML, taking an increasing amount of gr ound- truth data into account between subsequent iterations. In this study, we describe MBO -assisted laser reduction of graphene on quartz and polyimide. We discovered experimental conditions that lead to LIG patterns with Raman G/D ratios of at least a factor of four higher than those found in the literatur e. High Raman G/D ratios indicate a high degree of graphene formation in our LIG patterns, which restores the unique properties that are desirable in graphene -based devices. Additionally, we demonstrate the human interpretability of the model that gives insight into the substrate -dependent LIG processing conditions. The mechanisms which govern the formation of LIG patterns are discussed. Since the MBO approach presented here suggests promising configurations iteratively, we consider this work an initial step in the development of human- in-the-loop automated LIG patterning and characterization systems. 4 2. Experimental s etup and m ethods 2.1. Graphene oxide films Graphene oxide (GO) was synthesized from graphite using the improved Hummers’ method, as it avoids the generation of hazardous gases such as NO and NO 2 from nitrate that is used in the conventional Hummers’ method [29] . In this procedure, powder ed samples , ground and sieved to 20 µm , were mixed in concentrated H 2SO 4 and H 3PO 4 and placed in an ice bath. Then, KMnO 4 was added at a mixture temperature of 35 ˚C and increased further to 98˚C before termination with ultrapure water (Millipore) and H 2O2. The filtrate was then washed with HCl and subsequently with water repeatedly until a pH -level of about 6.5 was obtained. The GO inks were produced using 25 mg of the freeze -dried GO powder, which was diluted in 100 ml deionized water and ultra -sonicated with a cooling system. After the sample was centrifuged, the remaining supernatant was repeatedly diluted and ultra -sonicated until a 200 ml dilution was obtained. The GO inks were spray -coated onto a 1 cm2 quartz or polyimide substrate (Kapton HN 125 μm, Dupont) in multiple passes until a thickness of 1 μm was achieved . The GO thickness was verified by scratching the surface and measuring the height using an optical profilometer (VK -X1000, Keyence). 2.2. Laser- induced graphene Laser -induced graphene (LIG) spots were patterned by reducing GO films deposited on quartz and polyimide, and by carbonization of polyimides directly. We denote GO/Quartz, GO/Polyimide and Polyimide as samples GOQ, GOPI and PI, respectively. The patterning setup is shown in Fig. 1. The deposited GO films were placed in a sample chamber which allows patterning in an air, argon or nitrogen environment with pressures up to 1000 psi. LIG patterns were irradiated using a 532 nm diode -pumped solid- state continuous -wave (CW) laser (Sprout G -12W, Lighthouse Photonics). The laser beam was focused with a 50× microscope lens to a spot size of 20 μm on the sample surface. Irradiated beam spots were positioned sufficiently far apart from each other to ensure pri stine precursor material for each experiment. The sample area is about 1 cm2, allowing approximately 256, 25, and 25 patterns for samples GOQ, GOPI, and PI, respectively. Taking into account the sample preparation and repeated measurements, we set our expe rimental budget to the maximum number of successful patterns carried out for our most spatially limited material within a work week, 5 which is 70 runs for PI -based samples. The same number of runs were performed for quartz - based samples for consistency. Figure 1: Schematic of the laser patterning and in situ Raman measurement setup . The quality of the LIG spots was determined using a Raman spectrometer (Isoplane SCT320, Princeton Instruments). As the sample is irradiated, the laser beam is backscattered and filtered through a long -pass filter to increase the sensitivity of the signal. Using the same laser source for patterning and Raman spectroscopy, the identical spot was characterized in - situ. The Raman data for each spot were averaged over 10 measuremen ts with a collection time of 3 s at laser power <10 mW for each measurement. The Raman spectra were post - processed with a linear background subtraction to 0 and normalization of the maximum peak to 1. The G - and D -bands were fitted using Lorentzian functions and the ratio of their intensities computed as the ratio of the areas under the fitted functions. The G/D ratios indicate the degree of reduction of GO and is used as a proxy for electric conductivity. 2.3. Model -based Optimization The schemati c in Fig. 2 sh ow the procedure of the machine -learning -assisted LIG patterning based on the Bayesian model -based o ptimization ( MBO) algorithm , which is a sequential approach to optimize a black -box function f(x) [30] . In this study, we use the mlrMBO package [31] to model the parameter space, build the surrogate models and determine the most promising configuration for the next evaluation. We do not optimize parameters one- at-a-time [32] , but simultaneously to account for possible interactions between them . Here, the G/D ratio = f (x) is the objective function specific to our LIG patterns, and x is the process configuration consisting of a complete assignment of values to all parameters, i.e. CW -laser power, irradiation time, gas pressure and type of gas (see Table 1 for limits) ; for simplicity, we denote ‘ G/D ratio’ as ‘ratio’ throughout this study unless stated otherwise. The data {𝑥𝑥 𝑛𝑛,𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑛𝑛}𝑛𝑛=1𝑁𝑁 obtained from the past N LIG patterns and ratio measurements are used to build a surrogate model to predict the value of f (x) at different x, 6 substituting the expensive underlying process. As we have a mixture of continuous numerical and discrete categorical parameter s, we use a random forest with 500 trees as our surrogate model . Specifically, the predicted ratio for a point x corresponds to the mean of the predictions across all trees μ( x) and the uncertainty σ( x) is estimated as the sum of variances over all t raining points. The point -wise variance is computed as the average of the jackknife - after -bootstrap and infinitesimal jackknife variance estimates [33] , which effectively captures uncertainty due to finite size of the training data. Figure 2: Schematic of the model -based Bayesian optimization approach used in this study. We use random forests as the surrogate model as the input consists of both continuous numerical and discrete data. The acquisition function is Expected Improvement for maximizing the Raman G/D ratio of LIG. The objective of the MBO is to maximize the ratio within an experimenta l budget of 70 runs, which is the largest number of successful patterns performed for our most spatially limited material within a work week . We split the total budget into initial training data of 20 randomly selected configurations to fit the initial sur rogate model and 50 iterations of the MBO. The acquisition function proposes the next configuration to evaluate based on the predictions made by the surrogate model . We use the Expected Improvement (EI) acquisition function, which searches for the configur ation that most likely has a better objective value than the current best measured point [31] . More concretely, the EI takes both the mean μ( x) and uncertainty σ(x) into account to trade off exploration and exploitation – once the uncertainty becomes sufficiently small for a known good region, the larger uncertainty of unexplored areas of the parameter space will guide it towards that [31] . 7 Our MBO approach applies the focus search algorithm to identify the most promising config urations according to the acquisition function. It shifts from coarser to finer grids for sampling configurations to evaluate the surrogate model on, as even with a cheap surrogate model an exhaustive grid search is prohibitively expensive. The finer grids are centered around the areas where the best configurations were discovered in previous iterations. In this work, the shift is iterated 10 times and the best point over all iterations is chosen. Note, there is no guarantee that this process will find the globally optimal configuration, and no guarantees of convergence. This is because there is no closed functional form for the process we optimize and hence only an exhaustive evaluation of the entire parameter space would allow to identify the global optimu m and prove that it is the global optimum. As the evaluations of the process are very expensive (running an experiment), this is infeasible. To summarize, we begin the iterative parameter tuning process by building an initial surrogate model using random parameter values and selecting the next configuration to evaluate based on the results of the focus search and the acquisition function. Thus, the model proposes new LIG parameters for the next run and the resulting measurement is used to update the surroga te model for the next iteration. This process is iterated until we reach 50 iterations, our stopping criterion. To account for the randomness of the initial data, we performed three experimental runs for each material. The MBO code, dataset and metadata are made available in the supplementary information following the FAIR guiding principles [34]. Table 1 : Parameter space limits for MBO Parameters Lower limit Upper limit Instrument precision Number of possible values CW-laser power [W] 0.01 5.55 0.01 554 Irradiation time [s] 0.500 20.000 0.001 195000 Gas pressure [psi] 0 1000 10 100 Gas type Argon Nitrogen Air - 3 3. Results 3.1. Maximization of G/D ratio in LIG Raman spectroscopy is the traditional technique to characterize the structural and electronic properties of graphite . It relies on the inelastic scattering of laser photons after they interact with the vibrating molecules in the sample probe. For graphite, Raman experiments 8 extract characteristic information such as the de fects (D -band at ~1350 cm-1), in-plane vibrations of sp2 carbon atoms (G -band at ~1580 cm-1) and the stacking order of the carbon basal planes (2D -band at ~2700 cm-1) [35] . The Raman spectrum of GO and LIG can be distinguished most notably by the presence of the 2D -band in Fig. 3. Although the 2D peak can be fitted with only one Lorentzian peak at ~2700 cm-1, similar to graphene [35] , the full- width half -maximum (FWHM) is much larger ~60 cm-1 in our LIG , as reported by other s [10, 36]. In contrast, the 2D -band in bi layer or multilayer ed graphene is commonly fitted with more than one Lorentzian function [35] . Although the 2D -band was not explicitly optimized in this study, its characteristics are further discussed in Section 3.5. Finally , the G/D ratio indicates the degree of graphitization in the LIG patterns . Figure 3: Raman spectra of graphene oxide (GO, bottom) and laser -induced graphene (LIG, top). The spectra are offset for clarity. A series of three experimental campaigns where conducted for each material. The experimental setup considered parameters and objective were the same, but different batches of GO inks, w hich were manufactured using the same process, were used to deposit the GO films and different training data points were used. The aim of running a series of experiments was to show that our method to maximize the ratio works with slight variations in the precursor material, as two thin films of graphene oxide are not exactly the same even if they are manufactured using the same process, and to show that we can achieve improvements with different initial training data. Figure 4 shows the best -performing opt imization runs. For all materials, the random initial parameters for the initial surrogate model are depicted as scatter plots up to the 0th iteration. Note that the performance improves significantly as soon as we explore configurations that were suggested by the model -based parameter optimization approach. Generally, the 9 predictions become more accurate as the distance between the actual and predicted values become smaller with each iteration and as more data are added to the surrogate model. Significant performance improvements are observed in all three campaigns, as evident from the t-test results between the randomly sampled and MBO -optimized patterns (Table S1); the corresponding box plots comparing such differences are shown in (Fig. S4). LIG Patterning Parameters Ratio Ref. 10.6µm, 14µs, pulsed 2.25 [10] 10.6µm, 4.8W, CW 1.38 [36] 532nm, 100mW, CW 0.73 [60] 532nm, 12mW, CW 1.03 [55] Figure 4: Progress of G/D ratio optimization for various materials in campaign 1. The box plots show the statistics of the achieved ratios before (iteration <1) and during optimization (iteration >= 1); box plots for all campaigns are shown in the supplement (Fig. S4). The results achieved during the optimization are statistically significantly different from the results achieved before for each sample, as evident from t -tests (Table S1). The G/D ratios are improved by a factor of about 3, 4 and 2 in 9, 13 and 1 optimization iterations for GO/quartz (GOQ), GO/polyimide (GOPI) and polyimide (PI), respectively. The grey bounds show the uncertainty of the predictions. State -of-the-art literature values are tabled for comparison. 10 SEM of the LIG patterns are shown in Fig. 5. Before optimization, the random processing parameters led to LIG patterns with significant d amage at the center of the irradiated spot. After optimization, the irradiated pattern is reduced in diameter with no signs of damage . The porous structures commonly reported in LIG patterns on PI substrates are seen [10] . In contrast, no such structures are observed in GO on quartz substrates, suggesting that they are substrate dependent . Figure 5: SEM of the LIG patter ned (top) before and (bottom) after optimization for (a),(b) GOQ (c) ,(d) GOPI and ( e),(f) PI . Scale ba r is 100 µm. 3.2. Interpretation s of the surrogate model Model interpretation is of paramount importance in any ML study. T he importance of each parameter can be determined to gain insight into structural features that impact the LIG patterns. We apply partial dependence analysis to obtain post -hoc interpretations of our surrogate models . Similar to descriptor significance analysis, the partial dependence elucidates the effect of a single parameter on a property by marginalizing over all other parameters [37] . For a more robust interpretation of the model, we have trained the surrogate model by combining all the data points from all three campaigns. The partial dependence for the GOPI and PI samples show comparable behavior (Fig. S5). Considering these results, we have combined the data points from six experimental campaigns for both materials into one. The combinations resul t in a total of 210 and 420 data points for the LIG on substrates quartz (LIGQ) and polyimides (LIGPI), respectively. 11 3.3. Partial dependence of lasing parameters The partial dependence plots in Fig. 6 a – 6f visualize the marginal relationship between the conti nuous LIG process parameters and the prediction of the ratio. Note that the ratios in these plots are not measured but predictions made by the surrogate model trained on this data. Therefore, the ratios are not as high as the maximum measured data in the p revious section. The grey area denotes how the values vary with all possible values for other parameters i.e. irradiation time, gas pressure and type. For the type of gas, the actual measurements are shown as we cannot compute partial dependencies for cat egorical parameters. The partial dependence for laser power shows that t he optim al ratio in Fig. 6 b is found between 1.9 W – 2.5 W . This is in agreement with the threshold power of 2.4 W, at which conversion from PI to LIG occurs [10] . Note th at mos t studies on PI were done with µ m- pulsed 10.6µm lasers, while this study , to the best of our knowledge, is the first LIG investigation using a CW 532nm laser ; in PI, the absorbance at 10.6µm is twice that of 532nm [38]. Although we expect some therma l effects from the CW laser, the SEM micrographs in Fig. S 5 show that, apart from some striations and cracking of the GO thin film, no significant differences are observed in the morphology between LIG patterned on GOPI and PI . For LIGQ, the optical profilometer micrographs of patterns irradiated at high laser power show depressions in the center of the pattern with no visible melt zones near the edges , which is consistent with a quasi -phase explosion process observed by other groups [39, 40] . As the CW laser beam is absorbed, the vaporization of interlayer water builds up internal pressures between the GO sheets until the ma terial is ejected through rapid local heating. The ejecta re- deposit as clusters spattered outwardly in a 380 µm diameter from the irradiation spot ( Fig. S6b); these clusters are not visible in lower power regimes (Fig. S 6a). Interestingly, clusters were observed for LIGQ but not for LIGPI, which indicates that their origins are substrate - dependent . Additional micro -Raman mapping shows that the clusters in the outer region are not sp 2-hybridized (Fig S6c ), which suggests that the re -deposited material may be re - oxidized, as observed by others [41] . In contrast, high ratios at the edges of the depression center indicate that graphitization occurs following explosive ejection and re -deposition near the irradiated area. 12 Figure 6: Partial dependence plots of G/D ratio for different parameters for (a, b) laser power , (c, d) irradiation time , gas (e, f) pressure and (g, h) type of LIG on quartz (LIGQ) and on polyimides (LIGPI) , respectively . In both cases, the best ratios are found at lower irradiation times; between powers o f 1W -3W for LIGPI; and low and high pressures for LIG on quartz and polyimide substrate, respectively. The predicted ratios in Fig. 6 c and 6d show a decreasing trend with irradiation time. In contrast to pulsed laser systems often used in the LIG literatur e [42] , the uninterrupted CW laser results essentially in an infinite ‘‘pulse width”. The LIG mechanism for CW lasers 13 involves local heat deposition, which deoxygenates the irradiated surface and reorganizes sp3 to graphene -like sp2 structures [7] . However, prolonged exposure produce d cluster regions that increased in diameter with time (Fig. S7). We believe the extended transfer of energy to the surface induced continuous material ejection and eventually ablates the region within the beam diameter of the laser, causing structural dam age that disrupts sp2 carbon bonds, resulting in a drop of ratios. The surrogate model has learned exactly this effect, as seen in the drop of predicted ratios in Fig. 6 c. Discussions of the partial dependence of laser power and irradiation time on ratios by themselves are rare to find in the literature. This is partly due to the fact that laser power is by itself “time -averaged” and that strong interactions between the two parameters may elucidate LIG mechanisms via instantaneous energy transfer onto the m aterial [43] . Bischl et al . have quantified the post -hoc interpretability of models and show that higher interaction strengths results in less relia ble partial dependence analyses [37]. In our ANOVA tests in Tables S 2 and S3 we show that the interactions between power and irradiation time are insignificant, suggesting that our partial dependence analysis is reliable. 3.4. Local environment, plume and thermal effects A controlled surrounding environment can promote or inhibit chemical reactions, ablation, melting, etc. that are relevant to the LIG patterning process. Here, the LIG approaches in inert gas show the best performance, as seen in Fig. 6 g and 6h. The poor ra tios achieved in air may be due to oxidative burning in thick GO [44] , resulting in consistent preclusion of this gas type from the configurations proposed by the MBO . In Fig. 6 e, the ratio in LIGQ is predicted to increase towards lower pressures, implying that LIG patterning on quartz in vacuum atmosphere helps to reduce defects. This is in agreement with Sokol ov et al ., who demonstrated successful laser reduction of graphite oxide in vacuum or high purity nitrogen inert gas [41] ; our surrogate model has predicted best performance in both environmental conditions. The role of the inert gas in LIG may be related to a photothermal pathw ay based on the plume formation, i.e. ablated material being ejected normal to the surface plane that is shaped by pressure gradients [41, 45] . In vacuum, the plume expands freely. In inert gas, the expansion is restricted until the pressure within the plume equals the gas pressure. The plume acts as a shield to prevent beam penetration and that generates heat in the material, which, in 14 turn, carbonizes the film. Under these conditions, the low -pressure conditions that are optimal for high ratios in LIGQ suggests large plume expansions of high temperatures. In contrast, the optimum in LIGPI at high pressures confines the plume in a small area at lower temperatures than LIGQ . This is observed clearly by comparing the heat -affected zone s, i.e. darkened region outside the delineated circular boundary, as observed in the optical micrographs in Fig. S8, Moreover, the high pressure may reduce the PI substrate from swelling when ablated at increased fluence [45] . The local heating in LIGQ and LIGPI is evident from the Raman shift in Fig. 7. The incident laser is absorbed by the GO or PI surface, causing an increase in local temperature that was observed using in situ Ram an spectroscopy. The increase in temperature can be estimated by the shift of the G -band according to T = (ω − ω 0)/χ, where ω 0 is the frequency of the G -band when temperature T is extrapolated at 0 K, and the temperature coefficient χ is 1.5·10-2 cm−1/K and 2.2· 10-2 cm−1/K for graphene [46] and PI [47] , respectively. It is justified to use the temperature coefficient of graphene for GO here as the thermal relaxations that govern the laser reduction process is in the ps -scale [48] , which is a negligible fraction of the irradiation times used in th is study. For LIGQ, the G-band is downshifted by 14 cm-1 from the base position of the G -peak ~1585 cm-1, at low pressures (<250 psi) , increasing local temperature by ~ 933 K . For LIGPI, the G -band shift of 15 cm-1 at high pressures (>750 psi) increased the local temperature by ~ 750 K; the position of the G -band varies widely and overlaps in both pressure environment s, suggesting higher sensitivity to thermal effects , as expected from soft substrates such as PI. The differences in temperature increase show that the underlying substrate plays an important role as a heat sink in graphene formation, as observed by others [44] . As the thermal diffusivity of quartz is a factor of seven larger than PI [49], we believe that the net photothermal effect from the plume and the heat sink promotes LIG patterns to form. 15 Figure 7: Position of the G -band at low (<250 psi) and high (>750psi) pressures for (a) LIGQ and (b) LIGPI samples. In contrast to LIGPI, the Raman shift in LIGQ at high gas pressures indicate a larger local temperature increase. The quality of LIG depends on the net thermal effect and the underlying substrate. At high pressures for LIGPI, the localized heating is restricted, promoting carbonization. On quartz, high- quality LI G is patterned at low gas pressure that expands local heating. 3.5. 2D-band characteristics of LIG For a rigorous identification of graphene or LIG, the analysis of the G/D ratio alone may be insufficient as this parameter may be affected by the amorphization t rajectory described by Ferrari and Robertson [50] . Hence, the characteristics of the 2D -band , particularly with respect to G/D ratio that has been optimized in this study , will be discussed here. The analysis of the 2D /G parameter is presented in Fig. 8 shows that, a part from some outliers, the 2D -band is mostly present for both LIGQ and LIGPI with a larger 2D/G distribution for the latter, presumably due to the substrate -dependent thermal effects discussed in the previous section. The 2D-band i s an overt one of the D -band, but unlike the latter, it is not activated only at proximity to defects. As a result , the 2D -band can be present in LIG even at the absence of the D -band and it is not sensitive to defects . The 2D -band can be tuned to determine graphene layer thickness . For the monolayer graphene , the 2D -band is observed to be one symmetrical peak with a FWHM of ~30 cm-1 [51]. With the addition of successive layers, others report broadening and splitting of the 2D - band such , that the distinct band shape can differentiate between single and multilayer graphene for layer thickness of less than 4 layers [51] . Although the FWHM reported here generally decreases with increasing G/D ratio , the 2D-bands in our samples are typically 16 twice as broad as those reported for monolayer graphene , as seen in other LIG reports [10, 36]; they have a symmetrical profile ( Fig. 8 d) more closely related to 2D graphite with its layers randomly stacked along the c-axis [52] . Figure 8 : The 2D/G ratio for (a) LIGQ and (b) LIGPI as a function of G/D ratio, which was optimized using MBO. The color bar shows the F WHM(2D) in cm-1; FWHM(2D) and G/D ratios are negatively correlated (Spearman coefficient: (a) - 0.79 and (b) -0.64). (c) Boxplot for the 2D/G ratio distributions and (d) symmetrical 2D -band profile for both LIG substrates. It is worth pointing out that monolayer graphene can be identified by analyzing the ratio between the 2D - and the G-band. The 2D/G ratio s in Fig. 8a and 8b generally increase with G/D ratio up to G/D ~ 2 and plateau at median s of 0.32 and 0.38, respectively . The maximum 2D/G values presented for LIG here are comparable to those reported in the literatur e (see Table 2) . Note, in comparison to the G/D parameter, the 2D/G parameter for LIG is not widely reported and the values presented here are roughly estimated . These values range between a quarter to half of the values typically observed in defect -free single layer graphene materials (2D/G ~ 2) [53] . Interestingly, to the best of our knowledge, the D -band does not 17 diminish through the laser induction process and thus high- quality LIG are not defect -free. Therefore, it is reasonable to take the D -band into account during single -objective optimization of G/D ratios. Nevertheless, through additional analysis of the 2D/G ratio , FWHM (2D) and the 2D -band profile , we have deduce d high -quality LIG. 3.6. Related work The best ratios averaged over the three LIG campaigns are 5.3 6, 7.13 and 5.38 for GOQ, GOPI and PI, respectively, and were achieved within 50 optimization iterations . To note, while this has cost a week of experimentation, prior study without MBO -assistance required not only over a month’s work but also yielded suboptimal results [54] . While the base ratio was measured to be about 1 for unreduced GO films, the characteristic Raman G -, D- and 2D - features were not found in non- irradiated polyimides. A list of LIG parameter conditions in the literature with the corresponding ratios are shown in Table 2. To the best of our knowledge, other groups with comparable experimental setups, i.e. 532 nm CW laser patterning of GO thin films, demonstrate a maximu m ratio of 1.03 [55] . Hawes et al. have carbonized PI substrates using CW CO 2 laser and demonstrated a ratio of 1.38 [11] . While we acknowledge that a higher ratio of 2.25 has been observed in LIG on PI [10] , it is not trivial to have a direct comparison of our work with that in the community as the results of laser - material interactions are highl y dependent on the characteristics of the laser system e.g. pulse width and wave length. We therefore prioritize and compare studies using CW lasers , since the underlying mechanism for CW laser reduction is thermally mediated [7] . The introduction of automated parameter tuning with surrogate models improves performance in a small number of iterations with relatively little experimental ef fort by a factor of at least four . We emphasize that the advantage of automated parameter tuning used here is not limited specifically to our experimental setup but can be applied in other context s as well, e.g. pulsed laser systems. 18 Table 2: List of l aser-induced g raphene (L IG) parameter conditions from various precursor materials and their Raman G/D and 2D/G ratios Laser Type Patterning Parameters Precursor Material/Substrate G/D Ratio 2D/G Ratio Ref Pulsed 790nm, 120fs GO/glass 1.10 0.13 [56] Pulsed 10.6µm, 2.4W Hydrated graphite oxide film 1.28 - [57] Pulsed 355 nm, 20ns Graphite oxide film 1.15 0.05 [58] Pulsed 522nm, 500fs PI 0.86 0.06 [11] Pulsed 10.6µm, 14µs PI 2.25 0.95 [10] Pulsed 10.6µm, 14µs PI 1.03 0.46 [59] CW 10.6µm, 4.8W PI 1.38 0.48 [36] CW 532nm, 100mW Thin film GO/glass 0.73 - [60] CW 663nm, 80mW GO/Quartz 0.93 - [61] CW 788nm, 5mW Graphite oxide/polycar bonate 0.94 0.26 [62] CW 532nm, 12mW Thin film GO/SiO 2/Si 1.03 0.13 [55] CW 532nm, <5.5W, <1000 psi Thin film GO/quartz GO/PI PI 5.36 7.13 5.38 0.49 0.52 0.51 Our MBO approach 3.7. Discussion There are several avenues for future studies. First, despite our reliable model interpretations, we believe that our understanding of the underlying LIG processes has not necessarily improved. With a view to develop explainable artificial intelligence systems, leveraging these models would allow us to close the scientific loop and have AI work in tandem with human researchers to not only achieve better experimental results, but also to develop new theories. More explicitly, the MBO partial dependence analysis has guided us to find the physical domains to investigate further post -optimization , and from that we formulated a substrate- dependent mechanism in LIG formation based on thermal and plume effects in a local environment. F uture work could include ML approaches that considers the 19 discrepancies between simulation and experiments , which is expected to outperform random ML training approaches. The analysis of the partial dependencies we performed in this paper is a valuable first step to integrate phys ics within ML. Second, while we have achieved good ratios, the simultaneous optimization of other objectives, such as the electric and thermal properties, can be beneficial to fabricate electronic devices. A common problem in multi -objective optimizations is that one property of the material cannot be improved without degrading another. Instead, an optimal boundary, the so-called Pareto front, defines the trade -off between configurations that are optimal with respect to at least one property. Multi -objectiv e optimization has been demonstrated to be data exhaustive [63, 64] . Most materials discovery studies therefore focus on single -objective optimization, which can support multi -objective optimization by combining all objectives into one. Another practical challenge is that often measurements can only be taken at multiple levels. For example, Raman characterization of a material can happen immediately after a LIG patterning has been performed, but measurements of electric conductivity using a four - point probe can only be performed after an experimental campaign is finished as this requires removal of the sample from the reaction chamber. The integration of such multi -level measurements that occur at diff erent frequencies into a Bayesian optimization process would potentially allow further improvements. 4. Conclusion We have demonstrated the application of state -of-the-art automated parameter optimization techniques to laser -induced graphene (LIG) patterning and improved the quality of LIG patterns based on Raman G/D ratios reported in the literature by at least a factor of four. The best LIG ratios were achieved within 50 configurations suggested by the MBO . The initial surrogate models were trained with rand om parameter evaluations that are independent of the skills and experience of individual researchers – our method does not rely on skilled operators to “guide" it towards good results. Further, our system does also not require a background in machine learn ing, and the achieved improvements are demonstrated to be reliable and reproducible. We further demonstrate how the surrogate model can be used to improve our understanding of the underlying processes for LIG patterned on quartz and on polyimide. The parti al dependence analysis shows that the LIG mechanism is governed by substrate -dependent local thermal effects and plume formations. This conclusion is supported by optical characterizations. 20 Acknowledgements The authors gratefully acknowledge the financial support by of the University of Wyoming’s Engineering Initiative of the College of Engineering and Applied Science and School of Energy Resources. 21