{'fraction_non_alphanumeric': 0.03333333333333333, 'fraction_numerical': 0.01043453401943968, 'mean_word_length': 4.752898610311652, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
An Exactly Solvable Phase-Field Theory of Dislocation Dynamics, Strain Hardening and Hysteresis in Ductile Single Crystals

An exactly solvable phase-field theory of dislocation dynamics, strain hardening and hysteresis in ductile single crystals is developed. The theory accounts for: an arbitrary number and arrangement of dislocation lines over a slip plane; the long-range elastic interactions between dislocation lines; the core structure of the dislocations resulting from a piecewise quadratic Peierls potential; the interaction between the dislocations and an applied resolved shear stress field; and the irreversible interactions with short-range obstacles and lattice friction, resulting in hardening, path dependency and hysteresis. A chief advantage of the present theory is that it is analytically tractable, in the sense that the complexity of the calculations may be reduced, with the aid of closed form analytical solutions, to the determination of the value of the phase field at point-obstacle sites. In particular, no numerical grid is required in calculations. The phase-field representation enables complex geometrical and topological transitions in the dislocation ensemble, including dislocation loop nucleation, bow-out, pinching, and the formation of Orowan loops. The theory also permits the consideration of obstacles of varying strengths and dislocation line-energy anisotropy. The theory predicts a range of behaviors which are in qualitative agreement with observation, including: hardening and dislocation multiplication in single slip under monotonic loading; the Bauschinger effect under reverse loading; the fading memory effect, whereby reverse yielding gradually eliminates the influence of previous loading; the evolution of the dislocation density under cycling loading, leading to characteristic `butterfly' curves; and others.


Introduction
This paper is concerned with the formulation of an exactly solvable phase-field theory of dislocation dynamics, strain hardening and hysteresis in ductile single crystals. The present theory is an outgrowth of the phase-field model of crystallographic slip of Ortiz (1), which was restricted to a single dislocation loop moving through point obstacles of uniform strength and was based on a line-tension approximation. The present theory accounts for: an arbitrary number and arrangement of dislocation lines over a slip plane; the long-range elastic interactions between dislocation lines; the core structure of the dislocations described by means of a Peierls potential; the interaction between the dislocations and an applied resolved shear stress field; and the irreversible interactions with short-range obstacles and lattice friction, resulting in hardening, path dependency and hysteresis.
A chief advantage of the present theory is that it is analytically tractable, in the sense that the complexity of the calculations may be reduced, with the aid of closed form analytical solutions, to the determination of the value of the phase field at point-obstacle sites. In particular, no numerical grid is required in calculations. On the sole basis of the value of the phase field at the obstacle sites, the theory reconstructs, analytically and in closed form, the location of-possibly large numbers of-discrete dislocation lines, often in complex arrangements and undergoing intricate topological transitions, and endows each dislocation line with a well-defined core. Furthermore, the theory characterizes the equilibrium configurations of the dislocation ensemble directly, with the result that no short transients associated with artificial kinetics need to be resolved in calculations. By virtue of these attributes, the theory would appear to be advantageous in comparison to other recent phase-field models of dislocation dynamics (2; 3), which require the use of large computational grids and time integration over large numbers of time steps.
For definiteness, in this paper we restrict our attention to systems consisting of a dislocation ensemble moving within a single slip plane through a random array of discrete forest dislocations under the action of an applied shear stress. In the present theory, the dislocation ensemble populating the slip plane is represented by means of a scalar phase field. Specifically, the value of the phase field at a point of the slip plane simply records the extent of slip in quanta of Burgers vector. The phase field is, therefore, integer valued, and its value at a point may alternatively be regarded as recording the number of dislocations, with proper accounting for their sign, which have crossed the point. Thus, the individual dislocation lines may be identified with the lines over which the phase field jumps by one. The phase-field approach furnishes a simple yet efficient means of representing arbitrary dislocation geometries, possibly involving very large numbers of individual dislocations, as well as enabling the tracking of intricate topological transitions, including loop nucleation, pinching and the formation of Orowan loops.
The central objective of the theory is to characterize the evolution of the phase field. The theory is couched within a general variational framework for dissipative systems developed by Ortiz et al. (4;5;6;7;8), and accounts for energetic and kinetic effects. The energy terms contemplated by the theory include: the core energy of the dislocations, represented by a piecewise-quadratic Peierls potential (9); the longrange elastic interactions between primary dislocations and between primary and forest dislocations; and the energy of interaction with the applied resolved shear stress. The particular piecewise quadratic form of the Peierls potential adopted here lends itself to an effective analytical treatment based on the Fourier transform (9). The main effect of the Peierls potential is to endow the dislocations with a well-defined core energy and to coarse-grain the slip plane by suppressing wavelengths shorter than the lattice parameter.
At zero temperature and in the absence of irreversible processes, the equilibrium configurations of the dislocation ensemble follow directly from energy minimization. The resulting variational problem is strongly nonlinear owing to the all-important constraint that the phase field be an integer-valued function. Furthermore, the minimization problem is nonconvex owing to the multiwell structure of the Peierls potential. In addition, the variational problem is nonlocal, owing to the presence of long-range elastic interactions.
These attributes render the energy minimization problem mathematically non-trivial and confer its solutions a rich structure. Despite these difficulties, the choice of a piecewise quadratic Peierls potential lends the problem analytical tractability and, remarkably, a general solution to the problem can be found analytically in closed form.
In the present theory, hysteresis arises from the assumed irreversible short-range interactions between primary and forest dislocations and from lattice friction. The strength of some of these interactions has recently been investigated using atomistic and continuum models (10; 11; 12). Thus, we assume that the crossing of a forest-dislocation site by a primary dislocation always costs energy, regardless of the direction of the crossing. In the context of a phase-field representation, this is tantamount to assuming that any variation in the phase field at sites occupied by obstacles requires the supply of a certain amount of work, regardless of the sign of the variation. Physically, this work may be identified with the energy required to dissolve the dislocation-obstacle reaction product. Evidently, this form of interaction is irreversible and dissipative, and thus cannot be described by means of an energy function. It is possible, however, to develop an incremental variational framework which characterizes the evolution of the system by means of a sequence of minimization problems (4; 5; 6; 7; 8). The function to be minimized over a given step includes the incremental work of dissipation incurred as a result of obstacle crossings, and depends on the state of the system at the beginning of the step, which results in path-dependent and hysteretic behavior.
The theory predicts a range of behaviors which are in qualitative agreement with observation, including: hardening and dislocation multiplication in single slip under monotonic loading; the Bauschinger effect under reverse loading; the fading memory effect, whereby reverse yielding gradually eliminates the influence of previous loading. Simultaneously with the deformation and hardening characteristics of the system, the theory naturally predicts the evolution of the dislocation density. In particular, no independent equation of evolution for the dislocation density needs to be supplied.

Dislocation energies
In this section we derive explicit expressions for the energy of a dislocation ensemble contained within a slip plane in an elastic crystal. The dislocation ensemble is described by means of a phase field defined over the slip plane. In addition to the elasticity of the crystal, the slip plane is endowed with a piecewise-quadratic Peierls interplanar potential. The elastic interaction energy formulated in this section is a special case of a general class of energies for continuously distributed dislocation loops in isotropic elastic crystals derived by Ortiz and Xu (13). Extensions to anisotropic crystals have also been given by Xu (14).
We consider a crystal undergoing deformations characterized by a displacement field u i . Following Kröner (15), we begin by decomposing the displacement gradient in the additive form: where u i,j is the displacement gradient, or distortion field, and, here and subsequently, commas are used to denote partial differentiation. For simplicity, we shall specifically focus on crystallographic slip occurring on a single slip plane S. Under these conditions, the plastic distortion β p ij is supported on S and has the form where δ i is the displacement jump across S, m is the unit normal to S, and δ S is the Dirac distribution supported on S. The elastic distortion β e ij is assumed to be continuous. Neither the elastic nor the plastic distortion is required to be compatible, i. e., to be a gradient.
We shall additionally assume that the energy of the dislocation ensemble may be written as the sum of three terms: a core energy expressible in terms of a Peierls interplanar potential; the elastic interaction energy of the dislocations; and the energy of interaction with the applied stress field. These assumptions lead to the consideration of an energy functional of the form: Here, φ denotes the Peierls interplanar potential and the first term in (3) represents the misfit or core energy of the dislocations. In the second term, which represents the elastic interaction energy, c ijkl are the elastic moduli of the crystal. Finally, the third term represents the interaction energy between the dislocations and a self-equilibrated applied stress field resulting in a distribution of tractions t i over the slip plane. The aim now is to derive a phase-field representation of each of the terms in (3).

Elastic interaction energy
Suppose that the slip distribution δ i over S, or equivalently, β p ij , is prescribed. The corresponding elastic interaction energy may then be computed as follows. Insertion of (1) into (3) and minimization with respect to the displacement field yields the equilibrium equations The corresponding displacement field is: where G ki is the Green's function and (⋆) denotes the convolution operator. The distortion field then follows as and the elastic distortion as Finally, the elastic interaction energy E int is obtained by inserting (7) into the second term of (3). A convenient explicit expression for E int may be obtained by recourse to the Fourier transform. Using the convolution theorem, the Fourier transform of the elastic distortion follows from (7) aŝ where a superposed (ˆ) denotes the Fourier transform of a function, k i is the wavenumber vector, and the Fourier transform of the Green's function is determined by the relation: Finally, the interaction energy follows from an application of Parseval's identity, with the result: Here the symbol ( * ) denotes complex conjugation and we writê Mura (16) has shown that E int can also be expressed in terms of Nye's (17) dislocation density tensor. Indeed, it is readily verified from (8) and (9) that, if the plastic distortion β p ij is compatible, i. e., if it is a gradient, then the elastic distortions, and hence the elastic interaction energy, vanish identically.
In order to facilitate analysis we shall resort to several simplifying assumptions. Thus, for the special case of an isotropic crystal the Green's function reduces tô where µ is the shear modulus of the crystal, ν is its Poisson's ratio, and we write k = |k|. The interaction energy (10) then simplifies to where we write η i = k i /k. In addition, we shall adopt the constrained displacement hypothesis of Rice (18;19) according to which displacements, and the attendant shear resistance, take place predominantly in the direction of the dominant Burgers vector, i. e., where s = b/|b| is a constant unit vector and δ is a scalar function. For certain crystals, this conjecture has found support in atomistic calculations (20; 21; 22; 23). By virtue of these assumptions, (14) simplifies to where we have chosen axes such that the slip plane is parallel to the (x 1 , x 2 )-plane and the Burgers vector points in the direction of the x 1 -axis. Finally, for the particular case in which the slip distribution is confined to the plane x 3 = 0, the elastic interaction energy further reduces to which is the sought expression. By way of a simple illustrative and verification example we may consider the case of an edge dislocation pile-up, for which the slip distribution is of the form whence the elastic energy per unit length of dislocation follows from (17) as For a single Volterra dislocation at the origin,f (k 1 ) = 1/ik 1 and the elastic energy per unit length evaluates to where, in order to avoid logarithmic divergences, we have introduced lower and upper cutoff radii. We verify that (20) coincides with the well-known expression for the energy per unit length of a straight edge dislocation in an isotropic crystal (24). For a screw pile-up, a similar derivation yields: For a single screw dislocation at the origin the dislocation energy per unit length evaluates to which again coincides with the classical result (24).

Core energy
In the Peierls theory of the dislocation core, the interplanar potential φ(δ) is identified with the energy per unit area that results when two semi-infinite crystals are taken through a relative rigid displacement δ i . From symmetry considerations it follows that, in the absence of an applied field, the energy of the crystal attains minima when the displacement jump δ i is an integral multiple of a Burgers vector b i of the lattice, i. e., at where Z denotes the set of all integer numbers. These special slips determine the location of the wells of φ.
With a view to enabling the application of Fourier transform methods, we shall assume that φ is piecewise quadratic (9). In this model, the interplanar potential is taken to be of the form The moduli C jl may be determined by equating (24) with the energy per interatomic plane of a crystal undergoing a simple shear deformation of the form: where d is the interplanar distance. The result is: By way of example we may consider the case of crystallographic slip on a {111}-plane of an fcc crystal. Then a simple calculation gives d = a/ √ 3, with a the cubic lattice parameter, and In this expression c 11 , c 44 and c 12 are the three independent cubic elastic moduli.
If the slip is additionally constrained to take place in the direction of the Burgers vector, as in eq. (15), then the piecewise quadratic interplanar potential (24) reduces to the form where δ = δ · s and The function φ(δ) and its derivative, which gives the resolved shear stress as a function of slip, are shown in Fig. 1.

Computation of the energy from the phase field
At zero temperature, the stable configurations of the dislocation ensemble may be identified with the relative minima of the energy (3). We begin by attempting a characterization of the absolute mimimizers of (3), and defer the more ambitious program of understanding the relative energy minimizers as well. Inserting (17) and (28) into (3) leads to the energy: where is a normalized slip function, and In this expression, is the resolved shear stress field, and we write The normalized slip distribution ζ of interest now follows from the problem: In definition (30), the function ξ may be regarded as an integer-valued phase field defined on the slip plane S, with the property that ξ(x) equals the number of dislocations which have passed over the point x ∈ S. Thus, the quantized phase field ξ describes an ensemble of perfect, or Volterra, dislocations on S. In writing (30) and (35), X × Y denotes the space of slip distributions and integer-valued phase fields of finite energy, i. e., the linear space of functions ζ : R 2 → R and ξ : R 2 → Z such that E[ζ|ξ] < ∞, which is the physically relevant space of solutions.
Next we seek to express the energy of the dislocation ensemble directly in terms of ξ. To this end, we begin by inserting definition (30) into the variational problem (35), with the result: At this point, we may invert the order of minimization in (36), which results in the reduced minimum problem: where This minimization leads to a linear problem in ζ which may be solved analytically by an application of the Fourier transform. To this end, we apply Parseval's identity to obtain: Minimization with respect to ζ gives:ζ is the slip distribution for ξ = 0, i. e., the slip distribution due to the elasticity of the interatomic plane. Insertion of (40) into (38) gives the sought energy: where is independent of ξ and represents the elastic energy of the interatomic in the absence of dislocations, i. e., at ξ = 0. Thus, for instance, if s = τ = constant, (43) corresponds to an elastic energy density per unit area of (d/2µ)τ 2 , which indeed coincides with the elastic energy density of the interatomic plane. Therefore, the excess energy may be identified with the energy proper of the dislocation ensemble. By virtue of (42), this energy is now expressed entirely in terms of the phase field ξ. If is evident from the form of (42) that, in the absence of an applied field, the dislocation energy vanishes for uniform slip, as required by symmetry. Indeed, in this casê ξ = nδ D (k), where n ∈ Z and δ D (k) denotes the Dirac delta at k = 0, and, consequently, Kξ vanishes identically in the sense of distributions. As already noted, the phase field ξ in (42) takes integer values and, therefore, represents a distribution of Volterra or perfect dislocations. It is evident from the form of (42) that the Peierls potential has the effect of regularizing linear elasticity so as to render the energy of Volterra dislocations finite. The lattice parameter is retained in the energy functional (42) through the interplanar distance d. The effect of the core regularization resides in the factorφ appearing in both terms of the energy (42). This factor tends to suppress wavelengths on the scale of the lattice parameter or shorter. The same factor relates the integer-valued phase field ξ to the core-regularized slip distribution ζ, eq. (40). In this relation, the function ϕ d (x) acts as a mollifier: the integral of ϕ d (x) over the entire plane is one, or, equivalently,φ d (0) = 1; and ϕ d (x) defines a Dirac-delta sequence, i. e., in the sense of distributions. In real space eq. (40) takes the convolution form: which shows that the effect of the piecewise-quadratic Peierls interplanar potential is to smooth out the phase field ξ. A consequence of this smoothing is that the perfect dislocations described by ξ acquire a core of width commensurate with d.
It should be carefully noted that, despite the quadratic appearance of the energy functional (42), the attendant variational problem is strongly nonlinear owing to the all-important constraint that ξ be an integervalued function. Furthermore, the minimization problem is nonconvex owing to the lack of convexity of the set X. This nonlinear and nonconvex structure of the variational problem is inherited from the-similarly nonlinear and nonconvex-structure of the Peierls potential (28). In addition, the variational problem is nonlocal, owing to the presence of long-range elastic interactions. These attributes render the problem (37) mathematically non-trivial and confer its solutions a rich structure. Despite these difficulties, the choice of a piecewise quadratic Peierls potential lends the problem analytical tractability, as demonstrated subsequently.

Illustrative examples
As an illustration of the type of core structure predicted by the piecewise quadratic model, we may consider the case of a general straight dislocation. For simplicity, we restrict our attention to the isotropic case. It is convenient to introduce auxiliary orthonormal axes (x ′ 1 , x ′ 2 ) with x ′ 1 normal to the dislocation line. In this coordinate frame one has:ξ corresponding to uniform slip over the half-plane x ′ 1 < 0. Here, δ D denotes the Dirac-delta distribution, and i is the unit imaginary number. A straightforward calculation gives the slip distribution aŝ where is a characteristic core width, and θ is the angle subtended by the normal to the dislocation line and the Burgers vector. Thus, the case θ = 0 corresponds to an edge dislocation, whereas the case θ = π/2 corresponds to a screw dislocation. A plot of the core profile is shown in Fig. 2. In addition, the energy per unit dislocation length follows in the form where R is an upper cut-off radius. It is clear from this expression that c plays the role of a core cut-off radius. The width of the core depends on the orientation of the dislocation line relative to the Burgers vector, and it attains its minimum (maximum) value for edge (screw) dislocations. This example illustrates how the introduction of a Peierls interplanar potential renders the dislocation core structure, including the core width, well-defined. In particular, the high wave number divergence in the dislocation energy is eliminated. The specific core structure (49) is, of course, a result of the assumed piecewise quadratic form of the Peierls interplanar potential. It bears emphasis that here the interplanar potential is regarded simply as a convenient device for regularizing the equations of elasticity, and no attempt is made to model real dislocation core structures. However, it seems reasonable to expect that, for sufficiently well-spaced dislocations, the only macroscopically relevant core parameter is the core energy per unit length, and that the details of the core structure play a limited role as regards the overall energetics of the dislocation ensemble. The Peierls potential introduces a lengthscale d into the formulation of the order of the crystal lattice parameter. The behavior of dislocation loops may therefore be expected to differ sharply according as to whether their size is much larger or, contrariwise, much smaller, than d. These regimes are exhibited in the simple example of a circular dislocation loop of radius R. In this case, ξ is the characteristic function of the circle of radius R, e. g., centered at the origin, and its Fourier transform iŝ where J 1 is the Bessel function of the first kind, and we have written (x 1 , x 2 ) = (r cos θ, r sin θ) and (k 1 , k 2 ) = (k cos φ, k sin φ). Insertion of (52) into (42) gives (25), and we have set ν = 0 for simplicity. From the asymptotic behavior of the Meijer function it follows that As expected, for fixed d the energy of the loop tends to zero as R → 0. Conversely, for fixed R, the energy diverges as d → 0, which corresponds to the limiting case of a perfect Volterra loop in a linear elastic crystal. Once again, this divergence illustrates the breakdown of linear elasticity in the presence of perfect dislocations and the crucial role played by the core regularization in eliminating that breakdown. The form (55) for the energy of large loops is consistent with expressions derived from linear elasticity (24). However, it should be noted that the core width c = d/4 implied by (55) is predicted by the present theory, which stands in contrast to the theory of linear-elastic dislocations wherein the core cut-off radius is an ad-hoc and extraneous parameter. The issue of metastability and the existence of relative energy minimizers may be illustrated simply by the stable configurations of a dislocation dipole under the action of a resolved shear stress τ . In order to preclude uniform slip distributions, ξ may be required to decay to zero at infinity. Suppose that slip occurs on the interval −r/2 < x 1 < r/2, where r is dipole width. Then, the corresponding phase field iŝ For purposes of the present discussion it suffices to consider the asymptotic limit of r ≫ d. In this regime (42) yields the relation and From these two expressions we find: for some constant E 0 independent of r and c. Equilibrium now demands: The dipole energy at equilibrium is, therefore, where r is the solution of (60). The preceding analysis shows that, in the regime under consideration, there exist stable dipoles of arbitrary width r, with Eq. (60) supplying the resolved shear stress which equilibrates the dipole. It is evident from (61) that, for sufficiently large r, the dipole energy is positive and, therefore, greater than the dislocation energy corresponding to a phase field ξ = 0, which is zero. Hence, the dipole configuration is metastable and constitutes a relative minimum of the dislocation energy.

Energy-minimizing phase fields
The variational problem (37) which characterizes the stable equilibrium configurations of the slip plane may be regarded as a constrained minimization problem, in which the unknown field ξ is constrained to take integer values everywhere on the slip plane. This constraint is nonlinear and nonconvex, which renders the problem (37) mathematically nontrivial. In addition, the long-range elastic interactions between dislocations render the energy nonlocal. These difficulties notwithstanding, in this section we obtain a distinguished set of minimizers of the energy (42) analytically in closed form.

Constrained minimizers and energy-norm projection
A compelling characterization of the solutions of the constrained minmization problem (37) may be derived as follows. Let η be the unconstrained energy minimizer. Thus, η is the solution of the unconstrained minimization problem: In particular, η is allowed to take arbitrary real values. Then, the constrained phase field is the solution of the problem (e. g., (26)): is an energy norm corresponding to the inner product The solution of (63) may be expressed as where P X denotes the closest-point projection of Y onto X in the sense of norm (64). The solution procedure leading to the core-regularized slip distribution ζ may now be broken down into three steps. The first step consists of the solution of the unconstrained minimization problem (62), resulting in the real-valued unconstrained minimizer η. The second step consists of the solution of the constrained minimization problem (63), resulting in the projection of η onto its closest integer-valued phase-field ξ. The third step consists of the smoothing of ξ according to (40), leading to the introduction of smooth cores around all dislocation lines.
The exact evaluation of the projection (63) entails some difficulty. The corresponding stationarity condition consists of the requirement that the functional be stable with respect to a unit increment or decrement of ξ over an arbitrary region of the slip plane. Thus, if R ⊂ R 2 is one such region and χ R denotes its characteristic function, stationarity demands that the inequalities be simultaneously satisfied. After some straightforward manipulations (67-68) may be recast in the form These inequalities set bounds on the possible values of the resultant of the residual tractions acting over arbitrary regions of the slip plane. We proceed to proof the remarkable result that ξ(x) is the closest integer to η(x) pointwise on the slip plane, i. e., where P Z denotes the closest-point projection of R onto Z, i. e., ∀x ∈ R, P Z x is the integer closest to x. In order to verify the bounds (69), we choose circles of radius ǫ → 0 as the test regions R. Let B ǫ denote the circle of radius ǫ centered at the origin, and set where x ranges over all of R 2 . Given a function u : In order to elucidate the structure of this inner product as ǫ → 0, we may introduce the normalized variables: whereupon (72) becomes In the limit of ǫ → 0 we have, asymptotically, or, in terms of the original set of variables, In addition, distributionally as ǫ → 0 we have χ Bǫ ∼ πǫ 2 δ D , and (76) simplifies to This expression finally evaluates to whereas choosing u = χ R , with χ R as in (71), yields in agreement with (54). Inserting (79) and (80) into (69) finally gives the stationarity condition This result is remarkable in several notable respects. Firstly, given the nonlocal character of the energy, the local character of the projection (70) is unexpected. Secondly, as we have already noted, problem (38) is non-convex owing to the non-convexity of the admissible set X. It is therefore remarkable that the problem is amenable to an exact, and particularly simple, solution. Finally, it is interesting to note again the crucial role played by the core regularization of the energy which, provided that d > 0, results in finite energies (79) and (80).

Unconstrained equilibrium equations
Next we turn to the unconstrained minimization problem (62). The corresponding Euler-Lagrange equation is: It is clear from this equation that all functions of the formη ∝ δ D , corresponding to uniform slip distributions over the entire slip plane, are solutions of the homogeneous equation. Therefore, it follows from the Fredholm alternative theorem that a necessary condition for (82) to have solutions is thatŝ(0) = 0, i. e., This condition expresses the requirement that the resultant of the resolved shear-stress field be zero. If this condition is met, thenη follows from (82) aŝ where C is an arbitrary constant. An application of the inverse Fourier transform to this expression gives where the fundamental solution represents an unconstrained slip distribution which decays to zero at infinity and is in equilibrium with a resolved shear stress in the form of a Dirac delta applied at the origin.

The general solution
We proceed to collect all the solution steps outlined in the foregoing and to provide an explicit expression for the slip distribution ζ as a function of the applied field s. Thus, the solution ζ follows in three steps, namely: s → η → ξ → ζ. The first step yields the unconstrained phase field η, the second the integer-valued phase field ξ, and the third step returns the core-regularized slip distribution ζ. Gathering the relations (85), (66) and (40) corresponding to each of the three steps just enunciated gives, explicitly, where the projection P X is given by (70) and, as noted earlier, we require s to have zero mean, whereupon the solution is determined up to an arbitrary integer C ∈ Z. At this point, it is illuminating to revisit the question of uniqueness and metastability, especially in view of the existence of relative energy minimizers such as presented in Section 2.4. In this respect, it should be noted that the solution procedure presented in the foregoing returns a unique solution, modulo uniform slips. Indeed, the unconstrained equilibrium equations are linear and, within the constraints imposed by the Fredholm alternative, determine a unique unconstrained slip distribution η. The local truncation (70) in turn yields a unique integer-valued phase field ξ. Finally, smoothing as in (47) results in a unique coreregularized phase field ζ. It is therefore clear that the solution procedure outlined above selects a 'preferred' slip distribution among a vast array of competitors. Thus, for instance, in regions of the slip plane where the dislocation density is small, i. e., where the distance between neighboring dislocations is large compared to the core size d, metastable slip distributions may be obtained by inserting small loops in equilibrium with the local resolved shear stress within the intervening area between dislocations. Examples of metastable solutions of this type have been discussed in Section 2.4. It is intriguing that, in those examples, e. g., in the case of a plane slipping under the action of a uniform resolved shear stress, the unique preferred solution ζ = 0 delivers the minimum attainable energy, i. e., it is an absolute energy minimizer. An open mathematical question of some interest is whether the preferred solution (87) does indeed always deliver an absolute energy minimizer.

Averages and macroscopic variables
Many of the macroscopic quantities of interest pertaining to the behavior of a slip plane may be recovered by taking the appropriate averages. In order to render the operation of taking averages over the slip plane well defined, we may simply assume that the phase field is periodic with unit cell Ω, i. e., for some basis vectors {a 1 , a 2 }. Under these conditions, an application of Orowan's relation gives the macroscopic slip strain as or, in view of (40), alternatively as where is a reference slip strain. Here l is the slip-plane spacing, which is assumed known, and γ 0 is the slip strain which is attained when one loop sweeps over the entire slip plane. The dislocation line density per unit volume is related to the phase field as Indeed, the integral of |∇ξ| over a region of the slip plane simply measures the total length of all dislocation lines contained in that region. Finally, the global equilibrium condition (83) may equivalently be expressed in the form: which requires that the mean resolved shear stress be zero.

Irreversible processes and kinetics
The preceding developments have focused on the energetics of a dislocation ensemble in an elastic crystal endowed with a piecewise quadratic Peierls interplanar potential. In addition to this energetics, the dislocation ensemble may undergo inelastic interactions with the lattice, resulting in lattice friction; and short-range inelastic interactions with an assortment of obstacles, such as second-phase particles, forest dislocations, and others. We assume that these interactions are irreversible and, therefore, kinetic in nature. The essential assumption is that the crossing of an obstacle by a dislocation 'costs' a certain energy, regardless of the direction of crossing. The energy toll depends on the strength of the interaction and is dissipated, e. g., as heat, and irreversibly lost to the system.

Variational formulation
The irreversible dislocation-obstacle interactions may be built into the variational framework developed previously by recourse to time discretization (5; 4; 6; 7; 8). Thus, henceforth we consider a sequence of discrete times t 0 , t 1 , . . . , t n , t n+1 , . . . , presume the slip distribution ζ 0 at time t 0 to be known, and seek to compute the slip distributions ζ 1 , . . . , ζ n , ζ n+1 , . . . , at all subsequent times. The central problem is to determine the slip distribution ζ n+1 at time t n+1 given the solution ζ n at time t n and the applied field s n+1 at time t n+1 . To this end, following (5; 4; 6; 7; 8) we introduce the incremental work function: where the elastic energy E[ζ] is given by (30), and the second term represents the incremental work of dissipation. In this term, the field f (x) ≥ 0 represents the energy cost per unit area associated with the passage of one dislocation over the point x, i. e., with a transition of the form ζ(x) → ζ(x) ± 1. Thus, the field f (x) represents the distribution of obstacles over the slip plane and is assumed known. By the work and energy identity, it follows that W equals the total work supplied to the system during the interval [t n , t n+1 ].
Eq. (94) simply states that part of this work is invested in raising the energy of the crystal, whereas the remainder of the external work supplied is invested in overcoming the obstacle resistance. The updated slip distribution now follows from the minimum principle: It should be carefully noted that the work function (94) depends on the initial conditions ζ n for the time step, which allows for irreversibility and hysteresis.

Min-max formulation
In order to enable the application of the general solution procedure outlined in the foregoing, we proceed to rephrase problem (95) as a min-max problem (e. g., (27), § 36). We begin by noting the identity: Since f (x) ≥ 0, it follows from this identity that which may be viewed as a statement of the principle of maximum dissipation. The field g n+1 (x), which is required to lie within the bounds ±f (x), may be regarded as the reaction-force field exerted by the obstacles on the dislocation ensemble at time t n+1 . Inserting representation (97) into (95) and using standard properties of saddle points of concave-convex functions (e. g., (27), Theorem 36.3) leads to the problem For fixed g n+1 (x) the central minimization problem stated in (98) is now of the form (35) analyzed earlier, and the solution is given by (87). The problem which then remains is to optimize the field g n+1 (x) in order to maximize the dissipation at the obstacles. The optimization of g n+1 (x) may conveniently be reduced to a problem with linear constraints by the introduction of Lagrange multipliers (e. g., (27), § 28), which leads to the Lagrangian: The corresponding Kuhn-Tucker optimality conditions are (e. g., (27), Theorem 28.3): which must be satisfied simultaneously by the solution. In (100), the updated slip distribution ζ n+1 follows from the general solution (87) as and the global equilibrium condition (83) now becomes which requires that the total applied force be exactly equilibrated by the resultant of the obstacle reactions. It is interesting to note how the optimality conditions (100-103) give mathematical expression to the assumed frictional interaction between the phase and obstacle fields, including stick-slip behavior. Thus, the equalities (103) prevent λ + (x) and λ − (x) from being nonzero simultaneously. If, for instance, ζ n+1 (x) − ζ n (x) > 0, then by virtue of (102) and (100) it necessarily follows that λ − (x) = 0 and λ + (x) > 0, and, by the first of (103), g n+1 (x) = f (x). Conversely, if ζ n+1 (x) − ζ n (x) < 0, then (102) and (100) require that λ − (x) > 0 and λ + (x) = 0, and, by the second of (103), g n+1 (x) = −f (x). These two cases correspond to slip conditions at x, or yielding. If the conditions (101) are satisfied as strict inequalities, then it follows from (101) that, necessarily, λ + (x) = 0 and λ − (x) = 0, which, in view of (100), requires that ζ n+1 (x) − ζ n (x) = 0. This case corresponds to stick conditions, or elastic unloading. The theory is thus capable of describing the loading-unloading irreversibility characteristic of plastic materials.

Special case of short-range obstacles
Next we consider the special case of an obstacle field consisting of a uniform Peierls stress and a distribution of short-range obstacles. This case merits special attention owing to the fact that the problem can conveniently be reduced to the determination of the phase field at the obstacles, which greatly facilitates calculations.
For definiteness, we shall assume the resolved shear-stress field to be of the form where τ (t) is the applied resolved shear stress at time t and s 0 (x) is a self-equilibrated field representing the long-range elastic stresses induced by the obstacles. For instance, in calculations of forest hardening s 0 (x) may describe the action on the slip plane of the long-range stress field of the forest dislocations. In this case, in the vicinity of a secondary dislocation s 0 (x) tends asymptotically to the field of the osculating straight dislocation, and, therefore, has the form for some function A(θ) which depends on the geometry of the secondary dislocation. Here (r, θ) are polar coordinates centered at the secondary dislocation. The corresponding unconstrained phase field η 0 diverges logarithmically, i. e., for some function B(θ). The attendant phase field ξ 0 represents an infinite pile up of perfect dislocation loops encircling the secondary dislocation. The depth of this pile up is truncated and rendered finite when the core of the secondary dislocations is taken into account. Since the exact geometry of the secondary dislocations away from the primary slip plane is unknown within the present formulation, one possibility that immediately suggests itself is to model s 0 as a random field. The use of random fields to model barriers to dislocation slip, or other similar processes such as first-order phase transitions in magnetic systems, has been proposed by Hardikar et al. (28) and Sethna (29). For simplicity, we shall additionally assume that the obstacle distribution and, correspondingly, the phase field, is periodic with unit cell Ω. For the obstacle system under consideration, the obstacle-strength field may be taken to be of the form: where τ P ≥ 0 is the Peierls stress, x i and f i ≥ 0, i = 1, . . . , N , are the positions and strengths of the obstacles in Ω, respectively, and the function ψ d (x) represents the structure of the obstacles. As noted earlier, the detailed modeling of the structure of obstacles is beyond the scope of the present work. However, the treatment of short-range obstacles as points, corresponding to setting ψ d = δ D in (109), inevitably leads to logarithmic divergences of the phase field under the obstacles. In order to avoid these divergences the obstacles must be endowed with a finite core. Conveniently, however, for distributions of well-separated short-range obstacles the precise form of ψ d plays a role only in the immediate vicinity of the obstacles and is otherwise largely irrelevant. We shall simply require that the overlap function ψ d ⋆ ϕ d define a Dirac-delta sequence as d → 0, i. e., that in the sense of distributions. Thus, for points such that |x − x i | ≫ d the overlap function ψ d ⋆ ϕ d may be treated as a Dirac-delta. In addition, in order to avoid divergent integrals we shall require that the productψ d (k)φ d (k) decay as 1/k 2 as k → ∞. For forest obstacles, a simple choice consistent with these requirements is which simply states that the core structure of the forest obstacles is identical to the core structure of the primary dislocations. The problem to be solved now follows by inserting (109) into the Lagrangian (99). In order to facilitate the solution of this problem we make the ansatz that the obstacle reaction field is of the form: where the energy functional E[ξ] is given by (42). The Kuhn-Tucker optimality conditions (100-103) now reduce to: with i = 1, . . . , N . If γ n+1 = γ n , i. e., in the presence of macroscopic slip, eqs. (114 -117) simply require that where sgn(x) = x/|x| is the signum function. In addition, if |g n+1 0 | < b τ P , then it necessarily follows from (114-117) that γ n+1 = γ n , corresponding to elastic unloading. Thus, the Peierls stress has the effect of introducing an initial threshold, or yield point, for plastic activity.
In evaluating the optimality conditions (118 -121), the phase field at the obstacles is computed as which is a special case of the pointwise projection (70). In addition, the unconstrained phase field follows by inserting (112) into (85) and taking periodicity into account, with the result Specializing this expression at obstacle sites we obtain the linear system of equations which directly relates the unconstrained phase field at the obstacles and the obstacle reactions. In addition, the global equilibrium condition (105) further reduces to where we have made use of the identity s 0 = 0. Eq. (126) simply states that the applied resolved shear stress must be exactly equilibrated by lattice friction and the resultant of all obstacle reactions. In view of (110), the influence coefficients coupling pairs of distinct obstacles in (125) may be computed as In arriving at this expression we have made use of the assumption that the short-range obstacles are wellseparated, so that they effectively interact as point obstacles. The evaluation of the self-interaction coefficients G ii in (127) requires some care, as the obstacles can no longer be treated as points and their assumed structure must be taking into consideration. In this case we have where the introduction of the function ψ d describing the obstacle core is essential in order to avoid a divergent value of G ii . If, by way of example, we choose the obstacle structure (111), then a straightforward calculation gives: In general, this constant is well defined and scales as 1/bd provided that, as previously assumed, the product

Algorithmic implementation
The preceding relations provide a complete basis for updating the phase field incrementally. The precise algorithm employed in the calculations reported subsequently proceeds as follows: 1. Assembly and factorization. The calculations are started by computing the matrix G ij using formulae (127), (128) and (129). This matrix is symmetric and positive-definite and has dimension N × N . The matrix G ij is then factorized. It should be noted that, for a constant distribution of obstacles, the matrix G ij remains constant throughout the entire deformation process and its factorization may be performed once and for all at the start of the calculations.
2. Initialization. At time t 0 , the constant C 0 in (125) is set to zero. In addition, the unconstrained phase field η 0 induced by the residual tractions s 0 is set up by superposition of fields of the form (108), suitably truncated at the origin so as to account for the obstacle core. The sign of the superposed fields is alternated randomly in such a way that the mean value of s 0 is zero, as required.
3. Incremental update. Assume that the state of the slip plane is completely known at time t n , and that a new value C n+1 of the constant in (125) is prescribed at time t n+1 . Thus, the constant C is chosen as a convenient control or loading parameter. The state of the slip plane is updated by means of the following operations: (a) Stick predictor. Set the predictor unconstrained phase fieldη n+1 i = η n i , and compute the predictor reactions:g 4. Phase-field evaluation. Once the obstacle reactions g n+1 i are known, the unconstrained phase field η n+1 (x) may be computed from (124). At the obstacles, the unconstrained phase field follows directly as For x sufficiently distant from all obstacles, we may treat the obstacles as points, whereupon (124) simplifies to: Finally, the local projection (70) returns the phase field ξ n+1 (x), and (47) gives the core-regularized phase field ζ n+1 (x).

Application to the forest hardening mechanism
In this section, we apply the general framework developed in the foregoing to the forest hardening mechanism. To this end, we consider a slip plane traversed at random locations by fixed secondary or forest dislocations and acted upon by an applied resolved shear stress. The objective is to characterize the dislocation patterns which arise in response to the applied loading, and to determine the effective behavior of the system measured, e. g., in terms of the macroscopic slip strain and dislocation line density. In situations where the forest hardening mechanism is dominant, the kinetics of the primary dislocation ensemble is primarily governed by: the interaction between the dislocations and the applied resolved shear stress; the long-range elastic interactions between the primary dislocations; the line-tension effect resulting from the core structure of the dislocations; the interaction between the primary dislocations and the long-range elastic stress field of the forest dislocations; and the short-range interactions between primary and forest dislocations such as may result in jogs, the formation of junctions, and other reaction products. Thus, the forest hardening mechanism provides a convenient framework for illustrating the range of behaviors predicted by the theory.

Monotonic Loading
We begin by considering the base case of a single slip system containing a fixed concentration of point obstacles of uniform strength and deforming under the action of a monotonically increasing shear stress. We assume a periodicity, and the unit cell Ω is taken to be a square of dimension 100b. We randomly select 100 points within the periodic cell as the obstacle sites. To these obstacles we assign a uniform strengths f = 10µb 2 . For simplicity, we take the Peierls stress τ P = 0, and set Poisson's ratio ν = 0.3.  slip strain. Fig. 4 shows the corresponding evolution of the phase field ξ and the dislocation pattern. In this latter plot, the dislocation lines are identified with the lines across which the phase field ξ jumps by one, or equivalently, with the level contours of η + 1/2 at integer heights. The phase field itself counts the number of dislocation which have passed over a given point of the slip plane, with proper accounting of the sign of traversal. Two well-differentiated regimes, with a transition point at γ/γ 0 ∼ 10, are clearly discernible in Fig. 3: a first regime of 'micro-slip' dominated by dislocation bow-out; and a second regime characterized by generalized yielding. The saturation value of the resolved shear stress in this latter regime follows from global equilibrium as i. e., the saturation stress is set by the mean obstacle strength. During the micro-slip regime, the value of the phase field remains close to zero at the obstacles, which are essentially impenetrable to the dislocations. The dislocation loops move reversibly by bowing through the open spaces between obstacles, and by subsequently pinching behind the obstacles, leading to the formation of Orowan loops. As increasing numbers of dislocations bypass the obstacles by the bow-out mechanism, the number of Orowan loops surrounding the obstacles correspondingly increases, Fig. 4. The details of this process are illustrated in Fig. 5, which zooms on a particular pair of obstacles and shows the evolution of the dislocation/obstacle interactions with increasing stress. In this figure, individual dislocations move from left to right and are color coded for ease of tracking. Fig. 5a shows the initial configuration. In Fig. 5b, the dislocation has bypassed the obstacles and left Orowan loops in its trail. In subsequent frames, new dislocations arrive from the left, bow through the obstacles, and eventually bypass them, leaving behind additional Orowan loops. The ease with which the phase-field representation describes these geometrical and topological transitions is quite remarkable.
At larger applied stresses, the obstacles gradually yield and are overcome by the dislocations. Figs. 4de correspond to the transitional phase between the micro-slip and the generalized yielding regimes. At sufficiently large applied stresses close to the saturation stress τ 0 , all obstacles yield and are crossed by the dislocations. Under these conditions, the phase field ξ increases uniformly without change of shape over the entire slip plane. In this regime, the level contours of ξ remain ostensibly unchanged, even as their height rises steadily, and the dislocation pattern becomes frozen in place, Fig. 4e and 4f. The hardening rate correspondingly drops and the stress-strain curve saturates asymptotically, Fig. 3a. The presence of a saturation stage is consistent with observations of fcc single crystals undergoing single slip, e. g., during the easy glide or stage I or hardening.
The evolution of the dislocation line density (92) with slip strain is shown in Fig. 3b. The computed dislocation densities are normalized by the reference density which is the limiting or saturation density corresponding to an arrangement of parallel straight dislocations at intervals of b. As may be seen from the figure, the slip plane contains a nonzero dislocation density at zero slip strain. This initial density is induced by the long-range elastic-stress field s 0 of the secondary dislocations. The process of slip is accompanied by a steady increase in the dislocation line density (92). Dislocation multiplication and the proliferation of dislocation loops with increasing slip strain are clearly evident in Figs. 4a-f. In the micro-slip regime, the increase in macroscopic slip is accompanied by a steady supply of dislocations bowing between the obstacles. This results in an initial parabolic growth rate, ρ ∼ γ 2 , Fig. 3b. By contrast, as already noted when τ approaches the saturation stress τ 0 the dislocation pattern becomes frozen, e. g., Figs. 4e and f, and the dislocation density remains ostensibly constant. This parabolic growth regime and the subsequent saturation phase are indeed consistent with observation (30; 31). The ability of the theory to predict the evolution of the dislocation density, in addition to predicting hardening rates, is noteworthy. Fig. 6 shows the effect of a variable obstacle strength on the computed stress-strain curve and dislocation density evolution during monotonic loading. In order to investigate this effect we consider obstacle strengths of the form f = 10µb 2 (1 + ω), with the random variable ω distributed uniformly within the intervals [−0.2, 0.2] and [−0.5, 0.5]. When the obstacle strengths exhibit variability, the weakest obstacles tend to yield first and are overcome by dislocations, which are pinned at the stronger obstacles. Consequently, the transition between the micro-slip and generalized yielding regimes becomes more gradual. It is interesting to note that, since the average obstacle strength is the same in all cases, the saturation stress is not affected by the obstacle strength variability. This example illustrates the ability of the theory to account for the combined effect of obstacles of different species.

Cyclic Loading
Next we investigate the behavior predicted by the theory under load reversal and cyclic loading. As in the preceding example, we assume a periodicity and take the unit cell Ω to be a square of dimension 100b. We randomly select 100 points within the periodic cell as the obstacle sites and assign to these obstacles a uniform strengths f = 10µb 2 . The Poisson's ratio ν is set to 0.3, and, for simplicity, we take the Peierls stress τ P = 0. We consider the cases of partial unloading, where the system is reloaded after a small amount of reverse yielding; and a fully-reversed loading cycle. In both cases, the unloading and reloading events bring to the fore the irreversible and dissipative nature of the system, and necessitate full use of the incremental theory developed in Section 4. Figs. 7a and b show the predicted stress-strain curve and the evolution of the dislocation density in the cases of fully-reversed loading. The corresponding evolution of the dislocation ensemble is shown in Fig. 8 for the case of fully-reversed loading. The various matching sections of these curves are labeled for ease of reference. The interval O-a of the curves in Fig. 7 simply reproduces the monotonic response described in the preceding section. Upon load reversal at point a, the system initially unloads elastically, Fig. 7a, segment a-c. In this stage, the applied shear stress continues to push the dislocations against the obstacles in the loading direction. However, as the applied load is decreased the force exerted by the dislocations on the obstacles correspondingly decreases and, consequently, remains below the obstacle strength. Under these conditions, the dislocations remain pinned at the obstacles and simply recoil elastically. This process gives the linear unloading branch a-c in the stress-strain curve, Fig. 7a.
It should be carefully noted that, owing to the frictional sliding of the dislocations over the obstacles that occurs during the loading state, the interaction forces between the dislocations and the obstacles do not vanish when the applied shear stress is reduced to zero. Instead, the system is left in a state of permanent or residual deformation in which the dislocation/obstacle force system is self-equilibrated. The corresponding dislocation pattern for the fully-reversed case is shown in Fig. 8c. This remanent field is in analogy to the self-equilibrated residual stress fields which remain in elastic-plastic solids upon unloading. The combination of a reversed applied stress and the residual force field causes the system to yield prematurely in the reserve direction, Fig. 7a, point d, and thus the theory predicts the Bauschinger effect.
Once reverse yielding commences, the obstacles continue to oppose the motion of the dislocations, which now takes place in the reverse loading direction. In particular, the frictional forces at the obstacles switch sign relative to the loading phase. This process of switching is clearly apparent in Fig. 9, which shows the evolution of the phase field during the loading cycle. Thus, during the loading phase the obstacles pull down on the phase field, Fig. 9b-d, causing it to cusp downward at the obstacles, whereas during the reverse loading phase the obstacles pull up on the phase field, Fig. 9g-i, causing it to cusp upward. The stress-strain curve shown in Fig. 10 reveals that the system exhibits 'fading memory'. Thus, when the system is reloaded, the stress-strain curve gradually transitions towards the virgin loading curve, and, with sufficient reloading, the system eventually 'forgets' the unloading cycle. The extent of this loss of memory depends on the extent of reverse yielding. Consider, for instance, the case in which unloading is purely elastic and no reverse yielding occurs. This corresponds to reloading, from point b in Fig. 10. Because no reverse sliding over obstacles has occurred at this point, the dislocation/obstacle configuration remains undisturbed. Consequently, upon reloading all the obstacles yield simultaneously and the system exhibits an abrupt yield point coinciding exactly with the point of unloading, point a in Fig. 10. Thus, in this case no loss of memory occurs and the system exhibits return-point memory in the sense of Sethna et al (29). This type of behavior should now be compared with that corresponding to a limited amount of reverse loading. In this case, the dislocations slip over the obstacles in the reverse loading direction, and some of the details of the dislocation/obstacle configuration established at the unloading point d are lost in the process. When the system is reloaded, point e in Fig. 10, the obstacles yield anew gradually and not all at once, as in the case of a purely elastic unloading, and the system exhibits a gradual transition towards the virgin loading curve, Fig. 10, segment e-f. In the case of fully-reversed unloading, Fig. 7, the amount of reverse yielding, segment d-e, is large enough that the unloading point is effectively wiped out from the memory of the system, with the result that the reloading curve e-i is ostensibly identical to the unloading curve a-e. These trends are in good agreement with the experimental cyclic stress-strain data for structural steels reported by Ortiz and Popov (32), which was obtained from tests specially designed to exhibit the fading memory effect just described.
The evolution of the dislocation density during a loading cycle is of considerable interest as well, Fig. 7b . Thus, upon unloading the dislocation density decreases as a result of the elastic relaxation of the dislocation lines. The dislocation density bottoms out -but does not vanish entirely -upon the removal of the applied stress, Fig. 7     state, Fig. 8c. The dislocation density increases again during reverse loading, Fig. 7 b, segment c-d, and the cycle is repeated during reloading, Fig. 7b, segment c-d, segment e-i, giving rise to a dislocation density vs. slip strain curve in the form of a 'butterfly'. This type of behavior is indeed observed experimentally (33). It also arises in models of the stored energy of cold work (34), and is in analogy to the hysteretic loops exhibited by magnetic systems (35; 29; 36). Finally, we exercise the model over the range of periodic cell sizes 10 2 b − 10 5 b, while at the same time keeping the number of obstacles constant at 100. The corresponding obstacles density c = N/|Ω| consequently varies from 10 −2 /b 2 to 10 −8 /b 2 . Of primary interest here is to ascertain how the macroscopic behavior of the slip plane depends on the obstacle density. Our calculations show that the stress-strain curve, when expressed in terms of the normalized variables τ /τ 0 and γ/γ 0 is ostensibly independent of c. It should be carefully noted that, for obstacles of uniform strength, τ 0 = f c, and thus the saturation strength scales in direct proportion to the obstacle density. The variation of the dislocation density evolution with c is shown in Fig. 11. It is interesting to note from this figure that the saturation value of the dislocation density scales as 1/ √ c, but that otherwise the evolution of the dislocation density with macroscopic slip strain remains ostensibly identical in all cases. This example also serves to underscore the ability of the theory to effectively deal with large domains of analysis. Thus, owing to the absence of a computational grid the complexity of the calculations scales with the number of obstacles and is otherwise independent of the cell size, which enables the consideration of large cells sizes in the relevant range of observation.

Dislocation line-energy anisotropy
We conclude this section with a brief parametric study of the effect of dislocation line-energy anisotropy on the hardening characteristics of the system. In particular, we consider systems which are identical in every respect except for the value of Poisson's ratio. By virtue of eq. (51), different choices of Poisson's ratio result in different ratios between the energy per unit length of edge and screw segments. We recall that, for ν > 0, the minimum (maximum) energy per unit length is attained for pure screw (edge) segments. All remaining parameters are as in the preceding examples. In particular, the unit cell dimension is 100b, the number of obstacles is 100, the obstacle strength is f = 10µb 2 , Poisson's ratio ν is 0.3, and the Peierls stress τ P is set to 0. Fig. 12 shows the predicted stress-strain curves and dislocation density evolution for ν = 0.0, 0.3 and 0.5. Since the saturation stress τ 0 solely depends on obstacle strength and concentration, all three stressstrain curves saturate at the same level. By contrast, as ν is increased the dislocation line-energy of both edge and screw segments increases and the system correspondingly stiffens, with the result that lower values of the slip strain are attained, Fig. 12a. The rate of dislocation multiplication also decreases with increasing ν, Fig. 12b. Perhaps more interestingly, Fig. 13 shows the effect of the anisotropy on the dislocation geometry. Thus, as Poisson's ratio is increased screw segments become energetically favorable relative to edge segments, and the dislocation pattern exhibits a preponderance of elongated screw segments, Fig. 13c.

Summary and Concluding remarks
We have developed an analytically tractable phase-field theory of dislocation dynamics, strain hardening and hysteresis in ductile single crystals at low temperatures. The phase-field representation furnishes a simple and effective means of tracking the motion of large numbers of dislocations within discrete slip planes through random arrays of point obstacles under the action of an applied shear stress. The theory rests on a variational framework for dissipative systems and accounts for energetic and kinetic effects (4; 5; 6; 7; 8). The energetics accounted for in the theory include: the core energy of the dislocations, represented by a piecewise-quadratic Peierls potential (9); the long-range elastic interactions between primary dislocations and between primary and forest dislocations; and the energy of interaction with the applied resolved shear stress. The kinetics of the system stem from the assumed irreversible interactions between dislocations and obstacles, and from lattice friction, and result in hardening, path dependency, and hysteresis in the macroscopic behavior.
The theory predicts a range of behaviors which are in qualitative agreement with observation. Thus, in the base case of single slip under monotonic loading, the theory predicts saturation at a stress which depends on obstacle strength and density. It also predicts dislocation multiplication at an initial parabolic rate, during a 'micro-slip' regime, followed by saturation at a maximum dislocation density. The ease with which the phase-field representation allows for complex geometrical and topological transitions, as the dislocation ensemble percolates through the obstacles, is quite remarkable. The cyclic behavior predicted by the theory is particularly noteworthy. As required, the theory predicts the essential phenomena of loading/unloading irreversibility and hysteresis. But the theory also captures more subtle aspects of the cyclic behavior of crystals, such as the Bauschinger effect, consisting of the premature yielding of a slip system under reverse loading; the fading memory effect, whereby reverse yielding gradually eliminates the influence of previous loading; and the evolution of the dislocation density, leading to characteristic 'butterfly' curves also observed in magnetic systems (35; 29; 36).
There are other aspects of the theory which help to make contact with actual materials. For instance, it is possible within the theory to account for obstacles of different species, and to accord to them varying densities and strengths. This furnishes an avenue for building into the theory results from recent work concerned with dislocation-dislocation interactions (10; 11; 37; 12). However, in order to build additional physical realism into the theory a number of extensions and enhancements immediately suggest themselves, to wit: 1. It should be possible, within the general framework outlined in Section 2, to formulate a fully threedimensional theory of crystallographic slip accounting for the full complement of slip systems in a crystal class and, e. g., an infinite stack of uniformly-spaced slip planes within each slip system.

2.
A topological restriction of the phase-field representation of crystallographic slip, as presented here, is that, since the dislocation lines follow the level contours of a scalar phase field, they must necessarily trace closed loops or form networks of zero net Burgers vector. This topological constraint precludes the effective representation of cross slip, and therefore points to the need for non-trivial mathematical extensions of the theory.
3. In the formulation presented in this paper, the inelastic interaction between dislocations and the crystal lattice, and, between dislocations and obstacles, has been modeled as a rate-independent frictional interaction in the interest of analytical tractability. However, the incremental variational framework developed in (4; 5; 6; 7; 8) is general enough to enable consideration of more complex kinetic relations, e. g., accounting for the temperature and strain-rate dependence of the Peierls stress.
4. As presented here, the theory accounts for the edge/screw line-energy anisotropy predicted by linear elasticity. However, this type of anisotropy does not suffice to model the vastly different behavior of screw and edge segments observed in bcc crystals, which largely owes to differences in the mobility of both types of segments. Within the present theory, a proper accounting of this effect can be achieved simply by allowing the Peierls stress to depend on the gradient ∇ξ of the phase field, which describes the local orientation of the dislocation segments within the slip plane.
Finally, the fact that the theory permits the reconstruction of the full dislocation ensemble from the value of the phase field at the obstacle sites can hardly be overemphasized. Indeed, the geometry of each dislocation line is given by the theory analytically and in close form. In addition, the same analytical treatment endows each dislocation with a well-defined core, which effectively regularizes -and eliminates the divergences which afflict -linear elastic dislocation theory. The practical effect of this analytical tractability is the elimination of any need for the use of numerical grids to discretize the slip plane or, worse, the entire crystal. Since the theory targets the equilibrium configurations of the dislocation ensemble directly, no costly resolution of transients, often introduced as numerical artifacts in order to advance the solution in time and requiring the use of exceedingly small time steps, is required. These attributes render the complexity of the calculations required by the theory, in its present state of development, commensurate with the number of obstacles. This drastic reduction in complexity may open the way for embedding the theory within a finite-deformation formulation of single-crystal elastic-plastic behavior (e. g., (7; 38)), with a view to its use in large-scale finite element calculations of macroscopic samples.