Rémi Dussart
email: remi.dussart@univ-orleans.fr

Thomas Tillocher

Gaelle Antoun

Jack Nos

Philippe Lefaucheux

Aurélie Girard

Christophe Cardinaud

Plasma cryogenic etching processes: what are the mechanisms involved at very low temperature?

In 1988, the first cryoetching process of silicon was introduced by a Japanese team to form anisotropic profiles [1]. The substrate was cooled down to a very low temperature (of the order of -100°C) in order to freeze chemical reactions on vertical sidewalls while etching the structure bottom submitted to the ion bombardment. Some years later, it was shown that a passivation mechanism was actually responsible for the etching anisotropy [2]. Standard deep cryoetching of silicon consists in using a plasma mixture of SF 6 /O 2 to etch silicon and form a SiO x F y inhibiting film at the sidewalls. Many cryogenic etching experiments have been reported showing high aspect ratio structures and interesting behaviors by varying different parameters, but without completely explaining the physical and chemical mechanisms behind this particular process at low temperature [3][4][5]. The role of physisorption of SiF 4 in the cryogenic process was highlighted by Tinck et al. who performed molecular dynamic modelization at low temperature [6]. The SiF 4 physisorbtion was also investigated experimentally to better understand the involved mechanisms in the formation of the passivation layer [7]. A plasma of SiF 4 /O 2 regularly injected during the process can help in enforcing the passivation layer: this is the so-called StiGer process [8]. An example of deep etching of submicron trenches in silicon is given in Fig. 1.

Fig. 1: High aspect ratio trenches obtained using the STiGer process Decreasing the wafer temperature can also be useful to reduce the diffusion of reactive species at the surface. For example, when etching porous low K materials, the diffusion of reactive species can damage the material and modify the K value. Lowering the wafer temperature can help in maintaining the quality of the material as it was shown in [9,10]. More recently, a cryogenic process was used for atomic layer etching (ALE) of SiO 2 . Instead of initiating a plasma of C 4 F 8 at the beginning of each cycle, the gas was simply injected and C 4 F 8 molecules could physisorb at the SiO 2 surface. This step was followed by a plasma of argon to induce chemical reactions at the surface. An example is shown in Fig 2 where the evolution of the thickness of a SiO 2 layer measured by in-situ ellipsometry is shown versus time during 8 cycles for two different temperatures of -110°C and -120°C. At -110°C, no etching is observed whereas an etch per cycle of 0.3 nm is obtained at -120°C. Mass spectrometry analysis was also carried out to study the adsorption and desorption of different species during the process. The control of physisorbed species is crucial in this particular process: If the surface residence time of physisorbed molecules is too small, no etching occurs during the argon plasma step. One of the advantages of using this cryogenic process for ALE relies on the low contamination of the reactor walls, avoiding process drift. Finally, cryo ALE experiments have been investigated using SiF 4 /O 2 plasma as a modification step, before the Ar plasma step. By reducing the temperature, the regime can be tuned from deposition to etching. In-situ XPS analyses were performed at IMN at Nantes. It was shown that the stoichiometry oft he deposited layer is significantly modified when the substrate is cooled. This property could be used to selectively etch Si 3 N 4 layers while a deposition regime was obtained on Silicon.

All these different aspects of cryogenic processes will be presented more in detail to better understand the involved physical and chemical mechanisms.

Fig. 2 :

 2 Fig. 2 : Thickness evolution of the SiO 2 layer versus time during 8 cycles for two different temperatures.

Acknowlegment: This work was financially and scientifically supported by Tokyo Electron Ltd. It is also financed by ANR agency (contract 20-CE24-0014-02 -PSICryo). The authors would like to acknowledge CERTeM 5.0 platform, which provided most of the equipments used for this work and the CNRS-Réseau Plasmas Froids for the Optimist platform.


domain: info math

text:
Hubert Baty
email: hubert.baty@unistra.fr

