Phase-Field Modeling of Chemoelastic Binodal/Spinodal Relations and Solute Segregation to Defects in Binary Alloys

Microscopic phase-field chemomechanics (MPFCM) is employed in the current work to model solute segregation, dislocation-solute interaction, spinodal decomposition, and precipitate formation, at straight dislocations and configurations of these in a model binary solid alloy. In particular, (i) a single static edge dipole, (ii) arrays of static dipoles forming low-angle tilt (edge) and twist (screw) grain boundaries, as well as at (iii) a moving (gliding) edge dipole, are considered. In the first part of the work, MPFCM is formulated for such an alloy. Central here is the MPFCM model for the alloy free energy, which includes chemical, dislocation, and lattice (elastic), contributions. The solute concentration-dependence of the latter due to solute lattice misfit results in a strong elastic influence on the binodal (i.e., coexistence) and spinodal behavior of the alloy. In addition, MPFCM-based modeling of energy storage couples the thermodynamic forces driving (Cottrell and Suzuki) solute segregation, precipitate formation and dislocation glide. As implied by the simulation results for edge dislocation dipoles and their configurations, there is a competition between (i) Cottrell segregation to dislocations resulting in a uniform solute distribution along the line, and (ii) destabilization of this distribution due to low-dimensional spinodal decomposition when the segregated solute content at the line exceeds the spinodal value locally, i.e., at and along the dislocation line. Due to the completely different stress field of the screw dislocation configuration in the twist boundary, the segregated solute distribution is immediately unstable and decomposes into precipitates from the start.


Introduction
The dependence of the material properties on chemical composition, temperature and pressure (stress) is central to the phase relations, thermodynamics and behavior of many materials. In the case of engineering alloys, for example, the dependence of elastic or magnetic material properties on chemical composition can have a significant influence on alloy thermodynamics, phase relations, and mechanical behavior due for example to defects such as dislocations.
A classical example of the latter, going back to the seminal work of Cottrell (e.g., Cottrell and Bilby [1]; see also Hirth and Lothe [2]), is the composition dependence of lattice distortion due to solute misfit, resulting in a contribution of the stress field to the solute chemical potential, and so to the driving force for solute diffusion, for example to lattice defects such as dislocations (e.g., Kuzmina et al. [3], Kwiatkowski da Silva et al. [4], Kwiatkowski da Silva et al. [5], Zhou et al. [6]). More recent work has focused on further aspects and details of this type of chemoelastic coupling. In Ma et al. [7], solute segregation and which order pertains. The scalar product A · B of two arbitrary-order tensors A and B is defined by A · B := A ijk... B ijk... (contraction; sum on repeated indices). In particular, then, a · b = a i b i represents the scalar product of two vectors. Given this, A T b · c := b · Ac defines the transpose A T of A, sym A := 1 2 (A + A T ) its symmetric part. Let v and T be differentiable tensor fields. The curl of these can be defined by a · curl v := div(v × a) and (curl T) T a := curl(T T a) with respect to any constant vector a. Additional notation and relations will be introduced as needed.
In addition to the balance relations (1), the current model for binary alloy chemomechanics is based on the general constitutive form for the free energy density ψ := ε − θη. In this case, the molar form f of ψ is also determined, with υ m the (here assumed constant) alloy molar volume (units m 3 mol −1 ). Since attention is focused in this work on purely bulk behavior, the generalized no-flux boundary conditionṡ hold on the boundary of any region with unit normal n. In addition, the dependent constitutive relations apply, with δ x g := ∂ x g − div ∂ ∇x g the variational derivative of g. Given non-negative dislocation m a (units J −1 m 3 s −1 ) and solute m c (units J −1 mol m 2 s −1 ) mobilities, (3) and (4) are sufficient for non-negative δ.
Analogous to the purely chemical case (e.g., Cahn and Hilliard [25]), the chemomechanical free energy density ψ(F, c, φ, ∇c, ∇φ) = ψ ho (F, c, φ) + ψ gr (φ, ∇c, ∇φ) (5) in (2) is modeled here as the sum of "homogeneous" ψ ho and "gradient" ψ gr contributions. The former is determined in general by elastic (lattice) ψ el , dislocation ψ hd , and chemical ψ hc , parts, respectively. The latter ψ gr (φ, ∇c, ∇φ) = ψ gd (φ, ∇φ) + ψ gc (∇c) (7) consists in general of dislocation ψ gd and chemical ψ gc parts. The material properties determining ψ gd are also generally dependent on c, but this is neglected here for simplicity. The form of ψ gc is based for simplicity here on pairwise interaction and regular solution theory, in which case ψ gc is independent of c (e.g., Cahn and Hilliard [25], Lass et al. [27]). All of these contributions to ψ are discussed in more detail in what follows.

