C.PashartisM.J.vanSettenM.HoussaG.PourtoisDepartment of Physics and Astronomy, KU Leuven, Leuven, BelgiumImec, Leuven, Belgium

###### Abstract

Advancements in modern semiconductor devices increasingly depend on the utilization of amorphous materials and the reduction of material thickness, pushing the boundaries of their physical capabilities. The mechanical properties of these thin layers are critical in determining both the operational efficacy and mechanical integrity of these devices. Unlike bulk crystalline materials, whose calculation techniques are well-established, amorphous materials present a challenge due to the significant variation in atomic topology and their non-affine transformations under external strain. This study introduces a novel method for computing the elastic tensor of amorphous materials, applicable to both bulk and ultra-thin films in the linear elastic regime using Density Functional Theory. We exemplify this method with a-SiO_{2}, a commonly used dielectric. Our approach accounts for the structural disorder inherent in amorphous systems, which, while contributing to remarkable material properties, complicates traditional elastic tensor computation. We propose a solution involving the inability of atomic positions to relax under internal relaxation, near the boundaries of the computational unit cell, ensuring the affine transformations necessary for linear elasticity. This method’s efficacy is demonstrated through its alignment with classical Young’s modulus measurements, and has potential for broad application in fields such as Technology Computer Aided Design and stress analysis via Raman spectra. The revised technique for assessing the mechanical properties of amorphous materials opens new avenues for exploring their impact on device reliability and functionality.

###### keywords:

Density Functional Theory; Amorphous; Glasses; Bulk; Ultra-Thin Films; Surfaces; Elasticity; Stiffness; Elastic Tensor

## 1 Introduction

Amorphous materials are abundantly integrated into modern semiconductor technologies, assuming a diverse array of roles, from ‘passive’ functions like substrates or electrically insulating layers, to ‘active’ ones, involving charge carrier transport. For example, hydrogenated amorphous silicon (a-Si:H) plays a pivotal role in both photovoltaic solar cells and the control transistors of thin-film liquid crystal displays. This is attributed to its direct band gap, moderate electron mobility, and excellent optical absorption capabilities, see for example Refs. [1, 2]. The anticipated successor, amorphous InGaZnO_{4} (a-IGZO), employed in high-end display technologies, is reported to offer substantial advantages over traditional silicon-based thin-film transistor materials. Notably, its deposition temperature is compatible with polymeric substrates. These advantages encompass reduced rigidity, improved stability, and enhanced electron mobility, making it feasible to manufacture flexible thin-film transistors for large-area electronic applications. [3, 4]

Furthermore, in Si-based Complementary Metal-Oxide–Semiconductor (CMOS) technology, amorphous silicon dioxide and hafnium dioxide act as dielectric barriers to prevent charge carrier leakage.[5, 6] Amorphous metals such as Ti based alloys are also frequently employed in gate metals [7], and the effective work function of the device can be adjusted using thin layers of other metals, such as amorphous TiAlC [8, 9]. Additionally, amorphous bulk materials find extensive use as low-dielectric constant materials to reduce the parasitic capacitance of the CMOS addressing interconnection circuitry or as cost-effective fill materials. [10, 11]

As the transistor’s dimensional reduction continues to approach the physical limit, the electronic industry is being driven to transition to three-dimensional integration schemes to sustain performance gains (see Refs.[12, 13]). As a direct consequence the impact of the mechanical properties of material on both the device operability and performances are becoming of paramount interest. The convergence of requirements arising from the device’s intended operation and the decreasing thickness of functional material layers necessitates stringent control measures to ensure the structural integrity of these ultra-thin material layers and their associated interfaces. They are also introducing new (or enhanced) sources of stress or strain variability. The effect of strain engineering in layered materials is a well studied topic, with examples including the enhancement of the silicon charge carriers mobility in Si based nanooelectronics [14, 15, 16]. At the atomic level, the dynamics of interfacial effects (strain, defects, etc.) can strongly impact different aspects such as adhesion, charge carrier mobility, band alignment, or the pinning of the Fermi-level.[17, 18, 19, 20] Our earlier work illustrated that sub-5nm thick crystalline Ru layers exhibit an increased mechanical flexibility [21]. In general, a less stiff elastic tensor has different consequences in terms of device affected properties, which include adhesion of material layers [22] and enhanced diffusion [23, 24]. Amorphous materials, due to their intrinsic atomic disordered nature, offer the perspective to accommodate external strain sources differently than their crystalline counterpart would, and are potentially opening new engineering mitigation strategies. For this to happen, there is a need for understanding how the mechanical integrity of amorphous materials is being impacted by thinning. Unfortunately, measuring the elastic properties of a few nanometer thick film is a non-trivial task [25, 26, 27], but one that is quantifiable from ab initio. An important advantage of computing the elastic tensor through ab initio methods is its adaptability for use in a range of simulations and experiments, including applications such as Technology Computer-Aided Design (TCAD) and stress analysis with Raman spectroscopy [28]. In this paper, we report a methodology to compute the elastic tensor of both bulk and ultra-thin film amorphous materials from first-principles. We illustrate the practicality of our approach using amorphous silicon dioxide, a widely utilized dielectric material in Si CMOS technology.

The methodology for numerically evaluating the bulk elastic tensor from ab-initio for crystalline materials is widely known, with implementations such as ElaStic [29] and Pymatgen [30]. However, there is a large difference when treating an amorphous material model that is not accounted for with their ordered counterpart; the non-affine transformation under external strain and lack of periodicity. In addition, when compared to ordered crystalline systems, amorphous materials are atomically disordered materials; in which, the bond lengths, coordination numbers (CNs), and atomic densities vary wildly throughout the model used. The difficulty is then to build a representation that captures a statically relevant distribution of the atomic disorder, while relying on condensed matter physics formalism with periodic boundary conditions. This typically results in the use of large periodic unit cell models, which can often exceed over 100 atoms, in which atomic disorder is being built. Due to these differences in morphology with respect to ordered systems, different local energetic minima solutions can co-exist. Therefore, modeling amorphous systems requires sufficient sampling of atomic configurations large enough to provide statistically relevant results. In what follows, we give a brief review on the background of amorphous materials, followed by a discussion of the algorithm used to compute the elastic tensor from ab-initio, and finally the application to SiO_{2} bulk and ultra thin films.

