Junction formation rates, residence times, and the rate of plastic flow in FCC metals

During plastic flow in metals, dislocations from slip systems with different glide planes collide to form junctions. After being in-residence within the dislocation network for some period of time, these junctions then break, thereby liberating the attached dislocation lines. In this work we use random forest discrete dislocation dynamics simulations to quantify the junction formation rate and junction residence time as a function of stress for all junction types in face-centered cubic metals. We then relate these quantities to the dislocation link-length distribution, which is found to exhibit an exponential form. This enables us to quantify the mean junction strength and also the slip system interaction coefficients. Finally, using the link-length model we obtain a flow rule for our systems which is physics-based with all parameters determined from DDD simulations. The insights here provide a path forward for a dislocation network theory of plastic flow based on the link-length distribution.


Introduction
During plastic deformation of metals, activation of different slip systems leads to evolution of the dislocation network (Devincre et al. 2008;Kubin et al. 2008;Sills et al. 2018).An important aspect of dislocation network evolution is the formation of dislocation junctions.Junctions form when two dislocations from different slip systems intersect and react, as shown in Fig. 1.If the two slip systems have different glide planes, we often refer to the slip system with higher resolved shear stress as the primary slip system and the other as the forest slip system.Junctions are important because they inhibit the motion of dislocations and increase the flow stress of the dislocation network (Madec et al. 2002;Kubin et al. 2003;Sills et al. 2018).Numerous aspects of dislocation junctions have been considered in the past, including their likelihood of formation (Madec et al. 2003;Kubin et al. 2003;Capolungo 2011), their strength (Madec et al. 2003;Devincre et al. 2008;Devincre and Kubin 2010;Capolungo 2011;Bertin et al. 2014;Wu et al. 2013Wu et al. , 2016)), and their contribution to the flow stress and strain hardening (Bulatov et al. 2006;Devincre et al. 2008;Kubin et al. 2008;Bertin et al. 2014;Sills et al. 2018).However, the role of junction formation in dislocation network evolution is not as well understood.
Page 2 of 21 Zhang and Sills Journal of Materials Science: Materials Theory (2024) 8:11 One area that has been extensively studied is the influnce of junctions on the quasi-static strength of crystals.Perhaps the most popular approach focuses on the so-called generalized Taylor law, which says that the quasi-static flow stress of slip system i is (Devincre et al. 2006;Devincre et al. 2008;Kubin 2013) where µ is the shear modulus, b is the magnitude of the Burgers vector, a ij is the interac- tion coefficient for slip systems i and j, and ρ j is the dislocation density of the forest slip system j.In Eq. ( 1), all quantities are readily determined except a ij , which is dependent on the strength of and mean spacing between the dislocation junctions that form during plastic flow.A commonly utilized approach for determining a ij is the random forest model (Madec et al. 2002;Devincre et al. 2006).In this approach, discrete dislocation dynamics (DDD) simulations are performed using dislocations from just two slip systems, the primary and forest systems.A glide plane in the primary slip system is randomly populated with forest dislocation line segments giving forest density ρ f , defined as the number of forest intersections with the glide plane per unit area (see Appendix), and one or more lines from the primary system.An applied strain rate then loads the primary slip system, forcing the primary dislocations to glide through the forest.The critical stress at which plastic flow ensues can then be used to determine a ij through use of Eq. ( 1).Computation of a ij in this way has been performed for face-centered cubic (Devincre et al. 2006;Madec and Kubin 2017), body-centered cubic (Queyreau et al. 2009;Madec and Kubin 2017), hexagonal close-packed (Bertin et al. 2014), fluorite (Madec et al. 2023), and ice (Devincre 2013) crystals.
Of course, while powerful, the generalized Taylor law does not provide a complete picture of plasticity and the influence of dislocation junctions.For one, it is only valid under quasi-static loading conditions, where dislocation drag forces do not influence plastic flow (Akhondzadeh et al. 2020).In many settings, high strain rate loading is relevant, necessitating constitutive laws that go beyond the quasi-static limit.Secondly, Eq. ( 1) only dictates when plastic flow commences (e.g., when τ > τ i c ), but provides no information on what the rate of plastic flow is.In plasticity modeling, a flow rule must also be specified to prescribe the plastic shear strain rate γ i as a func- tion of the resolved shear stress τ i .For example, a power-law flow rule is commonly employed (Roters et al. 2010): (1) where γ0 and m are material parameters.Equation ( 2) is simply a phenomenological equation that captures the increase of plastic strain rate with stress and has no micromechanicistic basis.While other flow rules have a sounder physical basis, we argue that all existing flow rules lack a rigorous connection to the processes underlying evolution of the dislocation network, such as dislocation bowing and formation/breaking of dislocation junctions.Through our work here, we take a different approach to plasticity modeling which goes beyond the Taylor law and phenomenological flow rules.
The fundamental goal for this work is to quantify the influence of dislocation junctions on dislocation network evolution in FCC crystals, and thereby enable a unified theory for plasticity which does not require phenomenology and is valid beyond the quasi-static limit.By analyzing the physical operations associated with forest interactions in random forest DDD simulations (e.g., junction formation and breaking), the stress-dependence of junction formation and the junction residence time are found.We then show how the junction residence time, the quasi-static flow stress, and the average junction strength are interconnected.Finally, we integrate all of this information to construct a physicsbased flow rule with no arbitrary fitting constants which matches the DDD results.