Local Kinematics and Elastic Energy
Let be the lattice (elastic) local deformation, right Cauchy-Green deformation, and Green strain, respectively. The residual (i.e., zero-stress) local deformation in the current model is determined by contributions from solute misfit F S and dislocation motion F D on g systems a = 1, . . . , g. Here, N S represents the (infinitesimal, linear) distortion per unit solute concentration due to solute misfit. Restricting attention to dislocation glide, holds. Here, γ a := b a /d a , γ a := b a /d a , b a := |b a |, s a := b a /b a , b a is the Burgers vector, d a the lattice slip plane spacing, and n a the slip plane normal. In the current work, attention is restricted to the special case that s a · n b = 0 for a = b, resulting in the simplified form of the third relation in (9) based on (10). Assuming "small" lattice strain |E L | 1 (i.e., outside the dislocation core), the harmonic form for ψ el applies, with C el the elastic stiffness, assumed independent of c for simplicity. From (9) and (12) follow the forms for the concentration and phase-field derivatives of ψ el , where represent the Kirchhoff, "lattice" second PK, and Mandel, stresses, respectively. For later purposes, it will be useful to compare this non-linear elastic model with its linear counterpart. As usual, this is based in particular on the displacement gradient H := ∇u = F − I and linear lattice distortion H Ll := H − H Rl , with H Rl = H Sl + H Dl = c N S + ∑ g a=1 φ a γ a ⊗ n a corresponding to (9) with (10) in the non-linear case. Given these, the linear elastic energy density takes the form with E Ll := sym H Ll = E − E Rl , and E := sym H the linear strain. From (15) follow corresponding to (13) and (14) in the non-linear case, where T is the linear elastic stress. In contrast to the third relation in (13) for ∂ 2 c ψ el in the non-linear model, note that the third relation in (16) for ∂ 2 c ψ le in the linear case is independent of lattice strain. As discussed in more detail below, this has consequences for the system binodal and spinodal in the current chemomechanical context.

Dislocation Energy
In general, the homogeneous dislocation energy ψ hd is related to energy barriers resulting in lattice resistance to dislocation transitions and motion. Examples of this include the stacking fault energy with respect to {111} glide planes (of spacing d 111 = a 0 / √ 3) in the fcc case, or the energy related to (screw) core spreading in the bcc case.
Restricting attention to planar dislocation cores for simplicity, the dislocation core energy ψ gd is modeled here by the simple quadratic from ψ gd (φ, ∇φ) = ϕ gd |G D (φ, ∇φ)| 2 in terms of the gradient energy coefficient ϕ gd (units J m −1 ) and dislocation tensor G D := curl F D . As already mentioned above, ϕ gd will also depends on c in general, but this is neglected here for simplicity, and ∂ c ψ gd = 0. Given further (11) Adopting the initial conditions F D0 = I, ∇φ 0 = 0, and G D0 = 0, the linearized form G D (∇φ) = ∑ g a=1 γ a ⊗ ∇φ a × n a of G D holds. On this basis, the simplified model relations are employed in the sequel.