### 1.1 The Origins of Elasticity and Plasticity of Amorphous Systems

From a first-principles perspective, the modeling of amorphous materials requires accounting for a large number of atoms to capture any semblance of information that is representative of the macroscopic properties. It is hence computationally difficult to properly assess the mechanical properties. Not only does the atomic configuration change within the same stoichiometric configuration, but the long-range hom*ogeneity and topological heterogeneity can also radically vary. In some instances, the atomic packing density can be very dense, while in others, there could be a formation of sub-lattices within the model (see Ref.[31]). Since many amorphous solids exhibit brittleness, much of the field of mechanical properties of glass has been concerned with determining the fracture or stress yield point in which the material breaks or enters the plastic regime.[32] These analyses have been performed using a variety of molecular dynamic approaches by simply applying varying normal strain/stress to a material system. The procedure corresponds to deformations normal to the material surface which result in the Young’s moduli. However, when considering ultra-thin films or shear deformations of an amorphous system, there is a caveat.

Due to the disordered nature of bonding in amorphous materials, they are quite heavily prone to changes in atomic position - typically constructing a new topological phase of the initial system. In other words, plastic-like rearrangements of the system are not uncommon. The bond lengths between atoms, the CNs per element, and the topology of voids can all greatly vary depending on the phase and the approach used to build the system. Under external strain, this makes amorphous materials a more ductile alternative to their crystalline counterparts, and one that is capable of reordering and reorganizing under external influence. In fact, metallic glasses have been observed to atomically rearrange under deformation to accommodate strain.[33]

In crystalline solids, we often deal with elastic and plastic deformations, where the latter is translated by the formation of dislocation defects. The shear stress required to push a system into a plastic regime, where any deformation introduces irreversible changes in shape, is typically referred to as the yield stress. The shear stress is the limit in which the deformation rate ($\dot{\gamma}$) goes to zero [34],

$\sigma_{\gamma}=\lim_{\dot{\gamma}\to 0}\sigma(\dot{\gamma})\hskip 2.84544pt.$ | (1) |

In contrast, since there is no recognition of traditional order in amorphous solids, it is difficult to define what a plastic transition is. This is especially so given every amorphous solid of the same stoichiometric composition is unique or a locally stable solution to the configuration of the system. In other words, a plastic deformation realistically yields another amorphous phase of the system. Therefore, any deformation to an amorphous system, whether it be large or small in magnitude, can push the system into a new configuration. Shear stresses on amorphous materials is an ongoing topic of investigation and one which is grounded in Falk and Langer’s early work on dynamics of viscoplastic deformation in amorphous solids, where they found that the atomic flow occurs in localized zones or clusters of molecules that shift to release stress [35]. These zones then can cascade local changes in the system, further changing the mechanical properties. In their system, they found linear elasticity up to 2% of shear strain, after which the onset of the plastic regime occurs. Their modelling work was verified to be observed in slowly sheared systems by Maloney and Lemaitre.[36] Though both of these systems are two-dimensional, they provided valuable information into the very slow changes that occur at the microscopic level of molecules and atoms, and one which we can observe in our work herein. An extensive review is performed by Bonn et al. [34] over the current state of theoretical models of shear plasticity and the corresponding yield stress.

With regards to the elastic regime, the glass phase transition temperature ($T_{g}$) typically has no consistent correlation to the stiffness of the glass, even though its mechanical properties is highly dependent on the shear modulus or viscosity at the transition point. For example, SiO_{2} and other refractory glasses, have high transition temperatures but are also typically defined by their low atomic packing density.[37, 38] Metallic glasses on the other hand, have been shown to have their transition temperature correlated with their elastic moduli by $T_{g}=2.5E$ .[39] Even the Young’s modulus widely varies from 5GPa for glassy water to 365GPa for tungsten-based metallic glasses.[40, 41] The structural and electronic contribution to the mechanical properties have been suggested to be correlated with the average bond energy and density (or coordination number) [31]. Given the variation of interatomic bonding options available to many of the known amorphous solids, the topology or intrinsic structure also contributes. A comprehensive review of various amorphous systems and their elastic properties and values can be found in Ref.[38].

In the linear elastic regime, local changes still occur and compared to crystalline systems, there is a non-affine linear elastic response to deformation. The domination of shear-elastic fluctuations has been heavily reported in theoretical studies [42], suggesting that the bond restructuring and phase changes occur even in the elastic regime. In these models, the particles are typically studied using displacements, forces, model potentials (such as the Leonard-Jones), and frequencies of the sound waves. Though these abstractions don’t formally constitute a real system, they yield valuable information, such as the heterogeneity of the local stiffness. The level of disorder (when compared to ordered structures) at the local regime can affect the transverse degrees of freedom of the stiffness resulting in broadening of the energy spectra or the existence of the so-called boson peak. [43, 44] This broadening is very clearly observed in Raman Spectra of amorphous solids (see a-Si for example).

The information that even at the atomistic level, the elastic components can vary, in particular, the shear components, leads to the realization that amorphous materials are extremely good at reordering bonds to distribute external stresses. In addition, these local elastic variations do not exist at the crystalline level, leading to the conclusion that the elastic tensor of a sufficiently large amorphous material is more akin to a global average of the local variations. Unfortunately, determining the local stiffness is difficult from an ab-initio viewpoint as the forces between atoms is not an ingredient for minimization or one that can be gleaned to reconstruct the internal atomic stresses from theorems such as the Virial Theorem. Since we have established the crux of employing a mechanical analysis at the atomic level is due to both the capability of amorphous materials to reorder and construct new phases, these must be taken into account in any elastic tensor generation procedure.

## 2 The Amorphous Elastic Tensor Method

Both experimentally and theoretically (ab-initio, molecular dynamics, etc.), one can deform a section of an amorphous solid and measure the stress response of the system, giving the Young’s modulus if the stress and strain are parallel. To compute the full elastic tensor, deformations are made to the model unit cell to perturb the system and determine either the change in energy or the change in stress. The former requires the second derivative with respect to the strain applied. In contrast, the latter case is related linearly to the strain applied tends to be more precise,

