2.1.1.

Data-driven surrogate sub-modeling

Contributions to this topic have prevalently considered separable submodels as part of a superordinate model consisting of several components. Plasma chemistry presents a specific case detailed later in Sec. 2.2. PSI is similarly addressed separately in Sec. 2.3. A major aspect to be considered in the following concerns the description of electromagnetic fields, in particular in electrostatics description. In this case, the Poisson equation presents an elliptic partial differential equation (PDE) subject to time-dependent boundary conditions (BCs). This problem stands out because of its computational effort and difficulties in an efficient parallel solution in multiple spatial dimensions (as the information of the whole domain is required at every iteration). Early work by Zhang et al.²³ on two-dimensional (2D) Monte Carlo (MC) simulations of a PN junction (conceptually related to methods applied in LTP) has considered reference data from the solution of the inhomogeneous Poisson equation subject to Dirichlet BCs using a finite-difference (FD) method. A multi-scale ANN was trained using the internal and boundary values of the electric potential as ground truth. Good agreement paired with a more than 10 times speedup has been reported. In an attempt to transfer a similar concept to LTPs, Trieschmann et al.²⁴ implemented a one-dimensional (1D) self-supervised learning scheme using an FD discretization to establish a physics-based loss, as, for instance, proposed by Dissanayake and Phan-Thien.²⁵ The data-driven solver was further assessed within a particle-in-cell (PIC) simulation of two collisionless counterpropagating electron streams. A relative error of the order of ${10}^{-4}$ was reported for the electric potential. However, as the force acting on each particle requires the electric field, evaluation of the first derivative implied significant distortion of the long-term spatiotemporal dynamics. This finding has been corroborated by Aguilar et al.²⁶ who have reported similar findings using fully connected ANNs [a.k.a. multilayer perceptrons (MLPs)²⁷] and convolutional neural networks (CNNs).²⁸ A 2D solution to Poisson’s equation coupled to a particle balance equation for a double headed streamer simulation (in drift-diffusion approximation for electrons and otherwise no transport) has been pursued by Cheng et al.²⁹ (closely related to the work Özbay et al.³⁰) Therein, a comparison of different ANN architectures was presented, specifically a multi-scale network and a U-Net configuration.³¹ Whereas a significant speedup has been consistently reported, the accuracy of the predictions was observed to agree only for short term prediction ($<2\mathrm{ns}$). Note that a simple rectangular domain was considered with the addition of a standardization scheme to extend to variable domain sizes. Further, the loss was composed of the inner domain and Dirichlet boundary values, paired with a physics-based loss representing the Laplacian. Zeng et al.³² followed a different paradigm using physics-informed neural networks (PINNs).^33,34 Therein, physical coordinates are input to an ANN and the underlying PDE is implemented as a physics-informed loss function. In general, this loss function may include a set of possibly coupled PDEs to be solved within the domain, the BCs, and the initial conditions in case of transient phenomena.^33,34 Automatic differentiation is exploited to obtain the respective local error gradients with respect to the input variables.³⁵ It resembles a paradigm closely linked to the backpropagation algorithm²⁷ for training ANNs and enables a mesh-less representation of the PDE solution. In the referenced work, the inhomogeneous Poisson equation with rectangular and circular dielectric objects was solved in two and three dimensions, subject to time-dependent Dirichlet and/or Neumann BCs. Note that the authors have initially included time as a variable to consider time-dependent BCs (but later reverted to a time-invariant formulation), although Poisson’s equation essentially poses a static PDE problem. Hence, the error metric was composed of a domain loss representing the Laplacian and the loss corresponding to Dirichlet BCs. Whereas the application to dielectric barrier discharges (DBDs) was proposed as use case, an integration into a complete LTP model remains due. An approach dedicated to the requirements of complex LTP simulations has been put forward by Siffa et al.³⁶ based on a combination of U-Net,³¹ DenseNet,³⁷ and transformer^38,39 network architectures, called TransDenseUNet. On a 2D Cartesian mesh, the non-overlapping domain information was encoded into types (plasma, dielectric, BCs), corresponding physical values (dielectric constant, charge density, and boundary potential/normal derivative), and spatial coordinates. A feature scaling of the physical equation was proposed and all input/output quantities were normalized in prior training. The data set was composed of a few physical reference cases for intended applications and a substantial amount (10 k to 50 k) of randomly generated data samples. The training loss was combined with various contributions [including a gradient loss and a structural similarity index measure (SSIM)⁴⁰] to increase the accuracy and the smoothness of the electric potential prediction (and the associated electric field). The proposed data-driven model was demonstrated to successfully generalize on arbitrary geometries, BCs, and charge distributions—including realistic cases encountered in LTP processing. Also here, an implementation within a complete plasma model remains due.

With the goal of studying the complex interaction of electromagnetic waves with plasma, Desai et al.⁴¹ suggested a U-Net architecture to predict the 2D distribution of the root mean square electric field based on plasma density distributions. The ground truth data set was obtained from finite-difference time-domain simulations. Interestingly, while a simple regularized L2 loss was used, the quality of the prediction was assessed using the SSIM and the average percentage error. Whereas overall good agreement with reference simulations was achieved, an issue generally noticeable with many ANN surrogate models was observed. Specifically, the convergence toward a reliable prediction is subject to a spectral bias. Lower frequency features converge significantly faster than higher frequency features. Consequently, it may be computationally intractable to train until convergence is achieved for all relevant frequencies.⁴²

2.1.2.

Data-driven discharge surrogate modeling

An early example by Verma et al.⁴³ has considered data sets from 1D time-resolved fluid simulations of capacitively coupled radio-frequency discharges to establish surrogate models for integrated quantities, such as the plasma sheath voltages/currents, that capture reduced temporal modes of operation. In contrast to physical model based solutions, the authors have highlighted the equation-free representation of their long short-term memory [LSTM, cf. Sec. 3.1] surrogate model. The corresponding reduced order description and its latent space representation was suggested for parameter space exploration. This may otherwise be infeasible due to computational requirements. Another approach to capacitively coupled plasma surrogate modeling by Hamaguchi and co-workers has been devised based on 1D fluid and particle-in-cell with Monte Carlo collisions (PIC/MCC) simulations of an Ar discharge.⁴⁴ Instead of a grid-based approach for training and prediction (where the whole spatial domain is regressed on at once), MLPs were trained on phase-average plasma quantities sampled at varying positions within the discharge. The dynamics of an inductively coupled plasma with additional radio frequency substrate bias has been implemented based on hybrid plasma equipment model⁴⁵ simulations by Ko et al.⁴⁶ 2D axisymmetric simulations of an Ar etch plasma were performed with a fixed geometry (rectangular mesh) and varying process conditions (pressure, source power, bias power at 1 and 13.56 MHz). The devised data set was further used to train a multi-encoder/decoder ANN, which encodes the process variables and geometry into a latent representation to finally decode the corresponding phase-averaged discharge quantities. Note that although the geometry was fixed in the presented study, the proposed ANN architecture with dedicated geometry inputs (sources, material, and properties) may also be capable of taking varying geometrical features (e.g., chamber aspect ratio, wafer diameter) into account. Due to its low computational costs, the ML surrogate model allowed for a rapid exploration and optimization of a uniform and sufficient etch rate using the multi-objective particle swarm optimization algorithm.