Chemical Energy
For simplicity, attention is restricted here to disordered phases, crystalline regular solid solution theory, and pairwise interaction. In this case, hold for the homogeneous f hc part of the chemical energy per unit mole and its gradient part f gc , respectively. These depend on the molar ab bonding energy e ab (units J mol −1 ), the relative interaction energy w AB := 2e AB − (e AA + e BB ), and the energy modulus N gc (units J m 2 mol −1 ). Note that R in (19) 1 is the universal gas constant.

Driving Forces for Solute Flux, Chemomechanical Binodal and Spinodal
The current energy model (5)-(7) determines the forms for the (chemomechanical) chemical potential µ from the first relation in (4) and its spatial gradient, respectively. In particular, the latter represents thermodynamic force driving solute flux j via the third relation in (4), and so solute segregation. For example, the c dependence of f hd drives Suzuki [28] segregation (i.e., to stacking faults) in the fcc case.
Analogously, that of f el due to solute misfit lies behind Cottrell segregation. A second consequence of (5)- (7) are the relations (recall that µ = 0 implies µ A = µ B ) for the chemomechanical binodal and spinodal hypersurfaces, respectively, in (c, F, φ) space. In the current model, then, both the binodal and spinodal deviate from their purely chemical counterparts ∂ c f hc = 0 and ∂ 2 c f hc = 0 due to the solute concentration dependence of the dislocation f hd and elastic f el contributions to the energy of the binary alloy. Recall that the former is related to the stacking fault energy in the fcc case, or to the (screw) core spreading energy in the bcc case. In the binodal case, both the non-linear elastic relation in (13) for ∂ c f el , and linear elastic relation in (16) for ∂ c f le , predict a dependence of the chemomechanical binodal on the (non-linear, linear) lattice strain state via the stress. As shown by the third relation in (13) for ∂ 2 c f el and the third relation in (16) for ∂ 2 c f le , such a dependence is also predicted for the chemomechanical spinodal by non-linear elastic model, but not by the linear elastic one. These aspects of the current model are examined more closely in the following with the help of simulation.

Reduction to Cubic Symmetry
All analytical and simulation results to be discussed in what follows are for the case of cubic single crystals. In this case, the misfit distortion per unit solute concentration N S takes the cubic form where υ S is the scalar dilatation per unit concentration. Given (22), the non-linear elastic relations in (13) reduce to with I · A = A 11 + A 22 + A 33 and k el := I · C el I = 3(C 11 + 2C 12 ) for cubic symmetry. Analogously, the linear elastic relations (16) 1,3 reduce to via (16) 4 in the cubic case. Since I · H Dl = 0 for dislocation glide, note that I · E Ll = I · H − 3 ν S c. Since dislocation glide does contribute to I · E L , this is another difference between the non-linear and linear models. Additional simplifications in the cubic case include that for N gc in (19), where a 0 is the lattice spacing in the solvent, and κ gc the chemical gradient energy coefficient (units J m 2 mol −1 ). In this case, the chemical energy (19) reduces in essence to the cubic Cahn-Hilliard (CH) form [25]. To emphasize that we are working with the CH model for the chemical energy, let represent the chemical part of f in what follows based on e BB = e AA . Recall that the formulation of Cahn and Hilliard [25] is based on energy per atom, rather than energy per mole as in the current work.

Non-Dimensional Model Relations
Scaling is based as usual in particular on a typical length 0 (e.g., system size) and time t 0 . In what follows, g * := g/g 0 represents the scaled/non-dimensional form of any quantity g. In particular, ∇ * := 0 ∇ is the non-dimensional gradient operator. Given these, the CH chemical energy (per unit mole) (26) takes the form with a * 0 := a 0 / 0 and θ 0 := w AB /2R. Likewise, one obtains the scaled forṁ of the model field relations from (1) 1,2 and (4) 2,3 , where are typical timescales for solute diffusion and dislocation glide, respectively. Note that m c −1 0 µ 0 represents the solute diffusivity corresponding to m c . In the following, the typical length 0 is determined by the largest system/simulation cell size, e.g., L z = 160 a 0 in the simulations to be discussed below. For a typical fcc lattice constant a 0 = 4 × 10 −10 m, for example, this implies 0 ∼ 10 −7 m, which is adopted here. In addition, for the case of solute segregation to static dislocations, m a = 0, t a = ∞, and t c is the material timescale of interest. To facilitate investigation of solute interaction with moving dislocations, solute diffusion is assumed to be much faster than dislocation glide, i.e., t c t a . In all cases, final results are based onċ * = 0 andφ * a = 0 on the timescale t 0 .