$\vec{\sigma}=\textbf{C}\vec{\epsilon}\hskip 2.84544pt.$ | (2) |

Where $\vec{\sigma}$ is the Voigt representation of the Cauchy stress tensor, C is the fourth-order Elastic tensor, and $\vec{\epsilon}$ is the Voigt representation of the Green-Lagrangian strain tensor. The algorithm from pymatgen [30] produces at most six unique deformations with four varying magnitudes per deformation. The magnitudes for the strains normal to the surface are $[-1\%,~{}-0.5\%,~{}0.5\%,~{}1\%]$ whereas the shear strains follow $[-6\%,~{}-3\%,~{}3\%,~{}6\%]$. These magnitudes are typically sufficient for crystalline systems, but for amorphous systems, one must still ensure to remain in the elastic regime.

As has been reviewed earlier, shear deformations are the linchpin that drive large variations of atomic displacement in amorphous materials. For instance, Figure1a) illustrates the evolution of the morphology of an amorphous SiO_{2} model after the minimization of its atomic degrees of freedom with a shear deformation of 3%. A colour spectrum is provided to demonstrate the amount of normalized displacement occurring per atomic site, accompanied with the direction of displacement as a green vector. It is observed that independent sections of the material shift in differing directions. For example, in the top facet of the unit cell, there is a clockwise-like shift in atomic positions. The resulting elastic tensor for such a system yields varying results since a different applied strain would nudge the system into a different energetic minima. As a result, the elastic stability criteria of Born [45, 46] is typically not satisfied, leading to negative elastic tensor’s eigenvalues and hence to a non-physical result.

At first glance, one might assume that either i) the applied strain is too large, which perturbs the system out of the linear elastic regime into a plastic one, or ii) that due to amorphous nature of the system, the configuration undoes the applied strain. The consequences on the Voigt-stress components of the elastic tensor in a normally applied direction (Figure2a)). At this stage of the discussion, we are focusing on reproducing the linear elastic response of the material, therefore neglecting any non-linear, viscoelastic, or any other mechanical variant. In this regime, the stress-strain relationship must be strictly linear. This assumption is reasonable (but should be used on a per case basis) as for instance, many glasses, such as silica (a-SiO_{2}), silicates, and aluminates, have been reported to exhibit elastic properties.[37]

The linear elastic regime is clearly broken when observing $\sigma_{1}$ in Figure2a), the system is responding to the external strain with non-affine transformations yielding a new local structural minimum. Additionally, recall that the relative magnitudes of each curve do not matter, but the variation of the stress is used to evaluate the elastic tensor. In crystalline materials, 1% strain does not force the material into a plastic regime for the majority of pure ordered materials. Is the non-linearity due to anharmonicity or is it truly due to the plastic limit? Instead of determining either of these options, if case ii) is correct, then the application of elastic tensor methodology to ultra-thin films or surfaces can be applied to amorphous materials.

If, due to applying strain to the unit cell of an amorphous material, the system undoes the applied strain, then we cannot use the finite perturbation technique in the quasi-harmonic potential limit to determine the mechanics of the material. One of the key-defining features of an amorphous material is its inherent lack of translational symmetry or order, so any enforced deformation to a unit cell without symmetry is therefore capable of reverting to its original state. This is very much the exact same phenomena observed when developing the elastic tensor methodology for ultra-thin films. The solution is to artificially build a region of symmetry that maintains the applied strain. To this end, taking a small volume ($dV$) within the edges of the unit cell periodic boundaries (see Figure3), atoms can be selected to remain invariant under internal relaxation. The internal relaxation of such a scheme can be clearly observed in Figure1b), where the atoms along the unit cell boundary do not move. The remaining interior atoms still relax their positions, but to a much lesser degree than in Figure1a). The result is an elastic tensor which typically satisfies the Born criteria for stability. The constraint also yields in a linear response to the applied strain, with respect to stress, as observed comparing the non-held system in Figure2a) to that of the held system in Figure2b).

The process to calculate an elastic tensor for an amorphous system then simplifies to:

- 1.
Relax the constructed amorphous system from a Monte Carlo, molecular dynamics or other scheme

- 2.
Apply the same deformations for crystalline systems, but enforce a shell of volume near the surface to constrain atomic sites lying within

The distance from the cell edges should be converged such that the number of atoms selected is minimized and the elastic tensor components remain stable and consistent. In practice we found a good metric was to ensure that less than 20% of the atoms should be constrained and that a very thin shell in the neighbourhood of 0.05-0.25Å (see Supplementary Material) is sufficient to yield an elastically stable tensor. Of course, general convergence to Density Functional Theory (DFT) parameters should still be carried out. There is no guarantee such a scheme correctly calculates the stiffness of an amorphous material, but by applying an Ackland-Jones analysis [48] on the local bonding nature of the structure, we observed that by holding certain atoms the symmetry is maintained under strain. Without the technique, the atomic shift ensued could signify that the material is in a non-elastic regime. We instead suppose that that the changes are due to an artifact of the lack of inherent symmetry in the computational unit cell. To assess how valid holding boundary atoms is, a-SiO_{2} was constructed and the elastic tensor analysis was performed.

To compute the ultra-thin elastic tensor refer to our previous paper [21], where the position of the atoms at the surface are held still (which is given for free by the amorphous algorithm). Additionally, the stress tensor must be scaled accordingly using a standardized definition of the thickness of the film which accounts for the electron density extending past its surface.

## 3 Computational Details

### 3.1 Generation of Amorphous Structures

The a-SiO_{2} models were generated with a decorate and relax scheme (see Ref.[49]). The procedure begins with an atomic structure built from a random distribution of silicon, the atomic position is then compressed and subsequently decorated by oxygen atoms to reproduced coordination numbers similar to the crystalline phase. These models are then optimized using DFT. This methodology has been heavily employed in the study of IGZOs and in general leads to structures with reasonable radial distribution functions and band gaps [50].

A total of ten 96-atom structures and fifteen 200-atom structures were generated to sample the variations in configurations available at densities ranging from 2.3–2.9$g\cdot cm^{-3}$ (see Figure4). Two different sizes of unit cells were studied, a 96-atom cell (Si_{32}O_{64}) and variants of a 200-atom type cell (ex. Si_{66}O_{132}). The radial distribution functions of the 200-atom type models for the nearest neighbours and different bond types can be found in Figure5, which reproduce values observed in literature [51, 52, 53].