Methods
Here we use random forest simulations to study plastic flow in FCC crystals.Random forest simulations have the advantage that they are topologically simple (only binary junctions can form) and, as we will show below, achieve a steady-state flow condition.These features make them relatively easy to analyze for our purposes here, in comparison to more traditional strain hardening simulations (Bulatov et al. 2006;Bulatov and Cai 2006;Sills et al. 2018).However, random forest simulations have the disadvantage that they feature unphysical dislocation geometries, such as a forest comprised of a single slip system and forest dislocations which terminate inside the crystal (e.g., pinned end nodes).Our results below may be affected by these simplifications and should be further assessed in the future with other simulation techniques and experiments.
All DDD simulations were performed using ParaDiS with the GPU implementation of the subcycling time integration scheme (Arsenlis et al. 2007;Sills et al. 2016;Bertin et al. 2019).Simulation parameters are given in Table 1.Cross-slip was not allowed in our simulations.Forest models of all four types of junctions in FCC crystals (collinear, glissile, Hirth, and Lomer) were constructed with a constant forest density ρ f = 1.53 × 10 12 m −2 using forest dislocation line segments of random orientations and spatial positions.Each forest segment had a length of 0.7/ √ ρ f .Since the forest segment orientations are ran- dom, the forest density is related to the total density in the forest system by ρ f = f ρ , where f = 2 π (see Appendix).The three initial mobile dislocation lines were evenly separated in the primary slip plane and periodic boundary conditions were imposed in all directions in the 24 micron cubic simulation cell, resulting in a primary density of ρ p = 7.49 × 10 9 m −2 .An example simulation snapshot is shown in Fig. 2. Based on con- straints of the GPU ParaDiS implementation, the crystallographic coordinate system was used for the simulations, resulting in a (1 1 1) glide plane that is oblique to the periodic (2) , simulation cell.With this chosen configuration, the mobile dislocations never contacted each other during our simulations and also never moved fully through the periodic simulation cell.Note that our choice of forest and mobile dislocation densities is arbitrary.
In the Results section, we have non-dimensionalized our results in the hope that they may apply to other dislocation densities, but future work is necessary to confirm this.Our box size was chosen (with fixed forest density) based on convergence testing (e.g., results are insensitive to box size).Furthermore, our choice of forest dislocation length is arbitrary and may influence our results; we have not specifically determined its influence on our results.Specific stress tensors were applied to each junction case such that the resolved shear stresses on the forest slip system and glissile junctions were zero, the term resolved shear stress in the following text only refers to the resolved shear stress on the primary slip system.Information of the slip systems and stress tensors for each junction case is given in Table 2.This approach differs from other works using the random forest model where constant strain rate loading was used.Using constant stress loading instead allows us to systematically study the stress-dependence of plastic flow.
During each simulation, the details of all topological operations were saved and translated into physical operations including collisions between forest and primary dislocations, formation of junctions, and breaking of junctions.These physical operations were then subsequently analyzed in order to compute the forest collision rates, junction formation rates, and junction breaking rates.Additionally, the simulation time of each junction formation and breaking event was used to compute the residence time for each junction, defined as the difference between the breaking and formation time for a junction.
Two techniques were used to analyze plastic strain rates in our simulations.The total shear plastic strain rate from the motion of all dislocation segments was obtained via linear regression of the plastic strain vs. time curve.To obtain detailed information on the contributions to the shear strain rate from individual dislocation segments, the positions for each dislocation node at every time step were saved.These were then subsequently analyzed to compute the area sweep rate Ȧm for each dislocation segment m.By definition, the total shear plastic strain rate is γ = b m Ȧm 2V , where V is the volume of the simulation cell.