Numerical Solution of Initial-Boundary-Value Problems Based on MPFCM
This is based in particular on the "weak" form Ubachs et al. [29], Shanthraj et al. [30] f wgc (c,c, ∇c, φ) : of the gradient chemical energy in terms of the auxiliary fieldc and penalty parameter α.
In this context, the difference betweenc and c is minimized via minimization of the last two terms in (30) with respect toc. As usual, the corresponding Euler-Lagrange relation is necessary for this and provides a field relation forc. In the context of (30) and (31), note that µ = ∂ c f ho + α (c −c) approximates the first relation in (20) for the chemical or diffusion potential. Numerical solution of the independent field relations in (1) and (31) is carried out in a staggered fashion. Initial conditions here include uniform solute concentration in each case. Boundary conditions include zero external loading (stress control). Iteration proceeds untilċ = 0. To minimize the difference c −c, a large value α * = 10 8 of α * is employed in all simulations. Changing this value an order of magnitude either way has no influence on the simulation results.
As in the case of previous applications of MPFCM, e.g., to the modeling of dislocationsolute interaction in Ni-Al-Co in [23], the model is implemented as a module in the simulation software toolkit DAMASK. This is an open-source toolkit for the numerical solution of initial-boundary-value problems based on coupled field relations like (1) 1,2 and (4) 4 with (4) 1−3 . Numerical solution based on both finite-element and spectral methods is employed. For more information, the interested reader is referred to the DAMASK website https://damask.mpie.de (accessed on 1 April 2021).

Simulation Set-Up
Unless otherwise stated, all simulation cells are fully periodic and cubic with cell side vectors (L x i x , L y i y , L z i z ). In the case of fcc edge dislocations, for example, Dislocation simulations assume initially perfect edge dislocation dipoles with glide plane normal i y .