### 3.2 Density Functional Theory

All calculations presented were performed using the first-principle CP2K simulation package, using 3D periodic boundary conditions. A Fermi Dirac distribution for the occupancy of the valence band structure has been imposed with an electronic temperature of 1000K. The Goedecker-Teter-Hutter CP2K pseudopotential library [55, 56, 57], was used and combined with the Perdew–Burke-Ernzerhof (PBE) [58] exchange correlation functional. DZVP basis sets were used from the CP2K standard library.[59] The simulations were converged with respect to the elastic tensor components for the bulk conventional structure to achieve a precision of 10GPa. The bulk lattice parameters and elastic coefficients were found to be converged with a kinetic energy cutoff of 600Ry and a relative energy cutoff of 50Ry (for the Gaussian grid). The maximum force cutoff used for unit cell relaxation was $1\cdot 10^{-4}$bohr^{-1}Ha. A single k-point was used to perform the calculation since amorphous materials have no Bravais lattice due to lack of symmetry. Convergence of DFT parameters was performed on a single amorphous system. The unit cells were minimized to reach their optimum lattice vector orientation and position.

### 3.3 Additional Details

The distance from the unit cell edges for the boundary atoms was converged to select the minimum number of atoms, at a distance of 0.3Å from the edge. For example, in the SiO_{2} system with 198 atoms, on average, 11% of all the atoms in the system were constrained after the deformations were applied.

The ultra-thin SiO_{2} films were constructed with the relaxed amorphous bulk variants, oriented with the direction of thickness aligning with the Cartesian-z axis. They were given a vacuum height of 20Å. The films were allowed to relax while maintaining the same lattice vector directions but not magnitudes, using same tolerances as the bulk before performing the necessary deformations. The atoms at the edge were constrained, combined with the amorphous methodology of constraining atoms near the boundaries of the unit cell.

For the definition of the thickness of the film, at the time of the research, the more precise thickness fitting to the Hartree potential (see Ref.[21]) was not implemented, therefore these results are reported with the vacuum removed. Note that the difference is observed to be minimal. As a consequence, the atomic densities are then reported without thickness correction and are strictly the edge-to-edge or vacuum removed for the z-axis of the volume.

Note that for the orientation of the elastic tensor, both in the bulk and ultra-thin film structures, the first lattice vector is colinear with the Cartesian x-axis. While in the bulk, the topological hom*ogeneity implies that the lattice vector orientations will not significantly matter, assuming a well-distributed amorphous configuration and a large enough configuration space was sampled for stiffness. The selection of the third lattice vector in the direction of the Cartesian-z for ultra-thin films allows those structures and their elastic tensor components in that direction to be easily compared without worry of a change of coordinates.

Coordination numbers were analyzed using the pymatgen CrystalNN method, which utilizes Voronoi analysis to ascertain the neighbourhood of atoms to select for bond pairs.[60] A recent study by Pan et al. 2021 [61] suggests it to be one of the best frameworks consistent with experimental evidence of coordination numbers.

## 4 Results & Discussion

### 4.1 Bulk a-SiO_{2}

The mechanical properties of a-SiO_{2} are widely varied due to the plethora of methods in which to synthesize it, such as the deposition process, and varying temperatures or pressures in the growth in environment. Its characterization produces differing results, which is understood to be due to the differences in structural configuration natural to amorphous materials. Depending on the strength of the covalent bonds, the orientations, their number, the coordination of each atom, and the number and locations of voids, the mechanical properties can be expected to change. The variation in the density and its associated shear, bulk, and Young Hills’ moduli are shown in Figure6. For reference, densities close to 2.2 $g\cdot cm^{-3}$ are used in devices. The best-fit was computed for the data and tends towards the values previously computed with ab-initio from Ref.[62], without the elastic tensor. Experimental values of a-SiO_{2} have found the Young’s modulus to be 70GPa [37] for bulk silica and 76.6$\pm$7.2GPa [63] for nano wires, measured at room temperature. Extrapolated to 0K, the experimental Young’s modulus for bulk silica is near 65GPa (note the density from this reference is not reported), placing the trends of our results and those of other DFT calculations within good accuracy of experimental data.

A linear best-fit line was computed to demonstrate the general trend of the data, as it is too sparse to fully determine if the relation is linear or non-linear. The R^{2} values of each line of fit ranges between 0.55-0.56, but a $\chi_{red}^{2}$ cannot be computed due to the lack of a quantity of variance. Traditionally, the Young’s modulus and the atomic density are given on a log-log plot, typically called a material selection chart (see Ref.[64]), but this does not imply that any material is non-linear with respect to the density and modulus. For amorphous materials, there have been reports of linear scaling with the two properties, one example being Harms et al. 2003 [65], where different annealing temperatures generated amorphous structures of different densities and stiffness. Therefore, any best-of-fit lines reported in this amorphous work should be considered general trends and not necessarily validation of linear behaviour.

The values reported by Bondi et al. 2010 for a-SiO_{2} consist of up to 100-atom amorphous configurations where the moduli were computed using the standard technique of straining in a particular direction (or isotropically) and the subsequent stress response is measured. Namely, the Young’s modulus is typically computed as

$Y=\frac{\sigma_{xx}}{\epsilon}\hskip 2.84544pt,$ | (3) |

with $Y$ being the Young’s modulus, $\sigma_{xx}$ the stress in the xx component (or 11), and $\epsilon$ the magnitude of strain in the x direction.In contrast, the Young’s modulus for this work uses the components of the elastic tensor,

$Y=\frac{K_{VRH}\cdot G_{VRH}}{3K_{VRH}+G_{VRH}}\hskip 2.84544pt,$ | (4) |

where $K_{VRH}$ and $G_{VRH}$ are the Voigt-Reuss-Hill average bulk modulus and shear modulus. The shear components of the elastic tensor are included in this formulation, suggesting that the methodology applied to compute the tensor is within the variation of existing literature. At the higher densities, we observe a variation of approximately $\pm 20$GPa for the Young’s modulus and $\pm 10$ GPa for the bulk modulus. This range of stiffness is also observed by Bondi et al. [62] (with a lower density). They suggested, when compared to a-Si, that the addition of H or O to the system weakens the bonds, creating a more deformable system, while reducing the coordination numbers, demonstrating another key-defining attribute of less stiff amorphous systems.

