Crystal Plasticity Simulation of Magnesium and its Alloys: A Review of Recent Advances

: Slip and extension twinning are the dominant deformation mechanisms in Magnesium (Mg) and its alloys. Crystal plasticity is a powerful tool to study these deformation mechanisms. Different schemes have incorporated crystal plasticity models to capture different properties, which vary from the simple homogenization Taylor model to the full ‐ scale crystal plasticity finite element model. In the current study, a review of works available in the literature that addresses different properties of Mg and its alloys using crystal plasticity modes is presented. In addition to slip and twinning, detwinning is another deformation mechanism that is activated in Mg and its alloys. The different models that capture detwinning will also be addressed here. Finally, the recent experi ‐ mental frameworks, such as in ‐ situ neutron diffraction, 3D high energy synchrotron X ‐ ray tech ‐ niques, and digital image correlation under scanning electron microscopy (SEM ‐ DIC), which are incorporated along crystal plasticity models to investigate the properties of Mg and its alloys, are addressed. Future research directions towards improving the deformation response of Mg and its alloys are identified, which can lead to increased deployment of the lightest structural metal in en ‐ gineering applications.


Introduction
Mg alloys can significantly impact many industrial applications such as automobile, aerospace, and transportation, not to mention the energy sector, by effectively reducing the weight of the designed system. This can lead to a significant improvement in fuel economy and reducing emissions [1]. However, the current understanding of commercial Mg alloys and their physical metallurgy is less developed compared to other established structural metals and alloys. To overcome many of the challenges in using Mg alloys in these applications, different aspects of these alloys, including strength, formability, and fatigue resistance, should be enhanced. To fulfil any of these objectives, the underlying deformation mechanisms of Mg alloys should be thoroughly understood.
Two important deformation mechanisms of Mg alloys are plastic slip and extension twinning. Mg alloys have Hexagonal Closest Packed (HCP) crystal structure. In the case of slip, different modes can be activated, including basal 〈a〉 0001 〈112 0〉 , prismatic 〈a〉 101 0 〈112 0〉 , pyramidal 〈a〉 101 1 〈112 0〉 , and pyramidal 〈c a〉 1 1 22 〈1 1 23〉 . In the case of unalloyed Mg polycrystals, the basal mode has a very low critical resolved shear stress (CRSS) and becomes the dominant plasticity deformation mechanism in most of the loading conditions. In the case of Mg alloys, although the basal mode CRSS increases, it still becomes the lowest among the slip modes. Although the basal mode CRSS is very low, its corresponding resolved shear stress disappears during the tensile deformation along the c-axis. The extension twinning 101 2 〈1 011〉 becomes an important deformation mechanism of Mg alloys in this loading condition [2][3][4][5][6][7][8][9]. Similarly, the extension twinning is the controlling deformation mechanism during the compression loading parallel to the basal plane. The highly asymmetric response of Mg alloys for some special orientations can be attributed to the activation of extension twinning, which has a complex morphology and micromechanics [3,[10][11][12][13][14][15][16]. Although consideration of slip and twinning mechanisms is enough to describe the response of Mg alloys in most of the loading types, which are monotonic, this is not the case for the loading types, which involve unloading, such as cyclic loading. These loadings are especially important in many applications, such as fatigue. In the case of cyclic loading parallel to the basal plane, twinning occurs during compression loading, which results in the 86.3° reorientation of the basal pole. During subsequent reversed cyclic loading, i.e., tensile loading parallel to the basal plane, the size of the twinned regions shrinks, and in some cases, the twin fully disappears, a phenomenon that is commonly known as detwinning [7,[17][18][19][20][21][22]. This twinning and detwinning occur alternately during the cyclic loading of Mg alloys.
In-situ experiments have a key role in unraveling the underlying deformation mechanisms of Mg alloys. Three in-situ experimental schemes, that were used to investigate the Mg alloy response were Digital Image Correlation (DIC), Synchrotron X-ray techniques, and neutron diffraction. In-situ DIC is another experimental technique that is used to address the deformation mechanisms of Mg alloys. Different aspects of Mg alloys have been studied by in-situ DIC, including the effect of twin nucleation on strain field [23], twin morphology [24], and the effect of aging on the accumulation of microscale plasticity [25]. Two common in-situ Synchrotron X-ray techniques of three-dimensional X-ray diffraction (3DXRD) and high-energy diffraction microscopy (HEDM) were incorporated in the investigation of Mg alloy behavior. Aydıner et al. [26] investigated the response of the AZ31 Mg alloy using the 3DXRD technique and analyzed the evolution of stress during deformation. It was observed that the stress state inside the twinned region was different from the untwinned grain. Other studies of Mg alloys have incorporated the in-situ 3DXRD and HEDM to analyze the micromechanics of twinning [27], the effect of age hardening on the deformation behavior [28], evaluate the CRSS for deformation modes [29,30], twin growth [31], and detwinning [7,22,32]. Using in-situ neutron diffraction, it was shown that the twinning in extruded Mg-7.7 at.% Al alloy was manifested as the change in the intensity of diffraction peaks during loading [33]. The in-situ neutron diffraction experiment on the Mg alloy showed that the stress relaxation occurred in the twinned region compared to the untwinned region of grain [34]. In-situ neutron experiments were also used to address the deformation detwinning in Mg alloys [35][36][37][38][39]. Although these insitu experiments can provide invaluable information regarding the underlying deformation mechanisms of Mg alloys, they are very expensive, time-consuming, and complex. Potential simulation frameworks benefit the research community to complement these experimental observations. Namakian and Voyiadjis [40] incorporated atomistic simulations and proposed a new mechanism for 101 2 〈1 011〉 twinning in which it was postulated that partial stacking faults (PSFs) play a key role in the formation of 101 2 twins. A detailed crystallographic study was conducted for different HCP metals. Accordingly, PSFs were created by 1 ⁄ 3 〈101 0〉 displacement vectors on every other basal plane and created a faulting plane on 101 2 plane. Figure 1 shows the geometrical properties associated with the transformation of the untwinned parent into the twinned child unit cell. The model defines the properties of HCP crystals in a way that atoms can undergo spontaneous cooperative movements within these crystals. Accordingly, considerable stress relaxation occurs due to twin formation. This can be described using a deformation gradient, which shows that the macroscopic effect of the current mechanism is a simple shear process.
Namakian et al. [41] incorporated atomistic simulations to model both deformation mechanisms of slip and twinning in single crystal Mg along with the interrelationships on the loose π_1L and dense π_1D first-order crystallographic planes. The generalized stacking fault energy (GSFE) analysis was used to investigate the slip mechanisms to study the 〈c + a〉 dislocation's core structure, dissociation mechanism, and mobility. Namakian et al. [41] stated that the screw component of a dissociated pyramidal-I 〈a〉 dislocation had a key role in compression twinning (CTW) nucleation. They also showed that CTW can grow by activating pyramidal-I 〈a〉 slip on the preexisting twin boundaries. Figure 2 shows the minimum energy path (MEP) of four deformation mechanisms associated with π_1 plane, including homogeneous CTW nucleation, heterogeneous nucleation mechanism of pyramidal-I 〈c + a〉 slip, i.e., 1/3 1 1 23 101 1 slip to CTW 1 012 101 1 (or zonal twinning mechanism), heterogeneous nucleation mechanism of pyramidal-I 〈a〉 slip 1/3 1 21 0 101 1 to CTW 101 2 101 1 _ , and pyramidal-I 〈c + a〉 slip 1/ 3 1 1 23 101 1 using Liu et al. [42] embedded-atom method (EAM) and Wu et al. [43] modified embedded-atom method (MEAM) interatomic potentials. Namakian et al. [41] showed the importance of the interrelationships of non-basal slips and CTW on π_1 plane. They concluded that these interrelationships should be incorporated in simulation frameworks with larger length scales, including crystal plasticity, to accurately capture the mechanical response of Mg. The schematics of the transformation of the untwinned parent unit cell (blue segments and red and blue atoms) into the twinned children unit cell (cyan segments and magenta and cyan atoms for the first two layers close to the twin boundary) during 101 2 1 011 twinning (After Namakian and Voyiadjis [40]).

Figure 2.
Minimum energy paths of four deformation mechanisms associated with π_1 plane, i.e., homogeneous CTW nucleation, heterogeneous nucleation mechanism of pyramidal-I 〈c+a〉 slip, i.e., 1/3 [ The deformation mechanisms and response of metallic systems and alloys have been studied using frameworks at different length scales spanning from density functional theory (DFT) [44,45], atomistic simulations [11,40,41,[46][47][48][49][50][51][52][53][54][55][56], discrete dislocation dynamics [57,58], and continuum mechanics [59]. However, all these schemes except continuum mechanics suffer from length and time scale gaps to model a real-size sample or experiment. Crystal plasticity (CP) is a very powerful continuum mechanics framework to model the response of Mg alloys with both slip and twinning deformations. Various CP models have been developed to simulate slip and deformation in HCP polycrystals [60][61][62][63][64][65][66][67][68][69]. A combination of CP and finite elements, which is commonly known as CP finite element (CPFE), benefits the advantages of both frameworks. Two main categories of CPFE frameworks were developed to capture both slip and twinning deformations, which were developed by Staroselsky and Anand [63] using rate-independent formulation and Kalidindi [62] using rate-dependent formulation. In the case of Mg and its alloys, the rate-independent framework developed by Staroselsky and Anand [63] was used to model the slip along with twinning [66,68,69]. Kalidindi [62] introduced a modified plastic velocity gradient tensor by incorporating twinning. Various researchers have included twinning into their CPFE framework following the formulations introduced by Kalidindi [62], such as Abdolvand and his coworkers [65,[70][71][72] and Zhang and Joshi [67]. Qiao et al. [73] included the stress relaxation due to twinning into the CPFE formulation, a CPFE framework that can capture the stress relaxation effect for twinning. Hama et al. [74] incorporated twinning into the CPFE framework to model the anisotropy of Mg alloy sheets subjected to a twostep loading. Prasad et al. [75] simulated the ductile fracture in pure Mg by adding the twinning to plane strain CPFE.
The CP models were used in combination with other schemes to model twinning. One of these candidates was a viscoplastic self-consistent (VPSC), which was first used by Lebensohn and Tomé [76] to model the twinning that occurred in zirconium alloys. Agnew et al. [77] used the VPSC to model the slip and twinning in Mg alloys. Beyerlein and Tomé [12] simulated the twinning in pure zirconium using VPSC by introducing a probabilistic model for the twin nucleation. Kumar et al. [78] used a full-field elasto-viscoplastic model based on fast Fourier transformation (FFT), which included twinning to model the response of Mg. Lévesque et al. [79] introduced twinning into a Taylor-type CP model to simulate the response of AM30 and AZ31B Mg alloys. The multiscale framework of VPSC-FE was incorporated to simulate the slip and twinning deformations [80][81][82]. Accordingly, the response of each FE material point was defined using an RVE, which was modeled by VPSC. A similar multiscale scheme has been used by Ardeljan et al. [83] and Feather et al. [84] to simulate the response of Mg alloys in which an RVE modeled by a Taylor-type CP defines the constitutive model of an FE material point. A CP model has also been coupled with a phase field method to better describe the twin morphology in HCP metals [85,86].
Most of the studies have focused on the response of the Mg alloys during monotonic loading, which only requires the deformation slip and twinning. In the case of more complex strain paths such as cyclic loading, however, detwinning should also be defined in the CP model. A meso-scale composite grain (CG) model was the first detwinning model, which was included in the VPSC framework [87,88]. An empirical model defines the twin nucleation and propagation in the CG model. Guillemer et al. [89] used a simple phenomenological detwinning model along a self-consistent model to capture the cyclic behavior of extruded Mg. The first physically-based model, which includes both mechanisms, was the Twinning-Detwinning (TDT) model developed by Wang and his coworkers [90][91][92], which was used along the EVPSC framework. Qiao et al. [93] modeled the cyclic response of the ZK60A Mg alloy using the TDT/EVPSC scheme. The twinning and detwinning can be effectively handled by EVPSC framework by introducing new grains as twin nucleation and removing that twinned grain for complete detwinning, while the change in the volume of the twinned grain manifests the twin growth or shrinkage. However, in the case of studying the local variable information in a polycrystal, one cannot use the EVPSC model. One way to obtain the local information is using the CPFE framework to capture twinning and detwinning mechanisms. However, most of the CPFE models introduced to capture the twinning and detwinning mechanisms are phenomenological and neglect some of the main elements of the deformation twinning and detwinning. The biggest challenge in CPFE is that the introduction of twin nucleation, growth, shrinkage, and detwinning is not as straightforward as EVPSC. In the case of Mg and its alloys, the CPFE method has been used recently to model the cyclic behavior, including twinning and detwinning [22,[94][95][96][97]. A material point is treated in these models as two states of twinned or not twinned, and the reorientation criteria are commonly based on the most active twinning system, which is called the predominant twinning reorientation (PTR) scheme [61]. Before the reorientation, the stress inside the twinned region does not contribute to the stress state of the material point, while after reorientation, the stress state inside the parent grain is neglected. However, the stress state of the twinned region and untwinned region are not similar, which was shown using in-situ experiments [26,34]. Yaghoobi et al. [9] introduced a multiscale CPFE framework in which the deformation twinning and detwinning were captured using a physically-based TDT model proposed by Wang and his coworkers [90][91][92]. They used a Taylor-type homogenization scheme to model both twinned and untwinned regions, coexisting at a material point.
Traditional crystal plasticity models were developed largely without a connection to grain size and shape effects. In Magnesium alloys, it has been observed that the grain size strengthening follows the Hall-Petch effect [98]. The incorporation of grain size effect into constitutive models for single slip began in 1962 with Armstrong et al. [99], who modified the Hall-Petch equation to correspond to the flow stress on a slip system (the "micro-Hall-Petch relation"). The interrelationship between grain size and texture was not considered until 1983, when Weng [100] employed the mean grain size in the equation for slip system resistance through the micro-Hall-Petch relation. Sun and Sundararaghavan [101] presented a crystal plasticity model for grain size and shape effects where the slip length of each slip system at a material point was considered in a micro-Hall-Petch strengthening term for each slip system. Recent HREBSD experiments [102] in Mg-4Al alloy have revealed the interaction of slip systems and grain boundaries is dependent on not only the slip length but also strongly on the angle between the incoming and emerging slip planes at a grain boundary. In the future, there is a need to better incorporate such insights in crystal plasticity.
In the current work, a review of works available in the literature, which addresses different properties of Mg and its alloys using crystal plasticity models, is presented. In addition to slip and twinning, detwinning is another deformation mechanism, which is activated in Mg and its alloys. The different models that capture detwinning will also be addressed here. Finally, the recent experimental frameworks, such as in-situ neutron diffraction, 3D high energy synchrotron X-ray techniques, and digital image correlation under scanning electron microscopy (SEM-DIC), are incorporated along crystal plasticity models to investigate the properties of Mg and its alloys are addressed. Future research directions towards improving the deformation response of Mg and its alloys are identified, which can lead to increased deployment of the lightest structural metal in engineering applications.

Overview
The CP theory is developed according to the assumption that the plasticity in metals and alloys occurs as a result of slip-on prescribed slip systems. The finite deformation continuum mechanics is usually incorporated to describe the formulation. The deformation gradient tensor is decomposed as below: (1) where and are the elastic and plastic deformation gradients, respectively. Figure  3 shows the deformation described in Equation (1). Elastic distortion of the crystal lattice and plastic slip are two deformation mechanisms, which sustain the applied deformation in this formulation. Accordingly, the macroscopic velocity gradient can be decomposed as below: (2) where and are the elastic and plastic velocity gradients, respectively. The novel relation of CP is to link the macroscopic velocity gradient to micro deformation as below: ( 3) where is the shearing rate on slip system , is the number of slip systems, and is the Schmid tensor for the slip system , which can be defined as follows: where unit vectors and denote the slip direction and slip plane normal, respectively, in the deformed configuration. The resolved shear stress on the slip system can be obtained as follows: where is the Cauchy stress tensor and ⋅ operator denotes the standard inner product of tensors. The shearing rate on the slip system can be calculated depending on if the CP formulation is rate-independent or rate-dependent. In the case of a rate-independent model, the yield surface is defined for each system, i.e., , as follows: where is the slip resistance for slip system α. The values of is then obtained using the yield surface and considering the hardening model, which describes the evolution of slip resistance.
In the case of rate-dependent formulation, the evolution of shearing rate on the slip system can be defined as follows: where and are the material parameters denoting the reference shearing rate and rate sensitivity of the material. One should note that the kinematic hardening is neglected in both Equations (6) and (7).

Twinning
In the case of Mg and its alloys, extension twinning 101 2 〈1 011〉 is very important, along with the plastic slip ( Figure 4). The simplest modification to include twinning into the CP formulation is to consider twinning as pseudo-slip systems (See, e.g., Staroselsky and Anand [64]), and their contribution to the plastic velocity gradient tensor can be described as follows: (8) where is the number of twinning systems. Figure 5 shows the kinetics of slip and twinning defined in Equation (8). The relation between the plastic slip of twinning systems and their corresponding twin volume fraction can be described as follows: (9) where is the rate of change in twin volume fraction of twin pseudo-slip system and is the characteristic twin shear strain, which defines the amount of shear associated with twinning. The value of depends on the / ratio, which can be defined as follows [104]: where and are depicted in Figure 4. In the case of Mg and its alloys, the 0.129.  Although the twin pseudo-slip system's contribution to plastic velocity gradient tensor is similar to that of the slip systems, they are not precisely treated similarly due to the polar nature of twinning. In other words, only the tensile component of stress along the c-axis can trigger the extension twinning. This is reflected in the CP formulation for the twinning pseudo-slip systems as follows: Kalidindi [62] further enhanced the effect of twinning on crystal plasticity formulation. The work considers the contribution of stress in both untwinned and twinned region as follows: where is the reoriented volume fraction of grain according to the twin system . Kalidindi [62] then presented the modified macroscopic plastic velocity gradient tensor , which includes the contribution of multiple twinned systems as below: where , , and denote the Schmid tensors for the slip systems in parent grain, twin systems in parent grain, and slip systems in twinned children, respectively. denotes the number of slip systems in the twinned region. Equation (10) describes three mechanisms contributing to the plastic velocity gradient . The first term defines the contribution of slip inside the parent grain, the second term defines the contribution of twin systems in parent grain, and the third term is the summation of the contributions of slip systems in all twinned children.
In the twinning models, there is a key part that should be defined as the reorientation. Figure 6 shows the crystallography of the deformation twinning, which includes the untwinned region (parent grain) and twinned region (child). While Equation (9) defines the amount of twin volume required for a certain amount of shear strain, it does not define the reorientation in the model. Different criteria have been introduced for the reorientation. Van Houtte [60] was the first researcher to propose a method for reorientation. In this method, reorientation of the grain is decided depending on the volume fractions calculated at every incremental step along with using a random criterion. Tomé et al. [61] modified the model presented by Van Houtte [60] using the predominant twin reorientation scheme (PTR) in which the grain is reoriented according to its predominant twin variant as the integration of twin volume throughout the time satisfy specific criteria.

Stress Relaxation
The in-situ neutron diffraction experiment on the Mg alloy showed that the stress relaxation occurs in the twinned region compared to the untwinned region of grain [34]. Different crystal plasticity models were incorporated to capture the stress relaxation during twinning. The simplest model is to consider larger CRSSs for twin nucleation than growth [12,93]. Wu et al. [105] presented a comprehensive physically-based twin model implemented in the EVPSC framework. They presented this twin model, which was called 'TNPG,' to capture twin nucleation, propagation, and growth for magnesium alloys (Figure 7). In Figure 7a, the twinning process initiates with twin nucleation. It is assumed that the grain boundary is the preferred nucleation site. As the twin propagates according to Figure 7b, stress relaxation occurs. Different notation is used throughout the current review for twinning resistance as compared to the one used by Wu et al. [105], which used ̂ (Figure 7). The twinning resistance decreases during the propagation phase until it reaches its minimum value of when the twin volume is . Afterward, the twinning resistance increases due to the hardening occurring in the twin thickening process. Wu et al. [105] captured the stress relaxation by modifying the twin systems resistances as follows: where Γ is the accumulated shear of twinning system, which is calculated after its corresponding twin volume surpass (see Figure 7). They used this model to capture the stress-strain response of AZ31B Mg alloy obtained by Lou et al. [106]. They also compared the predicted twin volume versus the experimental results. Figure 8 shows that the model developed by Wu et al. [105] can successfully capture both the stress-strain and twin volume during uniaxial compression and tension. Qiao et al. [73] reviewed the avail-able scheme to model stress relaxation in the CPFE framework. They used the model developed by Wu et al. [105] for EVPSC and incorporated it in the CPFE model to capture the stress relaxation in Mg single crystal during twinning.

Detwinning
In the case of more complex strain paths, the detwinning mechanism should also be included in addition to twinning and slip modes. Various models address the detwinning along twinning and slip modes using crystal plasticity [22,[87][88][89][94][95][96][97]. Here, a physically-based Twinning-Detwinning (TDT) model is elaborated, which was developed by Wang and his coworkers [90][91][92] for the EVPSC framework and extended later on for the CPFE framework [9]. In this model, the twinning and detwinning mechanisms can be divided into 4 major operations as follows ( Figure 9):  Operations A: Twin nucleation and growth due to parent grain reduction.  The polarity of simple twinning model should be modified for these operations. In the parent grain, the crystallographic systems include slip systems, twin systems for Operation A, and twin systems for Operation C. In the case of the twin child, the crystallographic systems include slip systems, one twin system for Operation D, and one twin system for Operation B. Operation A can get activated in parent grain for when the resolved shear in the twin system is larger than the corresponding slip resistance, i.e., . Operation B is activated in the twin child, which is nucleated by the activation of twin variant, when 0 and .
Operation C initiates in the parent grain for the twin variant when 0 and | | . Operation D initiates inside the twin child, which is nucleated by the activation of twin variant, when . In the TDT/EVPSC model developed by Wang and his coworkers [90][91][92], the twinned child is manifested as a new grain, and the Operations A-D is reflected by the change in the volume of the parent and twin grains. Wang et al. [90] used the developed TDT model to capture the cyclic response of AZ31B Mg alloy for complex loading paths. They compared the model prediction versus the experimental data provided by Lou et al. [106] (Figure 10). They also showed the variation of twin volume and how the deformation twinning and detwinning governs the twin volume and stress-strain curve shape. Furthermore, they compared the prediction of TDT/EVPSC model versus the CG results reported by Proust et al. [88] for the experimental results of AZ31B Mg alloy subjected to the case of in-plane compression followed by through-thickness compression ( Figure 11). Although the CG model was able to capture some features of twinning-detwinning mechanisms, the TDT/EVPSC model considerably enhances the accuracy of the simulation results.  Unlike the EVPSC framework, the TDT operations are not trivial to apply to a material point in the case of CPFE framework. Yaghoobi et al. [9] used a multiscale framework in which the untwinned and twinned regions were considered within a sub-scale model at a material point and the stresses homogenized using a Taylor-type scheme ( Figure 12). This allows multiple variants to coexist at a material point, which avoids the need for the selection of predominant variants as in the PTR models. They implemented the model in an open-source CPFE software PRISMS-Plasticity [103]. The twin volume evolution of a material point for the twin variant is obtained as below: here and are the shear rate inside the parent corresponds to operations A and C, and ̅ and ̅ are the shear rates within the child nucleated due to the activation of the twin system. Yaghoobi et al. [9] calibrated the model using the uniaxial experimental data of extruded ZK60A Mg alloy presented by Wu [107] (Figure 13). The predictions of the modeled cyclic response of the ZK60A alloy along the extrusion direction were compared with the experimental data by Wu et al. [17] and Wu [107]. Yaghoobi et al. [9] further studied the cyclic response of ZK60A by comparing the predicted twin variation versus the normalized intensity of the 0002 diffraction peak along the longitudinal direction measured by Wu et al. [35] and Wu [107]. The result showed that not only can the multiscale TDT CPFE framework successfully capture the cyclic stress-strain response of the ZK60A (Figure 14), it can also capture the deformation twinning and detwinning during the cyclic loading ( Figure 15). Furthermore, Yaghoobi et al. [9] showed the evolution of predicted basal 0001 pole figures at strains of 1.2% (Figure 16). It was shown that the deformation twinning resulted in the reorientation of the basal pole, which intensifies in the basal 0002 peak in the loading direction at max compressive strains of 1.2% (Figure 16 (b) and (d)). On the other hand, the detwinning mechanism removed the increased basal 0002 peak intensity at a max tensile strain of 1.2% (Figure 16 (c) and (e)).  [17] and Wu [107] during the cyclic loading along the extrusion direction in ZK60A Mg alloy (After Yaghoobi et al. [9]). Figure 15. The variation of predicted twin volume versus the strain compared to the experimentally measured change in the normalized intensity of the 0002 diffraction peak along the longitudinal direction obtained by Wu et al. [35] and Wu [107] in ZK60A Mg alloy during the cyclic loading along the extrusion direction (After Yaghoobi et al. [9]).

Microscale Deformation Mechanisms
Githens et al. [25] used the in-situ DIC under scanning electron microscopy (SEM-DIC) to study the deformation mechanisms of WE43-T5 Mg alloy with weak basal texture subjected to uniaxial tension and compression along the rolling direction. The SEM-DIC technique provided the full-field strain maps. They first investigated the response of the WE43 alloy during uniaxial tension ( Figure 17). They showed that the heterogeneous pattern of strain does not vary after the yield point. They also studied the strain probability distribution, which confirms that the strain map does not vary after yielding ( 2.91% and 4.86%). However, this pattern is different from the one before yielding ( 1.23%). They simulated the sample using CPFE simulation open-source software of PRISMS-Plasticity [103] and compared the predicted strain map with the experimental observations, which agrees well ( Figure 18).  They not only verified the CPFE framework using the SEM-DIC results, but they also resolved and categorized the strain from individual slip traces for the first time in any Mg alloy, as shown in Figure 18 (a). They used this information to categorize the active slip/twinning modes with respect to their nominal Schmid factors in both SEM-DIC and CPFE results ( Figure 19). The results show that the basal modes sustain the most plastic deformation during uniaxial tension. The non-basal slip modes are less active. The least active is the extension twinning, which is because of the tested sample texture and loading direction. Githens et al. [25] also studied the distribution of active slip/twin modes during the compression loading ( Figure 20). The activity of extension twinning is considerably increased during compression compared to the tensile loading ( Figure 19). In addition, the non-basal slip modes contribution to the deformation is very limited during compression ( Figure 20).

Effects of Heat Treatment
Ganesan et al. [8] incorporated the in-situ DIC along the SEM technique to investigate the effect of heat treatment on the response of WE43 Mg alloy with the weak basal texture. They started with the WE43-T5 sample and applied different heat treatment conditions of the solution treated (ST) and aging times of 15 min, 2 h, 4 h, and 16 h (T6 condition). They investigated the effect of heat treatment on the macro responses of WE43 Mg (Figure 21). The results showed that the T5 condition had the highest yield stress, which can be attributed to its finer grain size. Within the heat-treated samples, the T6 condition had the maximum yield strength, which was due to the role of precipitates. Each of the samples was then subjected to the in-situ uniaxial tension loading along the rolling direction, where the SEM-DIC technique was used to gather the deformation maps. Figure 22 shows the normalized maximum principal strain maps at 3.23% in samples with different aging times. They showed that the strain map of the sample with T6 condition had fewer localized strains. In order to quantitatively investigate their observation, Ganesan et al. [8] studied the strain probability distribution of the tested samples ( Figure 23). They showed that it is ~3.84 more probable to find a localized strain of 4 Average applied strain in underaged samples compared to the T6 condition.  They further investigated the effect of heat treatment using the CPFE simulation. They validated the CPFE framework and showed that the simulation can capture the stress-strain response of samples with different heat treatment conditions ( Figure 24). They further validated the simulation by predicting the strain map of WE43-T6 at different tensile strains (Figure 25). The results showed that not only can the CPFE capture the macro responses, but it can also accurately model the local field maps. After validation of the CPFE framework, Ganesan et al. [8] used the CPFE simulation to investigate the deformation mechanisms involved in the response of T5 and T6 conditions. Accordingly, they compared the contribution of each deformation mode for each of these heat treatment conditions of T5 and T6 during uniaxial tension and compression ( Figure 26). They showed that the heat treatment condition alters the contribution of different deformation modes. A sample with a T5 condition has finer grains, which alter the CRSS of deformation modes. For instance, they obtained the ratio of ⁄ 1.02 in the T5 condition, while this ratio was 0.88 for the T6 condition. They stated that this was the underlying mechanism for more active pyramidal in the T5 condition compared to the activity of prismatic, while pyramidal slip mode was nearly inactive in the T6 condition. Finally, Ganesan et al. [8] used the CPFE simulation results to capture the Hall-Petch constant for different deformation modes in WE43 Mg alloy, which were 4.53, 11.33, 8.9, 11.18, 5.1 MPa mm ⁄ for different deformation modes of basal, prismatic, pyramidal , pyramidal c , and extension twinning, respectively.

Micromechanics of Twinning
Abdolvand et al. [27] used a combination of CPFE and in-situ 3DXRD experiments to investigate deformation twinning in the AZ31B Mg alloy. They applied the uniaxial tension on the sample with the texture favored for twinning during the loading and used the in-situ 3DXRD to study the microstructural evolution. They calibrated the CPFE model to capture the stress-strain response, as shown in Figure 27. They also studied the effect of twinning on texture. Accordingly, they measured the 0002 pole figures at a tensile strain of 0% and 0.4%. In the case of 0%, as shown in Figure 28 (a), basal poles were observed to be aligned almost parallel to the normal direction. However, in the case of 0.4%, as shown in Figure 28 (b), the deformation twinning reorients the 0002 poles, which were observed towards the rolling direction and transverse direction on the outer edges of the pole figure. Abdolvand et al. [27] also predicted the twin volume fraction using the CPFE simulation, which agrees well with the measured twin volume fraction (Figure 28 (c)). They further investigated the CPFE simulation results by studying the variation of predicted twin volume versus misorientation angle ( ) between normal to the basal plane of the grain and the loading direction predicted by CPFE at a tensile strain of 0.4% ( Figure  28 (d)). In the case of grains with 40°, no twinning was predicted by CPFE. Furthermore, the scatter of twin volume for grains with similar showed that global was not the only governing factor of twinning volume as the effect of the local neighborhood and load sharing considerably varied the predicted twin volume fraction.  Another important subject Abdolvand et al. [27] addressed was the stress inside the twinned children and parent grains using the in-situ 3DXRD results ( Figure 29). They selected 15 parent grains along 20 twinned children nucleated within those parent grains. Figure 29 (i) shows different components of the stress presented in the results summary. Different stress components were the normal stress , stress along the c-axis of the untwinned parent , and normal and shear stress according to the twinned child plane, i.e., and , respectively. The mean value of each stress was represented by a line. Abdolvand et al. [27] observed that the mean normal stress and almost followed the applied stress trend. Because of the texture (c-axis is almost aligned with the loading direction), the values of and were close. The values of in untwinned parents and twinned children grains were close, which were in the region of 23 30 MPa, which did not vary considerably during the loading. Interestingly, Abdolvand et al. [27] observed that the shear stress, which was resolved on the twinning plane, i.e., , was considerably different in untwinned parent and twinned children grains. They reported that the average value of in untwinned parent grains was 22 MPa higher than that of the twinned children grains. To further analyze the experimental results, Abdolvand et al. [27] focused on an untwinned parent grain with its two twinned children grains, as shown in Figure 29 (e-h), which were nucleated due to the activation of the twin variants 1 and 6. The twin variant 1 was nucleated at the normal stress of ~20 MPa, which was close to that of the parent, as shown in Figure 29 (a-d). The twin variant 6 appeared at the strain value of 0.4%, in which the shear stress is ~19 MPa in the untwinned parent grain and ~ 10 MPa for the twinned grain. The results showed that the shear stress, which was resolved on the twinning plane, i.e., , was smaller for both twin variants compared to the parent grains, which was in line with the results shown in Figure 29 (a-d).

Detwinning
The deformation slip and twinning mechanisms were enough to define the response of the Mg alloys during simple monotonic loadings. However, in the case of complex loading paths, such as cyclic loadings, this was not the case. One should also define the detwinning mechanisms in these cases, as described in Section 2.3. In-situ Synchrotron Xray techniques have contributed to investigate the detwinning mechanisms. In the case of Mg and its alloys, Murphy-Leonard et al. [7] incorporated the HEDM technique to investigate the twinning-detwinning in pure Mg. They applied cyclic loading with three different applied strains of 0.4%, 0.52%, and 0.75%. The sample was extruded, and loading was applied along the extrusion direction. The c-axis was mostly perpendicular to the extrusion direction. Accordingly, extension twinning was a favored mechanism during compression loading. To monitor the twin volume, Murphy-Leonard et al. [7] measured the X-Ray diffraction peaks of the basal 0 0 0 2 planes. In other words, initially, there was no basal 0 0 0 2 peak intensity along the loading direction. As a twinned child is nucleated and grows, a 86.3° reorientation of the basal pole occurs within a parent grain, and the peak intensity starts increasing. Accordingly, Murphy-Leonard et al. [7] introduced the formulation for the twin volume fraction based on the HEDM data as follows: where is the basal 0 0 0 2 peak intensity along the loading direction, and is the initial basal 0 0 0 2 peak intensity along the normal direction. Figure 30 shows the cyclic response of the sample at different loading cycles, including the stress-strain loops, basal 0 0 0 2 peak intensity along the loading direction, and basal 0 0 0 2 peak intensity along the normal direction. In cycle 1, Figure 30 (b) shows that the basal 0 0 0 2 peak intensity along the loading direction remained zero during the tensile part of the loading, which means there no twin nucleation during the initial tensile loading. The basal 0 0 0 2 peak intensity along the loading direction started increasing at the C-1 point in compression and kept increasing until point B, which was the maximum compression loading. At the same time, at point C-1, Figure 30 (c) shows that the basal 0 0 0 2 peak intensity along the normal direction started dropping. This was in fact, the reorientation of the c-axis along the normal direction towards the loading direction due to twinning.
As soon as unloading starts, the basal 0 0 0 2 peak intensity along the loading direction starts decreasing, while the basal 0 0 0 2 peak intensity along the normal direction starts increasing. In other words, the previously twinned children during compression loading were reoriented back to the parent grains. The twin exhaustion occurred at the point T-2 during the second loading cycle, as shown in Figure 30 (d-f). The variation of basal 0 0 0 2 peak intensity in the second cycle followed the same trend similar to the first cycle. However, the twinning initiates at smaller compressive strains C-2 compared to C-1 in the first cycle. In both cycles 1 and 2, the twin was fully exhausted in the tensile loading. However, as shown in Figure 30 (g-i), all the twinned children were not fully detwinned as the minimum twin volume fraction can be observed at the maximum tension strain A. This twin volume content was named residual twins. Murphy-Leonard et al. [7] showed that the residual twins increased by the number of cycles. One can compare the residual twins in loading cycle 500 at maximum tensile strain (point A) to that of loading cycle 200 to observe the increase in the residual twin content as the number of loading cycle increases. . The results for a single untwinned parent grain and two twin variants nucleated from that grain (e) , (f) , (g) , (h) . (i) The schematics of the stress components presented in (a-h) (After Abdolvand et al. [27]).
Murphy-Leonard et al. [7] also investigated the variation of twin initiation stress and detwinning exhaustion stress versus the number of cycles, as shown in Figure 31, in the case of cyclic loadings with the strain amplitude of 0.75% and 0.52%. The results showed that as the number of cycles increased, the magnitude of twinning initiation stress decreased while the magnitude of the twin exhaustion stress increased. They further investigated the results and analyzed the variation of twin intensity at maximum tensile strain, which can be inferred as the twin residuals, versus the number of cycles in the case of cyclic loadings with the strain amplitude of 0.75 % and 0.52 %, as shown in Figure 32. The results showed that before the loading cycle 100, there was no twinning intensity at maximum tensile strain independent of the strain amplitude, i.e., there were no residual twins. In the cases of loading cycles higher than 100, the residual twin content increased as the number of cycles increased. Furthermore, the residual twin content of samples subjected to the strain amplitude of 0.75% was larger than those subjected to the strain amplitude of 0.52%. In other words, the residual twin content increased by the strain amplitude. An open-source CPFE software PRISMS-Plasticity [103] was later used to capture the cyclic response of the sample, as shown in Figure 33. The results showed that the CPFE simulation can successfully capture the experimental cyclic stress-strain response.

Micromechanics of Twinning
Brown et al. [34] incorporated the in-situ neutron diffraction to investigate the evolution of internal strain and texture during extension twinning in the AZ31B Mg alloy. They applied different loading paths that include in-plane compression (IPC), in-plane tension (IPT), and through-thickness compression (TTC), where the c-axis of most of the grains were perpendicular to the loading direction. Among these three loading paths, the extension twinning was the governing deformation mechanism in the IPC, which was the main focus of their study [34]. They used the VPSC framework to simulate the sample response during in-plane compression. Figure 34 presents the stress-strain responses of AZ31B Mg alloy during TTC and IPC. The circles on the IPC curve represent the strains at which the in-situ neutron diffraction measurements were conducted. The VPSC results were also presented as the dashed line for the IPC loading. The stress-strain response of IPC loading yields at the stress value of ~60 MPa, which was followed by a plateau with a small hardening rate. This plateau is commonly attributed to the role of extension twinning [34]. The hardening rate starts increasing at the inflection point that occurs at the strain of ~5%. For strain values larger than ~8%, the stress during IPC becomes equal or larger than that of the TTC loading, and the rate of work hardening decreases again. Brown et al. [34] investigated the diffraction pattern of the AZ31B Mg alloy subjected to the IPC loading path in parallel and perpendicular detector banks at different stress values, as shown in Figure 35. Initially, the parallel detector banker does not detect any diffraction peaks from 0 0 0 2 , while the transverse detector bank shows very strong diffraction peaks from 0 0 0 2 . As the applied compression increases, the diffraction peaks intensity detected by the parallel detector banker increases, while the one detected by the transverse detector bank decreases. Brown et al. [34] stated that the increase in the diffraction peaks intensity of 0 0 0 2 detected by the parallel detector banker manifests the development of twinned grains. Accordingly, they measured the twinning variation and showed that the twinning initiates at the strain ~1%. It was observed that the twinning volume linearly increases with the applied strain up to a strain ~6%. The twinning is then saturated at the strain ~8%, with the maximum twinning volume fraction of 80%. Brown et al. [34] mentioned that the variation of twin volume fraction is almost proportional to the increase in the diffraction peaks intensity of 0 0 0 2 detected by the parallel detector banker, with the slope of 0.145 per unit plastic strain. The VPSC predicted faster twin formation, which saturates at the strain ~4%. The increased hardening rate at larger strain observed in Figure 34 can be attributed to the twin saturation and initiation of pyramidal slip in the twinned regions [34]. Figure 36 shows the variation of lattice strain parallel and perpendicular to the loading axis in untwinned parent grains and daughter twinned grains versus the applied load [34]. Brown et al. [34] obtained the lattice strain for a grain orientation contributing to a given ℎ diffraction peak as follows: where denotes the unstrained interatomic spacing, which can be approximated by the initial interplanar spacing. Figure 36 (a) shows the internal strain in the untwinned parent grain parallel and normal to the loading axis versus the applied load. The theoretical linear elastic relation is shown as dashed lines. Prior to the applied compressive stress of 70 MPa, the untwinned parent grains behave according to the linear elasticity. Afterward, the lattice strain becomes less than the elastic prediction as other grains with different orientations, including the daughter grains, which sustain stresses larger than average applied stress.   [34] reported an interesting observation in which the lattice strain of newly twinned daughters was 700 , which was considerably less than the untwinned parent and the surrounding grains. In other words, the daughter twinned grains were in a relaxed state compared to the surrounding grains. In the region of applied compressive loading of 70 MPa 220 MPa, the increase in lattice strain was much larger than the elastic prediction. The stress/lattice strain slope was reported as 32 GPa, which was smaller than the elastic slope of 48 GPa calculated for grains with basal pole parallel to the loading axis. This can be attributed to the fact that the daughter twinned grains were hard orientations and sustain more elastic strain compared to the untwinned parent and surrounding grains, which have soft orientations and deform inelastically [34]. For applied stress greater than 220 MPa, one can see an inflection point in the variation of lattice strain versus the applied loading as the variation of strain sustained by the daughter grain is less than the elastic prediction.
Brown et al. [34] calculated the CRSS required for twin formation using the in-situ neutron diffraction experimental data. They observed that the twin formation occurs at the stress in the region of 50 MPa 70 MPa. In the case of the tested sample, the maximum Schmid factor is 0.499, which leads to the twinning CRSS of 25 to 35 MPa for AZ31B Mg alloy. In the last step, Brown et al. [34] tried to derive how much strain is sustained by twinning. They derived the contribution of twin to the total plastic strain as follows: where the coefficient 0.32 is obtained based on the texture of tested AZ31B Mg alloy, is the characteristic twin shear strain defined in Equation (10), and is the twin volume fraction. They estimated that 61% of the plastic deformation is sustained by twinning mechanisms at the strain of 5%. Afterward, the twinning rate versus the strain decreases in which it saturates at the strain of 8%. At this strain, 42% of the plastic deformation is sustained by the twinning mechanism.

Conclusions and Future Works
In the present study, the crystal plasticity models, which have been incorporated for Mg and its alloys, are reviewed. These models include different deformation mechanisms such as plastic slip, twinning, and detwinning. The recent experimental frameworks, such as in-situ neutron diffraction, 3D high energy synchrotron X-ray techniques, and digital image correlation under scanning electron microscopy, which have been incorporated along crystal plasticity models to investigate the properties of Mg and its alloys, are also reviewed. Although many studies have tried to model the behavior of Mg and its alloys using crystal plasticity, which was often coupled by the advanced in-situ techniques, there are still some critical challenges that need to be addressed. A summary is presented for the future challenges as below:  Real-time crystal plasticity simulation coupled to in-situ experiments to guide identification of outliers that can in-turn improve crystal plasticity theories.  A general map to include the effect of alloying for a variety of Mg alloys using crystal plasticity models along with synchrotron X-ray techniques in a consistent framework.  Using machine learning techniques to learn the crystal plasticity models and generate surrogate models which can be used to design specific Mg alloy loading paths to achieve target properties.  Developing a crystal plasticity model with a physically based twinning and detwinning model, which include the correct isotropic and kinematic hardenings to capture the appropriate cyclic response of Mg alloy. This is extremely important in the prediction of fatigue simulation using crystal plasticity simulation.  Developing an integrated framework of crystal plasticity models and phase field simulation to better capture the twin morphology in Mg alloys.  The interaction of slip modes and twinning and detwinning mechanisms, which is typically reflected as latent hardening in crystal plasticity models.
 Improved modeling of slip/twin and grain boundary interactions: Effects of grain size in different Mg alloys using crystal plasticity models, specifically via micro-Hall Petch models whose parameters can be inferred through experiments [98,102] and including the effect of grain boundary on twin nucleation and growth in crystal plasticity models.  Integrating the crystal plasticity models of Mg and its alloys with the PRISMS-Fatigue framework [109] to investigate the effects of texture, grain morphology, sample size, multiaxial strain, and strain amplitude on their fatigue response.  Coupling the crystal plasticity models with phase-field simulations to address the effect of deformation mechanisms such as plastic slip and twinning on the dynamic recrystallization of Mg alloys.