Results
During each simulation, we observed an initial transient phase prior to reaching a steady-state plastic flow condition.During steady-state plastic flow, the plastic strain rate was approximately constant in time.Figure 3 shows a few examples of the plastic strain vs. time data for a few levels of stress.The brief initial transient period is visible, and during the subsequent steady-state plastic flow the plastic strain increases linearly with simulation time.The steady-state plastic strain rate is obtained by linear regression based on the data points of steady-state plastic flow.For all analyses below, we only consider data obtained during steady-state plastic flow.We present results in both dimensional and non-dimensional form.Specifically, we normalize residence times by the characteristic time scale B/µ , where B is the dislocation drag coefficient, and resolved shear stresses by the characteristic Taylor stress τ * ≡ µb √ ρ f .

Junction formation
First we analyze the rate of junction formation.As previous research has shown (Madec et al. 2003;Kubin et al. 2003), junctions only form under a subset of line intersection geometries.However, prior studies only studied the likelihood of junction formation under no applied stress.Here we seek to understand what fraction of forest collisions κ result in junction formation as a function of the applied stress level.
To determine κ , we analyze the forest collision, junction formation, and junction breaking events from our simulations.Figure 4 shows the number of each event type as a function of time for four different stress levels.Similar to the plastic strain rate data, we see an initial non-linear transient period, after which the rates of forest collision R C , junction formation R F , and junction breaking R B become constant.We extract these rates via linear regression fits in the steady-state region.Note that in the steadysteady regime, R F ≈ R B so that the number of junctions in the system at any point in time is approximately constant.The gap between the formation (green) and breaking (blue) curves in Fig. 4 indicates the number of junctions existing in the system during steady-state.As the stress increases, this gap shrinks, indicating that there are fewer and fewer junctions in the system, as expected.Note that during steady-state flow, none of the junctions are permanently stored, eventually every junction breaks.
By definition, the junction formation fraction is simply the ratio of the formation and collision rates, Using this definition, we plot κ as a function of the resolved shear stress for all four junction types in Fig. 5.In all cases, κ monotonically decreases with stress in a sig- moidal fashion.This indicates that as stress increases, junction formation becomes where 2) with erfc as the complementary error function.The param- eters τ µ and τ σ are the mean and standard deviation, respectively, for the normal distri- bution underlying the Q-function.κ ∞ is the asymptotic value of κ when τ → ∞ .Values for these parameters for each junction type are given in Table 3.The dislocation operation with the highest energy dissipation rate is preferred (Bulatov and Cai 2006).Junction formation occurs if the energy dissipation rate of continuous plastic flow is lower than the energy dissipation rate of junction formation.As we increase the resolved shear stress, the energy dissipation rate of continuous plastic flow increases and fewer and fewer junctions are preferred to form.In other words the stress-dependence of junction formation fraction of Eq. ( 4) reveals the distribution of critical stress of junction formation.We demonstrate below that the parameters τ µ correlates with the mean junction strength, indicating that the distribution function Eq. ( 4) derives from the statistical distribution of junction strengths in the dislocation network.

