Crystal Plasticity Modeling of Creep in Alloys with Lamellar Microstructures at the Example of Fully Lamellar TiAl

A crystal plasticity model of the creep behavior of alloys with lamellar microstructures is presented. The model is based on the additive decomposition of the plastic strain into a part that describes the instantaneous (i.e., high strain rate) plastic response due to loading above the yield point, and a part that captures the viscoplastic deformation at elevated temperatures. In order to reproduce the transition from the primary to the secondary creep stage in a physically meaningful way, the competition between work hardening and recovery is modeled in terms of the evolving dislocation density. The evolution model for the dislocation density is designed to account for the significantly different free path lengths of slip systems in lamellar microstructures depending on their orientation with respect to the lamella interface. The established model is applied to reproduce and critically discuss experimental findings on the creep behavior of polysynthetically twinned TiAl crystals. Although the presented crystal plasticity model is designed with the creep behavior of fully lamellar TiAl in mind, it is by no means limited to these specific alloys. The constitutive model and many of the discussed assumptions also apply to the creep behavior of other crystalline materials with lamellar microstructures.


INTRODUCTION
In high-temperature applications, the resistance to creep and thermomechanical fatigue are the key criteria for the choice of structural materials. After decades of intense research, the mechanisms behind the creep behavior of the well-established high-temperature structural materials like, for example, Ni-based super alloys are well understood, and as a result, their creep resistance has been optimized to a great extent. As the weight of components and the resulting mass forces are highly relevant in most high-temperature applications, there is, however, an ongoing search for lightweight alternatives to these commonly used alloys. In this context, intermetallic titanium aluminide (TiAl) alloys have frequently been discussed as promising candidates to replace the much heavier conventional alloys (Appel et al., 2011). Due to their beneficial combination of high specific strengths and good thermomechanical properties (Appel et al., 2011), TiAl alloys open up a significant weight-saving perspective, provided they are used to their full potential. This does, however, necessitate profound understanding of their high-temperature (creep) behavior ideally manifested in a comprehensive constitutive model.
In TiAl alloys, creep becomes relevant, starting from 650°C to 700°C Appel et al. (2011) and is influenced by many factors like composition, processing, or microstructure Appel et al. (2011), Zhang andDeevi (2003). When it comes to creep resistance, TiAl alloys with fully lamellar microstructures have consistently been reported to show the best properties with creep rates that are at least one order of magnitude smaller than in duplex microstructures Appel et al. (2011), Zhang and Deevi (2003), Maruyama et al. (1997).

Fully Lamellar Microstructures in TiAl
In TiAl alloys, fully lamellar microstructures mainly consist of two intermetallic phases, namely, the γ phase and to a minor fraction of the α 2 phase. Fully lamellar microstructures originate from a eutectoid phase transformation during which the globular parent grains become subdivided into thin straight lamellae (see Appel et al. (2011) for details). The lamella interfaces are strictly parallel so that these lamellar grains-the so-called colonies-are distinguished by the normal of their lamella interfaces.
Within each colony, there is a strict orientation relation between the lattices of the γ and the α 2 lamellae in which their close packed planes are arranged strictly coplanar Inui et al. (1992a), Appel et al. (2011). As illustrated in Figure 1, several orientations of the γ lattice fulfill the respective orientation relation with respect to the α 2 lattice. Therefore, the γ lamellae are further subdivided into the so-called domains of different crystal orientation Appel et al. (2011), Inui et al. (1992a). Figure 1 illustrates these coexisting types of microstructural interfaces with their corresponding spacing.
Due to the strict orientation relation between the lattices of the γ and the α 2 lamellae, it is possible to uniquely categorize their slip and twinning systems based on their slip/twinning direction s and the slip/twinning plane normal n (Lebensohn et al., 1998): • longitudinal (s lamella interface; n⊥ lamella interface), • mixed (s lamella interface; n⊥ / lamella interface), and • transversal (s ¤ lamella interface ; n⊥ / lamella interface). TABLE 1 | Slip and twinning systems in the tetragonal γ and hexagonal α 2 phase with morphological classification according to Lebensohn et al. (1998).
Although both phases exhibit a c a ratio ≠ 1, standard Miller index, respectively, Miller-Bravais index notation is used for a better readability. Throughout this article, the index α is used for slip systems, whereas the index β is used for twinning systems. Table 1 summarizes all slip and twinning systems in the γ and the α 2 lattice with their corresponding morphological classification. Figure 2 illustrates the deformation modes that correspond to the different morphological classes.

