{'fraction_non_alphanumeric': 0.04093979441997063, 'fraction_numerical': 0.012834067547723935, 'mean_word_length': 4.5188006482982175, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Physics-separating artificial neural networks for predicting initial stages of Al sputtering and thin film deposition in Ar plasma discharges

Simulations of Al thin film sputter depositions rely on accurate plasma and surface interaction models. Establishing the latter commonly requires a higher level of abstraction and means to dismiss the fundamental atomic fidelity. Previous works on sputtering processes addressed this issue by establishing machine learning surrogate models, which include a basic surface state (i.e., stoichiometry) as static input. In this work, an evolving surface state and defect structure are introduced to jointly describe sputtering and growth with physics-separating artificial neural networks. The data describing the plasma-surface interactions stem from hybrid reactive molecular dynamics/time-stamped force bias Monte Carlo simulations of Al neutrals and Ar$^+$ ions impinging onto Al(001) surfaces. It is demonstrated that the fundamental processes are comprehensively described by taking the surface state as well as defect structure into account. Hence, a machine learning plasma-surface interaction surrogate model is established that resolves the inherent kinetics with high physical fidelity. The resulting model is not restricted to input from modeling and simulation, but may similarly be applied to experimental input data.


I. INTRODUCTION
Physical vapor deposition is commonly used for the sputter deposition of thin films. Various processes and types of gas discharges may be exploited, whereas a direct current or radio frequency voltage of a few 100 V to kV are typically used to sustain a plasma. Inert or reactive gases may be used depending on the application and the desired thin film. Sputtering processes have in common that the strong electric field in the plasma sheaths accelerates ions (e.g., Ar + ) towards the electrodes. [1][2][3][4] At the target surface, the ion bombardment leads to a collision cascade and eventually emission of material. This is referred to as sputtering and be can be studied theoretically by performing transport of ions in matter (TRIM), Monte Carlo (MC), or molecular dynamics (MD) based simulations. After release of the sputtered atoms from the surface with the characteristic Sigmund-Thompson distribution [5][6][7][8][9][10][11][12][13], they are subject to collisional transport through the plasma [2,[14][15][16]. Collisions with electrons may ionize the otherwise neutral atoms. Transport phenomena can be resolved by either fluid or particle in cell/Monte Carlo collision (PIC/MCC) models coupled to test-particle methods (TPM) or direct simulation Monte Carlo (DSMC) models [14,[16][17][18][19][20]. When the sputtered atoms arrive at the opposing substrate surface, they contribute to thin film growth.
The different length and time scales of the solid and the gas/plasma states of matter render a joint simulation framework with preserved atomic fidelity impractical [22][23][24][25]. Hence, surface and plasma simulations need to be evaluated separately, whereas often a priori evaluated particle fluxes onto the surfaces and effective surface kinetics are considered (e.g. Berg model [21,26], Sigmund-Thompson theory [11][12][13]), look-up tables).
This issue has been addressed in previous works by proposing data-driven machine learning plasma-surface interaction (PSI) surrogate models that are intended to complement plasma simulations [22,[27][28][29][30][31][32][33]. In particular, a multi layer perceptron (MLP) was trained to generalize on data describing the sputtering of a Ti 0.5 Al 0.5 composite target, simulated with TRIDYN [7,22]. TRIDYN was also used in a subsequent study to generate data for Ti 1−x Al x composite targets. The stoichiometry x was introduced to the model as an additional degree of freedom and preliminary surface state. Moreover, the MLP was replaced by a β-variational autoencoder establishing a reduced latent space representations of the sput-tered particle energy angular distribution functions as well as a dedicated mapper network correlating them with the input (i.e., Ar + ion energy distribution functions) [7,27,34,35].
In this work, the concept of a machine learning surface surrogate model is advanced byamong other aspects -exploiting hybrid reactive molecular dynamics/time-stamped forcebiased Monte Carlo (RMD/tfMC) simulation data of sputtering and deposition of Al thin films with Ar working gas [36][37][38][39]. The manuscript is structured as follows: In Section II the setup depicting the surface interactions is introduced. In Section III, the utilized methods and the corresponding parameters are detailed. The results are collected and discussed in Section IV. Finally, conclusions are drawn in Section V.

II. SETUP
The considered plasma-surface interactions relate to a previous simulation study investigating the sputter deposition of Al thin films with Ar working gas, and specifically Ar + ion and Al neutral impingement on an Al(100) surface. The evaluation procedure and physical interpretation is detailed in [39]. In summary, hybrid reactive molecular dynamics/timestamped force-biased Monte Carlo (RMD/tfMC) simulations have been conducted to describe individual surface interactions [36][37][38]. The simulations have been performed using the open-source Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [40,41]. The crucial role of Ar atoms incorporated in the Al surface has been identified for the description of the properties of the deposited thin film (e.g., stress, defect structure). A rigorous treatment of the incorporated Ar atoms has been found to be of particular relevance in this regard. In the following, the key aspects of the simulation setup and the provided data are reviewed: Particle fluxes to the surface The particle fluxes incident from the plasma to the surfaces are defined as follows. Two particle doses consecutively impinge onto an Al(001) surface, each of 1.21 × 10 15 particles/cm 2 . The fluxes consist of five combinations of Al neutral to Ar + ion flux ratios Γ in Al /Γ in Ar + ∈ [0, 1]. 20 different Ar + ion energies E in Ar + are selected from a square root energy axis in the range of 3 to 300 eV. Al neutrals are assumed to be thermalized (k b T = 7 eV) [39]. The particle fluxes from the plasma to the surface are characterized by Surface defect structure The ongoing ion bombardment leads to surface reconstructions due to the accompanying collision cascades and due to point defects (e.g., vacancies, interstitials). The latter result from sputtering (emission) as well as forward sputtering (peening).
The defect structure after the 1st and 2nd dose of incident particles is quantified by a ring statistical connectivity profile [42][43][44]. It describes the systems' bond topology and consists of four properties: i) The number of rings per supercell divided by the number of nodes R C,r for rings with r nodes; ii) The fraction of nodes at the origin of at least one ring with r nodes P N,r ; iii) The fraction of nodes with r sized rings as shortest closed paths P min,r ; iv) The fraction of nodes with r sized rings as longest closed paths P max,r [39,[42][43][44][45]. The defect structure is characterized by the descriptors D = {R C,r , P N,r , P min,r , P max,r } for r ∈ [3,6]. Particle fluxes from the surface Finally, the fluxes of particles emitted from the surface due to reflection or sputtering are characterized by their outgoing energy distribution functions (EDFs) f . These are evaluated by collecting the individual particle trajectories. The details of how the outgoing particle information is processes are presented in Section III A.

III. METHODS
In the following we proceed with an extension of the strategies detailed in previous publications on machine learning surface surrogate models [22,27]. Specifically, after assembling the data set from hybrid RMD/tfMC simulations and defining the relevant descriptors, the data is pre-processed in terms of its data structure and normalization. In this work, data augmentation is added thereafter, as elaborated below. The data set is divided for training, validation, and testing and a physics-separating artificial neural network (ANN) is optimized and trained subsequently. Finally, the model may be assessed given a number of metrics specified below.
A. Data preparation, training and metrics a. Outgoing EDFs pre-processing The ensemble of outgoing particles are characterized by the individual particle properties (especially kinetic energies). Equivalently they can be depicted by their respective EDFs. The EDFs are obtained from phase-space information of the individual particle in two steps. (1) A kernel density estimator (KDE) with tophat kernel is used to obtain a continuous representation of the EDFs based on the particle energies.
KDE is favored instead of simple histogram binning providing a more accurate description despite the statistical challenges. Specifically, the simulation case with the best statistical description, that is the largest number of emitted (i.e., reflected, sputtered) particles, is b. Data normalization All data (input data x and output data y) is normalized to range from zero to one: x → (x − x min )/(x max − x min + ) and y → (y − y min )/(y max − y min + ). The min/max values are extracted from the training data set and stored for the final denormalization. For the EDFs, the maximum value per species is used as a reference.
Notably, this normalization is merely relevant for an optimal training procedure.
c. Data set size and splitting In total n data = 200 individual cases represent respective training samples from a combination of 2 × 5 × 20 varied parameters. The data set is shuffled and split into 80 % training, 10 % validation, and 10 % test set (i.e., n train data = 160, n val data = 20, n test data = 20) for a study of the hyperparameters (HPs) of the ANN and its subsequent training.
d. Data augmentation A virtual (artificial) extension of such a small data set by augmentation with "combined" data samples improves the training robustness in spite of providing no further information (but only a hypothesis) to the network [47][48][49]. In this work, a constrained version of mixup augmentation is proposed. Mixup augmentation is commonly used for image classification, but in this work found to facilitate training a more robust regression model [50]. Augmented inputx = λx i + (1 − λ)x j and corresponding output samplesŷ = λy i + (1 − λ)y j are computed for samples i and j, respectively. The factor λ is sampled from the beta distribution λ ∼ Beta(α, β), whereas α = β is considered as HP.x andŷ are recomputed each epoch [50]. However, sample pairs i and j are only mixed when k, l, and m specify the indices along the incident Ar + ion energy axis E in Ar + ,k , the Al neutral to Ar + ion flux ratio axis Γ in Al,l /Γ in Ar + ,l , and the 1st dose (m = 0) or the 2nd dose (m = 1) of impinging particles, respectively. The augmentation range r c is considered as HP. It allows for a continuous transition from no augmentation (i.e., r c → 0) to ordinary mixup augmentation (i.e., r c → ∞).
Mixing of the data can be interpreted as providing the hypothesis of a linear superposition of sample pairs with a maximum index radius r c to the network. This mixup combination is at most conducted once per sample pair. The training set is then replaced by all combinations of mixup training set sample pairs fulfilling this constraint. This means that for instance for r c → ∞ the number of augmented training set samples is squared (i.e., n train data,aug = n train data 2 ).
e. Training procedure and metrics The mean absolute errors (MAEs) of the model's predictions on each reference training batch are backpropagated through the network to update the network weights (degrees of freedom). The stochastic gradient descent algorithm adaptive moment estimation (Adam) is applied [51]. The batch size n batch is chosen to be as close as possible to equal 32 while fulfilling |n batch − n train data,aug %n batch | ≤ 4. This requirement ensures that all data samples contribute almost equally to the training progress.
The coefficient of determination R 2 is utilized as an additional metric to assess the network's prediction accuracy. A coefficient of determination equal to one resembles a fully explained residual variance.
The learning rate is set to 10 −3 for an initial simulated annealing phase (described in Section III B). Subsequently, it is reduced by a factor of 10 whenever the validation loss is not smaller than its previous minimum value over the course of 10 epochs. Early stopping ends the training phase when there is no further reduction of the validation loss in 25 epochs.
f. Hyperparameter study After training, the validation loss may be used to find a suitable set of HPs as detailed in Section III C. 10-fold Monte Carlo cross validation (MCCV) is applied to address the data splitting induced variance (potential bias) in model training and validation, as well as its test error in prediction. The last is used to asses the model's performance on unknown data (no feedback to the model). Once the final set of HPs is selected, the corresponding network is trained with the validation and augmented training set by performing a 100-fold MCCV. To obtain final training metrics, the MAE of the unnormalized predictions on the three (not augmented) subsets (i.e., training, validation, test) are evaluated to accurately asses the model's fitness.
g. Production run For the final production run, first, the validation set is merged with the test set (i.e., n val data = 40). Second, an ensemble of ANNs is trained with the final set of HPs (see Section III C) to address the data splitting induced variance (potential bias) in the model's training phase [52]. The mean values of the ANN ensemble predictions are finally output. Third, the machine learning model is used to predict the case study described in Section II. The considered input Ar + ion energy range (i.e., 3-300 eV) and incident Al neutral to Ar + ion flux ratios (i.e., 0 to 1) are sampled more continuously (i.e., 1000 steps along each axis). The predicted output surface state S s after the 1st dose is fed as input to the model for the 2nd dose of impinging particles.

B. Physics-separating artificial neural network
In the following, the physical phenomena or processes to be modeled are subdivided into smaller physically disjoint contributions, each of which described by their own individual ANN model. Once established, the combination of submodels form a physics-separating neural network (PSNN) for the intended application. One advantage of such a subdivision is that it allows for the combination of dependent but possibly unavailable sets of descriptors.
The model and its training can be set up to depend solely on available descriptors (e.g., measurable quantities), whereas unavailable descriptors (e.g., non measurable or dependent quantities based on theory) are hidden as internal system states. The PSNN applied in this work follows this idea. It consists of two connected regression machine learning models using conditional variational autoencoders (CVAEs) architectures.
A CVAE maps an input x to an output y [53,54]. The information flow is depicted schematically in Figure 1. Model predictions are marked by primes (e.g., y ).
First, the training phase is described (blue and orange elements in Figure 1). An encoder (blue) that is conditioned on x compresses y|x to a n ls dimensional latent space (ls) with The sampled latent space variables z ls |y, x project the output from the encoder (blue) to the input of the decoder (orange). The decoder is further conditioned on the input x. It finally transforms the inputs z ls |y, x and x and predicts a representation of the output y , trying to reconstruct y. This reconstruction is governed by a reconstruction loss, e.g., given by the mean absolute error (MAE) |y − y |. During the training phase, the network initially resembles an ordinary β-variational autoencoder (β-VAE) whose encoder and decoder are conditioned on the input x.
The data transformation through multivariate normal distributions structures the information from the ANN inputs (x and y) and distributes it over the latent space coordinates z ls . Initially, conditioned projections z ls |y, x are distributed with respective standard deviations and mean values. They are subsequently forced toward a common standard normal distribution due to the consideration of an additional Kullback-Leibler (KL) divergence, which pushes N (µ ls , σ ls ) toward N (0, 1). To achieve this, the KL loss is multiplied with a HP β [34,35]. Simulated annealing is applied to this loss contribution, gradually scaling it with a factor of 10 −3 to 1 over the course of ten epochs (i.e., per batch over a log-axis).
Notably, the encoder (blue elements in Figure 1) is dismissed post training for prediction (e.g., on the validation or test set). The latent space variables z ls are then sampled directly from a standard normal distribution. In the end, the decoder (orange elements in Figure 1) represents the final regression model [53,54].
a. PSNN training and inference The PSNN combines two CVAE decoders for inference as depicted schematically in Figure 2. PSI (i.e., ion bombardment induced sputtering, damage formation) and damage to surface state conversions (DSS; i.e., defect structure related changes of experimentally attainable surface properties) are modeled by the PSI decoder and the DSS decoder, respectively. The reasoning to separate the surface interaction physics in two ANNs is that the defect structure D is not experimentally accessible and, hence, may be hidden as an internal state in the PSNN. This decision is built on the hypothesis that D as well as S s are equivalently sufficient to describe the system state appropriately. Finally, the DSS decoder output S s may be fed to the PSI decoder as input for the next dose of impinging particles, establishing a recurrent connection and assuming that S s describes the system state sufficiently.
Notably, while prediction proceeds utilizing only the CVAE decoder networks, training is conducted with the complete CVAE networks. For the training it is moreover fundamental to combine both CVAEs to a single PSNN. Both CVAEs are then trained simultaneously, whereas the encoder networks are provided with the inputs and outputs as shown in Figure 1.
The MAEs of each output (f , D for PSI and A s , S s for DSS) are backpropagated through the PSNN all the way to its inputs. If reference values for the internal descriptors (i.e., D) would not be available, e.g., for hypothetical experimental data, the MAE for these descriptors may be set to zero, suspending error backpropagation.
b. PSI-CVAE Both sets of descriptors P and S s for the regression task as well as f and D for the autoencoding task are provided to the PSI-CVAE as input. It predicts the ion bombardment induced defect structure D as well as the EDFs of reflected and sputtered particles f . In its encoder, the EDFs f are passed through three convolutional layers (CLs) with kernel sizes of 9, 7, and 5, filter numbers n f (HP) and strides 2 to extract relevant features. These are then concatenated with D, P and S s . They altogether are fed to n hl (HP) hidden dense layers (DLs) with n npl (HP) nodes per layer resembling an MLP. The output layers of the encoder provide the mean values µ ls and standard deviation σ ls . The activation functions for any hidden and output layer are set as rectified linear unit (ReLU) and linear, respectively. The CVAEs are chosen to be symmetric and, hence, the decoder is built in the inverse order. CLs are replaced with upsampling transposed convolutional layers.
A last custom layer enforces the output y of each CVAE to lie between lower and upper bounds y lo and y hi , respectively. If the input to this custom layer is considered as x, the constraint is enforced in the two consecutive computations: i) y tmp = x − 2(x − y lo )(x < y lo ), ii) y = y tmp − 2(y tmp − y hi )(y tmp > y hi ). The lower bounds y lo for all output properties but the biaxial stress σxx+σyy 2 and shear stress τ xy are zero. The upper bounds y hi for the Ar concentration x Ar , the Ar gas porosity φ Ar , and the ring statistical properties P N,r , P min,r , P max,r } are 100 (relative percentage) and 1 (relative values), respectively. All remaining upper bounds are chosen to approach infinity. A detailed figure of the network structure can be found in the Appendix.
c. DSS-CVAE The DSS-CVAE is trained to predict the surface properties A s and surface state S s after the ion bombardment on the basis of the provided (or previously predicted) defect structure D. The structure of the PSI-CVAE is adopted for the DSS-CVAE, but no convolutional layers are utilized. In the encoder, D , A s , and S s are concatenated directly. The decoder mirrors the encoder network with the addition of the previously mentioned custom layer to enforce output bounds. For training, D is provided as input for the regression task whereas A s and S s are provided as input for the autoencoding task. A detailed figure of the DSS-CVAE network structure can be found in the Appendix. Notably, each CVAE is assigned with its own unique set of HPs.

