Influence of CH 2 content and network defects on the elastic properties of organosilicate glasses

Organosilicate glasses (OSGs) are used as low- k intermetal dielectrics for advanced integrated circuits. In this application, the material must fulﬁll two conﬂicting requirements: It has to have low density to reduce the dielectric constant while being sufﬁciently mechanically stable to withstand thermomechanical and other stresses during subsequent steps of integrated circuit manufacture. Recent experimental advances in improving the mechanical and electrical properties of these materials have not yet been systematically studied theoretically at the ab initio level due to the large model sizes necessary to realistically describe amorphous materials. In this paper we employ the density-functional-based tight-binding method to achieve an accurate description of OSG properties at different compositions. We analyze the inﬂuence of composition and local network defects on the density and bulk modulus of nonporous OSG. We ﬁnd that the dependence of density and that of mechanical stiffness on chemical composition are of different natures. This difference is traced to a transition between mechanisms of elastic deformation in silica glass and in silicon hydrocarbide, which is also the reason for the two materials’ different sensitivities to network defects.


I. INTRODUCTION
The ever-increasing complexity and clock frequencies of integrated circuits (ICs) lead to fundamental problems in chip design. One of these is the limited speed of signal propagation between switches. The inductance and capacitance of the metal interconnects increase the time necessary for a logic state, i.e., voltage level, to be communicated between the transistors of an IC. This signal propagation time is one of the factors limiting the maximum clock frequency at which an IC can operate, along with the switching time of the transistors. While the induction is mainly influenced by the design, the capacitance of the interconnects is also governed by the dielectric constant k of the insulating dielectric around the wire. For a long time, SiO 2 was the insulator of choice because of its favorable mechanical properties, its high breakdown field, and its ease of manufacture. However, technological advances in microelectronics have made the use of low-k materials necessary. One important class of such materials are organosilicate glasses (OSGs). These materials are often manufactured by annealing of siloxane precursors, such as octamethylcyclotetrasiloxane at moderate temperatures up to 600 K, or by plasma-enhanced chemical-vapor deposition. These processes lead to an imperfect topology of the covalent bond network with a large number of terminating methyl groups. Furthermore, porogens are often added to further increase the density and thus the k value of the resulting OSG. While these materials are already widely employed in the manufacture of complex logic devices, the origins of their properties are not yet well understood, especially at the nanometer scale. This lack of insight makes it difficult to improve the materials' properties in an effective, knowledgebased manner.
Some attempts have been made to understand the atomistic structure of OSG by using empirical interatomic force fields, but they suffer from limitations in the construction of their atomistic models. For instance, Yuan et al. 1 employed force fields that do not allow changes in the bonding topology, which limits the model search space. Previous ab initio studies of OSG have focused on spectroscopic properties of individual structural features. 2,3 In this study, we rely on simulation methods that do not force us to choose between size and topology constraints. We also use a different approach to building our models: instead of trying to replicate the manufacture process from siloxane precursors by assembling a model of the complex material by combining various building blocks, we regard amorphous silica as well as OSG as variations on amorphous silicon. In this picture, silica derives from an amorphous silicon network in which O atoms have been inserted into the Si-Si bonds. A dense OSG is then a modification of silica in which some of the -O-bridges are replaced by CH 2 groups. Network defects are modeled as vacancies in the Si superstructure with appropriate termination of the surrounding bond network. This ensures that proper coordination of Si atoms by either CH n or O is retained and no dangling bonds are introduced. Removing an O or CH 2 bridge and terminating the resulting dangling bonds by methyl groups is precluded for steric reasons.
In Sec. II we describe the computational methods, and in Sec. III we discuss the different structural models used for this work. In Sec. IV we present our results, which we discuss in Sec. V.

II. METHODS
We calculate total energy and atomic forces of our structural models using the density-functional-based tight-binding (DFTB) method. 4,5 This approximation of density functional theory 6,7 (DFT) allows calculations with an accuracy approaching that of generalized gradient approximation [8][9][10] (GGA) DFT at much reduced computational effort by using a minimal (Slater-type 11 ) basis set and pretabulated Hamiltonian and overlap matrix elements. We use the well-tested to a series of models strained around the equilibrium cell volume. In the above equation, U is the total energy, V is the cell volume, V 0 and U 0 denote the equilibrium volume and energy,B 0 is the bulk modulus, and B 0 is the pressure gradient of the bulk modulus. We average over models of the same composition.