Polysynthetically Twinned Crystals
Investigating the plastic deformation of fully lamellar microstructures and in particular the underlying micromechanics is a cumbersome task and often suffers from problems in separating the contributions of different microstructural boundaries or deformation mechanisms/ systems. In fact, most of the understanding of the micromechanics of fully lamellar TiAl has been gained from experiments with the so-called polysynthetically twinned (PST) crystals. PST crystals are macroscopic specimens that only contain lamellae of one specific orientation, that is, they basically represent a single colony without the influence of the neighboring colonies. First reported in Fujiwara et al. (1990), PST crystals helped to reveal many of the intricate details of the highly anisotropic, plastic deformation of the lamellar colonies in fully lamellar TiAl (see, e.g., Kishida et al. (1998), Umakoshi and Nakano (1992, 1993), Fujiwara et al. (1990, Inui et al. (1992b), Yao et al. (1995), Umeda et al. (1997), Uhlenhut (1999), Bartels and Uhlenhut (1998)) and served as the starting point for most micromechanical modeling efforts (see, e.g., Lebensohn et al. (1998), Zambaldi et al. (2011), Roos et al. (2004), Parteder et al. (1995), Schlögl and Fischer (1996, 1997b, 1997a, Grujicic and Batchu (2001), Werwer andCornec (2000, 2006)).

Creep Characteristics of Fully Lamellar TiAl
Throughout this article, creep is discussed in terms of the commonly observed three creep stages Kassner (2009), which are schematically shown in Figure 3.
During primary creep, the density of stored dislocations ρ dis increases. As a result of the corresponding work hardening, the creep rate in the primary creep stage decreases over time/plastic strain. With increasing dislocation density, the recovery/ annihilation of dislocations increases until it eventually catches up with the dislocation generation resulting in an effectively constant dislocation density. Once the dislocation density reaches its saturation value ρ dis sat , there is effectively no additional work hardening, and the creep rate remains constant. This stage is known as secondary or steady-state creep. After a certain time/creep strain, damage sets in and leads to a drastic increase of the strain rate and ultimately to failure of the specimen/component. This stage of accelerating the creep rate is termed tertiary creep.

Frequent Absence of Secondary Creep Stage
Most of a component's life is usually spent in the secondary creep stage. Therefore, the creep resistance of different alloys is often compared in terms of their steady-state creep rates _ ε ss where the creep resistance is said to be higher when _ ε ss is lower . However, in TiAl alloys, a pronounced secondary creep stage is often not observed Appel et al. (2011), Kassner (2009. Instead, primary creep is often directly followed by tertiary creep, that is, by an increasing creep rate Appel et al. (2011). As illustrated in Figure 4, this manifests in an inflection point (instead of a plateau) in the creep rate vs. creep strain plot. As this inflection point-per definition-appears before the saturated state is reached, the respective experiments are regularly compared/quantified in terms of their minimum creep rate _ ε min (i.e., the creep rate at the inflection point) instead of _ ε ss . As it will become clear in the course of this article, the frequent absence of a pronounced secondary creep stage has led to some misconceptions in the theoretical assessment and comparison of creep experiments with fully lamellar TiAl.

Gradually Changing Stress Exponent
The stress dependence of the creep behavior is commonly discussed in terms of the creep law which relates the steadystate creep rate _ ε ss to the applied stress σ and absolute temperature θ via an Arrhenius type equation Kassner (2009).
In this, A is a material constant and n is the stress exponent. Q c denotes the activation energy for creep and k B is the Boltzmann's constant. As Eq. 1 is used to describe the steady-state creep rate, it is implicitly assumed that the dislocation density has reached its saturation value ρ dis sat and that the strength of the material consequently does not change anymore. Therefore, Eq. 1 can conveniently be evaluated for a given stress and temperature without the need to consider any evolution equation(s). However, Eq. 1 does neither allow any assessment of the primary creep stage nor of the onset of creep damage and is thus only of limited use for predicting the creep behavior. For TiAl alloys, Eq. 1 has frequently been used to evaluate and discuss the stress dependence of their steady-state or minimum creep rates (see, e.g., Appel et al. (2011), Wegmann et al. (2000, Zhang and Deevi (2003), Marketz et al. (2003)). The experimentally determined stress exponents do, however, show a significant scatter and often are much higher than the values of 3 < n < 5 that are typically expected for dislocation creep Appel et al. (2011), Kassner (2009. In fact, the stress exponent itself has regularly been found to be stress dependent increasing from 1 < n < 2 at very low stresses to 7 < n < 9 at high stresses Zhang and Deevi (2003), Appel et al. (2011). This contradicts the whole idea of the existence of a single, constant stress exponent and questions the suitability of Eq. 1 to generally describe the creep behavior of fully lamellar alloys over a broad range of stresses Zhang and Deevi (2003), Appel et al. (2011). Therefore, the experimentally observed stress exponents are termed apparent stress exponents in the following.
The gradual change in the apparent stress exponent is most likely associated with the transition between different ratelimiting creep mechanisms. For very low stresses (and reasonably high temperatures), there is evidence that interface glide is the dominant creep mechanism Zhang and Deevi (2003), Appel et al. (2011) with apparent stress exponents of 1 < n < 2. At intermediate stresses, the apparent stress exponents of 3 < n < 5 fit the expected values for dislocation creep Kassner (2009). The cause of the apparent stress exponents of n > 5 that have frequently been observed for high stresses does, however, appear to be not yet clear.
A possible explanation for the unusually high apparent stress exponents may be found in the frequent absence of a secondary creep stage or more precisely in the fact that the respective experiments are compared in terms of their minimum creep rate. The direct transition from primary to tertiary creep as it is  Figure 4-per definition-appears before a saturation state has been reached. Else at least a short steadystate plateau would occur in the creep rate vs. strain graph. Consequently, the corresponding minimum creep rate is higher than the creep rate that would occur in the saturated state if damage would not have set in already. Thus, it is to be expected that the stress exponent of the minimum creep rates is higher than the maximum value n 5 for dislocation creep. Consequently, the apparent stress exponents of n > 5 that have been observed for high stresses are not related to another creep mechanism but most likely describe the stress dependence of the onset of tertiary creep which gradually reduces the extent of the secondary creep stage with increasing stress until it finally disappears completely.

Influence of Microstructural Interfaces
While it is well accepted that the strengthening by the different microstructural boundaries in fully lamellar TiAl can be described by a Hall-Petch type relation Appel et al. (2011), Umakoshi and Nakano (1993), Dimiduk et al. (1998), their influence on the creep resistance has only been reported qualitatively Zhang and Deevi (2003), Maruyama et al. (1997), Chatterjee et al. (2002), Marketz et al. (2003). In general, it has been stated that the creep resistance increases with decreasing lamella thickness λ L , while the colony size λ C appears to only affect the creep resistance for λ C < 100 μm Appel et al. (2011), Maruyama et al. (1997, Zhang and Deevi (2003). However, as there is evidence that for high temperatures and very low stresses the creep rate is determined by interface glide, the microstructural interfaces can-contrary to their overall beneficial effect-also become the carriers of creep deformation for certain specific load cases Maruyama et al. (1997), Appel et al. (2011). Experiments with PST crystals-as quasi prototypes of single isolated colonies-revealed a strong anisotropy in their creep behavior, that is, for a given stress, the minimum/steady-state creep rates vary with the angle between lamella plane and loading direction Parthasarathy et al. (2000), Wegmann et al. (2000), Asai et al. (2002), Kim et al. (2002). The anisotropy in the creep resistance of PST crystals shows characteristics that are similar to what has been observed for their yield stress. For the yield stress, there is a region of intermediate loading angles around 45°w.r.t. the lamella plane for which it is significantly lower than in the extreme loading angles near 0°and 90°Fujiwara et al. (1990), Uhlenhut (1999), Umeda et al. (1997), Yao et al. (1995), Inui et al. (1992b), Umakoshi and Nakano (1992). Analogously, the creep resistance of PST crystals has been found to be low for intermediate loading angles around 45°and significantly higher for loading angles near 0°and 90°Parthasarathy et al. (2000), Wegmann et al. (2000), Asai et al. (2002), Kim et al. (2002). This anisotropy is directly related to the influence of the lamella boundaries which are strong barriers to dislocation motion and twin propagation. In terms of the morphological classification, this means that mixed and transversal slip/ twinning systems are generally much stronger than longitudinal systems which are only strengthened by the weak and widely spaced domain boundaries.

Modeling Creep of Fully Lamellar
TiAl-State of the Art So far, remarkably few studies have dealt with modeling creep of fully lamellar TiAl alloys. Depending on the level of pragmatism, the reported models followed a more phenomenological approach by relating the minimum/steady-state creep rate to the applied stress without considering the details of primary creep and its transition to secondary creep Marketz et al. (2003), Zhang and Deevi (2003) or were designed to explicitly consider the underlying micromechanics Chatterjee et al. (2002), Estrin and Mecking (1984). To the best of the authors' knowledge, the onset of tertiary creep in TiAl alloys has not been modeled at all so far.
1.3.1 Marketz et al. (2003) In order to capture the anisotropy in the creep behavior of single colonies and the influence of the lamella boundaries, in Marketz et al. (2003), the material coefficient A of the classical model (i.e., Eq. 1) has been defined to be a function of the lamella thickness λ L and the loading angle φ. With constant stress exponents and activation energies taken from γ and α 2 single crystal creep experiments, _ ε c min and _ ε α 2 min were set up separately for both phases. The model was then calibrated against experimentally determined minimum creep rates for differently oriented PST crystals reported in Wegmann et al. (2000). From this, a polynomial fit for A(φ, λ L const.) was obtained. The dependence of A on λ L was fitted against experiments with polycolony sheet material, while reasonably assuming that the relation A(φ) of the lamella colonies within a polycolony microstructure is the same as for PST crystals. Subsequently, the resultant model was applied to an RVE of a polycolony microstructure and was successful in reproducing the corresponding minimum creep rates for different lamella thicknesses.
While certainly being a pragmatic and to an extent successful approach, the model from Marketz et al. (2003) still lacks some generality as it largely relies on phenomenological fits. In particular, the choice of a constant stress exponent of n 7 for creep in the γ lamellae is questionable as comprehensively described above. As the model from Marketz et al. (2003) is based on Eq. 1, it naturally suffers from the same restrictions, that is, it does not allow any statements on the extent of primary creep and the onset of tertiary creep.

Zhang and Deevi (2003)
In order to phenomenologically capture the observed increase in the stress exponent with increasing stress, Zhang and Deevi (2003) replaced the power law relation in Eq. 1 by the hyperbolic sine relation In this, σ 0 is a fitting constant that represents a back-stress for dislocation motion. With this model, the steady-state/minimum creep rates of a broad variety of experimental results have been fitted reasonably well. Only for the low stress regime, Eq. 2 failed to reproduce the experimental data. Thus, in this regime, the steady-state/minimum creep rates have been fitted using a variation of Eq. 1 with a stress exponent of n 2 in order to represent interface sliding. As opposed to the classical creep law, Eq. 2 relates the applied stress to a back-stress σ 0 . From the fitting procedure in Zhang and Deevi (2003), it has been found that σ 0 generally increases with decreasing λ L . This has qualitatively been associated to the strengthening effect by lamella interfaces.
From an engineering point of view, the model presented in Zhang and Deevi (2003) suffers from the same lack of generality as the model from Marketz et al. (2003). While successfully fitting the steady-state/minimum creep rates of a variety of experimental studies including the gradual change in the apparent stress exponent, the model from Zhang and Deevi (2003) does not allow any statements on the extent of the primary or the onset of the tertiary creep stage and is thus only of limited use for predictive purposes.

Chatterjee et al. (2002)
In order to overcome the limitations of classical creep laws to only model the constant creep rate in the secondary creep stage (or the minimum creep rate respectively), the model from Chatterjee et al. (2002) explicitly takes into account the competititon between work hardening and recovery, which causes the transition from primary to secondary creep stage.
In order to model the gradual change of the strain rate from primary to secondary creep stages, in Chatterjee et al. (2002), an evolution equation for the stored dislocation density _ ρ dis f (ρ dis , _ ε, θ) is combined with an equation that describes the strain rate dependence of the flow stress σ f (ρ dis , _ ε, θ) for a given ρ dis . In this, the evolution of ρ dis with plastic strain has been modeled in a generation-recovery format Estrin and Mecking (1984), Chatterjee et al. (2002).
with _ ρ dis gen denoting the rate at which dislocations are immobilized and _ ρ dis dynrec describing the dynamic recovery of dislocations. Furthermore, b is the magnitude of the Burgers vector and L corresponds to the free path length of dislocations. R dyn is the dynamic recovery coefficient.
Evolution equations in the form of Eq. 3 have been successfully applied in many dislocation density-based models (see, e.g., Anand et al. (2015), Evers et al. (2004), Beyerlein and Tomé (2008)). As their mutual interaction is one of the main causes for dislocations to get immobilized, the free path length L in Eq. 3 is commonly assumed to be correlated to the dislocation spacing L dis : with ζ being a fitting constant.
The hybrid one parameter model from Estrin and Mecking (1984) on which the model from Chatterjee et al. (2002) is largely based does, however, take into account that in most materials there are additional obstacles for dislocation motion such as point defects or grain boundaries all with their own spacing L. These different coexisting obstacles all increase the rate at which dislocations become immobilized, that is, contribute to _ ρ dis gen in Eq. 3. Consequently, the influence of the lamella interfaces on the dislocation storage was modeled in Chatterjee et al. (2002) by introducing an additional free path length L L with ξ again being a fitting constant. Following Estrin and Mecking (1984), has been introduced into Eq. 3 in order to model the combined effect of dislocation interaction and lamella interfaces on the dislocation storage.
In addition to Eq. 3, the kinetic equation has been set up in Chatterjee et al. (2002), where σ Y denotes the current yield strength that is modeled to evolve as a function of ρ dis via a classical Taylor hardening law: In this, _ ε 0 is a reference strain rate and n inst is the instantaneous strain rate sensitivity exponent. M denotes the Taylor factor, and μ denotes the shear modulus.
Combining Eqs. 3, 7, and 8, integrating the resultant equation and solving for σ const. 1 , yields an analytical relation between the creep rate and the creep strain which has subsequently been calibrated against experimental creep curves obtained from sheet material of fully lamellar TiAl.
As the model explicitly considers the evolution of the dislocation density toward the saturated state, it is capable of reproducing the transition from primary to secondary creep stage and thus is-from an engineering point of view-of much more relevance than the classical creep laws. Unfortunately, the model was calibrated in such a way that the saturation state coincides with the minimum creep rates that have been found in the respective experiments. As discussed earlier, this is not physical as the minimum creep rate occurs at the onset of tertiary creep and is thus not related to the saturated state. Furthermore, neither interface glide nor the tertiary creep stage is considered in Chatterjee et al. (2002), so that the model consequently fails to describe the respective deformation regimes. However, the model from Chatterjee et al. (2002) still represents an important step toward a comprehensive model of the creep behavior of fully lamellar TiAl, and many of the underlying assumptions will be adopted in the course of the present article.

Conclusion
In the context of the frequent absence of a secondary creep stage and the gradual change in the apparent stress exponent, classical creep models that are based on Eq. 1 are not well suited to describe the creep behavior of fully lamellar TiAl. While to a certain extent being successful in describing the stress dependence of the minimum creep rates by using different stress exponents for different stress regimes, these models miss some vital details of the creep behavior. Even if correctly predicting the minimum creep rate for a given stress, the corresponding creep strain is not known from models based on Eq. 1. As the minimum creep rate does, however, occur at the onset of tertiary creep and thus marks a deformation state the occurrence of which has to be avoided by any means in engineering design, such models are of little use for engineering purposes. Thus, a precise description of the primary creep stage is vital for modeling the creep behavior of fully lamellar TiAl.
Thus, some of the ideas of the previously reported models (especially Chatterjee et al. (2002), Estrin and Mecking (1984)) will be picked up in the following and extended to a crystal plasticity framework in order to model the creep behavior of fully lamellar TiAl in a microstructure sensitive way.

Base Model
In previous works, Schnabel and Bargmann (2017), Schnabel et al. (2019), a thermomechanically coupled, defect densitybased crystal plasticity model has been set up with the focus on precisely capturing the microstructure sensitivity of the yield and work hardening behavior of fully lamellar TiAl. This model incorporates the Hall-Petch strengthening effect of the different coexisting microstructural boundaries as well as the work hardening interactions between twins and dislocations in the γ phase. Furthermore, the evolution of the dislocation density in this model takes into account dynamic and static recovery.
As this model has been extensively discussed elsewhere (Schnabel and Bargmann (2017), Schnabel et al. (2019)), the governing equations are just briefly summarized in the following. The meaning of the used symbols can be found in Table 2.
Kinematics: Split of deformation gradient Plastic deformation gradient Plastic velocity gradient 2 Interaction coefficients between slip system α And non-coplanar (ncp) twinning systems β h β β′ Interaction coefficients between twinning system β And non-coplannar (ncp) twinning systems β′ C βα Interaction coefficients between twinning system β And dislocations on slip systems α Bulk modulus E Young's modulus ] Poisson's ratio α t Thermal expansion coefficient S 0 Absolute entropy density c p Heat capacity Q −κGradθ Heat flux vector κ Thermal conductivity r External heat supply per unit mass 2 cf. Rice (1971), Kalidindi (1998). Note that as opposed to the original model from Kalidindi (1998), no reorientation of the twinned region and subsequent slip within the twins is considered here. This assumption is reasonable as long as the twins remain very thin as it is the case in lamellar microstructures.
Elastic right Cauchy-Green tensor Stress measures: 2nd Piola-Kirchhoff stress (15) Flow/twinning rule: Twinning rate 4 CRSS/Work Hardening (WH): CRSS of twinning systems Hall-Petch for slip systems Hall-Petch for twinning systems WH of slip systems 5 WH of twinning systems 6 Defect density evolution: Thermomechanical coupling Helmholtz free energy 9 2nd Piola-Kirchhoff stress Temperature evolution As the α 2 phase only constitutes few Vol.-% (usually less than 10%) in commercial TiAl alloys, its contribution to the plastic 3 cf. Pierce et al. (1983). 4 cf. Kalidindi (1998 deformation is hard to assess from experiments with two-phase alloys. Therefore, the α 2 phase is modeled in a simplified way, by using the (temperature dependent) initial critical resolved shear stresses from single crystal experiments and neglecting Hall-Petch strengthening effects as described in Butzke and Bargmann (2015). The model that is described by Eqs. 9-28 has been successfully applied to investigate different aspects of the plastic deformation of fully lamellar TiAl alloys from room to operating temperature Schnabel and Bargmann (2017), Schnabel et al. (2019), Schnabel (2018). However, in order to capture the creep behavior of fully lamellar TiAl, this model needs to be extended as it will be shown in the following.

Model Extension to Creep
Most equations and assumptions of the base model are valid also for modeling creep deformation. It has to be noted, however, that in Equation (11), the kinematics of dislocation climb is neglected. This appears reasonable as long as the creep deformation is dominated by (climb-assisted) glide. In line with the arguments of the one-parameter model from Chatterjee et al. (2002) and Estrin and Mecking (1984), the work hardening behavior only depends on the dislocation densities and the twinned volume fractions, that is, the respective work hardening model Eqs. 18-23 is applicable without any changes. The same holds for the thermomechanical coupling and the resultant stress as well as the temperature evolution.
The only parts of the model that will be altered/enhanced in the following are the flow rule (Eq. 16 in the base model) and the dislocation density evolution model (Eq. 24 in the base model).

Shear Rates
In crystal plasticity, the kinematics of the plastic deformation of a crystal's lattice is described in terms of shear (rates) on its crystallographic slip and twinning planes; cf. Eq. 11. In this, the shear rates _ c α and c T _ g β are the continuum representation of the collective shear deformation of a crystal by all dislocations/ twins on its slip/twinning systems α/β. The actual motion of dislocations does, however, take place in a successive sequence of mobilization of pinned line segments which then glide on their slip plane until being pinned once again. Once mobilized, dislocations move with a velocity that is close to the speed of sound. Thus, the effective shear rate on a slip system is essentially determined by the rate at which dislocations are mobilized.
In order for a dislocation to become mobile, there has to be a sufficient driving force (stress and/or thermal activation) that supports the dislocation to overcome obstacles and/or the local stress field that hinders its motion. Depending on the combination of stress and temperature, different characteristic types of dislocation mobilization can become rate limiting. The intricate details of the discrete motion of dislocation line segments and their interactions with point defects and the surrounding dislocation network are naturally beyond the scope of continuum mechanics. However, it is still possible to incorporate different rate-limiting micromechanical effects into a crystal plasticity model by additively decomposing the shear rate _ c α on slip system α, that is, by expressing it as the sum of the shear rates that correspond to different mobilization mechanisms as proposed in Cassenti (2008, 2011).
In the context of creep, essentially two different strain rate limiting mechanisms of dislocation mobilization occur: • Instantaneous plasticity: If the applied (macroscopic) stress is above the (macroscopic) yield stress at the given temperature, the resolved shear stress on at least one slip system is sufficiently high to force the dislocations on that slip system to move through/past any obstacles or local stress fields. This (yield) deformation takes place with a-in the context of creep-very high strain rate and is thus referred to as instantaneous plasticity in the following. • Viscous glide: If the (macroscopic) stress is below the yield point, the resolved shear stress on neither of the slip systems is sufficient for the dislocations to break free from any pinning points or local stress fields. If the temperature is, however, sufficiently high, pinned dislocations may climb to parallel slip planes to bypass obstacles or escape local stress fields. Once a dislocation has climbed to a slip plane above the obstacle or with a (locally) relaxed stress field, it may glide even under the given resolved shear stress until getting pinned again. The strain rate of this climb-assisted glide is thus limited by the diffusion kinetics.
Considering these two rate limiting cases, the slip shear rates are written as the sum (cf. Cassenti (2008, 2011)) As twinning partial dislocations are bound to their glide plane, the corresponding twinning rates do not have a viscoplastic creep part and thus simply read

Instantaneous Plasticity
Instantaneous plasticity occurs whenever the resolved shear stress τ on a slip or twinning system approaches its current strength. This behavior is modeled by Eqs. 16 and 17 in the base model, which are consequently adapted here to model the instantaneous shear rates for slip and twinning Unless stated otherwise, _ c inst 0 0.001 and n inst 50 are assumed throughout this work.

Viscoplastic Response
The viscous shear rates on slip systems α are modeled by (cf. Cassenti (2008, 2011)) in order to represent the temperature-dependent climb behavior which is rate limiting for the viscoplastic response. In this, _ c visc 0 denotes a reference creep rate. The strain rate sensitivity exponents should generally obey n visc ≪ n inst with 3 ≤ n visc ≤ 5 for climb-assisted glide.

Defect Density Evolution
The work hardening relations, that is, Eqs. 22 and 23, describe the strength of slip and twinning systems as a function of their current surrounding microstructure. As this microstructure evolves with plastic strain, temperature, and time, an evolution model for the corresponding internal variables-namely, the dislocation densities ρ dis α and twinned volume fractions f β -is required.
In the following, the evolution of the dislocation density on a slip system α is modeled by the generation recovery format As opposed to the dislocation density evolution models in, for example, Chatterjee et al. (2002), Estrin and Mecking (1984), Evers et al. (2004), Beyerlein and Tomé (2008) which incorporate the accumulation of dislocations _ ρ dis α,gen and their dynamic (i.e., the deformation dependent) recovery _ ρ dis α,dynrec , Eq. 34 also considers the static (i.e., the time dependent) recovery behavior via _ ρ dis α,statrec . The dislocation accumulation term _ ρ dis α,gen is modeled as in Eq. 3, but in a slip system resolved form, In the following, some of the ideas from Estrin and Mecking (1984), Chatterjee et al. (2002) are adapted in order to make Eq. 35 microstructure sensitive. In the model from Chatterjee et al. (2002), Estrin and Mecking (1984), the free path length L was defined by combining the dislocation spacing and the obstacle/lamella interface spacing as shown in Eq. 6. While this is an inevitable (and by no means bad) assumption to combine both effects in an analytical, macroscopic model, crystal plasticity models allow the individual definition of the spacing of obstacles for every slip system (thus the subscript α in L α ). The free path length L α in Eq. 35 depends on the orientation of the respective slip system w.r.t. the lamella interfaces, i.e., on its morphological class. Slip systems whose slip plane crosses the lamella interface (mixed and transversal) are restricted by the lamella thickness λ L , while for slip systems with slip direction and slip plane parallel/co-planar to the lamella interfaces (longitudinal systems), the free path length is given by the domain size λ D . 10 Following Eq. 5, the free path length of longitudinal slip systems is thus defined as while the free path length of mixed and transversal systems reads However, the free path length is only defined by the domain or lamella boundaries as long as the dislocation density on the respective slip system is low. With ongoing plastic deformation, dislocations begin to pile up against these boundaries, steadily reducing the free path length. This effect is modeled in a continuum sense by adapting the dislocation spacing from Eq 4, that is, and evaluating the free path length L α in Eq. 35 via for each slip system α in each time step of the simulation. The stored dislocation density may also decrease by dynamic and static recovery (see Eq. 34). Dynamic recovery describes the recombination of dislocations that move on the same slip plane and thus only takes place for a non-zero slip shear rate _ c α . Static recovery, on the other hand, describes the recombination of dislocations on different slip planes by diffusion-assisted climb and is thus highly temperature dependent.
The dynamic recovery rate _ ρ dis α,dynrec is modeled following the well-established formulation from Estrin and Mecking (1984) As the recombination of dislocations can only take place between dislocations on the same slip system, the dislocation density ρ dis α was chosen to be the relevant measure for dynamic recovery on slip system α in Eq. 40. The static recovery behavior is modeled according to the base model (cf. Anand et al. (2015) and Roters et al. (2010), Ma et al. (2006)) The expression 〈•〉 in Eq. 41 takes a value of x for 〈x〉 > 0 and of 0 for 〈x〉 < 0 and thus ensures that static recovery only takes place when the dislocation density ρ dis α is above a minimum threshold 10 In the α 2 phase, obviously no domain boundaries exist, so that the corresponding free path length would be the colony size in case of polycolony microstructures. However, most of the statements presented here refer to the γ phase as the α 2 phase only constitutes several Vol.-% and is thus of minor interest. Throughout this work, the dislocation density is assumed to evolve starting from an initial value of ρ dis α,init 10 5 1 mm 2 for all slip systems α which corresponds to a total dislocation density of ρ dis 1.2 × 10 6 1 mm 2 . This is also assumed to mark the minimum dislocation density, that is, ρ dis α,min ρ dis α,init .

APPLICATION TO PST CRYSTALS
In the following, the modified thermomechanically coupled, defect density-based crystal plasticity model is applied to reproduce and discuss some aspects of the creep deformation of PST TiAl crystals.

Simulation Setup
As in previous studies Schnabel and Bargmann (2017), Schnabel et al. (2019), a simplified representative volume element (RVE) of a PST crystal is set up (see Figure 5).
In this, the thickness of the single modeled α 2 lamella is adjusted while fixing the RVE's length in x-direction in order to fit the reported α 2 volume fraction of the experiments to which the simulations are compared. The six different γ orientation variants are assumed to all have the same volume fraction. The orientation relation between the crystal lattices of the α 2 phase and the different γ orientation variants are set according to Figure 1 with the [110] direction of the c I M domains either in y-or z-direction according to the experimental characterization in each case. The geometry is meshed using 42 linear hexahedral finite elements. The RVE is subjected to periodic boundary conditions in order to mimic an infinite repetition of the RVE in all directions. 11, The periodic RVE is then embedded into a dummy element with negligible stiffness. The displacements of the periodic RVE's master nodes are coupled to the averaged deformation of the dummy element via a rotation relation as introduced in Werwer and Cornec (2000), Werwer (2005) and described in Schnabel (2018). This rotation relation allows the investigation of the RVE's constitutive response as a function of the loading angle φ (i.e., its anisotropy in the xz-plane) by applying a uniaxial deformation to the dummy element.

Thermoelasticity
Although the elastic deformation is usually not discussed in detail in the context of plasticity, the (thermo-) elastic model parameters still have to be chosen properly as they, e.g., occur in the work hardening relation via the shear modulus μ. The thermoelastic parameters (i.e., E(θ), ](θ), c a (θ), ρ 0 (θ), c p (θ), κ(θ), and α t (θ)) that are used throughout this article are extracted from the literature as indicated in Table 3.

Instantaneous Plasticity
As discussed earlier, the constitutive response is dominated either by instantaneous or by viscoplastic deformation depending on stress level, applied strain rate, and temperature. Instantaneous plasticity occurs whenever the resolved shear stress on any slip/ twinning system is above its current strength. Thus, instantaneous plasticity dominates the deformation behavior for experiments with applied stresses above the yield stress and constant strain rate experiments with strain rates that are above typical creep rates at the given temperature. Furthermore, the instantaneous plastic deformation depends on temperature only through the temperature dependence of the yield stress and the thermoelastic parameters.
Making use of this fact, it is possible to separately calibrate and evaluate the parts of the model that describe the instantaneous plasticity by only considering load scenarios that favor this deformation behavior.

Constant Strain Rate Room Temperature Experiment
Classical uniaxial tensile/compression tests at room temperature are naturally dominated by instantaneous elastic-plastic deformation. Correspondingly, such room temperature experiments are a good starting point to calibrate the respective parts of the discussed model.
At room temperature, static recovery and viscoplastic (creep) deformation are negligible. Consequently, _  11 Note that this simplifies the RVE to a columnar structure in y-direction.
Frontiers in Materials | www.frontiersin.org February 2021 | Volume 7 | Article 581187 been considered. However, the enhanced dislocation density evolution model requires the newly introduced model parameters (i.e., ξ D , ξ L , ζ, and R dyn ) to be identified, and also necessitates a (slight) adjustment of previously identified work hardening interaction coefficients. The parameters for the work hardening model (Eqs. 22 and 23) and the dislocation density evolution (Eq. 34) are calibrated against the room temperature compression tests with differently oriented PST crystals that are reported in Uhlenhut (1999). The corresponding calibration procedure has extensively been discussed before Schnabel and Bargmann (2017), Schnabel et al. (2019) and will therefore not be revisited here. A summary of the corresponding simulation results can, however, be found in the Supplementary Material, and some details about the identified parameter set are given in the following.
The work hardening due to the interaction of slip dislocations is modeled by classical Taylor hardening (cf. Eq. 22), and thus do not require any model parameters to be identified. The work hardening interaction between slip and twinning systems does, however, require the calibration of the corresponding interaction coefficients in Eqs. 22 and 23 (i.e., h αβ , h ββ′ and C βα ). In the most general case, there are different hardening interaction coefficients for all possible combinations of slip and twinning systems. However, when applying crystal plasticity models to specific materials, the slip/twinning systems can naturally be categorized into groups of systems with the same behavior (e.g., based on crystallography), which consequently can be described by the same set of model parameters. 3 | Thermoelastic model parameters used in this study. The reported experiments were carried out in the indicated temperature range. Where necessary, the linear approximations of the temperature dependences are extrapolated to the temperatures of interest. If no parameters for the single phases are available, the parameters of two phase alloys were used instead. The lattice constants a and the c a ratios of the unit cells are necessary to calculate the (temperature dependent) magnitude of the Burgers vectors and the corresponding slip/twinning directions s and slip/twinning plane normals n as described elsewhere Schnabel (2018   For fully lamellar TiAl, it has been shown by several studies that the number of model parameters of the crystal plasticity formulation can be drastically reduced by making use of the morphological classification (see, e.g., Lebensohn et al. (1998), Werwer and Cornec (2006), Zambaldi et al. (2011)). By assuming that the hardening interaction between all slip/twinning systems of one morphological class with all slip/twinning systems of another morphological class can be described by the same coefficient, the number of parameters reduces to the interaction matrix that is shown in Table 4. The calibration against the experiments from Uhlenhut (1999) (cf. Supplementary Material) moreover reveals that only 6 of the 16 remaining work hardening interaction coefficients are actually needed (i.e., have values ≠ 0) to reproduce the room temperature work hardening behavior. Having a closer look at the corresponding parameters in Table 4, it becomes clear that it is the Hall-Petch type hardening effect of emerging transversal twins on longitudinal and mixed slip/twinning systems that dominates the slip-twin work hardening interaction (cf. Schnabel and Bargmann (2017), Schnabel et al. (2019)). In addition, only longitudinal twins appear to be significantly hardened by an increasing number of slip dislocations (i.e., C βα ≠ 0 in Table 4).
The parameters of the dislocation density evolution (i.e., ξ D , ξ L , ζ, and R dyn ) were identified from the same simulations and take values that are comparable to what has been reported in Chatterjee et al. (2002). The results are given in Table 5.

Constant Strain Rate High Temperature Experiment
Although thermally activated processes become increasingly relevant with increasing temperature, uniaxial tensile/compression tests at elevated temperatures are still dominated by instantaneous elastic-plastic deformation as long as the strain rate is sufficiently high. Therefore, the viscoplastic (creep) deformation remains negligible and is thus again set to _ c visc α 0 in the following. As discussed above, instantaneous plastic deformation occurs as soon as the resolved shear stress on any deformation system reaches its current strength. In preparation for the next section in which the creep behavior of differently oriented PST crystals will be discussed at the example of the experiments from Parthasarathy et al. (2000), it is thus crucial to capture the temperature dependence of the model parameters that determine the current strength of deformation systems.
The work hardening as described by Eqs. 22 and 23 is not explicitly depending on temperature. Equations, 22 and 23 are only indirectly depending on temperature through the shear modulus μ and the magnitudes of the Burgers vectors b α/β . Therefore, the work hardening parameters identified from room temperature simulations/experiments (see Table 4) are applicable to high temperature problems without any changes. The same holds for the parameters of the dislocation generation and dynamic recovery model, that is, Eqs. 35 and 40 (see Table 5).
The yield stress of PST crystals and thus the initial slip/twinning system strengths τ Y α,0 and τ T β,0 (cf. Eqs. 18 and 19) are highly temperature dependent and exhibit a temperature anomaly as it frequently occurs in intermetallic alloys Appel et al. (2011), Veyssière (2001. This means that the yield stress counterintuitively increases with increasing temperature up to a certain peak strength.
Above the corresponding peak temperature, the yield stress decreases drastically with further increasing temperature.
In Butzke and Bargmann (2015), this yield stress anomaly of PST crystals has been modeled according to the experimental findings from Umakoshi and Nakano (1992) by making the Hall-Petch coefficients k HP i in Eqs. 20 and 21 temperature dependent. Following the same line of arguments, the critical stress concentration that needs to be exceeded by a dislocation pileup in order to overcome lamella respectively domain boundaries has been modeled as a function of temperature just recently in Ilyas and Kabir (2020). However, the temperature at which the anomalous peak occurs as well as the height of the peak varies from study to study Umakoshi and Nakano (1992), Inui et al. (1995), Parthasarathy et al. (2000). Thus, the models from Butzke and Bargmann (2015) and Ilyas and Kabir (2020) cannot generally be applied to any experimental study without recalibration.
Thus, the temperature-dependent initial critical resolved shear stresses τ Y α,0 and τ T β,0 (cf. Eqs. 18 and 19) are calibrated as a whole to match the experimental results from Parthasarathy et al. (2000) FIGURE 6 | Yield stress (σ 0.2 ) of differently oriented PST crystals. Comparison of simulation results to experiments from Parthasarathy et al. (2000).
TABLE 6 | Identified initial critical resolved shear stresses in the γ phase as a function of temperature for the PST crystals tested in Parthasarathy et al. (2000).

Initial critical resolved shear stresses in the γ phase
Parameter Unit ls (α 1-3) ms (α 4-6) ts (α 7-12) Temperature for different characteristic loading angles. As the evolution of the dislocation density does not play a significant role at the onset of yield (even at the σ 0.2 yield point), the static recovery behavior is neglected here (i.e., _ ρ dis statrec 0 is assumed again) for the sake of simplicity. Furthermore, the initial strengths of all slip/twinning systems of the same morphological class is assumed to be the same because of their identical free path lengths. Thus, only three initial critical resolved shear stresses remain to be identified as a function of temperature; one for longitudinal, one for mixed, and one for transversal slip/twinning systems each. Figure 6 shows the simulated yield stresses of PST crystals under the characteristic loading angles of 0°, 45°, and 90°over temperature in comparison to the experimental results from Parthasarathy et al. (2000). The corresponding initial critical resolved shear stresses for the γ phase are gathered in Table 6.
Although only covering the discrete temperatures from the creep experiments in Parthasarathy et al. (2000), Figure 6 clearly shows the yield stress anomaly of PST crystals. While for a loading angle of 90°, the peak respectively plateau temperature has obviously been already reached below 700°C (cf. Umakoshi and Nakano (1992), Inui et al. (1995)), the yield strength for loading angles of 0°is still increasing even at 815°C. For loading angles of 45°, the yield stress basically remains constant over temperature, as it has been found in other studies Umakoshi and Nakano (1992), Inui et al. (1995).

Viscoplastic Creep Deformation
In the following, the presented model is applied to reproduce the creep experiments from Parthasarathy et al. (2000). In Parthasarathy et al. (2000), PST crystals of the characteristic orientations 0°, 45°, and 90°were subjected to compressive stresses at temperatures of 700°C, 760°C, and 815°C.
As the dislocation density evolution (and thus the work hardening rate) has been modeled in a microstructure sensitive way (cf. Eqs. 36 and 37), the microstructural parameters λ L and λ D need to be set according to the microstructure characterization in the experiments. The lamella thickness has been reported to be λ L 0.37μm in Parthasarathy et al. (2000). The γ domain size has not been reported in Parthasarathy et al. (2000) and is thus chosen to be λ D 50λ L 18.5μm which corresponds to typically observed aspect ratios of the γ domains (cf., e.g., Nakano (1992, 1993)).
The instantaneous shear rates _ c inst α are modeled via a viscoplastic power law (cf. Eq. 31). While this formulation ensures a unique identification of active slip/twinning systems, it comes at the expense of introducing artificially high strain rate dependence if n inst in Eq. 31 is chosen too low. In order to minimize the unwanted contribution of _ c inst α to the viscoplastic deformation (i.e., in order to really separate the two deformation regimes), n inst 200 is chosen for the creep simulations.
With these assumptions, only the viscoplastic shear rates _ c visc α from Eq. 33 and the static recovery rate _ ρ dis statrec from Eq. 41 remain to be identified in order to reproduce the creep behavior of differently oriented PST crystals.
In the creep experiments from Parthasarathy et al. (2000), PST specimens have been subjected to a constant creep load until the strain rate vs. time plots reached a plateau. Once such a plateau has been reached, the load was increased. By repeating this stepwise stress increase procedure, it has been possible to extract the creep rates that correspond to the observed plateaus as a function of applied stress by using the same specimen. By applying different stress series to different specimens, it has further been shown in Parthasarathy et al. (2000) that the creep rates that correspond to the observed plateaus only depend on the applied stress irrespective of the load history. From this experimental procedure, a stress exponent of n 5 has been identified for all orientations and temperatures for which the PST crystals have been tested Parthasarathy et al. (2000). This is reasonable in the context of dislocation creep and suggests that a steady state (i.e., saturation of the dislocation density) has in fact been observed in the respective experiments as claimed by the authors. In consequence of the findings from Parthasarathy et al. (2000), the stress exponent in Eq. 33 is thus set to n visc 5 in the following.
With the stress exponent at hand, only the two Arrhenius type terms _ c visc 0 exp − Qc kBθ from Eq. 33 and R stat exp − QR kBθ in Eq. 41 remain to be identified in order to reproduce the creep behavior. As the creep behavior has been found to be insensitive to the α 2 volume fraction Appel et al. (2011) which is either way quite low, no viscoplastic deformation is considered in the α 2 phase. In the simulations, the creep load is applied via a linear ramp over a time of 1 s in order to facilitate numerical handling of the load jump. As soon as the simulations reached the creep strain that corresponds to the primary creep strain that has been reported in Parthasarathy et al. (2000), the load has been increased to the next level, again over a period of 1 s. If for a given creep load, the simulation did not show saturation behavior up to the experimentally determined primary creep strain, the load was kept on the respective level until a plateau appears in the strain rate vs. time plot. Figure 7 shows the steady-state creep rate vs. stress plot as extracted from the respective simulation results in comparison to the experimental results from Parthasarathy et al. (2000) for different lamella orientations and temperatures. Figure 8 shows the simulated strain rate vs. time plots that have been used to identify the steady-state creep rates in Figure 7 as well as the corresponding strain rate vs. strain plots.
The parameters that were applied to produce the simulation results from Figures 7, 8 are gathered in Table 7.

Discussion
In the following, the simulation results that have been shown in the previous section will be critically discussed and compared to FIGURE 8 | Simulation results for differently oriented PST crystals at the temperatures and stress levels from the experimental results in Parthasarathy et al. (2000). Symbols: experimentally determined steady-state strain rates at the corresponding strain, that is, at the end of primary creep stage. Solid lines: simulations. the experimental findings from Parthasarathy et al. (2000) and other literature.

Creep Simulations
The comparison of the creep simulations to the experimental results from Parthasarathy et al. (2000) allows two observations: (1) The presented model is capable of reproducing the steadystate creep rates of differently oriented PST crystals reasonably well as a function of stress, loading angle, and temperature for the identified set of model parameters (see Figure 7).
(2) For certain load cases and orientations, the primary creep strain is significantly overestimated by the model (see Figure 8).
Having a closer look on the load cases/orientations for which the model overestimates the primary creep strain (i.e., the strain at which the saturation state is reached), it becomes apparent that a significant part of the estimated primary creep strain accumulates within the 1 s in which the load is applied/ increased in the simulations. In all cases for which the primary creep strain is overestimated, the applied stress exceeds the current yield strength in the model; that is, the resolved shear stress on at least one slip/twinning system is FIGURE 9 | Simulation results for differently oriented PST crystals at the temperatures and stress levels from the experimental results in Parthasarathy et al. (2000). Second simulation series, that is, the creep strain that accumulates during the 1 s loading periods is subtracted. Symbols: experimentally determined steady-state strain rates at the corresponding strain, that is, at the end of primary creep stage. Solid lines: simulations.
Frontiers in Materials | www.frontiersin.org February 2021 | Volume 7 | Article 581187 higher than its current strength. Thus, the overestimated primary creep strain is closely related to the occurrence of instantaneous plasticity in the model. In order to identify/quantify the contribution of the instantaneous strain to the primary creep strain, an additional series of simulations is set up. For these simulations, the model set up as well as all identified model parameters remain unchanged. Only the strain that accumulates during loading, that is, during the 1 s linear load increase, is subtracted when plotting the simulation results. The simulation times at which the load is increased are adjusted accordingly so that the load increase is again taking place once the primary creep strain from the experiments in Parthasarathy et al. (2000) has been reached.
The corresponding strain rate vs. strain and strain rate vs. time plots are shown in Figure 9. The secondary creep rates of this second series of simulations essentially remain the same as in the first series of simulations (cf. Figure 7) and are therefore not plotted here again.
By the results in Figure 9, that is, by subtracting the instantaneous strain that accumulates during the time span in which the load is applied, most of the differences between the primary creep strain in the simulations and the experiments can be rationalized without a negative effect on the predicted steadystate creep rates and the simulation results for the load cases/ orientations for which the primary creep strain was correctly predicted in the first place.
Obviously, the occurrence of instantaneous plasticity in the model could also be reduced/suppressed by increasing the initial yield strength (and thus the creep stress to yield strength ratio) or by increasing the work hardening rate (and thus the creep stress to current strength ratio). The former approach is, however, not reasonable as the model has been calibrated to match the yield strengths of the PST crystals tested in the same study (i.e., Parthasarathy et al. (2000)) at the temperatures that have also been used for the creep experiments (see Figure 6). The (athermal) work hardening behavior on the other hand has been calibrated against results from a different experimental study (see Supplementary  Material). However, the drastic increase of the work hardening rate that would be necessary to sufficiently reduce the instantaneous plastic response for the load cases in which the model overestimates the primary creep strain would in turn negatively affect the simulation results for the majority of load cases for which the model estimates the primary creep strains reasonably well.
The fact that the discrepancy between the experimentally determined and the simulated primary creep strain can largely be eliminated by subtracting the instantaneous strain does suggest that there is systematic cause.
Although not frequently discussed in the literature, much of the information on the instantaneous deformation is probably lost in typical creep experiments due to apparative limitations (Estrin and Mecking, 1984). First, the instantaneous strain during loading accumulates with a strain rate that is several orders of magnitude higher than the viscoplastic creep rates and is thus not detectable with equipment that is optimized for measuring creep strain rates. This does, however, not explain why the instantaneous strain might not be captured in creep experiments. Instead, this may result from the typical loading procedure. The disturbance of the sensitive strain measurement during manually applying dead loads impedes a precise determination of the strain during loading. In the context of the long duration of creep measurements, it is thus a pragmatic approach to start the measurement only after the load has been applied, consequently losing all information on the instantaneous deformation.
The instantaneous strain that accumulates during the load application is the sum of elastic strain and instantaneous plastic strain. While losing the information of the elastic part of the instantaneous strain is justifiable because the elastic strains remain small and can be estimated analytically, not detecting potentially occurring instantaneous plastic strains can lead to dramatic misinterpretation of the experimental results and in worst case leads to a drastic overestimation of the time that is needed to reach the admissible creep strain for a given stress. As many experimental studies (including the study discussed here (i.e., Parthasarathy et al. (2000))) did, actually, not exclusively apply creep stresses below the yield point, this hypothesis even rationalizes the frequent occurrence of mechanical twins in creep experiments Appel et al. (2011), Wegmann et al. (2000 as twins do naturally evolve when the load exceeds the yield strength (i.e., in the instantaneous plastic deformation regime).
In the end, this hypothesis cannot be unambiguously proven without further experimental investigation. However, the occurrence of (undetected) instantaneous plastic strains would rationalize many of the observed discrepancies and is supported by the fact that even the (analytically determined) elastic strain due to a load increase would in many cases be higher than the primary creep strains for the reported load cases in Parthasarathy et al. (2000).

Activation Energy for Creep
Numerous experimental studies reported the activation energy for creep in polycolony fully lamellar TiAl, however, with a significant scatter between 320 kJ/mol < Q c < 450 kJ/mol Appel et al. (2011). In the modeling part of Zhang and Deevi (2003), it has been shown that the temperature dependence of the creep behavior of different experiments with polycolony fully lamellar TiAl can be reproduced for an activation energy of Q c 375 kJ/mol which is well in the middle of the reported values. While confirming the findings for polycolony TiAl by reporting an activation energy of Q c 382 kJ/mol, the authors of Parthasarathy et al. (2000) found that for PST crystals, the activation energy varies with the loading angle from Q c 532 kJ(/mol) at 0°over Q c 398 kJ/mol at 45°to Q c 432 kJ/mol at 90°. In combination with the anisotropy of the creep resistance which has consistently been reported to be lowest for PST crystals under loading angles of 45°Asai et al. (2002), Wegmann et al. (2000), Parthasarathy et al. (2000), Marketz et al. (2003), these (apparent) activation energies give a sound However, these activation energies are determined from Arrhenius graphs, that is, by plotting the minimum/steadystate creep rate against the (inverse) temperature. Thus, the reported activation energies can be called apparent in the sense that they represent the macroscopic result of the thermally activated micromechanisms that lead to creep. The plotted macroscopic minimum/steady-state creep rates are, however, related to the creep shear rates on the slip systems via the Schmid/Taylor factor so that the slope in the Arrhenius plot (and thus the activation energy) is different when resolved to the single slip systems. Thus, the experimentally determined apparent activation energies and especially their apparent anisotropy for PST crystals may be rationalized by taking the crystallography into account. While this is obviously not easily done or even impossible in experiments, the simulation results presented here can help to relate the apparent macroscopic activation energies to their slip system resolved representation.
Therefore, the Arrhenius term _ c visc 0 exp − Q c k B θ from Eq. 33 was purposely treated as a single parameter and separately identified for longitudinal, mixed, and transversal slip systems during the model calibration (cf. Supplementary Material).
A detailed discussion of the dislocation mobilization mechanisms that are reflected in the activation energy for creep of PST crystals is naturally beyond the scope of continuum mechanics. However, plotting the identified _ c visc 0 exp − Qc kBθ values for longitudinal, mixed, and transversal slip systems over the temperature (see Figure 10) still allows the following phenomenological observations: • All identified values _  Figure 10 can be fitted using activation energies between 375 kJ/mol and 425 kJ/mol is in good agreement with experimental findings.
As the predominantly activated deformation systems in PST crystals vary with the loading angle, it is possible to phenomenologically relate the identified viscoplastic model parameters (and thus the tendencies in the activation energies) for longitudinal, mixed, and transversal slip systems to the observed creep behavior of PST crystals.
For a loading angle of 45°, mainly longitudinal deformation systems become active (even under creep conditions) as they have the highest Schmid factors and the lowest strength. For a loading angle of 90°, only transversal deformation systems have Schmid factors ≠ 0 so that longitudinal and mixed slip systems do not contribute to the plastic deformation. For 0°, longitudinal deformation systems have a Schmid factor of 0. Thus, mainly mixed and transversal deformation systems contribute to plastic deformation.
At first glance, the identified _ c visc 0 exp − Q c k B θ values do not reflect the experimentally observed anisotropy in the creep resistance as longitudinal and mixed slip systems can be described by nearly the same parameters although the creep resistance of PST crystals under loading angles for which these deformation systems dominate (i.e., 45°and 0°) differs significantly. The same holds for the _ c visc 0 exp − Qc kBθ values of transversal slip systems which correspond to low activation energy and thus do not fit the high creep resistance that has been observed for PST crystals under loading angles of 90°.
However, in the current model, the viscoplastic creep behavior has been described by Eq. 33, that is, in a slip system resolved manner in which the viscoplastic shear rate on a slip system is determined by normalizing the classical creep model (i.e., Eq. 1) by the slip system strength. Thus, the creep resistance of longitudinal, mixed, and transversal deformation systems is in fact not given solely by the activation energy but also incorporates the slip system strengths and is related to the experimentally observed values via the Schmid factors as already mentioned above.
Thus, the identified _ c visc 0 exp − Q c k B θ values do not contradict the experimental findings but instead represent the slip system resolved, micromechanical interpretation of the creep law in Eq. 1. In fact, the activation energies for which the _ c visc 0 exp − Q c k B θ can be fitted (cf. Figure 10) may even be interpreted in the context of the yield stress temperature anomaly shown in Figure 6. While the yield stress for a loading angle of 90°(and correspondingly the critical resolved shear stresses of the predominantly activated transversal deformation systems) monotonically decreases, the yield stresses for 0°and 45°( and correspondingly the critical resolved shear stresses for mixed and longitudinal deformation systems) show a plateau or even increase in the range of tested temperatures (see Figure 6). Thus, longitudinal and mixed slip systems can be said to show a temperature anomaly in the indicated temperature range (700°C-815°C). Consequently, the activation energy that successfully fits the _ c visc 0 exp − Qc kBθ values of both the longitudinal and mixed slip systems (cf. Figure 10) can be interpreted to be related to the mechanisms that cause the temperature anomaly. In the line of these arguments, it is reasonable that the activation energy that fits the _ c visc 0 exp − Q c k B θ values for transversal slip systems is lower because the experiments in which transversal deformation systems are predominantly activated (i.e., for a loading angle of 90°) show no anomaly in the yield stress between 700°C and 815°C.

The Role of Static Recovery
As opposed to other alloys for which static recovery can often be neglected Estrin and Mecking (1984), TiAl alloys have been found to exhibit a fast static recovery behavior Paul and Appel (2003), Lindemann et al. (2006), Schnabel et al. (2019), Appel et al. (1999). In Paul and Appel (2003), Appel et al. (1999), the fast static recovery kinetics of TiAl has been associated to the annealing-out of dipoles and debris. In fact, it has been shown in Lindemann et al. (2006) that these defects even recover during constant strain rate deformation at sufficiently high temperatures. Thus, the recovery of these thermally unstable defects can be considered as instantaneous in the context of creep deformation so that these defects do not contribute to the work hardening during primary creep. A possibility to model the additional work hardening that is provided by dipoles (and consequently their separate, fast recovery) is introduced via a separate dislocation evolution term in Ma et al. (2006). However, in the present work, dipoles are not modeled separately as they do not contribute to the work hardening in primary creep. Still, the static recovery behavior is not neglected here as the thermally more stable defect structures may also recover although with a significantly slower recovery kinetics.
Interestingly, the recovery term R stat exp − Q R k B θ had to be chosen ≠ 0 only for transversal slip systems during the model calibration (see Supplementary Material). This shows that the static recovery behavior of the dislocation network is in fact negligible for longitudinal and mixed slip systems. A possible explanation for the necessity to choose R stat exp − QR k B θ ≠ 0 for transversal slip systems may be found in the dislocation pile up against the lamella boundaries which results in a high density of immobile dislocations on parallel slip planes potentially facilitating recombination by climb. However, this surely needs further investigation.

CONCLUSION AND OUTLOOK
The results that are shown in the present work lead to the following conclusions regarding the creep behavior of PST crystals and the presented model: • The climb-assisted glide model from this work is able to reproduce the creep experiments with differently oriented PST crystals from Parthasarathy et al. (2000) reasonably well for a viscoplastic stress exponent of n visco 5. • The results of the presented numerical study indicate that the instantaneous (plastic) strain that accumulates during application of the creep load may not be captured in creep experiments due to apparative limitations. • The apparent anisotropy of the activation energy for creep that has been observed in experiments with differently oriented PST crystals changes its characteristics when resolving the shear rates to the predominantly activated deformation systems and normalizing the results to the different strengths of deformation systems. • Differences in the identified activation energies for creep of deformation systems from different morphological classes can potentially be rationalized in the context of the yield stress temperature anomaly of PST crystals. • The role of the fast static recovery kinetics in creep of TiAl alloys needs further experimental investigation in order to backup or-if necessary-refine the model.
In consequence, the present study shows that interpreting/ modeling the steady-state/minimum creep rates of PST crystals using models in the form of Eq. 1 leads to misconceptions in both the stress exponents n and the activation energies Q c .
In its present form, the derived model is capable of reproducing the transition from primary to secondary creep stage for PST crystals of TiAl 12 in a microstructure sensitive way. However, the modeled processes of dislocation creep only describe the bulk creep behavior. For very low stresses, the present model will thus underestimate the secondary creep rate as under this conditions interface glide may become the dominant creep mechanism. In addition, the steady state is not reached for certain load cases in fully lamellar TiAl due to the early onset of damage/tertiary creep. While the presented evolution model for the dislocation density is not invalidated in the tertiary creep stage, the creep rate for the respective deformation state will still be underestimated as damage is not yet incorporated in the model. Thus, in a next step, the here presented model may be extended to account for damage in the bulk of the lamellae. The interface related effects, that is, interface glide and interface damage, cannot be reasonably captured in the context of crystal plasticity. In order to incorporate these interface related effects, the cohesive element formulation presented in Scheider et al. (2020) can be introduced to the RVE of the lamellar microstructure.

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

AUTHOR CONTRIBUTIONS
JS developed and implemented the model, conducted the simulations, and analyzed the resultant data. JS and IS discussed necessary model assumptions, simulation results, and the resultant conclusions. JS prepared the manuscript. IS did proof reading.

FUNDING
This work received internal funding by the VirMat project of the Helmholtz Association of German Research Centres and the IDEA project Modeling the effect of interfaces on creep of lamellar TiAl from a multiscale perspective of the Helmholtz-Zentrum Geesthacht.