Pattern formation in a minimal model of continuum dislocation plasticity

The spontaneous emergence of heterogeneous dislocation patterns is a conspicuous feature of plastic deformation and strain hardening of crystalline solids. Despite long-standing efforts in the materials science and physics of defect communities, there is no general consensus regarding the physical mechanism which leads to the formation of dislocation patterns. In order to establish the fundamental mechanism, we formulate an extremely simplified, minimal model to investigate the formation of patterns based on the continuum theory of fluxes of curved dislocations. We demonstrate that strain hardening as embodied in a Taylor-type dislocation density dependence of the flow stress, in conjunction with the structure of the kinematic equations that govern dislocation motion under the action of external stresses, is already sufficient for the formation of dislocation patterns that are consistent with the principle of similitude.


Introduction
Work hardening during plastic deformation of crystalline solids is associated with significant changes in dislocation microstructure. The increase in dislocation density on the specimen scale is accompanied by the quasi spontaneous emergence of regions of low dislocation density and clusters of high dislocation density which to a large extent persist upon unloading. These metastable structures are denoted as dislocation patterns. Despite a significant degree of morphological variation depending on slip geometry and loading mode (e.g. cell [1] or labyrinth structures [2], dislocation accumulation in veins [3] or walls [4]), these patterns are characterized by some fairly universal scaling relationships. These relationships are commonly referred to as 'law of similitude' or 'similitude principle'. They relate the characteristic length Λ of deformation-induced dislocation patterns to the applied stress τ ext at which they have formed, and to their average dislocation density ρ: The wavelength Λ is proportional to the dislocation spacing and inversely proportional to the applied stress, Λ = D/ √ ρ = DaGb/τ ext arXiv:1501.03527v2 [cond-mat.mtrl-sci] 16 Jan 2015 To resolve these controversies it would be desirable to develop computational models which allow to activate or de-activate physical micro-mechanisms at will, and thus to decide which features of dislocation motion are essential and which ones are incidental to pattern formation -an approach chosen for instance by Madec et al. [14]. Three-dimensional discrete dislocation dynamics (3D DDD) simulations provide a complete description of dislocation motion and interactions on the level of individual dislocation lines and might therefore offer a suitable methodological framework for such an investigation. Clearly developed dislocation patterns have, however, rarely been observed in 3D DDD simulations. The reason for this is simple: If we require that, in order to avoid boundary effects, the periodic simulation cell of a bulk 3D DDD simulation should at least encompass ≈ 5 pattern wavelengths in each direction, and using D = 10 as the ratio between pattern wavelength and dislocation spacing (or mesh length of the dislocation network), we find that in order to observe patterning, the size of the simulation should exceed the dislocation spacing by a factor of about 50 in each spatial direction. This is very much at the upper end of the system sizes attainable by state-ofthe art 3D DDD simulations. In fact, to our knowledge, there is only a single example of a 3D DDD simulation of this size in the published literature [15] and even there, the achieved strain of < 2% may have been insufficient to form fully developed dislocation patterns. In the study of Madec et al. [14], a tendency towards dislocation clustering concomitant with the formation of dislocation-depleted regions is reported, but the scale of these features is comparable to the simulated volume and hence no definite conclusions on pattern morphology can be drawn. Patterning can be efficiently simulated using two-dimensional discrete dislocation dynamics (2D DDD) where dislocations are envisaged as straight parallel lines, but such models require phenomenological rules to represent genuinely three-dimensional processes associated with segment curvature, such as dislocation multiplication and junction formation. Attempts have been made to formulate and parameterize such rules using 3D simulations ('2.5D DDD'), which results in patterns that are consistent with the similitude principle [16]. However, it is clear that this approach owing to its geometrical restrictions cannot give access to dislocation patterns in general deformation geometries.
We expect that the above mentioned problems of 3D DDD, which mostly relate to the huge computational cost of the method, can be overcome by resorting to density-based continuum dislocation dynamics (CDD) frameworks. A suitable CDD framework must be capable of representing the kinematics of curved dislocations in a geometrically correct manner and should allow to include the essential features of dislocation interactions. Several attempts in this direction have been reported in the literature. Some studies focus on subgrain formation (e.g. [17]) -a process associated with geometrically necessary dislocations only -and can therefore not capture the initial stages of dislocation patterning which are associated with the clustering of statistically stored dislocations of zero net Burgers vector. Other investigations focus on specific mechanisms such as the sweeping of narrow dislocation dipoles by curved dislocations as envisaged by Kratochvil et al. [18] which may only be relevant in specific deformation conditions such as fatigue.
In the present manuscript we undertake first steps towards a generic model that captures the coupled dynamics of statistically stored and geometrically necessary dislocations while accounting for the specific kinematics of curved dislocation lines. To this end we use the kinematic framework of Hochrainer's CDD theory to represent fluxes of curved dislocations in a continuum setting [19,20,21,22]. We use this framework to investigate the following question: Given the kinematics, what are the minimal ingredients required for the emergence of dislocation patterns? In Sect. 2 we formulate our kinematic CDD framework together with a minimal model of dislocation interactions and a specification of boundary conditions. Scaling the system reveals the structure of the equations and allows to analyze the critical conditions required for patterning in Sect. 3. We then use numerical simulation to investigate the properties of the emergent patterns (Sect. 4) and their evolution with stress during strain hardening. We conclude with a discussion of the physical mechanism of dislocation patterning and point out perspectives for future work.

