Cell structure formation in a two-dimensional density-based dislocation dynamics model

Cellular patterns formed by self-organization of dislocations are a most conspicuous feature of dislocation microstructure evolution during plastic deformation. To elucidate the physical mechanisms underlying dislocation cell structure formation, we use a minimal model for the evolution of dislocation densities under load. By considering only two slip systems in a plane strain setting, we arrive at a model which is easily amenable to analytical stability analysis and numerical simulation. We use this model to establish analytical stability criteria for cell structures to emerge, to investigate the dynamics of the patterning process and establish the mechanism of pattern wavelength selection. This analysis demonstrates an intimate relationship between hardening and cell structure formation, which appears as an almost inevitable corollary to strain hardening itself. Specific mechanisms such as cross slip, by contrast, turn out to be incidental to the formation of cellular patterns.

Plastic deformation by dislocation motion is generally associated with dislocation patterning, leading to formation of heterogeneous dislocation arrangements. If multiple slip systems are active, dislocations form cellular structures where dislocation depleted 'cell interiors' are surrounded by dislocation rich 'cell walls' (Szekely et al., 2002). Such cell structures show an almost universal scaling behavior ('law of similitude') which is independent of loading condition, material or temperature: the characteristic pattern wavelength λ is proportional to the mean dislocation spacing (mds) ρ −1/2 0 where ρ 0 is the spatially averaged dislocation density, and inversely proportional to the applied stress (Rudolph, 2005;Sauzay and Kubin, 2011): λ ∝ ρ −1/2 ∝ 1/τ ext . This behavior results directly from fundamental scaling invariance properties of dislocation systems as discussed by Zaiser and Sandfeld (2014). Recent investigations by Oudriss and Feaugas (2016) indicate an even stronger form of the similitude principle according to which the components (cell walls, cell interiors) of cell structures obey the similitude principle separately, such that the wall thickness λ w is related to the wall dislocation density by λ w = Cρ −1/2 w and the cell dislocation density to the cell size λ c = Cρ −1/2 w , in such a manner that the proportionality coefficients C are identical. We note in passing that, under very specific conditions which may be the exception rather than the rule (namely, deformation of fcc crystals with the loading axis oriented along a [100] direction) fractal cell patterns with a wide spectrum of length scales may emerge . However, even in these exceptional cases, the length scales defined by the upper and lower boundaries of the fractal scaling regime of cell sizes obey the "law of similitude" as shown by ; Zaiser (1998).
Numerous models have been proposed for dislocation cell structure formation. Early models often relied on visual analogies of dislocation patterns with other patterning phenomena and adopted equations drawn from other realms of science (e.g. spinodal decomposition (Holt, 1970) and chemical patterning as described by reaction-diffusion models (Walgraef and Aifantis, 1985)). These equations were adapted to dislocations in a manner that, seen with malevolent eyes, might be envisaged as a mere re-labeling exercise. It is not easy to see how, if at all, such models for the specifics of dislocation topology, dislocation motion and dislocation interactions -for instance, it is immediately evident that the fundamental mode of dislocation motion under stress is not diffusion but directed glide. In recent years, efforts have been made to match chemical patterning inspired models more closely to actual dislocation processes, by distinguishing slip systems (Pontes et al., 2006) and providing physically motivated reaction terms (Aoyagi et al., 2013). However, in all these models the problem remains that diffusion terms do not appropriately describe the glide of dislocations, which needs to be described by transport terms that are of a hydrodynamic rather than of a diffusion-like character.
Discrete dislocation dynamics (DDD) simulation provides a powerful alternative to phenomenological adhoc models. DDD simulations faithfully represent the kinematics and interactions of dislocations and should be well suited for modelling dislocation pattern formation. While existing simulations (Hussein and El-Awady, 2016;Madec et al., 2002) indicate that simulations of systems sufficiently large to allow for a quantitative investigation of pattern morphology alongside a reliable determination of pattern wavelengths may still be challenging, such limitations may be overcome with time simply due to the expected increase in available computing power.
However, from an epistemological point of view the ability to provide a more or less faithful in vitro simulation of a real process should not be confounded with understanding: a sufficiently complex may encompass, besides essential, a large amount of redundant features and it may not be easy to decide which features of the collective dynamics are at the core of a collective phenomenon such as dislocation cell structure formation, and which are incidental to it. Rather than pursuing accuracy in detail, our own modelling strategy therefore is heavily poised towards simplicity -while at the same time we make sure that the most essential kinematic features and the structure of the interactions are represented correctly. Mathematical simplicity of the model allows us to obtain some results in an analytical or semi-analytical manner, and makes the essential features of the dynamics obvious. To this end we rely on a most basic version of density based dislocation dynamics in multiple-slip conditions. We start from the model used by Groma et al. (2016); Wu et al. (2017) for analysing the conditions for pattern formation in single slip, and generalize this to symmetrical double slip along lines proposed by Yefimov and Van der Giessen (2005) and Limkumnerd and der Giessen (2008). This framework not only provides us with some degree of analytical tractability but also with a solid theoretical foundation: The equations we use have been derived from statistical averaging of the underlying discrete dynamics (Groma et al., 2003;Valdenaire et al., 2016) and can be related via variational calculus to the statistically averaged energy functional of the dislocation system Groma et al. (2016);Zaiser (2015). Moreover, predictions obtained with these equations for size-dependent deformation in small samples and/or constrained geometries have been shown to be in quantitative agreement with discrete dislocation dynamics simulations (Yefimov et al., 2004;Yefimov and Van der Giessen, 2005). This makes us confident that the mathematical framework we used indeed captures some of the essential features of dislocation dynamics under load.
We note that other, more complex versions of density-based continuum dislocation dynamics have been applied to the patterning problem. Some of these approaches consider geometrically necessary dislocations only (Chen et al., 2013;Limkumnerd and Sethna, 2008). However, during the early stages of deformation the dislocations in the cell walls have near-zero net Burgers vector: they are predominantly not geometrically necessary dislocations. Application of such models to early stages of cell structure formation is therefore possible only if the spatial resolution is well below the actual dislocation spacing such that Burgers vectors do not cancel out. If one makes this numerical effort the results can be impressive, see the work of Xia and El-Azab (2015) which reproduces the morphology of three-dimensional dislocation patterns in impressive detail. A more coarse grained model that allows for co-existence of dislocations of different Burgers vector in the elementary volume but nevertheless captures effects of three-dimensional curvature was proposed by Sandfeld and Zaiser (2015). An interesting work which is conceptually close to the present one was recently published by Grilli et al. (2018). These authors consider two models which allow for dislocations of different Burgers vector in the same elementary volume, which are described by a set of densities obeying transport equations and applied to labyrinth-like patterns emerging under cyclic loading. These works are conceputally more complex than the present one, as they consider three-dimensionally curved dislocations Sandfeld and Zaiser (2015), distinguish various orientations Grilli et al. (2018), and including essentially three-dimensional processes such as junction formation Grilli et al. (2018) and cross slip (Xia and El-Azab, 2015). While these approaches are interesting in their own right, we demonstrate in the present paper that the added complexity is actually not essential for cell structure formation or dislocation patterning as such. In the following we first briefly introduce the governing equations of our model and then provide a stability analysis that allows us to establish necessary conditions for cell pattern formation. We show the results of numerical simulations of the evolution equations and compare our findings to experimental data. Finally we provide a conclusion where we discuss implications of our findings in view of some commonly held ideas regarding the nature of dislocation patterns and the requirements for their formation.

