Slip-free multiplication and complexity of dislocation networks in FCC metals

During plastic deformation of crystalline solids, intricate networks of dislocation lines form and evolve. To capture dislocation density evolution, prominent theories of crystal plasticity assume that 1) multiplication is driven by slip in active slip systems and 2) pair-wise slip system interactions dominate network evolution. In this work, we analyze a massive database of over 100 discrete dislocation dynamics simulations (with cross-slip suppressed), and our findings bring both of these assumptions into question. We demonstrate that dislocation multiplication is commonly observed on slip systems with no applied stress and no plastic strain rate, a phenomenon we refer to as slip-free multiplication. We show that while the formation of glissile junctions provides one mechanism for slip-free multiplication, additional mechanisms which account for the influence of coplanar interactions are needed to fully explain the observations. Unlike glissile junction formation which results from a binary reaction between a pair of slip systems, these new multiplication mechanisms require higher order reactions that lead to complex network configurations. While these complex configurations have not been given much attention previously, they account for about 50% of the line intersections in our database.


Introduction
A fundamental goal in dislocation theory and physical metallurgy is to understand how and why dislocations multiply. Such multiplication affects numerous mechanical properties of crystalline solids, such as the strain hardening rate, fracture toughness, fatigue-life, and creep-life, to name a few. However, theories which predict the multiplication rate on the basis of fundamental dislocation processes have been elusive; existing models for dislocation multiplication are either phenomenological (i.e., not connected to the underlying physical mechanisms) (Mecking and Kocks 1981;Roters et al. 2010) or derived on the basis of assumed multiplication mechanisms, such as storage via junction formation Kubin et al. 2008a). In the pursuit of a physics-based model for the multiplication rate which minimizes constitutive assumptions, the authors recently compiled a database of discrete dislocation dynamics (DDD) simulations for more than 100 loading orientations in which dislocation cross-slip was suppressed (Akhondzadeh et al. 2020). While analyzing that database, it was discovered that dislocation multiplication frequently occurs on slip systems which experience zero applied shear stress (i.e., zero Schmid factor) and have a plastic strain rate of zero; we termed such multiplication slip-free multiplication. We subsequently sought to clarify the mechanisms underlying slip-free multiplication and assess its significance in greater detail. In so doing, we discovered that the dislocation networks in our simulations did not evolve in the simple manner assumed by predominant theories; these networks are tremendously complex and cannot be readily tied to common dislocation theory concepts, such as the formation of binary junctions. We further discovered that slip-free multiplication appears to be, in part, the result of these complex interactions. The goal of this work is to elucidate these discoveries and discuss their importance and implications.
The prevailing view regarding dislocation network evolution in face-centered cubic (FCC) metals is that the interactions within the network can be analyzed in terms of pairwise interactions between slip systems (Argon 2008; Franciosi et al. 1980;Hirth 1961). Each mobile dislocation belongs to one of the twelve slip systems of the 1 2 110 {111} type. Due to the symmetry of these slip systems, there are only six distinct types of interactions possible between every pair of slip systems; namely, self, coplanar, Hirth, Lomer, glissile and collinear (Kubin 2013). Of these, the interactions between non-coplanar (forest) slip systems (Hirth, Lomer, glissile, collinear) lead to the formation of binary junctions (e.g., see Fig. 6a). The stability or strength of the six slip system interactions affects the evolution of the network, thereby influencing the flow stress on each slip system (τ i ), the hardening rate ( ), as well as dislocation multiplication rate (ρ i ), where the subscript i denotes the slip system (Akhondzadeh et al. 2020).  and  used specialized simulations to compute the interaction coefficients (proportionality constant between the flow stress and the forest dislocation density) between non-coplanar slip systems and concluded that the collinear junction is the strongest type of interaction (Devincre et al. 2005), and hence is expected to play a dominant role in the yield strength. We computed interaction coefficients by analyzing our large-scale DDD simulation database and found that the glissile junction has the largest interaction coefficient (Akhondzadeh et al. 2020). We further observed that for a loading orientation to exhibit a high hardening rate ( > 100 MPa), the dominant slip system (slip system with the highest Schmid factor) needed to have either a collinear or glissile interaction with another active slip system.
While slip system interaction coefficients provide a measure of the strength contribution from a given interaction, they do not directly inform our understanding of dislocation multiplication mechanisms. Additional analyses or assumptions are necessary to understand multiplication. For example, Kubin et al. (2008a) developed a theory for dislocation multiplication based on the assumption that dislocations are "stored" when stable binary junctions form, invoking the notion of a dislocation mean free path . In a similar vein, Stricker and Weygand (2015) examined the role of glissile junctions during the evolution of dislocation structures. Glissile junctions are unique among the four junction types because the junction dislocation belongs to a 1 2 110 {111} slip system, so is able to bow out and contribute to plastic flow. Glissile junction formation can also be thought of as a multiplication mechanism, since the junction dislocation increases the dislocation density of the 1 2 110 {111} slip system to which the junction belongs; we refer to this multiplication mechanism as the direct glissile mechanism. Stricker and Weygand (2015) found that dislocation segments originating from glissile junction reactions contributed to 20-50% of the total plastic strain and dislocation density. Building on this work, Sudmanns et al. (2019) proposed a multiplication model which accounts for the direct glissile multiplication mechanism. While most of the previous studies were focused on non-coplanar (forest) slip system interactions, the interactions between dislocations on coplanar slip systems have received less attention, especially in terms of their contribution to multiplication. However, in our analysis below we find that coplanar reactions are associated with a key mechanism of slip-free multiplication. Hence, we expect that coplanar reactions also provide a major contribution to the overall multiplication rate.
Slip-free multiplication has been observed by several other research groups. The dislocation density increase in inactive slip systems has been observed in experiments (Higashida et al. 1986) and ultra-scale molecular dynamics simulations (Zepeda-Ruiz et al. 2020). In the DDD simulations of Weygand (2014) and Stricker and Weygand (2015), slip-free multiplication was attributed to the direct glissile mechanism. While the direct glissile mechanism plays an important role, we show below that it cannot explain all of the behaviors we observe in our simulation database. For example, in loading orientations such as [ 0 1 1] and [ 4 9 10] where slip-free multiplication is observed, an insufficient number of slip systems are active to produce glissile junctions on the slip-free system via a binary reaction. Instead, we show below that a combination of glissile and coplanar interactions are, at least in part, responsible for dislocation multiplication on certain slip-free slip systems.
In trying to better understand the dislocation processes which give rise to slip-free multiplication, we discovered that the networks from our simulations are highly complex. By complex, we mean that the topology of the network cannot be readily explained on the basis of binary collisions between pairs of slip systems. Instead, the nodes and links (dislocation lines which connect nodes) of the network are largely unrecognizable when viewed through the lens of binary reactions. Furthermore, we find that slip-free multiplication is intimately tied to this complexity. Below we demonstrate two slip-free multiplication mechanisms which involve "unconventional" network configurations. However, it is likely that more multiplication mechanisms are at play; the complexity of the dislocation networks makes it tremendously difficult to identify mechanisms. All together, these observations bring into question the traditional view that the dislocation network can be characterized in terms of slip system pair interactions, and motivate more in-depth analyses of the network and its evolution.
The remainder of the paper is organized as follows. In "Simulations setup" section, we present a summary of the simulation setup. In "Data analysis" section, we investigate how microstructural features of slip system i, are correlated with each other and with their values on other slip systems. We observe that the dislocation density ρ i values of two slip systems are correlated with each other as long as they are coplanar, even when their Schmid factors S i are not correlated with each other. A related finding is that there are many cases in which dislocation multiplication occurs on slip systems with very small Schmid factors and negligible shear strain rateγ i . In "Physical analysis of slip-free multiplication" section we examine the physical mechanisms of slip-free multiplication, showing that glissile and coplanar interactions both play an important role. In addition, we show that the actual sequence of events leading to dislocation multiplication can be more complex than what is assumed in the conventional view. In "Discussion" section, we discuss the physical consequences of slip-free multiplication, including a modification to the Kocks-Mecking model of dislocation multiplication. A summary is given in "Conclusion" section.