Continuum dislocation dynamics theory
As kinematic framework on which we base our numerical investigation we use the CDD theory formulated by Hochrainer and co-workers [20,21,22]. The field variables describing the dislocation microstructure on a given slip system with normal vector n and Burgers vector b = bs are the total density ρ t , the dislocation density vector , and the curvature density q t . The latter quantity can be envisaged as an averaged product of dislocation line density and line curvature; its integral V q t dV = (2π)N d (V ) over a volume V yields the number of dislocation loops contained in the volume, hence, q t can also be envisaged as a loop density. We define the GND density vector as = −(1/b)∇ n γ where γ is the plastic slip on the considered slip system, and ∇ n = ∇ · (I − n ⊗ n) is the Nabla operator in the slip plane where I is the rank-2 unit tensor. We assume that dislocation motion occurs by crystallographic glide with average (scalar) velocity v. Hence, Orowan's equation [23] gives the evolution of γ as ∂ t γ = ρ t bv.
We note that our definition of differs from the definition of the dislocation density vector κ in [20,21,22]; both formulations are related by = −n × κ where n is the slip plane normal vector. The modified definition leads to an equivalent, but more intuitive formulation of the kinematic equations which are now given by These equations can be considered lowest-order terms of an alignment tensor expansion [22]. We close this expansion by expressing the curvature density vector Q (1) and the rank-2 alignment tensor A (2) in terms of the variables ρ t , q t and .
The curvature density vector Q (1) in (3) plays very much the same role as the dislocation density vector in (2). As pointed out by Hochrainer et al. [21], Q (1) can be approximated as the product of the dislocation density vector and the mean curvature q t /ρ t , In the vicinity of a homogeneous and isotropic dislocation arrangement, this approximation becomes exact. The rank-2 alignment tensor A (2) has trace Tr A (2) = ρ t and for an isotropic dislocation arrangement, A (2) = 0.5ρ t I. The deviatoric part A (2) dev characterizes imbalances in the distribution of dislocation characters (e.g. predominance of edge vs. screw dislocations). For a weakly anisotropic arrangement containing a fraction of geometrically necessary excess dislocations, we follow the proposal of Monavari et al. [24]) and write where l ρ is the unit vector in direction of the GND vector (the slip gradient) and Φ is a function with the properties Φ(0) = 0 and Φ(1) = 1. The approximation for A (2) proposed by Hochrainer et al. [21] is tantamount to setting Φ(u) = √ u which appears inappropriate since it leads to a formulation which is non-analytic at u = 0, i.e., in the vicinity of a homogeneous and isotropic dislocation arrangement. We therefore follow the suggestion of Monavari et al. [24] who from a Maximum Entropy argument deduce Φ(u) ≈ u(1 + u)/2.

