Simple Algebraic Expressions for the Prediction and Control of High-Temperature Annealed Structures by Linear Perturbation Analysis

: The prediction and control of the transformation of void structures with high-temperature processing is a critical area in many engineering applications. In this work, focused on the void shape evolution of silicon, a novel algebraic model for the calculation of ﬁnal equilibrium structures from initial void cylindrical trenches, driven by surface diffusion, is introduced. This algebraic model provides a simple and fast way to calculate expressions to predict the ﬁnal geometrical characteristics, based on linear perturbation analysis. The obtained results are similar to most compared literature data, especially, to those in which a ﬁnal transformation is reached. Additionally, the model can be applied in any materials affected by the surface diffusion. With such a model, the calculation of void structure design points is greatly simpliﬁed not only in the semiconductors ﬁeld but in other engineering ﬁelds where surface diffusion phenomenon is studied.


Introduction
The evolution of void structures in different materials is gaining presence in different research fields during the last years. For this reason, investigating new expressions for the calculation of such structures is a crucial subject in the engineering field. Fast numerical estimations can provide the needed insight for the design step of engineering processes. In this work, the focus is put on the calculation of the void shape evolution process in silicon or SON (Silicon-On-Nothing) process, introduced and developed by Mitzushima et al. [1] and Sato et al. [2,3], which is a particular case of a surface diffusion process driven by surface energy minimization. SON (Silicon-On-Nothing) can be used to build semiconductor devices at a reduced cost in comparison to SOI (Silicon-On-Insulator) devices, while being interchangeable [4]. Besides some advantages in on/off transitions [4][5][6][7], SON devices show promising results for low-voltage devices and lower parasitic resistance/capacitance [5,6,8]. Further applications include the fabrication of SON-based solar cells, as explained and manufactured by Depauw et al. [9,10], along annealing possibilities and possible manufacturing defects associated with the used annealing gas.
Regarding the void structures, several techniques have been developed to create SON morphologies such as double He+/He+ and H+ ion co-implantation with an annealing step [11,12] or direct lateral etching [5,6,13,14]. However, in this work, the void shape evolution process is studied and analyzed for obtaining SON microstructures [1,3,15,16]. The first step of the void shape evolution process consists of the etching of trenches in a silicon wafer. After that, the wafer is annealed under a high vacuum (about 10 Torr) in a deoxidizing gas, e.g., hydrogen and high temperatures (1000-1100 • C). At those conditions, a geometrical evolution occurs due to surface diffusion driven by surface 2 of 19 energy minimization. The output consists of an empty space in silicon (ESS), which is usually a spherical form, and a silicon layer above it (SON or Silicon-On-Nothing layer), which separates the void from the outer surface. If the distance between trenches is very small in a matrix of trenches, the voids may coalesce. A schematic description for both cases is given in Figure 1. The advantage of the void shape evolution process is that it is simpler and more cost-effective than other techniques. However, the prediction and control of such transformations are difficult to perform due to technical and computational difficulties. While some simulation methods have been proposed for simulating the void shape evolution process such as the moving mesh method [17], phase field/level set methods [18][19][20], or direct solving the equations with simplifications [21][22][23][24][25], a more simple, computationally inexpensive approach for calculating final equilibrium geometries is also particularly useful for quick size estimations and design, independent of the annealing conditions. The aim of this work is to provide such a model focused on the non-coalesced case, as a starting point for more complex geometries.
The developed algebraic model is general and, while based on silicon data, should be suitable for any substance that undergoes the same surface diffusion process.