Simulations setup
In this work, we make use of our large DDD simulation database presented in Akhondzadeh et al. (2020) to investigate the mechanisms of dislocation multiplication. In what follows, a brief summary of the simulation setup is given. The DDD simulations were performed using the ParaDiS program (Arsenlis et al. 2007) with subcycling time integration algorithm (Sills et al. 2016) implemented for GPU computation (Bertin et al. 2019). Material properties for FCC single crystalline copper were used, with shear modulus μ = 54.6 GPa, Poisson's ratio ν = 0.324, and Burgers vector magnitude b = 0.255 nm. Glissile dislocations on the 1 2 110 {111} slip systems follow a linear mobility law with drag coefficient B = 1.56 × 10 −5 Pa·s. Cross slip was suppressed in all of the simulations (see Akhondzadeh et al. 2020).
The initial configuration was obtained by relaxing randomly introduced straight dislocation lines whose character angles are one of 0 • (screw dislocation), 60 • or 90 • (edge dislocation). Initial straight dislocations are placed randomly in a (15 μm) 3 cell subjected to periodic boundary conditions in all three directions. No artificial pinning points are introduced. After relaxation the initial dislocation density was ρ 0 ≈ 1.2 × 10 12 m −2 . The relaxed configuration, as shown in Fig. 1a, was then subjected to uniaxial tension at a constant strain rate ofε = 10 3 s −1 along 120 different crystallographic directions, sampled in the symmetry-irreducible stereographic triangle, as shown in the inset of Fig. 2a. The simulations were performed until the shear strain on the dominant slip system reached values in the range of γ d ≈ 1% to 5%, depending on the loading orientation and multiplication rate. An example of the predicted shear stress -shear strain curves corresponding to three different loading orientations are shown in Fig. 2a. From the shear stress -shear strain curves, the strain hardening rate ≡ dτ/dγ was extracted by fitting a straight line to the post-yield regime, as shown in Fig. 2a. The resulting strain hardening rates for 120 different loading orientations are shown in the inset of Fig. 2a.
In FCC crystals, plastic slip takes place on twelve 1 2 110 {111} slip systems as detailed in Table 1. In this work, we will use the Schmid and Boas (SB) notation to distinguish these slip systems. The slip systems are comprised of four different slip planes, denoted by A, B, C, and D, each containing three different Burgers vectors (slip directions). There are, however, only six unique Burgers vectors, denoted by b 1 , b 2 , ..., b 6 . In SB notation, each slip system is labeled with a letter-number pair (e.g., B4), corresponding to the plane and Burgers vector indices. The SB label for each slip system is given in Table 1. Within each slip plane, the three different slip systems are said to be coplanar with each other. The dislocation microstructure is a network formed by (physical) nodes and dislocation links 1 . Nodes are the locations where dislocations on different slip systems meet, and a link refers to a dislocation line connecting two nodes (Sills et al. 2018) (see Fig. 6). Each link is associated with a Burgers vector and a slip plane, hence, belongs to one of the twelve slip systems or is a non-glissile dislocation (e.g., Lomer or Hirth junctions). By analyzing the links and nodes of the network, we can gain insight into how it evolves. Figure 2 shows some of the data that can be obtained from the DDD simulations for the three loading orientations [001], [011] and [111], plotting the shear stress, total dislocation density ρ, and total number of dislocation links per unit volume N, as functions of the shear strain on the dominant slip system. The total densities ρ and N are the sums of the contributions from individual slip systems, i, denoted by ρ i and N i , respectively, plus contributions from sessile dislocations. The evolution of ρ i , N i , and other relevant microstructural parameters at the slip system level, including the accumulated plastic shear γ i , were recorded during all the simulations. Due to the statistical nature of DDD simulations, each of these variables fluctuates with time during the simulation, while a coarse-grained expression to model their evolution is expected to give smooth values as a function of time. Hence, some averaging is needed to reduce the fluctuations in the raw Fig. 2 a Shear stress-strain curves predicted by DDD simulations, along three high-symmetry loading orientations. Inset figure shows strain hardening rates from DDD simulations under constant strain ratė ε = 10 3 s −1 along 120 different loading directions. b Total dislocation density as a function of shear strain for the three loading orientations shown in (a). The windows of averaging, denoted by w r , r = 1, 2, . . . , 5 are illustrated on curve II and averaged values are shown by red dots. c Total number of link per unit volume for the three loading orientations shown in a Table 1 Slip systems in FCC crystals with their slip plane and Burgers vector b and their Schmid-Boas (SB) notation. The column G i lists glissile reactions between slip-system pairs that produce a dislocation on slip system i. DDD data. To this end, the raw data trajectory of each simulation is divided into 5 blocks as shown in Fig. 2b, and the time-averaged values for N i and ρ i are computed for each block.
To compute the plastic strain rateγ i for each slip system, we first fit the raw (γ i , t) data to a third order polynomial of time t. We then take an analytic derivative of the polynomial and evaluateγ i at the center time of each block. Unfortunately, the above procedure cannot be applied successfully to evaluate theṄ i andρ i , because of their large fluctuations during the simulation. Therefore, in "Correlation between microstructural features on different slip systems" section we use another measure to qualitatively compare dislocation multiplication from different simulations; we evaluate N i and ρ i at γ d = 0.5%, denoted byN i andρ i , respectively. Since all the simulations started from the same initial configuration,N i andρ i can be used as a proxy for the dislocation multiplication rate up to γ d = 0.5%.