Solving stiff ordinary differential equations using physics informed neural networks (PINNs): simple recipes to improve training of vanilla-PINNs

Physics informed neural networks (PINNs) are nowadays used as efficient machine learning methods for solving differential equations. However, vanilla-PINNs fail to learn complex problems as ones involving stiff ordinary differential equations (ODEs). This is the case of some initial value problems (IVPs) when the amount of training data is too small and/or the integration interval (for the variable like the time) is too large. We propose very simple recipes to improve the training process in cases where only prior knowledge at initial time of training data is known for IVPs.

For example, more physics can be easily embedded in the loss function in problems for which the total energy is conserved. A better definition of the training data loss taking into account all the initial conditions can be done. In a progressive learning approach, it is also possible to use a growing time interval with a moving grid (of collocation points) where the differential equation residual is minimized. These improvements are also shown to be efficient in PINNs modelling for solving boundary value problems (BVPs) as for the high Reynolds steady-state solution of advection-diffusion equation.

I. INTRODUCTION

The use of neural networks (NNs) to solve differential equations is revisited in a tutorial paper (see Baty & Baty 2023 and references therein). Basically, a classical NN can give a non linear approximation of the whole desired solution by using a dataset of known particular values. This is a supervised learning method which consists in finding a mapping function between given inputs values (for example the time variable) and their output values (for example the solution). This dataset is used to parameterize the NN such that it minimizes the error between solution predicted by the NN and true known solution from the dataset during a training procedure. The convergence is achieved by minimizing a loss function which expression is based on error estimate, using for example the mean squared error. Finding "good" parameters is achieved by solving an optimization problem using a gradient descent algorithm that relies on automatic differentiation to back-propagate gradients through the network (Baydin et al. 2018).

However, the amount of available known solution data (so called training data) is in general very small. These are typically the initial or boundary values in problems involving ordinary differential equations (ODEs). In these cases, NN cannot be used as it is a bad extrapolation tool. An approach called physics-informed neural networks (PINNs), has been thus proposed in order to tackle the limitations of classical NN (Raissi et al. 2017(Raissi et al. , 2019)). The basic idea is to provide additional information corresponding to the physics. The method consists in evaluating the solution at some other set of data points (called collocation points) at which the equation residual is minimized. A second loss function corresponding to the physics is thus defined and added to the previous one in the learning process. The training is penalized by this additional constraint and the space of available solutions is thus restricted, being partly driven by the original data and also partly driven by the physics.