Data-driven surrogate modeling applied to different atmospheric discharges has been proposed by Zhang et al.^47–49 A drift-diffusion model taking into account particle balance coupled to Poisson’s equation, and the electron energy conservation equation was used to simulate pulsed He and ${\mathrm{CO}}_2$ discharges. Corresponding data sets were established for varying source properties and subsequently used for setting up an MLP. A reduced set of input descriptors was selected, such as time and source properties (e.g., voltage, pulse rise rate). As output descriptors the discharge current density or the spatial distributions of species were chosen. With reference to these simulations (and experiments), it was demonstrated that a rather simplistic ANN even in the absence of complex learning schemes is capable of capturing the spatiotemporal discharge dynamics as well as the chemical reaction kinetics.

2.1.3.

Physics-informed machine learning

PINNs were further used to directly consider the equations governing the physics of LTP discharges. Similar to the previously mentioned PINN solution to Poisson’s equation, corresponding PDEs may be implemented in the error metric during training.^33,34 Automatic differentiation is again exploited to obtain desired partial derivatives with respect to the input variables.³⁵ Early work in LTP has been proposed by Kawaguchi et al.⁵⁰ who have considered a PINN solution to the stationary Boltzmann equation subject to ionization and electron attachment processes. An MLP has been trained to predict the electron velocity distribution function (EVDF) following an exponential ansatz. The solution was validated for two reference gas mixtures, whereas the dependence on a variation of the reduced electric field was used as a case study. Excellent agreement between the data-driven approach and MC simulations was reported. Several improvements to their initial scheme have been suggested by the authors. Specifically, a nondimensionalization of the equation combined with an improved ANN architecture was used.⁵¹ Last is a fully connected ANN, which includes two transformer layers for an improved training of the PINN.^38,52 Furthermore, a scheme was proposed to take into account variable reduced electric field values $E/N$.⁵³ In both studies, overall good agreement with MC simulations was demonstrated. For Ar gas eight processes were included: momentum transfer, six electronic excitations, and ionization. For ${\mathrm{SF}}_6$ gas nine processes were included: momentum transfer, vibrational and electronic excitations, electron attachment, and ionization. An extension to the time-dependent Boltzmann equation and assessment of the memory capacity of the ANN of 0.1% compared to mesh-based direct numerical simulation (DNS) calculations of EVDFs in three-dimensional (3D) velocity space was further argued.

In the frame of PINNs for LTP, Zhong et al.⁵⁴ initially investigated stationary and transient 1D arc discharges. A fully connected ANN was used to implement the stationary Elenbaas–Heller equation (with decay term) as well as the time-dependent Elenbaas–Heller equation (considering mass and energy conservation). Whereas the PDE was generally solved using PINNs, the thermodynamic, transport, and radiation properties of the plasma were (partially) taken into account via so-called coefficient subnets. These resemble small regression networks to interpolate the corresponding plasma properties. In a follow-up study,⁵⁵ the authors have introduced this framework more generally in the frame of coefficient-subnet physics-informed neural network and Runge–Kutta physics-informed neural network. For the first approach, a comparison with cubic spline interpolation was included. The last approach embeds PINNs in an implicit Runge–Kutta formalism. In a case study for three examples of Boltzmann equation (with ionization, electron attachment, elastic collision, and excitation), Poisson-drift-diffusion equations (for electrons and ions), and the previously addressed 1D arc discharge using the time-dependent Elenbaas–Heller equation, excellent agreement was noted. Zhong et al.⁵⁶ lastly revisited the data-driven 1D arc discharge simulation, introducing the concept of meta-learning.⁵⁷ Therein, a meta loss was aggregated over multiple tasks (e.g., parameters governing discharge conditions, time-dependence) to enable an efficient two-step optimization and, therefore, reduce the training time. For each training task, a regular PINN was trained. In contrast, the meta network was trained to provide good initial weights for a PINN solving a new task. Based on a variation of arc discharge parameters and corresponding meta-training, a substantial acceleration with a factor in the range of 1.1 to 6.9 was reported.

A different approach to PINNs following Raissi et al.³³ involves the discovery of PDEs from observations. A single PDE or a set of PDEs may be formulated as prototype equations. Therein, unknown coefficients ${\lambda }_i$ determine the exact form of the equations. The problem can be phrased such that training a PINN not only yields the desired solution to the PDE but also includes the coefficients ${\lambda }_i$ as the parameters of the PINN. This concept has been used by Kawaguchi et al.⁵⁸ to “discover” electron transport coefficients from electron swarm maps measured by scanning drift-tube experiments. The temporal and spatial evolution of the electron swarm was formulated in two PDEs with corresponding transport coefficients. A fully connected ANN with transformer layers was utilized similar to previous works by the authors.⁵¹ The tunable parameters are optimized along with the ANN, without assuming any analytical form of the electron swarm map. It was concluded that more accurate transport coefficients could be obtained, including higher-order dependencies if required.

Xiao et al.^59,60 considered a data-driven surrogate model based on recurrent neural networks (RNNs) and a reduced-order model (ROM) for model predictive control (MPC) of a plasma etch process. It is, therefore, addressed more specifically in Sec. 2.4.

2.2.1.

Extraction of cross sections and reaction rate coefficients

The use of an ANN as an optimization technique for treating the inverse problem of obtaining electron collision cross section from electron transport coefficients has been proposed by Morgan.⁶⁴ In that pioneering work, it was assessed that neural networks are indeed useful to determine the elastic momentum transfer cross sections of electrons in real gases, but the accuracy level is highly dependent on the quality of input cross section data that are used for training of the model. Stokes et al.^65,66 improved the idea by Morgan and implemented deep neural network (DNN) using cross sections from the LXCat database⁶⁷ associated with electron transport coefficients found by numerical solutions of the electron Boltzmann equation or by experimental measurements. In their works, the authors demonstrated that their automatic solution using ANN had an accuracy comparable to that of a human expert in determining cross sections of electrons in He or in the biomolecule tetrahydrofuran (THF). The authors also showed that the use of large amount of synthetic training data generated by using the real cross sections available from LXCat provide good results for the prediction of elastic momentum transfer and ionization cross sections of He and Ar.⁶⁶ However, the architecture of the ANN used in Refs. 65 and 66 had minor improvements over the architecture proposed by Morgan.⁶⁴ Jetly and Chaudhury⁶⁸ implemented an ANN, specifically CNN and densely connected convolutional network (DenseNet)³⁷ to derive electron impact cross sections from a solution of the inverse electron swarm problem. The networks were trained using elastic momentum transfer, ionization, and excitation cross sections for different gases available in the LXCat database⁶⁷ and their corresponding electron transport coefficients, which were calculated using the BOLSIG+ solver.⁶⁹ For the first time, the performance and accuracy in cross sections determination from different network architectures was assessed, and it was found that DenseNet predicts cross sections with significantly higher accuracy compared to other ANN models due to its ability to extract both long and short term features from the electron transport parameters. Figure 1 shows that the DenseNet model, as opposed to other ANN-based models (e.g., CNN), is able to predict specific cross sections features in energy, such as the resonance of the elastic momentum transfer cross sections of electrons in ${\mathrm{N}}_2$. Improvements of such deep learning models, such as the ability to extract multiple excitation cross sections and/or the use of sophisticated synthetic data generation schemes, may facilitate the use of ML-derived electron-impact cross sections in plasma chemistry models.

Fig. 1

Predicted total momentum transfer cross sections of electrons with various gas species. Figure reproduced from Ref. 68 licensed under CC BY 4.0. (a) Nitrogen (${\mathrm{N}}_2$), (b) argon (Ar), (c) helium (He), (d) fluorine (F), (e) methane, (f) oxygen (${\mathrm{O}}_2$), and (g) sulfur hexafluoride (${\mathrm{SF}}_6$).