Correlation between microstructural features on different slip systems
Having all the values of S i ,γ i , ρ i and N i for every slip system i from DDD simulations, we can now perform various analyses on the data, such as the identification of correlations between various features across different slip systems. Such an analysis helps to identify slip system interactions which are relevant to dislocation multiplication, especially slipfree multiplication which is the focus here. To this end, for each microstructural feature (e.g., ρ i ) on every slip system i, we form a column vector with the data from every simulation and every time (averaging) block and calculate the Pearson correlation between the vectors from different slip systems. The resulting correlation coefficient, r, for each pair of features is a scalar bounded between -1 and 1; the closer is the correlation coefficient to r = 1, the stronger is the correlation (both increase together), and values close to r = −1 indicate anti-correlation (one decreases as the other increases, and vice versa). r = 0 indicates no correlation. Figure 3a shows the Pearson correlation matrix between the plastic shear strain ratesγ i on slip systems i. The slip systems are indexed according to Table 1, so that the indices for  Table 1. All 120 simulations whose hardening rates are shown in the inset of Fig. 2a are used in the plots. Heavy black boxes denote sets of coplanar slip systems coplanar slip systems are consecutive. To emphasize this point, we use heavy black lines to denote groups of coplanar slip systems. Because all of our DDD simulations were performed with the tensile loading axes lying in a specific stereographic projection triangle as shown in the inset of Fig. 2a, the slip systems are not sampled symmetrically, leading to the observed asymmetries in the correlation matrices. As can be seen in Fig. 3, generally speaking the values ofγ i are not correlated with each other, with a few exceptions. For example, it is observed that values ofγ C1 andγ C3 are correlated with each other with the Pearson correlation coefficient of r (γ C1 ,γ C3 ) = 0.94. This correlation, however, is not a general feature of plasticity in FCC crystals, but rather specific to the stereographic triangle in which we have performed our DDD simulations. In this triangle, S C5 < 0.08 and S C1 = S C3 + S C5 2 for all loading orientations. Hence, the values of S C1 and S C3 are always very close to each other, causing the slip systems C1 and C3 to have similar plastic strain rates. Correlation (i.e., r > 0.7) is also observed between the following slip-system pairs: (γ A2 ,γ B2 ), (γ A2 ,γ D4 ) and (γ B2 ,γ D4 ). In all of these cases, corresponding pairs of S i values are also correlated with each other as shown in Fig. 3b. Comparing Fig. 3a and b, we conclude that the values ofγ i on different slip systems are generally independent from each other, except between slip systems whose Schmid factors are correlated. Contrary to the correlations among the shear strain ratesγ i , ρ i and N i clearly exhibit regions of high correlation where the Schmid factors are uncorrelated, as shown in Fig. 3c and d. Most notably ρ i and N i among coplanar systems are strongly correlated. This correlation indicates that dislocation multiplication may occur in a slip system with a low Schmid factor (i.e., slip-free multiplication) as long as one or more of its coplanar slip systems are active. These coplanar correlations are strongest on planes B, C, and D, and weakest on plane A. For the stereographic triangle used here, plane A is always the most favorably orientated plane with respect to the loading direction. We observed that in loading orientations where slip-free multiplication does not occur (see Fig. 5), always exactly one of the slip systems on plane A multiplies (e.g., the [ 1 2 3] loading orientation with A3 being the only slip system that is active and multiplies). This results in lower correlation on this plane. Finally, we note that the correlations observed for ρ i and N i are quite similar, indicating that the dislocation density and the link density evolve in a similar manner.
The fact that ρ i and N i values of coplanar systems are correlated, but shear strain rateṡ γ i are not, means that dislocation multiplication ratesρ i andṄ i cannot be simply proportional toγ i . For example, the commonly used Kocks-Mecking model (Mecking and Kocks 1981;) (see "Discussion" section) assumes thatρ i ∝γ i . Our data indicate that additional terms are necessary in the model to capture slip-free multiplication. Recently, Stricker and Weygand (2015) showed thatρ i depends not only on the shear rate of the slip system i, but also on the shear rate of other slip systems that can produce dislocations on slip system i through a glissile reaction (direct glissile multiplication). Alternatively, in our previous study (Akhondzadeh et al. 2020) we show that the evolution ofρ i in DDD can be described by an expression that depends not only onγ i but also onγ i andγ i , where i and i denote the two other slip systems coplanar with slip system i.
In order to obtain a better understanding of the interplay between coplanar slip systems and dislocation multiplication, we define S i max ≡ max (S i , S i , S i ) as the maximum Schmid factor of the three coplanar slip systems, and examine the correlation between ρ i and both S i and S i max . Examining correlations in this way will enable us to see the interaction between slip-driven multiplication (driven by slip due to large S i ) and slip-free multiplication (driven by S i max when S i ≈ 0). Figure 4 shows the values ofγ i , ρ i and N i aṫ γ d = 0.5% for all loading orientations, as functions of S i and S i max . Note that all the data points in these plots are from DDD simulations that share the same initial configuration with (approximately) the same dislocation density on each slip system. Accordingly, we can use the dislocation density of each slip systemρ i = ρ i (γ d = 0.5%) as a measure for the amount of multiplication which has occurred on that slip system (i.e.,ρ i is a proxy forρ i in this plot). Figure 4b and c clearly indicate two different regions where multiplication increases with either S i or S i max , designated by blue and red colors, respectively. In the blue region both dislocation multiplication and slip rateγ i increase with the Schmid factor S i . However, in the red region dislocation multiplication can be significant while the slip rateγ i is negligibly small. Hence, we refer to this type of multiplication as slipfree multiplication. Furthermore, in the red region dislocation multiplication increases 5% as a function of S i and S i max . For clarity, data points are also projected onto the bottom plane. Inset figures show the same data plotted as a function of S i − 0.5 S i max , which is equivalent to viewing the data along the boundary line between red and blue data points with increasing S i max , even when S i ≈ 0. The boundary between the two regions can be delineated by the line S i −0.5 S i max = 0, which according to the inset in Fig. 4b corresponds to the minimum values of dislocation multiplication. It can be shown that the red points in Fig. 4, for which S i < 0.5 S i max , always correspond to the slip system with the lowest Schmid factor among the three coplanar slip systems (see "Dependence of multiplication rate on Schmid factors" section for more discussions). The inset of Fig. 4b shows that the dislocation multiplication rate is appreciable only if |S i −0.5 S i max | > 0.1, which is denoted by a set of vertical lines in the insets; this may be considered as a rule-of-thumb for a necessary condition for dislocation multiplication. Note that this condition captures both slip-driven multiplication (blue data) and slip-free multiplication (red data), but is purely an empirical correlation revealed by our DDD database. We will discuss the meaning of this rule-of-thumb in "Discussion" section.
The data presented above clearly demonstrate that slip-free multiplication is a frequent phenomenon within our DDD database. The correlation analysis further indicates that coplanar slip system interactions are likely to play an important role in slip-free multiplication. Next, we quantify the fraction of the total dislocation density constituted by slip-free multiplication.

