Surrogate model for optimizing annealing duration of self-assembled membrane-cavity structures

We propose a scheme to establish a surrogate model for optimizing the annealing duration of the self-assembled membrane-cavity structures from hole patterned silicon wafers. Although it has been reported that the design space of post-annealing shape can be extended by increasing the dimensions of hole arrays, the annealing duration for large hole arrays has not been well examined. A two-dimensional axisymmetric phase-field model in commercial FEM software is employed to establish the surrogate model with respect to three variables (i.e., radius, aspect ratio (AR), and normalized spacing). The established surrogate model based on the neural network indicates that the hole radius dominantly affects annealing duration and the temperature elevation (i.e., acceleration of diffusion speed) is necessary to achieve the practical annealing duration when the hole radius is larger than 1 μm.


Introduction
Self-assembly is one of the effective methods to fabricate the membrane-cavity (e.g., Silicon-on-Nothing, SON [1][2][3] or Germanium-on-Nothing, GON [4]) structures for thin [5] or flexible solar cells [6] via transferring membranes with the thickness less than 1 µ m, microchannel resonator [7] or multiscale fluidic channel [8] in robust mono-crystalline structures, and pressure sensors monitoring the deflection of suspended membranes [9,10]. During high-temperature annealing where surface diffusion is accelerated, the hole patterns on wafers are covered with adjacent materials while the covered hole becomes shorter and a sphere [11]. At a fixed annealing temperature or diffusion rate, a longer annealing duration is required to fabricate the membrane-cavity structures from larger hole patterns. However, the detailed analysis of the effect of hole dimensions (i.e., radius, height, and spacing) on the annealing duration has rarely been reported.
In this work, we propose a scheme to establish a surrogate model for optimizing the annealing duration of the self-assembly process concerning the dimensions of hole arrays. To this end, the phase-field model in commercial software (COMSOL Multiphysics v5.5) is employed considering the higher accessibility compared to in-house code. After verifying the commercial model in a full three-dimensional simulation, a two-dimensional axisymmetric model for a single hole pattern is established considering the computational cost. Finally, the sensitivity of annealing duration to hole dimensions is analyzed using the established surrogate model from the neural network as well as simulation results.

Phase-field model
The phase-field model is employed to simulate the surface diffusion of silicon and the merge of silicon or cavities during high-temperature annealing. Phase-field model in COMSOL Multiphysics is in the form of the Cahn-Hilliard equation as Eq. (1) [12] Open Access *Correspondence: bongjae.lee@kaist.ac.kr; junhcullee@kaist.ac.kr where φ and ǫ are non-dimensional phase variables and interfacial thickness, respectively. Mixing energy density ( = 3ǫσ 2 3/2 = 6.5 × 10 −8 N) is obtained from its relation to surface tension ( σ ) and ǫ . Mobility ( γ ) is obtained as γ = D s σ ν� 2 k B T 1 = 9.6 × 10 −22 m 3 ·s/kg by comparing Eq. (1) with the Mullins-Herring equation that was employed to simulate the surface diffusion process [3,13]. The properties used in the present work are listed in Table 1 while temperature (T) is 1423 K. The volumetric number of atoms and atomic volume ( ) are acquired from pure silicon in a face-centered-cubic unit cell, whereas ν is calculated by multiplying the volumetric number of atoms and ǫ . The tabulated surface tension is obtained from (001) surfaces [3]. On the other hand, ǫ is set as 45 nm considering that the smaller ǫ increases simulation accuracy while computational cost increases [12]. It is also noticeable from Eq. (1) that the properties of materials only affect the processing speed. The adopted computational resource has 32 core CPU (Intel Xeon ® gold 6130) with 198 GB RAM.
To examine the feasibility of employing the phase-field model in COMSOL Multiphysics for the simulation of self-assembly, full three-dimensional phase-field models are established using the identical hole arrays with the literature [7]. The hole patterns are in the hexagonal array and have a diameter (D) of 1 µ m, aspect ratio (the ratio of height to diameter, AR = H /D ) of 5 (i.e., H = 5 µm), and spacing (distance to the adjacent hole, S) of 0.4 µ m (see Fig. 1a). The numbers of holes on the one side of arrays (n) are 2, 4, and 6 as shown in Fig. 1a-c, respectively. As shown in Fig. 1b, the area of the simulation domain is generated considering the in-plane configuration of hole arrays (i.e., D, S, and n) with the residual width of 3 µ m for each side, while the domain height is set to be 10 µ m with the 7-µm-thick hole-patterned-silicon (i.e., 3-µ m-thick vacuum or inert gas layer on silicon). The shape evolution is simulated with a transient solver (MUMPS solver) without any modifications and obtained every 50 s for simulation time (t) of 500 s with the time step size determined by the solver. Figure 1d-f shows the selfassembled structures at 500 s for n = 2 , 4, and 6. The obtained membrane-cavity structures after shape evolution indicate the feasibility to employ the phase-field model to simulate the self-assembly process in comparison with Ref. [7], which is not shown.
On the other hand, it takes about 2 h to simulate threedimensional geometry with 510,064 elements for a hole array ( n = 2 ) that has the minimum elements among the considered hole arrays. Therefore it is practically inappropriate to employ a three-dimensional model due to its high computational cost considering that dozens or more sampled data at a longer simulation time (t) are necessary for building a surrogate model. A two-dimensional axisymmetric model of a single hole is employed to achieve a reasonable computational cost for the data sampling to build a surrogate model. Please note that the two-dimensional model can represent the effect of hole dimensions on the annealing duration, although it may not provide detailed dimensions of the self-assembled membrane-cavity structures.