Linear Chemoelastic Binodal and Spinodal in Defect-Free Cubic Crystals
In this case, the dislocation contributions ψ hd and ψ gd to ψ are zero, and (5)-(7) reduce to ψ = ψ el + ψ hc + ψ gc , with the sum of the latter two given by (26) for f CH . Note that the homogeneous chemical part of (27) for f * CH determines the forms for the non-dimensional chemical binodal temperature θ * CHb , the non-dimensional chemical spinodal temperature θ * CHs , and the chemical spinodal points c CHs± , respectively. In the linear chemoelastic case and (24), note that (33) generalize to the chemoelastic forms for dilatation d(H) := I · H control, and those for (non-dimensional) hydrostatic stress σ * h (T) := 1 3 k −1 el I · T control. Clearly, (34) and (35) reduce to (33) for υ S → 0. Recall that (34), and in particular (35), are based on neglecting the dependence of C el and υ S on c.
Recall that lim c→ 1 2 θ * CHb (c) = 1 = lim c→ 1 2 θ * CHs (c) represents the so-called (here lower) critical (solid) solution or consolute temperature. In the linear chemoelastic case, lim c→ 1 2 θ * leb (d(H), c) diverges for υ S = 0 except at the "critical" dilatation d(H) = 3 2 υ S . Indeed, at this dilatation, the second term in (34) 1 vanishes, and lim c→ 1 2 θ * leb (c, 3 2 υ S ) = 1 − 1 2 v m k el ν 2 S /w AB holds. Note that this dilatation corresponds to I · T = 1 3 k el I · E Ll = 1 2 k el υ S (1 − 2c). In this case, lim c→ 1 2 σ * h (T)/ tanh −1 (1 − 2 c) = 1 6 ν S ; (34) 1 and (35) 1 are then consistent. Figure 1 displays θ * leb (c, 3 2 υ S ) and θ * les (c) for selected values of υ S . At least at the critical dilatation, then, the binodal region (miscibilty gap) and its spinodal counterpart decrease with increasing solute misfit. For the spinodal region, this is also shown in Figure 2. Clearly, there is a signficant effect of elasticity, and in particular of solute misfit, on the binodal and spinodal in the context of linear chemoelasticity.  In the non-linear chemoelastic case ψ = ψ el + ψ CH , the corresponding binodal and spinodal are determined by (23) 1 for ∂ c ψ el rather than by ∂ c ψ le from (24) 1 . This nonlinear form is employed in all simulations in the sequel. Further, these are based on the scaling choices for the driving forces in (28) 1,3 . The value w * AB := w AB /v m k el = 10 −3 of w * AB employed in Figures 1 and 2 is adopted as well in what follows. Further, we work with θ * = 0.5, and the typical values C 11 /ψ 0 = 1.5 × 10 −1 , C 12 /ψ 0 = 9 × 10 −2 , C 44 /ψ 0 = 8 × 10 −2 , and ν S = 2 × 10 −2 , for an fcc metal. In this case, note that ψ * el ∼ 10 −4 for |E L | ∼ 10 −2 . These and other typical non-dimensional values employed in the simulations to follow are listed in Table 1.  The last three non-dimensional parameter values are related to dislocation dissociation (i.e., ψ * hd 0 ) and core energy (i.e., ψ * gd 0 , ψ * gd 1 ) relevant to the case of solute interaction with a gliding dislocation and discussed in more detail in Section 5.4 below. Note that the enhanced solute mobility at the dislocation core is neglected in this work, i.e., solute mobility m c does not depend on dislocation order parameters φ. However, as it will be shown in the dynamic case, solute drag due to dislocation motion is automatically captured by the model.

Single Static Perfect Edge Dislocation
For simplicity, the simulation examples to be discussed in the following three subsections neglect the dislocation core energy ψ gd in (7). In addition, dislocations involved are assumed to be of perfect Peierls-Nabarro (PN)-type, in which case ψ hd is of Frenkel potential-type (e.g., Hirth and Lothe [2], Schoeck [31]). Then g = 1, and the planar dislocation phase field/disregistry φ 1 (x) = φ PN (γ 1 · x) is determined by the analytic PN arctanbased disregistry φ PN (e.g., Hirth and Lothe [2]). Note for example that γ 1 = √ 3/2 s 1 in the fcc case. Further, m a = 0 (t a = ∞) in the static case as mentioned above. Lastly, an initially uniform solute concentration with c(0) = 0.11 is assumed in all cases.
Results for segregation to a perfect PN edge dislocation for two system/simulation cell sizes are shown in Figure 3. . Snapshots of solute segregation to a perfect Peierls-Nabarro edge dislocation dipole for larger (L x , L y , L z ) = (80, 20, 160) a 0 (above) and smaller (L x , L y , L z ) = (80, 10, 160) a 0 (below) simulation cells. The cell orientation is as given in (32). Dislocation lines are shown in yellow, and the 40% solute concentration isosurface in green, and t 0 = 10 2 t c . Note that segregated solute is below (above) the left (right) monopole, corresponding to the region of positive hydrostatic stress. See text for discussion.
Initially, solute segregation to the dipoles results in a uniform solute distribution (with maximum concentration c = 0.97) along the lines in both systems (left two snapshots). Whereas this distribution is stable in the larger system (above) up to t * = 6.2, it decomposes into a single cylindrical precipitate at one of the monopoles in the smaller one (below). In contrast to the larger system (above), the smaller system (below) contains too little solute (about 5.4 × 10 5 solute atoms below, 1.1 × 10 6 atoms above) for segregation alone to stabilize the uniform solute distribution along the monopoles against spinodal decomposition and precipitation for t * > 2.1. From the point of view of statistical thermodynamics, the larger system is more grand-canonical-like, and the smaller more canonical-like, in its behavior.