Junction residence times
Next we turn our attention to the junction residence times.The residence time for each junction is the total simulation time that a junction is present in the simulation cell (e.g., the time over which it is "in-residence").Collecting together junction residence times at different stress levels, we present the mean junction residence time tr as a function of stress for each junction type in Fig. 6.In all cases, the mean residence time decays monotonically with increased stress.In fact, all data sets obey an inverse proportion law of the form where B r and τ c,∞ are fitting parameters, values for which are given in Table 3.The parameter τ c,∞ is the critical stress at which tr → ∞ .In other words, when τ < τ c,∞ no junction breaking events occur, meaning that the dislocation network is fully immobilized by junctions.When τ > τ c,∞ , junctions begin breaking and plastic flow ensues.Hence, we can recognize τ c,∞ as the quasi-static flow stress.To further assess the validity of this asymptotic extrapolation, we plot τ vs. 1/ tr in Fig. 6b.According to Eq. ( 5), the data should exhibit linear scaling when plotted in this way.We see that all junction types exhibit a linear trend as 1/ tr → 0 , indicating the correct asymptotic behavior. (4) Table 3 Fitting constants obtained for all junction types To establish a more granular understanding of the junction residence times, we present probability histograms for the residence time distributions at four different stress levels in Fig. 7, and the distributions of junction residence time are fitted by Ŵ-distributions.These distributions show a similar form at all stress levels, with the distribution shifting self-similarly to lower residence times as the stress increases.Note that we have no theoretical reason for selecting a Ŵ-distribution, we simply find that it fits the data reasonably well.

Link-length distributions
Next we seek to establish a connection between the arrangement and topology of the dislocation network, and the junction parameters we have identified above.To do so, we focus on the so-called link-length distribution of the dislocation network, n(l) (Lagneborg and Forsen 1973;Sills et al. 2018).A dislocation link, is a section of dislocation line which terminates on either end at physical nodes.Physical nodes are points in the dislocation network where either three or more dislocations intersect (e.g., the nodes of a junction) or a single dislocation changes its glide plane (e.g., due to crosssip or collinear junction formation).The link-length distribution is defined such that n(l)dl is the number density per unit volume of links of length l.Accordingly, the dislocation density is given by ρ = ∞ 0 ln(l)dl and the number of links per unit volume is where p(l) is the probability density function of link length distribution.The mean link length can be related to the forest density by the relation where φ is a dimensionless constant (Sills et al. 2018).The product φρ f can be thought of as the "effective forest density", because it quantifies what portion of the forest dislocations are actively involved in the dislocation network (e.g., connected to dislocation links).
While the previous work by Sills et al. considered the link-length distribution for the entire network, we can also instead define separate link-length distributions for each slip system.In our case here, we are interested in using the link-length distribution to establish a flow rule.Accordingly, we only need to include links from the primary slip system in the link-length distribution, since only primary links experience resolved shear stress in our simulations (e.g., forest links never generate net plastic strain).We have analyzed the link-length distribution for the primary slip system obtained in our random forest simulations as follows.For each stress level, all link lengths obtained after reaching steady-state were collected together in order to obtain an average link-length distribution from the whole simulation.The resulting distributions are shown in Fig. 8 at the same stress levels as Fig. 4. At all stress levels, the link-length distributions exhibit strong exponential character, except for at very large link lengths (same observation was made by Sills et al.).To compare these results with the exponential distribution, we note that N = ρ p / l in our simulations, leading to the form By fitting Eq. ( 8) to the each link length distribution, we can estimate the value of φ .The resulting φ values are shown in Fig. 9a as a function of stress for each junction type.We see distinct curves for each junction type, with φ ranging from 0 to 0.4.At the limit of τ → ∞ , φ → 0 because no junctions form at all so the network is comprised of one ( 6) single link of length l .Interestingly, when we instead plot φ as a function of κ in Fig. 9b, we see a strong correlation with the data from all junction types collapsing together.In particular, we see an appoximately quadratic relationship with fitting function which is shown as a solid line in Fig. 9b.

