Interest in hexagonal close-packed (hcp) magnesium alloys in structural applications has intensified over the last few decades due to their low density, high specific strength, abundance, and recyclability [1,2] . A common strategy to bolster the mechanical properties of Mg alloys, such as AZ91, is to form a dispersion of precipitates throughout the microstructure that blocks the glide of dislocations [3][4][5][6][7] . The precipitates that form in AZ91 alloy, for instance, are either thin lath-shaped basal or rod-shaped prismatic Mg 17 Al 12 precipitates with a body-centered cubic (BCC) crystal structure [8,9] . The habit planes of basal-and prismaticprecipitates are parallel to the basal and prismatic planes of the parent hcp phase, respectively. Plastic deformation of the AZ91 al-loy under relatively low strains and at room temperature is accommodated by both dislocation slip and deformation twinning [9,10] . The amount of twinning that develops affects the ductility, mechanical stability, strain-hardening response, and texture evolution [11][12][13][14][15][16] . Compared to precipitate-dislocation interactions, the interactions between the precipitates and twins have received less attention but are clearly of equal importance in understanding the deformation behavior of precipitate-strengthened Mg alloys.
Twin development generally consists of three stages: nucleation, propagation, and growth, and precipitates can have a distinct influence on each [17] . Precipitate/matrix interfaces have been suggested to serve as nucleation sites due to their high defect content and stress concentrations [18] . Consistent with this notion, some studies have reported greater numbers of twin lamellae but similar twin volume fractions with precipitates than without [19][20][21] . In propagation, precipitates act as obstacles, causing the moving twin-tip to fully arrest or skew from the twin habit plane to bypass the precipitate [17,22] . Where the twin-tip impinges the precipitate, and the ratio of the twin thickness to precipitate width have been found to influence the path taken by the twin-tip [23] . Specifically, in our previous work, we have shown that the interac-tion between an impinging precipitate and twin develops a highstress concentration on the opposite side of the precipitate. The intensity of stress localization depends on the relative sizes of twins and precipitates, and also the interaction location, such as the precipitate tip versus the center of the precipitate. This local stress favors the formation of a new twin on the opposite side of the precipitate.
In regards to twin growth, it is commonly reported that the presence of precipitates hinders the twin thickening [19,[24][25][26][27][28][29] . Several numerical studies aim to quantify the effect of precipitates on twin growth by simulating the following two different configurations: (i) twin boundary interaction with the precipitates, and (ii) twin tip interaction with the precipitates. For example, using the phase field model, Liu et al. simulated the migration of tensile twin boundaries in the presence of precipitates [30] . In their simulations, a perfect twin boundary that is not pinned by any obstacles starts to migrate uniformly under macroscopic loading and gets pinned locally by the precipitates. In order to continue the migration of twin boundaries, additional macroscopic loading is needed to de-pin it from the precipitates. Using these simulations, the authors have reported an effective critical resolved shear stress (CRSS) for twin growth in the precipitate-hardened material, which is significantly higher than in precipitate-free materials and increases linearly with the precipitate volume fraction. On the other hand, using crystal plasticity modeling, Siska et al. studied the interaction of a twin tip with the precipitates and its influence on twin thickening [22] . The authors found that the interaction of a twin with the hard precipitates alters the CRSS for slip activity, and so increases the CRSS for twin growth.
Evidently, the interactions between propagating twins and precipitates are complex, making it difficult to predict how twins may develop in the presence of arrays of precipitates. Clarifying how twins interact with precipitates and eventually bypass them or engulf them is desired in the design of precipitate-strengthened microstructures in Mg alloys, such as AZ91. To the best of our knowledge, none of the existing literature rationalize the experimentally observed large twin domains by investigating the growth of twins arrested by arrays of precipitates in AZ91 alloy.
In this work, the interactions between Mg 17 Al 12 plate-shaped basal precipitates and { 10 1 2 } tensile twins at different stages of twin development in AZ91 are studied. Specifically, the growth of a twin that is blocked by precipitates is studied. An elasto-viscoplastic fast Fourier transform (EVP-FFT) model with a dislocation-density (DD) based hardening law is employed to calculate the spatial distribution of internal micromechanical fields and dislocation densities corresponding to each interaction. The results indicate that twins interacting with precipitates grow asymmetrically, with one twin boundary remaining relatively parallel to its crystallographic twinning plane and the other forming a serrated structure that deviates from its crystallographic plane. We discuss how these findings can explain the twin morphologies near precipitates observed in transmission electron microscopy (TEM) observations of deformed AZ91 [31,32] . Last, we propose a twostep mechanism for the growth of an early-stage thin twin that is enabled by interactions between sequentially formed twins and a precipitate.
We employ a crystal-plasticity-based elasto-viscoplastic fast-Fourier-transform (EVP-FFT) model that was originally developed by Lebensohn et al. [33] and later advanced by Kumar et al. [34] for explicit twins. A similar modeling framework was recently used to study pathways taken by twins after the initial interaction with a precipitate [23,35] . Here, hardening via dislocation density evolution in the HCP slip families is included. Both the EVP-FFT formulation and dislocation density law are briefly reviewed below.
In EVP-FFT, the initial microstructure is represented by a periodic unit cell, discretized into regularly spaced voxels whose center point is located at x . The Cauchy stress tensor field, σ (x ) , at a time t + t is described by a time discretized generalized Hooke's law:
where C(x ) is the elastic stiffness tensor, and ε e (x ) is the elastic strain tensor. Under a small strain approximation, ε e (x ) is the difference between the total strain tensor, ε(x ) and the inelastic strain tensor. The inelastic strain is given by the sum of the plastic strain due slip, ε p (x ) , and the homogeneous twinning shear transformation strain in the twin domain, ε tw (x ) , i.e.,
The ˙ ε p (x ) is constitutively related to the stress at each material point through a sum over the N active slip systems:
In the above, the Schmid tensor, m s , is the symmetric part of ( b s n s ) , where b s and n s are the unit vectors along the slip direction and the slip plane normal, respectively, of slip system s .
The slip resistance τ s c (x ) is the current critical resolved shear stress needed to activate slip system s . The rate ˙ γ o is a normalization factor and n is the stress exponent. In the twin domain, the total plastic strain due to twinning is equal to the product of the characteristic twin shear and its volume [36] . For convergence, it is necessary to apply the twinning strain increment, ε tw , on the twin plane in the twinning shear direction with sufficiently fine increments (i.e., 20 0 0 steps in this case). Tensors ε tw and ε tw outside of the twin domain are zero.
The evolution of τ s c of every slip system, s , is calculated based on a dislocation density (DD) based hardening law [37] . Each slip system has its own relationship to strain rate and temperature. The τ s c is comprised of a sum of multiple factors given by:
where τ s 0 represents the initial slip resistance, or Peierls stress, of slip system s . This resistance τ s 0 is dependent on the temperature, initial dislocation density, and intrinsic material properties. In this work, it is assumed to be a scalar quantity and scales with temperature as:
Here, T is the temperature (in K), and B and C are the material parameters. The resistance τ s debris is associated with the resistance to slip provided by arrays sessile dislocation substructures and is given by:
where K is a material-independent constant and is assumed to be 0.086 [37][38][39] . μ is the temperature-dependent shear modulus, b s is the Burgers vector, and ρ debris is the total debris dislocation density. The generation of ρ debris is equal to the sum of debris contributed by each slip system and is written as:
where q is set equal to 4 and defines the rate coefficient for how debris grow from point defects. The density ρ s rem is the population of dislocations from each slip system that contribute to the formation of debris and f s is a term that describes the dislocationdislocation interactions that lead to debris formation and can be written as:
In the above, A s is a term that relates the temperature dependency of dislocation interactions. Finally, τ s f orest represents the resistance to slip due to a random arrangement of forest dislocations [40,41] and is related to the square root of the dislocation density, ρ s , of slip system s via Taylor's law, i.e.,
In the above, χ is the dislocation interaction coefficient and is set to 0.9. h ss is a latent hardening matrix for the interactions between dislocations from different slip modes. The stored forest dislocation densities evolve via a competition between a dislocation generation term, k s 1 , that describes the storage of dislocations by statistical trapping of gliding dislocations by forest dislocations and a thermally controlled annihilation term, k s 2 , for dynamic recovery of stored dislocations [42][43][44] :
The annihilation term, k s 2 is given by [37] :
where k , ˙ ε , ˙ ε 0 , D s , and g s is the Boltzmann's constant, strain rate at each voxel, reference strain rate (set to be 10 7 1/s), drag stress, and normalized activation energy of slip system s .
Before performing the explicit twin-precipitate interaction simulations within an individual crystal, we need to determine the DD-based hardening model parameters of the dominant slip modes in precipitate strengthened AZ91. Experimental uniaxial tensile stress-strain curves at 293 K and 463 K for a strain rate of 1 × 10 -4 1/s of this alloy from Dini et al. [45] are used to aid parameter estimation. The model finds a single set of material DD parameters that will enable the DD-EVP-FFT simulations of a polycrystal of the same alloy to achieve agreement with all measured curves. Fig. 1 a shows the simulation cell of the polycrystalline microstructure. It is discretized into 6 4 × 6 4 × 6 4 voxels and consists of 128 equiaxed grains with a uniform distribution of crystallographic orientations and approximately 7.5% volume fraction of Mg 17 Al 12 precipitates, colored in black [45] . Accordingly, every grain is represented by approximately 2048 voxels. The Mg 17 Al 12 precipitates have a body-centered cubic (bcc) crystal structure and share a Potter orientation relationship with the parent grain that it resides in, described by ( 0 0 01 ) HCP 2 °from (011) BCC , [ 2 11 0 ] HCP // [ 1 1 1 ] BCC , and ( 01 1 1 ) HCP //(110) BCC [8,9,17] .
Over the deformation period, the Mg 17 Al 12 precipitates remain elastic, while the AZ91 matrix deforms elastically and plastically.
The cubic elastic constants C 11 , C 12 , and C 44 of the Mg 17 Al 12 precipitates at different tem peratures are not available, and thus, values calculated from Density Functional Theory (DFT) were used corresponding to 86.8, 29.0, and 20.0 GPa, respectively, for both temperatures [46] . The five elastic coefficients C 11 , C 12, C 13, C 33, C 44 of AZ91 were approximated by the coefficients of pure Mg at 293 K and 463 K with values of 59.0, 25.5, 21.5, 61.0, 16.7 and 57.0, 25.4, 21.1, 58.0, 15.0 GPa, respectively [47] . The shear modulus, μ, used in Eqs. ( 6) and ( 9) , for AZ91 at 293 K and 463 K were set to equal to C 44 , which is 16.7 and 15.0 GPa, respectively [48] . For plastic slip in the AZ91 matrix, three slip modes were made available: prismatic a , basal a and pyramidal-I c + a slip. Twinning deformation was not taken into account in the parameterization since its contribution was reported to be negligible in deformation [49] . The simulation begins with the initial values of the slip resistance for prismatic, basal, and pyramidal-I c + a slip modes for AZ91 from Bhattacharyya et al. [10] . with self and latent hardening set to be 0.15 [50] .
Fig. 1 b compares the calculated stress-strain response (red solid lines) to the experimentally measured curves (symbols). The corresponding constitutive law parameters for all three slip modes are listed in Table 1 . These parameters are characterized for the matrix material in AZ91 outside of the precipitates and therefore are independent of precipitate size or distribution. The model operates on the mesoscale length scale and does not account for the individual interactions of single dislocations with precipitates. It is possible for an additional term, τ s precipitate , to be included in Eq. ( 4) to phenomenologically capture the effects of the precipitate, since, in reality, the precipitate spacing can affect the flow stress of dislocations and twins. Here, the effect of the precipitates are accounted for by explicitly distributing the precipitates throughout the microstructure to closely replicate the microstructures commonly reported in experiments [31,32,45] . Thus, the simulated results may be to closely tied to the size, shape, distribution, volume fraction and type of precipitate. For validation, we also compare the evolution of the total (gray), prismatic (red), basal (blue) and pyramidal-I (green) dislocation densities with respect to applied strain at 293 K and 463 K in Fig. 1 c andd, respectively, with the measured densities (designated by a square marker) by Dini et al. [49] . Notably, the predicted densities fall within the reasonable agreement of experimentally measured values at both temperatures.
In the calculations that remain, we consider twin/precipitate interactions within crystals. The reorientation and shear of a finite twin lamella are represented explicitly through the reorientation of a predefined twin domain and the incremental development of the characteristic twinning shear homogeneously within the twin domain [34] . For the { 10 1 2 } tensile twin, the voxels within the domain are reoriented 86.3 °about 1 2 1 0 and uniformly sheared by the characteristic twinning shear of 12.9%. The precipitates are also predefined domains with crystallographic and material properties of the commonly observed Mg 17 Al 12 plate-shaped precipitates that lie on the basal plane of the parent matrix. Fig. 2 depicts the simulation setup with two precipitates of the same size: precipitate-1 and precipitate-2 in orange and red, respectively. The twins and precipitates are embedded in an AZ91 matrix (light blue). The crystallographic orientation in the Bunge convention of the parent grain is chosen to be (0 °, 88.16 °, 0 °), and it is illustrated by the hexagonal inset in Fig. 2 . Following the frequently observed Potter orientation relationship [8,9,17] , the crystallographic orientation of the precipitates is assigned to be (39.23 °, 114 °, 63.43 °). Other orientation relationships and variations in orientation may exist between the matrix and precipitates and can be easily accounted for in the model. In the present work, the most
Dislocation density-based constitutive model parameters that are found by fitting to the experimentally obtained stress-strain response of AZ91 alloy under uniaxial tension at 293 K and 463 K.
3.231e-10 3.231e-10 6.0771e-10 ˙
5.0e + 03 5.0e + 03 9.0e + 03 commonly observed orientation relationship was considered for simplicity.
The simulation cell size is discretized into 3 × 550 × 550 voxels. During the early stages of propagation, twins can assume a lamellar shape and may eventually thicken to an elliptical shape under further loading after the twin tips have become arrested [51][52][53] . The present work aims to study precipitate-twin interactions during the propagation stages of twinning, before significant twin thickening has occurred. The twin thickness, t, and precipitate width in the Y-direction are both four voxels, appropriate for representing thin, newly propagating twins. In experiments, the aspect ratio of precipitates can be seen to vary between 1 -100 [32] . Thus, the precipitate length in the Z-direction is taken as 40 voxels that corresponds to the aspect ratio of ten. The twin-tip/precipitate interactions depicted in Fig. 2 is often seen in experimental observations [32] . Furthermore, the precipitate length was chosen in or-der to keep the twin impingement site far away from the edges since these events make up a majority of twin-precipitate interactions. Edge impingement can also occur in either the Z-direction or the X-direction; however, it is less likely to occur due to the relatively large surface area of the plate-shaped precipitates. In the event that edge impingement does occur, the resulting stress states are more complex [23] . In this case, the propagating twin may not be effectively arrested by the precipitates and, as a result, the twin growth mechanism is distinct from the what is reported in this work. Seven different precipitate spacings, d, were considered, ranging from 4 to 200 voxels, representative of what is commonly observed in deformed AZ91 [32] .
The parent matrix is surrounded by a buffer layer (dark blue) of 50 voxels thick in the Y and Z directions, with uniformly distributed crystallographic orientations. This buffer layer approximates the response of a polycrystalline medium and shields the parent matrix from overlapping fields due to the periodic boundary conditions. There is no buffer layer in the X-direction. The simulation uses periodic boundary conditions, which causes the microstructure to adopt a columnar microstructure. This microstructural simplification is appropriate to use when studying the interactions of twins with thicknesses that are much smaller than the dimensions of the precipitates in the X and Z-direction, such as twins during the early stages of propagation being studied here. Furthermore, the crystallography of the twin, with its twin plane normal and twin shear direction in-plane, permits study of a columnar structure, as shown in Fig. 2 since any gradients in the out-of-plane direction will be negligible.
We adopt a simplification of the intricate microstructures found within twin and precipitate-hardened AZ91 alloy. The real microstructure may contain many twins and precipitate with geometries that are more complex than those represented here. Consideration of other types of precipitates and significantly different distributions of precipitates would involve recalibrating the model. However, the idealized model scenario permits the isolation of the interactions between two sequentially formed twins without the additional complications from other precipitates and twins. The model setup describes situations that may occur just as a twin has finished propagating in the twinning shear direction and has been arrested at a precipitate. Previously, such calculations have shown that depending on the precipitate-1 thickness with respect to the twin-1 thickness, a second twin, twin-2, may form on the other side of precipitate-1 and may eventually become arrested by another precipitate, say, precipitate-2 [23] . When faced with an array of precipitates, twins may continue to successively impinge on one side of each precipitate and reform on the other side in order to continue propagating. However, it is often observed that twins can also thicken and eventually engulf precipitates as they propagate through an array of precipitates. The influence of precipitates and precipitate spacing on twin thickening after impingements have yet to be investigated. The calculations that follow aim to identify the likely growth mechanisms of twin-1 and twin-2 and the role of the inter-precipitate spacing that stands in their way.
Uniaxial compression strain is applied in fine strain increments to a twin-free, precipitate-containing crystal in the a -axisdirection of the parent grain. To determine the path taken by a twin approaching the precipitates, the twin plane resolved shear stress (TRSS), a relevant driving stress component for twin formation or twin boundary propagation, is examined. An initial 4% compressive strain is applied, comparable to the strains experienced in precipitate-hardened AZ91 in the work by Xie et al. [32] , where many newly propagated twins can be seen interacting with precipitates. A TRSS of approximately 30 MPa is produced in the parent matrix, appropriate for the formation of tensile twins [54][55][56] . The stresses are inhomogeneous and result from the mismatch in elastic and plastic properties between the elastically hard precipitates and elasto-visco-plastically deforming AZ91 matrix [34,57] . Plastic deformation is seen to concentrate in the matrix around the hard bcc precipitates, regardless of precipitate spacing. Fig. 3 shows the distribution of the average slip resistance values after the initial 4% compression for each slip mode in the parent matrix in a small region around the precipitates outlined by dashed lines in Fig. 2 . For all three slip modes, the slip resistance increases or hardens the most along the short sides of the precipitates and the least along the long sides of the precipitates. Basal slip experiences the most significant hardening. The lowest slip resistances values were 2.9 times greater than its initial CRSS (44.8 MPa). Larger resistance to plastic deformation around the precipitates can result in higher stresses and stress concentrations in the parent grain with the addition of approaching deformation twins. Fig. 4 a shows the TRSS field of the twin-free, precipitatecontaining crystal at 4% compression strain. Along the short sides of the precipitates, the TRSS is approximately -81 MPa and resists twinning. The TRSS values along the long sides of the precipitate/matrix interface, however, are positive and enhanced. The peak TRSS lies at the center, reaching up to 32 MPa, more than three times higher than the surrounding matrix ( ∼ 9 MPa). The positive twin stresses can aid twin development, whether nucleation of a new twin from the interface or the growth of twin that propagates into the vicinity of the precipitate. The former situation would support the experimentally observed increase in the twinning density with the precipitates [21,26] , and although important, we focus this study on the latter situation. If an incoming twin impinges the short ends of the precipitates, its stress field will interact with an anti-twinning contribution from the precipitate, possibly resisting the propagation of the twin. In contrast, it would encounter less resistance should it impinge the center of the long side of the precipitate.
To study the stress fields corresponding to precipitate/twin interactions, deformation twins are introduced to the crystal containing precipitates under a fixed macroscopic loading. We first con-sider the case that twin-1 propagates up to precipitate at its center, while under 4% compressive strain. Twin-1 is much finer in thickness than the precipitate length, l, in the Z-direction. Fig. 4 b shows the calculated TRSS field resulting from this interaction. Fig. 4 b shows the calculated TRSS field resulting from this interaction. We note that the TRSS fields are qualitatively similar to the fields calculated and reported in [23] without pre-straining. A TRSS concentration forms on the opposite side of the precipitate, in the region ω opp , reaching up to 155 MPa, about five times greater than that from the initial macroscopic loading alone ( ∼ 32 MPa). A new twin of the same variant could, therefore, form in ω opp . To determine how this new twin could expand while twin-1 remains arrested, Fig. 5 c andd present the TRSS field in between the two adjacent precipitates along two lines in the vicinity of where this new twin would be. Near ω opp , corresponding to the left side of Fig. 5 c, andd, the TRSS is greatly enhanced from ∼32 MPa before twin-1 impingement (gray points) to values between 135 and 150 MPa after twin-1 impingement (red points). The five-fold increase in driving force supports forming a second twin-2 of the same variant as the first.
Fig. 5 a and b closely examine the calculated TRSS profiles in the parent matrix along the upper and lower twin boundaries of twin-1, which are relevant for twin boundary migration during the twin growth stage. Along the entire upper and lower boundary of twin- 1, the TRSS is negative, indicating resistance against twin boundary migration and the need for more applied deformation for twin growth. Here, we classify the decrease in TRSS stresses resulting from the reaction of hard, elastic precipitate to the twin shear as backstresses [57] . The largest backstresses are experienced near the precipitate, corresponding to the right side of Fig. 5 a andb. Since the backstresses are much less at the upper twin-1 boundary than the lower twin-1 boundary at the twin/precipitate junction site, it is anticipated that twin-1 would expand asymmetrically, initiating from the upper boundary first.
Based on the foregoing observations, calculations are performed for the situation in Fig. 2 , in which twin-2 is formed where the peak TRSS had developed and extended between precipitate-1 and precipitate-2 and made with equal thickness as twin-1. In the simulation, this process occurred under the same macroscopic loading after the formation of twin-1. Reminiscent of twin-1, the leading tip of twin-2 causes a rise in TRSS of up to ∼ 155 MPa on the opposite side of precipitate-2, seen in Fig. 4 c. Severe backstresses result in intensely negative TRSS in the vicinity of twin-2 boundaries. Specifically, along the twin-2 boundaries, the magnitude of the TRSS is nearly -160 MPa where twin-2 impinges upon both precipitate-1 and precipitate-2 at an acute angle, that is, on the left of the upper boundary and right of the lower boundary, see blue profiles in Fig. 5 c andd. The lower twin-2 boundary near precipitate-1 develops slightly higher TRSS due to its closer alignment with the peak TRSS concentration generated by twin-1 [23] . It is reduced to -100 MPa near the center of the twin. At the same time, the backstresses along the boundaries of twin-1 are relieved by the formation of twin-2, shown by the blue markers in The augmented TRSS fields from twin-2 on twin-1 suggest a two-step mechanism for the growth of twins that are interacting with precipitates. First, twin-1 must propagate and impinge upon a precipitate, triggering the nucleation and propagation of twin-2 on the other side of the precipitate. Secondly, twin-2 encourages the growth of twin-1 along the upper boundary near the precipitate that they co-impinge upon. With the formation of twin-2, little to no additional applied stresses would be needed for the growth of twin-1. The mutual aid between twin-1 and twin-2 allows for the migration of twin-1 boundary around precipitate-1. This result is similar to the observation reported in Kumar et al. [58] , where the effect of twin-twin interactions across the shared grain boundary fostered their mutual growth.
Thus far, the precipitate spacing d has been fixed. The effects of d on twin-1 growth are investigated next. Fig. 6 shows the evolution of the TRSS along the twin-1 boundaries after each step in the interaction process: 1) the initial 4% a -axis compression (gray), 2) the formation of twin-1 (red) only, and 3) the formation of twin-1 and then twin-2 (blue). In all cases, twin-2 extends from the first precipitate to the second and has length d . The average TRSS in the two regions shown in Fig. 6 a are studied: at the twintip/precipitate junction (solid lines) and in the center of the twin lamellae (dashed lines). Each region is equal to the twin thickness. In the second step, when the twin first interacts with the precipitate, the TRSS fields around twin-1 under external load are insensi-Fig. 6. a) An enlarged view of the simulated twin-precipitate interaction region. The fields along the upper and lower boundaries of twin-1 near the center of twin-1 and near the precipitate-1 interaction site are outlined by black boxes to be considered for further detailed analysis. The variation in average TRSS taken along the (b) upper and (c) lower twin boundary of twin-1 near the center and the impingement site as a function of precipitate spacings, d, normalized by the twin thickness, t. The average TRSS values after an initial 4% compression (prior to the formation of twin-1), after the formation of twin-1 (prior to the formation of twin-2), and after the formation of twin-2 are shown in gray, red, and blue lines, respectively. Solid and dashed lines represent the averages taken near the impingement site and center of the twin-1, respectively. tive to the spacing d . It appears that the heterogeneous stress fields around precipitate-2 are too far away to influence the TRSS fields at twin-1, even for the closest spacing of d/t = 1. Fig. 6 b shows that for all d/t ranging from 1 to 50, the average TRSS along the upper boundary reduced from ∼ + 20 MPa (gray dashed lines) to ∼ -20 MPa in the center (red dashed lines) and from ∼ + 30 MPa (gray solid lines) to ∼ -50 MPa at the twin tip (red solid lines) after twin-1 impinges on precipitate-1. In Fig. 6 c, for the lower boundary of twin-1, we observe a reduction in TRSS to ∼ -20 in its center and to ∼ -130 MPa at its tip. Thus, we find that migration of the lower boundary near the twin tip is more limited than that of the upper boundary.
The spacing d has an important influence on the stress field once twin-2 has formed and extended to precipitate-2. For tight precipitate spacings, d = 1 t to 5 t , the average TRSS along the top boundary of twin-1 increases rapidly from ∼ -20 to 62 MPa near the tip (solid blue lines) and modestly from ∼ -20 to -7 MPa near the center (dashed blue lines) of twin-1, seen in Fig. 6 b. The larger d , the greater the enhancement by twin-2 on twin-1, with the largest effect near the tip. For d > 5t , the enhancement saturates with TRSS values reaching up to ∼ + 80 and -2 MPa at the tip and central regions, respectively. Similar trends are observed along the lower boundary of twin-1, however, the TRSS values re-main negative and less favorable for twin boundary migration. As shown in Fig. 5 , for all d/t cases, both upper and lower boundaries of twin-2 remain resistant to twin boundary migration and consequently are not shown here. Although the subsequently formed twin-2 may not thicken, it can support the migration of the upper boundary of twin-1 when d is large, enabling twin-1 to grow around the precipitate. Without twin-2, the ability of a twin-1 to thicken independently and bypass the precipitate is unlikely due to the large backstresses that develop when the twin becomes constrained at both tips.
To sum, larger precipitate spacings further enable twin growth. However, a saturation limit exists for d ≥ 5 t , beyond which significant enhancement of twin growth is not expected. The likely region for the initiation of twin growth is shifted from the center of the twin to the twin tip near the upper twin-precipitate junction. In principle, this boundary migration mechanism for circumventing precipitates can continue with each new precipitate encounter. After twin-2 impinges precipitate-2, another twin, twin-3, can form on the opposite side of precipitate-2, resulting in the relief of the backstresses on twin-2 and helping it to grow. Several assumptions have been made in this study, including the amount of external load prior to twinning and relative thickness of the twin to the precipitate, which may alter the calculated twin stresses. However, it is found that even without any external load and for twin thickness twice are large that considered here, the results are qualitatively similar, yielding saturation limit at d ≥ 5 t These results are presented in the Appendix. Other assumptions made within this study involve the relative precipitate size, shape, distribution, volume fraction and types of precipitates considered in this study. Changes to these factors may potentially influence the calculated stress fields.
In the following section, we analyze and present the model predicted dislocation density evolution near the precipitate-twin interaction site. The strain fields calculated in the simulations are associated with slip activity. Fig. 7 depicts the dislocation density fields that arise during the different stages of the simulation: after the pre-twin strain, after forming twin-1, and after forming twin-2. In Fig. 7 , the precipitates are colored in black. After pre-twin straining, a heterogeneous distribution of dislocations emerges around the precipitates, with dislocation densities much greater than the initial value of 1 . 5 × 10 12 m -2 for all slip modes. Prismatic activity in the parent matrix was the greatest, resulting in densities that grew by two orders of magnitude. Greater densities of prismatic dislocations developed near the short edges of the precipitates than compared to the long edges of the precipitates, as seen in Fig. 7 a. The opposite trends are seen for the pyramidal dislocation densities, resulting in more dislocations along the long edges of the precipitates than the short edges, seen in Fig. 7 c. The development of basal dislocations was concentrated at the corners of the precipitates.
After forming twin-1, prismatic dislocation densities remained relatively constant, and additional pyramidal dislocations were generated only in the twin near the impinging twin tip. The basal dislocation field experiences the greatest change, with increased densities all along the precipitate/matrix boundary and large concentrations of dislocation appearing inside the impinging twin tip and in ω obt . After twin-2 forms, again, prismatic dislocation densities do not change. Evidently, twinning does not significantly influence prismatic slip activity. Pyramidal dislocations were primarily generated at the twin tips of twin-2 and twin-1 tip that coimpinges on precipitate-1. The formation of twin-2 caused basal dislocations densities to increase significantly inside and around each tip of twin-2. Furthermore, additional basal dislocations can be observed in the twin-1 tip that co-impinges the precipitate and, in the region ω obt .
Next, a closer inspection of the regions near the center of the twin-1 and near the twin-1 and precipitate-1 interaction site, as outlined in black boxes in Fig. 6 a, is presented. For the following analysis, only the upper boundary of twin-1 is considered since it has favorable stresses to support its migration, as seen in Fig. 6 . The increase in average dislocation density for each slip mode, i , due to the formation of twin-2 is calculated as
ρ a f ter twin -2 s -ρ be f ore twin -2 s (12) where ρ s is the dislocation density of slip system, s , in the region either near the twin center or impingement site. N i is the number of slip systems in each mode: prismatic, basal, or pyramidal-I.
Fig. 8 shows the change in the total average stored dislocation densities per slip family, prismatic , basal , and pyramidal -I , that accrues due to the formation of twin-2 for d between 1t to 50t The results show that twin-2 does not generate additional dislocation densities of any kind near the center of the twin-1, indicated by hollow box symbols in Fig. 8 . However, near the twin-1 tip, the total stored basal dislocation density increases by ∼ 6 × 10 12 for even the smallest precipitate spacing considered, d = 1 t , and increases up to ∼ 7 . 2 × 10 12 for d = 50 t , shown by solid markers in Fig. 8 b. The total stored pyramidal I dislocations also increased with the formation of twin-2, however, to a much lesser degree. The increase in total stored pyramidal-type I dislocation densities ranges from ∼ 0 . 9 × 10 10 to 2 . 7 × 10 10 for the corresponding d ranges from 1 t to 50 t , seen in Figure 8c. Pyramidal c + a -type I densities steadily increased with increasing d . The development of twin-2 is observed to increase the local concen-tration of defects in ω obt . Overall, the interaction between twin-2 and precipitate-1 increases the favorable driving stresses for twin-1 growth (shown in Fig. 6 ) and increases dislocation structures in the region ω obt .
It is generally thought that twin growth in precipitate hardened Mg-alloys is difficult because the elastically hard precipitates can pin migrating twin boundaries. This notion mostly does not account the effect of precipitates on twin front propagation. Instead, most of these studies directly start with a fully propagated twin and the migration of the twin boundaries through twin thickening in the twin plane normal direction [34] . Here, a mechanism is identified that allows twins to thicken around an array of precipitates, beginning with a newly propagating twin. Typically, the twin is thought to be arrested; however, when the approaching twin is of comparable thickness as the precipitate, the proposed two-step mechanism provides an alternative path for twins to engulf precip-itates, distinct from the studies that only treat twin growth as the translation of the entire twin boundary in the twin-plane-normal direction [28,30,50,59] .
The analysis suggests that the two-step bypass mechanism could explain the experimentally observed asymmetric twins and serrated morphology of twin boundaries in the presence of precipitates [32] . When twin-1 expands around the precipitate, its growth is easier along one side. For all the precipitate spacings considered here, the upper boundary gains sufficient driving forces to migrate while the lower twin boundary does not. Twin growth would initiate where twin-1 and the precipitate come together at an obtuse angle, leading to the formation of serrated twin boundary structures, as schematically illustrated in Fig. 9 a. Such serrated twin boundary structures are observed in the postmortem TEM micrographs of heavily twinned AZ91 with dispersions of Mg 17 Al 12 basal-plane precipitates, see Fig. 9 b [23] . Serrated features are exhibited by white arrows, and relatively straight twin boundaries are highlighted by white dashed lines, consistent with our simulations. Upon further loading, as the twin domains can continue to expand and engulf the precipitates, the serrated features from the locally augmented stress fields from the twin-precipitate interactions may be washed out by the driving stresses for twin expansion. The two-step twin growth mechanism described here helps to explain how large twin domains develop through the interactions of propagating twins and precipitates.
There is also the possibility that a second twin of the same variant could be independently formed, and it interacts with the same precipitate as the first twin but on two opposing sides of the precipitate. If this were the case, then the site where the second twin interacts the precipitate could occur anywhere along the long side of the precipitate. Two twins independently interacting the same precipitate may help to mutually relieve their backstresses if the location of their interaction sites were positioned near each other and on opposite sides of the precipitate. Any misalignment of the two independent twins could cause preferential growth of one of the twin boundaries over the other depending on how the two twins are offset from each other.
An elasto-viscoplastic fast-Fourier transform (EVP-FFT) model has been used to calculate the micromechanical fields and dislocation content associated with the interaction of a twin lamella and a precipitate during the early stages of twin propagation and growth in presence of an array of precipitates in AZ91. Specifically, the growth of a precipitate-blocked twin is investigated. The proposed twin-growth mechanism may help to explain how large twin domains can develop despite frequent reports that twin growth is suppressed in precipitate hardened magnesium alloys. Initially, it is reasoned that an initial thin twin can impinge upon a precipitate and cause another twin to form on the other side of the precipitate. The interactions between the precipitate with the impinging twin and the subsequent twin changes the local stress fields and dislocation content. The key findings are summarized below:
• The growth of a precipitate-blocked twin can be realized in two-stages. First, the initial propagation and impingement of a twin on a precipitate cause a local stress concentration to develop on the opposite side of the precipitate. Secondly, the formation and propagation of another twin from the far side of the precipitate counteract the backstresses that develop from the impingement of the initial twin. The augmented TRSS fields along the boundaries of the initial twin support its growth without any additional external macroscopic loading. • The amount of relaxation provided by the newly formed second twin on the first twin is related to the length of the sec-ond twin, which, in precipitate-hardened AZ91, is mostly determined by the precipitate spacing, d , a microstructural feature that can be controlled through alloying and heat treatment. A precipitate spacing that is five or more times the twin thickness allows for nearly the greatest enhancement of the driving stresses for twin growth. • The propagation of twin-2 also increases basal a and pyramidal-I c + a dislocation densities near the twin-1 impingement site.