Occurrence and quantification of slip-free multiplication
In this section we analyze the loading orientations that exhibit slip-free multiplication, and show that this type of multiplication can constitute up to 25% of the total dislocation density. For this purpose, in the following we give a more precise criterion for the occurrence of slip-free multiplication. A slip system i is considered to undergo slip-free multiplication if the following conditions are satisfied.
where ρ i is the change of dislocation density on slip system i from the beginning to the end of the DDD simulation, and the subscript d indicates the dominant slip system. 3 There are 90 slip systems from 61 (out of 120 total) different loading orientations satisfying the three conditions in Eq. (1).
As stated in the previous section, in our simulations slip plane A always contains the slip system with highest resolved shear stress. We observed in our data (with only two exceptions) that if plane A does not show slip-free multiplication in a simulation, then the other planes do not exhibit slip-free multiplication either. Hence, in determining the loading orientations with slip-free multiplication, we only need to consider plane A. Figure 5a shows the minimum values of (S i − 0.5 S i max ) for i =A2, A3 and A6. According to the rule-of-thumb given above, this figure can be used to predict the loading orientations with slip-free multiplication. We show as dashed lines in Fig. 5a the boundaries separating loading orientations expected to exhibit slip-free multiplication from those that should not on the basis of this rule-of-thumb.
Next, in order to quantify the relative importance of slip-free multiplication, we denote byρ sf the sum of the dislocation density on those slip systems classified as undergoing slip-free multiplication according to Eq. (1), and plot its ratioρ sf /ρ to the total dislocation densityρ at γ d = 0.5% in Fig. 5b. This figure shows that the predictions of slip-free multiplication by the rule-of-thumb as shown in Fig. 5a are consistent with the DDD simulation data and a more precise criterion of slip-free multiplication. It can be seen that in simulations where slip-free multiplication occur, it accounts for 5-25% of the total density. For comparison, this ratio is approximately as high as the ratio of dislocation junctions (Lomer, Hirth, etc.) to the total density. The ratioρ sf /ρ was observed to gradually increase during the course of the simulation.
In the next sections, we discuss specific examples where slip-free multiplication is observed, why the direct glissile mechanism is insufficient to fully explain these observations, and propose two new slip-free multiplication mechanisms.

Analysis of different loading orientations
In order to understand the underlying mechanisms for slip-free multiplication, we first analyze the [ 0 0 1] loading orientation which hasρ sf /ρ = 0.24. In this case, there are 4 slip systems with zero Schmid factor and 8 slip systems with the maximum Schmid factor as shown in Table 1. Dislocations on the 8 active slip systems can collide and react with one another, e.g., to form glissile junctions. A schematic representation of the glissile  1)) as a fraction of the total dislocation densityρ at γ d = 0.5% interaction mechanism is shown in Fig. 6a. For each slip system i, there are two combinations of slip system pairs whose reaction produces a glissile junction which is a mobile dislocation on slip system i. These combinations are denoted by G i and they are listed in Table 1 for every slip system i. In each combination, one of the two reacting slip systems is coplanar with slip system i. For example, segments on the slip system A6 can be formed by reactions of the type A2+C3 as well as A3+B2. For the [0 0 1] loading orientation, dislocations on each of the slip-free slip systems can be produced via glissile reactions between active slip systems.
The situation is similar for the [ 1 1 1] loading orientation which hasρ sf /ρ = 0.24; slip systems A2, C5, and D4 have zero Schmid factor but still show dislocation multiplication. This is consistent with the rule-of-thumb in "Data analysis" section, because each of these three slip systems satisfies the condition that |S i − 0.5 S i max | = 0.13 > 0.1. Each of these slip systems can be produced via glissile reactions since all of the participating slip systems in G i are active (with the maximum Schmid factor). Note that while slip systems B2, B4, and B5 also have zero Schmid factor, their dislocation densities do not increase with strain. Considering B2 as an example, G i consists of B4+C5=B2 and B5+D4=B2. However, the four slip systems, B4, C5, B5, and D4, are inactive. Hence, we do not expect glissile junctions to form frequently on slip system B2. This is also consistent with the rule-ofthumb in "Data analysis" section, because the Schmid factors on all three slip systems on plane B are zero, so that |S i − 0.5 S i max | = 0 < 0.1. Hence, for [ 1 1 1] loading, the trends in the slip-free multiplication can be explained on the basis of the direct glissile mechanism.
On the other hand, the slip-free multiplication observed in [ 0 1 1] loading orientation cannot be explained by the direct glissile mechanism. Among the 8 slip systems with zero Schmid factors only A2 and B2 show dislocation multiplication withρ sf /ρ = 0.15. Taking A2 as an example, the possible glissile interactions which can result in a dislocation Fig. 6 a Direct glissile junction formation on slip system k from dislocations on slip systems i and j. Slip system i shares its slip plane with slip system k but not with slip system j, e.g., A2 + C3 = A6. Both end nodes for the new glissile link (k) are G-type. b Coplanar junction formation on slip system k from dislocations on slip systems i and j. All three slip systems, i, j, k are on the same slip plane, e.g., A2 + A3 = A6. Both end nodes for the new coplanar link (k) are P-type on A2 are A3+D6=A2 and A6+C3=A2. However, slip systems D6 and C3 do not multiply (their dislocation densities stay close to their initial values). In fact, the dislocation density on A2 greatly exceeds that on D6 or C3, even though all these slip systems have zero Schmid factors. Therefore, dislocation multiplication on A2 cannot be the (passive) reaction product between, say, A3 and D6, because the continuous production of glissile junctions requires additional dislocations in slip systems D6 or C3 (which do not exhibit multiplication). Even if a limited amount of dislocation content on A2 can be formed by the direct glissile reaction between A3 and D6, new A2 links cannot bow out to undergo slip-driven multiplication as the Schmid factor on A2 is zero, making the driving force for bowing zero. This suggests that there must be a mechanism other than direct glissile junction formation that causes dislocation multiplication on slip system A2. In order to uncover this mechanism, we continue the analysis for the [ 4 9 10] loading orientation, which is a simpler case than [ 0 1 1] because only one slip system exhibits slip-free multiplication withρ sf /ρ = 0.10; the complexity of network evolution during [ 0 1 1] loading made it difficult to identify mechanisms.
The plastic strain rate and dislocation density evolution on all slip systems for the [ 4 9 10] loading orientation are shown in Fig. 7a and b. It can be seen that slip systems A3 and A6 exhibit both high shear strain rate and high multiplication rate, because of their high Schmid factors (0.43 and 0.40, respectively). However, slip system A2 exhibits the behavior of slip-free multiplication; its dislocation multiplication rate (while lower than that of A3 and A6) far exceeds all the remaining slip systems, in spite of its very low Schmid factor (0.03). Similar to the [ 0 1 1] case, the copious dislocation multiplication on slip system A2 could not have been the result of direct glissile multiplication via glissile junction reactions A3+D6 or A6+C3, because there are far fewer dislocations on D6 or C3 than on A2.