Tilt Boundary
Consider next the case of solute segregation to, and precipitation at, an array of static PN edge dislocation dipoles of the type from the last subsection. To this end, a simulation cell of size (L x , L y , L z ) = (80, 60, 160)a 0 is employed. Three glide plane/dipole spacings of d = 6 a 0 , d = 10 a 0 and d = 20 a 0 result in low-angle grain boundaries (LAGBs) of tilt angle 6.77 • , 4.05 • , and 2.03 • , respectively. Results for segregation of an initially uniform solute distribution to these arrays are displayed in Figure 4. . Snapshots of solute segregation to three different low-angle tilt grain boundaries. As before, dislocation lines are displayed in yellow, and solute 40% concentration iso-surfaces in dark green. In (II), the 10% solute iso-surface is displayed in light green, and t 0 = 10 2 t c .
For case (I) with tilt angle 2.03 • , the separation between the dipoles is such that their interaction is relatively low, and segregation takes place to each as an essentially isolated dipole. In addition, the system size is sufficiently large for the uniform segregated solute distributioin to stabilize against spinodal decomposition and preciptate formation along the line, analogous to the behavior in Figure 3 (top row). For case (II) with tilt angle 6.77 • , the dipoles are sufficiently close to each other that their stress fields shield each other, leading to an effective reduction of stress field strength and less segregation. Indeed, the maximum solute concentration here is about 0.37, which is much lower than in cases (I) and (III).
For the intermediate case (III) of tilt angle 4.05 • , the situation is similar to that of the single static dislocation in Figure 3 (bottom row). Indeed, increasing the number of dislocations in the system at constant system size effectively reduces the system size per dislocation. Note that the dislocation density is 7.8 × 10 15 in (I), 2.6 × 10 16 in (II), and 1.6 × 10 16 1/m 2 in (III). Consequently, solute content limitation is again stronger, and segregation alone cannot stabilize the initially uniform solute distribution along the monopoles against spinodal decomposition. As seen starting at t * = 2.3 in Figure 4e-g, because of this, the uniform distribution along the lines becomes unstable and precipitate formation leads to solute depletion along the lines. The resulting precipitates have maximum solute concentration close to the bulk binodal (about 0.9).

Twist Boundary
Analogous simulation is carried out for the case of low angle twist boundary. These boundaries often result in a network of screw dislocations, in this case a hexagonal network. Note that to satisfy the periodicity, the simulation box is divided into four sections with four twist boundaries with opposite twist angles, resulting in zero sum (analogous to the tilt boundary dipole above). The snapshots of this simulation are shown in Figure 5. In contrast to (linear elastic) continuum dislocation theory, non-linear effects accounted for in MPFCM result in a non-zero hydrostatic stress field in screw cores driving segregation to these as well. Note that atomistic modeling based on hybrid Monte Carlo molecular dynamics (Sadigh et al. [16]), or on diffusive molecular dynamics (e.g., Dontsova et al. [32], Ponga and Sun [33]), also predict segregation to screw dislocations. Again, this is in contrast to continuum modeling based on linear elasticity.
Although much smaller than its edge counterpart, the hydrostatic screw stress field is sufficient to drive Cottrell segregation to these as well. As shown in Figure 5, this is also the case for screw configurations like a twist boundary. In this latter case, maximum positive hydrostatic stress, and so segregation, appears to be at junctions where four sections meet and the twist angle is reversed. Due to the completely different stress field of the screw configuration in the twist boundary, note that the segregated solute distribution is apparently immediately unstable and decomposes into precipitates from the start, in contrast to the single static edge dislocation ( Figure 3) and tilt boundary (Figure 4) cases.