Model Equations
We consider crystal deforming in plane strain where two orthogonal slip systems are active. System 1 has Burgers vector b 1 = be x and slip plane normal n 1 = e y , and system 2 with Burgers vector b 2 = be y and slip plane normal n 2 = e x . The shear strains on the two slip systems are denoted as γ 1 and γ 2 . The plastic distortion is then given by (1) We define the plastic strain ε pl and plastic rotation ω pl as the symmetric and anti-symmetric parts of β pl . These are given by where γ = γ 1 + γ 2 and ω = γ 1 − γ 2 . Both slip systems contain straight parallel edge dislocations gliding in the directions of the respective Burgers vectors. We assume that each system contains equal numbers of positive and negative dislocations with the corresponding dislocation densities denoted as ρ ± 1/2 where the upper label distinguishes positive and negative dislocations, and the lower label distinguishes the two slip systems. Positive dislocations move under the action of a positive resolved shear stress in the positive Burgers vector directions, and negative dislocations move under the same shear stress in the negative Burgers vector directions, v ± 1/2 = ±v ± 1/2 b 1/2 /b where v ± 1/2 are scalar velocities.
In the spirit of defining a minimal model, we neglect dislocation reactions (which are anyway not expected to occur for energetic reasons), dislocation multiplication and annihilation. The dislocation densities are thus conserved quantities which obey the continuity equations ∂ρ + The dislocation velocities for these four types of dislocations are assumed to be linearly proportional to respective, effective shear stresses T s i where the index i ∈ {1, 2} distinguishes the two slip systems and s ∈ {−1, 1} distinguishes the two signs of the dislocations: In these equations, M 0 is a dislocation mobility coefficient (inverse drag coefficient). A closed mathematical model is then specified by relating the effective shear stresses to the dislocation densities. In line with the single-slip model of Groma et al. (2016); Wu et al. (2017), we consider the effective driving stresses T s i to result from the combination of sign-dependent local driving stresses τ s,dr i and friction stresses τ s,f i : The driving stresses combine the resolved shear stress τ i in the respective slip system with corrections describing short-range dislocation interactions associated with the mutual arrangement of individual dislocations (dislocation correlations) according to We discuss the three stress contributions in this equation separately: 1. The resolved shear stress τ i arises from the superposition of stresses caused by external tractions and internal stresses associated with the plastic eigenstrains -in other words, it is found by solving a standard elastic-plastic problem. The considered slip geometry has the peculiarity that this stress is the same in both slip systems and equals the xy component of the stress tensor, τ 1 = τ 2 = σ xy . In our calculations, we consider a bulk system with periodic boundary conditions and calculate this stress using from the plastic strain γ using a Green's function formalism as in (Wu et al., 2017;Zaiser and Moretti, 2005): where τ ext is a spatially constant external stress arising from remote tractions acting on the infinite contour, and G is an interaction kernel function with the Fourier transform G is the shear modulus of the material, ν is Poisson's ratio, and k x and k y are components of the Fourier wave-vector with modulus k. 2. The 'back stresses' τ b i stem from the mutual correlation of dislocations of the same sign and counter-act their accumulation. For single slip on some slip system i the back stress is given by where D is a non-dimensional factor of the order of unity and ρ = ρ + + ρ − is the total dislocation density on the considered slip system. The local excess density κ i is given by the difference of positive and negative dislocation densities and relates to the slip gradient on the slip system i via For multiple slip situations as considered here, Limkumnerd and der Giessen (2008) use a statisticalmechanical model of the density cross correlation functions to derive instead of Eq. 10 the superposition relations where θ i j are the angles between the Burgers vectors (slip directions) of slip system pairs. For the geometry considered here, cos θ i j = δ i j and hence, Eqs. (10) and (12) are equivalent. 3. Accordingly, we assume the 'diffusion stresses' τ d i to be given by where A is another nondimensional factor of the order of unity. The terminology 'diffusion stresses' is used because this stress, if inserted via Eqs (7), (6), (5) into the transport equations Eq. (4), give rise to a diffusion-like contributions to the evolution of the total dislocation densities ρ i .
All three stress contributions can be derived from the energy functional of the dislocation system, as discussed in detail by Groma et al. (2016), hence, they are associated with stored energy contributions. It remains to specify the friction-like stresses τ s,f i . In generalization of the expression derived by Groma et al. (2016) for single slip, we assume these stresses in the form where the latent hardening matrix H i j describes slip system interactions. The dependency on the κ i accounts for the fact that excess dislocations cannot be pinned by dislocations of the same slip system (the net force on the excess cannot become zero), for more details see Wu et al. (2017).
For the present system, the resolved shear stresses induced by a dislocation in both slip systems are equal, hence, it is reasonable to set H ii = H i j = 1 leading to where ρ = ρ 1 + ρ 2 is the total dislocation density. The friction stresses are of a different nature from the driving stresses: they represent friction-like stresses that are associated with dissipated, not with stored energy contributions. While these stresses arise naturally from direct averaging of the dislocation interactions, they cannot be derived from an energy functional but need to be added 'by hand' to an energy-based formalism where they enter in terms of a non-trivial, nonlinear mobility function with a mobility threshold Groma et al. (2016). The functional form of these stresses is that of Taylor stresses; in physical terms, they represent the mutual trapping of dislocations into dipolar or multipolar configurations. Their dependency on the κ i reflects the fact that the presence of an excess of dislocations of one sign implies reduced pinning of the majority and enhanced pinning of the minority population. Assembling all stress contributions, we find that the four dislocation density species under consideration fulfill, under the assumption that the local effective stress is positive and the system is everywhere in the flowing phase, the respective continuity equations The strains γ i evolve according to Before we proceed to analyze the model equations, it is important to comment on the nature and meaning of the non-dimensional parameters A, D, and α which enter the model in addition to the physical constants G, b, ν, and the mobility coefficient M 0 . All three parameters A, D, and α characterize correlations in the positions of individual dislocations and can in principle be evaluated in terms of integrals over dislocation-dislocation correlation functions, see their derivations by Groma et al. (2003), Limkumnerd and der Giessen (2008), Groma et al. (2016) and Valdenaire et al. (2016). All these parameters are of the same order of magnitude as they characterize the spacings of close dislocations whose positions, owing to their mutual interactions, are strongly correlated. Specifically, α is proportional to the characteristic spacing of dislocations that have trapped each other into dipolar or multipolar configurations, measured in units of the typical spacing of dislocations of the same slip system in the surrounding of a given spatial point -of course, as such α is nothing but the well known Taylor factor. If the dislocation arrangement is thought of as an assembly of isolated dipoles of height h, then α = (8π(1 − ν)(h √ ρ), but in more general circumstances, this factor needs to be modified to account for the influence of dislocations surrounding the dipole. The parameters A and D have an analogous interpretation, but 'probe' different aspects of short-range interactions: While α mainly captures the trapping effect of dipole-like interactions, D characterizes the interactions between dislocations of the same sign in piled-up configurations, which cause a net stress if there is a gradient in the 'geometrically necessary' density κ. Finally, A which controls the 'diffusion stress' accounts for the fact that dipoles and multipoles have finite extension, such that dislocation density cannot localize down to arbitrary narrow scales. In summary, all three factors are proportional to spacings of individual dislocations, with α mainly characterizing the spacing of slip planes of adjacent dislocations, D spacing of dislocations of the same sign in piled up configurations, and A the extension of dipoles and multipoles in glide direction.
Understanding the physical nature of the constants α, D, A is also beneficial for the physical interpretation of the respective stress contributions. Breaking of dipoles and formation of new ones is a dissipative process that occurs as soon as the local stress exceeds the dipole breaking stress, hence, the associated stress contribution has friction-like characteristics. Piling up dislocations against an obstacle, by contrast, leads to storage of energy that can be recovered if the stress causing the pile up is removed or reversed, hence, the associated energy contribution enters an appropriately averaged internal energy functional. The same is true for the work expended in compressing or expanding dipolar and multipolar configurations. It is in line with these intuitive arguments that, upon formal statistical averaging of the elastic energy of a dislocation system as performed by Zaiser (2015), the resulting density based functional allows to recover through variational calculus both the 'back stress' and the 'diffusion stress' but not the 'friction stress' Groma et al. (2016).

Reference state
We consider pattern formation first in an analytical framework where we focus on infinitesimal perturbations of a spatially homogeneous reference state where ρ i,s = ρ 0 /4 ∀ {i, s} and γ i = γ 0 /2 ∀ i. At this stage we envisage loading by a temporally constant applied stress τ ext . Depending on the level of stress, two situations need to be distinguished: (i) If τ ext < αGb √ ρ 0 then all velocities in the reference state are zero, hence, γ = 0 is constant in space and time and ρ i,s = ρ 0 /4 is a stationary solution of the evolution equations that is stable with respect to infinitesimal perturbations. (ii) If τ ext > αGb √ ρ 0 we are in a flowing phase. In this case the dislocations move with homogeneous and stationary velocity v 0 = M 0 b(τ ext − αGb √ ρ 0 ) and the slip system strains increase linearly in time, ∂ t γ i =γ 0 /2 = ρ 0 bv 0 /2. The stability of this state is analyzed in the following.
In our analysis we have a choice of variables. Instead of the four densities ρ i,s we may use the total and excess dislocation densities on the two slip systems, ρ i = s ρ i,s and κ i = s sρ i,s . Furthermore, instead of the excess dislocation densities we may alternatively consider the slip variables γ i which relate to the former via κ i = b i ∇γ i /b 2 . This is the choice we make, i.e., we consider the problem in terms of the four variables ρ i , γ i , i ∈ {1, 2}.

Dimensionless scaling
In the following we switch to a dimensionless formulation which helps to see the influence of all model constituents more easily. Only the final results are stated here, detailed information and derivations can be found in (Sandfeld and Zaiser, 2015;Zaiser and Sandfeld, 2014). We define the scaling relations between quantities with physical units and their dimensionless counterparts (indicated by a tilde) as τ = C ττ (for stresses), ρ s = C ρρ s (for dislocation densities), x = C xx (for lengths), and γ = C γγ , with the scaling factors Furthermore, we scale velocities in units of C v = M 0 bC τ , which implies a scaling for time according to t = C tt with C t = C x /C v . In non-dimensional form the equations of motion become where we have dropped the tildes on all variables and introduced the notations ∇ i = (b i .∇/b),D = D/α, and A = A/α. The scaled stress kernel is given byT = T/(αρ 0 ) and has in scaled variables (k → k/ √ ρ 0 ) the Fourier transform

Linearized evolution equations
We now write down the equations of evolution for small perturbations δρ i , δγ i of our reference state ρ i,0 = 1/2, γ i,0 = γ 0 /2. In linear approximation these perturbations are given by Defining the state vector δq = [δρ 1 , δγ 1 , δρ 2 , δγ 2 ] and using the Fourier Ansatz δq = q(k) exp(ik.r), we write these equations in matrix form: We now first investigate two simple cases where the eigenvalues can be computed analytically in a straightforward manner.

Symmetrical case
We first study the eigenvectors and eigenvalues of this matrix for the symmetrical case k x = k y = k/ √ 2. The matrix M simplifies to The eigenvectors of this matrix have the structure q 1 = ±q 2 where q 1 = [δρ 1 , δγ 1 ], q 2 = [δρ 2 , δγ 2 ]. We first consider the "-" case. The matrix equation then reduces to M − .q 1 = Λ − q 1 where The eigenvalues fulfil the characteristic equation Since τ ext > 1 in the flowing phase and bothÃ andD are positive, both roots of this equation have negative real parts for all k and τ ext , hence, no instability can occur. In the "+" case we get M + .q 1 = Λ + q 1 where The eigenvalues then fulfil the characteristic equation An unstable wavelength band may in that case occur if 1 ≤ τ ext and 8τ ext (3/2 − τ ext ) > T 0Ã . This band is comprised between the wavelengths k = 0 and k

Fluctuations along the cube axis
Next we consider the case where the fluctuation wave vectors are aligned with the x axis, k x = k, k y = 0 (the opposite case is symmetry equivalent). The matrix M simplifies to The characteristic equation is obtained by setting the determinant of the matrix M − ΛI to zero. Expanding the determinant with respect to the last column gives the straightforward result Curves Λ(k) are shown in Figure 1, for fluctuations in the glide directions and along the slip system symmetry axis. The instability occurs for fluctuations aligned with the slip systems, the wavevector of maximum amplification corresponds, for the parameters given in the Figure,

Condition for instability: physical interpretation
Since instability occurs first in [10] directions, the condition for instability to occur is, in non-dimensional representation, simply given by τ ext < 5/4 or, in dimensional units, τ ext < (5/4)αµb √ ρ 0 . To understand the physical nature of this condition, we define the total (scalar) flux of dislocations on slip system i in the homogeneous reference state as The derivative of the total flux j i with respect to the slip system dislocation density ρ i is then given by For the present case where ρ i = ρ 0 /2 we thus find that the dislocation density derivative of the total dislocation flux turns negative when τ ext < 5/4αµb √ ρ 0 which is precisely our instability criterion. We are, hence, dealing with a variant of a basic instability that has long been studied in hydrodynamic models of traffic flow, see e.g. Gerlough and Huber (1975). Importantly, no other terms in the evolution equation but the flux term and the friction-like stresses -which represent the isotropic hardening due to dislocation density accumulation -are needed to observe this instability which is, hence, a quite generic feature of dislocation dynamics.
We now proceed to numerical analysis of the full, nonlinear equations in order to show how these findings are reflected in the developing patterns.

Numerical analysis
We have performed a numerical analysis of the evolution equations for two different types of initial conditions, namely (i) a spatially uncorrelated Gaussian white noise of small amplitude and (ii) a localized small perturbation in the origin of the coordinate system. We implement periodic boundary conditions in x and y for the stresses and for the dislocation fluxes on the two slip systems. For the stress evaluation we use a Finite Element framework with periodic displacement boundary conditions. As initial conditions we use ρ ± (r, t) = ρ 0 /2 + δρ ± (r, t) where 1 and we consider two types of perturbation δρ ± : (i) a Gaussian white noise of unit amplitude and (2) a localized Gaussian 'blob' of width l = ρ −1/2 0 located at the center of the simulation cell. The system is loaded by imposing a constant external stress and keeping it fixed throughout the simulation.
The time evolution of the Fourier coefficients of the emergent patterns is shown in Figure 2 for both cases. The emergent patterns are dominated by fluctuations with wave-vectors oriented along the symmetry equivalent [01] and [10] lattice directions. From the initial growth rates of the discrete Fourier modes ρ(k) we deduce growth factors defined as Λ(k) = ∆ ln ρ(k)/∆t. Comparison with the analytical predictions for fluctuations oriented along [10] and [11] lattice directions shows good agreement. The wavelengths of the fully developed patterns match closely (within 20%) the predictions of linear stability analysis for the wavelength of the mode with maximum amplification. At longer times, satellites appear at multiples of the dominant wavelength and the Fourier spectrum assumes a grid-like pattern, indicating a non-sinusoidal periodic pattern with long-range order. While the initial growth rates of Fourier components are similar for localized and distributed perturbations, the ordering tendency seems to be more pronounced if patterning starts from a single localized perturbation ( Figure  2, bottom). The mode of growth depends on the initial conditions, see Figure 3: in case of a spatially distributed noise the emergent patterns have a crossed stripe-like character. If we use a localized perturbation as initial condition, two perpendicular walls start growing from the perturbation and then the wall pattern spreads into a grid-like pattern. The characteristic wavelength of the emergent pattern is, however, independent of the growth mode.
An interesting question concerns the applicability (or not) of the well-known composite model to our simulation data. According to the composite model as originally formulated by Mughrabi (1983), long-range internal stresses associated with slip heterogeneities develop in such a manner as to homogenize deformation. Regions of enhanced dislocation density (cell walls) have a higher local flow stress, accordingly, plastic slip is reduced in these regions. In regions of reduced dislocation density, the flow stress is reduced and slip is enhanced. The compability requirements between both kinds of regions imply presence of geometrically necessary dislocations which, so the model, create long range internal stresses that offset the flow stress differences. Ultimately, in quasi-static deformation one expects the local stress to everywhere match the local flow (friction) stress such that deformation can then proceed in a compatible manner: . Note that this relation is expected to hold independent of the length scale of the pattern: The 'composite' of the original composite model is considered in the spirit of classical composite mechanics which does not know about size effects. The composite model has some important corollaries. For instance, it can be seen immediately that patterning does, in the composite model, always lead to softening (reduction of flow stress) in comparison with the homogeneous reference state: Evaluating the spatial averages .... and noting that the because of stress equilibrium τ(r) = τ ext , we find that in the patterned state because of the triangular inequality τ ext = αGb √ ρ < τ ext,0 = αGb √ ρ 0 where ρ 0 = ρ is the homogeneous reference density. This finding is supposed to hold independently of the morphology or of the length scale of the heterogeneous patterns (Zaiser (1998)). Looking at the strain patterns in our simulations we find that they match the expectations: Strain is increased in the cell interiors and decreased in the cell walls. If we look at the internal stress patterns in our simulations, however, a more complex behavior is found. The internal stresses do not exhibit a strict correlation with the plastic strain, or with the dislocation density, see Figure 5.
To quantify the deviation from the composite model, we note that according to the composite model, in non-dimensional variables we expect the local internal stresses and dislocation densities to obey the relation where the angular brackets denote spatial averages. Figure 6 shows that a positive correlation which however is significantly below the value expected according to the composite model, exists only during the initial stage of patterning. This correlation actually decreases as patterns are formed and ultimately drops to zero. For patterns emerging from a localized perturbation, there is an additional complication since the correlation oscillates as walls are formed sequentially. Either way, in the fully developed pattern there is no appreciable correlation between local stress and local dislocation density. This raises the intriguing question how the patterns can deform compatibly. The shortfall is made up by the length scale dependent stress contributions τ b i (r) and τ d i (r) which may be considered non-local, strain and dislocation density gradient dependent generalizations of the classical composite model. This points to a limitation of the composite model which assumes an entirely classical composite mechanics framework: If applied to patterns that are heterogeneous on the micrometer scale, where in other composite systems size effects start to become relevant, composite models which neglects non-local stress contributions might systematically under-estimate the flow stress of heterogeneous dislocation arrangements, see also the discussion of strain gradient effects in the composite model by Mughrabi (2001)

Relation to experimental observations
At first glance a plane-strain slip geometry with two perpendicularly intersecting slip systems as studied in the present idealized model seems unrealistic. However, a quite faithful realization of this situation can be found in early deformation stages of ionic solids with KCl crystal lattice structure. This structure consists of two interlaced fcc sub-lattices containing the K + and Cl − ions, respectively. If the crystal is subjected to a uni-axial stress state with the stress axis oriented along the [100] crystal lattice axis, deformation can take place on four symmetrically oriented slip systems which form two conjugate pairs, namely the (110) (101)[101] systems. We make the following observations: 1. The active slip systems are such that, for tension along a [100] lattice axis aligned with the x axis, the conjugate pairs of active slip systems produce plane strain states in the xy and xz planes, respectively. 2. The slip systems in a conjugate pair intersect at right angles. Their mutual interactions are comparatively weak (forming a junction produces, in line tension approximation, no net energy gain). By contrast, there are strong interactions between pairs of slip systems belonging to different conjugate pairs, leading to significant latent hardening. 3. As a consequence, during the early stage of deformation a symmetry breaking takes place where deformation is taken over by one conjugate pair of slip systems while the second pair becomes inactive (Schwerdtfeger et al., 2010). This situation quite faithfully matches the slip geometry assumed in our simulations.  Streb and Reppich (1973); the insert has been taken from the simulation shown in Figure 3 and scaled according to the average dislocation density in the experimental image.
can be estimated by noting that the dislocations forming a dipole stem from independent sources, hence, it will be of the order of (1/5) mds which, with a typical dislocation density of ρ = 3 × 10 11 m −2 , translates into a spacing of the dislocations in the dipoles of the order of about 0.35 µm, well above the atomic spacing. Hence, annihilation of dislocations is not expected to be a relevant process here. The walls are formed by the mutual trapping of dislocations into dipole-like configurations (friction stress). They are stabilized by two effects that mutually compensate each other: On the one hand, excess of dislocations of positive sign pushes against the wall from one side ('pile up stress') , on the other hand, the dislocations within a dipole push each other back ('diffusion stress'). As a consequence we see a wall consisting of polarized dipoles, with positive and negative dislocations accumulating on opposite sides of the wall. The width of the walls, the corresponding width of the cells and the dislocation spacings are all in good agreement with the experimental observations. This can be seen in Figure 7, right, where a piece of the simulated dislocation density pattern could, after re-scaling to the dislocation spacing in the experiment, be seamlessly pasted into the experimental image.
We also investigate whether our patterns match the similitude principle in the strong form proposed by Oudriss and Feaugas (2016). To this end we study one-dimensional density profiles taken along the slip directions and define, for a given profile, the wall dislocation density ρ w i of wall i as the dislocation density at the corresponding density maximum and the channel dislocation density ρ c i as the dislocation density in the corresponding density minimum. Left and right wall boundaries x l i and x r i are defined as the locations where the dislocation density takes the respective values (ρ w i − ρ c i )/2 and (ρ w i − ρ c i+1 )/2. The width of wall i is then evaluated as λ w i = x r i − x l i and the width of channel i as λ c i = x l i − x r i−1 . Figure 8 shows lengths λ c,w as well as pattern wave-lengths λagainst the corresponding densities ρ c,w for different values of the average density ρ 0 . As can be seen, the data are well represented by a common fit function λ c,w = C √ ρ c,w with C ≈ 6, in good agreement with the findings of Oudriss and Feaugas (2016). Also the overall relationship between pattern wavelengt λ = λ c + λ w and total dislocation density ρ 0 matches well the experimental data of Oudriss and Feaugas (2016) (full data points in Figure 8). We thus conclude that our model is consistent with the strong similitude principle as observed by Oudriss and Feaugas (2016).

Discussion and Conclusions
We have presented a very simple model of dislocation cell structure formation in a 2D setting with two perpendicularly intersecting slip systems. Despite its simplicity, the model can be considered a elementary representation of dislocation processes in a real system, namely a crystal with KCl lattice structure deformed uni-axially along a cube axis. We find formation of cellular dislocation patterns with a cell size of the order of about 10 mean dislocation spacings. The patterns obey the similitude principle: their wavelength is proportional to the dislocation spacing and inversely proportional to the stress at which they form. The simplicity of the 2D model, which can not account for dislocation multiplication, does not allow us to consider strain hardening. However, if we impose a higher overall dislocation density ρ 0 , then deformation requires an accordingly higher stress that scales in proportion with √ ρ 0 , and similitude is maintained. It is instructive to discuss our findings in relation to commonly held viewpoints on dislocation patterns: (i) It is an often expressed viewpoint (see e.g. Madec et al. (2002); Xia and El-Azab (2015) that cross slip is essential for dislocation cell structure formation. However, it is easy to see that in KCL structures, as in our simulations, this mechanism is irrelevant since there is only one (110) slip plane for each [110] slip vector, hence, there are no cross-slip planes. Nevertheless, formation of cellular dislocation patterns is observed regularly in these structures and our simulations -where cross slip is excluded by construction of the model -provide an excellent match to the observed cellular patterns. We therefore conclude that cross slip is, in the end, incidental to dislocation patterning. (ii) The composite model predicts that a patterned dislocation arrangement deforms at a stress that is strictly below the stress needed for deforming a homogeneous reference arrangement. This assumption is predicated upon a classical treatment of internal stresses that does not allow for strain gradient dependent effects. Even within the classical continuum mechanics framework, it is clear that dislocation patterns or strain patterns of general morphology in general produce internal stress patterns that do not directly match the strain/dislocation patterns as required by the composite model, compare our Figures 5 and 4. In fact, for the present slip geometry a match between stress and dislocation patterns would be possible only if the dislocation patterns would form with a [11] orientation which they do not. Deformation compatibility must therefore be ensured by other means that cannot be described by standard continuum mechanics. Such effects are also needed to understand pattern wavelength selection. In our model these effects are provided by the gradient dependent stress contributions τ b and τ d , in other models (Sandfeld and Zaiser, 2015) a similar role is played by curvature related terms. (iii) The only essential requirement for patterning in our model is that, for a given stress, the local dislocation flux is a decreasing function of local dislocation density. Many models of work hardening fulfill this requirement for a wide range of deformation parameters. We therefore conclude that, if dislocation density evolution is described by appropriate transport equations, patterning is an expected feature of dislocation dynamics. Our investigation can be easily generalized to a wide range of stress-velocity laws in order to provide guiding principles that allow to decide under which deformation conditions heterogeneous patterns may form. It thus provides an important complement to microstructure-based plasticity models as proposed e.g. by Castelluccio and McDowell (2017) which investigate the impact of self-organization of dislocations into mesoscale structures on the macroscale deformation behavior under complex loading paths.