Dislocation interactions and dislocation velocity
In a density-based model, dislocation interactions enter in different ways depending on their range. Interactions on scales above the characteristic dislocation spacing (in the following referred to as long-range interactions) can be represented in terms of eigenstresses associated with a heterogeneous distribution of the plastic strain ( [25,26]) or, equivalently, with the presence of geometrically necessary dislocations (non-vanishing dislocation density vector ). Interactions on scales below the characteristic dislocation spacing, on the other hand (in the following referred to as short-range interactions) cannot be described in this manner, for the obvious reason that the continuous dislocation density fields and the associated strain field cannot meaningfully represent heterogeneities below the scale of the single-dislocation spacing. These short-range interactions are described in the standard manner in terms of a friction-like yield stress, which we assume to depend on the local value of the dislocation density according to the Taylor relationship. Additional terms relating to line curvature or to short-range interactions between geometrically necessary dislocations might be considered (see the discussion of such terms by Zaiser & Sandfeld [27]) but, in the spirit of formulating a minimal model, here we exclusively use the Taylor relationship to describe the shortrange dislocation interactions which control the flow stress. The dislocation velocity is then computed from the glide component of the Peach-Koehler force as where B is the drag coefficient. In this expression, τ is the resolved shear stress in the slip system, computed by solving the elastic boundary value problem for a body containing the eigenstrain field γ. This stress captures both the effect of the externally applied boundary tractions and/or constraints, and the long-range dislocation interactions associated with a heterogeneous plastic strain or, equivalently, with the dislocation field .

Boundary conditions
We study two versions of the model, namely (1) a situation where long-range interactions are absent but short-range interactions are present, and (2) a situation with long-range interactions but no short-range ones. In both cases we envisage a layer of plastically active material of thickness h where deformation occurs in single slip (a slip zone or slip band). The layer is assumed parallel to the xy plane which is also the slip plane, and the slip direction is the x direction. For numerical simulations, periodic boundary conditions are assumed in the x and y directions.
In case (1) we assume that the plastically active layer forms a thin film loaded by homogeneous shear tractions acting on its unconstrained top and bottom surfaces. In case (2) we consider the layer to be embedded into an infinite block of elastic material which is loaded remotely, again in a spatially homogeneous manner. The stress field in the plastically active layer is then given by case (1), where τ ext is the shear stress due to the externally applied tractions and the Fourier transform of the kernel G is given by (see e.g. [28]) The convolution integral in (7) can from this relation be evaluated in standard manner by a double Fast Fourier Transform.

Dimensionless scaling
To analyze the behavior of the model it is convenient to switch to a non-dimensional representation which allows to identify the independent parameters. We define the scaling relations between quantities with physical units and their dimensionless counterparts (indicated by a tilde) as τ = C ττ (for stresses), ρ t = C ρρt (for dislocation densities) and x = C xx (for lengths), with the scaling factors where ρ 0 = ρ t (t = 0) is the mean initial dislocation density. The same scaling applies to the total dislocation density and the dislocation density vector. The scaling relation for the lengths implicates that lengths are measured in multiples of average dislocation spacings. For the velocity we find for where H 0 = αGb √ ρ 0 /τ ext . The scaling of velocity implies a scaling of time, since the scaling factor of velocity can be understood as the ratio between the characteristic length C x and a characteristic time C t . Furthermore, the fact that the dislocation curvature has the dimension of an inverse length suggests to scale q t in proportion with C ρ /C x . Finally, we observe that the plastic slip can be envisaged as the product of the Burgers vector and the area per unit volume swept by dislocations. These observations lead to the derivative scaling An important implication of the scaling relationships for stress, length and dislocation density is that any patterns that may emerge from a simulation based on the scaled equations automatically fulfill the similitude principle (cf. [27]) -provided that it can be demonstrated that the patterns do not depend on initial and boundary conditions (i.e., on system size and deformation history).
In the scaled variables the kinetic equations read: where the velocityṽ is given by (10) and the quantitiesQ (1) andÃ (2) are obtained by replacing ρ, q t and in (4) with their scaled counterparts. The scaled stress is given bỹ whereG wherek = k2 x +k 2 y . The dislocation density vector relates to the scaled strain bỹ q = −∇ nγ , and the evolution equation forγ is It is interesting to note that, in the scaled formulation, the model has only two non-dimensional parameters. The first, H 0 = aGb √ ρ 0 /τ ext , measures the fraction of the externally applied stress that is needed to overcome short-range dislocation interactions. The second parameter, G 0 = Gbhρ 0 /τ ext measures the relative magnitude of the longrange dislocation interactions compared to the external stress.