Topology and complexity of the dislocation networks
To reveal the multiplication mechanisms on A2 for the [ 4 9 10] loading orientation, we examine the properties of all dislocation links on this slip system. Figure 7c shows the number and length distribution for dislocation links on A2 as a function of their average character angle (i.e., the angle between their end-to-end vector and their Burgers vector, see inset) for a network taken at γ d = 5% which is shown in Fig. 1b. A peak is observed at the 60°orientation, which is the character angle of newly formed glissile junctions. This is because the end nodes of glissile junctions are confined to remain on the intersection line between the two parent slip planes, thus enforcing the end-to-end vector of the junction to stay at 60 • . Hence a peak at the 60 • orientation is indicative of glissile junction formation events. However, we shall see that the situation is more complex. In fact, most of the A2 links at the 60 • orientation are not in the configuration of a direct glissile reaction.
To gain further insight into the nature of the links within the network, we analyze the nodes which connect the links. Table 2 shows the number of A2 links categorized based on the types of their two end nodes. Nodes are categorized into seven distinct types: G refers to a glissile node connecting the A2 dislocation with dislocations on A3 and D6, or with dislocations on A6 and C3. A link resulting from a glissile junction, such as one formed via the direct glissile mechanism, must have two glissile nodes of type G as shown in Fig. 6a. P refers to a coplanar node connecting the A2 dislocation with dislocations on A3 and A6. L refers to a collinear node connecting the A2 dislocation with an out-of-plane dislocation having the same Burgers vector. If the out-of-plane dislocation is on slip system B2, the node is categorized as binary collinear (L1). The collinear node is categorized as L2 if the out-of-plane dislocation is a junction dislocation (e.g., of the Lomer type). Hirth (H) and Lomer (M) nodes are three armed nodes, where exactly two arms are glissile and their interaction produces Hirth and Lomer junctions, respectively. For example, interactions of A2 with C1 and C5 produces Hirth and Lomer interactions, respectively. All other nodes are categorized as type O for "other, " indicating that they are not the result of a binary reaction. Examples of O-type are nodes connected to three glissile, non-coplanar arms (e.g., A2, C5, D4) or nodes connected to four or more arms. We will introduce one possible mechanism for the formation of O-type nodes in "Assisted glissile mechanism" section. In Table 2 we show the node configurations for all links, and also separate out links with a character angle of 60 ± 2.5 • from other character angles. For brevity, we grouped Hirth and Lomer-type nodes together, as well as L1 and L2-types.
A number of observations can be made from Table 2. First, a freshly created glissile dislocation on A2 must have both nodes of the G type and be oriented at 60 • . There are only Table 2 Number of dislocation links on slip system A2 with different end node types at γ d = 5% for the [ 4 9 10] loading orientation. Columns and rows denote the two end nodes of each link. There are 1488 dislocation links on slip system A2. G, P, L=L1+L2, H, and M refer to glissile, coplanar, collinear, Hirth, and Lomer nodes, respectively, as shown in Figs. 6 and 9a. All other nodes are categorized as type O. 259 A2 links that are connected to a glissile node of the type A2+B5=B4 and A2+B4=B5 have been excluded, because they are considered as A2 links participating in a glissile reaction, instead of being the product of a dislocation reaction 6 such dislocation links, i.e., a very small fraction out of 1488 total links on A2. Second, the number of A2 links with both ends being coplanar nodes (P) is 32, which is more than four times the number of newly created glissile dislocations. Third, a large number of A2 links (101) have one node of coplanar type (P) and the other node of collinear type (L); we will show one formation mechanism for such links in the next section. Fourth, most of the A2 links have one or both nodes of the type O, meaning that they are not the results of any binary reactions between ordinary dislocations on the 12 slip systems; they are the results of (secondary or tertiary) reactions between products of previous dislocation reactions.
To further investigate this point, we analyze the node types distribution on all slip systems at γ d = 0.5% for all different loading orientations. Results for the contribution of each category to the total number of nodes in the microstrucutures are shown in Fig. 8. Note that unlike Table 2, this figure does not contain link information; this figure only characterizes the nodes, not their connectivity. Here we find that, as a general feature among all dislocation microstructures for different loading orientations, O-type nodes constitute the largest fraction of the total nodes among all categories. Combining O and L2 nodes, we find that higher-order nodes (i.e., nodes that cannot result from a binary collision) amount to approximately 50% of the nodes in the networks, indicating the complexity of dislocation microstructure. This complexity cannot be explained solely by the binary reaction between slip system pairs nor by second-order junctions defined as glissile junctions interacting with forest dislocations (Madec and Kubin 2008), because these mechanisms cannot produce O-type nodes. Thus, the abundance of O-type nodes as shown in Fig. 8 suggests that other mechanisms are at play. One example of such mechanisms is the assisted glissile mechanism, which was frequently observed in DDD simulations, and will be presented in "Assisted glissile mechanism" section. This mechanism involves a series of reactions between nodes within the microstructure, as opposed to collisions between links. For these reasons, we conclude that dislocation networks are more complex than assumed in previous theories.

