Quantifying the High-Temperature Separation Behavior of Lamellar Interfaces in γ-Titanium Aluminide Under Tensile Loading by Molecular Dynamics

γ-titanium aluminide (TiAl) alloys with fully lamellar microstructure possess excellent properties for high-temperature applications. Such fully lamellar microstructure has interfaces at different length scales. The separation behavior of the lamellae at these interfaces is crucial for the mechanical properties of the whole material. Unfortunately, quantifying it by experiments is difficult. Therefore, we use molecular dynamics (MD) simulations to this end. Specifically, we study the high-temperature separation behavior under tensile loading of the four different kinds of lamellar interfaces appearing in TiAl, namely, the γ / α 2 , γ / γ PT , γ / γ TT , and γ / γ RB interfaces. In our simulations, we use two different atomistic interface models, a defect-free (Type-1) model and a model with preexisting voids (Type-2). Clearly, the latter is more physical but studying the former also helps to understand the role of defects. Our simulation results show that among the four interfaces studied, the γ / α 2 interface possesses the highest yield strength, followed by the γ / γ PT , γ / γ TT , and γ / γ RB interfaces. For Type-1 models, our simulations reveal failure at the interface for all γ/γ interfaces but not for the γ / α 2 interface. By contrast, for Type-2 models, we observe for all the four interfaces failure at the interface. Our atomistic simulations provide important data to define the parameters of traction–separation laws and cohesive zone models, which can be used in the framework of continuum mechanical modeling of TiAl. Temperature-dependent model parameters were identified, and the complete traction–separation behavior was established, in which interface elasticity, interface plasticity, and interface damage could be distinguished. By carefully eliminating the contribution of bulk deformation from the interface behavior, we were able to quantify the contribution of interface plasticity and interface damage, which can also be related to the dislocation evolution and void nucleation in the atomistic simulations.