C. Hyperparameter study
The 11 HPs to be optimized are listed in Table I (n HP = 11). This is achieved by applying an anisotropic self-adaptive evolution strategy with intermediate recombination (µ/µ I , λ)-σSA-ES. HPs are represented by individual genes g and step sizes (mutation strengths) σ that altogether form a list that is referred to as chromosome. The effect of the step sizes (mutation strengths) are detailed later. Each individual in the population has got its own chromosome. The overall population size is specified by λ. The value of each gene g is initialized by sampling randomly in between g lo and g hi (integer types are enforced when necessary). The initialization range intervals for all genes are provided in Table I respectively [56][57][58]. Second, each gene g is updated following g → σU(−1, 1)(g hi −g lo )/2+g h .
The HP study is initiated with 70 individuals (λ = 70) and reduced by one per generation until eventually a population size of seven is kept constant. The number of parents is determined by µ = λ/7 (integer values are enforced) [59]. Hence, a (7/7 I , 70)-σSA-ES is gradually turned into a (1, 7)-σSA-ES over the course of 63 generations. The (1, 7)-σSA-ES is then continued for the remaining generations.

IV. RESULTS
First, the outcome of the HP study and the machine learning model's generalization capability is assessed. Second, the results of the production run with the final set of HPs  are discussed concerning the conservation of physical principles observed in the underlying case study [39]. This study focuses on only the most relevant physical properties.

