Model discovery for studies of surface morphological modifications based on Kuramoto-Sivashinsky dynamics

A wide range of observations in studies of surfaces exposed to ion beams can be explained and analyzed successfully by continuum models of the Kuramoto-Sivashinsky type. Despite certain progress in the theoretical understanding of the model parameters on the basis of atomistic models, much of the applications are based on phenomenological determination of several unknown quantities. In this work a numerical tool is discussed and investigated, which allows us to determine model coefficients and complex model structures from experimental findings. The method resembles known approaches in machine learning and data-driven reconstruction techniques. To keep the discussion on a fundamental level, numerical simulations are conducted by employing a scaled test model. The reconstruction technique is demonstrated for this model system and shows a high accuracy in recovering input parameters for situations without beam noise. As an application to an unknown system to be explored, the algorithm is then applied to a system with lognormal distributed ion bombardment. The impact of the beam fluctuations in the proposed model are discussed. Perspectives of the numerical algorithm for an analysis of experimental data are addressed.

Figure 1

Panel (a) shows the temporal evolution of the layer thickness $W$ for the cases with constant beam flux characterized by $\langle J/J_0\rangle =0.5$, 1.0, 1.5 and $\langle {\tilde{J}}^2\rangle =0$. The reference time is defined by ${\tau }_0=K_0/{\nu }_0^2$. The saturation level $\langle W\rangle \approx 4{\nu }_0/{\lambda }_0$ is almost the same for all cases. The exponential growth for small times $t$ corresponds to the beam flux dependence of the linear growth rate ${\gamma }_{*}$ of Eq. (6). To illustrate this panel (b) shows a semilogarithmic plot of the same time traces as in the top figure with fit functions (giving the straight lines) included. According to linear theory the fit functions were chosen as $f\left(t\right)=bexp\left({\gamma }_{*}t\right)$, where ${\gamma }_{*}={\langle J/J_0\rangle }^2/\left(4{\tau }_0\right)$. Fit parameters found are$b=9.6\times {10}^{-4}, b=1.5\times {10}^{-3}$, and $b=2.0\times {10}^{-3}$ for $\langle J/J_0\rangle =0.5$, 1.0, 1.5, respectively. Therefore, the straight lines in the semilogarithmic plot become steeper as $\langle J/J_0\rangle$ increases.

Figure 1

Panel (a) shows the temporal evolution of the layer thickness $W$ for the cases with constant beam flux characterized by $\langle J/J_0\rangle =0.5$, 1.0, 1.5 and $\langle {\tilde{J}}^2\rangle =0$. The reference time is defined by ${\tau }_0=K_0/{\nu }_0^2$. The saturation level $\langle W\rangle \approx 4{\nu }_0/{\lambda }_0$ is almost the same for all cases. The exponential growth for small times $t$ corresponds to the beam flux dependence of the linear growth rate ${\gamma }_{*}$ of Eq. (6). To illustrate this panel (b) shows a semilogarithmic plot of the same time traces as in the top figure with fit functions (giving the straight lines) included. According to linear theory the fit functions were chosen as $f\left(t\right)=bexp\left({\gamma }_{*}t\right)$, where ${\gamma }_{*}={\langle J/J_0\rangle }^2/\left(4{\tau }_0\right)$. Fit parameters found are$b=9.6\times {10}^{-4}, b=1.5\times {10}^{-3}$, and $b=2.0\times {10}^{-3}$ for $\langle J/J_0\rangle =0.5$, 1.0, 1.5, respectively. Therefore, the straight lines in the semilogarithmic plot become steeper as $\langle J/J_0\rangle$ increases.

Figure 2

Distribution functions $\text{PDF}\left(h\right)$ for the cases with constant beam flux characterized by $\langle J/J_0\rangle =0.5$, 1.0, 1.5 and $\langle {\tilde{J}}^2\rangle =0$ (dashed lines). Gaussian fits are also shown in the figure for all three cases based on the first and second moments computed for the time traces of $h$ (solid lines). The statistical analysis was performed for the saturated and statistically stationary state.