Single Gliding Dislocation
As a last application, consider solute segregation to, and interaction with, a single gliding dislocation. In contrast to the previous examples, the dislocation energies ψ hd in (6) and ψ gd in (18) now play a role. Focusing attention here on the fcc case, for a single edge dislocation, g = 2, (s 1 , n 1,2 , In the first of these, γ sf (c, φ 1 , φ 2 ) is the fcc stacking fault energy whose representation is given for example in Mianroodi et al. [23]. For simplicity, the c dependence of γ * sf (c, φ 1 , φ 2 )/d * is assumed linear with scaled slope of 3 times the value of ψ * hd 0 in Table 1, resulting in a negative driving force for Suzuki segregation to stacking faults.
The corresponding results are shown in Figure 6. Figure 6. Solute segregation to, and spinodal decomposition at, an fcc edge dislocation gliding under an external shear deformationF xy = 0.10 in a system/simulation cell with orientation (32). 40% solute concentration isosurface shown in green. Dislocation line visualization is based on the scalar field |∇φ 1 | + |∇φ 2 |, which varies between zero (blue) and maximal (red) in the dislocation core. Note that the external shear deformation is applied at a rate much faster than t −1 c . In addition, t c = 10 −2 t a , and t 0 = 10 2 t a . See text for discussion.
At the start (t * = 0.3), solute segregates to the initially perfect dislocation dipole. Each monopole of the dipole then dissociates into Shockley partials (t * = 3.3) and begins to glide (t * = 4.0). In the process, the initially uniform solute distribution along each monopole destabilizes, driven in part by negative Suzuki segregation. Due to the boundary conditions, this takes the form of complete solute depletion along the left dissociated monopole, precipitate formation at the right dissociated monopole, and solute depletion along the rest of this monopole (4.0 t * 6.6). The higher stacking fault energy inside the precipitate results in a reduction of the stacking fault width between the partial dislocations in the right monopole and a bend in the dislocation lines at the precipitate interface.
Under the current boundary conditions, the force due to spinodal instability driving precipitate formation is stronger that the Cottrell force on solutes due to misfit and the positive hydrostatic core stress field attracting them to the dislocation line. Note that solute distribution due to precipitation increases the spatial separation between solute in the precipitate and the dislocation core, resulting in a reduction of the Cottrell force on solute and no solute transport due to dislocation glide (t * 4.6). Lastly, as glide continues (t * > 8.3), the leading partials on both sides interact with their periodic images, resulting in partial annihilation of the leading partials (t * = 9.0).

Summary and Discussion
Microscopic phase-field chemomechanics (MPFCM) has been employed in the current work to model solute segregation, dislocation-solute interaction, spinodal decomposition, and precipitate formation, at straight dislocations and configurations of these in a model binary solid alloy. In particular, (i) a single static edge dipole, (ii) arrays of static dipoles forming low-angle tilt (edge) and twist (screw) grain boundaries, as well as at (iii) a moving (gliding) edge dipole, have been considered. MPFCM is formulated for such an alloy in the first part of the work. Central here is the MPFCM model for the alloy free energy, which includes solute, dislocation, and elastic lattice, contributions. Due to solute lattice misfit, the latter energy is concentration dependent, resulting in a strong elastic influence on the binodal and spinodal behavior of the alloy. In addition, MPFCM-based modeling of energy storage couples the thermodynamic forces driving (Cottrell and Suzuki) solute segregation, precipitate formation and dislocation glide. As implied by the simulation results for edge dislocation dipoles and their configurations, there is a competition between (i) Cottrell segregation to dislocations resulting in a uniform solute distribution along the line, and (ii) destabilization of this distribution due to low-dimensional spinodal decomposition when the segregated solute content at the line exceeds the spinodal value locally, i.e., at and along the dislocation line. Due to the completely different stress field of the screw dislocation configuration in the low-angle twist boundary, the segregated solute distribution is immediately unstable and decomposes into precipitates from the start.
Like in previous works based on linear elasticity, the dependence of the elastic energy on solute misfit in the current non-linear treatment is central to the influence of elasticity on (Cottrell) segregation, spinodal decomposition, and precipitate formation. As shown by the treatment in Section 5.1, in this case, the binodal and spinodal depend in a constitutive fashion on the strain or stress (i.e., in addition to the solute concentration) in the (linear) chemoelastic context. In particular, this dependence holds in the case of spatially homogeneous solute concentration, strain and stress, the latter satisfying mechanical equilibrium trivially. On the other hand, again in the chemoelastic context, spinodal decomposition represents a transition from spatially homogeneous to inhomogeneous solute concentration, strain and stress. This was realized by Cahn [10] in his ground-breaking work on the role of solute misfit in the spinodal behavior of defect-free metallic alloys. Under the assumption that spinodal decomposition takes place in mechanical equilibrium, he showed that this process is affected by dependence of the (equilibrium) elastic energy on solute misfit not accounted for in the CH model [25].
To discuss this in more detail, consider the split φ =φ +φ of any field φ into meanφ := φ (i.e., volume averaged, spatially constant) and spatially fluctuatingφ parts. In this context, Cahn [10] assumed (i) spatial inhomogeneity in one dimension (x), (ii) isotropic linear elasticity, (iii) no defects, (iv) E Sl = υ Sc I, (v) E = E xx i x ⊗ i x , and (vi) T = T xx i x ⊗ i x + T yy i x ⊗ i y + T zz i x ⊗ i z for the stress. Under these assumptions, mechanical equilibrium reduces to div T = T xx,x i x = 0. Choosing then E Ll = E − E Sl in such a way that T xx = 0, Cahn [10] obtainedψ le = υ 2 S E c 2 /(1 − ν), with E Young's modulus, and ν Poisson's ratio. More recently, Onuki [34,35] extended the treatment of Cahn [10] to multiple dimensions and a dependence of the isotropic elastic constants on solute concentration. As discussed for example by Binder and Fratzl [36], the original 1D treatment of the cubic anisotropic case by Cahn [11] has been extended by Khachaturyan [12] and others to 3D and general anisotropy with the help of the Green-functionbased formal solution (See also [37]; also used in PFM: e.g., [8,9,38].) of linear elastostatic mechanical equilibrium div T = 0. In particular, this yields the form (Here and in what follows, the operator * represents both convolution and linear mapping.) E Ll =Ē Ll + M le * Ẽ Rl for the equilibrium lattice strain with M le := Γ le C el − I, where Γ le is the linear elastostatic Lippmann-Schwinger operator (Γ le (k) A := sym (Ĝ le (k) A(k ⊗ k)), with G le the corresponding Green function (Ĝ −1 le (k) a := C el [a ⊗ k]k).) (e.g., [39]). In turn,ψ le = 1 2Ē Ll · C elĒLl + 1 2 Ẽ Rl · A le * Ẽ Rl follows withÂ le (k) :=M T le (k) C elMle (−k). For the current case and H Rl = c N S + ∑ g a=1 φ a N Da , then,ψ le is determined in particular by the term 1 2 c N S · A le * c N S , representing a 3D anisotropic generalization of the above result of Cahn [10]. In addition, one obtains the elastic contribution δ cψle = −N S ·T − N S · A le * Ẽ Rl to the solute chemical/diffusion potential µ satisfying mechanical equilibrium in the spatially inhomogeneous case, in contrast to its constitutive counterpart ∂ c ψ le = −N S · T from (16) 1 . This is likewise the case for the elastic contribution δ 2 cψle = N S · C el N S + N S · A le * N S to the chemoelastic spinodal. This is of course also true for all simulation results in the current work based on elastic non-linearity and MPFCM.
Further analogous generalizations of the treatment of "closed" and "open" solid solution chemoelasticity in Larché and Cahn [40] to (i) non-linear chemoelasticity and (ii) non-ideal (e.g., regular) solutions, are also possible and represent work in progress.
Author Contributions: All authors contributed equally to conceptualization, methodology, manuscript writing, review and editing; funding acquisition P.S., B.S.; modelling, execution of simulations, data analysis and writing the corresponding J.R.M. and B.S. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.