The use of PINNs to solve ODEs has been clearly illustrated in Paper 1 (Baty & Baty 2023). In particular, it has been shown that data knowledge representing only the initial conditions can be sufficient when the equations are weakly non linear. However, this not the case for strongly non linear stiff problems (as for Van Der Pol oscillator for example). As a consequence, the use of PINNs in the latter cases can be successful under the condition that a larger set of training data is used. The training procedure can be improved by adding another physical information like the energy conservation (i.e. an additional constraint) in problems for which it is effectively conserved, but such limitations called failure modes in the literature remain (Krishnapriyan et al. 2021). Note that we focus on ODEs in this study for the sake of simplification, but the problems discussed in this paper also concern the integration of partial differential equations (PDE's). These difficulties are even worse when one to integrate over a rather large time interval. Many improvements of these socalled vanilla-PINNs have been proposed in the literature, that are often based on self adaptive procedures (of the collocation points, activation function, etc.). However, none of these methods appears to be effective on all the equations of concern (Karniadakis et al. 2021, Xiang et al. 2022). In this work, we propose very simple recipes that can be easily implemented with the vanilla-PINNs in order to ameliorate the training procedure.

The improved results obtained in this study are illustrated by directly comparing to some benchmark results shown in Paper 1.

The paper is organized as follows. In Section 2, we give a short summary of the PINNs technique. Two recipes based on modifying the loss function in two ways are investigated in Section 3. In Section 4, we investigate a third recipe where the grid of collocation points is modified along with the progression of the training process. Finally, conclusions are drawn in Section 5.

II. PHYSICS-INFORMED NEURAL NETWORKS

A. The basics of neural networks for ODEs

We first consider the desired solution u(t) of an ODE (see below) with u θ u being the approximated solution at different t values, where θ is a set of model parameters. For IVPs, the variable t would generally be a time parameter and should be replaced by a spatial coordinate parameter for boundary value problems (BVPs). Using a classical neural network approximating the desired solution, we can write,

u θ (t) = (N L • N L-1 ... N 0 )(t), (1)
where the operator • denotes the composition and θ = {W l , b l } l=1,L represents the trainable parameters (with weight matrices and bias vectors) of the network (see Paper 1 for more details on the N l functions). The network architecture schematized in Figure 1, is organized in L + 1 layers with neurons connected in adjacent layers. A single input layer containing the input variables t is connected to L -1 hidden layers (two layers with four neurons in the schematic example of Figure 1), and finally to an output layer for the solution u θ . The goal is to calibrate its parameters θ such that u θ approximates the target solution u(t). An activation function is also necessary in order to to introduce non-linearity into the output of each neuron. In this work, the most commonly used hyperbolic tangent tanh function is chosen.

The optimization problem is based on the minimization of a loss fonction that can be expressed as,

L data (θ) = 1 N data N data i=1 u θ (t i ) -u i data 2 , (2)
where a set of N data data is assumed to be available for the known solution at different times t i that are called the training data (i = 1, N data ), which includes the initial and/or boundary conditions.

B. The basics of vanilla-PINNs for ODEs

Let us now introduce an ODE in the following residual form with imposed initial conditions and/or boundary conditions depending on the problem considered. The exact number of initial/boundary conditions necessary to solve the equation obviously depends on the order of the equation. Note that in case of an ODE with an order n higher or equal to two, an equivalent system of n equations can be also used (see below).

F t, u, du dt , d 2 u dt 2 , ... = 0, t ∈ [t 0 , T ] , (3)
In the original form, the basics of vanilla-PINNs is based on the use of a second loss function defined as

L F (θ) = 1 N c Nc i=1 |F[u θ (t i )]| 2 , (4)
that must be evaluated on a set of N c points generally called collocation points that are not necessarily coinciding with the training data points. Note that one can evaluate exactly the differential operators at the collocation points in L F and F by using automatic differentiation. This same technique of automatic differentiation is also used to compute derivatives with respect to the network weights (i.e. θ), that is necessary to implement the optimization procedure (see below). Note that contrary to the use of standard numerical schemes, the derivatives can be obtained at machine precision. In this work, we use Pytorch Python open source software libraries facilitating thus the latter operations.

A composite total loss function can be consequently formed as

L(θ) = ω data L data (θ) + ω F L F (θ), (5)
where an optimal choice of values for hyper-parameters (ω data , ω F ) allow to ameliorate the eventual unbalance between the partial losses during the training process. These weights can be user-specified or automatically tuned. In the present work, for simplicity we fix the ω data value to be constant and equal to unity, and the other weight parameters including ω F are determined with values varying from case to case. A gradient descent algorithm is used until convergence towards the minimum is obtained for a predefined accuracy (or a given maximum iteration number) as for the residual equation and loss function L F . Consequently, one gets

θ k+1 = θ k -η∇ θ L(θ k ), (6)
L E (θ) = 1 N c Nc i=1 |E(t i ) -E 0 | 2 , (7)
The total loss function is consequently modified by incorporating a new term weighted by a hyper-parameter ω E ,

L(θ) = ω data L data (θ) + ω F L F (θ) + ω E L E (θ). ( 8
)
It has been shown that the results are considerably ameliorated when compared to a case without this additional constraint. For example, this is illustrated in Figure 8 of Paper 1

for the harmonic oscillator problem.

B. Using an hybrid loss function for data

Nevertheless, for the following anharmonic equation,

d 2 u dt 2 + ω 2 0 u 3 = 0, (9)
two training data were necessary (see Figure 13 in Paper 1 and left panel of Figure 3 in present paper). This is not a complete surprise as two conditions are necessary to integrate such second order equation using analytic or classical numeric methods. Another option would be to solve an equivalent system of two first order differential equations, as done for

N eval i=1 u θ (t i ) -u i eval 2
, where the evaluation u θ (t i ) is done on N eval = 1000 points uniformly distributed within the whole time interval, and where u i eval is the expected exact solution at t = t i .

the non linear pendulum (Figure 12 in Paper 1). Indeed, in the latter case the NN is learning the solution and also the first order derivative on the whole time interval.

In this study, we propose another strategy. Indeed, we can choose the value of the first order time derivative evaluated at the first collocation point, in order to constrain it to converge towards its true initial value du dt (t = 0). Thus, the previous loss fonction L data that contains only the initial data value on y(t = 0) = y 0 can be modified to include a second term,

L data = |u θ (t = 0) -u 0 | 2 + ω d |u θ (t = 0) -u 0 | 2 , (10)
where u 0 = du dt (t = 0) and ω d is a new weight parameter not necessarily equal to one. This is an hybrid loss fonction as the first term involves the training data and the second term the first collocation point (also used to evaluate the loss fonctions on equation residual and total energy residual). This is illustrated in Figure 3 for the anharmonic oscillator with initial conditions u 0 = 1.5 and u 0 = 0, using also ω 0 = 15.5 with an integration over the time interval t ∈ [0, 1]. Indeed, we have obtained two PINN solutions choosing the following hyper parameters, η = 1.1 × 10 -3 , ω data = 1, ω F = 1 × 10 -5 , and ω E = 1 × 10 -6 for the two cases. The NN architecture is made of 5 hidden layers with 32 neurons per layer.

As explained above, in left panel of Figure 3 is plotted the converged PINN solution for the method without the improvement using two training data points (i.e. at t = 0, and t = t 1 > 0) and N c = 44 collocation points. In right panel, one can see the second PINN solution for the method using the improvement (i.e. with the hybrid data loss function) with only one training data point and N c = 32 collocation points. In the latter case, the chosen additional weight is ω d = 1 × 10 -3 . The corresponding histories of the total loss and mean squared error (M SE) that are plotted in Figure 4, show similar convergence during the training processes. Note that, this is important to also examine the M SE, as it is a direct measure of the error contrary to the total loss that is a composite function in PINNs.

We can conclude that this first recipe using the hybrid data loss not only allows to reduce the training data set to the sole initial condition on u (that is u 0 , as the initial derivative is imposed using the first collocation point), but it also allows to reduce the minimum number of collocation points.

IV. APPLYING A MOVING GRID WITH A GROWING COLLOCATION IN-TERVAL

In Paper 1, it has been shown that when the differential equation is particularly stiff, the use of vanilla-PINNs requires a minimum number of training data that is significantly higher than unity, and which is distributed over the time integration interval. This is indeed the case for Van Der Pol oscillator.

A. Initial value problem -Van Der Pol oscillator

The Van Der Pol oscillator equation is

d 2 u dt 2 + ω 2 0 u -ω 0 (1 -u 2 ) du dt = 0, ( 11
)
where ω 0 is a normalized angular velocity, and t ∈ [0, T ]. Finally is a parameter having a value which determines the amplitude of a limit cycle in the phase space, and consequently determines the stiffness of the equation. PINN solutions obtained in Paper 1 for = 0.33, 1, and 5, show that 3, 7, and 16 training data points respectively were necessary for convergence of the training process. One must note that, in this case the total energy is not conserved and therefore it is not possible to use the loss energy constraint.

We have found that vanilla-PINN, even when using the previous improvement strategy (with hybrid data loss), fails to train when the integration interval covers a few periods of the oscillator. This is however not the case when the integration is done on a restricted time interval that is typically smaller than one period. A strategy has been proposed to overcome PINNs failures in long time for PDE's, consisting in learning progressively sequence-to-sequence until the entire space-time solution is obtained (Krishnapriyan et al. 2021). In this strategy, the data set of collocation points is progressively increased in order to finally invade the whole integration domain. New points are added along with the progression of the training process, but they remain at fixed position in time once created.

Hence, a considerably high number of points is necessary at the end of the training process.

We propose a variant of the latter strategy which requires a moderate and fixed number with 10000 time steps uniformly distributed over the integration interval.

For completeness, we have also investigated a variant of the previous PINN modelling in which, instead of using the hybrid data loss function to take into account the initial time derivative condition, we solve an equivalent system of two equations that is,

     du 2 dt - ω 0 u 1 = 0, du 1 dt -ω 0 (1 - u 2 1 3 )u 1 + ω 0 u 2 = 0 (12)
The new expected variables are defined as u 1 = u, and u 2 = u -u 3 3 -u ω 0 (u being the time derivative of u). Contrary to another possible choice where u 2 = u , this choice of variables called Lienard's transform is known to lead to values of u 1 and u 2 that have similar magnitudes (see below). Note that the initial condition for u (t = 0) = 0 is translated now into u 2 (t = 0) = u 0 -u 3 0 /3, that is u 2 (t = 0) = 2/3 for the case investigated. The result of the corresponding PINN integration using our moving collocation grid procedure is illustrated in Figure 7. We have chosen ω 0 = 15 and integrated for t ∈ [0, 1.5]. We also take the initial condition u 0 = 1 and zero initial derivative. The NN has 4 hidden layers with 42 neurons per layer. The number of training data and collocation points are N data = 2 (one per solution variable) and N c = 90 respectively. The chosen weights are ω data = 1, and ω F = 9 × 10 -3 .

The learning rate is η = 7 × 10 -4 . Now, we do not need the hybrid data loss function, as the initial first order time derivative is imposed via the training data set for u 2 .

However, one must note that our simple recipies have their own limitations. Indeed, when the stiffness of the system increases too much, one is forced to find another strategy or an enlarged set of training data.

B. Boundary value problem -steady-state solution of 1D convection diffusion equation

In partial differential equations of high Reynolds convection-diffusion problems, imposed boundary conditions involve the formation of strongly localized boundary layers. This is a stiff system associated to BVP's which is considered to be numerically challenging. In this section, we only focus on the steady-state solution in one dimension. Hence, we consider the following equation,

a du dx -ν d 2 u dx 2 -1 = 0, ( 13
)
where a is a constant speed, and ν is a viscosity parameter. The spatial integration domain is for x ∈ [0, 1], x being a normalized spatial coordinate. The expected exact solution is,

u(x) = x - e -R(1-x) -e -R 1 -e -R , (14)
when the imposed boundary conditions are u(0) = u(1) = 0, with the definition of the Reynolds number R = a/ν.

In order to illustrate our PINN solution using the previously described moving grid strategy with a growing collocation interval, we have considered a high Reynolds case for R = 50 Finally, when the Reynolds number is higher, again our simple recipe shows its limitation and more sophisticated strategy is required.

V. CONCLUSION

In this work, we have presented a few simple recipes with the aim to improve drawbacks of the vanilla-PINNs for solving differential equations. This is indeed the case of stiff ODE's, requiring thus the use of a too large set of data representing the prior knowledge of the solution. In other words, the sole knowledge of the initial/boundary conditions is not sufficient. Modifying the total loss function in two ways, via an hybrid definition of the data loss, or adding another partial loss associated to the conservation of the total energy (when it is possible) is a first possible improvement. Another interesting idea consists in using a moving collocation grid with a growing interval for the evaluation of the equation residual. Compared to classical numerical integration methods, PINNs still fail in terms of robustness because of these failures inherent to the the training process. However, this technique using neural networks is relatively recent and many other improvements are probably expected in the future years. Nevertheless, the PINN formulation offers interesting advantages over classical methods. Indeed, it is a meshless method, and once trained the solution can be quasi-instantaneously generated.

FIG. 1 .

 1 FIG. 1. Schematic representation of a structure example for a standard NN. The input layer has one input variable (i.e. one neuron noted σ 0 1 ) representing for example a time coordinate. Three hidden layers with four neurons per layer are connected with the input and the output layer (one neuron noted σ 4 1 ), where the latter has a single variable (one neuron) representing the desired solution u θ .

FIG. 2 .

 2 FIG. 2. Schematic representation of the structure for a vanilla-PINN modelling an ODE solution. The previous NN architecture (see previous figure) is used to evaluate the residual of the ODE equation (via u θ and corresponding derivatives). The two partial loss functions (L data and L F ) are used to form a total loss function with associated weights (see text) that is finally minimized.

Figure 2 .

 2 Figure 2. Note that in this schematic figure, a single input neuron representing time or space coordinate for ODEs must be replaced by two neurons for a partial differential equation having spatio-temporal (x, t) dependences. Moreover, in cases where a set of n differential equations is considered, the output neuron must be replaced by n neurons associated with the n solution variables that need to be learned.

FIG. 4 .

 4 FIG. 4. Histories of the total loss function L and M SE corresponding to the two cases of the previous figure respectively. The M SE is evaluated using the standard expression, M SE = 1 N eval

FIG. 5 .FIG. 6 .

 56 FIG. 5. PINN solution for the Van Der Pol equation (see text), obtained using the moving grid procedure with a growing collocation interval. The different snapshots from top to bottom and left to right panels corresponds to the progression of the process as a function of the training step.The time at which the collocation points are defined are visible with the green circles.

  figure. Note that the exact solution is obtained by using a classical Runge-Kutta of order 4

FIG. 7 .

 7 FIG. 7. PINN solution of the Van der Pol equation solved as a two-dimensional system (see Equation 12) for the two variables u 1 and u 2 . The moving grid with a growing collocation domain procedure is applied. The different snapshots from top to bottom and left to right panels corresponds to the progression of the process as a function of the training step. The time at which the collocation points are defined are visible with the green circles.

FIG. 8 .

 8 FIG. 8. PINN solution of the steady-state 1D convection-diffusion equation (see Equation 13). The moving grid with a growing collocation interval procedure is applied. The different snapshots from top to bottom and left to right panels corresponds to the progression of the process as a function of the training step. The space coordinate at which the collocation points are defined are visible with the green circles.

corresponding to a = 1

 1 and ν = 0.02. For this BVP, the two boundary conditions are imposed via the training data set at x = 0 and x = 1. The results are plotted in Figures 8-9. Indeed, one can clearly see that a modest number of only N c = 27 collocationpoints is sufficient to obtain a good training process. We have checked that, without such improvement recipe, more than 40 collocation points are necessary. Moreover, a classical numerical integration method also requires an even significantly larger number of points in the spatial domain for uniform grid. For example, a finite-difference scheme requires a few points in order to resolve the quasi-singular layer which thickness is of order 1/R. For the PINN integration, we have chosen 3 hidden layers with 42 neurons per layer. The number of training data and collocation points are N data = 2 and N c = 27 respectively. The chosen weights are ω data = 1, and ω F = 5 × 10 -1 . The learning rate is η = 1 × 10 -3

  In this way, the training process is made more progressive. More sophisticated methods based on self-adaptive methods could be also developed (see McClenny & Braga-Neto 2023, Karniadakis et al. 2021, and Cuomo et al. 2022 for reviews).

ACKNOWLEDGMENTS

The author thanks Léo Baty (CERMICS, ENPC) for his help in Python programming. Some of the Python codes used to make the figures are available from the Github repository at https://github.com/hubertbaty/PINNS-EDO2.