Fast prediction of collision cross sections based on simulated cross sections has been initially suggested by Zhong.⁷⁰ While the author used density functional theory (DFT) calculations to devise a data set for electron impact ionization cross sections for molecules, it has been argued that the approach would be similarly applicable to experimentally measured cross sections. Based on a support vector machine (SVM) model with radial basis function kernel, the prediction of accurate cross sections has been verified. Moreover, the capability to generalize from small to large molecules has been argued, without the need for computationally costly DFT calculations. Following a similar paradigm,⁷¹ Harris and Nepomuceno recently developed a conceptually simple but effective ANN model to predict electron-molecule ionization cross sections from a database of measured cross sections. In particular, by exploiting the chemical configuration (i.e., number of atoms in a molecule) as input to the network, the authors have evaluated the capability of a generalized prediction of cross sections for additional molecules.

Direct determination of reaction rate coefficients and their temperature dependence is also important in plasma modeling. In fact, each reaction needs to be assigned with its own kinetic parameters and these are sometimes not available or subject of large uncertainties. Reiser et al.⁷² used ML to derive rate coefficients for radio-frequency atmospheric pressure plasmas operated with various gas mixtures of ${\mathrm{O}}_2$, ${\mathrm{CH}}_4$, and He. In their work, the authors used a genetic algorithm to derive rate coefficient values that bring the results of the integrated set of particle balance equations for different species in the plasma into good agreement with the experimental data. Even if this approach could be prone to biases induced by user experience and be limited by the availability of experimental data, it has been demonstrated that results from the genetic algorithm open up the possibility of learning more about the significance and correlation of individual reactions. Hence, this method is potentially interesting for post-processing and evaluating experimental data using simplified reaction schemes. Hanicinec et al.⁷³ explored the use of ML to supply unknown reaction rates in plasma chemistries, thus allowing complete chemistry sets to be generated without resorting to estimations by analogy or an educated guess. In their work, ML was used to provide rate coefficients of binary chemical reactions using three distinct optimized regression models: an SVM model, a random forest model, and a gradient-boosted trees model. The models were trained on kinetic data for binary heavy-species collisions at or near room-temperature extracted from the QBD,⁷⁴ KIDA,⁷⁵ NFRI,⁷⁶ and UfDA⁷⁷ databases. After removing duplicate reactions, the final data set consisted of 9470 reactions involving 1080 distinct species. As a sample use case, the ML results were used to augment the chemistry of a ${\mathrm{BCl}}_3/{\mathrm{H}}_2$ gas mixture. In this example, the reactions between the various ${\mathrm{BCl}}_x$ species and H were missing. Results of their work showed the importance of taking such reactions into consideration for which the ML algorithm likely gives better estimates than intuitive guesses while being faster than measurements or ab initio calculations.⁷³

2.2.2.

Derivation of reduced reaction mechanisms

Several methods have been developed to provide interpretation and obtain proper reduction of complex chemical kinetic sets. There are two main methods for reducing reaction mechanisms: reparametrization of the chemical state space and graph theory related algorithms. For the first category, we can refer to principle component analysis (PCA)⁷⁸ and intrinsic low-dimensional manifold (ILDM)⁷⁹ as two methodologies that are used for LTP kinetics reduction. In PCA, a manifold is generated from linearization of the chemical source terms. Reduction is obtained since the chemical state space that is described by a small number of parameters, the principal components (PCs), and balance equations for the PCs are solved instead of the full set of equations. Peerenboom et al.⁸⁰ applied PCA to a 0D state-to-state kinetic model of ${\mathrm{CO}}_2$ including vibrational levels of ${\mathrm{CO}}_2$ and CO. The a priori analysis showed that log-transformation and scaling are important for the manifold identification. Moreover, a significant speedup in central processing unit (CPU) time has been obtained due to the fact that PCA reduced not only the number of variables but also the stiffness of the equations. The ILDM method has been applied by Rehman et al.⁸¹ to Ar and ${\mathrm{H}}_2$ plasma kinetic schemes. It has been shown that the ILDM method automatically extracts the relevant information about fast and slow time scales in a chemical reaction system. After a short interval of time, the fast time scale processes will quickly move toward this low-dimensional manifold and the slow time scale processes will move tangential or along the manifold to finally reach the equilibrium point.⁸¹ Despite successful application in plasma chemical systems, PCA and ILDM also present some drawbacks. For example, in PCA, the manifolds are constructed based on the hypotheses of linearity between original and transformed variables that cannot be always satisfied for complex chemical networks. Moreover, ILDM can be problematic for chemical kinetics that have more than one equilibrium point or for more complex systems where the increase in the dimensionality of the manifolds and subsequently of the generated look-up tables can require higher computational resources.

Graph theory has been applied to chemical kinetics for visualization, topological statistics calculation of a chemical network, and determination of pathways from a sequence of reactions. Sakai et al.⁸² have been the first to apply it to complex plasma chemistry. In their work, they have constructed a network structure for methane and silane plasmas by representing the species as nodes and reactions as directed unweighted edges and they demonstrated that this network analysis allows one to estimate substantial effects of given species on other ones existing in the network. Then, the work was extended with betweenness and centrality indices to highlight not only the influentiality of each species on/to the others but also global structures in chemical networks.⁸³ In this respect, these works have opened up possibilities for simpler prediction of chemical reactions that one would like to understand and/or control. A similar approach for graph construction has been adopted by Sun et al.,⁸⁴ with a different choice of undirected edges weighted based on steady-state reaction rates. The reduction method has been applied to ${\mathrm{CO}}_2$ chemical kinetics for a gliding arc plasma. Murakami and Sakai⁸⁵ have used a graph-based approach for extracting a low-dimensional reaction set based on centrality indices. In particular, their approach was applied to $\mathrm{He}-{\mathrm{O}}_2$ and He-humid air plasma chemistry, and they have shown that a reduced chemical reaction set for these systems can be achieved with the help of topological centrality (closeness and betweenness) and scale-freeness. In particular, in Fig. 2, it is shown that a plasma-induced reaction network can be described with a scale-free network, where the number of degrees originating from a particular species follows a power-law distribution. The advantage of their approach is that it can be used as preprocessing to reduce complex plasma chemistry before multi-dimensional simulations. Mizui et al.⁸⁶ extended the analysis in Ref. 85 to transient states by constructing temporal graphs. Their classification method was used to deduce information from silane and methane plasma chemistries. In particular, quantitative insights into chemical mechanisms have been derived with axes of multi-centrality indices. Holmes et al.⁸⁷ used a directed graph representation of plasma mechanisms in which both species and reactions are identified as nodes, with the edges weighted based on the corresponding reaction rate coefficients. Furthermore, they have used a connectivity matrix to get qualitative information on fast and slow reaction pathways. With this approach, they have shown that the new approach can produce a clearer image for the interpretation of complex chemistry. Venturi et al.⁸⁸ presented a novel method for reducing plasma chemistry mechanisms by taking into account the uncertainties affecting the reaction parameters. As a difference with respect to previous works, their approach relies on an ensemble of weighted directed graphs having both species and reactions as nodes. Edges are weighted based on the reaction contribution to the particle balance equation. The methods were applied to a plasma mixture of 23 species including 15 vibrational states of the ${\mathrm{H}}_2$ molecule, three electronic states of the H atom, four ions, and electrons. As a result of their work, it has been shown that a graph ensemble significantly improves the methodology’s robustness while preserving the predictive accuracy. In general, graph theory related algorithms have been successfully applied to find patterns in large chemical data set. These algorithms could be used to identify and analyze patterns in chemical kinetic data prior running experiments or simulations.⁸⁸

