Machine intelligence supports the full chain of 2D dendrite synthesis

1 Machine intelligence supports the full chain of 2D dendrite synthesis Wenqiang Huang 1,2,3, † , Susu Fang 2, † , *, Xuhang Gu 3, † , Shen’ao Xue 1,2 , Huanhuan Xing 2 , Junjie Jiang 1,2 , Junying Zhang 1,2 , Shen Zhou 4 , Zheng Luo 2 , Jin Zhang 3 , Fangping Ouyang 1,5,6 *, and Shanshan Wang 2,3 * 1 School of Physics, Hunan Key Laboratory for Super-Microstructure and Ultrafast Process, and Hunan Key Laboratory of Nanophotonics and Devices, Central South University, Changsha 410083, China 2 College of Aerospace Science and Engineering, National University of Defense Technology, Changsha 410000, China 3 School of Advanced Materials, Guangdong Provincial Key Laboratory of Nano-Micro Materials Research, Peking University Shenzhen Graduate School, Shenzhen 518055, China 4 College of Science, National University of Defense Technology, Changsha 410000, China 5 School of Physics and Technology, and Xinjiang Key Laboratory of Solid-State Physics and Devices, Xinjiang University, Urumqi 830046, China 6 State Key Laboratory of Powder Metallurgy, and Powder Metallurgy Research Institute, Central South University, Changsha 410083, China † These authors contributed equally. *Corresponding authors: S.F. (fangsusugfkd@nudt.edu.cn); F.O. (oyfp@csu.edu.cn); S.W. (wangshanshan08@nudt.edu.cn) Exemplified by the chemical vapor deposition growth of two-dimensional dendrites, which has potential applications in catalysis and presents a parameter-intensive, data-scarce and reaction process-complex model problem, we devise a machine intelligence-empowered framework for the full chain support of material synthesis, encompassing rapid process optimization, accurate customized synthesis, and comprehensive mechanism deciphering. First, active learning is integrated into the experimental workflow, identifying an optimal recipe for the growth of highly-branched, electrocatalytically-active ReSe 2 dendrites through 60 experiments (4 iterations), which account for less than 1.3% of the numerous possible parameter combinations. Then, a prediction accuracy- guided data augmentation strategy is developed combined with a tree-based machine learning (ML) algorithm, unveiling a non-linear correlation between 5 process variables and fractal dimension (D F ) of ReSe 2 dendrites with only 9 experiment additions, which guides the synthesis of various user- defined D F . Finally, we construct a data–knowledge dual-driven mechanism model by integration of cross-scale characterizations, interpretable ML models, and domain knowledge in thermodynamics and kinetics, unraveling synergistic contributions of multiple process parameters to the product morphology. This work demonstrates the ML potential to transform the research paradigm and is adaptable to broader material synthesis. 2 Introduction Dendrites, a prevalent microstructure in solidification, play significant roles in capacitors 1 , rechargeable batteries 2–4 , and sensors 5 . Its two-dimensional (2D) counterparts are featured with a high ratio of surface atoms, exotic edge states, and superior mechanical flexibility, displaying unique prospects in catalysis 6,7 , non-linear optics 8,9 , and dielectric memristors 10 . Chemical vapor deposition (CVD), which has achieved great success in the structure-controlled synthesis of various 2D single crystals 11–16 , homo/heterostructures 17–20 , and superlattices 21–24 , is regarded as the most promising approach for the large- scale, high-quality preparation of 2D dendrites 25,26 . However, the CVD approach involves numerous synthesis parameters (e.g., growth temperature, precursor type/dosage, reaction time, and carrier gas), which constitute a vast search space with complex growth mechanisms lying behind, making the traditional trial-and-error method, such as the one factor at one time (OFAT), suffer from low efficiency, high cost, and poor transferability. Fortunately, machine learning (ML), as a mathematical tool adept at solving high-dimensional problems, is transforming the research logic of materials science 27–31 and holds great promise in the full- chain support to synthesis, specifically in achieving three long-sought goals: (i) Efficient process optimization. The ML approaches can find the optimal process corresponding to a material’s extreme properties in a huge parameter space with as few experiments as possible. (i) Customized synthesis on demand. Complex functional relationships between multiple process/component variables and material properties can be fitted via ML methods, so that researchers can easily obtain the recipe for the desired property values. (i) Multifactorial mechanism deciphering. The impact of multiple parameters on the product performance can be evaluated simultaneously and quantitatively via interpretable ML approaches (e.g., SHapley Additive exPlanations, correlation analysis), so that a panorama of growth mechanisms involving multifactor importance analysis and synergies between variables can be unveiled. Increasing efforts have been devoted to ML-assisted material synthesis in recent years. Notable studies include constructing ML-guided frameworks, such as active learning loops, supervised-learning workflows and transformer-based language models, to design all-natural plastic substitutes with diversified user- desired properties 32 , to discover multi-metallic perovskite oxides and alloys with record catalytic activities 33–35 , to optimize the manufacturing process for device efficiency improvement of solar cells 36 , and to customize nanomaterial synthesis with controllable morphology, density, and optical wavelengths 37– 41 , etc. However, several limitations remain in prior work. First, there is a lack of research cases where machine intelligence is utilized to empower the material growth from process optimization to customized synthesis to mechanism deciphering, especially for small data, which may weaken the academic community’s perception of ML’s capability in the full-chain support for material synthesis. Second, in terms of mechanism unraveling, the results given by the ML models are sometimes difficult to fully integrate with the material characterizations and expert knowledge. Most interpretations involve only single variables, with a dearth of quantitative analysis of multiple factors, hindering the establishment of a growth mechanism that is truly interpretable, understandable, and comprehensive from the dual perspective of human and machine intelligence. In this work, we select the CVD growth of 2D ReSe 2 dendrites as a model problem, in which ReSe 2 belongs to the low-symmetry triclinic crystal system with a more complex growth behavior that has rarely been studied compared with its high-symmetry counterparts (e.g., graphene, MoS 2 , etc.), thus providing a 3 superior platform to validate machine intelligence empowerment in the whole chain of material synthesis. First, an active learning loop is designed to find the optimal recipe (involving 5 process variables) for the synthesis of 2D dendrites with a high fractal dimension (D F ) of 1.71 using only one-week experiments. Second, a nonlinear mapping between multiple growth parameters and dendritic morphology of 2D ReSe 2 is constructed using a gradient boosting framework with a prediction accuracy-guided data augmentation strategy, yielding a predictive capability of R 2 =0.86 across the entire design space with a small budget of ˂ 70 experiments. Finally, a multifactorial growth mechanism is established by jointly applying cross-scale microscopic/spectroscopic characterizations, domain knowledge, and interpretable ML approaches, thus unraveling synergistic effects of the rhenium source temperature, concentration, and the substrate type on the 2D ReSe 2 dendrite growth, which aligns with both human cognition and machine intelligence conclusions. Results Overview of the ML framework Our ML framework comprises three modules: process optimization, customized synthesis, and mechanism deciphering, demonstrating the full-chain support of machine intelligence in the CVD synthesis of 2D ReSe 2 dendrites. Fig.1 Scheme of the ML framework depicting the full-chain support of machine intelligence in the CVD preparation of 2D ReSe 2 dendrites. Three modules are included, i.e., (a) process optimization, (b) customized synthesis, and (c) mechanism deciphering. Process optimization. This module is designed to quickly retrieve the optimal process conditions for the CVD growth of highly branched 2D ReSe 2 from a vast parameter space using as little experimental cost as possible (Fig. 1a). First, the input variables (x i ) and the target output (y) for the active learning model are determined so that an initial dataset for iterative optimization can be established. Five key process parameters are identified as x i based on empirical knowledge, which are the heating temperature of 4 the rhenium source (T Re ), the heating temperature of the selenium powder (T Se ), the concentration of the rhenium source (c Re ), the flow rate of the hydrogen gas (f H2 ), and the substate type (sub.). Their ranges are decided by a group of binary classification experiments (can grow vs. cannot grow) to ensure that the probability of growing 2D ReSe 2 within this sample space (involving 4752 condition possibilities) is high (Supplementary Tables 1-3). The fractal dimension (D F ), a metric quantifying geometric branching fineness, is introduced as 푦 to evaluate the dendritic degree of ReSe 2 crystals, whose value can be programmatically extracted from scanning electron microscopy (SEM) images using a standard box-counting method (Supplementary Fig. 1). Twenty different process conditions are selected via Latin hypercube sampling (LHS) in order to generate a near-random sample of the input variable values from a multidimensional distribution (Supplementary Figs. 2 and 3). Their corresponding average D F of the as-grown crystals are achieved by CVD experiments, preparing an initial dataset for the first round of model training. Second, a Bayesian optimization (BO) strategy is utilized for active learning with Gaussian process regression (GPR) as the surrogate model. This is mainly due to two reasons: i) there is a lack of any theoretical or empirical formulas to describe the explicit mathematical relationship between the CVD process variables and the crystal morphology, while BO-GPR excels at solving such a black-box function approximation problem; i) BO-GPR is inherently capable of quantifying uncertainty over predictions. This nature allows a good balance between exploration and exploitation when sampling in new iterative cycles via acquisition functions, thus reaching a global optimum with a small experimental budget 34,42,43 . Third, Max-value entropy search (MES) is employed as the acquisition function to suggest new sampling points in the next round after training the BO-GPR model with the initial dataset. This acquisition function leverages an information-theoretic perspective to maximize global information gain (accomplished by minimizing the information entropy of the objective function’s maximum value) with minimum sampling cost 44 . Finally, CVD experiments corresponding to the process conditions suggested by MES are conducted to achieve new D F which update the BO-GPR model iteratively in the active learning loop until the improvement of D F is close to saturation. Customized synthesis. This module aims to construct a mapping between the five process variables and the dendritic morphology of ReSe 2 , enabling easy determination of process conditions corresponding to various user-defined D F (Fig. 1b). Firstly, we select an ML model that is best suited to perform this regression task by comparing how well different models fit the data generated in Module 1. Various models with distinct structures are considered, including kernel-based, network-based, instance-based, and tree-based ones. Notably, the optimization purpose of Module 1 leads to an uneven distribution of the accumulated data at this stage, which is skewed toward the high D F regions. This deficiency results in a hardly satisfactory performance even for the best regression model, if solely using data achieved in Module 1 for training, since our target in Module 2 is to construct an accurate mapping between process variables and D F across a wide value range. Therefore, secondly, we develop a D F prediction accuracy-guided data augmentation strategy. Such discrepancy-based method supplements data for model tuning around the current data points with the poorest D F prediction accuracy, which optimizes the worst-fitting regions with minimum addition of experiments. It is worth noting that we define a PA score, the square of the difference between the model-predicted and the experimental D F values, to quantify the D F prediction accuracy of every sampling point, which displays high sensitivity to data points with large prediction errors (see 5 equation (2) and more details in Fig. 3). After a few rounds of data augmentation and model optimization, a global mapping with satisfactory prediction performance can be attained. Mechanism deciphering. This module is to unveil a multifactorial growth mechanism by integrating interpretable ML models, cross-scale spectroscopic/microscopic characterizations, and researchers’ domain knowledge in thermodynamics and kinetics (Fig. 1c). We employ a feature-based interpretability method, called SHAP analysis, to quantify the importance of multiple process variables to D F (SAHP values) and evaluate synergistic effects between features (SHAP interaction values) 45–47 . At the same time, a series of characterization techniques spanning from atomic to micrometer scales are performed, including Raman spectroscopy, annular dark-field scanning transmission electron microscopy (ADF-STEM), SEM, etc., to uncover morphology, orientations, and crystallinity of the as-grown 2D ReSe 2 . Finally, the quantitative explanation from the ML algorithms and the qualitative understanding from material characterization and expertise are integrated to establish a data–knowledge dual-driven growth mechanism, which is comprehensive, precise, and understandable from both perspectives of human and machine intelligence. Active learning-guided CVD process optimization The active learning loop achieves prominent D F improvement by a total of 60 experiments, which occupied <1.3% of the 4752 possible combinations in a grid search of the CVD parametric space. It begins with an initial dataset of 20 samples and undergoes 4 iterative cycles, each comprising 10 experiments. Fig. 2a is a box plot showcasing the progressive enhancement of D F as the iterative cycle increases. We chose the median (hollow squares in the middle of boxes) to reflect the central tendency of D F in each cycle, which is independent of data distribution and outliers compared with the mean. The D F medians rise from 1.36 (iteration 0, initial dataset) to 1.61 (iteration 4), which attains a 69.4% boost in the feasible range (D F ∈ (1,2), see Methods, “Range of D F ”) and is intuitively manifested by increased roughness and intricate branches of the crystal edges via SEM images (top insets). In the final iteration, 2D ReSe 2 crystals with D F up to 1.71 are obtained, exhibiting snowflake-like self-similar configurations with three-fold symmetry (inset, bottom right). Furthermore, the elevation rate of D F with iterations exhibits a rapid jump at the early optimization stage, yielding a 47.2% enhancement in iteration 1, and gradually reaches saturation with a dramatic decay of the rising speed by more than an order of magnitude to <4% in iteration 4, which provides a clear decision-making indication for terminating the active learning loop. Interestingly, distinct from the monotonic increase of the D F median, the distribution of D F with iterations shows an “accordion”-style alternating expansion and contraction, as revealed by the widths of boxes in Fig. 2a. This suggests that the monotonic rise of D F may be achieved through a trade-off between exploration (sampling in regions of high uncertainty to maximize information gain) and exploitation (sampling in regions with optimum predicted D F to seek the largest objective reward) via the MES acquisition function during the iterative BO-GPR model training. To uncover the optimization process, we drew 3D contour plots, which involve three key parameters (i.e., T Re , c Re , and T Se ) that contribute the most to D F as the three coordinate axes (Supplementary Fig. 4) and the predicted D F /uncertainty represented by the color of the contour map (Fig. 2b). When the BO-GPR model is trained using the initial 20 data points randomly selected by LHS, the temperature plane with T Re ≈ 620°C displays superior predicted D F with low uncertainty. It leads the algorithm to concentrate all the sampling points on this T Re plane in iteration 1 (exploitation strategy) but with a broad dispersion along the feasible range of T Se in order to balance exploration (red dots and double-headed arrow in iteration 0→1, 6 Supplementary Fig. 5). Since the algorithm adopts an exploitation strategy for the feature of T Re , which has a dominant contribution to the model output, it brings about a fast increase of D F with data distribution prominently narrowed compared with iteration 0 (magenta box in Fig. 2a). In iteration 2, the model chooses to scatter the sampling points in the c Re dimension and prefers regions of higher uncertainty (red dots and double-headed arrow in iteration 1→2, Supplementary Fig. 6). Such an exploration-dominated strategy results in not only an obvious decrease in the D F improvement rate but also an increase in the D F dispersion compared to iteration 1 (green box in Fig. 2a). After a large-scale exploration of T Se and c Re in the first two iterations, sampling in iteration 3 focuses on a small range of T Se ≈ 280°C and c Re ≈ 0.1 mol/L with relatively high predicted D F and uncertainty on the T Re ≈ 620°C plane (red circle in iteration 2→3, Supplementary Fig. 7). Such evolution of the sampling distribution from a large-range linear dispersion (iteration 2) to a small-range point-like concentration (iteration 3) leads to a drastically narrowed D F distribution (yellow box in Fig. 2a). In the last iteration, the sampling points remain concentrated in a small region with a slight adjustment of T Se and c Re on the T Re ≈ 620°C plane in order to maximize the predicted D F and minimize the uncertainty (red circle in iteration 3→4, Supplementary Fig. 8). This exploitation strategy leads to a minor enhancement of D F but with a remarkable widening of data distribution compared to iteration 3 (gray box in Fig. 2a), which signifies that BO is approaching its limit. Fig. 2. Active learning-guided optimization of CVD process parameters for 2D ReSe 2 dendrites. a, Box plot showing the improvement of D F with iterative cycles. Boxes: 10 th -90 th percentiles; whiskers: 3×interpercentile range; 7 hollow squares: medians. Insets (top) are false-color SEM images of 2D ReSe 2 dendrites corresponding to median D F in iterations 0 to 4. Inset (right bottom) displays a typical SEM image of the ReSe 2 dendrites with a maximum D F = 1.71±0.03 achieved in the final iteration (black star). b, Visualization of the BO process by exhibiting the three- dimensional CVD parameter-D F plot and the sampling strategies in every iteration s→(s+1) (s=0,1,2,3). CVD parameters: T Re , c Re , and T Se . Color bar (top panels): predicted average D F . Color bar (bottom panels): uncertainty of the predicted D F . The contour maps and gray dots represent the CVD parameter-D F relation and the sampling points in iteration s, while the red dots represent the sampling points in iteration (s+1). c, Evolution of five process parameters’ values across four iterations. The initial dataset is called iteration 0 with a size of 20 experiments, represented by blue dots. 10 data points are appended in each new iteration, represented by magenta, green, yellow, and gray dots, respectively. d, SHAP interaction value matrix in the logarithmic form based on the BO-GPR model trained by the initial 20 data points in iteration 0. e, Histogram of I scores for the five CVD process parameters in 4 iterations. f, Cathodic linear sweep voltammetry (LSV) curves of ReSe 2 with D F of 1.40±0.02, 1.59±0.02, and 1.71±0.03. g, The corresponding Tafel slopes. The visualization of the active learning process reveals that the algorithm employs different sampling strategies (exploration vs. exploitation) for different input features. Such feature-dependent sampling behavior can be unveiled more clearly in Fig. 2c. Parameters of sub. and T Re converge quickly to c-Al 2 O 3 and ≈620°C, respectively, as soon as iteration 0 finishes, which corresponds to a continuous employment of the exploitation-based sampling on these two features throughout the BO process. c Re experiences a large-scale sampling in the initial rounds, followed by a gradual convergence to ≈ 0.75 mol/L from iteration 3, indicating the application of an “exploration first, exploitation later” sampling strategy. Both T Se and f H2 exhibit non-convergent behavior as their values change round by round even until iteration 4, implying the dominance of exploration-based sampling in these two variable dimensions. This feature-related convergence behavior suggests that the independence level of different features’ impact on D F varies. Process parameters that converge quickly mean that their values corresponding to the optimum D F remain unchanged regardless of other parameters’ values, whereas the variables that are difficult to converge imply that their values for the optimum D F are heavily influenced by other features. We leveraged SHAP feature value analysis from collaborative game theory and proposed a measure, called I score, to quantify the independence level of each input variable’s impact on the model output. I score = |SHAP main effect x i | ∑ |SHAP interaction x i -x j | i≠j = |SHAP interaction x i -x i | ∑ |SHAP interaction x i -x j | i≠j (1) where SHAP main effect x i represents the SHAP main effect value of feature x i , indicating the contribution of x i to the model output when acting alone. It is equal to the SHAP interaction value of x i with itself, denoted as SHAP interaction x i -x i . SHAP interaction x i -x j represents the SHAP interaction value of two different features x i and x j , evaluating the synergistic contribution of x i and x j to the model output (i≠j) (See Methods, “Quantification of the independence level of feature’s impact on D F ”). The higher the I score is, the stronger the independence of this feature’s impact on D F becomes. The absolute values of SHAP interactions for the five CVD parameters after each iteration are displayed in logarithmically formed matrices (Fig. 2d and Supplementary Fig. 9) with the corresponding I scores calculated in Fig. 2e. The I scores of T Re and sub. remain the top two highest among the five parameters in all iterations, which indicates their strong independence level of impact on D F and is consistent with their fast convergence from iteration 1 in Fig. 2c. The I scores of c Re are situated in the middle in most iterations, while those of T Se and f H2 are at the bottom, in line with the moderate convergence rate of c Re and non-convergence behaviors of T Se and f H2 in active learning. It is noteworthy that the I score of f H2 reaches a relatively high value in iteration 0 with a 8 drastic decrease in the following three iterations, which coincides with the short-term stability of this feature’s value in iteration 1 followed by three rounds of continuous fluctuations in BO. Taking advantage of an abundance of nanoscale tips and unsaturated edge atoms, 2D ReSe 2 dendrites with high D F are expected to exhibit superior hydrogen evolution reaction (HER) performance, which was evaluated using a three-electrode system in 0.5 M H 2 SO 4 electrolyte (Supplementary Fig. 10). Cathodic linear sweep voltammetry (LSV) curves show that ReSe 2 with a D F of 1.71±0.03 exhibits an overpotential of only 195 mV at a current density of 10 mA cm -2 , representing a reduction of approximately 534 mV compared to ReSe 2 with a D F of 1.40 ± 0.02 (Fig. 2f). Its Tafel slope is 88.9 mV dec -1 , indicating a Volmer- Heyrovsky mechanism with the Heyrovsky reaction as the rate-determining step 48,49 (Fig. 2g and Supplementary Fig. 11). The increase in D F promotes the rate-determining step to transform from the Volmer reaction (H₃O⁺ + e⁻ → H₂O + H*) to the Heyrovsky reaction (H₃O⁺ + e⁻ → H₂O + H₂). The positive correlation between the HER performance of dendritic structures and their D F not only verifies D F as an effective descriptor to link morphological characteristics and catalytic properties but also suggests the efficacy of ML algorithms in optimizing material performance. Process parameter-D F model construction and interpretation After finding an optimal parameter combination for the CVD synthesis of highly branched ReSe 2 , we investigated the construction of an accurate relationship between the five process variables and D F (Module 2). Six representative ML models, which can be categorized into instance-based (k-nearest neighbors (KNN) regressor), network-based (multilayer perceptron (MLP)), kernel-based (GPR and support vector regressor (SVR)), and tree-based (extra trees (ET) regressor and eXtreme Gradient Boosting (XGBoost) regressor) algorithms, are selected for training using the 60 data points accumulated in Module 1. The coefficient of determination (R 2 ) and mean square error (MSE), two metrics widely used to quantify the goodness of fit of a regression line to actual data, are employed to evaluate the model’s performance (see Methods, “Model evaluation”). The best model should display both the maximal R 2 and the minimal MSE. Five repeated 5- fold cross-validations are conducted for each model to mitigate potential biases introduced by random data partitioning. Fig. 3a shows that two tree-based models outperform the other types of algorithms, where XGBoost wins the championship (Supplementary Fig. 12). It achieves an R 2 score of 0.74 and an MSE of 3.4×10 -3 , which are 8.8% higher and 22.7% lower than those of the suboptimal model (i.e., ET). However, due to the inhomogeneous distribution of data collected in Module 1, fitting for the process parameter-D F curve is still not satisfactory across the entire feasible range of the parametric space, even utilizing the best XGBoost regressor. To improve the model predictive capability with minimal experiment addition, we develop a D F prediction accuracy-guided data augmentation strategy, where a metric (called PA score), expressed as the square of the difference between the predicted and the experimental D F values, is introduced to quantify the prediction error of each data point by the current model: PA score = 1 n ∑ (y pred. − y exp. ) 2 n i=1 (2) where y exp. and y pred. represent the experimentally measured and model-predicted D F values, respectively. n represents the number of times the 5-fold cross-validation experiments are conducted. In this work, we performed 5 repeated 5-fold cross-validations for each data point to achieve its average prediction error (i.e., n = 5). To conduct data augmentation, we sorted the PA scores of all data points in Module 1 in descending order, chose the top three (red dots, left panel in Fig. 3b, Supplementary Table 4), and supplemented experiments in the vicinity of these three points in the process parameter space (green dots, 9 middle panel, Supplementary Table 5). This means the CVD process parameters of these new experiments are similar to those of the red dots, but not the same. Results of the supplemented experiments were incorporated into the original dataset for model retraining, which are expected to acquire significant improvement in the original worst-fit regions (right panel). After circulating this data augmentation process 3 times with only 9 experiments added in total (Supplementary Tables 6-9), R 2 and MSE of the XGBoost model are optimized step-by-step to 0.86 and 2.0×10 -3 , which are 16.2% higher and 41.2% lower than before the data supplementation (Figs. 3c-e). In comparison, if 9 data points are added at once rather than divided into three iterative rounds, the model’s MSE only decreases by 14.7% with an 8.1% enhancement of R 2 (Supplementary Table 10 and Supplementary Fig. 13). The global mapping of process parameter-D F is displayed in Fig. 5 with a detailed analysis of the multifactorial growth mechanism. Fig. 3. Model construction and interpretation to reveal the process parameter-D F relationship. a, Performances of the KNN, MLP, GPR, SVR, ET, and XGBoost models, evaluated by R 2 and MSE using the 60 data points collected in Module 1 for training. b, Schematic illustrating the principle of the D F prediction accuracy-guided data augmentation strategy. c, Prediction performance of the XGBoost model before data augmentation. Yellow circles mark the three data points having the largest D F prediction errors, which help locate the regions for data supplementation in the first iterative round. d, Evolution of R 2 and MSE values of the XGBoost model with the increase of the supplementary data points. e, Prediction performance of the XGBoost model after data augmentation with 9 experiments added. f, Feature importance of five CVD process parameters on D F determined by normalized SHAP values. g, Beeswarm plot of SHAP values. h, SHAP interaction value matrix in the logarithmic form based on the XGBoost model after data augmentation. The inset is a histogram displaying the I scores of the five process parameters. While the process parameter-D F mapping regressed by the XGBoost model provides efficient guidance for the synthesis of 2D ReSe 2 with various user-defined D F (Supplementary Fig. 14 and Supplementary 10 Table 11), the black-box nature of this ML model poses a challenge to understanding the internal rationale behind predictions. To improve the model’s interpretability and uncover the process-D F correlations for better mechanism deciphering (Module 3), we employed SHAP analysis, which can explain predictions of any ML model by calculating the importance of the input features and is versatile to provide both local explanations (per data point) and global feature importance. Fig. 3f shows that T Re is the most influential process parameter to D F , whose contribution accounts for nearly 50%. c Re ranks second with a dramatic SHAP value decrease to 26.6%. The remaining three parameters, T Se , sub., and f H2 , account for less than a quarter of the overall contribution. This suggests that T Re and c Re play dominant roles in tailoring the morphology of 2D ReSe 2 dendrites and should be the focus of the domain knowledge-incorporated mechanism elucidation in the next section. It is worth noting that feature importance is not in alignment with the convergence behavior of different features in active learning. The input variable of sub., which displays the fastest convergence rate in BO (Fig. 2c), has the second-lowest feature importance among the five parameters. This demonstrates the necessity of our proposed new metric, I score, which assesses the contribution independence of each feature to D F , to explain the discrepancy in convergence behaviors of different features. To get insights into the detailed impact of each feature on the model’s output, we drew a beeswarm plot of the SHAP values (Fig. 3g), where five input features are arranged in a row and each dot represents a sample in the dataset with its color indicating the value of the feature. The horizontal position of each dot is determined by the SHAP value, indicating the magnitude and direction of a feature’s impact on D F . A positive value refers to a positive impact, and vice versa. Taking T Re as an example, its widest distribution of SHAP values from −0.206 to +0.079 compared to the other inputs implies the strongest influence of this feature on D F . Both high and low T Re values make a negative contribution to D F , while only intermediate ranges enhance D F . Such a non-monotonic impact on the model’s output is also observed in 푐 Re . In addition, the c-Al₂O₃ substrate (represented by a high feature value) demonstrates obvious D F -promoting capability compared to MgO. Finally, the SHAP interaction matrix shows that the main effect of the features (i.e., the interaction effect of features with themselves) rather than the interaction effect between different features plays a dominant role to D F , as manifested by larger values of the diagonal elements than non-diagonal elements (Fig. 3h). It indicates the absence of strong entanglement effects of these CVD process parameters on D F . Besides, The I score ranking of the five features based on the XGBoost model after data augmentation is consistent with that obtained in active learning based on the GPR model, where T Re and sub. are the top two, verifying their strong independence level of impact on D F regardless of the algorithm type and sample size. Multiscale characterizations Apart from extracting the apparent correlations between process parameters and D F using interpretable ML algorithms, we also carried out multiscale microscopic and spectroscopic characterizations to unveil the morphology, orientation, and crystallinity of 2D ReSe 2 dendrites, which facilitates the proposal of a comprehensive and reliable growth mechanism (Module 3). Three key findings were discovered. First, 2D ReSe 2 crystals evolve from isotropic circles to anisotropic dendrites as T Re increases from 580 to 660°C (Fig. 4a). It is noteworthy that the fractal dimension of ReSe 2 grown at T Re = 620°C (D F =1.7) is higher than that at T Re = 660°C (D F =1.5), which is consistent with the results implied by SHAP 11 analysis (Fig. 3g) and further verifies the non-monotonic change of D F with T Re . Second, ReSe 2 dendrites exhibit oriented growth strongly correlated with the substrate symmetry. For ReSe 2 grown on c-Al 2 O 3 with three-fold rotational symmetry (C 3v ), two preferential orientations twisted by 60° are present, which is revealed by the angle concentration at 0° and 60° when measuring the smallest misorientation angle (α 1 ) between c-Al 2 O 3 [112 ̅ 0] and the primary branches of ReSe 2 (Fig. 4b, top). Moreover, the minimum angles between the primary branches (α 2 ) and the minimum angles between the primary and secondary branches (α 3 ) are concentrated at 120° and 60°, respectively, showing orientation preference that corresponds to C 3v symmetry of c-Al 2 O 3 (Fig. 4b, bottom). When the substrate changes to MgO(001), ReSe 2 dendrites are oriented along MgO[110] with the inter-branch angles (β 2 and β 3 ) concentrated at 90° (Fig. 4c, Supplementary Fig. 15), matching well with the four-fold rotational symmetry (C 4v ) of this substrate surface. Third, the CVD-grown 2D ReSe 2 dendrites are polycrystalline. We measured the angle-dependent Raman intensities of the vibrational model at 162 cm -1 in the parallel polarization configuration at different locations on one ReSe 2 crystal, where disordered patterns were achieved (Fig. 4d, blue and magenta lines), showing stark contrast to the well-oriented four-lobed patterns captured from an exfoliated single- crystalline sample (Fig. 4d, gray and brown lines). The ADF-STEM image shows that the edges of ReSe 2 dendrites are rugged, consisting of abundant ribbon-like structures that can be deflected by 60° (yellow markers, Fig. 4e, top), which demonstrates the edge roughness and the polycrystalline nature of ReSe 2 dendrites at the nanometer scale. When zooming in at a ribbon structure, well-organized Re4 chains are observed with negligible defects at the atomic scale, indicating a high level of local crystal quality (Fig. 4e, bottom). Fig. 4. Cross-scale characterizations of 2D ReSe 2 dendrites. a, SEM images showing the shape evolution of 2D ReSe 2 grown at T Re =580, 620, and 660°C, respectively. The other CVD process parameters are fixed at c Re = 0.10 mol/L, sub.= c-Al 2 O 3 , T Se = 250°C, f H2 = 2 sccm. Scale bars (from left to right): 10 μm, 5 μm, 5 μm. b, Histogram showing the angle distribution of α 1 , α 2 , and α 3 of 2D ReSe 2 grown on c-Al 2 O 3 . α 1 : the minimum angle between c- Al 2 O 3 [112 ̅ 0] and the primary branch of ReSe 2. α 2 : the minimum angle between the primary branches of ReSe 2 . α 3 : the minimum angle between the primary and secondary branches of ReSe 2 . Insets schematically illustrate α 1 , α 2 , and 12 α 3 . c, Histogram showing the angle distribution of β 1 , β 2 , and β 3 of 2D ReSe 2 grown on MgO(001). β 1 : the minimum angle between MgO[110] and the primary branch of ReSe 2. β 2 : the minimum angle between the primary branches of ReSe 2 . β 3 : the minimum angle between the primary and secondary branches of ReSe 2 . Insets schematically illustrate β 1 , β 2 , and β 3 . d, Angle-resolved Raman intensities of the vibrational model of ReSe 2 at 162 cm -1 captured at two different locations on a CVD-grown ReSe 2 dendritic crystal (blue and magenta) and two different locations on a single-crystalline exfoliated ReSe 2 (gray and brown). e, Medium- and high-magnification ADF-STEM images of ReSe 2 dendrites. The yellow arrows indicate the b-axis of ReSe 2 (top), while the rhombus markers represent the Re4 chains (bottom). Scale bars: 100 nm (top), 1 nm (bottom). Multifactorial mechanism deciphering We integrated quantitative information given by machine intelligence with qualitative understanding from material characterizations to establish a data–knowledge dual-driven growth model of 2D ReSe 2 . Figure 5 (middle panel) displays a non-monotonic, non-linear correlation between the CVD process parameters and the ReSe 2 crystal morphology on the c-Al 2 O 3 substrate, which is fitted by the XGBoost model using 69 sets of experimental data accumulated in Fig. 3. It involves the top two important features (T Re and c Re ), identified by SHAP analysis, as the two x-coordinate axes, and D F as the y-axis. Combining results from ML and characterizations, key phenomena of 2D ReSe 2 grown on c-Al 2 O 3 can be summarized: (i) When T Re is below 600°C, 2D ReSe 2 domains exhibit isotropic, circular, non-dendritic, polycrystalline structures; (i) When T Re rises above 600°C, 2D ReSe 2 polycrystals transform to anisotropic, dendritic morphology with branch orientations strongly correlated with the substrate symmetry; (i) At T Re > 600°C, different combinations of T Re and c Re lead to a wide-range modulation of D F from 1.4 to 1.7 (labels (i) to (iv) in Fig. 5, middle panel); (iv) T Re has the strongest impact (~50%) on D F followed by c Re (~27%). By incorporating domain knowledge, we proposed a multifactorial growth mechanism to explain these phenomena. When T Re is below 600°C, 2D ReSe 2 is thermodynamically controlled by an attachment- limited growth, meaning that the adsorption and bonding of the precursor species to the domain edges is the rate-dominant step (Fig. 5, left panel). In this case, the surface energy of the system determines the preferential adsorption sites of the active species. Compared with crystal tips and edges, the attachment of precursor molecules to concavities reduces the surface energy the most. Therefore, any protrusions at domain edges will be quickly smoothed out under this mechanism, prohibiting the formation of dendritic structures thermodynamically. In addition, previous literature has reported that Re4 chains are prone to deflect integer multiples of ～60° on c-Al 2 O 3 during crystal growth due to the capability of forming highly coherent grain boundaries with low energy 50,51 . This not only generates polycrystalline structures of ReSe 2 but also leads to a circular crystal shape due to isotropic growth rates in all lateral directions independent of either the substrate surface symmetry or the intrinsic crystal structure of the triclinic ReSe 2 . When T Re exceeds 600°C, ReSe 2 swifts to kinetically controlled diffusion-limited growth (Fig. 5, right panel), where migration of precursor species on the substrate is the rate-dominant step. Our previous work has applied density functional theory calculations to unveil that the diffusion energy barrier of Re atoms on c-Al 2 O 3 displays prominent anisotropy dependent on the lattice orientation of c-Al 2 O 3 52 . Therefore, the macroscopic shape of polycrystalline dendrites and the preferential growth orientations of branches exhibit obvious association with the C 3v symmetry of c-Al 2 O 3 . It is noteworthy that, under diffusion-limited synthesis, the growth rate of 2D ReSe 2 at a certain temperature is affected by the concentration of precursor species absorbed on the substrate surface (commonly refers to Re concentration, because Se is in excess). 13 In our space-confined CVD setup constructed by two substrates placed face-to-face, Re precursor is spin- coated onto the lower substrate, which can be vaporized to adsorb onto the upper c-Al 2 O 3 surface via heating for ReSe 2 dendrite growth (Supplementary Fig. 16). Under this setup, the Re concentration adsorbed on the upper substrate surface is jointly influenced by T Re and c Re , because T Re dictates both the heating temperature of the bottom substrate coated with Re source and the temperature of the upper substrate at which adsorption equilibrium occurs, while c Re impacts the supply amount of the Re source from the bottom substrate. Therefore, different combinations of T Re and c Re can vary D F of dendrites. Fig. 5. Multifactorial growth mechanism deciphering of 2D ReSe 2 dendrites. Middle panel: Correlation between the most important CVD parameters (T Re and c Re ) and D F of 2D ReSe 2 crystals on c-Al 2 O 3 fitted by the XGBoost model. Left panel: Attachment-limited growth when T Re is below 600°C, where the formation of dendrites is thermodynamically prohibited due to a higher growth rate at the concavity region of the branch bottom (v c ) than at the branch tip (v t ). Right panel: Diffusion-limited growth when T Re rises above 600°C, where a negative concentration gradient at the crystal edges results in a higher growth rate at the branch tip than at the branch bottom (v t > v c ), allowing the formation of dendrites. v t-i, v t-i , v t-i , and v t-iv represent the growth rates of branch tips in cases (i) to (iv), respectively, which vary due to different combinations of T Re and c Re . Take cases (i) and (i) corresponding to low T Re with high c Re and low T Re with low c Re as an example to compare their products’ D F (Fig. 5, bottom right). Due to their same T Re , the equilibrium Re adsorption concentration at the upper substrate surface (defined as c equ. ) is solely determined by c Re , thereby resulting in a higher c equ. of case (i) than that of case (i). Furthermore, the heating temperature of the upper substrate, which equals T Re , affects the migration ability of Re atoms on the substrate surface, thus determining the thickness of the depletion zone at the ReSe 2 crystal edge. The lower the T Re temperature is, the narrower the depletion zone thickness becomes 53,54 . Here, the depletion zone is the region with a negative concentration gradient of precursor species, extending from the edge of the ReSe₂ crystal to the location where the Re adsorption concentration reaches c equ. . When a branch is accidentally generated at the edge of a ReSe 2 domain due to energy perturbation and inserted into the depletion layer, the Re concentration at the tip of the branch is higher than at its base due to the presence of a negative concentration gradient within the depletion layer, which makes the branch not to be smoothed out and disappear, but rather to grow longer and longer. This explains why dendrites can be produced in diffusion-limited growth. Since case (i) has a higher c equ. than case (i) and both cases adopt a low T Re , a sharper negative Re concentration gradient is produced in the depletion layer under case (i) than under case (i). Therefore, the branch tip under case (i) 14 is situated at a higher Re concentration than under case (i), which gives ReSe 2 dendrites in case (i) a larger D F . The D F variations by other combinations of T Re and c Re can be demonstrated in a similar way. Discussion In this work, we display full-chain support of ML approaches to material growth under small datasets, including efficient recipe optimization, customized synthesis, and multifactorial mechanism deciphering, exemplified by the CVD growth of 2D ReSe 2 dendrites. First, the BO-GPR algorithm is performed, which improves D F medians by 69.4% and achieves HER-active dendrites with a high D F of 1.71 using 4 iterations (60 experiments), occupying less than 1.3% of the 4752 possible combinations in a grid search. Second, a prediction accuracy-guided data augmentation strategy is developed, which utilizes only 9 additional experiments combined with an XGBoost model to realize the establishment of a nonlinear mapping between 5 process variables and dendritic morphology of 2D ReSe 2 . The model’s predictive capability reaches R 2 =0.86 across the entire feasible range and formulates data-informed guidelines to synthesize ReSe 2 dendrites with various user-defined D F . Finally, we jointly apply cross-scale microscopic/spectroscopic characterizations, expert knowledge in thermodynamics and kinetics, and interpretable ML methods (SHAP feature analysis) to construct a data–knowledge dual-driven growth mechanism of 2D ReSe 2 dendrites, which unravels independent and synergistic contributions of multiple CVD parameters (T Re , c Re and sub.) to D F . In addition, this work proposes two measures from the perspective of machine intelligence, namely I score and PA score. The former evaluates the independence of different features in contributing to the model output, so as to better understand the convergence behavior of input features during active learning process. The latter can guide the sampling strategy to improve the model’s predictive capability with minimum experimental budget. Our approach can be generalized to the preparation of diversified materials and, more importantly, establishes a stylized ML-empowered framework for rapid material development and precise recipe-property relation establishment, which may revolutionize the research logic of material synthetic methodology. Methods Preparation and transfer of 2D ReSe 2 (i) Synthesis of 2D ReSe 2 . The growth of 2D ReSe 2 was performed with a spatially confined CVD system 16,18 . The spin-coated substrate served as the bottom chip, with another c-Al 2 O 3 or MgO substrate placed face-to-face on top of it to form a space-confined reaction microchamber, which was positioned at the center of the high-temperature zone of the furnace (580-680°C). A 0.025-0.15 mol/L solution of ammonium perrhenate (NH 4 ReO 4 , 99%, RHAWN) was applied as a rhenium source and spin-coated onto the plasma-treated substrate surface at 5000 rpm for 30 seconds. The temperature of this substrate is the rhenium source temperature (T Re ). c-Al 2 O 3 or MgO was treated with oxygen plasma to create a hydrophilic surface. For selenium (Se) source, 80 mg of Se powder (99.99%, RHAWN) was put in a ceramic boat at the center of the temperature zone of the furnace (220-300°C), and this temperature is defined as the T se . Before 2D ReSe 2 synthesis, the CVD system was flushed with argon (Ar) gas flow for 10 min to remove the air from the tube. Subsequently, a mixed gas of Ar/H 2 was introduced, with Ar and H 2 gas flow rates (f H2 ) of 100 sccm and 1~4 sccm, respectively. The CVD system was maintained under the above-mentioned 15 conditions for 20 minutes to ensure the growth of ReSe 2 on the target substrate. After the growth period, the furnace was naturally cooled to room temperature under a continuous flow of Ar at ~500 sccm. (i) Transfer of 2D ReSe 2 . A thin film of poly methyl methacrylate (PMMA) was spin-coated onto the ReSe 2 /substrate surface. The coated substrates were baked at 80 °C for 3 minutes to strengthen the interaction between the PMMA film and the ReSe 2 and then floated on a 2 mol/L potassium hydroxide (KOH) solution to separate the PMMA/ReSe 2 film from the substrate. The detached PMMA/ReSe 2 film was transferred to deionized water three times to remove the residual etchants. Then, the film was transferred to a TEM grid or a glass carbon electrode, dried naturally in air, and baked on a hotplate at 100°C for 15 min. The sample was immersed in acetone for 8 h to remove the PMMA. Characterizations (i) Raman spectroscopy. Raman spectroscopy was employed to study structural properties, utilizing a Renishaw inVia confocal Raman system with an excitation wavelength of 532 nm. The laser power was kept below 100 μW, and the system was calibrated with the Raman peak of silicon at 520.7 cm -1 . The experimentally measured data of A 1g -like mode at 162 cm⁻¹ with parallel polarization x̅ (y )x or configurations are plotted as the blue and red dots, respectively. Under parallel polarization configuration, the angular dependence of the Raman intensity for this mode follows equations 55,56 ： I A g ∥ ∝ | a | 2 cos 4 θ+ | d | 2 sin 4 θ+4 | b | 2 sin 2 θcos 2 θ +2 | a || d | cosφ ad sin 2 θcos 2 θ(3) +4 | a || b | cosφ ab sinθcos 3 θ +4 | d || b | cosφ db sin 3 θcosθ Here, I A g ∥ is the Raman intensity for the A g -like (162 cm -1 ) mode in parallel polarization configurations. | a | , | b | and | d | are the value of Raman tensor element a, b, and d. θ is the angle between the incident laser polarization and the Re-chain direction (b-axis). φ a d =φ a -φ d , φ a b =φ a -φ b , and φ db =φ d -φ b . (i) ADF-STEM imaging. The measurements were conducted at room temperature on aberration-corrected Titan Cubed Themis G2 300. STEM was operated under an accelerating voltage of 300 kV. The experimental setup included a condenser lens with an aperture diameter of 50 m. The convergence semi- angle was set to 21.3 mrad, and the collection angle ranged from 39 to 200 mrad. For each frame, the dwell time was set to 2 μs per pixel. Imaging was performed at a pixel size of 0.012 nm and a beam current of 30 pA. Electrochemical measurements HER measurements for 2D ReSe 2 with different D F were performed at room temperature using a CHI 760E electrochemical workstation. The experimental setup consisted of a glassy carbon disk as the working electrode, a platinum plate as the counter electrode, and a saturated Ag/AgCl electrode as the reference, forming a standard three-electrode configuration. Linear sweep voltammetry (LSV) was conducted at a scan rate of 5 mV s −1 in a 0.5 M H 2 SO 4 electrolyte. Tafel slopes were obtained from the linear portions of the LSV curves using the Tafel equation, with current densities normalized based on the material’s exposed surface area. Potentials were referenced to the reversible hydrogen electrode (RHE) using the relation E (vs RHE) = E (vs Ag/AgCl) + 0.197 V + 0.0592 pH. 16 ML models (i) Bayesian optimization. BO operates by constructing a surrogate model to approximate the objective function and employs an acquisition function to guide subsequent experimental sampling, thereby identifying the global optimum with minimal experimental iterations. We set up the BO framework through the package Emukit in Python 36,57 . The GPR model was selected as the surrogate model 43 . GPR infers the distribution of the objective function across the parameter space using a limited number of data points and provides uncertainty estimates for predictions. Its mathematical formulation is as follows: F(x)~N (μ(x), K(x, x’)) (4) where μ(x) is the mean function and K(x, x’) is the covariance function (kernel function). The acquisition function recommends the experimental sampling points of the next iteration based on the surrogate model’s predictions, balancing exploration (sampling regions with high uncertainty) and exploitation (regions predicted to yield near-optimal values). We employ the MES acquisition function, which selects sampling points by maximizing the information gain associated with the optimal value of the objective function 44 . Its mathematical formulation is as follows: MES(x)=H(p(f(x)))-E p(f(x * )) [H (p(f(x)|f(x * )))] (5) while f(x) and f(x * )represent the predicted value and the predicted maximum value of the GPR, respectively. It is worth noting that all data were normalized prior to being fed into the model. (i) Model evaluation. Prior to feature analysis, ML models were evaluated using five repetitions of 5- repeated 5-fold cross-validation. Two indicators are selected to evaluate the model's performance, which are R 2 and MSE. R 2 is a statistical measure that provides information about the goodness of fit of a model. In the context of regression, it is a statistical measure of how well the regression line approximates the actual data. R 2 is (−∞,1], where a value closer to 1 indicates better model performance. MSE measures the average squared difference between predicted and actual values, which disproportionately weights larger deviations, resulting in greater sensitivity to outlier values. R 2 and MSE were utilized to assess model performance, with their calculations shown as follows: R 2 = 1- ∑ (y pred. –y exp. ) 2 n i=1 / ∑ (y pred. – y mean ) 2 n i=1 (6) MSE = 1 n ∑ (y pred. -y exp. ) 2 n i=1 (7) where y pred. , y exp. and y mean are the predicted D F value, the experimental D F value and the average of all the measured values, respectively. (i) XGBoost model. As an advanced implementation of gradient-boosted trees, XGBoost distinguishes itself through its efficient regularization techniques and parallel computing capabilities 58 . Its core principle involves iteratively constructing multiple decision tree models and integrating them into a robust ensemble framework to enhance predictive accuracy. Its mathematical formulation can be expressed as: y pred. = ∑ f m (x i ) M m=1 (8) where f m denotes the base estimator, y pred. donates the predicted D F value, and x i donates the growth parameters value. m=1,2,..., M. 17 The objective function comprises a loss function and a regularization term, mathematically expressed as: Obj(θ) = ∑ L N i=1 (y pred. , y exp. )+ ∑ Ω(f m ) 푀 푚=1 (9) where L(y pred. , y exp. ) represents the MSE loss function, and ∑ Ω(f m ) 푀 푚=1 is the regularization term designed to control model complexity and mitigate overfitting. i = 1, 2,..., N. Model interpretability (i) Quantification of feature importance. We employed SHAP method to quantitatively assess the contribution of each input feature to predictions in the XGBoost model, thereby enhancing ML model’s interpretability 46 . Feature importance values are assigned by SHAP through calculating the marginal contributions of all possible feature subsets for each data instance. Global feature importance is derived by aggregating the mean absolute SHAP values across the entire dataset. SHAP also enables local interpretability by explaining predictions for individual instances through SHAP values. Features with larger mean absolute SHAP values exhibit higher importance and vice versa. (i) Quantification of the independence level of feature’s impact on D F . SHAP values can be decomposed into the independent effects of features (main effects, SHAP main effect x i ) and the mutual influences between features (interaction effects, SHAP interaction x i -x j , i≠j). For any model prediction f(x), SHAP guarantees the following additive decomposition: f(x) = φ 0 +SHAP main effect x i +SHAP interaction x i -x j (10) where φ 0 is the baseline value, i.e., the model prediction when no features are considered. Therefore, using Equation (1), which is the ratio of the main effect to the interaction effect, to quantify the independence level of a feature’s impact on D F is reasonable and easy to understand. Range of D F . In this paper, we focus on the linear fractal dimension of dendrites, where D F ∈ (1,2). Specifically, a value of 1 corresponds to a closed, smooth curve (exhibiting minimal geometric complexity). A value approaching 2 indicates a structure with non-zero area, which strictly speaking no longer qualifies as a “line” but instead approximates a plane. The D F medians rose from 1.36 (iteration 0, initial dataset) to 1.61 (iteration 4). The improvement amount can be calculated using the normalization principle expressed as: (1.61-1.36)/(1.36-1)×100%= 69.4%. Data availability Relevant data supporting the key findings of this study are available within the article and the Supplementary Information file. All raw data generated during the current study are available from the corresponding authors upon request. Code availability All data supporting this study and the code are available on GitHub (https://github.com/csuhwq0421/ML4Dendrites). References 1. Zhang, M. et al. Ultrahigh energy storage in high-entropy ceramic capacitors with polymorphic relaxor phase. Science. 384, 185–189 (2024). 18 2. Li, L. et al. Self-heating–induced healing of lithium dendrites. Science. 359, 1513–1516 (2018). 3. Ma, M. C., Li, G., Chen, X., Archer, L. A. & Wan, J. Suppression of dendrite growth by cross-flow in microfluidics. Sci. Adv. 7, eabf694119 (2021). 4. Ning, Z. et al. Dendrite initiation and propagation in lithium metal solid-state batteries. Nature 618, 287–293 (2023). 5. Xiao, B. et al. Preparation of hierarchical WO 3 dendrites and their applications in NO 2 sensing. Ceram. Int. 43, 8183–8189 (2017). 6. Ji, Q. et al. Morphological engineering of CVD-grown transition metal dichalcogenides for efficient electrochemical hydrogen evolution. Adv. Mater. 28, 6207–6212 (2016). 7. Shao, G. et al. Shape-engineered synthesis of atomically thin 1T-SnS 2 catalyzed by potassium halides. ACS Nano 13, 8265–8274 (2019). 8. Chen, J. et al. Homoepitaxial growth of large-scale highly organized transition metal dichalcogenide patterns. Adv. Mater. 30, 1704674 (2018). 9. Wang, B. et al. Fractal growth of 2D NbSe 2 for broadband nonlinear optical limiting. Adv. Funct. Mater. 34, 2401490 (2024). 10. Kang, Y. et al. Conductive dendrite engineering of single-crystalline two-dimensional dielectric memristors. Innovation 6, 100885 (2025). 11. Wang, M. et al. Single-crystal, large-area, fold-free monolayer graphene. Nature 596, 519–524 (2021). 12. Liu, L. et al. Uniform nucleation and epitaxy of bilayer molybdenum disulfide on sapphire. Nature 605, 69–75 (2022). 13. Wang, L. et al. Bevel-edge epitaxy of ferroelectric rhombohedral boron nitride single crystal. Nature 629, 74–79 (2024). 14. Moon, D. et al. Hypotaxy of wafer-scale single-crystal transition metal dichalcogenides. Nature 638, 957–964 (2025). 15. Wang, J. et al. Dual-coupling-guided epitaxial growth of wafer-scale single-crystal WS 2 monolayer on vicinal a-plane sapphire. Nat. Nanotechnol. 17, 33–38 (2022). 16. Jiang, J. et al. Chirality-transferred epitaxy of circular polarization-sensitive ReS 2 monolayer single crystals. Nat. Commun. 16, 7119 (2025). 17. Wang, S. et al. Atomic-scale studies of overlapping grain boundaries between parallel and quasi- parallel grains in low-symmetry monolayer ReS 2 . Matter 3, 2108–2123 (2020). 18. Chen, Y. et al. Constructing slip stacking diversity in van der waals homobilayers. Adv. Mater. 36, 2404734 (2024). 19. Sahoo, P. K., Memaran, S., Xin, Y., Balicas, L. & Gutiérrez, H. R. One-pot growth of two- dimensional lateral heterostructures via sequential edge-epitaxy. Nature 553, 63–67 (2018). 20. Li, W. et al. Robust growth of 2D transition metal dichalcogenide vertical heterostructures via ammonium-assisted CVD strategy. Adv. Mater. 36, 2408367 (2024). 21. Xie, S. et al. Coherent, atomically thin transition-metal dichalcogenide superlattices with engineered strain. Science. 359, 1131–1136 (2018). 22. Zhao, Y. et al. Supertwisted spirals of layered materials enabled by growth on non-Euclidean surfaces. Science. 370, 442–445 (2020). 23. Xu, M. et al. Reconfiguring nucleation for CVD growth of twisted bilayer MoS 2 with a wide range of twist angles. Nat. Commun. 15, 562 (2024). 24. Fortin-Deschênes, M., Watanabe, K., Taniguchi, T. & Xia, F. Van der Waals epitaxy of tunable moirés enabled by alloying. Nat. Mater. 23, 339–346 (2024). 25. Zhang, Y. et al. Monolayer MoS 2 dendrites on a symmetry-disparate SrTiO 3 (001) substrate: formation mechanism and interface interaction. Adv. Funct. Mater. 26, 3299–3305 (2016). 19 26. Li, J. et al. Fractal-theory-based control of the shape and quality of CVD-grown 2D materials. Adv. Mater. 31, 1902431 (2019). 27. Dan, J. et al. Learning motifs and their hierarchies in atomic resolution microscopy. Sci. Adv. 8, eabk1005 (2022). 28. Zhao, Q. et al. A machine learning–based framework for mapping hydrogen at the atomic scale. Proc. Natl. Acad. Sci. 121, e2410968121 (2024). 29. Huang, W. et al. Auto-resolving the atomic structure at van der Waals interfaces using a generative model. Nat. Commun. 16, 2927 (2025). 30. Luo, Z. et al. Exploring structure diversity in atomic resolution microscopy with graph. Adv. Mater. 37, 2417478 (2025). 31. Crozier, P. A. et al. Visualizing nanoparticle surface dynamics and instabilities enabled by deep denoising. Science. 387, 949–954 (2025). 32. Chen, T. et al. Machine intelligence-accelerated discovery of all-natural plastic substitutes. Nat. Nanotechnol. 19, 782–791 (2024). 33. Moon, J. et al. Active learning guides discovery of a champion four-metal perovskite oxide for oxygen evolution electrocatalysis. Nat. Mater. 23, 108–115 (2024). 34. Suvarna, M. et al. Active learning streamlines development of high performance catalysts for higher alcohol synthesis. Nat. Commun. 15, 5844 (2024). 35. Kim, M. et al. Searching for an optimal multi-metallic alloy catalyst by active learning combined with experiments. Adv. Mater. 34, 2108900 (2022). 36. Liu, Z. et al. Machine learning with knowledge constraints for process optimization of open-air perovskite solar cell manufacturing. Joule 6, 834–849 (2022). 37. Lin, X. et al. Machine learning-assisted dual-atom sites design with interpretable descriptors unifying electrocatalytic reactions. Nat. Commun. 15, 8169 (2024). 38. Tang, B. et al. Machine learning-guided synthesis of advanced inorganic materials. Mater. Today 41, 72–80 (2020). 39. Zhang, J. et al. Toward controlled synthesis of 2D crystals by CVD: learning from the real-time crystal morphology evolutions. Nano Lett. 24, 2465–2472 (2024). 40. Guo, H. et al. Machine learning-guided realization of full-color high-quantum-yield carbon quantum dots. Nat. Commun. 15, 4843 (2024). 41. Li, Y. et al. Transforming the synthesis of carbon nanotubes with machine learning models and automation. Matter 8, 101913 (2025). 42. Shields, B. J. et al. Bayesian reaction optimization as a tool for chemical synthesis. Nature 590, 89– 96 (2021). 43. Seeger, M. Gaussian processes for machine learning. Int. J. Neural Syst. 14, 69–106 (2004). 44. Wang, Z. & Jegelka, S. Max-value entropy search for efficient Bayesian optimization. In Proc. of the 34th International Conference on Machine Learning, 7, 5530–5543 (2017). 45. Lundberg, S. M. & Lee, S. I. A unified approach to interpreting model predictions. In Proc. of the 31st International Conference on Neural Information Processing Systems, 4766–4775 (2017). 46. Lundberg, S. M. et al. From local explanations to global understanding with explainable AI for trees. Nat. Mach. Intell. 2, 56–67 (2020). 47. Muschalik, M. et al. shapiq: shapley interactions for machine learning. In Proc. of the 38th International Conference on Neural Information Processing Systems. 4141, 130324-130357 (2024). 48. Xie, J. et al. Controllable disorder engineering in oxygen-incorporated MoS 2 ultrathin nanosheets for efficient hydrogen evolution. J. Am. Chem. Soc. 135, 17881–17888 (2013). 49. Wang, X. et al. Single-atom vacancy defect to trigger high-efficiency hydrogen evolution of MoS 2 . J. Am. Chem. Soc. 142, 4298–4308 (2020). 20 50. Li, X. et al. Nanoassembly growth model for subdomain and grain boundary formation in 1T′ layered ReS 2 . Adv. Funct. Mater. 29, 1906385 (2019). 51. Chen, X. et al. Diverse spin-polarized in-gap states at grain boundaries of rhenium dichalcogenides induced by unsaturated Re-Re bonding. ACS Mater. Lett. 3, 1513–1520 (2021). 52. Hu, P. et al. Lateral and vertical morphology engineering of low-symmetry, weakly-coupled 2D ReS 2 . Adv. Funct. Mater. 33, 2210502 (2023). 53. Xu, J., Srolovitz, D. J. & Ho, D. The adatom doncentration profile: a paradigm for understanding two-dimensional MoS 2 morphological evolution in chemical vapor dposition growth. ACS Nano 15, 6839–6848 (2021). 54. Zhuang, J. et al. Morphology evolution of graphene during chemical vapor deposition growth: a phase-field theory simulation. J. Phys. Chem. C 123, 9902–9908 (2019). 55. Choi, Y. et al. Complete determination of the crystallographic orientation of ReX 2 (X = S, Se) by polarized Raman spectroscopy. Nanoscale Horizons 5, 308–315 (2020). 56. Zhao, H. et al. Interlayer interactions in anisotropic atomically thin rhenium diselenide. Nano Res. 8, 3651–3661 (2015). 57. Paleyes, A., Mahsereci, M. & Lawrence, W. Emukit: A Python toolkit for decision making under uncertainty. In Proc. of the 22nd Python in Science Conference. 68-75 (2023) 58. Chen, T. & Guestrin, C. XGBoost: A scalable tree boosting system. In Proc. of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. 785–794 (2016). Acknowledgments S.W. acknowledges support from the National Natural Science Foundation of China (52222201, 52172032, 22494464013), Major Fundamental Research Project of Hunan Province (2025JC0005), Shenzhen Science and Technology Innovation Commission Project (KQTD20221101115627004), and National University of Defense Technology (ZZCX-ZZGC-01-07). F.O. acknowledges support from the Key Project of the Natural Science Program of Xinjiang Uygur Autonomous Region(Grant No. 2023D01D03) and the computing resources at the High Performance Computing Center of Central South University. Author contributions S.W. initiated the project and generated experimental protocols. W.H. wrote the code. S.W. prepared the STEM samples and captured images. X.G. grew 2D crystals. S.X. and S.F. prepared the samples and conducted measurements. H.X., J.J., S.Z., Z.L., J.Z. and Y. O. discussed the work and gave suggestions. All authors contributed to the data analysis, manuscript writing, and revision of the manuscript. Conflict of interest The authors declare no competing interests.