A. Hyperparameter study
For the evolution strategy of HP optimization, the MAE and R 2 evolved for 200 and averaged over the course of 40 generations are shown in Figure 3. The initial deviations between the errors on the validation as well as the test set are removed after 60 generations.
The gap to the errors on the training set is reduced within 140 generations. Both findings indicate an improved generalization. Afterwards, the evolution strategy converges. The final values of all HPs can be found in Table II, whereas the evolution is depicted in Figure 4. The final MAE on the training, the validation, and the test set are presented in Table II for all considered properties. The reference values include statistical noise as well as physical fluctuations (i.e., cycle of Ar incorporation, clustering, outgassing) and, hence, do not resemble any kind of ground truth [39]. This provides an inherent lower limit for the MAE.
For instance, the shear stress is predicted correctly to equal zero for any incident particle flux by the PSNN. The provided MAE of the shear stress τ xy in Table II describes solely the statistical noise. Hence, the resultant errors on the unknown test set are considered to be sufficiently small when comparing them to the underlying data in terms of magnitude, fluctuation, and noise [39]. Reflection of incident Ar + ions and outgassing of incorporated Ar atoms lead to an overall higher Ar flux leaving the surface. Thus, a higher MAE is  observed for outgoing EDFs of Ar f Ar than for Al f Al .