Fig. 2

Power-law distributions of the number of degrees for (a) He + humid air plasma, (b) $\mathrm{He}+{\mathrm{O}}_2$ plasma, 20-species model, and (c) $\mathrm{He}+{\mathrm{O}}_2$ plasma, 12-species model. Figure reproduced from Ref. 85 licensed under CC BY 4.0.

2.2.3.

Chemical engineering

In many applications, the plasma community is interested in knowing “how we can” generate a specific concentration of species. In this respect, modeling and simulations are powerful tools to allow one to actively modulate plasmas and/or discharges for production of chemically reactive species. In this section, advances in ML for chemical engineering of plasma systems are explored. Zhu et al.⁸⁹ developed a DL framework called DeePlasKin for non-equilibrium plasma systems that can be described by a global chemistry model. The model is based on a predictor-corrector approach detailed in Sec. 2.4. For a given kinetic scheme and pre-defined temporal evolution of target species densities, it allows one to extract the temporal profile of the reduced electric field ($E/N$) and all other species densities. Pan et al.⁹⁰ investigated a multi-layer feed-forward neural network for plasma kinetic modeling. The deep learning model has been applied to kinetic schemes for ${\mathrm{CH}}_4-\mathrm{Ar}$ pulsed discharge and a ${\mathrm{N}}_2-{\mathrm{H}}_2$ plasma catalysis model for ammonia production. The magnitude of $E/N$ and time were chosen as the input data of the DNN, whereas the output was selected as the target densities of each species calculated at every instance of time by inquiring the kinetic model. The results showed that the network can replace complex kinetic models with large reactions set while significantly improving the computational efficiency. Wang et al.⁴⁹ proposed an MLP to describe the discharge characteristics and plasma chemistry of ${\mathrm{CO}}_2$ pulsed DBD discharges at atmospheric pressure. The MLP was trained using data obtained from fluid simulations. As a result of their work, the trained model takes a few seconds to predict various features of ${\mathrm{CO}}_2$ pulsed discharges, which significantly improves the computational efficiency with respect to the 30 h CPU time that are required for the fluid model. Overall, these works have shown the great potential of ML for optimization and design of plasma sources in practical applications. In fact, the proposed models have been able not only to automatically extract plasma parameters that are difficult to measure in experiments but also possess high computational efficiency, which provides a promising approach for multi-dimensional plasma modeling.

2.3.1.

Surface chemical reactions

In addition to the species and surface chemical reactions inherent to ordinary gas-phase processing (e.g., chemical vapor deposition, atomic layer deposition, catalysis), plasmas introduce radicals, excited species, and ions and hereby enable a manifold of new reactions at the surface. This circumstance is employed for instance for the plasma enhanced atomic layer deposition (PEALD) of metal oxides. It consists of two self-limiting half cycles with purging cycles in between, out of which the second half cycle, the plasma cycle, is dedicated to the deposition of O atoms, molecules, and radicals. Ding et al.⁹² modeled PEALD of ${\mathrm{HfO}}_2$ thin films with a kinetic MC scheme. Reaction path ways and energy barriers were obtained by combining an intensive literature research with DFT calculations. Such simulations were used to assemble a data set to train a Bayesian regularized artificial neural network (BRANN).⁹³ The trained network was used for the prediction of the half cycle process times, which are essential for industrial considerations, as a function of total precursor pressure and surface temperature. The BRANN was constructed by combing an MLP network architecture with a Bayesian regularization, which penalizes large weights to achieve a better generalization. In a follow-up work,⁹⁴ the authors pair their kinetic MC model with a computational fluid dynamics (CFD) process simulation. An RNN was then correspondingly trained to predict the evolution of the surface profile as well as gas-phase kinetics. Recurrent links establish internal states, which allow for the prediction of the temporal dynamics. The issue of vanishing as well as exploding gradients during the training of the RNN was addressed by including an LSTM (cf. Sec. 3.1).

Plasma enhanced atomic layer etching (PEALE) of copper has been suggested by Sheil et al.⁹⁵ to include a half cycle of plasma induced oxidation followed by a half cycle of exposing the surface to etchant molecules (e.g., HCOOCH), which solely removes the just formed oxide layer. While the latter step is naturally self-limiting, the former is not. This finding has been seconded by Xia and Sautet,⁹⁶ who performed MD simulations for a corresponding plasma oxidation of copper. A high-dimensional neural network potential was trained by considering the ReaxFF potential^97,98 and DFT calculations for the initial and on-the-fly data generation, respectively. The ML potential was demonstrated to overcome the limitations of classical as well as reactive interaction potentials imposed by their functional form, which are argued to not allow for accurate descriptions of very short range interatomic distances.

For the particle fluxes from the plasma, not only radicals but molecules with high vibrational excitation are assumed to play an integral role for plasma catalysis. This is assessed although their contribution and cannot be experimentally clarified due to an overlay of effects, for instance, as induced by high electric fields. In contrast, simulations are well suited to provide experimentally unattainable insight. Wan et al.⁹⁹ conducted DFT simulations to investigate the role of electric field-dipole interactions for the ${\mathrm{CH}}_3$ synthesis catalyzed with Ru surfaces. The resource intensive computations were accelerated by approximately five orders of magnitudes by generalizing the aggregated data with ANNs. Relevant atom configurations were processed with graph neural networks (GNNs),^100,101 specifically graph convolutional networks,¹⁰² whose network structure consists of nodes (i.e., atoms) and edges (i.e., configurational similarities). Multiple shared MLPs were used to combine the output with information on the electric field, adsorbate energy, and surface energy to predict all quantities of interest (i.e., adsorption energies, dipole moments, and polarizabilities of reactive intermediates). Kedalo et al.¹⁰³ performed MD simulations to show that the activation of highly excited ${\mathrm{N}}_2$ molecules on Ru surfaces during ${\mathrm{NH}}_3$ synthesis can be approximated by the Fridman–Macheret $\alpha$ model.¹⁰⁴ The interatomic interactions were described by an ML potential, utilizing smooth overlap of the atomic positions (SOAP) descriptors,¹⁰⁵ to enable ab initio accuracy but with increased computational efficiency. Bal and Neyts¹⁰⁶ demonstrated that the Fridman–Macheret $\alpha$ model is insufficient to describe the vibrational nonequilibrium (i.e., lowering of the free energy barrier) during the catalytic dissociation of ${\mathrm{H}}_2$ and ${\mathrm{CH}}_4$ on Ni surfaces at low temperatures of the order of 750 K. The multitudes of potential reaction pathways in between two metastable states were simplified by applying a harmonic linear discriminant analysis as a supervised dimensionality reduction algorithm.¹⁰⁷ The states themselves were studied by conducting short MD simulations with enhanced sampling methods (i.e., variationally enhanced sampling,¹⁰⁸ metadynamics^109,110) to introduce the vibrational nonequilibrium.

Siron et al.¹¹¹ proposed a high-throughput workflow for the study of amorphous material (i.e., a-C, a-Si, $\mathrm{a}\text{-}{\mathrm{SiO}}_2$, $\mathrm{a}\text{-}{\mathrm{Al}}_2{\mathrm{O}}_3$) surface reactions relevant to dry plasma etching. Surface sites were characterized with SOAP descriptors and clustered with a Bayesian Gaussian mixture model to effectively reduce sites to be studied for reaction kinetics, specifically etching.