Two-dimensional axisymmetric model
As the reference case, a two-dimensional axisymmetric model of a single hole is built with the radius of 0.5 µ m and AR of 3 while the silicon width beside the hole is set to be 0.35 µ m to keep a cavity from crossing the simulation domain (Fig. 2a). The properties of the two-dimensional axisymmetric model are identical to those of the three-dimensional model, with the exception of ǫ , which is reduced to 40 nm considering the minimized computational cost, and thus becomes 5.8 × 10 −8 N, while γ remains unchanged ( 9.6 × 10 −22 m 3 ·s/kg). By reducing one dimension and modeling only one hole, the computational cost of the two-dimensional axisymmetric model is dramatically reduced to 12 s for a simulation time of 200 s following the reduced number of elements to 1734.
On the other hand, the aspect ratio of the cavity is employed to determine the end of the self-assembly process for a single hole. The moment when the cavity becomes a sphere is considered as the end of shape evolution considering the volume preservation of a cavity [11]. In the current case, after 200 s of transient shape evolution, the cavity becomes a half-circle (i.e., a sphere in the three-dimensional domain), as shown in Fig. 2b. For the following parametric study, the required annealing duration ( t a ) for each case is obtained when the ratio of the cavity width to height ( r 2 /D 1 ) becomes larger than 0.49, as shown in Fig. 3. In addition, the simulation time of 200 s is designated as the reference time ( t r ) for the parametric study.

Parametric study procedure and condition
The parametric study is automated and conducted using the interconnection between COMSOL and MAT-LAB as shown in Fig. 4a. The 200 sampling points are obtained using the Latin-Hypercube (LH) sampling method [14] within the ranges in Table 2 (see Fig. 4b). Hole dimensions are represented as the radius ( r hole ), AR, and normalized spacing ( S/S min ). Especially, spacing is normalized by the minimum spacing ( S min = r sphere − r hole ) to prevent a cavity from crossing the simulation domain considering the equivalent radius of the sphere [ r sphere = ( 3 2 r 3 hole AR) 1/3 ] which have the same volume as the initial hole. As shown in Fig. 4c, geometry and mesh are built using the identical setting for each case. The transient results are obtained every 100 s for t = 30 × t r = 6000 s, which is expected to be equivalent to 5 h in the experiment considering that the experimental annealing duration was 10 min for the hole  array with the similar dimensions in Ref. [11]. After each simulation terminates at 30 × t r , the required annealing duration (i.e., simulation duration, t a ) according to hole dimensions is acquired when the cavity becomes a sphere ( r 2 /D 1 > 0.49 ) (see Fig. 4d-e). If a cavity does not become a sphere ( r 2 /D 1 < 0.49 until 30 × t r ), t a is set to be zero to mark cases where extended annealing duration is required. Again, the criterion of 30 × t r , which is expected to be 5 h in the experiment, is already too long for the practical fabrication process considering that the experimental annealing duration is generally less than 30 min in previous works on SON applications [7,9,10]. It takes about 10 min to obtain the required annealing duration for each case. Figure 5a shows the normalized annealing durations ( t a /t r ) according to r hole and AR. It is found that 30 × t r is not enough to terminate the shape evolution process or to obtain a spherical cavity when r hole is larger than 1 µ m (i.e., D > 2 µ m) regardless of the other variables. In other words, enlargement of the r hole by only two times results in more than 30 times of t a . On the other hand, the required annealing duration also increases with respect to AR, although the effect of AR on t a is less dramatic than that of r hole . Therefore, only increasing the annealing duration is not sufficient to extend the design space of membrane-cavity structures, and the acceleration of the diffusion process is necessary by elevating the annealing temperature or decreasing the pressure [4,15,16].

Parametric study and sensitivity analysis results
The obtained results that have the annealing duration of less than 30 × t r are then trained using a neural network with two layers and ten nodes. The obtained sensitivity map in Fig. 5b indicates that the annealing duration is more sensitive in the order of r hole , AR, and S/S min . The sensitivity order is related to the shape evolution process. As shown in Fig. 2b, the diffusion of material over the hole and the shrinkage of the cavity are the main processes of the shape evolution, and thus the t a is more sensitive to r hole and AR than S/S min . On the other hand, the shrinkage of the cavity occurs at two sides (i.e., the top and bottom sides of the hole) once the shape evolution occurs, while the diffusion of the material over the hole occurs in one direction (see Fig. 2b). Therefore, the effective  influence of AR on the t a seems weaker compared to r hole . It is also found that the accompanied variation of r hole and AR is dominant to determine t a . Although the present work did not consider the merge of cavities, the proposed scheme will provide possibilities to optimize the annealing duration according to hole dimensions.

Conclusion
In this work, the scheme to build a surrogate model of annealing duration is proposed using commercial FEM software and a neural network. The sensitivity analysis of annealing duration to three variables (i.e., radius, aspect ratio, and normalized spacing) of hole arrays shows that  hole radius is the dominant variable and a higher diffusion rate is necessary for the practical annealing duration when the hole radius is larger than 1 µ m. The presented scheme for establishing a surrogate model will provide a deeper explanation of the self-assembly process as well as optimized annealing duration with respect to hole dimensions. It is also expected that the surrogate model based on a full three-dimensional model will serve as the design tool to fabricate membrane-cavity structures with the desired dimensions. Moreover, the suggested surrogate modeling scheme can be employed for other annealing processes, such as removing surface defects on semiconductor wafers and thermal oxidation or recrystallization of semiconductor material.