Flow rule
Having now established a relationship between junction formation and the structure of the dislocation network (i.e., how κ and φ are related), we are able to formulate a physics-based flow rule.This is accomplished as follows.First, we introduce the notion of a link area sweep rate function, Ȧ(τ , l) , which quantifies the rate at which a link of length l sweeps out area under resolved shear stress τ .Since plastic strain fun- damentally results from the area swept out during dislocation motion, Ȧ(τ , l) provides a physics-based linkage between the dislocation flow rule γ (τ ) and the link-length dis- tribution description of the dislocation network.Specifically, the total plastic strain is the superposition of the swept area from all links in the dislocation network, giving a plastic strain rate expression of the form: In order to make use of Eq. ( 10), we must establish a functional form for Ȧ(τ , l) .For our DDD simulations here, we use an isotropic, linear mobility law, meaning that the dislocation velocity v is proportional to τ/B , where B is the isotropic drag coefficient.The area sweep rate should scale as vl.However, if the stress on a link is not great enough to break either junction which anchors it, the area sweep rate is zero.Combining these ideas, we introduce the area sweep rate function where τ b (l) is the link-length-dependent breaking stress for the junctions bounding the link, which can also be thought of as a back-stress since it is being subtracted from the applied stress.Basic scaling arguments give that the breaking stress of an isolated junction should follow the form (Kubin 2013) where β ij is a dimensionless constant which depends on the type of junction (similar to a ij ) and l is the length of the link connected to the junction.β ij can be thought of as the non-dimensional "intrinsic" strength of a junction which is independent of the link lengths connected to the junction.
To proceed further we must determine β ij , which we can do by asking the following question: for the dislocation network, what is the average stress at which a junction will break?Each junction is bounded by two links in the primary slip system.According to Eq. ( 12), the longer of the two links will precipitate junction failure, since the junction strength scales as 1/l.Thus, the average junction strength is dictated by the expectation of the maximum link length when two links are randomly sampled from the dislocation network, E[max{l i , l j }] (assuming that the link lengths connected to a junction are uncorrelated).For an exponential distribution, it turns out that the result for this calculation is exceedingly simple: E[max{l i , l j }] = 3 2 l (Ross 2019).Hence, for our simulations where the exponential distribution is always satisfied, we expect that the average critical stress for breaking junctions is The junction breaking stress is equal to the flow stress at the quasi-static limit, τ c,∞ , since junctions are the only resistance to plastic flow in our systems.In other words, ( 10) we can equate Eq. ( 13) with the generalized Taylor relation for the quasi-static flow stress, Eq. ( 1); doing so while making use of Eq. ( 7) results in the relationship where φ c is the value of the φ parameter in the quasi-static limit (obtained by determin- ing κ at τ c,∞ using Eq. ( 4) and then making use of Eq. ( 9)).To determine a ij , we use the quasi-static flow stresses obtained from our asymptotic analysis of the junction residence times, τ c,∞ , in conjunction with the generalized Taylor law.The resulting values are given in Table 3.Finally, we obtain β ij for each junction type using Eq. ( 14), and these values are also given in Table 3.Now we are able to compute the stress-dependence of the plastic strain rate using Eq. ( 10) by substituting in Eq. ( 11).We numerically evaluate the integral in Eq. ( 10) over a range of stresses using scipy.integrate.quad in Python, and making use of Eqs. ( 4), ( 8), and ( 9).In Fig. 10, we compare the predictions of Eq. ( 10) with DDD results by plotting γ as a function of τ for all four junction types.The link-length-based model captures the form of the data well, in addition to the differences between the different junction types.This fit is captured without any free parameters in the model, all parameters were obtained from prior analyses, demonstrating the self-consistency of our approach and model formulation.At very low stresses, the plastic strain rate is zero.As the stress increases, the plastic strain rate gradually increases.At very high stresses the plastic strain rate converges to a linear dependence on stress, indicating dislocations are exhibiting free-flight dynamics with no resistance besides the drag force Bv.Note that the plastic strain rate curves in Fig. 10 are the result of both the Ȧ function and the evolu- tion of the link-length distribution with increasing stress through the φ parameter and ( 14) its junction-specific dependence on κ .To further validate our analyses, we plot the Ȧ values for all links in our simulations for each junction type as a function of τ l in Fig. 11, along with the fit from Eq. ( 11).All data at different stresses and link lengths collapse together into a line, as expected.In Fig. 10, the shape of the flow rule is similar for all junction types near the critical stress.To explore this further, we plot in Fig. 12 the plastic strain rate data from Fig. 10 with each dataset shifted by its respective τ c,∞ value, i.e., plotting γ vs. τ − τ c,∞ .At low stresses near τ − τ c,∞ = 0 all of the datasets collapse together, indicating that the behav- ior at the onset of plastic flow is universal to all junction types.Note that at high stresses the datasets diverge, since the dislocation networks evolve differently with stress for each junction type.