#### 4.1.1 Analysis of Bonds and Coordination

In an effort to determine the source of the changing elastic properties for a given density, a coordination analysis was performed on each structure’s relaxed state. The variation of the average coordination number with the Young’s modulus in Figure7 suggests that as the density is increased, the coordination number typically increases. This is observed from the colours becoming lighter as the density increases. Except for one or two coordination pairs, the majority of the pairs are Si-O. The result coincides with the literature that suggests that given any atom, the more bonds it has (higher CN), then the more resistant it should be to external stresses, thereby stiffening the material system. [31] Much like a simple spring system or ball and chain model of atoms, the more connections added to an atom, the more difficult it is to perturb it. The result does assume that the chemical environment, or strength of the bonds, remains approximately the same throughout each configuration computed. The change in average CN is rather small, on the order of 0.15 at the ends of the spectra. Within a small unit of density, the average CN does indeed change non-uniformly in some cases, but the overall trend towards higher stiffness and increased average CN remains.

To determine whether bond direction enforces stiffness in SiO_{2}, the CN result was further post-processed to yield the pseudo-bonds or vectors to the nearest neighbours, each vector was unique to avoid double counting of the same pair (only one of Si-O or O-Si). After taking the element-wise absolute values of the vector, it can then be projected onto the Cartesian axis to determine the portion of the bond in the $\hat{x},\hat{y},\hat{z}$ directions. Enforcing that each component is positive, is necessary to ensure that each vector adds constructively and not destructively which would yield cancellation. If these vectors are normalized to ignore bond length bias, then averaged, the resulting gradient in Figure8 and 9 is observed (colour scale). In particular, the averaged and standard deviation normalized bond direction components are shown with the $c_{11},c_{22},c_{33}$ (xx, yy, zz) elastic tensor components over the atomic densities. We note that in a sufficiently statistically sampled and large enough amorphous material such as SiO_{2}, the mechanical properties are expected to be isotropic. Since the $c_{11}$ and $c_{22}$ components have similar best-fit trends, but the $c_{33}$ component is different, this highlights the difficulty in sampling enough topologies and varied densities to properly determine the mechanical properties of amorphous materials.

In Figure8, there does not appear to be correlation over the spectrum of densities computed with regards to the elastic tensor component and the pseudo-bond in the same direction. If the average bond strength remains similar across each structure over every atom, one would expect the respective component of the tensor to be stiffer with more bonds along the same direction. To observe this, for any small width of atomic density, the data below the best linear fit should be darker in colour than those above. This suggests that either insufficient structures have been used to properly gather the necessary statistics, or that the direction of the bonds has no effect on the macroscopic stiffness of the amorphous material. The former being the obvious choice since amorphous materials restructure and compensate for external sources of stress.

The variation in average bond direction along each component is very small among the calculated structures, inferring that the amorphous configurations is well disordered. The standard deviation of the bond direction within each structure is shown in Figure9. The values suggest standard deviations on the order of 2% of the projection of the pseudo-bond components, which is substantially small, further confirming that the amorphous structures are sufficiently disordered. The lack of variation does indeed suggest that the average pseudo-bond directions for each structure do not affect the final macroscopic stiffness of the material. To further support this observation, the average bond length variation over the densities computed is $10^{-3}$Å, which is within numerical precision. Meaning that variations for the Si-O bonds within a given structure are extremely small or negligible.

### 4.2 Ultra-thin a-SiO_{2}

### 4.3 Results & Discussion of Ultra-Thin Amorphous Silicon Dioxide

What has been discerned thus far is the dependency of SiO_{2} on both the bond strength and the coordination number. Extending such a study to ultra-thin silicon oxide is valuable due to its prevalence in CMOS transistors, where they are typically on the order of sub-2nm. Decoupling the dielectric layer to study individually as a free-standing film enables us to study its individual mechanical properties, which is presented here. At the time of this research, it is not clear if mechanical properties of ultra-thin SiO_{2} have been computed or derived experimentally.

Combining the approach of ultra-thin film and the amorphous elastic tensor methodologies, the elastic components of the different amorphous configurations (from their corresponding bulk structure) are calculated in Figure10. The films after cell optimization range from 1-1.6nm in thickness (surface to surface distance). The statistical averages are shown as the large dots, with the values of the many configurations given by the smaller, lighter dots. A guiding line is drawn between these averages to demonstrate the expected decrease for any configuration, on average, of a given tensor component. It is evident that, like the metallic and crystalline ultra-thin film systems, the most affected component in magnitude is that of $c_{33}$ of 50%- the one related to stress or strain applied normally to the surface. However, both the $c_{13}$ (43%) and the $c_{23}$ (48%) components demonstrate large percentage deviations relative to the bulk average, as observed in Table1. Which are responsible for the in-plane stiffness due to a normally applied stress. Any deviations of the tensor components with smaller magnitudes should be considered carefully, as the precision of the calculation is on the order of 10GPa, representing a 20% potential error. Note that due to selecting the thickness to be the surface-to-surface distance, the values of the elastic tensor will migrate to a larger thickness (rightward on the figure) and decrease in stiffness by the equivalent amount of thickness added. This is due to the stress being inversely proportional to any elastic tensor component.

Avg. Bulk (GPa) | Avg. Film (GPa) | % $\Delta$ | |
---|---|---|---|

$c_{11}$ | 158 | 134 | -15 |

$c_{22}$ | 161 | 129 | -20 |

$c_{33}$ | 162 | 81 | -50 |

$c_{12}$ | 50 | 38 | -24 |

$c_{13}$ | 50 | 28 | -43 |

$c_{23}$ | 49 | 26 | -48 |

$c_{44}$ | 50 | 34 | -31 |

$c_{55}$ | 48 | 36 | -26 |

$c_{66}$ | 52 | 42 | -19 |