Linear stability analysis
To analyze stability of the model, we envisage a spatially homogeneous reference state characterized by the valuesρ t = 1,˜ = 0,qt =qt 0 ,γ =γ 0 , with dislocation velocity v 0 = 1 − H 0 , and study the evolution of small fluctuations around this reference state. It turns out to be convenient to use the setρ t ,qt,γ as independent variables. We define fluctuations around the homogeneous state bỹ The operator∆ n =∇ n ·∇ n is the 2D Laplacian in the slip plane, and the corresponding term in the expression for the curvature density can be envisaged as the curvature of the geometrically necessary dislocations. In linear approximation, we find that The linearized evolution equations for δρ t , δq t , and δγ are then given by Eigenvalues of this matrix are Λ = Λ(q t 0 , H 0 , G 0 ,k). We first consider the hypothetical case where long-range stresses are appreciable but short-range stresses are absent (H 0 = 0, G 0 > 0). From the Routh-Hurwitz criterion it follows in this case that for any values ofk, all Eigenvalues have negative real part. The system is thus absolutely stable and no patterning can occur. More interesting is the opposite case, G 0 = 0, H 0 > 0. In this case the Routh-Hurwitz criterion yields the condition 1 > H 0 > 6/7 for eigenvalues with positive real part to occur, with the largest positive values occurring for H = 13/14. Within this unstable regime, the critical eigenvalue branch is real for allk > 0, indicating that fluctuations on all scales become undamped while only the homogeneous statek = 0 remains marginally stable. Finally, if both H 0 > 0 and G 0 > 0, the instability condition remains unchanged, however, in this case long-wavelength fluctuations are damped and only fluctuations above a critical value k c can become unstable. The regime of unstable k vectors is given bỹ An example showing the critical eigenvalue branch as a function of the wave-vector is shown in Fig. 1. In all cases where unstable modes occur, the unstable mode is characterized by antiphase fluctuations of γ and ρ t (regions of high dislocation density correspond to regions of low strain and vice versa). The same is true for the fluctuations of q t and ρ t , hence, the dislocation curvature (and thus the multiplication of dislocations) decreases in regions of high dislocation density and increases in regions of low density. This is consistent with the well known fact that, in heterogeneous dislocation patterns, multiplication occurs by the expansion of segments in the dislocation-depleted regions (channels or cell interiors).

Simulation of pattern formation
The analysis of the previous section can establish conditions for a homogeneous dislocation arrangement to become unstable. However, it cannot give access to the morphology of the emerging patterns which depends on the non-linear dynamics of the system. Moreover, the instability is in linear approximation not characterized by any dominant wavelength, hence, also the characteristic length of the patterns remains undetermined by the linear stability analysis. To access pattern morphology and pattern scale (if any), we resort to numerical simulation of the model.

Numerical methods
For the numerical solution of the scaled CDD evolution equations (12)- (14), (17) we use a Galerkin finite element scheme together with an implicit time integration. The spatial discretization consists of approximately 50 triangular elements in each spatial direction. In order to properly represent first and second derivatives within the finite element scheme we use third order Lagrangian shape functions. For reasons of numerical stability a small viscous term was added to each of the evolution equations, e.g. in (12) we replace the divergence term∇ n · (ṽ˜ ) with∇ n · (ṽ˜ + ∇˜ ), where is chosen sufficiently small such that the physical behavior of the equations is not affected but numerical oscillations are suppressed. For each time step we update the velocity field (i.e. the Taylor term) and then evolve the dislocation fields. One of the benefits of the dimensionless scaling is that we can solve systems with widely different parameters based on the same finite element mesh; the small but artificial viscous terms for stabilizing the numerical scheme also have the same influence on the results for widely different systems, which improves our control over the numerics.