Discussion
Through our analysis above, we connected the formation of junctions with the evolution of the dislocation network, and then in turn with the plastic strain rate.In the end, we were able to explain the plastic strain rate in terms of basic dislocation processes (junction breaking and dislocation gliding) with no free parameters.As a further demonstration of the self-consistency of our analysis, it is interesting to compare the two critical stress quantities we obtained in our analysis.During analysis of the junction formation fractions, we extracted what we argued was the mean junction strength, τ µ .By analyzing the junction residence times, we also extracted the quasi- static flow stress, τ c,∞ .In Fig. 13 we plot these two strength quantities with respect to each other for each of the junction cases.The dashed line shows that the mean junction strength τ µ scales in proportion to the the quasi-static flow stress τ c,∞ , as expected since the mean junction strength controls the flow stress.This scaling  3 demonstrates the self-consistency of our analysis and also bolsters our interpretation of the junction formation fraction as deriving from the junction strength distribution.
Another interesting connection between dislocation junctions and the dislocation network is the relationship between κ and φ shown in Fig. 9b.In the original work where the φ parameter was proposed by Sills et al. (2018), it was noted that φ evolved with strain but no mechanistic reason was given.Our analysis here reveals that evolution of φ is connected to the probability of junction formation κ .κ and φ are both geometric factors that describe mean length scales in the system.Specifically, 1/ √ κρ f is the mean junction spacing in the glide plane and 1/ φρ f is the mean link length of the network.One may expect that these two lengths should scale in proportion with each other.However, Fig. 9b indicates that as the mean junction spacing decreases ( κ increases), the mean link length decreases ( φ increases) at a faster rate.The ori- gin of this nonlinear scaling is the junction strength distribution and how it is sampled by the applied stress, as we show in the Appendix.If the mean junction strength was independent of applied stress, we would see that instead κ increased more rapidly than φ .However, as the stress increases, so too does the mean junction strength in the network since weaker junctions are "screened out" by the stress.This causes κ to increase more rapidly than φ.
Our results indicate that the junction residence time follows a complex probability distribution which is sensitive to applied stress and junction type.In order for a junction to break, the links connected to the junction must sweep out sufficient area so that the plastic work done exceeds the energy reduction associated with formation of the junction and the energy increase due to elongation of the links as the area is swept (Kubin 2013, I. Duan, R.B. Sills, Crossed-state bowing and the strength of binary dislocation junctions.In-Preparation).The rate at which the dislocations bow out and sweep area thus dictates the junction residence time.This is a complex process which is sensitive to the lengths of the links bounding the junction, the level of the applied stress, and also the strength of the junction itself (which is controlled by the orientations and Burgers vectors of the lines forming the junction).Additional research is necessary to understand the specific form of the residence time distribution function observed in Fig. 7. Despite this complexity, we find that to a good approximation the mean residence time obeys the inverse proportion law in Eq. ( 5), especially in the asymptotic limit of the quasi-static flow stress.
One drawback of the random forest model is that the forest system must be free of resolved shear stress.Otherwise the dislocation forest will multiply during the simulation, making the results difficult (or impossible) to analyze (Devincre et al. 2006).On the other hand, under multi-slip loading the forest system will experience a nonnegligible resolved shear stress.It is well established that the resolved shear stress on the forest system affects the junction strength (Dupuy and Fivel 2002;Wu et al. 2016; I. Duan, R.B. Sills, Crossed-state bowing and the strength of binary dislocation junctions.In-Preparation).Hence, since our results above are largely a consequence of the strength of junctions, we expect that they will also be sensitive to the stress on the forest system.We note that this issue also applies to computation of interaction coefficients for the generalized Taylor relation.Future efforts should seek to better understand the influence of the stress on the forest system.
Another shortcoming of the random forest model is that it assumes that the dislocation network remains topologically simple, comprised solely of dislocation links bounded by binary junctions.However, in reality the network is topologically more complex than this (Akhondzadeh et al. 2021) for numerous reasons, one of which is the formation of ternary junctions (Bulatov et al. 2006;Madec and Kubin 2008).A related issue is the assumption that collinear "junctions" remain intact during plastic flow.Unlike other junction types, collinear junctions do not have a physical junction dislocation that ties the junction nodes together topologically.This means that if the network evolves considerably after a collinear junction has formed, the junction nodes may "lose each other" making it impossible for the junction to "break" by recombination of the junction nodes.Further research is necessary to clarify this possibility and its impact on plastic flow.
Finally, we compare the a ij values here with those from the literature.Comparing our values in Table 3 with those from the literature (Devincre et al. 2006;Kubin et al. 2008;Madec and Kubin 2017), our values are on the order of 2× larger.One explanation is the sensitivity of a ij to the core energy model and the dislocation core radius value used in the simulations.Together the core energy and core radius dictate the dislocation line energy, which directly affects a ij ; a larger line energy means that the dislocation is more resistant to bow out, leading to larger a ij values.As an example, we present in the Appendix results using a core radius of r c = 4b compared to the value of r c = 0.1b used in Table 3.The interaction coefficients with r c = 4b are about half of the interac- tion coefficients with r c = 0.1b , which demonstrates the fact that the value of interaction coefficient is sensitive to the selection of core radius.In some sense, any value for a ij can be obtained by changing the core energy or core radius.Clearly this indicates that a ij coefficients are not as well-defined as often implied in the prior literature.Another difference between our results and those from the literature is that we use an extrapolation to the quasi-static limit (of vanishing strain-rate) to determine a ij , whereas prior studies used results from finite strain-rate simulations.