2.3.2.

Physical sputtering

The bombardment of surfaces with energetic ions spawns collision cascades in the subsurface region, which eventually lead to persistent radiation damage and the emission of surface material atoms from the surface. The latter is referred to as sputtering and is of major importance for the sputter deposition of thin films as well as etching of high aspect ratio features.

Hamedani et al.¹¹² used an ML interaction potential to conduct an MD study on the radiation damage and sputtering of Si. The utilized ML potential, specifically the Gaussian approximation potential (GAP),¹¹³ was originally trained to equilibrium properties and, hence, complemented with a DFT repulsive potential (i.e., DMol¹¹⁴) for short range interactions to resolve the evolution of the collision cascade accurately. The simulations provided sputtering yields that were in good agreement with experimental reference values. A greater number of isolated small defect clusters was observed than provided by MD simulations based on traditional interaction potentials (i.e., Stillinger–Weber,¹¹⁵ Tersoff¹¹⁶). Moreover, a new defect structure was revealed that is three split interstitials surrounding a central vacancy. This structure and the likelihood for its formation was further analyzed in a follow-up study.¹¹⁷

Sputtering yields can also be predicted by utilizing the Yamamura model.¹¹⁸ Phadke et al.¹¹⁹ shown that its model parameter $Q$, which relates to linear scaling, has a simple positive correlation with the sputter threshold energy. This result was obtained by applying a Bayesian Markov chain MC algorithm to TRIDYN¹²⁰ simulations of ${\mathrm{Ar}}^{+}$ and ${\mathrm{Ne}}^{+}$ ions bombarding 70 elemental targets. Eckstein and Preuss¹²¹ applied Bayesian probability theory to improve the accuracy of the Yamamura model’s sputtering yield predictions by empirical means. Preuss et al.¹²² conducted SDTrimSP^120,123 simulations of ${\mathrm{D}}^{+}$ ions sputtering Fe targets to revise the model’s internal parameter (i.e., surface binding energy). This was achieved by to applying a reduced-order spectral expansion¹²⁴ uncertainty quantification based on the Bayesian framework.

Krüger et al.¹²⁵ conceptualized a PSI ML surrogate model for the ${\mathrm{Ar}}^{+}$ ion bombardment of a TiAl composite target as indicated in Fig. 3. The ML model, specifically an MLP, was meant to map the incident ions’ energy distributions to the energy angular distributions (EADs) of sputtered particles. It was trained with noisy data from TRIDYN simulations (poor statistics), but demonstrated to be capable of predicting smooth EADs, which were almost congruent to reference TRIDYN simulations with high statistical quality. Gergs et al.¹²⁶ continued those efforts by considering TiAl composite targets with varying compositions, introducing the stoichiometry as a basic system state to the model. A similar data aggregation scheme based on TRIDYN simulations was executed. The ML model was replaced by a combination of a convolutional $\beta$-variational autoencoder ($\beta$-VAE)^127–129 and CNN, defining a reduced dimensional (latent) space for the sputtered particles’ EADs and mapping of the incident ion energy distributions as well as TiAl target composition to this particular latent space, respectively. Predictions and reference TRIDYN simulations agreed with each other, even though the assembled PSI ML surrogate model had merely 0.39% of the MLP’s¹²⁵ degrees of freedom and, hence, a significantly reduced likelihood for overfitting. Kim et al.¹³⁰ pursued a similar but more advanced approach for the sputtering of fluorinated Si surfaces due to ${\mathrm{Ar}}^{+}$ ion bombardment. MD instead of TRIDYN simulations were conducted to obtain information on the emitted (i.e., sputtered, reflected) particles’ EADs as a function of the surface F coverage as well as ${\mathrm{Ar}}^{+}$ ions’ incident angle (but not kinetic energy). The aggregated EADs were used for a deep learning-based ROM (DL-ROM), comparing the performance of different ML architectures [i.e., autoencoder²⁷ (AE), denoising AE,¹³¹ VAE, conditional VAE (CVAE)¹³²] with each other. The CVAE was found to best its competitors. The PSI ML surrogate model was completed by introducing an additional MLP, which mapped the input information (i.e., surface F coverage, ${\mathrm{Ar}}^{+}$ ions’ incident angle) to the latent space established by the DL-ROM, enabling resource efficient predictions with atomic, MD fidelity. Gergs et al.¹³³ performed hybrid MD/time-stamped force-bias Monte Carlo (tfMC)^134–136 simulations to setup a PSI ML surrogate model but treated sputtering of Al targets due to ${\mathrm{Ar}}^{+}$ ion bombardment and growth of Al thin films at the substrate as generalized wall interactions. The incident particle flux composition as well as ${\mathrm{Ar}}^{+}$ ion energy was varied to setup the data set.¹³⁷ Two CVAEs were stacked to setup a physics-separating artificial neural network (PSNN), separating the ion bombardment induced damage from the damage to surface state (e.g., composition, mass density) translation. This network architecture has been argued to allow for training with experimental and simulation data in spite of an inherent discrepancy regarding the information accessibility (damage cannot be measured by similar means). The predictions mitigated the noise but maintained key physical features (e.g., threshold for Ar clustering, reflection probability) for small particle doses. In a follow-up work,¹³⁸ the authors have shared a methodology that enabled predictions for particle emission as well as growth during the reactive sputter deposition of AlN thin films in $\mathrm{Ar}/{\mathrm{N}}_2$ plasmas on experimental timescales, that is, minutes to hours. First, the MD interaction potential was revised with an evolution strategy implemented in the genetic algorithm-based reactive force field optimizer method (GARFfield)¹³⁹ to guarantee accurate PSI simulations.¹⁴⁰ Second, a high-throughput, randomized data generation scheme was pursued to efficiently populate the relevant parameter space with hybrid MD/tfMC simulations.¹³⁸ Third, two CVAE were trained to predict either the PSIs or the diffusion processes. They were combined to form a PSNN¹³³ as PSI ML surrogate model, which can readily be used to complement either plasma simulations or experimental diagnostics. The latter scenario was considered to successfully validate the model. Process times of 45 min were predicted within merely 34 graphics processing unit (GPU) hours. In contrast, hybrid MD/tfMC simulations were argued to take more than $\sim 8\text{million}$ CPU years to finish on a comparable case study.

Fig. 3

Conceptual diagram of the ANN plasma-surface interface model and the data flow for training set generation and run-time model evaluation. Figure reproduced from Ref. 125 licensed under CC BY 4.0.

2.4.1.

Model predictive control

In a series of contributions, Woelfel et al. developed MPC schemes for radio-frequency reactive magnetron sputtering fundamentally based on the Berg model.¹⁴¹ Following the identification of an ROM taking into account the nonlinear process dynamics due to surface poisoning,¹⁴² they have devised a controller for process operation at the transition between poisoned and clean surface modes using the previous reduced model, combined with measurements from multipole resonance probe experiments.^143,144 An improved reduced model and novel control scheme (also based on the Berg model) has been introduced later, taking into account the plasma state and its dynamical behavior.¹⁴⁵

The control problem of delivering a specific thermals dose during application of a kHz atmospheric pressure plasma jet (APPJ) operated in helium has been detailed by Mesbah et al. in a sequence of publications discussed in the following. Subsequent to model identification of an ROM (from experimental measurements), a feedback control strategy for 2D spatial thermal dose delivery was proposed.^146,147 The substrate temperature as obtained from infrared thermal imaging was controlled based on the He flow and lateral shuttle position over a surface. The hierarchical feedback control strategy enabled reliable dose delivery despite abrupt disturbances that were introduced. An extension of the approach was suggested as treatment protocol, where the control of the surface temperature included the treatment time in an optimal control problem. One aspect was the adaptation of the control scheme to converge to the optimal solution for the real APPJ system.¹⁴⁸