INTRODUCTION
Rising concern for environmental protection demands the quest for lightweight high-temperature materials to improve fuel efficiency in civil aviation. In recent times, titanium aluminide (TiAl) alloys (Dimiduk, 1999) are evolving as a potential material for rotating components operating at elevated temperatures, especially turbine blades. Compared to other high-temperature materials, alloys based on γ-TiAl (Appel et al., 2000) possess a high specific modulus and melting point, a low density, and excellent corrosion resistance. In particular, the specific strength of γ-TiAl is higher than that of superalloys at elevated temperatures, making it an attractive engineering material (Clemens and Mayer, 2013). γ-TiAl can exhibit a duplex, equiaxed, nearly lamellar (NL), or fully lamellar (FL) microstructure. Among them, the FL microstructure exhibits superior creep resistance and fracture toughness (Bewlay et al., 2016). It is characterized by thin lamellae, each of which belongs to one of the two constituent intermetallic phases, the γ phase (ordered face-centered tetragonal TiAl) and the α 2 phase (ordered hexagonal close-packed Ti 3 Al). In the FL microstructure, a lamella generally consists of multiple domains within which an ordered crystallographic structure is preserved. On the other hand, multiple lamellae form the socalled (grain-shaped) colonies, leading to interfaces at different length scales, namely, colony level, domain level, and lamellar level. Two types of interfaces in FL TiAl substantially affect its high-temperature deformation properties: the colony boundaries on the larger scale and the lamellar interfaces on the lowest scale, which are either of the type (c/c or c/α 2 ).
Continuum mechanical modeling of the micromechanics of materials has attracted considerable attention over the last decades. Appel et al. (2016) recently reviewed the different modeling approaches to TiAl proposed so far. Crystal plasticity models can help to account for specific deformation mechanisms (slip and twinning) and for describing the anisotropic material behavior in TiAl. Recently, such models (Schnabel and Bargmann, 2017;Ji et al., 2018;Schnabel and Scheider, 2020) were used to predict the colony boundary strengthening coefficient as a function of lamella thickness and to map the deformation in lamellar TiAl. Microstructureinformed multi-scale models bear promise to accelerate the alloy development in γ-TiAl. To increase their accuracy, they should also include the role of interfaces and alloying elements. Obtaining the parameters required for continuum mechanical interface models from experiments is difficult so that molecular dynamics (MD) simulations appear to be the most promising tool to this end.
In MD simulations of TiAl, the different intermetallic phases (c and α 2 ) are modeled with their respective crystallographic systems (i.e., γ: face-centered tetragonal and α 2 : hexagonal closepacked) and appropriate orientations as well as with suitable interatomic interaction potentials (Zope and Mishin, 2003;Kim et al., 2016). In recent time, several MD studies have been devoted to deepen our understanding of the deformation behavior of γ-TiAl (Zhou et al., 2004;Xie et al., 2015;Kanani et al., 2016;Wu et al., 2016;Li et al., 2017;Feng et al., 2018;Hui et al., 2018;Cao et al., 2019;Ding et al., 2019;Feng et al., 2019;Li et al., 2019;Li et al., 2020) in combination with experiments. Kanani et al. (2016) performed MD shear simulations of a distinct c/c interface, observing different deformation mechanisms and strong in-plane anisotropy of shear strength. Wu et al. (2016) studied the influence of surface defects on the mechanical behavior of the γ phase under tensile deformation. They (Wu et al., 2016) observed for the defect-free bulk sample brittle fracture, but plastic deformation before fracture for samples with surface defects normal to the loading direction. In an extensive study using first principles, Ji et al. (2018) obtained a deformation map for plastic deformation modes using the orientation-dependent effective energy barrier. Li et al. (2019) studied the homogeneous tensile deformation behavior of the c/α 2 interface at 0.01 K and found an influence of boundary condition on the observed failure mode with plastic deformation for free boundaries and fracture for constrained boundaries. These observations can also be interpreted by the Schmid factor analysis of the corresponding slip system. A recent MD study (Ding et al., 2019) has been devoted to the effect of grain size and temperature, however, for polycrystalline TiAl. A decreasing Young's modulus and average flow stress was observed with increasing temperature. Together, these studies can help to understand experimental findings about room temperature ductility. However, more atomistic insights are still necessary to understand the thermomechanical deformation behavior of a single lamellar interface in γ-TiAl at elevated service temperatures. Our understanding of the combined microscopic effects and mechanisms occurring in plasticity and fracture of TiAl still remains limited (Appel et al., 2016), even though some small-scale testing has been conducted revealing the failure mechanisms in the microstructure (Werwer et al., 2007;Barbi et al., 2012).
The finite element method is frequently used in mesoscale and macroscale models of TiAl. Such models should include the deformation behavior of bulk constituent phases as well as different interfaces (c/α 2 and c/c). In γ-TiAl, the anisotropy arising from the complex microstructure can be captured by crystal plasticity models (Schnabel and Bargmann, 2017;Schnabel and Scheider, 2020). For modeling the interfaces between colonies, domains, and lamellae, cohesive zone models are promising tools. These phenomenological models initially developed to handle discontinuities due to crack propagation and fracture (Hillerborg et al., 1976;Li and Chandra, 2003;Scheider, 2018) can be embedded in finite element models by the so-called interface elements to also describe the material separation behavior at both macroscale and micro-scale interfaces (Scheider, 2009a;Simonovski and Cizelj, 2015;Scheider, 2018). For lamellar γ-TiAl, micromechanical lamellae separation has been studied in Werwer and Cornec (2000), Werwer et al. (2007), and Wei et al. (2009). Cohesive zone models rely on a traction-separation law (TS law). On the macroscale, its parameters can be identified by appropriate mechanical tests (Cornec et al., 2003). However, debonding processes are difficult to observe with in-situ experiments on the nanoscale of atomistic lattice boundaries. Classical MD simulations are typically more suitable on that scale. The identification of TS laws from atomistic simulations instead of experimental tests has started in the early 2000s and has remained an active field of research since then (Yamakov et al., 2006;Zhou et al., 2008;Zhou et al., 2009;Krull and Yuan, 2011;Möller et al., 2018). For high-temperature alloys, such as γ-TiAl, in particular, temperature dependency has to be included in such studies. An important aspect when the TS law is derived from micromechanical or atomistic simulations is that usually the damage in a finite volume (called the damage zone) is taken as a basis for the separation behavior. The elimination of the "bulk" response from the damage zone (i.e., solely considering the interface damage) has been addressed here not only based on a method already used by the authors (see (Scheider, 2009b, Scheider, 2009a) but also in a similar fashion employed by others (see, e.g., (Paggi and Wriggers, 2011)).
In this article, we perform MD simulations of four single lamellar interfaces (one c/α 2 and three c/c in different orientations) in γ-TiAl under tensile load from room temperature to elevated service temperatures. This way, we first obtain general atomistic insights into the separation behavior of single lamellar interfaces. Second, we obtain atomistically informed traction-separation (TS) curves for each material interface, depending on temperature, strain rate, and preexisting defects. These are used to formulate an atomistically informed cohesive zone model, which can be used in future finite element mesoscale simulations of TiAl.

Molecular Dynamics Model
At the micro-or nanoscale, materials can be modeled as manyparticle systems. MD is a numerical technique to study the evolution of such systems in time by computing the trajectories of the individual particles. This is accomplished by numerical integration of Newton's equation of motion. If the particles represent atoms, MD simulations provide an atomic resolution (Cai et al., 2012). Herein, all MD simulations were performed using the open source package LAMMPS large-scale atomic/molecular massively parallel simulator (Plimpton, 1995) distributed by the Sandia National Laboratories. For metallic materials, interatomic interactions in MD are often modeled in a (semi)empirical manner using the so-called embedded-atom method (EAM) of Daw and Baskes (1984). Following this approach, we define the total potential energy (U Total ) of the many-particle (N atoms) system as Here, r ij and ϕ ij are the distance and interaction potential between the ith and jth atom. N indicates the total count of atoms in the system. F i denotes the embedding energy of the ith atom due to the cumulated host electron density j ≠ i ρ j (r ij ) provided by the surrounding atoms. The pair interaction potential and electron density are summed over the neighboring j atoms defined by the cutoff. For γ-TiAl, we used the EAM potential from Zope and Mishin (2003), which was validated by a comparison of experiments and ab initio simulation data and found to yield, for example, accurate lattice constants, cohesive energy, and elastic constants. In particular, its ability to predict the lattice thermal expansion makes this potential suitable for MD simulations at high temperatures. Therefore, it has also been used by Xie et al. (2015), Kiselev and Zhirov (2014) for studying the synthesis of intermetallic γ-TiAl from Al and Ti plates at 1400 K and the solidification of α 2 -Ti 3 Al from the liquid phase at 2000 K. Remarkably, this potential correctly predicts the α 2 -Ti 3 Al equilibrium structure although the relevant information was not included in the fitting database. For these reasons, we chose this EAM potential of Zope and Mishin (2003), whose developers underline its suitability for large-scale atomistic simulations of plastic deformation and fracture in γ-TiAl even at high temperatures.

Cohesive Zone Model
In a cohesive zone model, interfaces are idealized as two surfaces initially attached to each other without any material volume in between. Tractions on the interface are defined by a vector field t and may cause a relative displacement of the two surfaces forming the interface, which can be modeled by the so-called separation vector EuF u + − u − with u + , u − being the displacements of the two surfaces forming the interface, respectively. The relation between traction vector and separation vector can be formulated by the so-called TS law t t(EuF), which can be understood as the mechanical constitutive laws of the interface itself. A host of different ways to define such laws has been proposed (see, e.g., the reviews in Park andPaulino (2011), Scheider (2018)). Here, a recently developed model is used, which splits the separation additively into a reversible (elastic) and an irreversible (plastic) part EuF EuF el + EuF pl .
(2) Similar to bulk plasticity models, the interface traction depends only on the elastic part of the separation EuF el . A large variety of TS laws is in principle possible. Herein, we focus on the mechanical behavior of the interlamellar interfaces in TiAl normal to the interfaces under normal loading only. Therefore, we only consider the normal components of the traction and the separation vector. The normal components of the full, elastic, and plastic separation vector are denoted by Eu n F, Eu el n F, and Eu pl n F. Analogously, the normal component of the traction vector is denoted by t n . The simplest relation of separation and traction is a linear one, which is also very common for cohesive zone models since the elastic separation is usually small. Since the traction is the derivative of the elastic potential (per unit area) W, we start with the following second order polynomial Here, W 0 describes the elastic energy in the absence of any damage, and the damage parameter d ranges from d 0 for the damage-free state to d 1 at failure. The material parameter C n defines the interface stiffness. This relation can easily be generalized to three dimensions. The traction vector follows from the above elastic energy as t n zW zEu el n F (1 − d)C n Eu el n F.
Following the classical setting of damage and plasticity theory (Maugin, 1992), we define the so-called fictitious undamaged stresst n t n /(1 − d). The plastic regime of our traction separation law is assumed to be defined by a yield function Here, t 0 denotes an initial yield stress and D(d) is reduction due to damage. The results of our MD simulations indicated that plastic hardening of the interface is negligible, and it is thus not represented in Eq. 5. This also implies that the Helmholtz free energy (per unit area) of the interface is W from Eq. 3, which is thus used in Eq. 9 to compute the thermodynamic driving force of the single internal variable d of our model. The flow rule associated to the yield function g provides the plastic flow rate normal to the interface with _ λ being the plastic multiplier. The accumulated plastic separation α can be defined as The damage-related reduction of yield stress D in Eq. 5 is defined herein heuristically as with two material parameters p and S. The evolution of the internal variable d can be defined in a thermodynamically consistent manner (see, e.g., Maugin, 1992) via its thermodynamic driving force which by standard thermodynamic arguments from classical plasticity theory directly leads to the evolution equation Here, the initial value of d is zero, and α nucl represents an additional model parameter below which no damage at all is assumed to occur. Moreover, it can be shown that where α fail denotes the plastic separation at failure, that is, at d 1. Note that with Eq. 11, α fail can be used to replace effectively the model parameter S. Thus, our cohesive zone model can be written using only the following five material parameters: C n for characterizing the elastic properties of the interface and t 0 , p, α fail , and α nucl for characterizing plasticity and damage.
While the cohesive zone model has been implemented into the commercial finite element code ABAQUS ® , a python implementation of the above equations is sufficient for the calculations presented here. The python class for the cohesive zone model is available at Scheider (2020).

Molecular Dynamics Simulation-Preprocessing
In γ-TiAl processing, the target microstructure is often achieved through the controlled cooling of the high-temperature α phase, which results in the precipitates of γ-plates through the reactions α → α/c → α 2 /c and α → α 2 → α 2 /c (Kim and Dimiduk, 1991). Such nucleation of the γ phase from the parent α or α 2 phase leads to four different lamellar interfaces. In addition to the c/α 2 interface, there exist three different γ/γ interfaces in the FL microstructure (Jin and Gray, 1997;Appel et al., 2011;Kanani et al., 2016), which represent different rotational variants. Herein, we studied all the four different interfaces. To this end, we constructed atomistic simulation models (refer Figure 1) whose two halves were a combination of intermetallic constituent phases representing two different lamallae. The interface between both phases represented the respective interlamellar interface. Based on the general setting illustrated in Figure 1, we designed the simulation domains in the following way. The face-centered tetragonal crystallographic structure with a 0 3.998 Å and c/a 0 1.047 (Zope and Mishin, 2003) of the γ-TiAl phase was represented by 96 atomic layers along the X-direction [1-1 0] and Y-direction [1 1-2], respectively, and by 48 atomic layers in the Z-direction [1 1 1]. Likewise, the hexagonal close-packed crystallographic structure with a 0 5.7884 Å and c/a 0 0.820 of the α 2 -Ti 3 Al phase was represented by 96 atomic layers in the X-direction [1-1 0], 64 atomic layers in the Y-direction [1 1 0], and 48 atomic layers in the Z-direction [0 0 1]. Here, the respective number of atomic layers was chosen to ensure periodicity in all directions. The c/α 2 interface was aligned with the close-packed plane and direction, that is, {111} c (0001)α 2 (Appel et al., 2011). The three rotational variants of the γ phase were achieved by a rotation around 〈111〉 c by 60°, 120°, and 180°, respectively, resulting in a pseudo twin (PT), rotational boundary (RB), and true twin (TT) configuration, denoted as c/c PT , c/c RB , and c/c TT interfaces, respectively.
In our study, emphasis was placed on understanding the separation behavior of the single interlamellar interface subject to tensile loading. Accordingly, we consider a defect-free interface model designated as Type-1 (refer Figure 1) to assess the theoretical strength and deformation behavior. Further, we considered the single interface models with a defect, which we refer to as the Type-2 model (refer Figure 1), to study the influence of the preexisting defect at the interface. Here, a line defect, that is, dislocation, was avoided on purpose as we were interested in neither the dislocation dynamics nor the interaction Frontiers in Materials | www.frontiersin.org January 2021 |Volume 7 | Article 602567 of existing dislocation with chemical species. Moreover, the simplified model of single dislocation poses several limitations: stability of the dislocation at elevated temperature, the bias introduced by the initial dislocation type, and unrealistic starting dislocation density. To circumvent these limitations, we chose a void type defect instead. Further, the preexisting void can be associated with the crack at the interface, which is often considered in continuum fracture studies to understand the crack propagation under different loading conditions. From the atomistic standpoint, such defect may act as dislocation sources during deformation, overcoming any modeling bias, as we discussed before. Accordingly, the Type-2 model includes a penny-shaped void at the interface with a radius of 1.6 nm and a height of 1 nm (refer Figure 1). The center of this void was chosen to coincide with the center of the simulation domain so that the void covers two atomic layers in each of the two intermetallic phases in the simulation model. The influence of the initial void shape and size on the interface separation behavior falls beyond this study's scope. In test simulations, we found that a simulation domain of the size 13.7 nm × 8 nm × 22 nm was a good trade-off between computational cost and the requirement to reduce boundary effects. Using such domains, Type-1 models  consist of 147,456 atoms in total and Type-2 models of 147,002 (after removing atoms belonging to the void region).
In Figure 2, the first row shows the constructed atomistic models of the four interlamellar interfaces in TiAl after an initial energy minimization. The second row displays in this initial configuration for each atom the so-called centrosymmetry parameter (CSP), which measures the local lattice disorder around the atom, for example, due to defects, surfaces, or interfaces (Stukowski, 2009). For an atom whose direct neighborhood forms a perfect lattice, the CSP is zero. By contrast, the CSP takes on high values if the atom constitutes the core of a dislocation or is located next to an interface or surface. The CSP was computed using OVITO (Stukowski, 2009). Here, blue atoms represent the γ phase, green atoms represent either an interface or the α 2 phase, and red atoms represent the boundary of the simulation domain. Note that at the c/c RB interface, the original stacking order of the γ-phase is preserved and thus no elevated CSP was observed.
The preprocessing of the simulation samples was performed in several steps. In the first step, Type-1 (defect-free) interface models were energetically minimized using a conjugate gradient algorithm in LAMMPS large-scale atomic/molecular massively parallel simulator (Plimpton, 1995). For Type-2 models, now the atoms in the void regions were removed and energetically minimized. For both Type-1 and Type-2 models, the second step was assigning velocities to all the atoms using a Gaussian distribution to initialize the start temperature. Subsequently, the samples were heated up to the target temperature and equilibrated at constant temperature and pressure for 40 ps in case of Type-1 models and 100 ps in case of Type-2 models. Finally, in the third step, the samples were equilibrated at constant volume and target temperature for 100 ps so that the kinetic energy is well distributed within the model. The Nosé-Hoover thermostat was used for all the thermal equilibrations. A time step size of 1fs and periodic boundary conditions in all directions were applied in the preprocessing simulations.

Molecular Dynamics Simulation-Deformation
After the preprocessing as described in section 2.3, the deformation according to the applied strain rate was simulated. Three different regions were distinguished for imposing boundary conditions: a loading region, a fixed region, and a mobile region in between. These regions are depicted in Figure 3. Here, we choose a boundary-driven deformation technique to perform the tension test. Accordingly, we selected a few layers of atoms (spanning 2 nm) far away from the interface at the boundary of the simulation box in the Z-direction. The loading region's thickness is carefully chosen to ensure the smooth transfer of displacement rate into the system without any boundary artifact. In the loading region, we prescribe similar to Zhou et al. (2008), Gupta et al. (2016)  1 · 10 9 s −1 (both cases were studied in separate simulations) in the Z-direction was applied to all atoms in the loading region. It is a well-known fact that MD simulations possess severe time scale limitation and inherently suffer high strain rate in comparison to the experiment. Following an investigation conducted by Wu et al. (2016), who showed that the effect of the strain rate is reduced for rates smaller than _ ε glo ZZ 10 9 s −1 in γ-TiAl, we chose the above mentioned strain rates in this study. In the fixed region, the same procedure was used but the displacement rate was set to zero, such that atoms in this region remained stationary. The positions and velocities of the atoms in the loading region were updated based on the assigned constant displacement rates, whereas for the atoms in the mobile region, they were gained in the usual way by numerical integration in time with a time step size of 2 fs. A constant target temperature and periodic boundary conditions in X-and Y-directions were maintained during the deformation simulations. A nonperiodic boundary condition was applied in the Z-direction. Thus, the loading can be described as isothermal uniaxial straining. Figure 3 highlights a special zone within the mobile region, in the following denoted as the local observation region. It stretches over 3.5 nm in both directions from the interlamellar interface, that is, over in total h loc 7 nm. The role of this region is to observe the local separation behavior; hence it does not affect the simulation but supports the postprocessing of the results.
The local atomic stress tensor S i (units: stress times volume) of each atom i is calculated by at each discrete time step in the simulation according to Thompson et al. (2009). Here, ⊗ denotes a tensor product and N is the set of atoms in the model. The first term on the righthand side represents a kinetic, thermal contribution with m i being the mass of the ith atom and v i its velocity vector. The second term refers to the pair interactions between the ith atom and its surrounding atoms with r i and r j the position vectors of the ith and jth atoms, respectively, and F ij is the interaction force vector between the ith and the jth atom. For the transition to the continuum scale, the atomic stresses have to be divided by the current subdomain volume V box , giving where the index i runs over all atoms in the subdomain within the simulation model considered.
For mapping atomic stresses to the continuum level, we used two types of such averaged stresses, a global stress, S glo , and a local stress, S loc . For S glo , we averaged across the whole simulation domain in the deformed configuration (initially, V box 22 × 13.7 × 8 nm 3 ), which provides a quantity similar to the pressure tensor in the MD simulation model. For S loc , only the atomic stresses in the local observation region were averaged (initially, V box 7 × 13.7 × 8 nm 3 ). The ZZ component of this stress tensor, S loc ZZ , was assumed to correspond to the component of the interface traction vector normal to the interface (t n , as introduced in subsection 2.2) in the cohesive zone model (Dandekar and Shin, 2011). Therefore, whenever S loc ZZ is considered, the term (normal) traction will be used. S loc ZZ and S glo ZZ are expected to be equal since both correspond to the same total loading imposed by the boundary conditions at upper and lower boundaries. Nevertheless, we use these two different quantities in the following to match the variables that we use for discussing the deformation. For the deformation, it is favorable to distinguish between the one in the local observation region and the global one because the latter may, especially in the inelastic regime, also contain some bulk effects which we would like to eliminate from our discussion of interface mechanics.
For relating the atomistic to the continuum level, not only stresses but also strains and the interface separation have to be discussed. The global MD-based strain in the Z-direction ε glo ZZ is the change of the total simulation box size in the Z-direction divided by its initial length in the Z-direction. The MD-based (normal) interface separation in the local observation region is the difference between the average displacement in the Z-direction of the atoms in the layer above and below the local observation region (refer Figure 3), that is, Both the MD-based tractions and the interface separation were computed in our simulations every 100 fs.

Identification of Damage Specific Interface Separation
The separation Eu loc n F defined by the difference between u loc + Z and u loc − Z according to Eq. 14 depends on the size of the local observation region. The reason is that the deformation in this region is not only governed specifically by interface effects but also governed to some extent by deformation of the bulk material, in particular its elastic strain, which occurs independently of the interface and which should not be modeled as part of the deformation resulting from the interface. To eliminate such bulk effects, we follow the ideas used already in Scheider (2009b, a). That is, we first assume that the interface-specific effects such as damage, which we want to capture by our cohesive zone model, occur only in the neighborhood of the interface, that is, in the local observation region. Outside this region, we assume the deformation to be equal to the bulk deformation that would result from a global loading S glo also in the absence of the interface. This bulk deformation is thus assumed to correspond to a strain Now, to isolate the interface-related part of the deformation in the local observation region, we have to subtract from the overall deformation Eu loc n F of this region the part attributed to the general bulk deformation, giving With the interface-related separation Eu n F and an associated traction S loc ZZ , the parameters of the cohesive zone model from section 2.2 can be identified by fitting Eqs 4-10 to the results of the MD simulations.
Moreover, the total dislocation density in the MD simulations was calculated using the dislocation analysis tool provided in the visualization software OVITO (Stukowski, 2009;Stukowski and Albe, 2010) and used as a measure of plasticity on the continuum level. The starting point of micro-plasticity in a continuum mechanical sense can then be identified by a significant increase in dislocation density.

RESULTS AND DISCUSSION
In the following, we first present in section 3.1 the results of the MD simulations for different interlamellar interface models, considering the effect of temperature, strain rate, and preexisting defects. In subsequent section 3.2, the calculation of TS parameters for cohesive zone models from the results of the MD simulations is pointed out. Finally, in section 3.3, we discuss the results of this procedure in more detail.

Molecular Dynamics Simulation Results
In the following, the separation behavior of c/α 2 and c/c interfaces under normal loading in MD simulations is compared for the Type-1 and Type-2 interface models for Frontiers in Materials | www.frontiersin.org target temperatures T ∈ {300K, 500K, 700K, 900K} and for the two different strain rates _ ε glo ZZ 10 9 /s and _ ε glo ZZ 10 10 /s imposed by the velocities enforced on the atoms in the loading region of the simulation model (refer Figure 3). Figure 4 shows the global stress normal to the interface S glo ZZ for Type-1 models of the four single interlamellar interfaces according to Eq. 13 vs. the strain of the simulation box. One can observe a monotonous relation between stress and strain for all temperatures, interfaces, and strain rates up to a point, denoted as yield stress, which is characterized by significant nucleation of dislocations. In Figure 4, this point is marked by symbols (triangles for _ ε glo ZZ 10 10 /s and circles for _ ε glo ZZ 10 9 /s). Briefly after this point, we observe either immediately or after a short regime of plastic hardening a drop of stress. The yield stresses for all interface types and temperatures are summarized in Table 1. Generally, yield stress decreases with increasing temperature. At any given temperature, the c/α 2 interface was found to exhibit the highest yield strength, followed by the c/c PT and c/c TT interfaces and finally by the (weakest) c/c RB interface. The observed yield stresses ranged between 5.57 and 12.67 GPa. Of course, one has to keep in mind that these values are theoretical values for ideal, perfect interfaces.

Type-1 Interface
Some interface models (c/c PT and c/c TT ) show a sharp drop in stress after reaching the peak value at room temperature for low strain rate (refer Figure 4), suggesting a brittle-like fracture failure characterized by interface separation. Oscillations in the stress-strain curves are observed for all the interfaces after failure. Especially, for the higher strain rate (studied here),  the vibration effect is pronounced, visible by substantial oscillations of the stress-strain curve, which is understandable noting that these oscillations are dynamic effects increasing thereby the more the conditions are different from the quasistatic case. In general, one can observe that a reduction of strain rate typically reduces the yield stress for the cases studied herein (see Table 1). This observation reaffirms the reported strain rate sensitivity (Spearot et al., 2009), that is, reduction in strain rate reduces the flow/nucleation stresses, typically for materials with reduced defect activation volumes in the nanoscale regime, for example, nano-twinned copper (Zhu et al., 2007). Especially when such small-volume materials are deformed to the elastic limit, the dislocation nucleates result from strain localization. It is worth emphasizing that the reduced activation volumes impact the nucleation stresses, making them susceptible to strain rate and temperature sensitivity (Zhu et al., 2008). Interestingly, MD simulations of a pure γ-phase show such a strain rate sensitivity by the authors Wu et al. (2016), who also rationalized the observation to some of the above sources (Zhu et al., 2007(Zhu et al., , 2008Spearot et al., 2009). Similar to our observation in c/α 2 interfaces, the γ-phase is the softer phase carrying the most inelastic deformation (Wu et al., 2016). Thus, the above reported strain rate sensitivity also appears suitable to explain our results. For both applied strain rates, we observe an increase in hardening with a rise in temperature, signaling that the dislocation activities giving rise to hardening are thermally driven. Furthermore, the observed strength values are in good agreement with previous studies, especially for c/α 2 .  Type-2 Interface Figure 5 shows the global stress-strain curves for the Type-2 interface models, that is, with a preexisting void. The presence of the void type defect reduces the interface strength due to high local stress concentrations around the void and hence also the global load under which interatomic bonds break. Consequently, in comparison to the Type-1 models, the yield stress is reduced for all the interfaces and the high-stress regions near the void become the dominant dislocation sources. Such readily available dislocations lead to a more pronounced hardening after the yield stress has been reached compared to the Type-1 models. As in Type-1 models, the yield stress decreases also in Type-2 models with temperature (see Table 2). And again, it is highest for the c/α 2 interface, followed by c/c PT , c/c TT , and finally the (weakest) c/c RB interface at all temperatures studied. Generally, for low strain rate simulations, the void in Type-2 models reduces yield stress for all the four interfaces by around 35% for higher temperatures and more than 40% for 300 K. In Figure 5, all the Type-2 interface models show similar trends at the strain rate _ ε glo ZZ 10 10 /s. That is, the tensile stress increases further after the onset of dislocation nucleation (triangles in Figure 5) followed by a smooth drop indicating a rather ductile fracture process. This interpretation, where naturally substantial energy dissipation during the fracture process can be expected from its ductile nature, is also supported by the absence of major vibrations in the stress-strain curve after separation (in contrast to Type-1 models, refer Figure 4). For the strain rate _ ε glo ZZ 10 9 /s, after reaching the yield stress, typically a short dip in stress occurs followed by hardening. Whether the drop is more pronounced or the hardening (or both of them even level each other out) depends on the dislocation and void evolution. The void growth reduces the stress (softening), while the nucleation of dislocations leads to hardening. The Type-2 models illustrated in Figure 5 differ in two important aspects from the previously shown Type-1 models. First, the strain values where dislocation nucleation starts (yield strain, indicated by triangles and circles) are nearly independent of temperature for most of the interfaces studied. The reason may be the dominant role of the preexisting void for dislocation nucleation, which reduces the effect of temperature. Second, the hardening from dislocation-driven plasticity decreases with increasing temperature due to the pronounced dislocation mobility. More information on the dislocation evolution is shown in Figures 6 and will be discussed in the subsequent subsections. The observed high-temperature deformation behavior, along with the underlying dislocation-mediated mechanisms, is discussed in more detail in this article's Supplementary Material.

Traction-Separation Behavior of the Interfaces
As discussed above, TS curves, which form the core of cohesive zone models, can be extracted from MD simulations by fitting the parameters in the equations of the model to the Eu n F-S loc ZZ curves observed in MD simulations. Since the investigation of strain rate effects is (at least in this high strain rate regime) of minor practical relevance, this process is discussed only for the strain rate _ ε glo ZZ 10 9 s −1 in this section.

Type-1 Interface Models
First of all, it is important to note that not all configurations fail at the interface, namely, for the c/α 2 interface, fracture occurs in the interior of the γ-phase at most of the temperatures. In this article's Supplementary Material, the evolution of the simulation domain during the imposed deformation is illustrated (refer Supplementary Figures S1, S2), which clearly shows the location and the evolution of voids and complete failure for all configurations at T 900K. Hence, in case of Type-1 models, we focused on c/c interfaces only, which in general show interface fracture. 1 For the c/α 2 interface, it is in case of a Type-1 interface model neither possible nor of interest to identify a full TS curve. Rather it is reasonable to model this interface on the continuum level as an ideal connection between a γ-phase and an α 2 -phase and have an additional continuum damage model for the bulk γ-TiAl crystal, which is, however, outside the scope of this investigation. The atomistically computed local TS curves are shown in Figure 6 for the different c/c interfaces and temperatures. Together with the TS curves, the evolution of the dislocation density is also shown, indicating an onset of plastic deformation close to the peak traction. After this peak, the interface traction is decreasing, indicating the degradation of the load-bearing ability of the sample. The substantial differences between the TS curves considered in this section and the figures (refer Figures 4, 5) shown in the previous section is a consequence of the fact that in the previous section, we considered the overall load and deformation of the global simulation domain. In contrast, in this section, we focus on the local interface separation, according to Eq. 16 with effects of the bulk deformation already eliminated. Interestingly, eliminating these effects significantly reduces the size of the elastic regime. One may also notice that the separation is not monotonically increasing in the majority of TS curves. The reason for this behavior may be that the dislocation evolution is not homogeneous throughout the structure; in the first stage after the peak stress, dislocations may sometimes move more in the bulk region, sometimes more in the local region. The respective other region is then unloaded for a short time. However, this only happens in the first 10% of stress reduction and does not change the big picture substantially. The peak atomic tractions in the loading direction decrease with increasing temperature (Figure 6), similarly as in the global stress-strain curves discussed in the previous section. The regime directly after the stress drop depends on temperature as well, as discussed in the Supplementary Material for the global stress-strain behavior.
The shape of the TS law has been fitted to the MD results shown in Figure 6, assuming a nonhardening interface for simplicity due to the brittle fracture observed in Figure 6. Anyway, for the initially undamaged interface, damage starts almost together with plasticity. Eventually, the following parameters had to be identified: the elastic parameter C n from Eq. 4, the initial yield strength t 0 from Eq. 5, and the three damage parameters α nucl , α fail , and p from Eq. 11. Parameter fitting was automated using the optimization module of Scipy (2020). The results for the temperature T 900 K are shown in  The agreement between MD simulations and the fitted cohesive zone models depicted in Figure 8 is in each case very good. Remarkably, this holds not only for the TS curves themselves but also for the normalized dislocation density in the MD simulations and the corresponding damage-corrected plastic separation (1 − d)α. This indicates that our model is not purely phenomenological but appears to capture at least to some extent the physics. The fitted parameters for the cohesive zone model are given in Table 3 for a temperature of 900 K. The values for the other temperatures are listed in Supplementary Table S1 (refer Supplementary Material). It is worth mentioning that the parameters depend almost linearly on temperature enabling an easy transfer of our results to other temperatures in the same range. Figure 7 shows the TS curves for the Type-2 interface models with the preexisting void for different temperatures as observed in MD simulations. These curves generally exhibit three regimes: in the first one, traction increases up to a certain peak value. In the second one, the traction remains on a high level with some fluctuation (which is not unusual for MD simulations), followed by a third regime where it finally decreases. While the elastic behavior and the general temperature dependence are similar to the Type-1 models, also characteristic differences are observed.

Type-2 Interface Models
First of all, in contrast to Type-1 interface models, failure is always observed at the interface itself, also for the c/α 2 interface, so that this interface can also be addressed by a cohesive zone model. Second, the peak traction in Type-2 interface models is generally by around 35-50% lower than for the defect-free Type-1 models. Third, in Type-2 models, dislocation density increases more gradually after the onset of plasticity, and instead of the drop of the traction directly after the onset of yielding, it remains high for some separation. This indicates a significant regime of interface plasticity without damage evolution.
As for the Type-1 models, the parameters of a cohesive zone model were fitted to the results of MD simulations for the Type-2 models (see Figure 9 for T 900 K). Again a very good agreement between the results of the MD simulations and the fitted cohesive zone model is observed. In particular, the damagefree plasticity region with a constant traction corresponds extraordinarily well with the increasing branch of the dislocation density (black dash-dotted line), which is also reflected by the damage-corrected plastic separation (1 − d)α (blue dash-dotted curve). The parameter values fitted at T 900K are given in

Discussion
From the parameters identified for the different temperatures and interface types (Tables 3, 4; Supplementary Tables S1, S2 in Supplementary Material), it is possible to calculate the separation energies Γ 0 based on with Eu n F fail as the separation at which failure of the interface is observed and t n S loc ZZ . The resulting values of Γ 0 for all temperatures are given in Table 5. It turns out that the energies vary slightly for the different interfaces, and in general, they have a decreasing trend with temperature. Remarkably, the prexisiting void in the Type-2 interface models seems to decrease the interface strength by around 30% (see above), but for the interface separation energy, no such clear trend is found in our simulations.
A quantitative comparison of the cohesive zone model parameters with values from the literature is difficult due to the sparse data available. However, there are a few articles dealing with the strength of lamellar interfaces in γ-TiAl. A thorough combined experimental and numerical study was proposed by Werwer et al. (2007), where precracked PST (polysynthetically twinned) crystal specimens made of lamellar γ-TiAl were tested in different orientations at room temperature. Cohesive model parameters (i.e., cohesive strength and separation energy) for a TS law were determined by fitting fracture simulations to the experimental data. A very low cohesive strength of about 100 MPa was found, but the identified separation energies (Γ 0 3 . . . 6 J/m 2 ) are comparable with our findings (ref . Table 5A,5B).
An explanation why the strength observed in Werwer et al. (2007) is so much (roughly by a factor of 50) lower than the one determined herein (t 0 in Table 4) may be given on the basis of a recent study of Möller et al. (2018), who performed MD simulations of a crack tip in a brittle tungsten metal to calculate TS law parameters. They concluded that for a simulation box size comparable to the one considered in our study, the stresses are still within the K-dominated field. That is, the opening stress ahead of the crack tip is proportional to the square root of the distance, and a TS law with constant values may not be suitable to describe the interface strength properly. By contrast, they concluded that the separation energy may be largely independent of the size of the simulated domain already for a domain size comparable to the one considered herein. In another investigation by Krull and Yuan (2011) on a similar but more ductile material, aluminum, the authors did not consider the size effect and also found unrealistically high values for the cohesive strength of more than 10 GPa.
For the computational results presented in this article, this means that the separation energies of the interface are most likely suitable for utilization in a continuum level finite element simulation where the interface may be discretized by interface elements. By contrast, for the cohesive strength, one should probably use much lower values than reported herein if the length of the cohesive finite elements used on the continuum scale is larger than the size of the K-dominated zone.

SUMMARY AND CONCLUSION
To summarize, the high-temperature deformation behavior of the Type-1 (defect-free) and Type-2 (preexisting void) single lamellar interface models (one c/α 2 and three c/c) has been studied here under tensile load using MD simulations. The results of these simulations were used to fit the parameters of the traction separation laws of cohesive zone models, which can then be used for continuum level models, for example, using the finite element methods. From this study, we can draw a number of conclusions on the MD simulations alone and some further conclusions on the derivation of TS curves from the atomistic simulations. In our simulations, the c/α 2 interface was found to be the strongest, followed by the c/c PT , c/c TT , and c/c RB interfaces. Particularly, the c/α 2 interface for the ideal (defect-free) structure was found to be stronger than the adjacent bulk material, and failure never occurred at the interface itself, but first in the bulk material. Our simulations allowed for the first time a quantification of the reduction of the yield strength as temperature increases for such single interlamellar interfaces in TiAl.
In general, we observed in our simulations a quasi-brittle behavior with a sudden failure for the defect-free interface models, while a pronounced plasticity was observed for the more realistic interface models with a preexisting void. Beside the failure behavior, the evolution of the dislocation density was quantified by our simulations and it turned out that a dislocation burst at a high stress occurs for the defect-free case, while the dislocation density evolves gradually at lower stresses for interfaces with preexisting voids.
The evaluation of the damage behavior of the interface itself (with the effect of the bulk deformation having been eliminated by Eq. 16) revealed a TS behavior of the interface without significant elastic separation. Rather, the main contribution of the overall interface separation until failure came from interface plasticity and interface damage, which could be clearly distinguished in the TS law. Without elimination of the bulk deformation in our evaluation, a very different picture would have resulted where the TS curves would have exhibited a significant elastic regime. It is important to underline this because to the authors' best knowledge, comparable studies of other groups where traction-separation laws were deduced from MD simulations did not eliminate the effect of bulk deformation from their calculations as we did above. Our study indeed suggests that doing so makes a substantial difference and would therefore be worthwhile to be taken into consideration also by other groups in future work to characterize the physical properties of interfaces more accurately. Matching cohesive zone models with the results of MD simulations, plasticity could be correlated with an increase in dislocation density in the MD simulations, which underlines that the differentiation between elastic and plastic separation in our cohesive zone models is not simply a mathematical trick to cover a wider range of functional relations but that it captures indeed a specific form of physics that can be clearly observed on the atomic scale. Particularly, the damage-corrected plastic separation (1 − d)α parameter compares qualitatively well with the dislocation densities from the MD simulations. In fact, this study is to the authors' best knowledge the first one that deduces not only TS curves from MD simulations but also captures the inelastic deformation of the interface in a physically apparently meaningful way.
We also studied the effect of initial defects at the interface. They were found to decrease the interface strength as expected. However, for the interface separation energy, now such clear trend is observed in our MD simulations.
A comparison with results of other groups reported in the literature suggests that our MD simulations are suitable for identifying the interface separation energy correctly. By contrast, our results for the interface strength as observed in MD simulations should be taken with caution. They seem to overestimate the interface strength compared to experimental observations, possibly due to too small a simulation domain typical for nanoscale regime often covered in MD. This problem has also been encountered by others who tried to identify interface strengths by MD simulations, for example, for aluminum. Apparently, further research is required to overcome this key problem in the multi-scale simulation of materials.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
HG is responsible for all MD simulations; IS is responsible for all cohesive zone modeling issues; CC supervised the modeling work and contributed to its interpretation as well as the preparation of the manuscript.