Conclusions
We have used random forest DDD simulations to extract fundamental quantities that dictate dislocation network evolution during plastic flow, including the junction formation fraction, the junction residence time, and the area sweep rate per dislocation link.By specifying the stress state in the simulations rather than the strain rate, we were able to quantify the stress-dependence of each of these quantities.The stress-dependence of the junction formation fraction κ was shown to closely follow a Q-function, which we argue derives from the junction strength distribution.The fitting constants for the Q-function can thus be interpreted as the mean and standard deviation of the junction strength distribution.As stress increases, fewer and fewer junctions are strong enough to form, leading to a decaying κ function.The mean junction residence time tr was shown to follow an inverse proportion law.Extrapolating this law enabled us to extract the critical stress in the quasi-static limit, τ c,∞ .We then analyzed the junction link-length distributions, demonstrating that they followed an exponential distribution defined by parameter φ , which was shown to be a function of κ .Using this critical stress in conjunction with the junction strength scaling law, we quantified the average non-dimensional junction strengths β ij and interaction coefficients a ij .These quantities were used to define a stress τ b which governed the area sweep rate Ȧ(l, τ ) for each link.This area sweep rate was finally used to obtain the flow rule for each junction type, γ (τ ).
From beginning to end, our analysis connects the behaviors of individual junctions ( κ and t r ) to the emergent behavior of the whole system ( γ (τ ) ).This connection is enabled by our network (link-length distribution) interpretation of the dislocation ensembles.For example, if we had instead only focused on the dislocation densities rather than the link-length distributions it would have been difficult to attain a physics-based, parameter-free flow rule.However the picture is still incomplete because our analysis does not give a quantitative theory for predicting the multiplication rate of the dislocation The 4/3 exponent obtained in Eq. ( 23) is not consistent with the exponent of ∼ 2 obtained from our DDD data, indicating some of the details in this derivation need refinement.However unlike with Eq. ( 18), Eq. ( 23) exhibits the correct trend that φ increases more rapidly than κ .Furthermore the quantity ≈ 2.8 based on our data of Lomer junctions and simulation conditions, which is in rough agreement with the value of 1.2 we obtained in Fig. 9b.