Initial values
Initial values for the dislocation fields should be consistent with representing a system of curved and connected lines in a coarse-grained manner. To this end we proceed in analogy with the coarse graining of a stochastic point process into a superposition of Gaussian 'blobs': we consider a system of discrete dislocation loops of radius r and 'smear out' each loop perpendicular to its line tangent by using a Gaussian standard distribution function (s √ 2π) −1 exp(−ξ 2 /2s 2 ) where the standard deviation s characterizes the width of the distribution and ξ is the distance perpendicular to the loop. After having obtained ρ l t for one loop in this way we determine the corresponding dislocation density vector l by multiplication of ρ l t with the radial unit vector, and the curvature density by division of ρ t by r. Initial values for the fields are then obtained by randomly distributing such loops and summing up the loop fields ρ l t , l , and q l t . An example of a resulting initial structure is shown in Fig. 2. An alternative method to define initial conditions is to define random scalar fields γ, q t and ρ t as superpositions of randomly located Gaussian 'blobs' and determine as the gradient of γ, which automatically fulfils the solenoidality condition divα of the Kröner-Nye tensor. To make sure that our results do not depend on the construction of initial conditions, both methods have been used in the following simulations.

Parameter values
We consider the system which is defined by the parameters given in table 1. The scaled quantities were obtained from the physical values according to Eqs. (9), (11).