Glissile+coplanar multiplication mechanism
We propose that the majority of A2 links for the [ 4 9 10] loading orientation are created by reactions between coplanar slip systems, i.e., A3 + A6 = A2, as opposed to glissile Fig. 8 Distribution of dislocation node types from network configurations at γ d = 0.5%. In each microstructure the total number of nodes is in the range of 3500-9000. See text for node type definitions reactions A3 + D6 = A2 and A6 + C3 = A2. Such a proposition may appear unfounded at first, because when A2 links are produced via direct collisions between A3 and A6 links the resulting links must have both end nodes of type P as shown in Fig. 6b, and such P-P links are only a small fraction of the data in Table 2. However, as we shall see in details below, we have identified another two-step mechanism, which we refer to as the glissile+coplanar mechanism, which we believe is key and provides the necessary process to explain the variety of link types we observe. The glissile+coplanar mechanism produces an A2 link connected to a coplanar node (P) and a collinear (L) node. Because P-type nodes are quite mobile in the slip plane, they are able to encounter forest dislocations intersecting the slip plane, causing subsequent reactions that converts the nodes into other types. This mechanism can for instance lead to the creation of glissile (G) type node, thus providing an explanation for the origin of the A2 links of the type G-L, that are not at 60 • orientation (see Table 2).
The glissile+coplanar mechanism is illustrated in Fig. 9a and operates as follows. Consider a dislocation on slip system A3 which first encounters a forest dislocation with Burgers vector b 2 (which could be a dislocation on slip system B2 or a Lomer junction with the same Burgers vector). The reaction between these two dislocations will produce a dislocation on slip system A6, which is on the same plane as the A3 dislocation. (If the forest dislocation is on slip system B2, then this reaction is the glissile reaction A3 + B2 = A6.) Because slip system A6 has a large Schmid factor under [ 4 9 10] loading, the A6 dislocation readily bows out and moves on the slip plane. The A6 dislocation can then react with the same A3 dislocation which also bows out on the same plane (A3 + A6 = A2). The result is an A2 link with a P-type node and an L-type node, as illustrated in Fig. 9a. An example from DDD is shown in Fig. 9b. Examining Table 2, we see 103 links with P and L nodes, suggesting this mechanism is important for slip-free multiplication. Alternatively, A6 may encounter another A3 dislocation on the same plane, and form an A2 link with both nodes of P-type, an example of which is shown in Fig. 9c. There are 31 A2 links of this type in Table 2. Additional support for the role of glissile+coplanar mechanism for the slip-free multiplication is provided in Appendix B through performing specialized simulations.
As a result of this process, we also believe that the A2 links that are not 60°o riented are more likely the result of reactions involving A2 links created by the coplanar mechanism (A3 + A6 = A2) than those from the direct glissile mechanism (A3 + D6 = A2 or A6 + C3 = A2). This is because firstly, a coplanar (P) node is actively dragged by the A3 and A6 links, hence it is free to move in any direction in the slip plane, greatly enhancing its chance to collide with a forest dislocation leading to change of both node type and orientation. Secondly, it is unlikely that an A2 link formed by the glissile mechanism changes its character angle through interaction with forest dislocations, because dislocations on A2 (as well as all the dislocations on planes B, C and D) have very low Schmid factor, and are thus unlikely to bow out and react. We note that the average length of A2 links with an orientation of 60°is about 75% of the average length of A2 links with other orientations. This suggests that the A2 links formed by the direct glissile mechanism remain short compared with those with coplanar nodes that move around on the slip plane without any line constraint. Fig. 9 a Schematic of an A2 segment formed by a glissile+coplanar mechanism, following the glissile junction formation shown in Fig. 6a, with i=A3, j=B2, k=A6. Both A3 and A6 are active and can bow out to form a new A2 link. G, P and L refer to the type of link end nodes (see text). b Two consecutive DDD snapshots showing the glissile+coplanar mechanism. Only dislocations on one plane of type A and dislocations on slip system B2 are shown. c Two consecutive DDD snapshots showing the coplanar (A3+A6=A2) mechanism. Only the dislocations on a single plane of type A are shown. In b and c, the scale bar is in units of b = 0.255 nm, and the coloring scheme is consistent with a

Assisted glissile mechanism
In the previous sections we presented the direct glissile and glissile+coplanar mechanisms as two mechanisms that can lead to slip-free multiplication, with the latter able to operate even in loading orientations where the former is not possible. In this section, we discuss yet another slip-free multiplication mechanism which we have frequently observed in our DDD simulations, and show how it creates a node connecting three glissile dislocation links. As for the preceding sections, this mechanism provides additional support for the complex evolution of the dislocation network during plastic deformation.
As an example, let us consider the [ 1 1 1] loading orientation, in which slip systems A2, C5, D4 have zero Schmid factor and experience slip-free multiplication. Dislocations on each of the three slip systems can be produced by direct glissile reactions between dislocations on active slip systems, for example: A3+D6=A2, A3+C1=C5, C1+D6=D4. Note that slip systems i=A3, j=C1, and k=D6 all have the maximum Schmid factor. Figure 10a shows the slip planes and Burgers vectors labeled using SB notation. Consider the situation shown in Fig. 10b where a dislocation on slip system i has reacted with a dislocation on slip system j (A3+C1=C5), and also with a dislocation on slip system k (A3+D6=A2), but the two dislocations on j and k have not yet collided. In our DDD simulations, we find that in such a situation the glissile junctions A2 and C5 often extend along their line direction until they meet at the apex of the tetrahedron, creating a 4-arm node connecting dislocations A2, C5, D6, and C1, as shown in Fig. 10c. The new 4-arm node is unstable and quickly splits into two nodes that are connected to each other with a link on the D4 slip system, as shown in Fig. 10d. Since the dislocation link on D4, which could alternatively result from a direct glissile reaction between slip systems j and k, is instead formed through the interaction with slip system i, we call this mechanism the assisted Fig. 10 a The slip planes (A, C, D) and Burgers vectors (b 1 , · · · , b 6 ) in the SB notation. Note that the A, B, C, D labels in the SB notation do not match those in the Thompson tetrahedron notation. b A dislocation on slip system i=D6 involved in two glissile reactions with slip systems j=A3 and k=C1: D6+A3=A2 and D6+C1=D4. G refers to the type of the nodes. c The two glissile junctions A2 and C5 meet at the apex of the tetrahedron, forming a 4-arm node. This node is O-type. d The 4-arm node splits into two 3-arm nodes connected by a dislocation link on slip system D4. The node on the apex of the tetrahedron is connected to three 60 • dislocations on slip systems A2, D4 and C5 glissile mechanism. An interesting consequence of the assisted glissile mechanism is that the node at the apex of the tetrahedron is connected to three 60 • dislocations on three slip systems (A2, D4, C5) all having zero Schmid factor. Furthermore, all three dislocation links (A2, D4, C5) have one end node of G-type and the other node of O-type.
The assisted glissile mechanism can also operate for other combinations of (i, j, k) slip systems as long as they belong to three different slip planes. For example, in the case of [ 4 9 10] loading orientation, the combination of i =C1, j =D6 and k =A3 would result in an A2 link that is connected to A3 and D6 links on one end, while being connected to C5 and D4 links on the other. This A2 link has a character angle of 60 • and is characterized as G-O type. Hence, the assisted glissile mechanism is an example of a series of reactions by which a glissile node can transform into a more complex configuration (e.g., those in Table 2 categorized as G-O links). Note that the assisted glissile mechanism still enforces the 60 • orientation of the glissile links. This supports the picture of the dislocation network primarily consisting of nodal structures that are not the result of binary reactions between oridinary dislocations on the 12 slip systems.