III. MODELS
The basis of all of our models is a perfectly coordinated 64atom supercell model of amorphous silicon, which was generated by careful simulated annealing using the environmentdependent interatomic potential (EDIP) potential. 14-16 Ten different models of amorphous SiO 2 were generated by expanding the a-Si model and adding O atoms with Si-O-Si angles of 150 • at random positions around the axes of each former Si-Si bond. The resulting initial models were then relaxed in atomic positions and unit cell dimensions using DFTB. Since all ten initial models relaxed into the same final configuration, we conclude that we have found the global minimum configuration for the present network topology. The resulting 192-atom amorphous silica cell is shown in Fig. 1(a).
As can be seen from the data in Table I, the DFTBcalculated properties correspond very well to previous ab initio calculations by van Ginhoven et al. 17 108.92(9.01) Bulk modulus B 0 (GPa) 47.4 for the radial distribution function of our silica model. No static property varies significantly from previous volume-optimized generalized gradient approximation density functional theory (GGA-DFT) calculations. The bond length and bond angle distributions as well as further peaks of the radial pair distribution function also correspond well to previous theoretical results. 17,22 The results from both simulation works also show very good agreement with experimental values for fused silica. We conclude that our 192-atom amorphous silica model represents silica glass very well. The calculated elastic properties agree reasonably well with experiment; the bulk modulus deviates by ∼30%, which is a well-known behavior of GGA calculations. We proceed from silica glass to organosilicates SiO 2(1−x) (CH 2 ) 2x by randomly replacing a fraction x of O atoms from the glass network with CH 2 groups, followed by a reoptimization of cell vector lengths and atomic positions. We generate five different models for each O replacement fraction: x = 0.1, 0.2, 0.3, 0.4, 0.5, and 0.7, plus one model of silicon hydrocarbide (x = 1). Figure 1(b) shows one of the x = 0.2 models. Table II shows our results for Si(CH 2 ) 2 . The density decrease with respect to silica is caused by two factors in equal measure: the CH 2 group is 2 amu or 12.5% lighter than an O atom, and the Si-C bonds are ∼0.2Å or 12% longer than the Si-O bonds.
To examine the influence of network defects, we remove one or two single Si atoms (labeled V Si ) or a Si-O-Si bridge (V 2Si ) from the network and replace the dangling neighboring O atoms with CH 3 groups. For each of these defect types, we generated three silica models with the Si vacancies in different positions, ensuring a uniform distribution in the models containing two separate Si vacancies (2V Si ). We then generate a full set of composition variants for each defective model, following the same procedure as for the fully coordinated networks. These single Si vacancies cannot be regarded as pores, as the terminating methyl groups fill the whole volume of the void arising from removing SiO 4 . Only the V 2Si and larger defects actually introduce empty volume into the system.
It should be noted that the defects in the network connectivity introduced in this way do not lead to the formation of occupied electronic gap states. Therefore any errors by defect-defect interaction that could be regarded as spurious do not influence our results.

A. CH 2 content dependence
The density plot in Fig. 2(a) shows that the densities of the OSG vary almost linearly with the O to CH 2 replacement fraction. The trend suggests a small nonlinearity that can be approximated by adding a third-order polynomial with nodes at x = 0, x = 0.5, and x = 1 to a linear transition between the silica and silicon hydrocarbide densities: We find the optimal parameters for this curve to be ρ SiO 2 = 2.20 g/cm 3 , ρ Si(CH 2 ) 2 = 1.67 g/cm 3 , and γ = 1.02. The  Tables I and II) indicate that our model of the composition dependence is reasonable.
The bulk moduli calculated under the constraint of cubic cell dimensions, plotted in Fig. 2(b), show a more complex behavior. Their values vary much more at each data point than the densities. More important, however, is the interpolation of bulk moduli for SiO 2(1−x) (CH 2 ) 2x composites. A linear regression does not yield a good description of the transition from silica to silicon hydrocarbide at all. As Fig. 2(b) shows, the bulk moduli of the mixed SiO 2 -Si(CH 2 ) 2 phases are consistently lower than a linear mixture of the form B 0 = (1 − x)B SiO 2 0 + xB Si(CH 2 ) 2 0 would predict. Adding a parabolic nonlinearity that only acts on the mixed phases with roots at x = 0 and x = 1 leads to a reasonably good description of the bulk moduli: To elucidate the origin of this nonlinear lowering of stiffness, we examine the elastic deformation mechanisms of silica and silicon hydrocarbide: specifically, we compare the behavior of bond angles and bond lengths between a-SiO 2 , Fig. 3, and a-Si(CH 2 ) 2 , Fig. 4. Comparing the behavior of the bond angles, Figs. 3(a) and 4(a), we find that the mean bond angles in a-Si(CH 2 ) 2 do not change significantly under isotropic strain. Comparing the bond length changes with hydrostatic strain, Figs. 3(b) and 4(b), we see that the C-Si bonds in a-Si(CH 2 ) 2 are compressed much more than the Si-O bonds in a-SiO 2 , the relative bond length compression being more than one order of magnitude larger in Si(CH 2 ) 2 than in