The optimization of a plasma etch process has been investigated by Xiao et al. using a MPC scheme based on a physical multi-scale model pairing a 2D fluid description of the discharge and a kinetic MC simulation of the surface kinetics. Initially, based on the simulated fluxes at the substrate surface, an RNN surrogate model was set up. Another one was implemented for the surface description, taking into account the stochasticity of the MC model.⁵⁹ Subsequently, a similar physical modeling approach has been proposed for a hierarchical surrogate model strategy.⁶⁰ Using similar 2D plasma fluid simulation results as for the previous approach, an ROM was devised via proper orthogonal decomposition and Galerkin’s method.¹⁴⁹ This reduced model was argued to capture the intrinsic system dynamics on a low-dimensional representation. It thereafter was used for training another RNN surrogate model for rapid prediction of the plasma discharge dynamics during argon-chlorine etching of silicon. In both works by Xiao et al.,^59,60 MPC schemes have further been developed to control the etching depth and the feature bottom roughness of the process.

2.4.2.

Plasma state estimation

With the goal to estimate the internal state of an LTP system in real time despite limited access to observable quantities, a physically constrained extended Kalman filter (PC-EKF) has been devised by Greve et al.¹⁵⁰ First, the general principle including physical constraints (e.g., prohibiting negative quantities) was assessed at the example of the Lorenz system. The PC-EKF has been shown to closely resemble its chaotic system dynamics. Thereafter, the dynamics of a Hall effect thruster described by a global ionization model were considered. Based on particle balance equations for electrons and ions, the electron temperature was dynamically estimated based on the measured discharge currents. Successful recovery of the electron temperature was demonstrated by the investigation of four different discharge modes subject to varying types of ionization oscillations. The proposed scheme has been argued to benefit from improved robustness with the ability to estimate low-frequency dynamics in the presence of possibly not well resolved high-frequency dynamics (in the model or experiments). In a follow-up study also by Greve et al.,¹⁵¹ the PC-EKF approach was generalized to more complex Ar and $\mathrm{Ar}/{\mathrm{O}}_2$ chemical reaction sets and global models of ICP sources. Particle balance equations as well as a power balance equation for electrons were taken into account. In an extensive case study, the robustness of the scheme has been demonstrated. Among other cases, a real time estimate of one or multiple unknown state variables based on corresponding limited measurements (e.g., electron absorbed input power based on argon ion density) has been shown, including pulsed operation.

2.4.3.

Data-driven process recipe design and optimization

In contrast to real time process control and estimation, a study by Kanarik et al.¹⁵² has investigated ML aided process recipe design. A virtual process of a radio-frequency etch plasma with fluorocarbon and ${\mathrm{O}}_2$ chemistry coupled to a feature profile simulator to estimate etching and critical dimension (CD) has been used for reference. A schematic of the virtual process is depicted in Fig. 4. A comparison of junior and senior human experts against three computer algorithms to design a process recipe was evaluated, focusing on a minimization of the accumulated cost-to-target for achieving design goals. Their GPR model was found to greatly outperform the other algorithms (MCMC sampling, tree-structured Parzen estimator). Notably, a threshold was observed, determining the ideal handover from human to computer algorithm design for the specific challenge. This is linked to the number of trials (and hence data samples) required for achieving the specified target. The proposed human first–computer last strategy has been argued to reliably reduce the cost-to-target, whereas human experts and data-driven computer algorithms excel at different tasks (high versus low-dimensional parameter space exploration).

Fig. 4

Schematic of the virtual process used for the process recipe design challenge. Figure reproduced from Ref. 152 licensed under CC BY 4.0.

2.4.4.

Data-driven parameter space exploration and optimization

Data-driven parameter space exploration in LTP has been predominantly pursued for atmospheric plasmas. As highlighted in Sec. 2.2, Zhu et al.⁸⁹ investigated a global chemistry model (with the possibility to interface with a Boltzmann solver) for a DBD as a function of the electric field. A discrete time-stepping predictor/corrector scheme has been suggested: (a) a reconstruction step to obtain a guess of the next reduced electric field $E/N$ value for a selected target density (at the next time step). (b) A prediction step using an ordinary differential equation (ODE) solver of the physical model and the estimated $E/N$ to obtain the future system state. (c) A correction step where the procedure is reevaluated if the difference between the outcome and the target exceeds a given threshold. The model enabled tailoring of the LTP system state evolution, allowing for efficient parameter space exploration and control of the plasma chemistry.

A more general approach to parameter space exploration has been suggested by Shao et al.²² based on GPR. The authors have trained an MLP as a surrogate model based on experimental data. Following the paradigm of Bayesian optimization (BO), this data-driven model is then exploited to set up a GPR model with an expected improvement acquisition function. For a considered use case of nitrogen fixation experiments of an atmospheric DC pin-to-pin glow discharge, for example, the energy cost of ${\mathrm{NO}}_x$ generation was targeted for minimization. An active learning scheme based on this BO was established for exploration of a multi-dimensional parameter space.

3.1.1.

Artificial neural network topologies

Several approaches to the modeling of dynamical systems using ANNs have been proposed. Early concepts have evolved around the class of RNNs.¹⁵³ Recurrent network connections establish feedback loops between neurons, whereas the previous output of a neuron subject to a time delay is fed back to the same or another neuron. Feed-forward RNNs only maintain local recurrent connections per layer; fully recurrent network architectures establish a complete coupling between all neurons. Similar to ordinary fully connected ANNs, recurrent connections may be multiplied with learnable weight parameters and passed through nonlinear activation functions. Training of an RNN involves the backpropagation of errors through time, passing them on through the evolution of the system.¹⁵³ This might be challenging due to vanishing or exploding gradients and a large computational effort. The RNN may finally resemble the system dynamics, through its (possibly many) fitted weight parameters. Examples in LTP processing have already been addressed.^59,60

An important variant of RNNs are LSTM networks suggested by Schmidhuber et al.^154,155 While similar in terms of recurrent feedback, their success is largely enabled by memory cells and gate units instead of basic neurons. Each memory cell with a self-connection includes gate units to control the information flow into and out of the cell, as well as of the recurrent feedback. As such, this may be attributed to the concepts of learning and forgetting¹⁵⁵ and allows for a substantially improved learning efficiency and robustness. LSTM networks have been applied in fluid dynamics for the prediction of the temporal dynamics of turbulent flow through channels based on DNS data.¹⁵⁶ They have further been used for the prediction of disruption in magnetic confinement fusion plasmas.^157,158

Moreover, in a study on a hierarchical approach to multiscale data-driven modeling, LSTM networks have been included for comparison.¹⁵⁹ Another fundamentally related variant of RNNs that was included in this study are echo-state networks (ESN).^153,160 Therein, an RNN network reservoir is initially set up. A number of inputs are fed into the RNN reservoir, but a number of outputs are extracted from the reservoir. In contrast to ordinary RNN networks, which are comprised of learnable weight parameters at every neural node, only a subset of output weight parameters is adjusted during training. Consequently, a manifold of complex nonlinear system dynamics may inherently be evolved by the network, whereas only a small number of output weight parameters capture the desired system dynamics to be described by the data-driven model. To this end, a physics-informed ESN architecture has been investigated for the description of chaotic systems.¹⁶¹ For the chaotic dynamics of a Lorenz system and a Charney–DeVore system, robust training was achieved with an extended prediction time horizon of approximately two Lyapunov times.