In ultra-thin film systems (see Ref.[21]), dangling bonds at the surface contribute to the change in stiffness; increased strain at the surface was observed relative to where the surface would be if it were fully coordinated. Such an analysis cannot be performed easily for amorphous thin films, since the atomic configuration is widely varied, and the atoms drastically reorder themselves to accommodate for dangling bonds. Studying the change in average coordination of the top-most and bottom-most atoms within 1Å from the average coordination of the entire film suggests that indeed this is also the case in amorphous ultra-thin films between 1-1.6nm. The change in coordination number was found to be on average -0.9 with an associated standard deviation of 0.3, meaning that compared to the entirety of the ultra-thin film, there are on average atoms with one less bond at the surface. Thus it can be considered that the surface has dangling bonds. Additionally, the films were found to have expanded on average 8% times (standard deviation of 6%) in thickness when compared to its initial starting position from the fully coordinated equivalent. Combined, the lowering of the coordination at the surface and the stretching of the system in the direction of vacuum, suggest that even in amorphous ultra-thin films (compared to crystalline ones [21]), the change in stiffness is due to the dangling bonds and to the relaxed degrees of freedom along the film confined direction. The introduction of the vacuum breaks the hom*ogeneity found in-plane throughout the bulk amorphous state, creating a more deformable system out-of-plane. We expect that at larger thicknesses, the surface effects will be overwritten by the contribution from the centre of the film body, outweighing the effect of the dangling bonds and surface relaxation, thereby increasing the stiffness of the system towards the bulk values.

## 5 Conclusion

As more amorphous materials are being included in semiconductor device stacks, particularly so as miniaturization continues, the need for multiscale properties increases to properly simulate and predict these systems. This paper has addressed the need to calculate the elastic tensor of amorphous solids, required to understand the effect of strain in various subjects such as thermally induced strain, strain-enhanced diffusion, or adhesion of layers. We have reviewed the origins and discussed the plastic and elastic regimes of amorphous materials. Given that a material can be described in the linear elastic limit, a methodology to compute the elastic tensor has been shown to be in agreement with literature for a-SiO_{2}. The methodology can be extended to any amorphous system in the linear elastic regime, but the efficacy should still be evaluated. Due to the nature of amorphous materials, small changes in external strain often result in a change the topological phase, often leading to plastic transformations. Since the ab-initio methodology requires internal relaxation after externally applied strain, this often leads to a sudden change in topology, breaking the affine transformation requirement of the linear elastic regime. To accommodate this issue, a thin shell of atoms around the boundaries of a computational unit cell are selected to remain stationary during internal atomic relaxation. It was found that for a small statistical sample of bulk a-SiO_{2}, the average coordinate number combined with the atomic density is correlated with the Young’s Hill modulus. Additionally, there is no nearest neighbour directional component to the corresponding direction of elastic tensor component (ex. $c_{11}$ with the Cartesian X direction). For ultra-thin film a-SiO_{2} constructed from the bulk counterpart, it was confirmed that similarly to crystalline materials, the dangling bond or electron density at the surface contributes to a reduced stiffness in the direction perpendicular to the surface.

## 6 Data availability

Pre-density functional theory relaxed unit cells of the bulk a-SiO_{2} have been deposited in the Cambridge Crystallographic Data Centre (CCDC) online database under the deposition numbers of 2313547-2313596, which can also be found searching for this publication.

## 7 Acknowledgements

The resources and services used in this work were partly provided by the VSC (Flemish Supercomputer Center), funded by the Research Foundation - Flanders (FWO) and the Flemish Government. The various members of our group for discussions and aid over the course of this research.

## References

- [1]T.Matsui, H.Sai, K.Saito, M.Kondo, High-efficiency thin-film silicon solarcells with improved light-soaking stability, Progress in Photovoltaics:Research and Applications 21(6) (2013) 1363–1369.
- [2]M.Stuckelberger, R.Biron, N.Wyrsch, F.-J. Haug, C.Ballif, Progress in solarcells from hydrogenated amorphous silicon, Renewable and Sustainable EnergyReviews 76 (2017) 1497–1523.
- [3]K.Nomura, H.Ohta, A.Takagi, T.Kamiya, M.Hirano, H.Hosono,Room-temperature fabrication of transparent flexible thin-film transistorsusing amorphous oxide semiconductors, nature 432(7016) (2004) 488–492.
- [4]T.Arai, T.Sasaoka, 49.1: invited paper: emergent oxide tft technologies fornext-generation am-oled displays, in: SID Symposium Digest of TechnicalPapers, Vol.42, Wiley Online Library, 2011, pp. 710–713.
- [5]M.Houssa, High k gate dielectrics, CRC Press, 2003.
- [6]B.Wang, W.Huang, L.Chi, M.Al-Hashimi, T.J. Marks, A.Facchetti, High-kgate dielectrics for emerging flexible and stretchable electronics, Chemicalreviews 118(11) (2018) 5690–5754.
- [7]T.Matsukawa, Y.Liu, W.Mizubayashi, J.Tsukada, H.Yamauchi, K.Endo,Y.Ishikawa, H.Ota, S.Migita, Y.Morita, etal., Suppressing v t and g mvariability of finfets using amorphous metal gates for 14 nm and beyond, in:2012 International Electron Devices Meeting, IEEE, 2012, pp. 8–2.
- [8]Y.Kamimuta, K.Iwamoto, Y.Nunoshige, A.Hirano, W.Mizubayashi, Y.Watanabe,S.Migita, A.Ogawa, H.Ota, T.Nabatame, etal., Comprehensive study of v fbshift in high-k cmos-dipole formation, fermi-level pinning and oxygen vacancyeffect, in: 2007 IEEE International Electron Devices Meeting, IEEE, 2007, pp.341–344.
- [9]H.Arimura, L.-Å. Ragnarsson, Y.Oniki, J.Franco, A.Vandooren, S.Brus,A.Leonhardt, P.Sippola, T.Ivanova, G.A. Verni, etal., Dipole-first gatestack as a scalable and thermal budget flexible multi-vt solution fornanosheet/cfet devices, in: 2021 IEEE International Electron Devices Meeting(IEDM), IEEE, 2021, pp. 13–5.
- [10]C.Wu, Y.Li, M.R. Baklanov, K.Croes, Electrical reliability challenges ofadvanced low-k dielectrics, ECS Journal of Solid State Science and Technology4(1) (2014) N3065.
- [11]A.Grill, S.M. Gates, T.E. Ryan, S.V. Nguyen, D.Priyadarshini, Progress inthe development and understanding of advanced low k and ultralow kdielectrics for very large-scale integrated interconnects—state of the art,Applied Physics Reviews 1(1) (2014).
- [12]S.Subramanian, M.Hosseini, T.Chiarella, S.Sarkar, P.Schuddinck, B.T.Chan, D.Radisic, G.Mannaert, A.Hikavyy, E.Rosseel, etal., Firstmonolithic integration of 3d complementary fet (CFET) on 300mm wafers, in:2020 IEEE Symposium on VLSI Technology, IEEE, 2020, pp. 1–2.
- [13]B.Vincent, J.Boemmels, J.Ryckaert, J.Ervin, A benchmark study ofcomplementary-field effect transistor (CFET) process integration optionsdone by virtual fabrication, IEEE Journal of the Electron Devices Society 8(2020) 668–673.
- [14]S.Gupta, W.T. Navaraj, L.Lorenzelli, R.Dahiya, Ultra-thin chips forhigh-performance flexible electronics, npj Flexible Electronics 2(1) (2018)1–17.
- [15]M.Wang, S.Cai, C.Pan, C.Wang, X.Lian, Y.Zhuo, K.Xu, T.Cao, X.Pan,B.Wang, etal., Robust memristors based on layered two-dimensionalmaterials, Nature Electronics 1(2) (2018) 130–136.
- [16]H.Jiang, L.Zheng, Z.Liu, X.Wang, Two-dimensional materials: From mechanicalproperties to flexible mechanical sensors, InfoMat 2(6) (2020) 1077–1094.
- [17]E.H. Poindexter, P.J. Caplan, Characterization of si/sio2 interface defectsby electron spin resonance, Progress in surface science 14(3) (1983)201–294.
- [18]C.R. Helms, E.H. Poindexter, The silicon-silicon dioxide system: Itsmicrostructure and imperfections, Reports on Progress in Physics 57(8)(1994) 791.
- [19]E.H. Poindexter, Mos interface states: overview and physicochemicalperspective, Semiconductor Science and Technology 4(12) (1989) 961.
- [20]M.Houssa, E.Chagarov, A.Kummel, Surface defects and passivation of ge andiii–v interfaces, Mrs Bulletin 34(7) (2009) 504–513.
- [21]C.Pashartis, M.van Setten, M.Houssa, G.Pourtois,Onthe elastic tensors of ultra-thin films: A study of ruthenium, AppliedSurface Science 592 (2022) 153194.doi:https://doi.org/10.1016/j.apsusc.2022.153194.