B. Local network defects
To test the influence of local network defects we repeat the simulated stress-strain experiments on the V Si , 2V Si , and V 2Si defect models described in Sec. III. All defect models lead to a noticeable decrease in density and bulk modulus compared to the dense phases. Figure 5 and Table  III show that the densities decrease with increasing mass of the removed atoms, as does the density nonlinearity. However, the influence of vacancy distribution and size differs in silica and in silicon hydrocarbide. In silica, the density reduction of a single defect is almost identical to the mass reduction of the unit cell, indicating only very small changes in the unit cell volume, regardless of the defect size. A second V Si , however, leads to a considerably larger density reduction than what is caused by the mass defect alone. A higher density of missing network links allows the silica network to expand. In contrast, in Si(CH 2 ) 2 , the density reduction induced by missing network links is always much larger than the mass defect. Here, even the lowest density of network defects that can be represented by our model allows for a considerable relaxation of the cell volume. This indicates that a dense network structure of silicon hydrocarbide is under considerable internal stress.
The size of the defect (V Si vs V 2Si ) does not have a noticeable influence on the bulk modulus of silica, while there is a definite dependence on the number of defect centers, as shown in Fig. 6. In Si(CH 2 ) 2 , introducing a single network defect leads to a considerable loss in stiffness, while a second defect that does not introduce a void into the cell has no significant influence. However, introducing an additional empty volume to the cell leads to a further considerable reduction in bulk modulus.

A. Density
Positive values of the nonlinearity factor γ in Eq. (2) mean that near pure silica and silicon hydrocarbide the behavior of the total system is governed by the majority compound. Small inclusions of the respective minority material cannot exert their full effect on the volume of the supercell. Their effect is in part compensated for by the buildup of internal strain. Figure 5 shows clearly that this effect is strongest in the dense material. The presence of network defects hinders the surrounding material from exerting its full constraining force on the inclusions. This is evident in the reduction of the nonlinearity with increasing defect size (cf. Table III). The results also show that the smallest density of defects we can model already has a very large impact on the density nonlinearity.
The different responses of OSGs of different composition to the defects shown in Table III can be understood in the context of the different deformation mechanisms of silica and silicon hydrocarbide. Internal stress can be effectively accommodated in silica by the flexible Si-O-Si angles [cf. Figs. 3(a) and 4(a)]. In Si(CH 2 ) 2 all stress must be accommodated by straining Si-C bonds, which is energetically expensive. When a Si(CH 2 ) 4 unit is removed, four Si-C bonds that could transmit considerable tensile forces are removed with it, allowing the whole material to expand around the vacancy. Statistical analysis of the Si-C bond distributions, shown in Fig. 7, reveals that the dense silicon hydrocarbide is subject to considerable internal strain, which is relieved by the introduction of even one Si vacancy. While the mean bond length of 1.89Å and the standard deviation of 0.02Å remain unchanged, the skewness of the bond length distributions jumps from 0.16 to 0.39 upon introduction of one V Si unit and exceeds 0.44 when two Si atoms are removed. This clearly indicates that a significant number of Si-C bonds can relax to shorter distances in a defective network. In the dense network, many of the Si-C bonds are stretched due to the overly strong constraints by the very stiff Si and C sp 3 tetrahedra.