Following the physics-informed learning paradigm,³³ the number and variants of approaches exceeds the frame of what can be covered here. The interested reader is referred to an extensive review on PINNs in scientific ML by Cuomo et al. (which also partially covers RNNs and its variants).¹⁶² It may be noted that the technical implementation of PINNs is subject to ongoing research. This is partially due to the limited performance and convergence characteristics of PINNs in general and of ordinary ANNs such as an MLP in particular.^52,163 Another approach to physics-informed ML extending this idea is based on deep neural operator networks (DeepONets).¹⁶⁴ With the reasoning that ANNs can be considered universal operator approximators (similar to the property that ANNs are universal function approximators), these networks comprise a trunk network and a (single or stack of) branch network(s). The branch part takes as input discrete sensor samples $u(x_i)$ of a function $u(x)$. The trunk part takes as input the quantity $y$. The output of both sub-networks is combined to provide a functional representation $G(u)(y)$. Following an offline training for a variety of examples, Lu et al.¹⁶⁴ demonstrated online inference performance significantly improving on other reference ANNs (among others LSTM networks). Different variants and extension taking advantage of Fourier and Laplace transformations as nonlinear kernels have further been proposed.^165–167

Transformer architectures have been originally developed in NLP.³⁸ The examples previously addressed have fundamentally explored similar architectures without substantial conceptual modifications (e.g., either by using standalone transformer layers^51,53,58 or similar to an image-to-image translation task³⁶). Classical transformer networks process input information through embeddings that transform the input data into a corresponding vector space representation. Recently, an extension to the concept of transformer networks has been suggested by Geneva and Zabaras for surrogate modeling of physical systems.¹⁶⁸ Therein, alternative embeddings have been suggested based on Koopman dynamics in a two step learning procedure. Initially, embedding ANNs are trained in an encoder–decoder structure to approximate the Koopman based dynamics in a latent space representation. Thereafter, a transformer network is trained using the pretrained Koopman embeddings. The obtained data-driven surrogate models have been demonstrated to effectively capture the chaotic system dynamics of a Lorenz system, the 2D fluid dynamics of a transient flow, and 3D reaction-diffusion dynamics. In comparison, the Koopman based transformer models have been demonstrated to outperform alternative surrogate models of such dynamical physical systems (including other embeddings and LSTM networks).

3.1.2.

Reduced order modeling and model identification

ROMs have previously been used in LTP modeling and simulation as exemplified above. The particular notion of sparse identification of nonlinear dynamics (SINDy)^169,170 has been demonstrated to provide interpretable and generalizable nonlinear differential equations of reduced order, which effectively captures the dynamics of the full model. This is achieved by identifying the relevant system dynamics from a library of potential functional forms by means of a sparse regression. It has been applied to space and astrophysical plasmas based on PIC simulations by Alves and Fiuza demonstrating the potential of recovering the fundamental dynamics of the high-fidelity kinetic plasma simulations.¹⁷¹ It may be straightforwardly applicable to LTP modeling and simulation data to capture their nonlinear dynamics, albeit a complex chemistry may necessitate further adaptation. Sparse regression has similarly been applied by Kaptanoglu et al. in the context of magnetic confinement fusion plasmas for the optimization of magnetic fields in stellarator fusion plasmas.^172,173 Moreover, MPC based on such ROMs has been suggested with limited availability of training data.¹⁷⁴ An aspect that has been applied using time-dependent simulations of the tokamak plasma boundary.^173,175 It may be specifically relevant also in the context of LTP processing with its diverse requirements for precise process control.

3.1.3.

Learning schemes

In contrast to classical supervised learning schemes, alternative approaches may be relevant also due to limited availability of training data. In a Bayesian setting, provided with an estimate of the prediction uncertainty, active learning has been exploited for efficient parameter space exploration.¹⁷⁶ While such an uncertainty metric is straightforward in the context of GPR modeling,²² similar techniques have been adapted to MD simulations and active learning of corresponding ML surrogate models.¹⁷⁷

Transfer learning enables an alternative route to ML in the small-data limit.¹⁷⁸ The fundamental idea relates to the property of generalization that data-driven models pursue. Given an initial problem task with abundant data, an ML model may be robustly trained. The pretrained model may be subsequently retrained on a related, yet different small-data task, whereas just a minor subset of weight parameters is allowed to be adjusted. Hence, the robustness and generalized prediction of the initial model may be successfully transferred to the more challenging learning task. This approach has been successfully applied for instance in the design optimization in inertial confinement fusion.¹⁷⁹

3.1.4.

Explainable AI

As surveyed by Gilpin et al.,¹⁸⁰ a multitude of explainable AI (XAI) concepts have been proposed for the interpretation and explanation of data-driven models. As an example, local interpretable model-agnostic explanations (LIME) have been suggested for identifying interpretable models that locally approximate the nonlinear behavior of any underlying model.¹⁸¹ While this is possibly an underexplored aspect in the context of LTP, an example where LIME relates to the data-driven classification and interpretation of OES experiments of plasma in aqueous solution has been proposed.¹⁸² Further works by Lundberg et al.¹⁸³ have included LIME and other methods in an extension Shapley additive explanations. It is similarly based on local explanation models but defines additive feature attribution methods to unify several of the included interpretable methods. There is arguably no mathematical proof that physics-based and physics-informed ML and corresponding models capture exactly the desired dynamics (and only those). LTP modeling and simulation has got its inherent foundation in a rigorous mathematical description and numerical solution. XAI methods may not only foster advances in data-driven approaches to the field but also strongly contribute to the general acceptance in LTP science and technology.

3.2.1.

Data-driven surrogate sub-modeling

Describing the evolution of dynamical systems (or sub-systems) accurately is an omnipresent task in scientific computing. The time integration schemes of differential operators have been correspondingly considered in different fields (e.g., pairing Runge–Kutta integration schemes with ANNs^{55,159,184,185}). For the time integration of multiscale physical systems, Liu et al.¹⁵⁹ suggested a hierarchical time-steppers (HiTSs) framework. The approach builds on the decomposition of the system dynamics into a hierarchy of prediction steps that can be arranged and interpolated in time. It has been shown to outperform other methods based on LSTM, ESN, and clockwork RNNs for the Lorenz system, 2D fluid flow past a cylinder, the Kuramoto–Sivashinsky (KS) equation, as well as music and video sequences. An illustrative example comparison of the prediction output is shown in Fig. 5. LTP modeling and simulation is facing essentially alike problems, suggesting it may similarly benefit from application of the HiTSs method (and/or other methods referenced) in the description of their multiscale dynamics.

Fig. 5

Output of different ANN architectures for the prediction of different training sequences (Kuramoto–Sivashinsky equation, music, flow past a cylinder, video). Figure reproduced from Ref. 159 licensed under CC BY 4.0.

An approach using the concept of dynamic mode decomposition (DMD) to obtain low-dimensional linear dynamics of given data has been applied to magnetic confinement fusion plasma simulations by De Pascuale et al.¹⁸⁶ The obtained dynamic modes are evolved individually solving a linear ODE. The constructed linear time advance operator based on DMD has been successfully applied for the analysis of several 1D diffusion transport problems, as well as for the dynamics of 1D profiles at the divertor targets of a tokamak plasma. This scheme could be similarly applied to complex LTP processes.

