Determining the Parameters of Importance of a Graphene Synthesis Process Using Design-of-Experiments Method

A systematic method to identify key factors that control the synthesis of Physical Vapor Deposition (PVD)-based graphene on copper is necessary for engineering graphene growth. The statistical design-of-experiments method is employed and demonstrated in this work in order to fulfill the necessity. Full-factorial design-of-experiments are performed to examine the significance of the main effects and the extent of the interactions of the controlling factors, which are responsible for the number of layers and the quality of the grown graphene. We found that a thinner amorphous carbon layer and a higher annealing temperature are suitable for the growth of mono-layer/few-layer graphene with low defects, while the effect of annealing time has a trade-off and needs to be optimized further. On the other hand, the same treatment, but with larger annealing times will result in multi-layer graphene and low defects. The results obtained from the analysis of the design-of-experiments are verified experimentally with Raman characterization.


Introduction
Graphene is a planar sheet of sp 2 bonded carbon atoms with single atomic thickness [1].It is increasingly popular due to its splendid properties.Excellent intrinsic mobility of about 200,000 cm 2 ¨V´1 ¨s´1 [2], current density tolerance of about 10 8 A/cm 2 [3] and thermal conductivity in the order of 5 ˆ10 3 W/m¨K [4] are among the few important properties of graphene that make it a promising candidate for a vast number of applications, such as ULSI (Ultra-Large Scale Integration) interconnects, transparent conductive coatings, photodetectors, optical modulator, sensors, metrology, etc., and many others.
Specifically, graphene with different numbers of layers have different potential applications.Single-layer or mono-layer graphene (SLG) is used as the top electrodes of semi-transparent organic solar cells fabricated by Liu et al. [5], as the transparent conductive anode in OLEDs (Organic Light-Emitting Diodes) investigated by Zhu et al. [6], as field-effect transistors studied by Li et al. [7], etc. Eda et al. [8], Kim et al. [9], Wang et al. [10] and Yan et al. [11] reported the applications of few-layer graphene (FLG) in thin-film transistors (TFTs), capacitive touch panels, touch screens and micro/nanometer-scale heat spreaders, respectively.On the other hand, multi-layer or many-layer graphene (MLG) films are used as thermal interface materials in photovoltaic (PV) solar cells, optoelectronic, photonic devices and systems [12].MLG films are also used as superconductors reported by Ballestar et al. [13] and as fillers in polymer composite reported by Chrissafis and Bikiaris [14] and Jia et al. [15], etc.
Several graphene growth methods have been reported since its discovery in 2004 by Novoselov et al. [1].Beginning from a scotch-tape technique [1] in order to mechanically exfoliate graphene flakes, many growth techniques have been reported [16][17][18], such as graphite sonication [19], epitaxial growth [20], graphene oxide reduction [21], etc.The chemical vapor deposition (CVD) technique produces large area graphene with high quality and low mass production, which otherwise is a big challenge for other growth methods [16,22,23].The CVD method uses a variety of carbon sources, including gaseous [7,[24][25][26] and liquid [27,28] precursors.Use of amorphous carbon as a solid source for graphene synthesis on Cu was reported by Ji et al. [29].However, their method involved graphene growth on a several micron-thick copper foil, and the graphene so obtained needs to be transferred onto the required substrate, which might introduce defects.
Recently, we demonstrated the feasibility of the crystallization of amorphous carbon (a-C) below sputtered copper (Cu) thin film, which acts as a catalyst to obtain graphene on the top surface of the Cu film experimentally [30].This graphene synthesis method is PVD (Physical Vapor Deposition) based, which is novel as both a carbon source and a metal catalyst are deposited using the sputtering method.Unlike other methods, such as conventional CVD (Chemical Vapor Deposition), where much of the gaseous carbon source is wasted and the synthesis temperature is very high, this PVD-based method provides low cost and direct transformation of pure carbon species underneath the copper film into graphene on the top of copper film at temperatures acceptable for the semiconductor industry, with better control over the number of layers and the quality of the grown graphene.Therefore, it can potentially be useful for enhancing the conductivity and reliability performance of ULSI copper interconnects.
In order to optimize the above-mentioned PVD-based process, the key process parameters (also known as factors) must be identified.One-factor-at-a-time (OFAT) experimentation is the most common method to identify the key factors, in which only one factor or variable is varied at a time while keeping others fixed [31].However, statistical approaches, such as the Design-of-Experiments (DoE) method, are more efficient in order to study two or more factors that vary simultaneously, while keeping the number of experimental runs to a minimum.Furthermore, DoE helps in obtaining the necessary information, especially the interaction effect of the various parameters [32], which is not possible with the OFAT method.DoE has been used as a statistical tool for several multi-disciplinary engineering applications, such as controlling the process parameters for various nano-suspensions [33] and nano-carriers [34][35][36].Furthermore, the optimization of biodiesel production reported by Vicente et al. [37] and surface roughness prediction reported by Choudhury and El-Baradie [38] also use this powerful tool.Wirtz et al. [25] implemented DoE in order to identify combinations of parameters that are suitable for graphene growth on copper using ethene.However, the implementation methodology of the latter work is not clear, and the reasons for the identified key parameters being important to the process are not explained in their work.
This work attempts to describe the implementation methodology of DoE for our PVD-based graphene synthesis method.The physical reasons for the identified key factors in affecting a specific outcome of the grown graphene film are also explained to demonstrate the consistency of the DoE method and the underlying physics of our PVD-based graphene synthesis method.

Design-of-Experiments
The full-factorial Design-of-Experiments (DoE) [32] method is implemented in order to study the effect of the levels of process parameters or factors necessary to facilitate PVD-based graphene growth on copper (Cu) using amorphous carbon (a-C) as the solid carbon source.There are many factors that can affect the growth of graphene on Cu in this process, namely the annealing temperature, annealing time, a-C and Cu layer thicknesses, gas mixture, gas flow rate and pressure values.
From our previous experiences with the process, we found that varying some parameters, such as gas mixture, gas pressure, flow rate and Cu thickness, will affect the possibility of graphene growth, and thus, they are held constant to values that can facilitate graphene growth in this work.This is because our focus here is on the values of parameters (levels of factors) that will affect the number of layers (single, few or multi-layered) of graphene and the quality (corresponding to the amount of defects) of graphene growth.Thus, we limit our scope to the study of only three parameters in this work, namely annealing temperature, a-C layer thickness and annealing time.
The thickness of the a-C layer is further determined by the combination of various parameters in the sputtering process of a-C layer, including the substrate temperature, the argon (Ar) flow rate, base pressure values, RF/DC (Radio-Frequency/Direct Current) power and deposition time.In order to achieve variable a-C layer thickness, deposition time is varied, while the other parameters involved in the sputtering process are held constant.
Each of the parameters or factors examined have two levels, as shown in Table 1.The factors are named as "A", "B" and "C", and the levels are coded as "´1" or low level and "+1" or high level.

Graphene Synthesis
The graphene synthesis process comprises sample preparation using the PVD method and subsequent annealing.
The samples are prepared by depositing amorphous carbon (a-C) thin film with different thicknesses on Si/SiO 2 (300 nm) substrate, followed by 800 nm-thick copper (Cu) (99.99%) film deposition.Supplementary Figure S1 shows a detailed schematic diagram of sample preparation.The choice of the thickness of the a-C layer is based on the levels depicted in Table 1.
The deposition is carried out using RF and DC sputtering for a-C and Cu respectively at a stable pressure of 3 mTorr, in the presence of Ar gas with a flow rate of 30 sccm.The substrate temperature during the sputtering process is maintained at 250 ˝C.
All of the samples are subsequently annealed in a hydrogen (H 2 ) environment with a flow rate of 50 sccm at a low pressure of 1 Torr.The annealing temperature and time are varied as per the levels chosen for DoE, as shown in Table 1.After annealing, the samples are cooled down, during which the H 2 flow rate is decreased to 30 sccm, and Ar gas is introduced at a flow rate of 20 sccm.Supplementary Figure S2 shows schematic diagram of annealing apparatus.

Characterization
The characterization of the annealed samples is performed using Raman scattering spectra, and the equipment used is PTT RAMaker Micro Raman/PL/TR-PL Spectrometer (Protrustech Co., LTD, Tainan, Taiwan) with a confocal Raman microscope system (laser excitation wavelength = 473 nm; laser spot-size = 0.5 µm).
Raman characterization is an important tool to determine the presence of graphene by observing 3 signature peaks, namely I D , I G and I 2D .It is well known that the I G peak is associated with the doubly-degenerate phonon mode at the Brillouin zone center, which comes from a normal first order Raman scattering process in graphene and originates from the interplanar sp 2 -bonded C-C stretching vibrations [39].The I 2D and I D peaks originate from a second-order process, involving two iTO phonons near the K point for the 2D band and one iTO phonon and one defect in the case of the D band, which can be often used to evaluate the number of graphene layers and the grain size of graphene, respectively [39].In general practice, I 2D /I G and I D /I G peak intensity ratios are used as metrics for evaluating the number of layers and quality of graphene, respectively [39][40][41].Hence, we use these metrics as the responses for DoE analysis.

Factorial Experiment Design and Experimental Results
The analysis is done on the design based on full-factorial DoE [42], which is a two-level factorial analysis, and all three factors are taken into consideration.The total number of runs performed is n ˆLF = 2 ˆ23 = 16, in a random order, where n is the number of replicates, L is the number of levels and F is the number of factors.All eight combinations of factors "A", "B" and "C" are chosen as shown in Table 2.The values of I 2D /I G and I D /I G peak intensity ratios (obtained from Raman characterization) for each of the eight sets of factor levels and each of the two replicates are collected, averaged and recorded in Table 2. Figure 1 shows the Raman Spectra for the experimental runs mentioned in Table 2. degenerate phonon mode at the Brillouin zone center, which comes from a normal first order Raman scattering process in graphene and originates from the interplanar sp 2 -bonded C-C stretching vibrations [39].The I2D and ID peaks originate from a second-order process, involving two iTO phonons near the K point for the 2D band and one iTO phonon and one defect in the case of the D band, which can be often used to evaluate the number of graphene layers and the grain size of graphene, respectively [39].In general practice, I2D/IG and ID/IG peak intensity ratios are used as metrics for evaluating the number of layers and quality of graphene, respectively [39][40][41].Hence, we use these metrics as the responses for DoE analysis.

Factorial Experiment Design and Experimental Results
The analysis is done on the design based on full-factorial DoE [42], which is a two-level factorial analysis, and all three factors are taken into consideration.The total number of runs performed is n × L F = 2 × 2 3 = 16, in a random order, where n is the number of replicates, L is the number of levels and F is the number of factors.All eight combinations of factors "A", "B" and "C" are chosen as shown in Table 2.The values of I2D/IG and ID/IG peak intensity ratios (obtained from Raman characterization) for each of the eight sets of factor levels and each of the two replicates are collected, averaged and recorded in Table 2. Figure 1 shows the Raman Spectra for the experimental runs mentioned in Table 2.

DoE Analysis Results
The data recorded in Table 2 are fed to Minitab software (trial package, version 17, Minitab Inc., Philadelphia, PA, USA, 2016) for DoE analysis.The analysis results of factors "A", "B" and "C" for two responses I and II, namely "I D /I G peak intensity ratio" and "I 2D /I G peak intensity ratio", respectively, are discussed as follows.
Figure 2 shows the main effect plots of factors "A", "B" and "C" for the two responses.Main effect plots represent the average change in output that follows from a change in the level [43].If a main effect plot is not horizontal, then different levels of the factor affect the response differently, and the steeper the slope of the plotted line, the greater the magnitude of the main effect that exists; whereas, if the plotted line is horizontal (parallel to the x-axis), there is no main effect present.In such a case, each level of the factor affects the response in the same way, and the response mean is the same across all factor levels [44].
The data recorded in Table 2 are fed to Minitab software (trial package, version 17, Minitab Inc., Philadelphia, PA, USA, 2016) for DoE analysis.The analysis results of factors "A", "B" and "C" for two responses I and II, namely "ID/IG peak intensity ratio" and "I2D/IG peak intensity ratio", respectively, are discussed as follows.
Figure 2 shows the main effect plots of factors "A", "B" and "C" for the two responses.Main effect plots represent the average change in output that follows from a change in the level [43].If a main effect plot is not horizontal, then different levels of the factor affect the response differently, and the steeper the slope of the plotted line, the greater the magnitude of the main effect that exists; whereas, if the plotted line is horizontal (parallel to the x-axis), there is no main effect present.In such a case, each level of the factor affects the response in the same way, and the response mean is the same across all factor levels [44].Figure 2 indicates that for both the responses of ID/IG and I2D/IG peak intensity ratios, the plotted line of all three factors "A", "B" and "C" are not parallel to the x-axis, and thus, they seem to affect the respective responses.
The interaction plots of the factors "A", "B" and "C" for the two responses in Table 2 are shown in Figure 3.An interaction plot indicates the dissimilarity in the response variation between the levels of one factor at different levels of other factors [45].In other words, the interaction between two factors shows that the effect of one factor on the response depends on the level chosen for its counterpart.Parallel plots are an indication of non-significant interaction, whereas non-parallel plots show interaction between the factors [46].Figure 2 indicates that for both the responses of I D /I G and I 2D /I G peak intensity ratios, the plotted line of all three factors "A", "B" and "C" are not parallel to the x-axis, and thus, they seem to affect the respective responses.
The interaction plots of the factors "A", "B" and "C" for the two responses in Table 2 are shown in Figure 3.An interaction plot indicates the dissimilarity in the response variation between the levels of one factor at different levels of other factors [45].In other words, the interaction between two factors shows that the effect of one factor on the response depends on the level chosen for its counterpart.Parallel plots are an indication of non-significant interaction, whereas non-parallel plots show interaction between the factors [46].
Figure 3 shows that all of the factors "A", "B" and "C" interact, producing interaction terms "AB" and "AC" for Response I and interaction terms "AC" and "C" for Response II, as their plots are non-parallel.In other words, the responses at different levels of annealing temperature, a-C layer thickness and annealing time chosen for our PVD-based graphene synthesis are dependent on the levels of each other.
The observed change in the two responses due to the factors individually and interactively could be due to the noise occurring in the experiments, and only statistically-significant changes can be considered as the main or interaction effects.Furthermore, the presence of a significant interaction term can mask the significance of the main effects.Hence, both the test for statistical significance and DoE plots need to be analyzed concurrently [32].Figure 3 shows that all of the factors "A", "B" and "C" interact, producing interaction terms "AB" and "AC" for Response I and interaction terms "AC" and "C" for Response II, as their plots are non-parallel.In other words, the responses at different levels of annealing temperature, a-C layer thickness and annealing time chosen for our PVD-based graphene synthesis are dependent on the levels of each other.
The observed change in the two responses due to the factors individually and interactively could be due to the noise occurring in the experiments, and only statistically-significant changes can be considered as the main or interaction effects.Furthermore, the presence of a significant interaction term can mask the significance of the main effects.Hence, both the test for statistical significance and DoE plots need to be analyzed concurrently [32].
In order to test for statistical significance, Analysis of Variance (ANOVA) statistics is employed, and the analysis results are presented in Table 3 and Table 4 for the two responses, respectively.
Table 3 indicates that the factors' interaction "AB" is the only statistically-significant term for Response I at the 10% significance level as the "p" value is less than 0.1.The presence of the main effect of the factors "A", "B", "C" and the factors' interaction term "AC" for Response I, as depicted by the plots in Figure 2 and Figure 3, respectively, are all found to be statistically insignificant.  4 indicates that the main effect of the factors "A" and "B" along with the factors' interaction term "BC" are all statistically-significance for the response II at the 10% significant level, as the "p" In order to test for statistical significance, Analysis of Variance (ANOVA) statistics is employed, and the analysis results are presented in Tables 3 and 4 for the two responses, respectively.Table 3 indicates that the factors' interaction "AB" is the only statistically-significant term for Response I at the 10% significance level as the "p" value is less than 0.1.The presence of the main effect of the factors "A", "B", "C" and the factors' interaction term "AC" for Response I, as depicted by the plots in Figures 2 and 3, respectively, are all found to be statistically insignificant.
Table 4 indicates that the main effect of the factors "A" and "B" along with the factors' interaction term "BC" are all statistically-significance for the response II at the 10% significant level, as the "p" value is less than 0.1.The presence of the main effect of the factor "C" and the factors' interaction term "AC" for Response II, as depicted by the plots in Figures 2 and 3, respectively, are all statistically insignificant.The interaction term "BC" is the most significant in affecting Response II, followed by the factor "A" and the factor "B".

Physical Reasons of the Statistically-Significant Terms
From the above ANOVA analysis, we have found the statistically-significant terms for Responses I and II, respectively.Let us now examine the physical reasons for these identified factors in affecting the two responses.
1.The presence of the significant interaction term "AB" mentioned in Table 3 can be explained with the aid of the following governing chemical reactions for graphene growth as given by Vlassiouk et al. [47].These reactions are shown in Scheme 1 in Figure 4 and are modified for our experiments since a-C is used as the carbon source instead of methane.

Physical Reasons of the Statistically-Significant Terms
From the above ANOVA analysis, we have found the statistically-significant terms for Responses I and II, respectively.Let us now examine the physical reasons for these identified factors in affecting the two responses.
1.The presence of the significant interaction term "AB" mentioned in Table 3 can be explained with the aid of the following governing chemical reactions for graphene growth as given by Vlassiouk et al. [47].These reactions are shown in Scheme 1 in Figure 4 and are modified for our experiments since a-C is used as the carbon source instead of methane.The formation of hydrocarbons responsible for graphene formation by the reaction of a-C and H2 at elevated temperatures as shown in Reaction (1) has been confirmed by Ji et al. [29], as well.
The rise in temperature causes thermo-mechanical stress [48] at the a-C/Cu interface due to the difference in the coefficient of thermal expansion for a-C and Cu.As a result, this thermo-mechanical stress causes the sample to warp downwards.This can be seen from the equivalent von-Mises stress analysis obtained from our Finite Element Analysis (FEA) simulation performed for the experimental test sample, as depicted in Figure 5.This FEA is performed using thermal and static-structural analysis in ANSYS 16.2 Academic Version (ANSYS, Inc., Canonsburg, PA, USA, 2015).The material properties used for FEA [49][50][51][52][53][54] are mentioned in Table S1 of the Supplementary Materials.Here, the resultant stress is given by Equation ( 1) and is popularly known as Stoney's equation [55,56].The formation of hydrocarbons responsible for graphene formation by the reaction of a-C and H 2 at elevated temperatures as shown in Reaction (1) has been confirmed by Ji et al. [29], as well.
The rise in temperature causes thermo-mechanical stress [48] at the a-C/Cu interface due to the difference in the coefficient of thermal expansion for a-C and Cu.As a result, this thermo-mechanical stress causes the sample to warp downwards.This can be seen from the equivalent von-Mises stress analysis obtained from our Finite Element Analysis (FEA) simulation performed for the experimental test sample, as depicted in Figure 5.This FEA is performed using thermal and static-structural analysis in ANSYS 16.2 Academic Version (ANSYS, Inc., Canonsburg, PA, USA, 2015).The material properties used for FEA [49][50][51][52][53][54] are mentioned in Table S1 of the Supplementary Materials.Here, the resultant stress is given by Equation ( 1) and is popularly known as Stoney's equation [55,56].
where ∆(1/R) = (1/R ´1/R 0 ), R 0 and R are the radius of the curvature of the samples before and after the thermo-mechanical stresses are developed, respectively, and "t" denotes the thickness of the layer in the sample with subscripts "s" and "f " referring to the substrate and film, respectively.Equation (1) shows that the film thickness and thermo-mechanical stress have an inverse relationship [55,56]; thus, the total stress will be higher for a thinner a-C/Cu interface, as shown in Figure 6.
where Δ(1/R) = (1/R − 1/R0), R0 and R are the radius of the curvature of the samples before and after the thermo-mechanical stresses are developed, respectively, and "t" denotes the thickness of the layer in the sample with subscripts "s" and "f" referring to the substrate and film, respectively.
Equation (1) shows that the film thickness and thermo-mechanical stress have an inverse relationship [55,56]; thus, the total stress will be higher for a thinner a-C/Cu interface, as shown in Figure 6.High stress at the interface and thermo-mechanical stress causing the downward warpage enhance carbon diffusion through Cu grain boundaries, resulting in the carbon reaching the surface of the Cu film.The von-Mises stress in the Cu film is tensile in nature (the values are believed to be high enough that they can either cause cracking in Cu film and/or enhancement of H2 diffusion in Cu film through its grain boundaries [57,58]).The enhancement of H2 diffusion will cause more carbon radical generation as per Reaction (1) in Figure 4. Furthermore, the a-C layer is found to be under compressive stress, which causes a "push" in the perpendicular direction towards the surface.Both of these mechanisms will enhance carbon diffusion through Cu boundaries.This renders the carbon radical generation Reaction (1) and graphene formation Reaction (2) favorable, and the defective sites, especially those at the edges, are repaired by the graphene formation, causing ID/IG ratio to fall.A where Δ(1/R) = (1/R − 1/R0), R0 and R are the radius of the curvature of the samples before and after the thermo-mechanical stresses are developed, respectively, and "t" denotes the thickness of the layer in the sample with subscripts "s" and "f" referring to the substrate and film, respectively.
Equation (1) shows that the film thickness and thermo-mechanical stress have an inverse relationship [55,56]; thus, the total stress will be higher for a thinner a-C/Cu interface, as shown in Figure 6.High stress at the interface and thermo-mechanical stress causing the downward warpage enhance carbon diffusion through Cu grain boundaries, resulting in the carbon reaching the surface of the Cu film.The von-Mises stress in the Cu film is tensile in nature (the values are believed to be high enough that they can either cause cracking in Cu film and/or enhancement of H2 diffusion in Cu film through its grain boundaries [57,58]).The enhancement of H2 diffusion will cause more carbon radical generation as per Reaction (1) in Figure 4. Furthermore, the a-C layer is found to be under compressive stress, which causes a "push" in the perpendicular direction towards the surface.Both of these mechanisms will enhance carbon diffusion through Cu boundaries.This renders the carbon radical generation Reaction (1) and graphene formation Reaction (2) favorable, and the defective sites, especially those at the edges, are repaired by the graphene formation, causing ID/IG ratio to fall.A High stress at the interface and thermo-mechanical stress causing the downward warpage enhance carbon diffusion through Cu grain boundaries, resulting in the carbon reaching the surface of the Cu film.The von-Mises stress in the Cu film is tensile in nature (the values are believed to be high enough that they can either cause cracking in Cu film and/or enhancement of H 2 diffusion in Cu film through its grain boundaries [57,58]).The enhancement of H 2 diffusion will cause more carbon radical generation as per Reaction (1) in Figure 4. Furthermore, the a-C layer is found to be under compressive stress, which causes a "push" in the perpendicular direction towards the surface.Both of these mechanisms will enhance carbon diffusion through Cu boundaries.This renders the carbon radical generation Reaction (1) and graphene formation Reaction (2) favorable, and the defective sites, especially those at the edges, are repaired by the graphene formation, causing I D /I G ratio to fall.A higher temperature will cause the stress to be higher, and it will also increase the rate of Reactions ( 1) and ( 2), thus rendering a larger reduction of the defects in the grown graphene.
The major difference between the diffusion of the solid carbon source and its gaseous counterpart is that in our method, the solid carbon source generates the carbon radical from beneath the Cu layer, and it diffuses through the Cu grain boundaries to stabilize on the surface as graphene.The thickness of the carbon layer becomes an important factor in such a case and influences the growth.In the gaseous carbon source, the carbon radicals are generated on the Cu surface, and there is no stress-enhanced generation of carbon radicals.
On the other hand, for a thicker a-C layer, the thermo-mechanical stress is lower, and Reaction (3) is more dominant.Under such circumstances, the rise in temperature will cause more etching of the synthesized graphene, resulting in an increase in the defective sites in the graphene layer due to the etching, and causing the I D /I G ratio to rise.
2. The presence of the significant interaction term "BC" mentioned in Table 4 can again be explained by the stress-related growth mechanism mentioned above.As mentioned earlier, for a thin a-C layer, the formation of graphene is enhanced, and the annealing time required is short, and vice versa for a thicker a-C layer.If the annealing time is sufficiently long, the response of the graphene growth will be independent of the a-C layer, as by then, even the case of the thick a-C layer would have graphene grown already.When the annealing time becomes short, there will be very little or even no graphene grown for a thick a-C layer in contrast to the case of a thin a-C layer.Hence, we can see the interaction of the factors "B" and "C" in Figure 3, which is statistically significant, as well.
The reasons for the statistical significance of the factors "A", i.e., the annealing temperature, and "B", i.e., the a-C layer thickness, as mentioned in Table 4, can be easily understood from the above explanation of the stress-related graphene growth mechanism and governing reactions shown in Figure 4. Factor "B" controls the amount of carbon supply, which is related to the variable stress at different levels of the factor "B", whereas the factor "A" influences the rates of Reaction ( 1), ( 2) and ( 3), as shown in Figure 4.

Optimal Level of DoE Factors
From the study of DoE, one can see that the three controlling factors do have a significant impact on the graphene growth and its quality.In order to find an optimal level of different factors for good quality graphene growth, i.e., low I D /I G peak intensity ratio (Response I) and different numbers of graphene layers, i.e., variable I 2D /I G peak intensity ratio (Response II), contour plot analysis is performed based on the results from DoE, as shown in Figures 7 and 8, respectively.Figure 8 depicts the contour and surface plots for Response II and shows the required levels of factors affecting the number of synthesized graphene layers.They suggest that a high level of the factor "A" and a low level of the factors "B" and "C" are the required levels for producing graphene with a high value of the I2D/IG peak ratio (>1.5), which is highly desirable for single-layer or few-layer graphene [40,41].In other words, a treatment combination of higher annealing temperature, thinner a-C layer and smaller annealing time will result in a higher I2D/IG peak intensity ratio, and hence, mono-layer/few-layer graphene will be synthesized.
On the other hand, if multi-layer graphene is desired, then there is a large region for the choice of levels in each of the contour plots shown in Figure 8.These regions represent the I2D/IG peak intensity ratio values below one, which is a signature for multi-layer graphene [40,41].
The trends suggested by Figures 7 and 8 are instrumental in bringing out a clear relation between the choice of levels for the factors and the quality/number of graphene layers, and they are summarized in Table 5.
For the case of mono-/few-layer graphene, the annealing time is required to be optimized as it affects the two responses in opposite directions, as shown in Table 5.Long annealing time favors low defects, but smaller annealing time favors mono-/few-layer graphene synthesis.Hence, the response surface method [37,45] must be employed for optimizing growth parameters, which we will present in the future work.These plots show the required levels of factors in affecting the quality of graphene.They indicate that a high level of the factor "A", a low level of the factor "B" and a high level of the factor "C" are the required levels for producing graphene with a low value of the I D /I G peak intensity ratio, which is highly desirable for good quality graphene [40,41].In other words, a combination of longer annealing time, a thinner a-C layer and a higher annealing temperature will result in a lower I D /I G peak intensity ratio and, hence, graphene with lesser defects.

Respective
Figure 8 depicts the contour and surface plots for Response II and shows the required levels of factors affecting the number of synthesized graphene layers.They suggest that a high level of the factor "A" and a low level of the factors "B" and "C" are the required levels for producing graphene with a high value of the I 2D /I G peak intensity ratio (>1.5), which is highly desirable for single-layer or few-layer graphene [40,41].In other words, a treatment combination of higher annealing temperature, thinner a-C layer and smaller annealing time will result in a higher I 2D /I G peak intensity ratio, and hence, mono-layer/few-layer graphene will be synthesized.
On the other hand, if multi-layer graphene is desired, then there is a large region for the choice of levels in each of the contour plots shown in Figure 8.These regions represent the I 2D /I G peak intensity ratio values below one, which is a signature for multi-layer graphene [40,41].
The trends suggested by Figures 7 and 8 are instrumental in bringing out a clear relation between the choice of levels for the factors and the quality/number of graphene layers, and they are summarized in Table 5.
Multi-/many-layer graphene with low defects Low Response I and low However, the case for multi-layer graphene with low defects is not straight-forward.The region for multi-layer graphene represented in Figure 8 is large, but the region for low defects in Figure 7 is smaller.In order to have a clearer picture, Figures 7 and 8 are superimposed, as shown in Figure 9.The coinciding regions (depicted by dotted ovals) in Figure 9 indicate that a high level of the factors "A" and "C" and a low level of the factor "B" are the required levels for many-layer graphene with low defects.These regions suggest that a treatment combination of higher annealing temperature, thinner a-C layer and longer annealing time will result in a low I2D/IG peak intensity ratio (close to one or less) and a low ID/IG peak intensity ratio; hence, multi-layer graphene [40,41] with less defects will be synthesized.For the case of mono-/few-layer graphene, the annealing time is required to be optimized as it affects the two responses in opposite directions, as shown in Table 5.Long annealing time favors low defects, but smaller annealing time favors mono-/few-layer graphene synthesis.Hence, the response surface method [37,45] must be employed for optimizing growth parameters, which we will present in the future work.However, the case for multi-layer graphene with low defects is not straight-forward.The region for multi-layer graphene represented in Figure 8 is large, but the region for low defects in Figure 7 is smaller.In order to have a clearer picture, Figures 7 and 8 are superimposed, as shown in Figure 9.

Experimental Verification
The coinciding regions (depicted by dotted ovals) in Figure 9 indicate that a high level of the factors "A" and "C" and a low level of the factor "B" are the required levels for many-layer graphene with low defects.These regions suggest that a treatment combination of higher annealing temperature, thinner a-C layer and longer annealing time will result in a low I 2D /I G peak intensity ratio (close to one or less) and a low I D /I G peak intensity ratio; hence, multi-layer graphene [40,41] with less defects will be synthesized.

Experimental Verification
In order to further validate the above statistical results besides the physical explanations, experiments are performed.Experimental samples with an a-C layer thickness of 12 nm are annealed at 1020 ˝C, and the annealing time is varied across multiple levels from 5 to 50 min.
Figure 10 shows the plots of the I D /I G and I 2D /I G peak intensity ratios vs. annealing time, for samples with thin a-C layer, annealed at a high annealing temperature.In order to further validate the above statistical results besides the physical explanations, experiments are performed.Experimental samples with an a-C layer thickness of 12 nm are annealed at 1020 °C, and the annealing time is varied across multiple levels from 5 to 50 min.
Figure 10 shows the plots of the ID/IG and I2D/IG peak intensity ratios vs. annealing time, for samples with thin a-C layer, annealed at a high annealing temperature.The plots in Figure 10 clearly show that the I2D/IG peak intensity ratio is larger for shorter annealing times (depicted by a dotted circle), whereas the ID/IG peak intensity ratio is lower for longer annealing times (depicted by a dashed oval), when the experiments are conducted at a high annealing temperature for samples with a thin a-C layer.This result is in clear agreement with Table 5.Therefore, in order to achieve a large I2D/IG peak intensity ratio and a small ID/IG peak intensity ratio simultaneously, the annealing time has to be optimized.
On the other hand, longer annealing times with the same parameters for the annealing temperature and a-C thickness will result in multi-layer graphene with low defects (depicted by dotted rectangles in Figure 10), in agreement with Table 5.Therefore, we can see that the type of The plots in Figure 10 clearly show that the I 2D /I G peak intensity ratio is larger for shorter annealing times (depicted by a dotted circle), whereas the I D /I G peak intensity ratio is lower for longer annealing times (depicted by a dashed oval), when the experiments are conducted at a high annealing temperature for samples with a thin a-C layer.This result is in clear agreement with Table 5.
Therefore, in order to achieve a large I 2D /I G peak intensity ratio and a small I D /I G peak intensity ratio simultaneously, the annealing time has to be optimized.
On the other hand, longer annealing times with the same parameters for the annealing temperature and a-C thickness will result in multi-layer graphene with low defects (depicted by dotted rectangles in Figure 10), in agreement with Table 5.Therefore, we can see that the type of desired graphene with respect to the quality and number of layers can be adjusted by choosing suitable levels obtained from the DoE analysis.

Conclusions
DoE has been successfully implemented to identify important factors that control the graphene synthesis using a-C as the carbon source and Cu as the catalyst.Full-factorial DoE analysis reveals that the interaction of annealing temperature and a-C layer thickness is statistically significant for the I D /I G peak intensity ratio and the quality of grown graphene.On the other hand, annealing temperature, a-C layer thickness along with the interaction of the a-C layer thickness and annealing time are statistically significant for the I 2D /I G peak intensity ratio and the number of layers of graphene.The physical reasons for the factors being significant for different aspects of the grown graphene are also explained, showing the consistency of the DoE method and the underlying process physics and chemistry.
Contour plots suggest that a higher annealing temperature and a thinner a-C layer are favorable for producing graphene with a smaller value of the I D /I G peak intensity ratio and a larger value of the I 2D /I G peak intensity ratio.However, the annealing time has to be optimized in this case, in order to synthesize mono-layer/few-layer graphene with low defects.On the other hand, higher annealing temperature, thinner a-C layer and longer annealing time will result in multi-layer graphene with low defects.This is verified by our experimental data.
In short, the application of DoE analysis on this type of graphene synthesis not only helps to determine the key controlling factors, it can also indicate a suitable combination of the levels for different factors for obtaining graphene of different qualities and different numbers of layers.

Supplementary Materials:
The following are available online at www.mdpi.com/2076-3417/6/7/204/s1. Figure S1: Sample preparation using PVD (Physical Vapor Deposition) method; 2 sets of samples each having 36 nm and 12 nm thick a-C (amorphous carbon) thin films are prepared.Figure S2: Thermal Annealing after sample preparation; Annealing is carried out at 2 levels of temperatures that are 820 ˝C /1020 ˝C and 2 levels of duration that are 10 min/50 min.Table S1: Thermo-mechanical properties of the materials used for Finite Element Analysis (FEA).

Figure 2 .
Figure 2. Main effect plots of the factors "A", "B" and "C" for Response I (a) and Response II (b).All of the factors seem to show the main effect on both responses.

Figure 2 .
Figure 2. Main effect plots of the factors "A", "B" and "C" for Response I (a) and Response II (b).All of the factors seem to show the main effect on both responses.

Figure 3 .
Figure 3. Interaction plots of the factors "A", "B" and "C" for Response I (a) and Response II (b); solid circles indicate the interaction points.Factor interaction terms "AB" and "AC" are visible for Response I, while factor interaction terms "AC" and "BC" are visible for Response II.

Figure 3 .
Figure 3. Interaction plots of the factors "A", "B" and "C" for Response I (a) and Response II (b);solid circles indicate the interaction points.Factor interaction terms "AB" and "AC" are visible for Response I, while factor interaction terms "AC" and "BC" are visible for Response II.

Figure 4 .
Figure 4. Scheme 1: Graphene growth mechanism's governing reactions; here C represents the carbon atom, Gr represents graphene, CH*/C* represent hydrogen bonded-carbon radical/carbon radical, H* represents the hydrogen radical and ∆H represents heat (permission to modify the scheme has been obtained).

Figure 4 .
Figure 4. Scheme 1: Graphene growth mechanism's governing reactions; here C represents the carbon atom, Gr represents graphene, CH*/C* represent hydrogen bonded-carbon radical/carbon radical, H* represents the hydrogen radical and ∆H represents heat (permission to modify the scheme has been obtained).

Figure 5 .
Figure 5. ANSYS simulation of thermo-mechanical stress (equivalent von-Mises stress) at the a-C and Cu interface at 1020 °C.

Figure 6 .
Figure 6.ANSYS simulation results of maximum principal stress distributions in the three samples with different a-C layer thickness at 1020 °C.The area of the maximum stress distribution is the largest in the sample with the thinnest a-C layer.

Figure 5 .
Figure 5. ANSYS simulation of thermo-mechanical stress (equivalent von-Mises stress) at the a-C and Cu interface at 1020 ˝C.

Figure 5 .
Figure 5. ANSYS simulation of thermo-mechanical stress (equivalent von-Mises stress) at the a-C and Cu interface at 1020 °C.

Figure 6 .
Figure 6.ANSYS simulation results of maximum principal stress distributions in the three samples with different a-C layer thickness at 1020 °C.The area of the maximum stress distribution is the largest in the sample with the thinnest a-C layer.

Figure 6 .
Figure 6.ANSYS simulation results of maximum principal stress distributions in the three samples with different a-C layer thickness at 1020 ˝C.The area of the maximum stress distribution is the largest in the sample with the thinnest a-C layer.

Figure 7 .
Figure 7. Contour plots for Response I with respect to the factors "A", "B" and "C"; insets show the surface plots.High levels of the factor "A" and the factor "C" and a low level of the factor "B" indicate a low value of the ID/IG peak intensity ratio.

Figure 7
Figure 7 depicts the contour and surface plots for Response I, wherein any two factors are displayed on the x-and y-scales, and the response (on the z scale) variable is represented by a contour and smooth surface, respectively.These plots show the required levels of factors in affecting the quality of graphene.They indicate

Figure 7 .
Figure 7. Contour plots for Response I with respect to the factors "A", "B" and "C"; insets show the surface plots.High levels of the factor "A" and the factor "C" and a low level of the factor "B" indicate a low value of the I D /I G peak intensity ratio.

Figure 8 .
Figure 8. Contour plots for Response II with respect to the factors "A", "B" and "C"; the insets show the surface plots.The area within dotted boundaries represents multi-layer graphene.

Figure 8 .
Figure 8. Contour plots for Response II with respect to the factors "A", "B" and "C"; the insets show the surface plots.The area within dotted boundaries represents multi-layer graphene.

Figure 7
Figure 7 depicts the contour and surface plots for Response I, wherein any two factors are displayed on the xand y-scales, and the response (on the z scale) variable is represented by a contour and smooth surface, respectively.These plots show the required levels of factors in affecting the quality of graphene.They indicate that a high level of the factor "A", a low level of the factor "B" and a high level of the factor "C" are the required levels for producing graphene with a low value of the I D /I G peak intensity ratio, which is highly desirable for good quality graphene[40,41].In other words, a combination of longer annealing time, a thinner a-C layer and a higher annealing temperature will result in a lower I D /I G peak intensity ratio and, hence, graphene with lesser defects.Figure8depicts the contour and surface plots for Response II and shows the required levels of factors affecting the number of synthesized graphene layers.They suggest that a high level of the factor "A" and a low level of the factors "B" and "C" are the required levels for producing graphene with a high value of the I 2D /I G peak intensity ratio (>1.5), which is highly desirable for single-layer or few-layer graphene[40,41].In other words, a treatment combination of higher annealing temperature, thinner a-C layer and smaller annealing time will result in a higher I 2D /I G peak intensity ratio, and hence, mono-layer/few-layer graphene will be synthesized.On the other hand, if multi-layer graphene is desired, then there is a large region for the choice of levels in each of the contour plots shown in Figure8.These regions represent the I 2D /I G peak intensity ratio values below one, which is a signature for multi-layer graphene[40,41].The trends suggested by Figures7 and 8are instrumental in bringing out a clear relation between the choice of levels for the factors and the quality/number of graphene layers, and they are summarized in Table5.

Figure 9 .
Figure 9. Superimposition of contour plots for Responses I and II, as shown in Figures 7 and 8, respectively; dotted ovals show the region of coincidence.

Figure 9 .
Figure 9. Superimposition of contour plots for Responses I and II, as shown in Figures7 and 8, respectively; dotted ovals show the region of coincidence.

Figure 10 .
Figure10.Plots of the ID/IG and I2D/IG peak intensity ratios vs. annealing time, for samples with a thin a-C layer annealed at 1020 °C.The dotted circle shows the region with a larger I2D/IG peak intensity ratio at a shorter annealing times.Dotted rectangles show the region with a smaller I2D/IG peak intensity ratio, and the dashed oval shows the region with a smaller ID/IG peak intensity ratio, both at longer annealing times.

Figure 10 .
Figure 10.Plots of the I D /I G and I 2D /I G peak intensity ratios vs. annealing time, for samples with a thin a-C layer annealed at 1020 ˝C.The dotted circle shows the region with a larger I 2D /I G peak intensity ratio at a shorter annealing times.Dotted rectangles show the region with a smaller I 2D /I G peak intensity ratio, and the dashed oval shows the region with a smaller I D /I G peak intensity ratio, both at longer annealing times.

Table 3 .
Analysis of variance for ID/IG using adjusted SS (Sum of Squares) for the tests.The text marked in red represents the statistically-significant term.

Table 3 .
Analysis of variance for I D /I G using adjusted SS (Sum of Squares) for the tests.The text marked in red represents the statistically-significant term.

Table 4 .
Analysis of variance for I 2D /I G using adjusted SS for the tests.The texts marked in red represent statistically-significant terms.

Table 4 .
Analysis of variance for I2D/IG using adjusted SS for the tests.The texts marked in red represent statistically-significant terms.

Table 5 .
Relation between the choice of levels and the quality/number of graphene layers.

Table 5 .
Relation between the choice of levels and the quality/number of graphene layers.