Dependence of multiplication rate on Schmid factors
In "Correlation between microstructural features on different slip systems" section we showed that a necessary condition for appreciable dislocation multiplication on slip system i to occur is |S i −0.5 S i max | > 0.1. Here we interpret the meaning of this rule-of-thumb, which was obtained empirically based on the data in Fig. 4. For simplicity, we limit the discussions in this section to the slip plane which contains the dominant slip system. For other slip planes, discussions are still valid as long as S i max of that plane is sufficiently high to cause dislocations multiplication. Recall that S i max is the maximum Schmid factor among the three slip systems coplanar with i, which we denote as i and i . According to the rule-of-thumb, multiplication is greatest on slip system i if either S i ≈ S i max or S i ≈ 0. Furthermore, we find that multiplication is minimized when S i ≈ 0.5 S i max . To help interpret the meaning of this expression, consider the schematic shown in Fig. 11, which shows the glide plane with the relative orientations of the Burgers vectors for slip systems Fig. 11 Schematics of relative orientation of applied shear stress, τ = t − (t · n) n, where n and t are unit vectors denoting the normal to plane containing slip systems i, i and i , and the loading orientation, respectively. a S i = S i max corresponding to rightmost boundary of blue data points in Fig. 4b; in this case system i experiences slip-driven (i.e. regular) multiplication. b S i − 0.5S i max = 0 corresponding to the boundary between blue and red data points in Fig. 4b; in this case multiplication on system i is minimum. c S i = 0 and S i = S i = S i max corresponding to leftmost boundary of red data points in Fig. 4b; in this case slip-free multiplication on system i can occur i, i , and i , and the applied shear stress in three cases: (a) S i = S i max , (b) S i = 0.5 S i max , and (c) S i = 0. These three cases correspond to (a) the rightmost boundary line of blue data points, (b) the boundary line between red and blue data points and (c) the leftmost boundary line of red data points, in Fig. 4b. Clearly when S i ≈ S i max , slip system i is most active so that it will experience slip-driven multiplication. On the other hand, when S i ≈ 0.5 S i max , Fig. 11b shows that either i or i has the largest Schmid factor and will experience slip-driven multiplication, but neither i nor the other inactive slip system is expected to experience multiplication according to Fig. 4b. Finally, in the S i ≈ 0 case which corresponds to a large fraction of slip-free multiplication as shown by Fig. 4b,  Fig. 11c shows that S i max ≈ S i ≈ S i . Accordingly, we expect that slip systems i and i will be equally active. As a result, coplanar collisions between i and i are likely, which results in the generation of dislocations on slip system i. Based on this reasoning, we see that a large fraction of slip-free multiplication may simply be the result of coplanar collisions between equally driven slip systems.
The scenario shown in Fig. 11c occurs in high symmetric loading orientations, i.e, orientations close to the three corners of the stereographic triangle, as well as orientations that are close to the [ 0 1 1] −[ 1 1 1] edge (e.g., [ 4 9 10]), as shown in Fig. 5a. Slip systems in the other loading orientations, e.g. those closer to the center of the stereographic triangle, exhibit the scenario of either Fig. 11a or b. Hence, the distinction between multiplication types based on Schmid factors as shown in Fig. 11 is consistent with the fraction of dislocation density on slip systems showing slip-free multiplication as shown in Fig. 5.

Modification to Kocks-Mecking model of dislocation multiplication
The Kocks-Mecking model (Mecking and Kocks 1981) for dislocation multiplication expresses the rate of change of the total dislocation density ρ with shear strain γ in the following form: where c 1 and c 2 are dimensionless constants and b is the magnitude of the Burgers vector.
It was subsequently generalized to describe dislocation multiplication on individual slip systems ). In our recent study (Akhondzadeh et al. 2020), we found that a correction term needs to be added to the generalized Kocks-Mecking relation, leading to the following expression for the dislocation multiplication rate on slip system i: where the first term is the generalized Kocks-Mecking expression (Devincre et al. 2008) and the second term is the correction to account for slip-free multiplication (i.e., contribution toρ i even whenγ i = 0). In Eq.
(3) a ij are interaction coefficients between slip systems i and j, andc 1 ,c 2 andc 3 are dimensionless parameters that depend on the loading orientation.
While the correction term in Eq.
(3) was motivated by the need to describe the dislocation density data from a large set of DDD simulations (Akhondzadeh et al. 2020), the findings in this work provide a mechanistic explanation for the coplanar form of the correction term (involving i and i ). If, on the other hand, one assumes that the slip-free multiplication is mainly the result of the direct glissile mechanism, then the modified Kocks-Mecking relation would take the following form (Stricker and Weygand 2015), where m and k are the slip systems that can react to form a glissile junction on slip system i, according to Table 1. In Eq. (4), the repeated index j = 1, ..., 12 is summed over, but the indices i, m and k are not. However, we can see that for loading orientations such as [ 0 1 1] and [ 4 9 10], Eq. (4) is not consistent with the DDD data, because there is no slip system pair that can form glissile junctions on slip system i (e.g., A2) with both systems active in the pair (i.e., with appreciableγ m ,γ k , and ρ m , ρ k ). In comparison, Eq. (3) is able to account for slip-free multiplication in these orientations. Furthermore, Eq.
(3) implies that slip-free multiplication occurs only when both of the other two slip systems on the same slip plane are active. This condition for which slip system can exhibit slip-free multiplication is indeed confirmed in our DDD simulations along all except one of the loading orientations. 4 As a result, when slip-free multiplication occurs in one slip system, then it is always observed that multiplication occurs in all three coplanar slip systems. This discussion demonstrates the value provided by mechanistic insight available through DDD when developing coarse-grained theories for dislocation dynamics. In "Glissile+coplanar multiplication mechanism" section we showed how the glis-sile+coplanar mechanism brings dislocations on slip systems i and i to the same plane to react and form dislocations on slip system i. For i = A2 in the [ 4 9 10] (as well as [ 0 1 1]) orientation, this requires the presence of forest dislocations with Burgers vector b 2 , such as dislocations on slip system B2. However, dislocations on B2 are not consumed in the glissile+coplanar reactions described in "Glissile+coplanar multiplication mechanism" section. Instead, they act as catalysts to bring dislocations on i =A3 and i =A6 slip systems to the same plane to react. Therefore, in Eq. (3) the correction term only involves the plastic shear strain rate on the coplanar slip systems,γ i andγ i , although the coefficientc 3 may depend on the dislocation density on slip system B2. Additional research is necessary to continue to optimize the functional form of these terms and determine strategies for obtaining values for the relevant parameters.