The solution to Newton’s equations for particle based simulations has been explored by Kadupitiya et al. using LSTM networks.¹⁸⁷ The accurate applicability has been investigated for few-particle systems subject to various interaction potentials (e.g., Lennard–Jones). Energy-conserving dynamics up to $4000\times$ larger timesteps compared to a Verlet integration have been demonstrated. Following the concept of GNNs,^100,101 which embed the interaction dynamics of nodes (e.g., particles) in their learnable edge weights (e.g., interaction), the integration of ODEs in a Hamiltonian framework have been discussed by Sanchez-Gonzalez et al.^188,189 The proposed Hamiltonian ODE graph network trained on Runge–Kutta scheme data generalizes the integrators to accurately represent the Hamiltonian dynamics. In the context of many-particle simulations, Li et al.¹⁹⁰ proposed a GNN accelerated MD framework that predicts the interaction forces for the system, given its momentary state (e.g., atom positions and types). It has been argued that the scheme can be easily integrated into existing numerical integrators. Due to the consideration of multi-particle interactions, this scheme may be specifically relevant in PIC simulations should Coulomb interactions be important.

3.2.2.

Data-driven discharge surrogate modeling

Many of the ML methods that may be transferred to data driven surrogate modeling of LTPs were discussed in Sec. 3.1. This relates specifically to the representation of dynamical systems through recurrent ANN models. A few examples from plasma science and adjacent fields were covered. However, it goes beyond the scope of this paper to cover related works from fluid dynamics. The reader is referred to corresponding comprehensive surveys by Brunton and Vinuesa et al.^191,192

Concerning data-driven process control and MPC, it should be noted that the plasma information based VM methodology for LTP etching processes previously introduced may be adapted based on physical or data-driven model data.^19,20,193 Despite a possible discrepancy between models and reality, in particular systematic deviations may be straightforwardly taken into account in a data-driven approach by means of discrepancy learning. A learnable translation layer may suffice to capture and combine the relevant information by means of data-fusion. This may enable data-driven MPC with limited observables (cf. also PC-EKF^150,151). In addition, MPC schemes based on LSTM networks that have been successfully demonstrated in magnetic confinement fusion plasmas subject to inherent complex nonlinear dynamics by Degrave et al.¹⁹⁴ may be promising candidates for adaptation in LTP processing.

3.3.1.

Reduction of complex chemical networks

The adoption of detailed mechanisms for chemical kinetics often poses two types of challenges. First, the number of degrees of freedom is large; and second, the dynamics is characterized by a wide range of time scales. This has motivated the development of several techniques for reducing the complexity of such kinetic models where only a few variables are considered in the simplified model. Grassi et al.¹⁹⁵ used AEs to reduce complex chemical data sets and decrease the computational time associated with numerical solutions of ODEs for a large number of species and reactions. In particular, an encoder is used to reduce the dimensionality of the chemical state space into a latent (lower-dimensional) space. Analogously to the chemical network in the original space, which is represented by the system of ODEs, the chemical network in the latent space can be evolved in time with a different set of ODEs, which has got a significantly smaller number of dimensions. Hence, a decoder is used to transform the variable from the latent space back to the original space. The proposed model enables a compression of chemical networks composed of 224 reactions and 29 species into a network with 12 reactions and 5 species that are evolved in time using a standard ODE solver. Their work has demonstrated that a 65-times speed-up can be achieved. Tang and Turk¹⁹⁶ proposed a combination of AE and neural ODE to model the temporal evolution of chemical kinetics in a reduced subspace. The architecture of the network is similar to the one by Grassi et al..¹⁹⁵ In their work, a 10-fold speed-up compared to commonly used astro-chemistry solver for a 9-species primordial network has been achieved while maintaining 1% accuracy across a wide range of density and temperature. As described in Sec. 2.2, it is well known that, due to the presence of fast and slow dynamics, the chemical systems are characterized by low-dimensional manifolds in the concentration space. Chiavazzo et al.¹⁹⁷ applied the diffusion Map (DMAP) approach for constructing reduced kinetics models for combustion applications. The DMAP approach can be considered as a non-linear counterpart of PCA that can be used for searching a low-dimensional embedding of a high-dimensional chemical data set.

3.3.2.

Fast and accurate models of complex chemistries

Stiff equations for chemical kinetics are computationally expensive to solve numerically. The efficacy of ML for applications that need to solve systems of ODEs or PDEs has been explored in many areas, such as CFD and combustion. Goswami et al.¹⁹⁸ developed a neural operator-based surrogate model (DeepONet) to efficiently solve stiff chemical kinetics. The model was applied to a chemical reaction system for $\mathrm{CO}/{\mathrm{H}}_2$ burning of syngas, which contains 11 species and 21 reactions, and a temporally developing planar $\mathrm{CO}/{\mathrm{H}}_2$ jet flame (turbulent flame) using the same syngas mechanism. It has been shown that the DeepONet, once trained, can accurately integrate the thermochemical state for arbitrarily large time advancements, leading to significant computational gains compared to stiff integration schemes. Moreover, both results obtained with DeepONet and AE-based DeepONet are computationally very efficient compared to conventional CFD solvers. One of the advantages of this approach is the possibility to be included in a CFD solver for solving such stiff chemically reacting problems, which can drastically reduce the computational cost. Campoli et al.¹⁹⁹ performed state-to-state numerical simulations of high speed reacting gas flow using data-driven ML regression models. The models were applied to a system of equations for a 1D reacting flow of a five-component air mixture including a total of 122 excited states of ${\mathrm{N}}_2$, ${\mathrm{O}}_2$, NO, N, and O species with detailed vibrational kinetics. Several methods have been compared for estimation of the state-to-state relaxation rate terms, such as kernel ridge (KR), SVM, k-nearest neighbor, Gaussian processes, decision tree (DT), random forest, extremely randomized trees, gradient boosting, histogram-based gradient boosting, and MLP. From their comparison, it has been found that KR reports the lowest error levels while the SVM reports the highest. Comparable error levels are reported by the remaining algorithms; however, there are noticeable differences in the prediction time, where DT is found to be the fastest algorithm. Finally, the solution of the state-to-state Euler system of equations was inferred by means of a DNN by-passing the use of the solver. As a result of this study, it has been demonstrated that embedding of ML algorithms into ODE solvers offers a speed-up of several orders of magnitude. Chemical reactor networks (CRNs) have been investigated by Savarese et al.²⁰⁰ for performing faster CFD simulations of combustion systems. In their work, a novel automatic data-driven method for the design of CRN models has been proposed. The method is based on a combination of unsupervised clustering and graph scanning algorithms. Results of their work show that the novel methodology is capable of extracting equivalent CRN models of realistic combustion devices from CFD data. Furthermore, the unsupervised clustering algorithm was able to automatically detect important regions within the combustor domain, such as inlets, inflame, and post-flame reactors. Overall, the approach is promising as it has the potential to drive toward the development of ROMs without a huge amount of user-based knowledge. Mao et al.²⁰¹ developed a novel CFD solver that combines the multi-functionalities of the OpenFOAM library²⁰² with the ML framework Torch²⁰³ and chemical kinetics program Cantera.²⁰⁴ In their solver, DNNs are used to accelerate chemistry solvers and improve the simulation efficiency. As a result, a speed-up of two orders of magnitude is achieved in a simple hydrogen ignition case when performed on a medium-end GPU. This work highlights the potential of the integration of ML techniques within current libraries and computing architecture for detailed CFD simulations of complex hydrocarbon fuels.