title: Ab initio investigation of the impact of structural motifs on the fracture toughness of silicon
title: Ab initio Untersuchung zur Auswirkung von struktuellen Motiven auf die Bruchfestigkeit von Silizium
Creators: name: Müller, Tobias
Publisher: name: Friedrich-Alexander-Universität Erlangen-Nürnberg (FAU)
Year: 2025

description: In the 1920s, Griffith introduced his continuum approach for fracture. Since then, great efforts have been made to describe crack propagation more precisely within continuum mechanics. Nevertheless, it is evident that the atomic structure of the crack tip itself plays a crucial role in crack advancement, thus making it necessary to extend the concept of Griffith to the discrete atomic level. This enables the possibility to investigate essential and material specific processes, such as local bond rearrangements and path-dependent activation barriers. In general, crack opening is a complex chemical process on a multi-dimensional potential energy surface with many local minima and saddle points. The complexity increases even further when parameters such as temperature, pressure, or chemical environments are included. In this study, we discuss first density-functional theory geometry optimizations for crack tips in silicon to investigate brittle fracture mechanics. Our results demonstrate that cracks can propagate through a multitude of local energy minimum configurations connected by a variety of energy barriers. We highlight the complex nature of crack advancement and the challenge of modeling and identifying the correct fracture pathway, even for a material with a simple crystal structure like silicon. One advantage of using silicon is that it is well-described in literature. It is purely brittle, and extensive data on its fracture properties from experiments and simulations is available. The lattice parameters, bulk elastic constants, and surface energies of the relaxed surfaces were determined using calculations done in the generalized-gradient approximation with the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional, using a plane wave basis set, Vanderbilt ultrasoft pseudopotentials, and the semi-empirical Grimme D3 dispersion correction for including van der Waals forces. This data is needed to determine the analytical solution of the long-range atomic displacement field around the crack tip, which is used to displace the atomic lattice accordingly. The finite cutouts from the displaced atomic lattice are referred to as "pacman" structures. The boundaries are saturated with hydrogen atoms, outer atoms are fixed, and the inner part, including the crack surfaces, are relaxed. Thus, the displacement field serves as a boundary condition for the ab initio calculations. In order to increase the external load, two different approaches are tested. The first approach involves taking cutouts of the pacman structures at different loads, or rather at different stress intensity factors K, and relaxing them individually. The second approach starts at a certain relaxed pacman structure and re-scales the atomic positions linearly in consecutive steps with increasing or decreasing load, with each re-scaling step involving a geometry optimization. As fracture criterion, the length of the first bond closest to the crack tip is chosen. Both approaches lead to different critical loads for fracture. The re-scaling approach, in particular, shows significantly higher critical K-values, since structures are often trapped in stable local energy minima with disadvantageous dimer transition states. The initial calculations indicate a complex potential energy surface with several possible transition minima for crack propagation. Thus, extensive attention is given to the different possible crack tip relaxations and resulting structures of the crack tip. Furthermore, a strong relation between the surface relaxation and the formed crack tip structure is observed. After probing for different local minima and crack shapes, the focus was shifted to determining possible activation barriers for crack propagation. According to statistical mechanics and for subcritical loads, the reaction path should follow the lowest energy minima. The change in energy is defined relative to the non-relaxed pacman structure and is referred to as center energy. Unfortunately, the center energies of the different crack tip structures are not size independent and exhibit deviating relative energies and system size behavior. This size inconsistency is caused by the mismatch between the originally chosen "mathematical" center of the displacement field and the true "physical" crack tip center after structure optimization. In other words, the structure optimization was carried out with incorrect boundary conditions that correspond to the fixed outer atoms of the cutout, as they were derived from a pacman structure with a wrong displacement field center. Thus, the boundary conditions have to be adjusted until the true physical displacement field aligns with the analytical one for the far-field region where linear elasticity theory applies. This is illustrated with a thorough analysis of the change in displacements caused by the relaxation of the initial pacman structure. The results of the study show that by shifting the initial displacement field position of the analytical solution within the finite cutout, a significant improvement for the transition to the far-field solution can be observed. The mean mismatch of the displacement field can be plotted against the displacement field center position, which shows a parabolic relation with the minimum at the position where the erroneous contributions from the boundary become minimal. To correct the outer boundary, an iterative correction scheme was implemented. This scheme performs a least square fit between the displacement fields of the relaxed structure and non-relaxed pacman structures with arbitrary displacement field centers. With the new boundary conditions from the selected pacman, the next relaxation calculation is performed. This process is repeated until no better match is found. This flexible boundary regime has the advantage that the position of the displacement field within the finite cutout can be used as a direct criterion for crack propagation. This suggests that both the boundary conditions and the geometry of the crack tip play a crucial role in determining the activation barrier for crack propagation. By taking into account the effects of boundary conditions and crack tip geometry, it is possible to predict the critical loads for fracture and understand the underlying mechanisms of crack propagation more accurately. Overall, the flexible boundary regime provides a powerful tool for studying fracture in complex materials and structures, and has the potential to improve our understanding of the fundamental processes underlying fracture in brittle materials.
description: In den 1920er Jahren stellte Griffith seinen Kontinuumsansatz für Brüche vor. Seitdem wurden große Anstrengungen unternommen, um die Rissausbreitung innerhalb der Kontinuumsmechanik genauer zu beschreiben. Es ist jedoch offensichtlich, dass die atomare Struktur der Rissspitze selbst eine entscheidende Rolle bei der Rissausbreitung spielt. Daher ist es notwendig, das Konzept von Griffith auf die diskrete atomare Ebene zu erweitern. Dadurch besteht die Möglichkeit, wesentliche materialspezifische Prozesse wie lokale Bindungsumlagerungen und pfadabhängige Aktivierungsbarrieren zu untersuchen. Im Allgemeinen ist die Rissöffnung ein komplexer chemischer Prozess auf einer mehrdimensionalen potentiellen Energieoberfläche mit vielen lokalen Minima und Sattelpunkten. Die Komplexität nimmt noch zu, wenn Parameter wie Temperatur, Druck oder chemische Umgebungen einbezogen werden. In diesem Zusammenhang diskutieren wir erste Dichtefunktional-basierte Geometrieoptimierungen für Rissspitzen in Silizium um die Sprödbruchmechanik zu untersuchen. Wir zeigen in unserer Arbeit, dass sich Risse über eine Vielzahl lokaler Energieminima-Konfigurationen ausbreiten können, die über verschiedene Energiebarrieren miteinander verbunden sind. Dabei gehen wir insbesondere auf die komplexe Natur der Rissausbreitung und die Schwierigkeit bei der Modellierung ein, sowie auf das Auffinden des korrekten Bruchweges, selbst für ein Material mit einer so einfachen Kristallstruktur wie Silizium. Silizium eignet sich besonders für unsere Untersuchungen, da es in der Literatur gut beschrieben ist. Es verhält sich rein spröde und es sind umfangreiche Daten über seine Brucheigenschaften aus Experimenten und Simulationen verfügbar. Wir haben die Gitterparameter, elastischen Konstanten und Oberflächenenergien der relaxierten Oberflächen bestimmt. Diese Berechnungen wurden in Generalized-Gradient-Annäherung mit dem Perdew-Burke-Ernzerhof (PBE)-Austauschkorrelationsfunktional durchgeführt unter Verwendung eines Basissatzes aus ebenen Wellen, ultrasoften Vanderbilt-Pseudopotentialen und der semi-empirischen Grimme-D3-Dispersionskorrektur zur Berücksichtigung von van der Waals-Kräften. Diese Daten werden benötigt, um die analytische Lösung des langreichweitigen atomaren Verschiebungsfeldes um die Rissspitze zu bestimmen, die das Atomgitter entsprechend verschiebt. Die endlichen Ausschnitte aus dem verschobenen Atomgitter werden als "Pacman"-Strukturen bezeichnet. Die äußeren Si-Atome sind mit Wasserstoffatomen gesättigt und fixiert und der innere Teil, einschließlich der Rissflächen, wird relaxiert. Das Verschiebungsfeld wird als Randbedingung für die ab-initio-Rechnungen verwendet. Es wurden zwei verschiedene Ansätze getestet um die externe Zugkraft zu erhöhen. Der erste Ansatz verwendet Ausschnitte der Pacman-Strukturen bei verschiedenen Zuglasten bzw. verschiedenen Spannungsintensitätsfaktoren K und lässt diese einzeln relaxieren. Der zweite Ansatz geht von einer bestimmten relaxierten Pacman-Struktur aus und skaliert linear die Atompositionen in aufeinanderfolgenden Schritten in denen die Zuglast erhöht oder verkleinert wird. Jeder Skalierungsschritt beinhaltet eine Geometrieoptimierung. Als Bruchkriterium wird die Länge der ersten Bindung gewählt, die der Rissspitze am nächsten liegt. Beide Ansätze führen zu unterschiedlichen kritischen Lasten für den Bruch. Insbesondere der Reskalierungsansatz zeigt deutlich höhere kritische K-Werte, da Strukturen oft in stabilen lokalen Energieminima mit ungünstigen Dimerübergangszuständen gefangen sind. Erste Berechnungen zeigen eine komplexe potentielle Energieoberfläche mit vielen möglichen Übergangsminima für die Rissausbreitung. Dabei wird den verschiedenen möglichen Rissspitzenrelaxationen und den daraus resultierenden Strukturen der Rissspitze besonders große Aufmerksamkeit geschenkt. Es wurde beobachtet, dass die Oberflächenrelaxierung einen starken Einfluss auf die Struktur der gebildeten Rissspitze hat. Nach der Untersuchung verschiedener lokaler Minima und Risskonfigurationen wurde der Fokus darauf gelegt, mögliche Aktivierungsbarrieren für die Rissausbreitung zu bestimmen. Gemäß der statistischen Mechanik sollte bei unterkritischen Belastungen der Reaktionsweg entlang der niedrigsten Energieminima erfolgen. Die Energieänderung ist dabei relativ zur nicht relaxierten Pacman-Struktur definiert und als Zentrumsenergie bezeichnet. Allerdings sind die Zentrumsenergien der verschiedenen Rissspitzenstrukturen nicht unabhängig von der Größe des gewälten Pacman-Ausschnitts und zeigen unterschiedliche relative Energien und Systemverhalten. Dies liegt daran, dass das ursprünglich gewählte "mathematische" Zentrum des Verschiebungsfeldes von dem tatsächlichen "physikalischen" Rissspitzenzentrum nach der Geometrieoptimierung abweicht. Das bedeutet, dass die Geometrieoptimierung mit falschen Randbedingungen durchgeführt wurde, die den festgehaltenen äußeren Atomen während der Rechnung entsprechen, da sie von einer Pacman-Struktur mit einem abweichenden Verschiebungsfeld abgeleitet wurden. Daher müssen die Randbedingungen angepasst werden, bis das tatsächliche physikalische Verschiebungsfeld mit dem analytischen des Fernfeldes, in welchen lineare Elastizität vorausgesetzt werden kann, übereinstimmt. Dies wird veranschaulicht, indem eine gründlichen Analyse der durch die Relaxation verursachten Änderungen an den anfänglichen Pacman-Strukturen und den damit einhergehenden Änderungen des ursprünglichen Verschiebungsfeldes durchgeführt wird. Eine deutliche Verbesserung des Übergangs zur Fernfeldlösung wird durch eine Anpassung des Ursprunges des anfänglichen Verschiebungsfeldes der analytischen Lösung innerhalb des endlichen Ausschnittes beobachtet. Die akkumulative Änderung der Verschiebungen kann gegen die Position des Ursprungs des Verschiebungsfeldes aufgetragen werden. Es zeigt sich eine parabolische Beziehung mit einem Minimum an der Stelle, an der die fehlerhaften Beiträge minimal werden. Es wurde ein iteratives Korrekturschema implementiert, das die äußeren Randbedingungen korrigiert, indem die minimalen Fehlerquadrate zwischen den Verschiebungsfeldern der relaxierten Struktur und den nicht relaxierten Pacman-Strukturen mit beliebigen Verschiebungsfeldzentren bestimmt werden. Mit den neuen Randbedingungen der ausgewählten Pacman-Struktur wird die nächste Relaxationsberechnung durchgeführt. Dies wird so lange wiederholt, bis keine bessere Übereinstimmung mehr gefunden wird. Diese flexible Randbedingung hat außerdem den Vorteil, dass die Position des Verschiebungsfeldes innerhalb des endlichen Ausschnitts als direktes Kriterium für die Rissausbreitung verwendet werden kann. Dies legt nahe, dass sowohl der Randbereich als auch die Geometrie der Rissspitze eine entscheidende Rolle bei der Bestimmung der Aktivierungsbarriere für Rissausbreitung spielt. Durch Berücksichtigung der Auswirkungen von Randbereich und Rissspitzengeometrie ist es möglich, die kritischen Werte für den Bruch zu bestimmen und die zugrunde liegenden Mechanismen der Rissausbreitung genauer zu verstehen. Insgesamt bietet die flexible Anpassung des Randbereiches die Möglichkeit, Brüche in komplexen spröden Materialien und Strukturen zu untersuchen und unser Verständnis der grundlegenden Prozesse zu verbessern.