URL https://www.sciencedirect.com/science/article/pii/S0169433222007553 - [22]K.Sankaran, J.Swerts, R.Carpenter, S.Couet, K.Garello, R.F. Evans,S.Rao, W.Kim, S.Kundu, D.Crotti, etal., Evidence of magnetostrictiveeffects on stt-mram performance by atomistic and spin modeling, in: 2018 IEEEInternational Electron Devices Meeting (IEDM), IEEE, 2018, pp. 40–5.
- [23]N.Cowern, P.Zalm, P.Vander Sluis, D.Gravesteijn, W.DeBoer, Diffusion instrained si (ge), Physical Review Letters 72(16) (1994) 2585.
- [24]B.Yildiz, “stretching” the energy landscape of oxides—effects onelectrocatalysis and diffusion, Mrs Bulletin 39(2) (2014) 147–156.
- [25]X.Deng, V.R. Joseph, W.Mai, Z.L. Wang, C.J. Wu, Statistical approach toquantifying the elastic deformation of nanomaterials, Proceedings of theNational Academy of Sciences 106(29) (2009) 11845–11850.
- [26]C.A. Amo, A.P. Perrino, A.F. Payam, R.Garcia, Mapping elastic properties ofheterogeneous materials in liquid with angstrom-scale resolution, ACS nano11(9) (2017) 8650–8659.
- [27]S.Zak, C.Trost, P.Kreiml, M.Cordill, Accurate measurement of thin filmmechanical properties using nanoindentation, Journal of Materials Research37(7) (2022) 1373–1389.
- [28]E.Anastassakis, A.Pinczuk, E.Burstein, F.Pollak, M.Cardona, Effect ofstatic uniaxial stress on the raman spectrum of silicon, Solid StateCommunications 88(11-12) (1993) 1053–1058.
- [29]R.Golesorkhtabar, P.Pavone, J.Spitaler, P.Puschnig, C.Draxl,Elastic: A tool forcalculating second-order elastic constants from first principles, ComputerPhysics Communications 184 (2013) 1861–1873.doi:10.1016/j.cpc.2013.03.010.