Sensitivity of a ij to core radius
Another set of DDD simulations were performed using a core radius of r c = 4b (compared to r c = 0.1b for results given above), and the core energy parameter is still 0 GPa (no core energy).The dislocation densities of forest and primary slip systems are the same as the previous DDD simulations.The fitting constants of these new simulations are listed in Table 4, where the flow stress of the Hirth type is slightly below zero, indicating it is too small for us to compute with the uncertainty in our analysis.The interaction coefficients with r c = 4b in Table 4 are about half of the interaction coefficients with r c = 0.1b in Table 3, which dem- onstrates the fact that the value of interaction coefficient is sensitive to the selection of core radius.
Geometric factor f converting ρ to ρ f Now consider two ways of computing dislocation densities.The forest density ρ f is given by where N is the number of dislocations penetrating the area S of the primary glide plane.And the total density ρ is given by where L is the total length of dislocation lines in the volume of V.These two densities can be related by introducing a geometric factor f: For computing the value of f, consider a forest dislocation segment penetrating the primary slip plane, the line direction of the forest dislocation will form an angle θ with the direction perpendicular to the primary slip plane.Assuming θ obeys a uniform distribu- tion in the range of [0, π 2 ] (since the forest dislocations are randomly oriented), the geo- metric factor f can be computed by averaging cos θ in the range of [0, π 2 ]:

Fig. 3
Fig. 3 Plastic strain versus time for several stress levels with Lomer junctions

Fig. 6 aFig. 7
Fig. 6 a Mean junction residence time tr as a function of resolved shear stress τ for each junction type.Markers indicate results from DDD simulations, solid lines are fits to Eq. (5).b Transformed plot showing that τ vs. 1/ tr is linear at small stresses, confirming that the residence time asymptotically scales as tr ∝ 1/(τ − τ c,∞ ) where l is the mean link length.It was recently shown by Sills et al. using DDD simulations of uniaxial deformation that the link-length distribution obeys an exponential distribution of the form(Sills et al. 2018):

Fig. 10
Fig. 10 Plastic strain rates as a function of resolved shear stress for a collinear, b glissile, c Hirth, and d Lomer junctions, showing DDD data and predictions from Eq. (10)

Fig. 11
Fig. 11 Area sweep rates Ȧ as a function of τ l for a collinear, b glissile, c Hirth, and d Lomer junctions.Solid lines are Eq.(11) with the β ij values from Table 3

Fig. 12
Fig. 12 Plastic strain rate data from Fig. 10 shifted by critical stress τ c,∞ for each respective junction type

Table 1
Simulation parameters used for DDD simulations

Table 2
Slip systems for each junction type and stress tensors used for loading

Table 4
Fitting constants obtained for all junction types with a core radius of r c = 4b