Importance and robustness of slip-free multiplication
In this work, we have attempted to elucidate the importance, pervasiveness, and origins of slip-free multiplication. We find that slip-free multiplication constitutes a significant fraction of the total dislocation density under multi-slip conditions (e.g., loading axes near the corners of the stererographic triangle, see Fig. 5b). Given that both the flow stress and the dislocation line tension are influenced by the total dislocation density (Anderson et al. 2017;Cai and Nix 2016), slip-free multiplication is expected to play a major role in metal plasticity. Furthermore, we expect that slip-free multiplication is operative in polycrystaline solids, since the multiaxial stress state within each grain typically activates multiple slip systems (Argon 2008;Kocks 1970). With that said, we note that all of our results here correspond to uniaxial loading; additional research is necessary to understand slip-free multiplication in more complex loading scenarios. For example, slip-free multiplication could be important in scenarios where the loading axis is changed (e.g., latent hardening experiments), so that a previously inactive and slip-free slip system becomes activated.
In our analysis above, we concluded that coplanar reactions are an important contributor to slip-free multiplication. At first glance, we may therefore conclude that slip-free multiplication via coplanar reactions would be rather limited, since such reactions require dislocations to be on the same (or very close) atomic planes. However, the glissile+coplanar and assisted glissile mechanisms demonstrate how other interactions within the network can "bring" dislocations together on the same glide plane, thereby enabling coplanar reactions. These mechanisms provide examples where it is incorrect to view slip system interactions simply on the basis of slip system pairs; Both the glis-sile+coplanar mechanism and the assisted glissile mechanisms involve three or more slip systems.
Finally, our analysis uncovers the role of the coplanar interaction as a key mediator for dislocation multiplication, thereby complementing the conventional view in which coplanar reactions are limited to dipolar, elastic interactions between dislocations on separate, parallel slip planes Kubin et al. 2008a;2008b).

Conclusion
We analyzed the DDD simulation results of 120 loading orientations in the absence of cross-slip, to investigate the relationship between dislocation multiplication and plastic shear rate,γ i . We observed that while the plastic strain ratesγ i on slip systems are generally not correlated with one another, dislocation densities ρ i on coplanar slip systems are correlated. Specifically, some slip systems multiply under zero shear strain rate, a phenomenon we refer to as slip-free multiplication. We observed that all the slip systems i which exhibit dislocation multiplication satisfy the rule-of-thumb |S i − 0.5S i max | > 0.1.
By analyzing the case of [ 4 9 10] loading orientation in detail, we demonstrated that the direct glissile mechanism is unable to explain slip-free multiplication for certain loading orientations. We then proposed the glissile+coplanar mechanism, which could explain the observed slip-free multiplication. The mechanisms justify a correction term to the Kocks-Mecking model for dislocation multiplication that depends on the plastic strain rate on coplanar slip systems. By analyzing the node connectivity of links on slip systems showing slip-free multiplication, we observed that most of these links are not the result of a simple reaction between ordinary dislocations on the 12 slip systems. Instead, they are the results of a series of dislocation reactions, such as the assisted glissile mechanism, leading to a dislocation network with a complex connectivity. All together, our results motivate the need for a deeper analysis of dislocation network topology which can act as a guide for physics-based constitutive theory development. Table 3 lists all the input simulation parameters used in creating the DDD database as explained in Akhondzadeh et al. (2020).

Appendix B Confirming the necessity of coplanar reaction in the [4 9 10] orientation
According to the coplanar mechanism proposed in "Glissile+coplanar multiplication mechanism" section for the [ 4 9 10] orientation, the multiplication of dislocations on slip system A2 requires an out-of-plane (forest) dislocation with Burgers vector b 2 , which brings two dislocations on A3 and A6 slip systems to the same plane in order to react and form A2. In the initial configuration for a DDD simulation, such forest dislocations should mostly consist of dislocations on slip system B2. On the other hand, if the direct glissile mechanism (A3 + D6 = A2 and A6 + C3 = A2) is responsible for the slip-free multiplication, we would expect that the presence of slip system B2 is not necessary and A2 would still multiply if only slip systems A3, A6, C3 and D6 are present.
In order to test the role of the coplanar mechanism for slip-free multiplication, we performed a DDD simulation along the [ 4 9 10] loading orientation starting from an unrelaxed initial configuration which only contains dislocations on the A2, A3, A6, C3 and D6 slip systems. To enable a fair comparison, we also repeated the same simulation along [ 4 9 10] as shown in Fig. 7 where all 12 slip systems are present, but starting from an unrelaxed initial configuration. Figure 12 shows the evolution of dislocation density ρ i on A2 slip system in these simulations. It can be seen that dislocations on A2 do not multiply, if the initial configuration consists only of those slip systems (C3 and D6) Simulation cell size 15μm Fig. 12 Evolution of ρ i for i =A2 from four DDD simulation along the [ 4 9 10] loading orientation with an unrelaxed initial configuration that only contains a subset of slip systems as shown in the legend involved in the direct glissile mechanism. However, by adding dislocations on the slip system B2 to this initial configuration, A2 multiplies at a rate of approximately 20% of the normal simulation. We confirmed that adding a different sixth slip system other than B2, would result in minimal multiplication of A2 slip system. We take this behavior as an indication of the essential role that the glissile+coplanar mechanism plays in the slipfree multiplication. It is interesting to note that by removing slip systems C3 and D6 from the initial configuration (that is, removing the possibility of A2 formation via glissile junction formation), A2 slip system shows multiplication of more than 60% of a normal simulation. One possible explanation for such behavior is that inclusion of more slip systems, increases the density of dislocation with b 2 Burgers vector (other than dislocations on B2 system), hence, increasing the rate of glissile+coplanar mechanism. For example, reaction of slip systems C5 and D4 result in a Lomer junction with Burgers vector b 2 . Another explanation could be that there are other mechanisms at play, responsible for slip-free multiplication on A2 which involves a series of reactions between different slip systems.