URL http://dx.doi.org/10.1016/j.cpc.2013.03.010 - [30]M.deJong, W.Chen, T.Angsten, A.Jain, R.Notestine, A.Gamst, M.Sluiter,C.KrishnaAnde, S.Van DerZwaag, J.J. Plata, etal., Charting the completeelastic properties of inorganic crystalline compounds, Scientific data 2(1)(2015) 1–13.
- [31]L.Wondraczek, J.C. Mauro, J.Eckert, U.Kühn, J.Horbach, J.Deubener,T.Rouxel, Towards ultrastrong glasses (2011).
- [32]A.Furukawa, H.Tanaka, Inhom*ogeneous flow and fracture of glassy materials,Nature materials 8(7) (2009) 601–609.
- [33]G.Wang, N.Mattern, J.Bednarčík, R.Li, B.Zhang, J.Eckert,Correlation between elastic structural behavior and yield strength ofmetallic glasses, Acta materialia 60(6-7) (2012) 3074–3083.
- [34]D.Bonn, M.M. Denn, L.Berthier, T.Divoux, S.Manneville, Yield stressmaterials in soft condensed matter, Reviews of Modern Physics 89(3) (2017)035005.
- [35]M.L. Falk, J.S. Langer, Dynamics of viscoplastic deformation in amorphoussolids, Physical Review E 57(6) (1998) 7192.
- [36]C.E. Maloney, A.Lemaitre, Amorphous systems in athermal, quasistatic shear,Physical Review E 74(1) (2006) 016118.
- [37]T.Rouxel, Elastic properties and short-to medium-range order in glasses,Journal of the American Ceramic Society 90(10) (2007) 3019–3039.
- [38]W.H. Wang, The elastic properties, elastic models and elastic perspectives ofmetallic glasses, Progress in Materials Science 57(3) (2012) 487–656.
- [39]W.Wang, Elastic moduli and behaviors of metallic glasses, Journal ofNon-Crystalline Solids 351(16-17) (2005) 1481–1485.
- [40]A.Migliori, J.Sarrao, W.M. Visscher, T.Bell, M.Lei, Z.Fisk, R.G.Leisure, Resonant ultrasound spectroscopic techniques for measurement of theelastic moduli of solids, Physica B: Condensed Matter 183(1-2) (1993) 1–24.
- [41]M.Ohtsuki, R.Tamura, S.Takeuchi, S.Yoda, T.Ohmura, Hard metallic glass oftungsten-based alloy, Applied physics letters 84(24) (2004) 4911–4913.
- [42]A.Marruzzo, W.Schirmacher, A.Fratalocchi, G.Ruocco, Heterogeneous shearelasticity of glasses: the origin of the boson peak, Scientific reports 3(1)(2013) 1–7.
- [43]S.Taraskin, Y.Loh, G.Natarajan, S.Elliott, Origin of the boson peak insystems with lattice disorder, Physical review letters 86(7) (2001) 1255.
- [44]A.Chumakov, G.Monaco, A.Monaco, W.Crichton, A.Bosak, R.Rüffer,A.Meyer, F.Kargl, L.Comez, D.Fioretto, etal., Equivalence of the bosonpeak in glasses to the transverse acoustic van hove singularity in crystals,Physical review letters 106(22) (2011) 225501.
- [45]M.Born, On the stability of crystal lattices. i, Mathematical Proceedings ofthe Cambridge Philosophical Society 36(2) (1940) 160–172.doi:10.1017/S0305004100017138.
- [46]F.Mouhat, F.-X. Coudert, Necessary and sufficient elastic stability conditionsin various crystal systems, Physical review B 90(22) (2014) 224104.
- [47]A.Stukowski, Visualization and analysis of atomistic simulation data withovito–the open visualization tool, Modelling and simulation in materialsscience and engineering 18(1) (2009) 015012.
- [48]G.Ackland, A.Jones, Applications of local crystal structure measures inexperiment and simulation, Physical Review B 73(5) (2006) 054104.
- [49]D.Drabold, Topics in the theory of amorphous materials, The European PhysicalJournal B 68(1) (2009) 1–21.
- [50]M.J. van Setten, H.F. Dekkers, C.Pashartis, A.Chasin, A.Belmonte,R.Delhougne, G.S. Kar, G.Pourtois, Complex amorphous oxides: propertyprediction from high throughput dft and ai for new material search, MaterialsAdvances 3(23) (2022) 8413–8427.
- [51]R.Temkin, An analysis of the radial distribution function of siox, Journal ofNon-Crystalline Solids 17(2) (1975) 215–230.
- [52]W.Ching, Theory of amorphous si o 2 and si o x. i. atomic structural models,Physical Review B 26(12) (1982) 6610.
- [53]S.Susman, K.Volin, D.Price, M.Grimsditch, J.Rino, R.Kalia, P.Vashishta,G.Gwanmesia, Y.Wang, R.Liebermann, Intermediate-range order in permanentlydensified vitreous sio 2: A neutron-diffraction and molecular-dynamics study,Physical Review B 43(1) (1991) 1194.
- [54]K.Momma, F.Izumi, Vesta 3 for three-dimensional visualization of crystal,volumetric and morphology data, Journal of applied crystallography 44(6)(2011) 1272–1276.
- [55]S.Goedecker, M.Teter, J.Hutter, Separable dual-space gaussianpseudopotentials, Physical Review B 54(3) (1996) 1703.
- [56]C.Hartwigsen, S.Gœdecker, J.Hutter, Relativistic separable dual-spacegaussian pseudopotentials from H to Rn, Physical Review B 58(7) (1998)3641.
- [57]M.Krack, Pseudopotentials for H to Kr optimized for gradient-correctedexchange-correlation functionals, Theoretical Chemistry Accounts 114(1)(2005) 145–152.
- [58]J.P. Perdew, K.Burke, M.Ernzerhof, Generalized gradient approximation madesimple, Physical review letters 77(18) (1996) 3865.
- [59]J.VandeVondele, J.Hutter, Gaussian basis sets for accurate calculations onmolecular systems in gas and condensed phases, The Journal of chemicalphysics 127(11) (2007) 114105.
- [60]S.P. Ong, W.D. Richards, A.Jain, G.Hautier, M.Kocher, S.Cholia,D.Gunter, V.L. Chevrier, K.A. Persson, G.Ceder, Python materials genomics(pymatgen): A robust, open-source python library for materials analysis,Computational Materials Science 68 (2013) 314–319.
- [61]H.Pan, A.M. Ganose, M.Horton, M.Aykol, K.A. Persson, N.E. Zimmermann,A.Jain, Benchmarking coordination number prediction algorithms on inorganiccrystal structures, Inorganic chemistry 60(3) (2021) 1590–1603.
- [62]R.J. Bondi, S.Lee, G.S. Hwang, First-principles study of the mechanical andoptical properties of amorphous hydrogenated silicon and silicon-rich siliconoxide, Physical Review B 81(19) (2010) 195207.
- [63]H.Ni, X.Li, H.Gao, Elastic modulus of amorphous sio 2 nanowires, AppliedPhysics Letters 88(4) (2006) 043108.
- [64]H.Shercliff, M.Ashby,Elasticstructures in design, in: Reference Module in Materials Science andMaterials Engineering, Elsevier, 2016.doi:https://doi.org/10.1016/B978-0-12-803581-8.02944-1.

URL https://www.sciencedirect.com/science/article/pii/B9780128035818029441 - [65]U.Harms, O.Jin, R.Schwarz, Effects of plastic deformation on the elasticmodulus and density of bulk amorphous pd40ni10cu30p20, Journal ofnon-crystalline solids 317(1-2) (2003) 200–205.