B. Production run
The first outputs of the PSNN ensemble are the EDFs f of reflected as well as sputtered Al and Ar atoms. An example for an Al EDF is shown in Figure 5 for an incident Ar + ion energy and Al neutral to Ar + ion flux ratio of 197 eV and 0.42, respectively. The noise is mitigated and the distribution function matches the physical expectation of a continuous representation. The peak at approximately 5 eV relates to the surface binding energy and is found to be more pronounced than in the displayed reference data.
The Al EDFs integrated over the energy axis yield the Γ out Al /Γ in Ar + flux ratios. Prediction and reference data are presented in Figure 6 (a). Fluctuations are reduced significantly but persist weakly. Above a threshold of approximately E in Ar + ≈ 50 eV an expected monotonously increasing dependence is observed. The threshold is commonly referred to as sputtering threshold and also relates to the surface binding energy. A raise of the Γ out Al /Γ in Ar + flux ratio with increasing incident Γ in Al /Γ in Ar + flux ratios is predicted by the PSNN, independent of the Ar + ion energy. This is related to the reflection of Al atoms, which in the hybrid RMD/tfMC case study are found to be Ar + ion energy independent too. The hybrid RMD/tfMC and PSNN simulations predict the mean reflection ratio for incident Al neutrals to equal 9.97 % and 9.57 %, respectively (assuming E in Ar + = 3 eV to rule out sputtering of Al surface atoms and separate the reflection of incident Al atoms). The reflection and sputtering of Al atoms are predicted accurately [39].
In the case study that provided the reference data for training the machine learning model, an ongoing cycle of Ar incorporation, clustering, and outgassing was revealed [39]. This cycle was found to cause correlated fluctuations for the Γ out Ar /Γ in Ar + flux ratio and other surface properties (e.g., Ar concentration x Ar ). The reference and predicted Γ out Ar /Γ in Ar + flux ratio are shown in Figure 6  on the incident Γ in Al /Γ in Ar + flux ratio is observed. Higher ion energies cause a deeper and more likely Ar implantation. The majority of incorporated Ar atoms may cluster but do not yet outgas. After the 2nd dose, a saturation of the surface slab and outgassing events lead to the predicted increase of Γ out Ar /Γ in Ar + . Its dependence on the incident Γ in Al /Γ in Ar + flux ratio, which is observed now, is reasoned twofold: i) A more pure Ar + ion flux (i.e., Γ in Al /Γ in Ar + → 0) consists of a higher number of impinging Ar + ions leading to faster saturation of the surface slab; ii) An eased (hindered) outgassing for a capping surface layer that becomes thinner (thicker) over time due to a number of sputtered Al atoms that exceeds (falls below) the number of built-in Al atoms per bombarding Ar + ion, Γ in Al /Γ in Ar + < Γ out Al /Γ in Ar + Γ in Al /Γ in Ar + > Γ out Al /Γ in Ar + . This may be identified from Figure 6 (a) and Figure 6 (b) [39].
The ion bombardment induced defect structure is quantified by means of the ring statistical connectivity profile (see Section II) [42][43][44]. Out of which R C,r was found to be the most relevant property and is shown in Figure 7 for r ∈ [3,6]. Rings with three or four members are sufficient to describe the ideal face-centered-cubic (fcc) Al lattice structure.
Rings with five nodes are related to interstitials and surface reconstructions, whereas rings with six members relate to vacancies and Ar atoms occupying vacant Al sites, respectively [39]. Hence, R C,3 and R C,4 (R C,5 and R C,6 ) are found to monotonously decrease (increase) for higher Ar + ion energies. An effect of the incident flux ratio Γ in Al /Γ in Ar + is predicted for Ar + ion energy E in Ar + greater than approximately 100 eV. This threshold was observed also in the underlying RMD case study and is due to the minimal kinetic energy required to cause the ongoing cycle of Ar incorporation, clustering, and outgassing [39]. Persistent fluctuations are attributed to the varying states within this cycle after the 1st and 2nd dose of impinging particles. Hence, the herein considered properties are intrinsically correlated, which is thus assumed to be accepted by the PSNN. The remaining properties of the connectivity profile (i.e., P N,r , P min,r , P max,r [42][43][44]) are in excellent agreement with the reference data (see Table II and Section IV A).
The predicted defect structure is used to further predict the corresponding surface state.
The mass density (not shown) and Ar concentration x Ar (shown in Figure 8 (a)) entail the same threshold at approximately 100 eV for the Ar cycle. Their dependence on the incident Γ in Al /Γ in Ar + flux ratio for ion energies greater than approximately 100 eV is similarly observed. A more pure Ar + ion flux (i.e., Γ in Al /Γ in Ar + → 0) leads to a higher number of impinging Ar + ions, saturating the surface slab more effectively [39]. The predicted biaxial and shear stresses do not hold such a dependency on the incident Γ in Al /Γ in Ar + flux ratio for any herein considered Ar + ion energy (consistent with [39]). Notably, fluctuations are mitigated for all predicted surface state properties (i.e., ρ , x Ar , σ xx +σ yy 2 , τ xy ) which after the 1st dose are fed as input to the PSNN for the 2nd dose of impinging particles.
The machine learning model is found to reliably describe the Ar concentration x Ar . The Ar clustering is assessed next by comparing the reference with the predicted Ar gas porosity φ Ar in Figure 8   leads to an increased volumetric partitioning due to the additional Ar + ions bombarding the surface, forming larger clusters [39]. The predictions' agreement with the reference data may be assessed more quantitatively by considering Table II in Section IV A. This work aims further advance the development of a data-driven plasma surface interaction model [27]. It is proposed that a common description of the target and the substrate plasma-surface interactions may be attained, which builds on data obtained from hybrid RMD/tfMC simulations. This description entails a comprehensive system state (i.e., surface state S s , defect structure D). It is established through a physics-separating artificial neural network that conceptually allows to consistently combine mixed input data with different descriptive power. In the future, a combination of simulation data (e.g., S s , D) with experimental data (e.g., only S s ) may be adopted. The PSNN is optimized by applying an anisotropic self-adaptive evolution strategy with intermediate recombination.
The PSNN is found to accurately predict the surface state S s and its properties A s taking only the ion bombardment induced defect structure D into account. Moreover, it describes accurately the first stages of a transient surface state evolution (i.e., S s → S s after the 1st dose → S s after the 2nd dose). This indicates altogether that the assumption (hypothesis) of a complete system state description by either the defect structure D or surface state S s to be justified [46].
The proposed machine learning model and procedure effectively adopts the physical fluctuations for an ongoing cycle of Ar incorporation, clustering, and outgassing. At the same time it reduces the statistical variations of the atomic system that previously hindered the interpretation of the hybrid RMD/tfMC case study. Observations, physical principles, and conclusions revealed in the referenced case study are found to be preserved by the PSNN [39]. Hence, the effect of the plasma (i.e., P = {E Ar + , Γ in Al /Γ in Ar + }) on the relevant surface properties (i.e., D, S s , A s ) can be revealed more clearly.
The PSNN has learned to generalize on a limited data set consisting of only 200 samples. The proposed constrained mixup augmentation is assumed to address this issue by establishing a more robust basis for training the model. Such an effect has been shown for other methods that virtually (artifically) extend the database [47][48][49]. It can be perceived as a hypothesis of linear superposition that is provided to the machine learning model. This hypothesis is accepted in this work for training set sample pairs along the Ar + ion energy, the Al neutral to Ar + ion flux ratio, and the dose axis with an indices radius equal less than 3.4.
In summary, the impingement of up to 15280 particles was simulated to generate a small data set that consists of only 200 samples [39]. The system trajectories are found to be predicted accurately for a few ( 2) steps (particle doses). This is enabled by the external recurrent connection and the corresponding feedback of the surface state S s . The PSNN model representation is assumed to not allow for a continued evolution ( 2 particle doses).
This issue is related to the herein considered four dimensional surface state and the 16 dimensional defect structure (probed just 200 times). The transient system state requires extrapolation, which may only be remedied by increasing the number of impinging particles with immense computational costs. This aspect is even more challenging when considering reactive systems (e.g., the sputter deposition of AlN thin films) due to the eventually required application of variable charge models (e.g., QTE, ACKS2) in the RMD or hybrid RMD/tfMC simulations [61][62][63].

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.