B. Elastic properties
The nonlinearity behavior in the elastic properties differs qualitatively from that of the density. While in the densities the "majority phases" govern the overall behavior, the bulk modulus is always closer to the properties of the weaker material, in this case silica. The question arises of why the surrounding majority phase does not lower the influence of small inclusions, as is the case for the density.
The comparison between Figs. 3 and 4 clearly shows that the origin of stiffness (or flexibility) differs between silica and silicon hydrocarbide. The flexible -O-bond angles open a new deformation mechanism that is not available in the rigid Si(CH 2 ) 2 network. Therefore, even small additions of such flexible links allow the network to accommodate strain while circumventing the unfavorable straining of Si-C bonds.
The presence of vacancy defects in the amorphous networks has a rather small effect on their overall properties. We attribute this mainly to the fact that the voids formed by single Si vacancies are almost completely filled by the terminating CH 3 groups, which hinders significant compression of the pore volumes. As can be seen in Fig. 6, the influence of the defects is much larger in Si(CH 2 ) 2 than in SiO 2 . This can be understood from the fact that a Si vacancy does not only remove bonds to be strained but also acts as a hinge, similar to an oxygen in the silica network but with almost zero energy cost for bending. The stronger influence of the V 2Si defect in silicon hydrocarbide can be traced to the ability of the material to relax into the pore upon compression. In silica, this opportunity is always present, due to the larger interstitial volumes. These are mostly filled by the hydrogens of the methyl groups in Si(CH 2 ) 2 .

C. Comparison to finite-temperature molecular dynamics
For a better comparison with experiment, we also performed finite-temperature molecular dynamics (MD) simulations on the same dense phase models, using the COMPASS 23 force field. 24 The results from these calculations are in good agreement with the static ab initio results. The trends in the composition dependence of density and elastic properties are confirmed in both simulations. However, two distinct differences are observed: silica is found to be stiffer in the MD simulations, while Si(CH 2 ) 2 is found to be more flexible, and the composition dependence shows a distinct density peak and a nadir of stiffness in the force-field simulations at low hydrocarbon admixture. These differences can be traced in part to the finite-temperature effects present in the MD simulations. The pure materials have different thermal expansion coefficients and different temperature dependence in their elastic constants. In addition to this, there is a disagreement between the DFTB results and COMPASS on the stiffness of Si-(CH 2 )-Si angles: while they appear extremely stiff in our DFTB calculations [cf. Fig. 4(a)], they are somewhat flexible in the force field. This accounts for some of the reduction in the stiffness of the Si(CH 2 ) 2 phase and also explains part of the initial loss of stiffness upon addition of small amounts of CH 2 , as additional deformation pathways are added to the network. At this time, it is unclear which of the two descriptions is more accurate. Experience shows that standard GGA tends to underestimate bulk moduli. Since the employed DFTB parametrization is fitted to GGA reference calculations, 12 the GGA reference results cannot be used to decide on this rather delicate detail.

VI. SUMMARY AND OUTLOOK
We have performed a comprehensive ab initio study of the influence of chemical composition and network defects on the density and bulk modulus of idealized, dense organosilicate glasses. Our results show that the transition between different mechanisms of deformation between Si-O-Si and Si-C-Si bonded materials leads to interesting nonlinearities of the mechanical properties with respect to C/O fraction in this class 054204-6 of materials. We also find that in O-rich OSG evenly distributed defects allow the greatest reduction of the total density while retaining mechanical stiffness. The properties of C-rich OSG are much more sensitive to the defect types we examined. All these findings can be explained by the different behavior of stiff Si-(CH 2 )-Si and flexible Si-O-Si bridges. Our results also suggest that the model sizes accessible to ab initio methods are not sufficient to examine the behavior of the elastic modulus for pore sizes larger than point defects. Our findings of internal stress in perfect network silicon hydrocarbide indicate that under realistic growth conditions a certain concentration of nonbridging sites will arise naturally to reduce the buildup of this strain.
The motivation to reduce the density of OSGs is to decrease their k value, which can be achieved by reducing the total electron density and the density of flexible dipoles. Other quantities, for example the strength of the flexible dipoles, which is connected to the bond polarity, also influence the dielectric constant. Therefore it would be interesting to calculate the dielectric constants directly, depending on composition and topology. The multitude of individual models generated in the course of this study can serve as the basis for such calculations.
Our results show the limitations of what can be achieved by means of even very effective ab initio methods. Larger models that allow lower pore concentrations as well as larger pore sizes need to be examined. In this direction, initial molecular dynamics simulations using the COMPASS 23 force field have been performed. 24 The results of those simulations and the results presented here largely agree. The disagreement on the behavior of -(CH 2 )-requires further careful validation against simulations at a higher level of theory, which will be the subject of future work.