Simulation results
We first look at statistically averaged quantities (volume averages) before we then analyze the system in terms of microstructural patterning.
The mean dislocation density ρ t and loop number N d are shown in Fig. 3 as functions of time. We observe that initially the dislocation density quickly increases until it saturates at a value that is about four times the initial value. At the same time the total number of dislocation loops N d stays exactly constant. This is expected since loop creation by Frank-Reaed sources and loop merging by mutual annihilation of segments from different loops are not accounted for in the present simple model. Accordingly, q t is a conserved quantity -a fact which is immediately obvious from (3).  Table 1; all fields are represented in a projection on the slip plane, instead of q t we show the more intuitive curvature q t /ρ t which is the inverse of the local dislocation curvature radius, The fact that N d stays perfectly constant in our simulation indicates that the numerical scheme performs well. In order to understand the saturation of the dislocation density, we note that the system can become stationary only by reducing the dislocation flux ρ t v -and thus the strain rate -everywhere to zero. This is possible if either the Taylor stress balances the external stress, or if the local dislocation density vanishes. This observation agrees with the fact that, in the dislocation-rich regions, the local dislocation density approaches the saturation density given by ρ sat = [τ ext /(aGb)] 2 = 2.99 · 10 14 m −2 where the Taylor stress equals the applied external stress. The dynamics of the system is further illustrated in the right graph of Fig. 3 which shows the average dislocation flux (the total strain rate) and the mean, maximum and minimum dislocation velocity. Initially, we see an increase in dislocation flux which is driven by the increase in dislocation density. This is followed by a decrease in dislocation flux caused by the decrease in dislocation velocity due to strain hardening. Even in the final state, there remain regions of non-zero dislocation velocity which, however, no longer produce strain since they are depleted of dislocations. Fig. 2, bottom, shows the CDD field variables in the stationary state; the supplementary movie additionally shows the time evolution of the CDD fields. Fig. 4 shows the concomitant evolution of the fluctuations of strain, dislocation density and dislocation velocity as characterized by the respective coefficients of variation (COV, defined as the ratio between variance and mean of the respective variable). In the initial period (cf. movie), the system is stable since the applied stress is significantly higher than the Taylor stress (H < H c = 6/7). Accordingly, fluctuations of strain or dislocation density are initially smoothed out as seen from Fig. 4. However, the increase in dislocation density drives the system towards the critical state which is reached, for the simulated parameters, at a density of about 2.2 × 10 14 m −2 . This point in time is marked by an arrow in Fig. 4. The ensuing instability is characterized by an increase in the fluctuations of strain and dislocation density, and an accelerated growth of the fluctuations in dislocation velocity. Once the critical state is reached, persistent dislocation density 'blobs' emerge from the fluctuating density distribution while the channels between these blobs become depleted of dislocations. In these channels, dislocations move rapidly (the maximum dislocation velocity increases, see Fig. 3) and curved line segments/loops expand and create additional line length which is deposited in the blobs. The emergent pattern is characterized by the following features (see Fig. 2
(ii) Dislocations that form the boundaries of these dense patches are aligned to each other and hence become GNDs.
(iii) The 'channels' between the dislocation rich regions exhibit the largest plastic strain.
Regions with high density, on the other hand, exhibit a reduced plastic strain.
(iv) The curvature is reduced inside the dislocation-dense patches but is high in the dislocation depleted 'channels'. The dislocation dense patches act like hard inclusions which constrain dislocation motion. Particularly high curvature is observed in spots where "corners" of the dense patches force the dislocations to assume strongly bent configurations.
(v) The emerging pattern morphology does not depend on the initial values: a smaller initial density is compensated by a longer simulation time (more multiplication due to loop expansion) leading to the same final density and morphology. Simulations with different type of initial conditions (random fields vs. smeared-out loops) also produce the same final morphology.
(vi) The emergent pattern morphology does not depend on the numerical discretization (as long as the mesh is small enough to be able to represent the patterns shape); the pattern borders in general do not coincide with the finite element faces.
(vii) Patterns formed in different simulation runs with the same shear stress exhibit approximately the same characteristic length.
The last observation is an auspicious feature of all our simulations -regardless of the chosen initial values, system size or material parameter. We analyze this behavior in the following section in terms of the 'law of similitude'.
Pattern evolution under changing stress and the principle of similitude In Fig. 5(a) we show density patterns in simulations where we increase the external stress in a step-wise manner, such that each stress step is followed by a relaxation step until a quasi-stationary configuration is reached before the next stress increment. We record the saturation density ρ sat after each relaxation step. The pattern morphology after relaxation remains essentially unchanged during the step sequence, however, the pattern size decreases with increasing stress while the dislocation density increases. To quantify this, we have determined the characteristic length of the patterns. This was done by use of the open source image analysis software ImageJ [29], which allows to determine the total number n, area fraction A i and perimeter S i of each density feature i of the density distribution in a semi-automatic way. The characteristic length is then obtained as the average value Λ = 4A i /S i of the characteristic sizes of all relevant 'density blobs'. These data are plotted in Fig. 5(b) together with the saturation density ρ sat . Fig. 5(a) shows examples of the pattern morphology for the largest and smallest  the measured pattern wave length by Λ = DaGb/τ ext with D = 9.5 we obtain good agreement for all data; thus the patterns are consistent with the similitude principle.

Discussion and Conclusions
We have formulated a model that describes the emergence of heterogeneous dislocation patterns that are consistent with the similitude principle. Linear stability analysis has demonstrated that the only requirement for the emergence of such patterns is the presence of short-range dislocation interactions as described by the Taylor relationship. Long-range internal stresses, on the other hand, are irrelevant for the occurrence of patterning -although they may influence the pattern morphology by suppressing long wavelengths.
The critical condition formulated in Sect. 3, 6/7 < H < 1, implies a lower and an upper boundary to the range of stresses where patterning can occur. If τ ext is too small (τ ext < τ f ), neither plastic flow nor dislocation patterning take place. This may appear trivial, but it formally demonstrates that the present model is inconsistent with the idea of quasi-equilibrium patterning due to energy minimization in the absence of externally driven plastic flow. If τ ext is too large (τ ext > (7/6)τ f ), i.e, if the deformation rate is too high or, more generally speaking, the rate-dependent contribution to the flow stress becomes too large, then patterning is also suppressed. This observation simply means that the flow stress must be predominantly (to at least 6/7) attributed to dislocation interactions, and not to interactions with the crystal lattice (rate-dependent drag component). Qualitatively this is consistent with the fact that heterogeneous dislocation patterns are almost generically observed in fcc metals, but not in bcc metals at low temperatures where the flow stress is controlled by the interaction of dislocations with the lattice.
It is interesting to compare the critical condition found in the present case with the results of Zaiser & Sandfeld [27] who carried out a similar analysis for a system of straight parallel dislocations. In that case, the instability criterion has a similar structure. In fact, the mechanism underlying the instability is in both cases the same: The dislocation density evolution equations have the structure of conservation equations. Therefore, in regions of reduced dislocation velocity, the dislocation density increases and vice versa. Due to strain hardening, a local increase in dislocation density in turn reduces the local dislocation velocity, leading to instability. We thus find that patterning is a direct consequence of the short-range dislocation interactions which control strain hardening and of the kinematics of dislocations. No other ingredients are needed.
What about dislocation multiplication? The present model naturally accounts for the increase of dislocation line length due to loop expansion. If we compare the resulting instability condition, 6/7 < H < 1, with the counterpart derived for straight parallel dislocations by Zaiser & Sandfeld [27], which can be written as 2/3 < H < 1, we find that in the latter case the unstable band is more wide. The difference is due to the fact that in the present model which considers curved dislocations, regions of low dislocation density develop increased curvature and thus dislocation multiplication localizes in the dislocation-depleted regions -a fact which is well known from in-situ observations of dislocation motion in heterogeneous dislocation patterns. At the same time, dislocation multiplication is reduced in dislocation dense regions and this effect tends to counter-act the instability.
Simulations of the model with increasing stress demonstrate that the patterns evolve in a manner that is consistent with the principle of similitude, i.e., the pattern wavelength is proportional to the mean dislocation spacing in the fully developed pattern, with a proportionality constant R ≈ 10, and decreases in inverse proportion with the stress. All these features are consistent with observations. At the same time, it is important to emphasize the very substantial limitations of the presented model. With this study we do not aim at realistically reproducing dislocation patterns observed in any real deformation experiment. Our aim is instead very much on the conceptual level: We want to present a minimal model which can be considered the skeleton of any realistic patterning model, while leaving it to future work to put in the realistic details. More realistic models will need to account for: • Boundary conditions typical of real deformation experiments. We chose rather artificial boundary conditions to represent, for the purpose of analysis, the possibility of a situation without long-range internal stresses (case 1), or with only long-range stresses (case 2). In reality, of course, both features will always be present.
• Though heterogeneous dislocation structures do form even in single slip, situations with multiple active slip systems are much more important. Accordingly, one needs to account for latent hardening terms coupling the dislocation densities on different slip systems in the expression for the short-range interaction stress, and for a more complex structure of the long-range internal stress field.
• The present model accounts for strain hardening in terms of the dislocation density increase that arises from loop expansion. However, in present form it does not represent a realistic hardening model: the reason is that according to (3) the number of dislocation loops is a conserved quantity. This leads to a parabolashaped hardening curve, τ f ∝ √ γ, which is not found in any real material. The reason is simply that in real materials the loop number increases by the action of Frank-Read sources, whereas at large strains and/or dislocation densities also the merging of loops by annihilation of encountering segments needs to be taken into account. All these features, while needed for obtaining realistic hardening curves, are however incidental to the patterning mechanism.
• In a realistic model, additional effects such as the influence of dislocation curvature/line tension on the dislocation velocity may need to be taken into account in an explicit manner.
As a consequence of these limitations, the blob-like pattern morphology observed in the present model differs from the morphology of dislocation patterns commonly observed in real samples. Other models proposed in the literature (see Sect. 1) produce much nicer patterns -but they cannot be shown to be in agreement with the similitude principle or to be based upon a realistic description of dislocation interactions. We are therefore convinced that the present approach, once the necessary model adjustments have been made, can serve as a robust framework for future investigation of dislocation patterning phenomena in a wide range of deformation settings.