{'fraction_non_alphanumeric': 0.03629982902906506, 'fraction_numerical': 0.013929631962566363, 'mean_word_length': 4.625797306874557, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Effect of a Micro-scale Dislocation Pileup on the Atomic-Scale Multi-variant Phase Transformation and Twinning

In this paper, we perform concurrent atomistic-continuum (CAC) simulations to (i) characterize the internal stress induced by the microscale dislocation pileup at an atomically structured interface; (ii) decompose this stress into two parts, one of which is from the dislocations behind the pileup tip according to the Eshelby model and the other is from the dislocations at the pileup tip according to a super-dislocation model; and (iii) assess how such internal stresses contribute to the atomic-scale phase transformations (PTs), reverse PTs, and twinning. The main novelty of this work is to unify the atomistic description of the interface and the coarse-grained (CG) description of the lagging dislocations away from the interface within one single framework. Our major findings are: (a) the interface dynamically responds to a pileup by forming steps/ledges, the height of which is proportional to the number of dislocations arriving at the interface; (b) when the pre-sheared sample is compressed, a direct square-to-hexagonal PT occurs ahead of the pileup tip and eventually grows into a wedge shape; (c) upon a further increase of the loading, part of the newly formed hexagonal phase transforms back to the square phase. The square product phase resulting from this reverse PT forms a twin with respect to the initial square phase. All phase boundaries (PBs) and twin boundaries (TBs) are stationary and correspond to zero thermodynamic Eshelby driving forces; and (d) the stress intensity induced by a pileup consisting of 16 dislocations reduces the stress required for initiating a PT by a factor of 5.5, comparing with that in the sample containing no dislocations. This work is the first characterization of the behavior of PTs/twinning resulting from the reaction between a microscale dislocation slip and an atomically structured interface.

Here, a focus is placed on metals in which the viscoplastic deformation can be conditioned by the interplay between all three deformation modes: dislocation motion, martensitic PTs, and twinning. In such scenarios, the coupling between these deformation modes will dictate the material microstructure evolution, such as PTs accompanied by the dislocation slip [12,13], and the PTs promoted by the dislocation-or twin-mediated plastic flow [8,14]. Several attempts on the study of the complex interaction between the dislocation-, twinning-, or PT-mediated plasticity include: (i) a simultaneous treatment of PTs, dislocations, and twinning within a continuum sharp interphase approach from the theoretical point of view in [15][16][17] and the numerical point of view in [18,19]; (ii) a smeared continuum description of the interaction of plasticity and PTs under high pressure; [20][21][22] (iii) continuum-level computational analysis of the role of interfaces and twinning on the PTs through crystal plasticity finite element (CPFE) simulations, where twinning is found to promote the α to ω PTs but suppress the reverse PTs in zirconium [8]; (iv) phase field approach (PFA) consideration of discrete dislocations and PTs at the nanoscale in [23][24][25][26][27] and also at the microscale in [7,28]; and (v) atomistic simulations of interaction between PTs and plasticity in [3,[29][30][31]. For a more comprehensive summary of the research progress in this field, a few representative review articles on this topic can be found in [32][33][34][35] and in [11,[36][37][38] from the materials and mechanics point of view, respectively. It should be noted that, in the above literature, twinning in martensite is often considered as one part of the multi-variant martensitic PTs.
In this work, we aim to probe the mechanisms underlying the interaction between dislocation slip, PTs, reverse PTs, and twinning. A co-operation of these three deformation modes is pertinent to materials subjected to extreme mechanical environments (e.g., high pressure, high torsion, shock loads, etc.) [11,[39][40][41][42][43]. Especially, the occurrence of dislocations and twinning may have a potential to decrease the critical pressure required for initiating PT. For example: (i) an irreversible PT from rhombohedral to superhard cubic boron nitride (BN) was obtained in [44] at 5.6 GPa under plastic straining but at 55 GPa under hydrostatic loading. A plastic straininduced PT from hexagonal to superhard wurtzitic BN can be induced under a pressure as low as 6.7GPa in a rotational diamond anvil cell (RDAC) while under hydrostatic loading it was not observed even at 52.8 GPa [45]; (ii) the plastic strain-induced PTs from graphite to hexagonal and cubic diamonds were observed in RDAC at 0.4 and 0.7 GPa, respectively, while under hydrostatic conditions they occurred at 20 GPa and 70 GPa, respectively [46]. Thus, the PT pressure reduction induced by a plastic straining may reach up to two orders of magnitude. Other than contributing to the PT pressure reduction, one more surprising effect is that the plastic straining under pressure may lead to new (hidden) phases that were not or could not be reached under hydrostatic conditions [11, 39-42, 47, 48].
In order to understand the dislocation-PT interaction under high pressure, the fundamental difference between the plastic strain-induced PTs and the pressure-or stressinduced PTs was introduced [37,49]. In details, the stress-induced PT refers to the PT under stresses below the yield strength. It starts at the pre-existing defects. In contrast, the strain-induced PT occurs at the defects generated during the deformation. Particularly, the dislocation pileup at the grain boundaries (GBs) is considered as the strongest internal stress concentrator responsible for the strain-induced PT [37,49]. As stated above, the plastic strain-induced PTs can occur at pressures one to two orders of magnitude lower than that under hydrostatic conditions. It thus requires a completely different thermodynamic/kinetic description and experimental characterization. The analytical solutions in [37,49] derived from a model by treating the dislocation pileup as a super-dislocation shows that, because all the stress components at a pileup tip are proportional to the number of the dislocations in a pileup, the PT pressure can be reduced by a factor of 10 and more. This conclusion was elaborated within a nanoscale PFA of coupling the evolution of discrete dislocations with PTs in [25][26][27][50][51][52] and a microscale (or scale-free) PFA simulating the simultaneous shear (dislocation or twinning) and PT in a polycrystalline aggregate under compression and shear [7,28].
Despite the success in providing explanation for the plastic strain-induced PT promotion, the above analytical treatments and nano-/micro-scale phase field simulations have a clear limitation: the GBs at which the dislocations pile up are modeled as fixed (with respect to material) lines, which cannot resolve the actual processes of interaction between dislocations and GBs. In order to capture more details about the slip-interface reaction and the subsequent structure changes, an atomistic resolution at the slip-interface intersection is necessary. Otherwise, the dislocation pileup's contribution to the subsequent PTs or twinning near the interface might be either under-or over-estimated. For example, our previous molecular dynamics (MD) simulations demonstrated that shuffle screw dislocations transmit through tilt GBs in silicon (Si) [53]. They simply cannot pile up and thus may not contribute much to a PT in Si. By contrast, 60 o shuffle dislocations pile up at GBs and significantly promote the amorphization [54]. Also, quite often the dislocations are located in the high-pressure phase and piled up against an immobile PB, causing PT within the remaining low-pressure phase. Such a situation has not been resolved at the atomic level yet due to difficulties of avoiding the interface motion.
Clearly, it remains a challenge using single-scale techniques/methodologies to address the full complexity associated with the microscale dislocation pileup and the subsequent structure changes because: (a) a dislocation pileup may be tens of micrometer in length and introduce a high internal stress field spanning tens of microns away from the slip-interface intersection [55,56]; (b) the PT at the tip of a dislocation pileup compete with slip transmission, twinning [8,[57][58][59], crack initiation [60][61][62], which all originate from the atomic scale and need to be resolved at a high resolution. Here we present a concurrent multiscale computational framework to address this challenge with a focus on understanding how a microscale dislocation pileup interacts with an atomically structured interface and how it contributes to the subsequent direct PTs, reverse PTs, and then twinning. It bridges the length scale gap between atomistic and continuum, and provides us with a platform for understanding the interactions between dislocation slip, PTs, and twinning from the bottom up.

II. METHODOLOGY
A concurrent atomistic-continuum (CAC) approach [63][64][65][66][67][68][69][70][71][72][73][74][75][76] built upon a formulation [77][78][79][80][81][82][83] that unifies the atomistic and continuum description of materials within one framework is deployed here. This formulation [80][81][82][83] is a generalization of the Irving-Kirkwood procedure [84][85][86][87] in statistical mechanics. It views the solid material as a collection of lattice cells continuously distributed in space, within each of which a group of discrete atoms is embedded. The continuum-level physical quantities, including mass density, linear momentum density, energy density, momentum flux (also referred as stress in continuum mechanics), and energy flux, are then defined in terms of the atomic positions, velocities, and interatomic forces through Dirac or Gaussian distribution functions [79,88]. An introduction of these physical quantities into the classical Newtonian mechanics leads to a series of equations, i.e., mass conservation, momentum balance, and energy conservation equations, which can govern the behavior of materials by considering them as a collection of atoms [80]. These equations are partial differential equations in the same form as the balance equations in classical continuum mechanics but with atomistic information being built-in. Thus, those equations can be solved using numerical techniques, such as finitedifference or finite element (FE), which are commonly used for solving the equations in continuum mechanics. The CAC simulation tool is an FE implementation [73,74] of this formulation. Unlike the FE model in classical continuum mechanics, which considers the material as a collection of mass points without any internal structures, the material body in CAC is discretized into finite number of elements, each of which is a collection of the lattice cells with the atomistic information (crystal structure, slip planes, cleavage planes, interstitial sites, etc.) being embedded. Furnished with such atomistic information, comparing with the other concurrent multiscale models, CAC has several unique features: (i) it does not need additional constitutive rules and can be directly driven by the traditional interatomic potentials or the machine learning-based potentials to be trained from ab initio calculation data; (ii) it does not need any special treatments to describe discontinuity, such as dislocations, PTs, or cracking, at the continuum level but with their atomistic nature being retained; (iii) it also does not need any sophisticated rules for passing the dislocation-mediated plastic flow from the atomistic domain to the continuum domain and vice versa [63][64][65]. At a fraction of the cost of full MD, CAC simulations have been performed to study: (a) the dislocation nucleation and growth in Cu, Al, Ni, and Si [66]; (b) the Si-I → Si-II PT in Si [74]; and (c) the dislocation transmission across a twin boundary (TB) in bi-crystalline metals [65]. It has a predictive capability at approximately the same level as that of MD but demands significantly less computational resources. This thus provides us with an ideal platform for understanding the microscale dislocation slip and the subsequent PTs/twinning nearby an atomically structured interface.
Here, as a first demonstration of its applicability in solving such problems, a two-dimensional (2D) two-phase material (the hexagonal and square phases co-exist with an incoherent interface in between) under compression and shear is selected as a model system because: (i) the 2D set-up enables a direct comparison of the simulation results with the analytical solutions under plane strain conditions for validation purpose; (ii) if desired, the interface structure in such systems can be manipulated to model the interfaces in a few recently developed highperformance metallic composites, such as Ti/Al [89], Mg/Nb [90], Cu/Nb [91], and among several others. (iii) the dislocation-induced local stresses in many realistic 2D materials (crystalline solids consisting of a single layer of atoms), such as one single layer of graphene sheet, boron nitride [92][93][94][95], or colloidal crystals [96,97], may have significantly contributed to the PTs in them but has not been fully understood yet.
The rest of this paper is structured as follows. In Section III, we briefly introduce the interatomic potential, the crystal structure, the PT variants, the boundary conditions, and also the loading strategies. The atomistic together with the multiscale computational analysis on the dislocation pileup formation process, the pileup-induced stress accumulation, as well as its role in the subsequent PTs, reverse PTs, and twinning are then presented in Section IV. Thereafter, we conclude this paper with a summary of our major findings as well as a brief discussion of future research in Section V.

A. Interatomic Potential and PT Variants
For the chosen material system, the atomic interaction in both square and hexagonal lattices is described using a modified Lennard-Jones (L-J) potential by Lee [98] as shown in Eq. (1): which was originally proposed to study the bcc-to-hcp PT in iron. The modified functional form of the L-J potential in Eq. (1) is the combination of a traditional L-J (12-6) term and an inverse Gaussian term. The 12-6 term leads to a stable hexagonal crystal structure at zero stress. The addition of an inverse Gaussian into the L-J potential gives rise to a square lattice. Here the parameters H, r mh , and σ h in Eq. (1) were chosen to stabilize: (i) an incoherent interface between the square and the hexagonal phase, which acts as an obstacle to the dislocation motion; and (ii) the core structure of a dislocation in the hexagonal phase. The parameters satisfying these two conditions are listed in Table I. In the absence of chemical heterogeneity, the atomic structure of materials with certain crystallographic configuration, such as fcc, bcc, hcp, diamond, zinc blend, and so on, can be determined by one single parameter: the potential energy per atom. The calculation of the potential energy as a function of the interatomic separation, i.e., the E − r relation, is often used to provide researchers with important information, such as: (i) the (a) The potential energy surface for 2D square and hexagonal lattices with the atomic interactions in them being described by a modified Lennard-Jones potential in Eq. (1); (b) two crystallographically equivalent variants (F + t and F − t ) of hexagonal phase (solid lines) and their orientation with respect to a square phase (dash lines).
r at which E becomes minimal is the equilibrium distance, r 0 , between atoms in the ground state; (ii) the pressure-induced PT from one phase to another phase occurs along the common tangent line of the E −r curves of those two phases. According to the potential energy functional form in Eq. (1) and its parameters in Table I, the calculated E − r curve is presented in Fig. 1a. Three key observations from Fig. 1a are: (1) this potential energy landscape has a minima at 0.595Å and 1.525Å for the hexagonal and square phases, respectively; (2) the square phase should be the ground state because it has a lower potential energy minimum than that of the hexagonal phase; and (3) the potential energy difference between hexagonal and square phases under zero stress is ∆Ψ = 0.066 eV. Although the square-to-hexagonal PT occurs under non-zero stresses, our simulations confirm that the hexagonal phase itself is metastable and can exist under zero stress. In this way, the chosen potential indeed enables us to stabilize two different lattices within one model. Figure 1b shows the unit cell of the square (dash lines) and the hexagonal (solid lines) phases with their lattice parameters being noted as a sq and a hex , respectively. In a 2D set-up, the transformation deformation gradient tensor F t , the stress and strain tensors can be all written as 2×2 matrices in a cartesian system with orthogonal unit basis vectors. Clearly, there are two crystallographically equivalent variants for the hexagonal phase due to the lattice symmetry. These two variants are in a twinned configuration with respect to each other. According to the lattice cell set-up of these two phases in Fig. 1b, the transformation deformation gradient associated with those two variants resulting from the square-to-hexagonal PTs are: The transformation deformation gradient can be then decomposed into orthogonal and symmetric parts through Here the superscript T means transposition. As such, the symmetric (rotation-free) transformation deformation gradient for these two variants can be calculated as

B. The Computer Model Setup
In order to determine the role of a dislocation pileup in the subsequent structure changes near the interface, the CAC model for a two-phase material system with square and hexagonal lattices co-existing is constructed (Fig. 2a). Clearly, when this system is subjected to loading, the deformation behavior near the square/hexagonal interface is critical and may dictate the material's overall microstructure evolution. The material domain near the interface is thus resolved at a fully atomistic resolution (Fig. 2a). In contrast, away from the interface, a CG description of the hexagonal phase using FE is deployed (Fig. 2b), which has much less degrees of freedom (DOF) than that of a fully atomistic model. To study the slipinterface reaction, similar to our previous work [63][64][65], the FEs in the CG domain are carefully aligned with their boundaries along a dislocation slip plane in the hexagonal lattice. Here, each FE contains sixty-four atoms and can slide with each other along the element boundaries. In this way, the dislocation-mediated slip far away from the interface can be accommodated in the CG domain at a fraction of the cost of MD. More importantly, since each FE in Fig. 2b is a collection of lattice cells, the forces acting on the FE nodes are calculated by converting the interatomic forces into the internal force density using a Gauss quadrature scheme elucidated in [73,74]. As such, the constitutive rule for the atomistic and CG domains in the CAC model is the same, i.e., the interatomic potential in Eq. (1). This differs from other continuum approaches which accommodate dislocations through the deployment of either a contact model in [7,28], a Heaviside step function in [99], or an additional DOF in CPFE [100,101].
Here we have carefully chosen the crystallographic orientation of both the hexagonal and square lattices to construct the incoherent interface ( Fig. 2b) such that: (a) the initially introduced dislocation slip is perpendicular to the interface, and (b) a dislocation transmission across the interface is suppressed, which has been justified through a detailed Schmid factor analysis and geometric compatibility factor analysis in [102]. In this way, a large number of dislocations can be piled up at the in-terface and a high internal stress is generated ahead of the slip-interface intersection, whose contributions to the subsequent PTs can be then well quantified without the need of considering transmission, cross-slip, and among several other complexities. Although the interface here satisfying the above two conditions may not be the one with an energy minimum, in practice, the majority of interfaces found in materials will not be at a minimum energy state anyway due to the presence of defects.
A queue of dislocations with a uniform separation of 12 nm in between are initially introduced into the hexagonal lattice, which is now in a CG description (Fig. 2b). As demonstrated in our previous work [63][64][65], because the CG domain in a CAC model can accommodate dislocations without smearing out their atomistic natures, the introduction of a dislocation into the CG domain is realized here through displacing FE nodes according to the analytical solution for the displacement field around an edge dislocation derived from the theory of elasticity. Different from traditional FE models in which the neighboring elements are connected by sharing the FE nodes, the FEs in CAC models across the slip plane are disconnected (See Fig. 2b). Thus, the displacement jump induced by a dislocation along the slip plane will be allowed. It should be also noted that the dislocation displacement fields, strain and stresses resulting from present CAC simulations will deviate from the linear isotropic elasticitybased solution. Instead, the actual atomic-level information, including crystal anisotropy, finite strain, nonlinear elasticity, and the interface structure heterogeneity will be all considered. The reason that we initially introduce dislocations through displacing the FE nodes according to the elasticity-based solution is to accelerate the convergence of our calculations by reducing the time required for equilibrating the system.
Due to the deployment of a CG description in the region away from the interface, the hexagonal lattice dimension, noted as L x−hex in Fig. 2a, in the CAC model can be at the micrometer level and even above. It enables us to initially introduce tens of dislocations into the model. This goes beyond the reach of a traditional MD model because MD usually has a limited length scale at nanometers, which can only accommodate a few dislocations in the pileup [103]. In the present simulations, on the other side of interface ahead of the pileup tip, the crystal structure is in a square lattice with a dimension of L x−sq along x direction. The sample dimension along y direction is chosen as L y = m hex a hex = n sq a sq , where m hex and n sq are the numbers of the hexagonal and square lattice cells along y direction, respectively. In addition to microscale CAC models, a series of nanoscale MD simulations are also performed using LAMMPS [104]. In all these simulations, the displacement along z direction of the sample is constrained to achieve a plane strain condition. The initial crystallographic orientations of square and hexagonal lattices are indicated in Fig. 2b. The dimensions and the number of the degree-of-freedom (DOF) in both nanoscale MD and microscale CAC models are listed in Table II. The left and right ends of the sample along x direction are free. A homogeneous shear within the xy plane is imposed until a desired shear stress is achieved. Thereafter, the top and bottom ends of the sample along y direction are fixed for equilibrating the dislocation configuration in the pileup. When a PT needs to be activated, a displacement controlled compressive loading along y direction will be then applied on the two ends of the sample.

A. Dislocation Pileup Process
The computational models containing the built-in dislocations are firstly relaxed for a duration of 20 ps with a time step of 1 fs to achieve a stable configuration. Thereafter, a shear stress, noted as τ ap , is imposed on the whole samples to drive the dislocations towards the interface. A certain level of τ ap is realized through applying a homogeneous shear strain while monitoring the resulting shear stress on the fly. If the shear stress does not match the desired τ ap , the shear stress is adjusted by increasing or decreasing the applied shear strain. When τ ap arrives at the desired value, the top and bottom boundaries of the sample are constrained and not allowed to move along both x and y directions any more.
An internal stress concentration, noted as τ xy , ahead of the slip-interface intersection is generated due to the large number of dislocations' arrival at the interface. Fig.  3a and Fig. 3b show the contour of the internal stress field obtained from CAC and MD simulations of 16 and 8 dislocations piling up at the interface under a shear of τ ap = 6 × 10 7 N/m and τ ap = 2 × 10 7 N/m, respectively. It should be noted that the contour here shows the magnitude of τ xy − τ ap , rather than the absolute value of τ xy . It is seen that the level of the internal stress concentration ahead of a pileup tip in Fig. 3a is significantly higher than that in Fig. 3b, because more dislocations have participated in the formation of a pileup in CAC (Fig. 3a) than that in MD (Fig. 3b).
With an atomistic resolution at the interface, CAC provides us with an opportunity of examining the atomiclevel structure evolution at the slip-interface intersection. Figure 4 presents the snapshots showing the atomic structure configuration (red: square lattice; blue: hexagonal lattice; green: defects) when different number of dislocations arrive at the interface. It is seen that, when τ ap is at a level of 1.3 × 10 7 N/m, one dislocation is behind the pileup tip and two dislocations arrive at the interface, a step with a height of 5.12Å, i.e., the magnitude of two Burgers vector, 2b, is formed (Fig. 4a). This step can be approximately considered as a single super-dislocation with a Burgers vector of N a b, where N a is the number of the dislocations arriving at the interface. Upon a further increase up to τ ap = 7.1 × 10 7 N/m, 9 dislocations arrive at the interface and eventually, leading to the formation of a step with a height of 9b (Fig. 4b). Under a shear stress of τ ap = 10.6 × 10 7 N/m, 15 dislocations arrive at the interface, Fig. 4c show the corresponding atomic arrangements. At this stage, approximately, the step at the interface has a height of around 15b.
As designed, the interface under consideration here indeed blocks the motion of dislocations without allowing any transmission. Also, as shown in [102], the dislocation configuration behind the pileup tip in Fig. 3 actually matches an analytical solution from Hirth and Lothe very well. More importantly, the step-like super-dislocation with a complex local structure is formed at the slip-interface intersection. Due to the local atomic structure reconstruction, this super-dislocation carries a sophisticated core structure rather than a simple additive of multiple dislocation cores. Similar features were observed in our previous full MD simulations of dislocation pileup against a tilt grain boundary in silicon, in which the amorphization has been found ahead of a slipinterface intersection [54]. It should be also noted that the presence of a thermal bath would likely promote the re-arrangement of the atoms at the interface and may promote a dislocation transmission. Without including any thermal fluctuations, the present simulations should be considered as extreme cases for the slip-interface reaction. A resolution of the thermal-induced atomic structure reconfiguration during the slip-interface reaction at finite temperature is beyond the scope here and will be studied in our future work.

B. Pileup-Induced Stress Intensity Factors
To characterize the stress profile ahead of the pileup tip, we construct a series of finite-sized volume elements ahead of the slip-interface intersection (see the inset picture in Fig. 5). Each volume element is at a resolution of 5Å× 5Å. The stress associated with each volume element can be then calculated using the Virial formula. Figure 5 presents the CAC-and MD simulation-predicted shear stress distributions ahead of the dislocation pileup tip. Two major observations are: (a) both CAC and MD simulations predict an obvious stress concentration at the pileup tip and its rise upon an increase of the applied shear, τ ap ; (b) this stress concentration decays away from pileup tip but spans a longer range (100nm, Fig. 5a) in CAC than that (60nm, Fig. 5b) in MD because more dislocations have been piled up in the CAC model. It suggests that, similar to what has been observed in experiments [55], the participation of a large number of dislocations in a slip may produce a long-range internal stress field, which can span tens of microns away from the pileup tip if tens of dislocations are piled up.
The obtained stress profile ahead of the pileup tip can be then fit into an equation as follows: In this equation, the dislocation pileup-induced internal stresses is, for the first time, decomposed into two parts. pileup tip; r is the distance between the pileup tip and the stress measurement site; K b is the stress intensity factor induced by the dislocations behind the pileup tip; z is the correction coefficient due to finite strain, nonlinear anisotropic elasticity, finite step size, and two-phase material.
A determination of those parameters is accomplished through fitting Eq. (4) into the stress profiles from two sets of simulations: (i) in simulation set-1, for a sample containing N dislocations, we apply a stress until all dislocations arrive at the interface and form a step, i.e., N a = N and N b = 0. Here, N b stands for the number of the dislocations behind the pileup tip (see inset picture of Fig. 6 for the physical meaning of N , N a , and N b ). When N a = N and N b = 0, the first term on the right side of Eq. (4) simply goes to zero. A fitting of simulation-based stress profile into Eq. (4) leads to a determination of the value of z, τ 0 , and r 0 ; (ii) in simulation set-2, a stress lower than that in set-1 is imposed on the same sample containing N dislocations. In this situation, N b = 0 and N a = N − N b . As such, the first term on the right side of Eq. (4) is non-zero. The K b in it can be then deter-mined by fitting the stress profile from simulation set-2 into Eq. (4). It should be noted that, at this step, the values of z, τ 0 , and r 0 fitted from simulation set-1 are used. In the fitting process, we only consider the stress profile up to r = 67nm where the internal stress does not relax due to the presence of free surfaces. Two different sample sizes are considered here. One is of a dimension of 490 nm × 84 nm and contains 8 dislocations. The other is of a dimension of 1.7 µm ×138 nm and contains 16 dislocations.
The τ ap − N a and τ ap − N b relations obtained from the MD together with CAC simulation of the dislocation pileup process are presented in Fig. 6a. It is seen that, in all simulations, when a pileup happens upon an increase of τ ap , N a increases and N b decreases. Fig. 6b presents the simulation-predicted relationship between K b and N b . Clearly, K b linearly increases upon an increase of N b . In addition to K b , another parameter that we have introduced into Eq. (4) is the correction coefficient z, which is needed here because: (1) the classical super-dislocation model only considers the defect embedded in an isotropic medium with a constant shear modulus of µ while we have considered the super-dislocation at a PB in a two-phase anisotropic material; (2) during the dislocation pileup process, both physical and geometric nonlinearities develop; (3) additivity of stress field from each dislocation is violated. Especially, when more than 8 dislocations arrive at the PB, the slope of z − N a relation is significantly larger than that when N a is less than 8 (Fig. 6c). This in turn, implies that the stress intensity factor upper bends to a very high level when tens of dislocations participate the formation of a pileup at the microscale, comparing with that when only several dislocations are involved in the formation of a pileup at the nanoscale. Such a a highly non-linear z − N a relation can be fit into an equation as: It should be noted that the data in Fig. 6c exhibits a trend of deviating from Eq. (5) when N b is larger than 2, especially when N = 16, N a = 12, N b = 4 or when N = 16, N a = 10, N b = 6. We believe this deviation is caused by a fundamental difference between a superdislocation model and the computer set-up for non-zero N b . A super-dislocation model, i.e., the second term on the right side of Eq. (4), is only exact when N b = 0, i.e, all the dislocations arriving at the interface participate the formation of a step. When N b is larger than 2, one can't rely on the super-dislocation model to estimate the internal stresses ahead of a slip-interface intersection any more. To confirm this assertion, we perform two additional CAC simulations. In one simulation, N = 12, N a = 12, and N b = 0. By contrast, in another simulation, N = 10, N a = 10, and N b = 0. Comparing with the data from simulations with N = 16, N a = 12, and N b = 4 or N = 16, N a = 10, and N b = 6, the results (pink circles) from these two new sets of simulations with N b = 0 obviously fits into the Eq. (5) with a significantly less error. Other than z and K b , another key parameter in Eq. (4) is r 0 , which indicates the location of the maximum stress concentration. This location obviously changes during the pileup process. r 0 is thus introduced in Eq. (4) to capture the such local structure evolution at the pileup tip. As expected, the values of r 0 obtained from fitting simulation data into Eq. (4) are found to be linearly proportional to the number of the dislocations arriving at the interface (Fig. 6d) with the slope equals to b/2, b as the length of the Burgers vector. This means that r 0 is located in the middle of a step produced by dislocations at the interface.
C. Pileup-assisted PTs and Twinning

The two-variant hexagonal phase formation
Thereafter the dislocation pileup formation, a compression along y direction is applied through a displacement-controlled boundary condition. Together with the pre-shear τ ap , the top and bottom boundary layers of the sample are held as rigid and forced to move towards each other at 0.2 m/s. This allows us to determine the role of a dislocation pileup in the subsequent structure changes, such as PTs, reverse PTs, and twinning, if there would be any.
Firstly, the dislocation pileup-induced strain localization together with the local atomic structure changes are analyzed by calculating the atomic-level deformation gradient tensor, F , the Lagrangian strain tensor E = 1/2(F T · F − I), and also the coordination number using OVITO [105]. The time sequences of the snapshots showing the atomic structure evolution during the PT process from CAC simulation (7 dislocations have arrived at the interface and another 9 of them are behind the pileup tip) are displayed in Fig. 7. The atoms in the first row of Fig. 7 are colored in a shear strain component of E xy and the atoms in the second row are colored according to their coordination number with a cutoff of 3.92Å. The coordination number analysis in Fig. 7 clearly shows that, at an early stage of the compression (from t = 12 ps to t = 15 ps), there exists a direct PT, i.e., from square (red) to hexagonal (blue), ahead of the dislocation pileup tip. During such direct PTs (square-to-hexagonal), the PT occurs through simultaneously activating two variants of the hexagonal phase as shown in Fig. 1b. When these two transformed domains are in contact with each other, a twin boundary (TB) is formed. In details, at t = 12 ps, two atomic layers ahead of the pileup tip in the square phase are sheared to form a stacking fault, a precursor of the hexagonal phase. With a further increase of the loading, the stacking faulted atomic layers transforms to the hexagonal phase at t = 14 ps. The shear strain associated with the newly formed hexagonal phase above the slip plane is in an opposite sign (blue) com- paring with that (orange) below this plane (see the first row of Fig. 7 from t = 14 ps and t = 17 ps). This corresponds to the two variants in the newly formed hexagonal phases, which are twinned with respect to each other and separated by a TB in between. Upon a further increase of the loading, the newly formed hexagonal phase grows into a wedge shape. The interface between the product (hexagonal) and the parent (square) phases consist of horizontal coherent regions separated by steps, avoiding the high-energy irrational interfaces.
Here the CAC simulation results, especially the material microstructure evolution ahead of the dislocation pileup, are also compared with that from MD simulations. In details, Fig. 8a shows the snapshots of atomic structure evolution from MD simulations of the deformation behavior in a sample under a pre-shear of τ ap = 2 × 10 7 N/m together with a compressive strain of ε y being increased from 0.054 to 0.080. In parallel, Fig. 8b presents the results from a coordinate number analysis of the CAC simulations under τ ap = 4 × 10 7 N/m together with ε y increasing from 0.042 to 0.063. Fig. 8c shows the results from CAC simulations of the dislocation pileup-assisted PT in a micron-sized sample under τ ap = 6 × 10 7 N/m with ε y being increased from 0.028 to 0.039. In Fig. 8a-c, the atoms in hexagonal and square phases are colored in blue and red, respectively. The atoms in any other configuration different from hexagonal and square lattice are colored in green. They represent either the PBs that separate the hexagonal and square phase, or the stacking faults/defects in the hexagonal phases. Similar atomic structure evolution in the same region from MD (Fig. 8a) and CAC (Fig. 8b and c) simulations is also analyzed in terms of atomic potential energy (PE) as shown in Fig. 9, where the atoms in hexagonal and square phases are colored in blue and orange, respectively, while those atoms with PE values deviating from the hexagonal or square lattice are in green.
Our major observations from Fig. 8 and Fig. 9 are: (1) the square-to-hexagonal PT occurs at ε y = 0.028 when τ ap = 6 × 10 7 N/m is applied on a micron-sized CAC model where tens of dislocations participate the pileup, while ε y = 0.054 when τ ap = 2 × 10 7 N/m is applied on nanoscale MD models where only a few dislocations participate the pileup; (2) in both CAC and MD simulations, the PT starts at the dislocation pileup tip. The newly nucleated hexagonal phase then grows into a wedge shape. Such a wedge-shape in the product phase resulting from a dislocation-assisted PT has also been observed in experiments [106], although in a different material; (3) the width of the interface between the parent and the product phases can be approximately estimated to be about 1 nm in both CAC and MD simulations; This value can be used for calibrating the PFA models in [50][51][52]; (4) the stacking fault in the product hexagonal phase only appears in the microscale CAC simulations at τ ap = 6 × 10 7 N/m, suggesting an increased microstructure complexity when more dislocations are piled up.

Twinning in square phase via a reverse PT
Interestingly, at a later stage of the deformation, the coordination number analysis in the second row of Fig.  7 when t = 15 ps and t = 17 ps shows that, at the slip- interface intersection, the newly formed hexagonal phases ahead of the pileup tip transform back to the square phases. This is referred as a reverse PT, i.e., hexagonalto-square PT. The finish of the square-to-hexagonal PT and then the onset of a reverse PT result from the large back stresses due to transformation shear and an atomic structure reconfiguration at the slip-interface intersection, which has largely released the local internal stresses. The newly formed square phase resulting from such a reverse PT together with the original square phases form a twinned structure. Although the twinning formation through a direct PT was reported in literature [107][108][109][110], our finding here is believed to be the first direct prediction of twinning formation during a monotonous loading through a dislocation pileup-induced PT and then a reverse PT. In other words, twinning does not occur directly in the square phase due to a high energy barrier there, but occurs via a direct-to-reverse PT instead. Analytically, when the elastic strains are neglected, the transformation deformation gradients, U i and U j , in the two regions across a twin boundary (TB) satisfy Here, γ t is the twinning shear, m is the twinning direction within the twinning plane which has a unit normal of n. Eq. (6) is used here to check whether the newly formed square phase resulting from the reverse PT is indeed twinned with respect to the original square phase. For the observed reverse PT, i represents an original square phase, and j is for the domain in which the material is in a square phase resulting from the hexagonal-to-square PT. Then Q i = U i = U j = I, and Eq. (6) is simplified as which has at least two solutions as following: γ + = 0.25; or To compare our simulation results against the above analytical solutions, the atomic-level deformation gradient components, F xx , F xy , F yx , and F yy , from CAC simulations at t=17 ps are calculated using the initially undeformed atomic configuration as a reference and shown in Fig. 10a. It is seen that the simulation-based F in the region where the reverse PT occurs agrees with the above analytical solution of Q very well. Therefore, the square phase resulting from such a reverse PT is indeed twinned with respect to the original square phase. In more details, for the region where the reverse PT occurs, the atomic configuration at the dislocation pileup tip is shown in Fig. 10b. Here all the atoms are displayed and colored by their deformation gradient, F xx , F xy , F yx , and F yy , respectively. Based on the detailed crystallographic orientation analysis (Fig. 10c), the TB formation process is also sketched in Fig. 10d for clarifications. During compression, accompanied by the direct and reverse PTs, the local stresses ahead of the dislocation pileup evolves. To characterize the full complexity of such a complex internal stress evolution, we measure the local stresses and again decompose it into two parts: in which K b and K a are noted as the intensity factors based on the Eshelby and super-dislocation model, respectively. During the process of the direct and reverse PT, a series of stress profiles ahead of the pileup tip are produced. A fitting of each stress profile into the above equation leads to a determination of the instantaneous K b and K a . A time evolution of K b and K a together with the snapshots of the atomic structure at the pileup tip are presented in Fig. 11. Obviously, before the occurrence of PT, both K b and K a are at a relatively high level. And the stress intensity, K a , induced by the dislocation at the pileup tip is larger than that, K b . Both K b and K a decrease when PTs occur. In details, at t = 13 ps, with the growth of the hexagonal phase, K a decreases from 165 to 80 while K b decreases from 100 to 75 (all in a unit of 10 7 N/m 1/2 ). It suggests that a direct PT largely relaxes the local stress concentration. Thereafter, the step at the slip-interface intersection is found to disassociate into multiple mini-steps. At this stage, the reverse PT starts and those mini-steps eventually emerge as a single step at t = 15 ps. Up to this moment, the stress intensity becomes low. At t = 14 ps, there is a cross-over of the stress intensity induced by the dislocations behind the pileup tip: K b becomes larger than that, K a , induced by the dislocations at the pileup tip.

Thermodynamic driving force at the interfaces
In addition to the atomic-level deformation gradient tensor as shown in Fig. 10, for the snapshot at t = 17 ps in Fig. 7, three components of the Lagrange strain tensor, E xx , E yy and E xy , associated with each atom ahead of the slip-interface intersection are also extracted and presented in Fig. 12. In order to quantify how the local stress contributes to the subsequent PTs, an atomic-level local strain together with a local stress-based thermodynamic driving force analysis (similar to that for nanoscale [51,52] and scale-free [7,28,111] PFAs) are also performed here. Based on the atomic-level local stresses, the thermodynamic Eshelby-type driving force for PT and twinning, X is defined as: Here, P = Jσ · F −1T is the first Piola-Kirchhoff stress expressed in terms of the Cauchy stress σ with J = detF . ∆Ψ is potential energy difference between the stress-free phases from both sides of each interface (which is zero for twin boundaries). F t1 and F t2 are the transformation deformation gradients in the two material domains across an interface. In particular, for the TBs in the hexagonal phases, they are F + t and F − t in Eq. (2), respectively. For the TBs in the square phases ahead of the pileup tip, they are I and Q + or Q + in Eq. (8) and Eq. (9). For the PBs between hexagonal and square phases, they are either F + t or F − t and the transformation deformation gradient in the square phase, I, Q + or Q − . The condition of X = 0 corresponds to the local thermodynamic equilibrium of the relevant interface under consideration.
The black lines in Fig. 12c correspond to the phase interface equilibrium condition of X = 0 with the deformation gradients in it being taken from the parent and product phases across the interface. In contrast, the white lines correspond to the twin interface equilibrium condition of X = 0 with the deformation gradients in it being taken from the two twinning variants across the TB. Figure 12c shows that, during the processes of direct PTs, reverse PTs, and twinning, all interfaces (including initial vertical phase interface) correspond to the continuum phase equilibrium condition of X = 0. Thus, these interfaces are stationary. Although this has been shown for cubic-tetragonal phase interface in nanoscale [51,52] and scale-free [7,28,111] PFAs at the continuum level, it is observed for the first time here at the atomic scale with a variety of different twinning interfaces being included.

D. The Role of the Pileup-Induced Local Stresses in the Subsequent PTs
In order to quantify how the dislocation pileup-induced local stresses contribute to the subsequent PT, we per-form a series of analysis on how the square phase responds to compression through plotting the (σ y -ε y ) curves (Figs. 13a-b). Here, ε y is the Lagrangian strain along y-direction. The stress, σ y , is an average of the true stresses along y-direction acting on the atoms located in a domain (100 nm × 22 nm in CAC and 60 nm × 15 nm in MD, respectively) ahead of the dislocation pileup tip. Results in Figs. 13a-b show that: (i) prior to ε y = 0.02, both CAC and MD simulations predict a linear stress-strain relation with the same modulus of E. At larger strains, the material behaviour is highly nonlinear though; (ii) thereafter, a plateau appears on all the stress-strain curves; (iii) a square-to-hexagonal PT is found to start at ε y where the plateau appears. This corresponds well to the atomic configuration evolution in Fig. 8 and Fig. 9; (iv) the critical compressive stress required for the occurrence of PT can be thus identified as σ c in Figs. 13a-b; and (v) σ c largely decreases when the pre-applied shear is increased from τ ap = 2 × 10 7 N/m to τ ap = 8 × 10 7 N/m. The considerable reduction of σ c with the increase of τ ap , i.e., from σ c = 16 × 10 7 N/m at τ ap = 2 × 10 7 N/m to σ c = 4 × 10 7 N/m at τ ap = 8 × 10 7 N/m, can be attributed to the large number of dislocations accumulated at the interface, which has introduced a high local stress assisting the PT at the pileup tip.
Furthermore, Fig. 14 presents the data from both CAC and MD simulations for quantitatively relating the critical stress, σ c , at which the PT starts, with the previously determined stress intensity factor K b and K a . tion by a factor of 5.5; and (ii) the σ c − K a relation can be divided into to two stages. Upon an increase of K b and K a , σ c in Stage-I reduces considerably faster than it does in Stage-II. It should be also noted that, for the material system under consideration here, σ y is not the only stress component that contributes to the PT. Our recent simulations show that all the stress components, i.e., σ x , σ y , and τ xy , all contribute to PTs, although the magnitude of σ x in the present simulations is found to be considerably small. A crystal lattice instability criterion based on all components of the stress tensor will be reported in a separate paper.

V. CONCLUDING REMARKS
To summarize, here we present an atomistic-tomicroscale computational analysis of the interplay between dislocation slip and PT/twinning in two-phase materials under compression and shear. One main novelty of the CAC approach deployed here is its capability in bridging the relevant length scales by resolving the structure changes near an interface at the atomic scale while the lagging dislocations away from the interface in a coarse-grained atomistic description. It thus expands the MD-simulation-based predictive capability from the nanoscale to the micrometer level. The main findings of this study are summarized as below: (1) The micron-sized CAC model can accommodate up to 16 dislocations in one slip at a modest computational cost. These dislocations may be blocked by the obstacles (an incoherent interface in this work) and form a pileup spanning a range of several micrometers (1.2 µm in the present model). In contrast, the nanoscale MD model using the same computational resource can only accommodate up to 8 dislocations in a pileup, the equilibrium configuration of which under certain shear stress only spans tens of nanometers; (2) The internal stress intensity induced by a dislocation pileup can be decomposed into two parts. One (noted as K b ) follows the Eshelby model and is contributed by the dislocations behind the pileup tip, the other (noted as K a ) follows the super-dislocation model and is caused by the step formed at the slip-interface intersection. Either the Eshelby or the super-dislocation model alone will largely underestimate the dislocation pileup-induced internal stress concentration. The local internal stress intensity is linearly proportional to the number of the dislocations in a pileup at the nanoscale, but upper bends to a very high level when tens of dislocations participate the pileup at the microscale.
(3) When the pre-sheared material sample is subjected to a compression, PTs and twinning ahead of the dislocation pileup where the local stresses concentrate occur through a two-step process: (a) the square phase transforms to the hexagonal phase. During such a squareto-hexagonal PT (also referred as the direct PT), two variants of the hexagonal phase simultaneously nucleate, grow, contact, and form a twin boundary in between; (b) at a later stage, a portion of the newly formed hexagonal phase transforms back to the square phases.The newly formed square phase resulting from this reverse PT is found to be twinned with respect to the original square lattices and relax the local internal stresses. This is the first direct observation of twinning formation from a dislocation slip-assisted PTs and reverse PTs during monotonous loading.
(4) The Eshelby driving forces at all the newly formed PBs or TBs are zero, confirming that all of them are in the thermodynamic equilibrium. Although such analysis was performed previously in the nano-and microscale PFAs at the continuum level, it was performed for the first time for atomistic modeling. This can be an evidence of the consistency between the atomistic and continuum thermodynamic treatments in our model.
(5) The direct PT largely relaxes the local stress concentration due to the transformation strain. Thereafter, the step at the slip-interface intersection is found to disassociate into multiple mini-steps. The reverse PT then starts and the mini-steps eventually emerge as one step.
(6) The critical compressive stress, σ c , required for ini-tiating the PT, largely decreases with the increase of the internal stress intensity factor. Our simulations suggest a possibility of reducing σ c for PT by a factor of 5.5 and even higher through injecting a considerably large number of dislocations into a micron-sized material sample, as observed in experiments for various material systems.
These findings highlight: (a) the insufficiency of only using nanoscale MD simulations to interpret the experimental observations on the slip-interface reaction, which may have involved hundreds of µm-long dislocations participating in the formation of a pileup; and (b) the possibility of using CAC to predict how the microscale dislocation-mediated plastic flow reacts with the buried interfaces in a variety of multi-phase materials, such as fcc/bcc, fcc/hcp, bcc/hcp metallic composites, Ti-/Zr-/high entropy alloys, among others, when subjected to a severe deformation. In such scenarios, the CAC simulation tool may provide researchers with suitable, if not the best, vehicle for simultanesouly considering the microscale plasticity together with the atomic-scale interface structure relaxation. Nevertheless, this framework is considered to be still at a preliminary stage because: (i) Comparing with MD, the length scale of the present CAC models is indeed one step closer but is not the same level as that in experiments yet. In particular, as far as the phase growth after the slip-interface reactions are concerned, a further length scaling up to microns through a CG description of the PTs in the material domain far away from the interface is needed. Otherwise, the growth of new phases from CAC simulations is confined at the nanoscale; (ii) In addition to the local stress, another factor that plays a vital role in the PT process is the thermal-induced atomic fluctuations. The implementation of a finite temperature algorithm into CAC for capturing the finite-temperature effects on dislocations, PTs, and their interactions is needed, especially when the correlation between phonon instability and PTs becomes a concern; (iii) The material system under consideration here is oversimplified in terms of crystal structure, chemistry, interatomic potential, and microstructure. A transfer of the present model or the gained knowledge for understanding the dislocations-interface reactions and subsequent structural changes in realistic multi-phase materials is not trivial. It demands: (1) the design of new finite elements for simultaneously accommodating complex dislocation activities; (2) the implementation of more sophisticated or machine learning-based interatomic potential to be trained from ab initio data for capturing more phase variants, as well as (3) the incorporation of realistic structures and chemistry at the GBs, PBs, and other interfaces in the materials. An expansion of CAC along those three directions and its applications in predicting the slip-interface reactions in real materials are being attempted and will be reported in the future.