Figure 3

Snapshots of the morphological structure for the cases with constant beam flux characterized by (a) $\langle J/J_0\rangle =0.5$, (b) 1.0, (c) 1.5 and $\langle {\tilde{J}}^2\rangle =0$. The domain size is chosen as $L=100K_0^{1/2}/{\nu }_0^{1/2}$. The change in the size of the patterns follows the scaling $l_{*}\sim {\nu }^{-1/2}\sim {\langle J/J_0\rangle }^{-1/2}$.

Figure 4

Shown are the squared components $|c_{ij}{|}^2$ of the time-averaged squared gradient $\overline{{\left|\mathbf{\nabla }h\right|}^2}$ in $m-n$ wave-number plane for the cases with constant beam flux characterized by (a) $\langle J/J_0\rangle =0.5$, (b) 1.0, (c) 1.5 and $\langle {\tilde{J}}^2\rangle =0$. The Fourier decomposition is given by Eq. (23). The change in structural size displayed in Fig. 3 is reflected by circular patterns of different radii.

Figure 5

Temporal evolution of the squared gradient $\overline{{\left|\mathbf{\nabla }h\right|}^2}$ for the cases with constant beam flux characterized by $\langle J/J_0\rangle =0.5$, 1.0, 1.5 and $\langle {\tilde{J}}^2\rangle =0$. The different saturation levels are well separated at $\overline{{\left|\mathbf{\nabla }h\right|}^2}=$ 2.5, 4.9 and $7.2{\nu }_0^3{K}_0^{-1}{\lambda }_0^{-2}$ for $\langle J/J_0\rangle =0.5$, 1.0, 1.5, respectively.

Figure 6

Shown are the 199 different results for the reconstruction of the coefficient $\lambda$ obtained by evaluation of subsequent pairs of 200 snapshots of the simulations all of them with $\langle J/J_0\rangle =1$ but with different increments ${\tau }_{\mathrm{\Delta }}/\mathrm{\Delta }t=1$, 10, 100 between the two respective snapshots used for the reconstruction method. The calculated values are represented by their relative deviation from the known input parameter ${\lambda }_0$ in percent.

Figure 7

Shown are the distribution functions (a) $\text{PDF}(\nu /{\nu }_0)$, (b) $\text{PDF}(\lambda /{\lambda }_0)$, and (c) $\text{PDF}(K/K_0)$ for the reconstruction using different increments ${\tau }_{\mathrm{\Delta }}/\mathrm{\Delta }t=1$, 10, 100 and based on 200 snapshots of the simulation with $\langle J/J_0\rangle =1$. Gaussian fits $G(\nu /{\nu }_0), G(\lambda /{\lambda }_0)$, and $G(K/K_0)$ are also shown in the figure based on moments of $\nu , \lambda$, and $K$ computed for the series of reconstruction results. The narrowest distribution occurs in all cases (a), (b) and (c) for the parameter ${\tau }_{\mathrm{\Delta }}/\mathrm{\Delta }t=10$, the widest for ${\tau }_{\mathrm{\Delta }}/\mathrm{\Delta }t=1$.

Figure 8

Plots of the relative error of the averaged values (a) ${\mu }_{\nu }$, (b) ${\mu }_{\lambda }$, and (c) ${\mu }_K$ obtained by averaging the reconstructed coefficients over all 199 reconstructions. For all cases $\langle J/J_0\rangle =0.5$, 1.0, 1.5 the relative error can be less than 1% as long as the incremental time ${\tau }_{\mathrm{\Delta }}$ is smaller than $5\times {10}^{-3}{\tau }_{*}$.

Figure 9