Materials and Methods
In this section, a simple novel analytical model to predict the evolution of silicon microstructures is introduced. The first necessary condition is to define the kinetic surface diffusion model. The migration of the surface atoms is provoked by surface energy minimization, which is the driving force of the surface diffusion phenomena [15,[26][27][28]. At constant volume, the surface energy would depend on the surface shape, defined by the curvature of the surface. Thus, a change of curvature is what causes the migration and the development of the empty space in silicon. In the following lines, this relation is demonstrated.

Phenomena Formulation
The void shape evolution evolves by the modification of the surface free energy. With S being the surface area and γ the surface energy density, the total free energy of a surface, G, can be calculated as: As the shape of the surface changes, the stored free energy on the surface is modified. The chemical potential, µ, is the thermodynamic variable that describes the discrete free energy change: The chemical potential is understood as a change in the associated free energy when an atom is transferred from an infinite crystal (reference system) to any arbitrary surface. When an isotropic energy surface receives an atom [21,29], the corresponding chemical potential is correlated to the curvature (in absence of external electric or stresses field) as follows: In Equation (3), K is the curvature (sum of the two main curvatures), and Ω is the atomic volume.
The kinetics driven only by surface diffusion phenomenon needs a driving force to allow the surface evolution [15,16,21,29,30]. This is theoretically achieved by the variation of the surface chemical potential as follows: Assuming an atomic Boltzmann distribution, applying the Nernst-Einstein relation, and knowing that k B is the Boltzmann constant, D S the surface diffusion coefficient, and T the temperature; the drift velocity of the particles on the surface can be calculated as [30]: The atomic current density is achieved by multiplying the atoms velocity by the number of atoms on the surface (isotropic case), X S : The divergence of the atomic current density is the number of atoms that diffuses out from the corresponding surface per unit area and time. If v n is the normal velocity of the surface, v n /Ω would be also the number of atoms leaving out the surface per unit area and unit time. Operating with this equality, the normal velocity of the surface yields: Equation (7) shows the kinetics of the process on every point of the surface depending on the curvature, which is the basis of subsequent calculations. The divergences and Laplacians are applied directly on the surface direction due to the surface nature of the curvature.

First Order Perturbation Analysis
The determination of the void shape evolution characteristic values is crucial for the understanding of the underlying transformation mechanisms. A first order perturbation analysis is performed to find such values, which are fundamental for constructing simplified transformation expressions from the initial geometry parameters.
The introduced analysis of first order perturbations is similar to that of Rayleigh [31,32] for instabilities of cylindrical geometries and jets. While Rayleigh analysis was focused on fluid mechanics, the studied phenomenon in this article is based on surface diffusion kinetics. Additionally, the analysis is inspired by other comparable perturbation analysis studies [2,25,[33][34][35][36][37]. However, in contrast to previous works, not only the critical wavelength solution is used to comprehend the geometrical transformation and build an analytical model, but the minimum wavelength able to drive a transformation also plays an important role. The general conditions of the perturbation analysis are presented below.

•
The initial geometry is an unperturbed, infinitely long void cylinder, and the surface normal points inwards, which is important for the curvature calculation and its sign.

•
The analyzed transformation is from an initial cylinder to a set of spheres. The evolution of the geometry is driven by an increasing sinusoidal perturbation of wavelength λ in the longitudinal direction. The evolution is schematically shown in Figure 2.

•
The analysis is performed on the first seconds of the transformation. This is already provides if the initial sinusoidal disturbance shrinks or grows.
diffusion kinetics. Additionally, the analysis is inspired by other comparable perturbation analysis studies [2,25,[33][34][35][36][37]. However, in contrast to previous works, not only the critical wavelength solution is used to comprehend the geometrical transformation and build an analytical model, but the minimum wavelength able to drive a transformation also plays an important role. The general conditions of the perturbation analysis are presented below.

•
The initial geometry is an unperturbed, infinitely long void cylinder, and the surface normal points inwards, which is important for the curvature calculation and its sign.

•
The analyzed transformation is from an initial cylinder to a set of spheres. The evolution of the geometry is driven by an increasing sinusoidal perturbation of wavelength λ in the longitudinal direction. The evolution is schematically shown in Figure  2.

•
The analysis is performed on the first seconds of the transformation. This is already provides if the initial sinusoidal disturbance shrinks or grows. Figure 2. Schematic representation of a void shape evolution from a cylinder of radius ρ to two void spheres of spherical radius RS.
The first step of the perturbation analysis is to find a suitable expression for the evolution of the radius, ρ, of the trench. The expression of a linear (n = 1) Fourier series even function [38][39][40] is selected due to its simplicity and suitability for the following operations. A general radius, ρ, would oscillate during the perturbation along the longitudinal axis, z, with the following function: where a0 and a1 are time-dependent functions, and a0(t) can be calculated from the mean value of the cylindrical radius. Under the mentioned sinusoidal perturbation, the average value of the radius should be similar to the initial unperturbed radius ρ0, especially considering the assumed very small initial perturbations. The average value of this perturbation of wavelength λ is: Consequently, the term a0 is then: According to the equation for the surface motion driven by surface diffusion (see Equation (7)), the movement of the surface is caused by a change in curvature. The curvature is the sum of the two main curvature values. The first main curvature is calculated from the angular direction, −1/ρ, which takes a negative value due to the orientation of Figure 2. Schematic representation of a void shape evolution from a cylinder of radius ρ to two void spheres of spherical radius R S . The first step of the perturbation analysis is to find a suitable expression for the evolution of the radius, ρ, of the trench. The expression of a linear (n = 1) Fourier series even function [38][39][40] is selected due to its simplicity and suitability for the following operations. A general radius, ρ, would oscillate during the perturbation along the longitudinal axis, z, with the following function: ρ(z, t) = (a 0 (t))/2 + a 1 (t)cos(2πz/λ), where a 0 and a 1 are time-dependent functions, and a 0 (t) can be calculated from the mean value of the cylindrical radius. Under the mentioned sinusoidal perturbation, the average value of the radius should be similar to the initial unperturbed radius ρ 0 , especially considering the assumed very small initial perturbations. The average value of this perturbation of wavelength λ is: Consequently, the term a 0 is then: According to the equation for the surface motion driven by surface diffusion (see Equation (7)), the movement of the surface is caused by a change in curvature. The curvature is the sum of the two main curvature values. The first main curvature is calculated from the angular direction, −1/ρ, which takes a negative value due to the orientation of the normal vector. The second main curvature corresponds to the axial direction which is calculated with the second derivative of the radius function due to the small slopes approximation (∂ρ/∂z << 1, [21,33]), as only the first moments of the perturbation are analyzed. The curvature, K, is then: Math. Comput. Appl. 2021, 26, 43

of 19
The sign criterion is chosen depending on the pointing direction of the corresponding normal vector. Concave surfaces exhibit a negative curvature. On the other hand, convex surfaces have a positive curvature. Substituting the expression in Equation (8) into Equation (11), the curvature yields: Equation (12) must be further derived to obtain the surface diffusion equation. In order to simplify the output, an approximation of the function using the Maclaurin-Taylor series is performed [36,[41][42][43]. In this case, it is convenient to use the first order series centered on the parameter a 1 (t): Under isotropic, isobaric, and isothermal conditions, all terms of the surface diffusion expression (see Equation (7)) can be assumed to be constant, except for the curvature. The parameter B is used as the product of all these terms. After differentiating the resulting curvature terms (see Equation (14)), the normal velocity of the surface is: Due to the small slopes approximation at the beginning of the motion, the normal velocity of the trench surface would be approximately the change in radius over time and equal to the velocity of Equation (15): The first derivative of parameter a 1 (t), solving Equation (16), can be expressed as: The solution of this simple ODE yields: where τ is a growth constant. Although the initial value of a 1 is unknown, the critical transformation wavelength can be obtained by examining Equation (18). The critical wavelength can be defined as the dominant wavelength over all the other potential wavelengths at which the transformation can be initiated. The highest exponential value in Equation (18) causes the fastest perturbation growth. Therefore, the minimum characteristic growth constant, τ, is the value that would lead to the fastest perturbation growth: The second derivative of the growth constant produces, with critical wavelength transformation, a positive value of 4ρ 2 /π 2 . Thus, the calculated value is a minimum and produces the fastest growing sinusoidal perturbation. On the other hand, other wavelengths, apart from the critical ones, can also produce a perturbation growth. The range of wavelengths that can lead a transformation is essential to obtain further transformation characteristic values. For a successful void evolution, the growth constant τ should be positive according to the surface diffusion dynamics: Another condition must be considered: surface energy minimization. The total surface area of the final geometry must be smaller than the initial cylindrical structure. If the volume is preserved, the radius of the final sphere would be (see Figure 2): Accordingly, the minimum required wavelength for surface minimization would be: Concerning the minimum wavelength, the condition of Equation (22) is less restrictive than the calculated condition of Equation (20). Thus, Equation (22) is not a sufficient condition for the geometrical evolution.
The last perturbation analysis step is to prove whether there exists a critical value at which the minimization of surface energy is maximal. The reduction in surface energy of the void shape evolution transformation is expressed as: This shows that the reduction in surface energy increases with increasing wavelength. There are no further critical values between the minimum and infinite wavelengths. Therefore, the critical wavelength calculated in Equation (19) is maintained as the preferred wavelength for the transformation.
The values obtained in this section are valid for infinitely long continuous cylindrical structures. However, in the real world, there are differences in comparison with this proposed ideal case.
The first difference would be the effect of the crystal orientation of the wafer and the crystal orientation obtained after the etching process. It is known that the transfer of atoms of the void shape evolution in silicon occurs between discrete crystallographic facets. For example, in the work of Sudoh et al. [16], the evolution of the trenches etched on a Si(001) substrate was carried out by the evolution and combination of {100}, {110}, and {113} discrete facets. Furthermore, the crystal orientations have different surface energies, ranging from 1.21 to 1.72 J/m 2 [44][45][46][47], depending on the source of data. As a consequence, a more complete description would take into account the spatial variation of the surface energy per facet and would require a discrete method to describe every possible facet in the system. This would need a modification and much added complexity of the surface diffusion Equation (7). Furthermore, the use of a discrete method such as the Wulff-shape-based method form by Kitayama et al. [48] for alumina would also be required. Unfortunately, that would neglect the advantages of using a continuous method, such as the perturbation analysis, to reduce the variables involved in the suggested model of this work. For this reason, while there is a dependency of the evolution on initial crystal orientation and created facets after the etching process, the presented model assumes a continuous and isotropic surface energy approach.
The second difference is the finite nature of the real structures. This causes a greater effect as this determines the pinch point of the separation of the internal void to the outer atmosphere along with the number of empty spaces in silicon that can be created. Consequently, further corrections must be performed to account for the top and bottom effects of the initial geometry. This is discussed and corrected in the following subsection.

Analytical Model of the Void Shape Evolution for Finite Structures
The adaptation of the infinitely large cylinder to a finite cylinder must be carried out by observing the literature experimental values and their reported real evolution patterns. The starting point of the analysis is the most probable wavelength values that create a full geometrical evolution, i.e., between the minimum and the critical one: The transformation wavelength must lie between the minimum and the critical (most probable) one. A cylindrical trench depth, which completely transforms into a spherical void, would have a depth equal to the transformation wavelength. If the depth of the trench has, at least, the same size as the minimum wavelength, the geometry would evolve to a complete sphere. As the depth size increases, the wavelength would also accordingly increase due to higher energy minimization until reaching the critical one, λ C , as determined by the kinetics.
According to literature observations (see Figure 3), the existence of the critical wavelength has been confirmed by the experiments of Sato et al. [2]. However, as opposed to the analyzed geometry, rectangular base trenches were used instead circular ones. This should not make a big difference as, as shown in their experiments, the rectangular base was transformed into a circular base, after some annealing time, with the same internal volume. To check the critical wavelength, an experiment was performed with the top of the geometry blocked by an SiN layer [2]. As a consequence, the evolution starts directly from below. The base was 0.14 × 0.24 µm, which corresponds to an equivalent cylindrical base radius of about 0.10 µm. The critical wavelength of such equivalent cylinder would be 0.89 µm, which is very similar to the value obtained experimentally (between 0.94 and 1.03 µm). Thus, the rectangular bases could be considered as circular ones with the same base area within a reasonable error margin. With regard to the longitudinal transformation of the work from Sato et al. [2], it can be observed that, when the critical wavelength is reached, the equivalent trench depth is transformed into a complete spherical void. Thus, the critical wavelength is the maximum length, in the absence of other limitations, which allows a complete void transformation. The difference of infinite and one-sided transformation of trenches, as well as the verification of the critical radius, was also tested and analyzed by simulations and experimental results of different etched silicon millefeuilles by Hernández et al. [49] and Garín et al. [50].

Analytical Model of the Void Shape Evolution for Finite Structures
The adaptation of the infinitely large cylinder to a finite cylinder must be carried out by observing the literature experimental values and their reported real evolution patterns. The starting point of the analysis is the most probable wavelength values that create a full geometrical evolution, i.e., between the minimum and the critical one: The transformation wavelength must lie between the minimum and the critical (most probable) one. A cylindrical trench depth, which completely transforms into a spherical void, would have a depth equal to the transformation wavelength. If the depth of the trench has, at least, the same size as the minimum wavelength, the geometry would evolve to a complete sphere. As the depth size increases, the wavelength would also accordingly increase due to higher energy minimization until reaching the critical one, λC, as determined by the kinetics.
According to literature observations (see Figure 3), the existence of the critical wavelength has been confirmed by the experiments of Sato et al. [2]. However, as opposed to the analyzed geometry, rectangular base trenches were used instead circular ones. This should not make a big difference as, as shown in their experiments, the rectangular base was transformed into a circular base, after some annealing time, with the same internal volume. To check the critical wavelength, an experiment was performed with the top of the geometry blocked by an SiN layer [2]. As a consequence, the evolution starts directly from below. The base was 0.14 × 0.24 μm, which corresponds to an equivalent cylindrical base radius of about 0.10 μm. The critical wavelength of such equivalent cylinder would be 0.89 μm, which is very similar to the value obtained experimentally (between 0.94 and 1.03 μm). Thus, the rectangular bases could be considered as circular ones with the same base area within a reasonable error margin. With regard to the longitudinal transformation of the work from Sato et al. [2], it can be observed that, when the critical wavelength is reached, the equivalent trench depth is transformed into a complete spherical void. Thus, the critical wavelength is the maximum length, in the absence of other limitations, which allows a complete void transformation. The difference of infinite and onesided transformation of trenches, as well as the verification of the critical radius, was also tested and analyzed by simulations and experimental results of different etched silicon millefeuilles by Hernández et al. [49] and Garín et al. [50].  The description of the process above is valid when the upper surface is passivated, for instance, with SiN. On the contrary, when there is no passivation, the void evolution is different. It has been observed [16] that the faster conversion occurs at the top of the cylinder where the void is closed and separated from the outer atmosphere. The intersection points at the top between the cylinder and the top surface are points of very large positive curvature (convex). The difference in curvature is higher between the top intersection points and the cylinder walls (cylinder walls have a negative curvature value) than at the bottom, where the intersection points between cylinder walls and the bottom surface of the trench have very large negative values (concave). Thus, the upper part evolves faster than the bottom. Formally, these intersections are points of theoretical infinite curvature and are, therefore, the source of the initial transformation. These differences in transformation speed are crucial for the definition of the presented analytical model.
Concerning the overall evolution, the starting of the void shape evolution would be the already-mentioned points of infinite curvature. After starting the transformation, the perturbation is propagated in the axial direction of the trench until achieving the condition of Equation (24). Starting from the top of the trench, the trench closes when the perturbation has grown until contacting the other side of the cylindrical walls. From the bottom infinite curvature points, a transformation would also occur in the opposite direction. The geometry starts to round up as the perturbation propagates to the center of the trench. The propagation of the top part of the trench has an opposite direction to the one occurring at the bottom, leading to different transformations and a kinetic competition.
Before continuing with the analysis, the boundary conditions must be defined. For this case, a regular matrix of trenches is selected, considering each trench as the unit structure which periodically repeats in all directions (see Figure 4).
for instance, with SiN. On the contrary, when there is no passivation, the void evolution is different. It has been observed [16] that the faster conversion occurs at the top of the cylinder where the void is closed and separated from the outer atmosphere. The intersection points at the top between the cylinder and the top surface are points of very large positive curvature (convex). The difference in curvature is higher between the top intersection points and the cylinder walls (cylinder walls have a negative curvature value) than at the bottom, where the intersection points between cylinder walls and the bottom surface of the trench have very large negative values (concave). Thus, the upper part evolves faster than the bottom. Formally, these intersections are points of theoretical infinite curvature and are, therefore, the source of the initial transformation. These differences in transformation speed are crucial for the definition of the presented analytical model.
Concerning the overall evolution, the starting of the void shape evolution would be the already-mentioned points of infinite curvature. After starting the transformation, the perturbation is propagated in the axial direction of the trench until achieving the condition of Equation (24). Starting from the top of the trench, the trench closes when the perturbation has grown until contacting the other side of the cylindrical walls. From the bottom infinite curvature points, a transformation would also occur in the opposite direction. The geometry starts to round up as the perturbation propagates to the center of the trench. The propagation of the top part of the trench has an opposite direction to the one occurring at the bottom, leading to different transformations and a kinetic competition.
Before continuing with the analysis, the boundary conditions must be defined. For this case, a regular matrix of trenches is selected, considering each trench as the unit structure which periodically repeats in all directions (see Figure 4). Following the described experimental observations and the defined boundary conditions, an analytical model based on modified perturbation analysis is constructed. The assumed steps from a cylindrical trench to a spherical void are: 1. The first transformation process is obtained through the competition between the effects at the top and at the bottom of the trench. 2. The top transformation has the priority. However, a balance between both extremes must be maintained, and this is used for obtaining the corresponding expressions (see Figure 5). 3. Once the wavelength at the top has been calculated, the pinch-off occurs at half the maximum wavelength, having the top surface as the reference (see Figure 5). 4. The void enclosed below the closing point transforms into one or several final spheres, depending on the aspect ratio of the void. The simplest transformation is Following the described experimental observations and the defined boundary conditions, an analytical model based on modified perturbation analysis is constructed. The assumed steps from a cylindrical trench to a spherical void are: 1.
The first transformation process is obtained through the competition between the effects at the top and at the bottom of the trench.

2.
The top transformation has the priority. However, a balance between both extremes must be maintained, and this is used for obtaining the corresponding expressions (see Figure 5).

3.
Once the wavelength at the top has been calculated, the pinch-off occurs at half the maximum wavelength, having the top surface as the reference (see Figure 5).

4.
The void enclosed below the closing point transforms into one or several final spheres, depending on the aspect ratio of the void. The simplest transformation is shown in Figure 5, where the complete enclosed void volume transforms into a spherical void geometry with the same volume.

5.
The volume left over the closing point transforms into a step, i.e., a difference of surface level with respect to the initial surface (see Figure 5). 6.
The only considered final equilibrium geometries are structures with constant curvatures: spheres and planes. No intermediate states can be predicted without a proper process simulation. 7.
Only the formation of one trench is considered, and no coalescence is assumed. This simplifies the expression of the size of the final void. The distance of the neigh-boring trenches is sufficient to have a void spherical transformation dominated by the aspect ratio instead of the distance to other trenches. However, the latter effect is not completely neglected, as the distance between trenches affects the step and SON-layer sizes. The described geometrical parameters are shown in Figure 5. R C is the cylindrical radius, L C is the depth, and D t is the distance between trenches. The proposed model gives different final equilibrium geometries depending on the initial aspect ratio of the trench L C /R C .
The presented calculations are based on cylindrical trenches. For rectangular trenches, the approach in [7] is assumed. If the rectangle has two sides l 1 × l 2 , a radius R C,eq of the circle with the same base area is employed: R C,eq = (l 1 l 2 /π) 1/2 (25) The first case to be addressed is for a depth lower than the minimum wavelength: Following the scheme of Figure 5, to have a void transformation, the perturbation wavelength that starts at the top of the trench λ top should have enough depth to exist as it has transformation priority according to the observed literature data. However, in this case, the depth is not sufficiently deep to achieve the transformation condition of Equation (24), and thus, there is no spherical formation. The equilibrium geometry at the end of the void shape evolution process would be just a flat plane without any internal void sphere, and that plane would be at a lower position than the initial substrate surface, i.e., a step is formed.
To calculate the step, it is assumed that the volume of the initial trench is transformed into the volume of the step. As a regular pattern of trenches was assumed (see Figure 4) and being A, the area of that unitary trench surface area (D t × D t ), the step size of surface vertical displacement of the substrate surface would be: That is the only geometrical parameter that can be calculated in this case with the lack of formation of any empty space in silicon. One example can be found in [3].
The second case of study is for 2πR C < L C < 4πR C . This case is represented in Figure 6. In contraposition to the previous case, the depth is sufficiently long to allow a geometrical evolution as it is longer than the minimum wavelength. As previously assumed, the transformation at the top has priority due to the observed kinetics and evolves at the minimum transformation wavelength. This is not higher due to the transformation occurring at the same time at the bottom in competition with the perturbation starting at the top. Once that the top part is closed, with an assumed closing point at half the wavelength point, starting from the top surface, the rest of the resting void volume automatically transforms into a sphere due to energy minimization. The proposed transformation is shown in Figure 6. Regarding the calculation of the final equilibrium geometries, the transformation at the top of the trench occurs at a wavelength of 2πR C . The closing point, which separates the outer atmosphere from the internal void, would be then at a distance of (λ top /2) with respect to the initial substrate surface (see Figure 6): As in the previous case, and for simplification, the step size would be calculated from volume preservation between the void volume of the trench over the closing point and the proportional part of the unitary trench area: Under these assumptions, the remaining void volume, which stays below the closing point, transforms into a void sphere with the time. Comparing the volumes of the initial trench below the closing point and the volume of a sphere, the spherical radius would be: The term (L C − πR C ) represents the part of the initial trench that stays below the closing point. The last geometrical feature would be the so-called Silicon-On-Nothing (SON) layer or T SON (see Figure 6). It is assumed that the void is situated at the middle point of the void volume enclosed below the resting point as similarly observed in the literature [51]. The SON layer has two contributions: (i) the distance between the step level and the closing point; and (ii) the transformation of the spherical void, i.e., the difference between half of the length of the volume below the closing point and the final spherical radius. Taking all the terms together from Equations (27), (29) and (30), the SON layer is given by the following expression (already reordered): The literature's experimental data about this transformation can be found in [3,16,51]. The third case of study is 4πR C < L C < (2 1/2 + 1)2πR C.
In distinction with the second case, the transformation wavelengths are different. The void shape evolution occurs, in finite trenches, as a competition of different wavelengths starting from the corner points at the top and at the bottom of the trench. Depending on the kinetics, one side is faster and more predominant than the other. In this third case study, the trench is sufficiently deep to have, at least, a minimum wavelength transformation at both sides: 2 × 2πR C . Due to the faster transformation of the upper part, the system tends to take the highest possible wavelength for the upper transformation, while staying in the minimum at the bottom part. In contraposition to the previous study case, there is evolution at the same time at the top and at the bottom. After the pinch-off of the trench, the remaining void volume below the closing point transforms into a sphere, as the remaining empty space is not long enough for making the transformation to two different voids (see Figure 7). Regarding the calculation of the equilibrium geometry parameters, the procedure is similar to the previous case. First of all, the closing point must be located at half of the top transformation wavelength (see Figure 7): The step size is calculated from the void volume over the closing point: The spherical internal void has a radius which equals a sphere whose volume is the remaining volume below the closing point: As in the previous case, the center of the final equilibrium sphere is assumed to be aligned to the central point of the resting void below the closing point. The SON layer, after taking into account all the contributions, is then: The literature data about this transformation can be found in [3]. The fourth case of study is (2 1/2 + 1)2πR C < L C < (3·2 1/2 + 2) πR C . This case is dominant when the total depth of the trench is equal to the sum of the critical and minimum wavelength or higher. The top of the trench transforms following the critical wavelength (kinetically faster than the other options from the perturbation analysis), 2 1/2 ·2πR C , while the bottom part starts to evolve with a wavelength equivalent to the rest of the trench length, with the maximum limit of the critical wavelength. However, this limit is not a determining factor for forming one void or more, as the upper part of the trench would close before having a complete bottom transformation. The determining factor for the upper depth range before achieving two empty spaces in silicon is that, after the pinch-off, the length of the remaining void volume is large enough to form two empty spaces. When approaching the limit where more than one empty space can be formed, the depth of the trench should be longer than two times the critical wavelength. At those conditions, the top part and the bottom part start to evolve with the critical wavelength. After the pinch-off, as the bottom part would be already transforming with the critical wavelength, the only necessary condition to form a second empty space is that, just below the closing point, the remaining length that it is not being transformed by the bottom wavelength, is equal to or higher than the minimum transformation wavelength. If this is not fulfilled, the complete rest of void below the closing point would transform into an equilibrium sphere with the same volume due to the lack of space to create two separate entities. Thus, the limit of this fourth case would be a trench depth whose remaining void, after the pinch-off λ' bottom , is as large as the sum of the critical and minimum wavelength (see Figure 8).
The step size is then: The spherical radius would be calculated by performing the equality of the resting void volume below the closing point and the spherical void: The center of the void below the closing point is assumed to be also the central point of the final spherical void, and, as in the previous case, the Silicon-On-Nothing layer is: The literature data about this transformation can be found in [3]. The described equations can be used as a direct calculation of final output geometries from initial trenches or, the opposite, with a non-linear system solver of algebraic equations, to obtain initial trenches from a desired final design void shape in silicon.

Results and Discussion
In this section, the results of the proposed simple model to obtain final equilibrium spherical geometries are shown. They are compared to examples found in the literature. In some cases, the sizes of initial and final structures are given, but there exist also examples where the sizes were estimated. The estimated values from pictures/diagrams are given with the symbol "~". Initial values are given in Table 1 for rectangular/square based  trenches and Table 2 for circular base trenches. The model case is dependent on the aspect ratio as explained in the methods subsection.  The literature values of Tables 1 and 2 are then used as an input for the given expressions in the methods subsection, according to the corresponding aspect ratio. The resulting final equilibrium values are expressed in Tables 3 and 4 for rectangular and circular base trenches, respectively. The absolute and relative errors are also given. It is important to notice that the uncertainty in the estimated values from the publications would lead to relatively significant relative errors due to size range and the lack of access to high-resolution pictures or scales.   To further illustrate the calculation process, the example number 6 is detailed: 1. l 1 = l 2 = 0.35 µm is obtained from Figure 4 of Sato et al. [3]. In that diagram, the value "a" is given which, according to Figure 2, is the side of the base squared side.

2.
The length of the trench would be estimated from the parameter "b" of Figure 4 [3], as the x-axis yields the aspect ratio "b/a". In this case, it is estimated to be L C = 1.98 µm.

3.
The next step is to calculate which of the proposed models must be used. That comes from the aspect ratio calculation L C /R C , which yields about 8.14. This value lies between 2π (6.28) and 4π (12.57) and thus is the model case 2. In other words, Equations (27)-(31) must be used for the calculation of the final parameters. 4.
The last unknown parameter is A. A is the surface area of the repeated pattern of the trench as shown in the colored area in Figure 4 of this work. The top trench area is a squared area as shown in Figure 2 of Sato et al. [3]. The information provided by Figures 2 and 3 and the work of Mizushima et al. [1] suggest that the distance between trenches "a + c" would be about 2a (which is equivalent to D t of Table 1) and the area A would be 2a × 2a or 2l 1 × 2l 1 = 0.49 µm 2 .

5.
As the base of the trench is a square, the equivalent cylindrical radius is given by Equation (25) and equals R C,eq = 0.1975 µm 6.
T STEP is calculated with Equation (29) and the values of R C and A. The obtained value equals T STEP = 0.16 µm 7.
R S is calculated with Equation (30) and the values R C and L C . The obtained value equals R S = 0.34 µm.

8.
T SON is calculated with Equation (31) and the values L C , R C , and A. The obtained value equals T SON = 0.80 µm.

9.
Calculated data are then added in Table 3 and compared to the results of Sato et al. [3], Figure 4, which, after an annealing time of 10 min at 1100 • C and 10 Torr, yields the diameter of the final equilibrium spherical void. The estimated measured radius can be inferred from such Figure 4 and is R S = 0.34 µm.
In case of using SEM data, the values are estimated from the scale given in the corresponding figure. The procedure and steps would follow the same structure. For this work, a simple Microsoft Excel file was employed for the calculations. However, programs with similar functionality should be sufficient for such simple equations.
With respect to the results of the rectangular base trenches, the estimations of the model are quite accurate. First of all, case 1 should not form any empty space in silicon, which is predicted also by the proposed algebraic model. The errors in the other different geometrical features are between 0.00 and 0.05 µm, which is a very reasonable result, considering the difficulty in obtaining the literature values. Regarding the spherical radius, an average absolute error of 0.014 µm (average relative error of 5%) is calculated, which is a quite accurate value for such a simple model. In the case of the Silicon-On-Nothing layer, an average of 0.015 µm (2.89% of average relative error) was calculated. The particularity of these results is that they are taken from the same literature source, under the same annealing conditions, and the model can equally reproduce them.
The results in Table 4 are less accurate. While the spherical radius is similar in experimental and calculated values, within a 0.01-0.02 µm (2-4% of relative error) margin, the other geometrical parameters are not so precisely predicted. The Silicon-On-Nothing layer absolute error is between 0.08 and 0.28 µm (7-33% of relative error), and the only step size which could be compared has an absolute error of 0.17 µm, which is quite large. These two geometries, while they are cylindrical like the model geometrical assumption, were not in equilibrium state when taking the data. For this reason, a higher deviation is expected. In the void volume, the difference is not so high as an average radius is calculated (they are actually ellipsoids), but it may produce a higher discrepancy in other geometrical parameters as the thicknesses of the different layers are not fixed yet. SON and step sizes are also typically more dependent on the most external side of the real trench pattern, which the model cannot completely consider.
In the Figure 9, the different results are represented and sorted by literature source. The lines are displayed for clarity and to distinguish between experimental and analytical values. It can be observed that the trend of the calculated SON layer for the less precise results is similar to the experimental taken values.  Table 3; Table 4 with respect to the proposed analytical model.
For such a simple model, the results for equilibrium geometries are within reasonable margins in most of the tested examples. The highest error is registered when compared to geometries not completely in equilibrium as the equilibrium values were completely estimated.
Regarding the aspect ratio of the initial trench and its transformation, Table 5 shows a summary of the different aspect ratios where one or more than one empty spaces (ESS) in silicon are created. Concerning the aspect ratio analysis, different results are observed. The first comparison is done with respect to the starting aspect ratio where, at least, one empty space is formed. As observed in Table 5, most of other literature values coincide with the analytical model of this work. The second greatest distance is found with the work of Grau et al. [17], where a different method (moving mesh with a virtual curvature function) was employed for the simulations and half a critical wavelength was the formation limit. In that work, a lower formation limit can be also observed for the formation of more than one empty space in silicon. The greatest distance is found in Garín et al. [50], which seems like an outlier. The differences may be larger due to the use of silicon millefeuilles with very long structures, changing their upper/bottom geometry configurations in comparison to the structures studied in this work. On the other hand, the analytical model gives a similar aspect ratio relation for the creation of more than one empty space in comparison to the work of Sato et al. [3]. Additionally, their experimental results are the most coincident with the values given in this work (see Figure 9), which are based on the same process.
Müller et al. [51] also studied the aspect ratio of the formation of void shape evolution structures in silicon. Similar aspect ratio limits in comparison with this work were found, between 2.76 and 8.28. Those are the results which better resemble the presented analytical values. However, Nichols et al. [52] predicted a theoretical aspect ratio of 7.2 for the formation of more than one ESS, but their approach and linear perturbation analysis were similar to the basis of the moving mesh method and virtual curvature calculation presented in Grau et al. [17]. Other similar literature values are given by McLean et al. [53,54], which were based on the evolution of cylindrical liquid lead inclusions in an aluminum matrix under a heating treatment; therefore, the surface orientation was different (top and bottom parts are equal as opposite to this work). In the case of Garín et al. [50], while there are some differences in the process, due to the study of large silicon millefeuilles, the formation limit of more than one empty of silicon is similar to this work. Finally, Stapley and Beevers [55], who studied the void transformation under an Ostwald ripening process, gave a quite different value for the aspect ratio limit, which would be caused by different phenomena.
After analyzing the data availability, the model gives results within a reasonable margin, taking into account the differences in literature results. The model is very easy to implement, and it is very fast. However, at this point, it is limited to the formation of just one spherical internal void without coalescence.

Conclusions
In this work, a novel simple predicting algebraic model is proposed for the calculation of final equilibrium structures after the annealing process of etched void trenches in silicon. The presented model gives the geometrical parameters of the process output, i.e., Silicon-On-Nothing layer, the spherical radius, and step size simply by using the initial trench geometry data such as the depth and radius of the trench, in addition to the distance between trenches.
The proposed calculations are inexpensive and can serve as reference design points for desired structures as evidenced in the comparison to the experimental literature values. The model could be used to extract information of the final structure from initial data possibilities or, by contrast, to solve the non-linear system of equations to obtain possible initial geometries from a desired final empty space.
The model is universal and should be able to equally predict similar void geometries in other materials undergoing a surface diffusion transformation, limited, in the current state, by the formation of one empty space without any coalescence.