Plots of the relative standard deviation (a) ${\sigma }_{\nu }$, (b) ${\sigma }_{\lambda }$, and (c) ${\sigma }_K$ obtained by averaging the reconstructed coefficients over all 199 reconstructions. For all cases $\langle J/J_0\rangle =0.5$, 1.0, 1.5 the curves are almost identical and a pronounced minimum occurs at around ${\tau }_{\mathrm{\Delta }}\sim 5\times {10}^{-3}{\tau }_{*}$.

Figure 10

Temporal evolution of (a) the layer thickness $W$ and (b) the squared gradient $\overline{{\left|\mathbf{\nabla }h\right|}^2}$ for the cases with a fluctuating beam flux characterized by $\langle J/J_0\rangle =1.0, \langle {\tilde{J}}^2/J_0^2\rangle =0.09$, and $\tau /{\tau }_0={10}^{-4}, {10}^{-3}, {10}^{-2}, {10}^{-1}, {10}^0, {10}^1$. For all values of $\tau /{\tau }_0$ in the statistically stationary phase the temporal average of the layer thickness is $\langle W\rangle \approx 4{\nu }_0/{\lambda }_0$. The respective mean values of the squared gradient $\langle \overline{{\left|\nabla h\right|}^2}\rangle$ lie in a range between 4.3 and $4.9{\nu }_0^3{K}_0^{-1}{\lambda }_0^{-2}$. The scattering of $\overline{{\left|\mathbf{\nabla }h\right|}^2}$ decreases when the scaled correlation time $\tau /{\tau }_0$ is reduced.

Figure 11

Reconstruction results for the cases with fluctuating beam flux characterized by $\langle J/J_0\rangle =1.0, \langle {\tilde{J}}^2/J_0^2\rangle =0.09$, and $\tau /{\tau }_0={10}^{-4}, {10}^{-3}, {10}^{-2}, {10}^{-1}, {10}^0, {10}^1$. Shown are the distribution functions (a) $\text{PDF}(\nu /{\nu }_0)$ and (b) $\text{PDF}(\lambda /{\lambda }_0)$ for an increment ${\tau }_{\mathrm{\Delta }}/\mathrm{\Delta }t=100$ and different correlation lengths $\tau /\mathrm{\Delta }t={10}^{-1}, {10}^0, {10}^1, {10}^2, {10}^3, {10}^4$. The dotted lines represent the raw numerical data and the solid line the lognormal fits $H(\nu /{\nu }_0)$ and $H(\lambda /{\lambda }_0)$. The evaluation is done for 2000 snapshots. The lognormal distribution of the beam flux is also indicated by the thick red curve.

Figure 12

Plots of the relative error of the averaged values (a) ${\mu }_{\nu }$ and (b) ${\mu }_{\lambda }$ obtained by averaging the reconstructed coefficients over 1999 reconstructions. Fluctuating beams have been considered with $\langle J/J_0\rangle =1.0, \langle {\tilde{J}}^2/J_0^2\rangle =0.09$, and $\tau /{\tau }_0={10}^{-4}, {10}^{-3}, {10}^{-2}, {10}^{-1}, {10}^0, {10}^1$.

Figure 13

Plots of the relative standard deviation (a) ${\sigma }_{\nu }$ and (b) ${\sigma }_{\lambda }$ obtained by averaging the reconstructed coefficients over 1999 reconstructions. Fluctuating beams have been considered with $\langle J/J_0\rangle =1.0, \langle {\tilde{J}}^2/J_0^2\rangle =0.09$ and $\tau /{\tau }_0={10}^{-4}, {10}^{-3}, {10}^{-2}, {10}^{-1}, {10}^0, {10}^1$. The reference value denoted by “ref” in the plots marks the standard deviation of the imposed beam fluctuations ${\langle {\tilde{J}}^2/J_0^2\rangle }^{1/2}=0.3$, which represents the limit of long correlation times $\tau /{\tau }_0\gg 1$, i.e., for cases when the surface dynamics is fast enough to adapt the slow beam fluctuation dynamics.