Applying diffusion coupled deformation theory, we investigate how the elastic properties of a solid body are modified when forced to keep its chemical potential aligned with that of its melt. The theory is implemented at the classical level of continuum mechanics, treating materials as simple continua defined by uniform constitutive relations. A phase boundary is a sharp dividing surface separating two continua in mechanical and chemical equilibrium. We closely follow the continuum theory of the swelling of elastomers (gels) but now applied to a simple two phase one-component system. The liquid is modeled by a local free energy density defining a chemical potential and hydrostatic pressure as usual. The model is extended to a solid by adding a non-linear shear elastic energy term with an effective modulus depending on density. Imposing chemomechanical equilibrium with the liquid reservoir reduces the bulk modulus of the solid to zero. The shear modulus remains finite. The stability of the hyper-compressible solid is investigated in a thought experiment. A mechanical load is applied to a rectangular bar under the constraint of fixed lateral dimensions. The linear elastic modulus for axial loading is evaluated and found to be larger than zero, implying that the bar, despite the zero bulk modulus, can support a weight placed on its upper surface. The weight is stabilized by the induced shear stress. The density dependence of the shear modulus is found to be a second order effect reducing the density of the stressed solid (chemostriction).

The thermodynamic state of a simple homogeneous liquid in thermal equilibrium is determined by the temperature and one more thermodynamic field variable, which can be density, pressure, or chemical potential. The value of the total free energy can be found by simply scaling the corresponding values of a reference system in volume or mass. Separate specification of mass and volume is not necessary, the two quantities being essentially equivalent. This equivalence is also reflected in the density. An increase in mass can always be compensated by a proportional increase in volume, leaving the density invariant. This invariance is a key principle of the thermodynamics of uniform liquids, as formalized by Gibbs. (Most of the material relevant to this paper can be found in the work of Hansen and McDonald.1) Density is a natural state variable for the description of an isothermal mechanical process. Alternatively, pressure can also be used.

However, complications arise for solids. An isochoric deformation of a homogeneous solid body, preserving volume and density, can still change its energy. The reason is the shear rigidity of solids. This is the subject of the classical theory of elasticity, which predates the Gibbsian thermodynamics of liquids. The state variables in this theory are strain and temperature. A variation in density is accounted for by an isotropic stretch of a reference state. (Out of the large number of textbooks on classical elasticity, we have selected for citation the recent work by Lubarda and Lubarda.2) Density is not treated as an independent degree of freedom and certainly does not contain sufficient information for the specification of the state of a solid (indeed, strain can be inhomogeneous even if density is homogeneous). There is therefore a disparity in the description of simple uniform liquids and solids. This is not a problem at interfaces between materially closed systems. The boundary conditions of continuum mechanics can take care of this. The question is how to match these different descriptions at the interface between liquid and solid phases if the liquid and solid are in chemical equilibrium, such as is the case if the liquid is the melt of the solid. This is the question addressed in the present contribution.

The study of chemical equilibrium between continuum models of liquid phases is an established field with a long history. This was the main driving force behind the development of the density functional theory (DFT) of inhomogeneous liquids.1,3 The discussion in the present paper is, however, strictly limited to boundaries between uniform phases as described in the framework of Gibbsian thermodynamics. Gibbs extended his scheme by representing the interface as a sharp dividing surface with its own energy scaling with the surface area of the interface.1,4 However, while these surface properties cannot be missing from a realistic description of phase boundaries, interface energy and surface tension will not be considered here. The focus is on how to generalize the dividing surface scheme for an interface between liquid and solid, bridging the difference in thermodynamic continuum theory for the two aggregation states.

This first of all requires an expression for the chemical potential of the solid accounting for the dependence on shear strain. The treatment of the balance of forces at the interface will also have to be refined now, involving the Cauchy stress tensor rather than only the hydrostatic pressure. Equal chemical potential no longer implies equal pressure (even in the absence of surface tension). Instead, chemical and mechanical equilibrium must be treated as separate but coupled boundary conditions. The theory to accomplish this generalization already exists. It is the diffusion coupled deformation (DCD) theory originally developed for the description of diffusion induced stress in solid solutions. The pioneers in this field are Larché and Cahn.5 Another influential theoretical contribution is the 1996 paper by Gurtin.6 For an update on the experimental evidence for Larché–Cahn theory, see Ref. 7. The theory was further extended for the study of the swelling of polymeric gels.8–14 A more recent application is in lithium ion battery research, in particular, for the understanding of the mechanics of charging and discharging of silicon cathodes.15–18 

The investigation in this paper follows the diffusion coupled elasticity theory as presented in the book by Gurtin, Fried, and Anand (GFA)19 almost to the letter (including the notation). There is, however, one difference. The diffusive motion of a metal species in alloys is defined relative to a lattice of sites, which can be occupied by atoms of different species or even vacancies. Similarly, in a gel, the network of elastic polymers in the dry state serves as the reference to quantify the migratory displacement of the solvent. However, the solid in our application is a simple one-component system. There is no external lattice providing a reference. This is not a fundamental objection. As pointed out by Baek and Srinivasa, an equivalent separation can be achieved by a decomposition of the mass current in a convective and diffusive flow.8 The continuum mechanical basis of the theory remains valid. Differentiation between a migratory component (solvent) and an elastic but immobile component (elastomer) is made by imposing certain “swelling” constraints on the balance laws.8,13,14 Here, we carry on without constraints proceeding to the next stage, the design of a constitutive model of the simple solid.

The approach taken in this paper is somewhat unusual from the perspective of a physical chemist. DCD theory is an application of the methods of continuum thermodynamics.19,20 The principle equations in continuum thermodynamics are balance equations involving time derivatives of densities, entropy, and fluxes. Dissipation is built in at a fundamental level (see also Ref. 21). It may seem that using non-equilibrium equations to derive thermal equilibrium values is an unnecessary detour. Variational methods based on free energy minimization are more direct. However, the dynamical route has certain advantages even under equilibrium conditions. System specific constitutive relations are kept separate from generally valid field equations. This guards against inconsistencies, which are more of a risk for “multiphysics” theories, which introduce further state variables in addition to strain.19,22 Mass density in DCD theory can be regarded as such an extended field variable. Moreover, the distinction between density and deformation, as proposed here, is based on a decomposition of the mass current, which is a dynamical quantity. Here, it is perhaps helpful to draw a parallel to the Lorentz continuum theory of polarization for which the foundation is also a dynamical theory, the Maxwell–Lorentz equations.20,23 In fact, as suggested by Baek and Srinivasa, mass density can be treated as the divergence of mass polarization, replacing density as the fundamental degree of freedom.

Mass exchange across a liquid–solid interface also requires flexibility in constitutive modeling, allowing for a natural transition between liquid and solid. The free energy density of simple liquids is due to the short range steric repulsion and long range cohesive attraction characteristic of atomic fluids. These interactions, while treated in a local approximation, are poorly described by linearization of isotropic expansion as used in elastic theory. This suggests to retain the liquid-like local free energy density in the solid, with a different parameterization, and add a shear energy term depending on isochoric deformation. For an isotropic solid, such a model contains just one parameter more determining the strength of the shear elastic energy term. The continuum solid in our model is missing a reference lattice. Setting the shear modulus to zero, therefore, directly turns the solid into a liquid. Unfortunately, this makes the model intrinsically non-linear, forcing us to use the demanding machinery of non-linear continuum mechanics.19,24 Non-linearity is also the price to pay if we want to allow for a density dependence of the shear elastic interaction. This leads, as we will show, to chemostriction, creating a pressure gap between solid and liquid. While small, this effect is of interest as a manifestation of nonlinear elasticity.

Finally, a word about the motivation for the work reported here: While chemomechanical equilibrium at liquid–solid interfaces is of physical interest, what prompted the author to look into diffusion coupled deformation was a theoretical problem encountered in a previous study on electromechanics of dielectric fluids.25 The aim of that investigation was a systematic derivation of the stress tensor of a simple dielectric fluid polarized by an external electric field and to compare the corresponding force density to the force density obtained from density functional treatment of the same free energy functional. The two expressions for this observable quantity, generally known as the Korteweg–Helmholtz force, were found to be identical. More specifically, it was shown that the density obtained by optimization of the deformation for a fixed but arbitrary material density in reference space satisfies the variational DFT equation for density in the current (deformed) space. It is not at all clear whether this is a general principle. The results of the present study suggest that at least for inhomogeneous systems material density and deformation are independent degrees of freedom, which must be computed separately by solving a set of coupled equations. We will return to this issue in the conclusion.

With the fundamental questions mentioned above in mind, the present paper is a joint theoretical and applied investigation. The emphasis is on the thermomechanical methodology. While almost completely gathered from the textbook of GFA, which has become a standard reference in the continuum mechanics literature, this approach may not be familiar to everybody in the physical chemistry community. The derivations in the theoretical sections are therefore more explicit and lengthy than strictly necessary claiming the larger part of this paper. Only Sec. VII is devoted to the application of a supported bar of solid in a bath of its melt. This application is meant more as a thought experiment. The simple continuum model developed here is too primitive for the investigation of realistic solid–liquid phase boundaries. We conclude with a discussion in Sec. VIII, where we also comment on the connection to the DFT of fluids and solids.

The diffusion coupled deformation methodology used in this work is patterned after the theory for the absorption (or release) of solvents by polymeric gels developed by the continuum mechanics community.8–14 (A recent review from a more physical chemistry perspective can be found in Ref. 21.) The solvent and the polymeric network are clearly distinct interacting components of a composite system, partly liquid and partly solid. The dry polymer is usually taken as the reference for both the elastic deformation of the polymer and the migration of the solvent. However, it was pointed out by Baek and Srinivasa that the balance equations can be formally developed, initially assuming that the gel is a one-component system with certain constraints imposed at a later stage8 (see also Refs. 13 and 14). We have adopted their approach leaving out the constraints. The argument presented in this section is a justification for treating the deformation of a one-component system and the density in the reference frame as separate independent degrees of freedom.

Using a notation common in continuum mechanics, the positions of the “material” particles are indicated by vectors x. The body is deformable, and therefore, the mass density ρ(x) can be inhomogeneous. Moreover, the body is in motion with the particles following trajectories that vary smoothly with their location x. As a result, the local density ρ = ρ(x,t) also depends on time. Let Pt be a spatial region that convects with the body; then,
(1)
is the total mass contained in Pt. Since our system consists of only one component, the quantity ρ can be formally set equal to the number density, as in DFT, although we will often refer to ρ as a mass density, such as already in the next sentence. While the material content of Pt, in principle, evolves in time, MPt is, in fact, a constant because mass is conserved by convective flow. It is this property that is modified by diffusion, and we will, therefore, elaborate the conservation of mass in more mathematical detail.
Following GFA, we introduce the notation
(2)
for the time derivative of the convecting integral of a scalar function φ. Applying the Reynolds transport relation, this can be written as19 
(3)
where v is the material (Lagrangian) velocity. φ̇ is the corresponding material time derivative,
(4)
where φ′ = ∂ϕ/∂t is the partial time derivative of φ at fixed x. The material velocity v is a central quantity in continuum mechanics and is best understood in terms of a formal definition relative to a reference frame. Indicating the positions of the particles in the reference by X (capital x), the motion of the body can be represented as a smooth function χ that assigns to each material point X and time t a spatial point
(5)
The second identity is an alternative notation. v can now be written as the time derivative
(6)
The partial time derivation is evaluated, not at a given position x in physical space (also called current space), but at fixed material position X.
The introduction of a reference space and motion map [Eq. (5)] gives the concept of a convecting volume Pt a more precise meaning. If P is some region in reference space,
(7)
The boundary surface of Pt is carried along by the same motion χt as the interior points. The total mass is constant, and therefore, in the notation of Eq. (2),
(8)
Applying the Reynolds relation [Eq. (3)] and then substituting for the material time derivative using Eq. (4), the integral equation [Eq. (8)] can be expanded as
(9)
Since Pt is arbitrary, we are allowed to localize this equation to
(10)
which is the well-known local expression of mass conservation, often referred to as point-wise mass balance compared to the part-wise mass balance of Eq. (9).
All of this is, of course, part of the standard kinematics setup of continuum mechanics (see the book of GFA and also Ref. 24) and reiterated here only to distinguish diffusive motion from convective motion. Driven by the velocity field v, the material particles cannot cross the moving boundaries of Pt, but diffusive flow can. A convecting volume can now acquire or loose mass. Defining a diffusive current h, the mass balance equation [Eq. (9)] is adjusted to
(11)
where we have used the conventional notation V for the boundary surface of a volume V. da is the elementary area of the surface, and n is the outward normal. The sign of the surface integral ensures that a positive h would lead to an influx of the mass, given the convention for the orientation of the normal n. We have also added a supply term, which can be interpreted as insertion of mass from an outside reservoir. Equation (11) can be regarded as an extended part-wise mass balance, which using the divergence theorem can again be localized to
(12)
where current j is a superposition of convective and diffusive motion (see Fig. 1),
(13)
Remember that ρ in this presentation is really a number density, and hence, the currents in Eq. (13) are not momentum densities but particle fluxes.
FIG. 1.

Separation of the total current j = ρv + h in convective mass flow ρv (solid arrows) and diffusive current h (dashed arrows), as written in Eq. (13). On the right is a deformed convecting volume (solid contour) with, on the left, its shape in material (reference) space (dashed contour). The curved arrows indicate the map χt of Eq. (5), moving points X from reference space to the positions x in the deformed body at time t. hR is the material diffusive current obtained by pulling the spatial current h back according to Eq. (24). While convective motion cannot cross the boundaries of the body, diffusive current h can. There is therefore residual diffusive current in reference space but no convective flow.

FIG. 1.

Separation of the total current j = ρv + h in convective mass flow ρv (solid arrows) and diffusive current h (dashed arrows), as written in Eq. (13). On the right is a deformed convecting volume (solid contour) with, on the left, its shape in material (reference) space (dashed contour). The curved arrows indicate the map χt of Eq. (5), moving points X from reference space to the positions x in the deformed body at time t. hR is the material diffusive current obtained by pulling the spatial current h back according to Eq. (24). While convective motion cannot cross the boundaries of the body, diffusive current h can. There is therefore residual diffusive current in reference space but no convective flow.

Close modal
For an additional insight into the separation between convective current and diffusive current introduced in Sec. II A, it is instructive to reformulate the mass balance [Eq. (12)] in terms of density and currents in the reference frame specifying the motion [Eq. (5)]. These transformations are controlled by the deformation gradient19,24
(14)
where χ is the placement map of Eq. (5). Note that the nabla (∇) operator indicates differentiation with respect to coordinates X in the reference space to be distinguished from the grad operator [see, e.g., Eq. (4)] standing for differentiation with respect to coordinates x in physical space containing the deformed body. The two gradients are related by a linear transformation,
(15)
where FT is the transpose of the deformation gradient [Eq. (14)] (FijT=Fji). In Eq. (15) and all following equations, repeated indices are summed over (Einstein convention) unless stated otherwise.
We further require that the determinant of F is definite positive,
(16)
where J is the Jacobian expressing an infinitesimal material volume element dvR in terms of the corresponding infinitesimal volume element dv in current space,
(17)
Therefore, if we want the referential density ρR to satisfy the integral identity
(18)
with Eq. (7) relating the finite volume Pt in current space to the volume P in reference space, we must set
(19)
The inverse of this relation, often referred to as the Euler continuity equation, is a multiplicative decomposition of the deformed space density
(20)
which can be compared to the additive decomposition [Eq. (13)] of the current. Differentiating the integral density [Eq. (18)] to time gives
(21)
where we could pull the time derivative of the material integral inside the integration because P, unlike Pt, is static. In the absence of diffusion [Eq. (8)], the material density ρR(X), while possibly varying with position X in reference space, is independent of time (ρ̇R=0).
To reformulate the diffusion adapted mass balance law [Eq. (12)] in terms of referential variables, we next transform the diffusive current h. In analogy with Eq. (18), we require that the total flux through a moving surface is invariant. Thus, defining a convecting surface St=χt(S), where S is some surface in reference space, we look for an expression of the material diffusive current hR such that
(22)
The transformation rule for oriented surface elements as derived in the book of GFA is
(23)
where we have used the notation FT=F1T. Therefore, to satisfy Eq. (22), we must have
(24)
Equation (24), changing a deformed space flux density to reference space, is known as a Piola transform. The magnitude h of the supply term [Eq. (11)] is a density and is transformed similar to ρ [Eq. (19)], and therefore,
(25)
Substituting Eqs. (21), (24), and (25) into Eq. (11), we obtain
(26)
Localizing and applying the divergence theorem in reference space, we find that
(27)
where the capitalized operator Div is the divergence in reference space coordinates X. The deformation flow has been eliminated. Only the diffusive motion remains, suggesting that material space is the natural realm for DFT.

The decomposition [Eq. (13)] of the mass current j in hydrodynamic flow ρv and diffusive flux h is rather suspect from the perspective of Newtonian hydrodynamics. There is no mass diffusion in one-component fluids (h = 0). The only diffusion process is momentum diffusion, which is manifested as viscosity. This changes when the system is solidified. The mobility of atoms is now hindered. Rigidity sets in, and solids respond to the applied shear stress not by viscous flow but by elastic strain. However, the mobility of atoms in not fully arrested. The thermal concentration of vacancies still allows atoms to change places and migrate over finite distances, given enough time. Under these conditions, the distinction between convective motion (the deformation of the lattice) and particle diffusion is physical even in a one-component system. This was formalized in the seminal paper by Martin, Parodi, and Pershan (MPP) on crystal hydrodynamics.26 MPP argued that vacancy diffusion must be considered the third Goldstone mode complementing the two shear modes to the three modes required by the Goldstone theorem applied to the translational symmetry breaking in crystals. This implies that vacancy fluctuations are not only long time but also long range. This analysis was further elaborated in Ref. 27 and is the subject of an important chapter in the textbook by Chaikin and Lubensky.28 

In addition to point defects, crystal mechanics must account for further defects not displayed by liquids. These are dislocations and other topological defects.28 Topological defects are also observed in mesophases such as liquid crystals and are characteristics of rigidity. Their dynamics can again differentiate between diffusive motion and convected motion. Restructuring of the near singular distortion inside dislocation cores is diffusive and is the rate limiting process for dislocation mobility. The relaxation of the long range elastic deformation induced by the dislocation is faster and can proceed by convective motion. This issue has been addressed in detail in recent studies of dislocation motion using the phase field crystal (pfc) method.29,30 The pfc method is a microscopic approach based on a modulation of the periodic structure of atomic density in crystals and is also used for the modeling of solidification and grain boundaries.31–35 These results suggest that separation between diffusion and convective motion is linked to the microstructure of density and could, therefore, also occur in systems without long range order such as fluids confined by hard walls.

These observations regarding flow in rigid systems are meant to be examples of systems where the decomposition of current [Eq. (13)] has a physical basis. However, in the work presented here, the separation in convective and diffusive currents is treated as virtual. Virtual currents are tools used in continuum thermodynamics to derive equilibrium constitutive relations. This formalism is briefly summarized in Sec. III. It is then applied to obtain thermodynamic derivative expressions for local chemical potential and pressure. In regular (“classical”) inhomogeneous fluids, changes in hydrostatic pressure p and local chemical potential μ are tied together by the relation
(28)
where x in the notation of Sec. II A is the position in current space [see Eq. (5)]. Equation (28) can also be expressed in the gradient form as
(29)
where we have suppressed the dependence on x. Under equilibrium conditions (no currents), liquids are uniform unless perturbed by an external potential ϕ(x). The gradient of ϕ directly enters the Gibbs equation for chemical equilibrium
(30)
and the Euler equation for mechanical equilibrium
(31)
Eliminating gradϕ leads to Eq. (29). Chemical equilibrium in fluids implies mechanical equilibrium and vice versa. For finite mass flow, Eq. (29) is an approximation known as local thermodynamic equilibrium.
Adopting a more thermodynamic perspective, Eq. (28) can be viewed as a local form of the isothermal Gibbs–Duhem relation. The global Gibbs–Duhem equation is a linear relation between the differentials of the thermodynamic fields in a one-component bulk system,
(32)
where V is the volume, N is the number of particles, and S is the entropy. T is the temperature. Dividing by V and setting dT = 0 gives dp = ρdμ similar to Eq. (28). However, as explained in every textbook on thermodynamics, there is more to Eq. (32). In bulk systems, pressure can equated to a volume derivative of the Helmholtz free energy F,
(33)
Similarly for chemical potential, we have
(34)
The global Gibbs–Duhem relation [Eq. (32)] is therefore more than only an equilibrium condition. It is a profound statement about the extensivity of the thermodynamic limit of bulk systems.

Can the local Gibbs–Duhem equation [Eq. (28)] be given a similar interpretation in terms of some kind of scaling relation? To answer this question, we first need a local equivalent of the thermodynamic derivative relations [Eqs. (33) and (34)]. This is the purpose of the Lagrangian density factorization [Eq. (20)]. It enables us to treat the local particle number represented by ρR and volume represented by J as independent field variables. The flux decomposition [Eq. (13)] is needed to define the corresponding adjoint thermodynamic forces using the thermomechanical scheme for diffusion coupled deformation, as developed by GFA. This procedure is outlined in Sec. III. Application to the simple liquid in Sec. IV reproduces, as expected, Eq. (28), which in the process is elevated to a constitutive relation for a simple fluid. This result is valid whether the partition of current in a diffusive and convective component is physical or not.

For solids, the status of a local Gibbs–Duhem relation is more involved. The Euler equation [Eq. (31)] is replaced by the Cauchy equation,
(35)
where T is the Cauchy stress tensor. The chemical equilibrium equation [Eq. (30)] remains unchanged. As a result, instead of Eq. (29), we now have
(36)
which can be interpreted as a generalization of the local Gibbs–Duhem relation [Eq. (29)]. However, while mechanical equilibrium again implies chemical equilibrium, the reverse is no longer true. Satisfying Eq. (30) still leaves the transverse component of the stress tensor unresolved.

Transverse stress is evaluated by solving the Cauchy equation taking account of the boundary conditions imposed by mechanical surface loading. In fact, body forces can be ignored for typical elasticity problems.2 The transverse stress is all that matters. It will also turn out the key point in the application presented in Sec. VII. Moreover, substituting for p in Eq. (29), the mechanical pressure p = −trT/3 is, in general, not correct either. The divergence of a (symmetric) tensor is not equal to the gradient of its trace. Ultimately, this is all again a consequence of rigidity. However, note that these complications not only apply to solids but also can affect inhomogeneous liquids. While liquids yield to shear stress, the equilibrium stress tensor at liquid interfaces can still have a deviatoric component. This has led to some controversy in the DFT of liquids.4 For a far reaching discussion, we refer to the papers by Percus and co-workers36,37 (see also the work of Lovett and Baus38,39).

The analysis in this paper remains strictly within the confines of continuum mechanics. Continuum thermodynamics is a self-contained theory, assuming a set of axioms (see the book of GFA). Derivations can, in principle, be carried out independently of an atomistic picture. A microscopic theory is needed at the application stage to fix constitutive relations for specific systems. The alternative is to obtain these relations from the experiment. Of course, a foundation in statistical mechanics is the ultimate goal in condensed matter science. This applies, in particular, to the separation in a convective and diffusive current for one-component systems introduced in Eq. (13). Fuchs and co-workers made significant progress toward this end for the case of crystal hydrodynamics.40,41 This concerns the problem raised by MPP macroscopic theory mentioned earlier.26,27 The question for microscopic theory is how to discriminate between the density in a crystal with vacancies and lattice deformation. Häring et al. designed a deterministic coarse graining procedure expressing these quantities in terms of instantaneous atom position. The advantage of defining hydrodynamic field variables by using a deterministic coarse graining scheme is that the fluctuations of these quantities can be analyzed by the methods of statistical mechanics and related to microscopic correlation functions.41 The coarse graining scheme of Ref. 40 is based on reciprocal space methods. An up-to-date review on real space deterministic coarse methods can be found in Ref. 42.

Continuum thermodynamics is a non-equilibrium theory in which fluxes play a central role. This is how we justified in Sec. II the inclusion of material density in the set of state variables describing the mechanical evolution of a single component system. The dynamical nature of continuum thermodynamics will again be used to derive equilibrium equations for these variables. Alternatively, the chemomechanical equilibrium conditions used for the application in Sec. VII could also be obtained using (free) energy-based variational methods. This is, in fact, the procedure followed by Baek and Srinivasa.8,13 However, the DCD theory of GFA is deeply rooted in non-equilibrium continuum thermodynamics, exploiting the constraints set by entropy production (see below). We will adhere to this philosophy and practice mainly because the constitutive theory for the solid proposed in Sec. V B is using another idea from the book of GFA, which is again formulated in terms of non-equilibrium continuum thermodynamics concepts.

A decisive step in the transformation of continuum thermodynamics into a practical tool for calculations was made by Coleman and Noll who showed how the second law can be used to derive formal state equations, given a set of state variables.43 This powerful procedure has become very popular in the continuum mechanics literature since it was proposed in 1963. However, departing in spirit and implementation from established methods in DFT, the Coleman–Noll (CN) procedure has received little attention from the physical chemistry community. Since the present paper is largely aimed at this audience, some of the key equations underlying the CN procedure are summarized below. Our main source is again GFA. A compact and accessible introduction can be found in the book by Kovetz20 who applied the CN procedure in a derivation of the Maxwell stress tensor.

We begin with the equation for the free energy of the subsystem Pt,
(37)
where ψm is the free energy per unit of mass (specific free energy density). If the system undergoes an isothermal process, the time derivative of Ψ(Pt) is constrained by the inequality
(38)
The quantities on the right-hand side are energy flows from the environment surrounding Pt into the system contained in Pt. W(Pt) is the power expended by mechanical forces,
(39)
where T is the Cauchy stress tensor and b is a volume force density accounting for long range interactions. Recall that, in addition to external forces, the traction Tn and volume force b also account for interactions of Pt within the body outside Pt (principle of Euler cuts). T and b satisfy the force balance (Cauchy) equation
(40)
where inertial effects are formally subsumed in b (see the book of GFA). T(Pt) is the energy flow due to diffusive transport
(41)
where h and h are the diffusive flux and supply term appearing in the mass balance [Eq. (11)].

The free energy imbalance [Eq. (38)], also known as the Clausius–Duhem equation, is the result of the combination of relations for part-wise energy balance and entropy production, assuming that the temperature is homogeneous and constant. Free energy imbalance plays an absolutely central role in the thermodynamics of deformable systems. The inequality is applied together with force balance equations such as the Cauchy equation [Eq. (40)], which are assumed, given, or are derived using the principle of virtual power.19,22 Here, the force balance equation is taken as a proven premise. A derivation using the virtual power balance can be found in Ref. 10.

The Coleman–Noll (CN) procedure, for our purpose, is best carried out in material variables. Recasting Eq. (37) in the referential form, we have
(42)
where ψR is the free energy per unit of reference volume. Recalling Eqs. (17) and (20), it is obtained from the specific free energy by multiplying with the reference density,
(43)
The mechanical energy flow of Eq. (39) in material formulation is written as
(44)
where TR is the first Piola stress tensor and bR is the referential transform of the volume force density. TR is obtained from the Cauchy stress tensor according to
(45)
Equation (45) is another example of a Piola transform [cf. Eq. (24)]. The transformation rule for b is simply
(46)
The counterpart of the spatial Cauchy equation [Eq. (40)] then becomes
(47)
Similarly, we write the chemical energy flux [Eq. (41)] as
(48)
with the material diffusive flux hR and supply term hR given in Eq. (24) and Eq. (25), respectively.
The next step is converting the surface integrals in Eqs. (44) and (45) to volume integrals over P. Applying the divergence theorem followed by the chain rule, the first term in Eq. (44) for the expenditure of mechanical power can be rewritten as
(49)
A: B is the notation for the matrix contraction,
(50)
Adding the second term in Eq. (44) containing bR gives
(51)
Replacing χ̇ by the time derivative of the deformation gradient using Eq. (14),
(52)
and substituting the force balance equation [Eq. (47)], we end up with
(53)
Rather unexpectedly, the volume force term due to bR cancels out.
Repeating a similar procedure for (minus) the surface term in the transport work rate [Eq. (41)], again formulated in material variables,
(54)
we see that the dependence on the contribution of the mass supply in Eq. (45) is eliminated by the equivalent term in mass balance [Eq. (27)]. The result is a chemical work rate
(55)
coupling the time derivative of density and its flux to the chemical potential and its gradient. Finally, because differentiation with respect to time can be interchanged with volume integration in material space,
(56)
the free energy imbalance [Eq. (38)] can be localized, yielding
(57)
Equation (57) is the key result in the diffusion–deformation coupling theory derived by GFA. They also present a variant for the deformed body in current space. The reference space formulation is, however, more suitable for our purpose.
The free energy imbalance [Eq. (57)] will now be applied to obtain expressions of the stress tensor and chemical potential in terms of derivatives of the reference free energy density ψR. This is the CN procedure at work. Adopting ρR and F as primitive variables, the functional form of the reference free energy density ψR defined in Eq. (43) is written as
(58)
with similar expressions for the Piola stress, chemical potential, and material mass flux,
(59)
Note that the chemical potential gradient enters as a further state variable, determining the diffusive flux. As explained by GFA, this dependence is required to satisfy the dissipation inequality. This is, however, not of direct concern here. Focusing on equilibrium quantities, the kinetics of transport will not be considered.
Expanding the time derivative of the free energy density in the velocities of the chosen state variables,
(60)
and substituting in the free energy imbalance [Eq. (57)], we obtain
(61)
CN then requires that this inequality holds not only for the actual velocities of Ḟ and ρṘ but also for all (admissible) virtual velocities. This is only possible if the affine forces (the prefactors) vanish,
(62)
(63)
In addition, the diffusive flux must satisfy a transport inequality
(64)
for all (F, ρR, ∇μ).

Equations (62) and (63) are what one would expect from a direct variational calculation carried out in reference space using ψR as free energy density. Equation (64) goes beyond what free energy minimization can produce. It is a kinetic relation reminiscent of Fick’s law for diffusion. In the present investigation, only the thermodynamic derivative relations [Eqs. (62) and (63)] are needed. For further discussion about kinetics, the reader is referred to the book of GFA.

Liquids have no memory. A Lagrangian detour via reference space should be an unnecessary exercise. Moreover, if we decide to take this route, we expect the choice of material density to be arbitrary. In this section, we verify that this is, indeed, the case. This derivation has been included not only for pedagogical reasons but also because this issue is going to be relevant later on for the constitutive definition of a simple solid at the melting line. We want the solid to resemble a liquid as much as possible in order to limit the number of order parameters to a minimum.

The free energy of a “body” of liquid in the local density approximation (LDA) is an integral over a local free energy density f(ρ),
(65)
The free energy density f(ρ) is exclusively a function of the density ρ in the spatial (current) frame. Translating in a free energy density ψR per unit volume in reference space, we once again using Eqs. (17) and (20) have
(66)
To work out the constitutive equation [Eq. (63)] for the chemical potential, we apply the chain rule. The deformation matrix F is now considered an independent degree of freedom and kept fixed,
(67)
J = det F is determined by F only, and therefore, the second term vanishes. Similarly, ∂ρ/∂ρR = J−1. Therefore, on account of Eq. (63), for the chemical potential, we have
(68)
reproducing the DFT equation for the chemical potential. Indeed, as GFA pointed out, μ is not a density and should be invariant under pull back to reference space. Equation (68) is a central identity for linking DFT and continuum mechanics.
For the evaluation of the LDA Piola stress tensor [Eq. (62)], we recall the Jacobi rule for derivatives of determinants,
(69)
which is then used to determine the tensorial deformation gradient derivative of the density
(70)
Substituting into Eq. (62)
(71)
and using Eqs. (69) and (70), we find
(72)
TR is a stress tensor in reference space. To change this to an expression for the Cauchy stress T in current space, we apply the inverse of Eq. (45) defining the Piola stress tensor TR in terms of the Cauchy stress tensor,
(73)
Applying this transformation to Eq. (72), for the Cauchy stress tensor, we find
(74)
where
(75)
is the familiar thermodynamic pressure of the simple liquid specified by the local free energy density f(ρ). Note that all explicit dependence on material quantities has been eliminated.
Returning to Sec. II C, we are now ready to deliver on the promise made there and formally verify the equivalence between chemical equilibrium and mechanical equilibrium for the LDA fluid. We will do this starting from the continuum mechanics end by substituting the constitutive relation for LDA stress [Eqs. (74) and (75)] in the force balance law [Eq. (40)]. We assume, as in Eq. (31), that the body force b can be derived from a one-particle external potential ϕ(x),
(76)
where b is balanced by the divergence of the Cauchy stress [Eq. (74)], which is evaluated using the identity
(77)
Inserting Eq. (75) yields
(78)
Substituting in the Cauchy equation with Eq. (76), we find
(79)
The next step is using the LDA constitutive relation [Eq. (68)] for the chemical potential, giving
(80)
Dividing by ρ, we have recovered the chemical equilibrium condition [Eq. (30)]. Similarly, eliminating the factor ρ in Eq. (79) and integrating results in the regular DFT equation for the equilibrium density,
(81)
where μ̄ is an integration constant to be identified with the chemical potential of the external particle reservoir. Reverting the right-hand side of Eq. (80) to a pressure gradient using the force balance gives the local formulation [Eq. (29)] of the Gibbs–Duhem relation. Finally, let us emphasize here again the more restrictive status of the local compared to the global Gibbs–Duhem relation. The local version is a constitutive property, not a thermodynamic principle.

It may seem that this elaborate continuum thermodynamics argument takes a lot of effort to reestablish a relation that was essentially derived in a couple of lines in Sec. II C. Note, however, that the treatment of the chemical potential μ is more general. μ was defined in Eq. (41) as the energy transferred in diffusive transport. This eventually led via Eq. (55) to the local free energy imbalance [Eq. (57)]. Particle transport remains active even in inhomogeneous closed systems without outside supply of particles [h = 0 in Eq. (41)]. The continuum thermodynamics definition of chemical potential is, therefore, more general than the role assigned to it in DFT. We regard this as a significant advantage of the continuum thermodynamics view, ultimately justifying the definition of chemical potential for elastic solids in Sec. VI.

How can the diffusion coupled continuum mechanics scheme for simple liquids outlined in Secs. III and IV be generalized to simple solids? Now, the mathematical machinery of non-linear continuum mechanics can no longer be avoided. The first step is replacing F by the right Cauchy–Green tensor,
(82)
Later, we will also need the left Cauchy–Green tensor,
(83)
The advantage of using C and B as state variables rather than F is that the frame indifference of the stored energy density is automatically built in. Frame indifference is not much of an issue for liquids. For elasticity of solids, it is considered crucial as discussed at length in the literature, including the main theoretical reference used in the present study.19,24 The reader is referred to these textbooks for further explanation. The constitutive relations [Eqs. (58) and (59)] will have to be reformulated with C and ρR as independent state variables,
(84)
where TRR is the second Piola stress tensor related to the first Piola stress tensor TR as
(85)
Unlike TR that changes when observers change frame, TRR is frame indifferent. With this change in state variables, Eq. (60) is modified to
(86)
which is to be substituted in the fundamental free energy imbalance relation of Eq. (57). Next, using the stress power relation (see the book of GFA)
(87)
we obtain
(88)
We are now ready to apply the Coleman–Noll procedure, as we did in Sec. III, yielding state equations for the second Piola stress tensor and chemical potential,
(89)
(90)
Using the inverse of the transformation of Eq. (85), for the (first) Piola stress tensor, we find
(91)
Equations (91) and (90) are an equivalent reformulation of the state equations [Eqs. (62) and (63)] valid for both simple liquids and solids.
This section is the mathematical pivot of our model for a simple solid. After switching from F to C, we need further refinement of the description of strain. Solids can sustain shear isochoric deformation. To explicitly account for the distinction between volumetric and isochoric deformation, we follow GFA and decompose the deformation gradient F in a volumetric factor Fv and isochoric factor Fi,
(92)
where Fv is an isotropic tensor proportional to the unit tensor 1,
(93)
and
(94)
With this definition, the determinant of the isochoric factor is unity,
(95)
Note that Eq. (95) is not a constraint but an intentional property of the kinematic variable Fi, which in this respect can be compared to a generalized coordinate in particle mechanics. For liquids in equilibrium, the multiplicative decomposition [Eq. (92)] is not necessary, and in the theory for the simple liquid in Secs. III and IV, we could continue to use the full unresolved and unconstrained deformation gradient F. The only F derived property that mattered in the end is its determinant J, leading to a fully isotropic hydrostatic stress tensor [Eq. (74)].
The decomposition of the deformation gradient Eq. (92) is reflected in a similar factorization of the Cauchy–Green tensor [Eq. (82)],
(96)
Multiplicative decomposition is applied in the continuum mechanics of anelasticity to separate plastic from elastic deformation. Equation (96) however, is not meant to describe an anelastic process. It refines the state variables describing purely elastic deformation without introducing additional degrees of freedom. The thermomechanics of Sec. III should remain valid, while the constitutive formalism of Sec. V A will have to be modified. This is the procedure applied by GFA, which is adopted here. The reference free energy density from which the constitutive relations are derived is now written as an explicit function of Ci, J and the reference density ρR,
(97)
The thermochemical principles used to derive equations of state in Sec. V A remain the same, but the mathematics is more involved. We begin by defining two auxiliary forces, which will be given a more physical interpretation later. The derivative with respect to C is split into a derivative with respect to the isochoric Cauchy–Green tensor,
(98)
and the determinant J,
(99)
where pJ will be referred to as the volumetric pressure.
As mentioned before, the constitutive theory for stress is based on Eq. (89) for the second Piola stress tensor. GFA showed that this expression can be split into two terms involving the forces of Eqs. (98) and (99),
(100)
This derivation is not straightforward. The author has not seen it presented anywhere else than in the textbook of GFA. In fact, it was a bit of a discovery finding it there. It is therefore briefly outlined in the  Appendix. Next, Eq. (100) is transformed into an expression for the first Piola stress tensor applying Eq. (85),
(101)
and from there on to the Cauchy stress tensor using Eq. (73),
(102)
GFA then continued their derivation by proving that (see the  Appendix)
(103)
The volumetric pressure pJ of Eq. (99) can therefore be identified with the actual mechanical pressure p. Introducing the notation dev A for the deviatoric component of a tensor A,
(104)
The Cauchy stress tensor can be written as
(105)
with
(106)
The corresponding expression for the chemical potential can be taken over from Eq. (90) without further manipulations,
(107)
This leaves us with the still rather mysterious tensor Πi, which will turn out to have a simple physical meaning related to the shear stress for an isotropic elastic solid, as will become clear in Sec. VI A.
A simple liquid can be turned into a simple solid by adding an elastic term eΛ to the local free energy density,
(108)
with eΛ being a function of spatial (deformed) density ρ and the right isochoric Cauchy–Green tensor defined in Eq. (96),
(109)
where Λ(ρ) is a density dependent (i.e., non-linear) shear elastic coefficient. In the limit of small deformation, Λ can be identified with the shear Lamé coefficient of linear elasticity theory, as will be verified in Sec. VI C.

Shear stress is the signature of the breakdown of mobility at the microscopic level. Failing a detailed representation of atomic position in our model, the elastic coefficient Λ plays the role of a continuum order parameter for the transition from a liquid to a simple solid. This is why Λ has been made dependent on the spatial density ρ and not on the material density ρR in the reference frame. In this way, the shear elastic response can be triggered by the modest difference in density between the liquid and solid phases. The free energy density of Eqs. (108) and (109) can therefore be thought of as defining a Landau-type free energy function for solidification. However, in a continuum mechanics framework, this distinctly non-standard construction raises a number of technical questions. In particular, shear strain is now an isochoric deformation of a reference space attached to the current deformed density. Further discussion will be deferred until after the presentation of the application in Sec. VII.

To obtain the constitutive relations for the isotropic solid, the free energy density [Eq. (108)] is subjected to a similar treatment as the liquid free energy density in Sec. IV A. First, the spatial free energy density is converted to a free energy density in material space,
(110)
For the local energy density fR, there is no difference with Eq. (66) for the liquid,
(111)
As intended, ΛR is of the same form,
(112)
with the prefactor J having the same effect (see below).
Equation (107) for the chemical potential is easy to work out. The reference density derivative of the isotropic first term has already been dealt with in Sec. IV. The ρR derivative of the second anisotropic term has an extra factor TrCi3, which is, however, independent of density. We can therefore directly copy from Eq. (68) to find
(113)
where μ̄ is the regular local density chemical potential,
(114)
and we have introduced the coefficient
(115)
characterizing the density derivative of the shear Lamé modulus. The shear elasticity adds a finite contribution to the chemical potential for strained systems but only when the elastic coefficient varies with density (κΛ ≠ 0). Where it concerns the effect on the chemical potential, a density dependent shear Lamé modulus plays the same role as a density dependent dielectric susceptibility in electrostriction, suggesting to refer to κΛ as the chemostriction coefficient.
The evaluation of the volumetric pressure of Eq. (99) involves taking a derivative to J. Inserting the model reference free energy density [Eq. (110)] gives
(116)
The partial derivative is carried out at fixed reference density ρR and shear strain Ci. Substituting Eq. (111) and applying the chain rule, we can rewrite the first J derivative in Eq. (116),
(117)
On account of Eq. (19), we can convert this into a derivative to current density ρ,
(118)
which is of the same form as what we found in Sec. IV A for the pressure of the simple liquid. Consistent with Eq. (114) for the local density chemical potential, p̄ will be referred to as the local density pressure. The J derivative of ΛR of Eq. (112) can be treated by a similar procedure, yielding
(119)
where κΛ is the chemostriction coefficient defined in Eq. (115). Combining, we obtain the expression for the total mechanical pressure pJ,
(120)
The suffix J of pJ can be dropped but will often be retained for clarity. Note that in the limit of linear elasticity (κΛ = 0), the mechanical pressure [Eq. (120)] remains sensitive to shear strain, while the chemical potential [Eq. (113)] is still determined by the local density contribution only. This contrast between chemical potential and pressure was already pointed out by Gurtin in Ref. 6.
The deviatoric Cauchy stress tensor T0 is the interesting mechanical property special to solids. This is the next quantity to evaluate. Equation (106) for T0 turns out to be rather simple for the isotropic solid [Eq. (108)]. The energy derivative is contained in the auxiliary stress-like tensor Πi defined in Eq. (98). This is a tensorial derivative with respect to the isochoric Cauchy–Green tensor Ci, where Ci appears as its trace in the energy term [Eq. (109)]. Fortunately, the derivative of the trace of a symmetric tensor is simple. Using Eq. (A2), which is valid for any symmetric tensor, we obtain
(121)
Substituting in Eq. (106) gives
(122)
Recalling the definition of the left Cauchy–Green tensor [Eq. (83)] and noting that tr Ci = tr Bi, we find that the deviatoric Cauchy stress tensor of a hyperelastic isotropic simple solid is proportional to the corresponding deviatoric component of Bi,
(123)
where we have used Eq. (112) to change back to the spatial form of the Lamé modulus.

Comparing to Eq. (74) for the liquid, we see that, in addition to a modification of the isotropic pressure term, now there is also an anisotropic term proportional to the shear Lamé modulus of the elastic model. This is, of course, a familiar result in non-linear elasticity theory (see the work of Ogden24). However, this expression is commonly given in terms of the full left Cauchy–Green tensor B. The fact that we can replace B by Bi is a convenient feature of the volumetric–isochoric decomposition of Eq. (92). The other important point to note is that the deviatoric stress of Eq. (123) and corresponding isotropic (total) pressure of Eq. (120) were derived in conjunction with a consistent expression for the chemical potential. Equation (113) makes it possible to define chemomechanical coupling coefficients in the form of the chemistry strain tensor introduced in Sec. VI B.

In thermodynamics, each independent state variable is associated with a thermodynamic force. In a hyperelastic (Green elastic) material, all thermodynamic forces can be obtained as derivatives of the stored energy density ψR in the reference frame with respect to the chosen set of state variables. This was implemented in Sec. V B for the volumetric/isochoric separation of the deformation gradient matrix. The corresponding thermomechanical response functions can therefore be formulated as second derivatives with respect to these state variables. Starting with the mechanical pressure p = pJ, this quantity, normally defined as a derivative of free energy with respect to volume, was in Eq. (99) identified with (minus) the derivative of ψR with respect to the determinant J of the deformation gradient matrix. Accordingly, we define a bulk modulus as the second derivative with respect to J,
(124)
Allowing for the variation of the reference density ρR, the chemical potential could be identified with the derivative with respect to ρR [Eq. (107)], leading GFA to define a chemistry modulus
(125)
As argued in Sec. II, reference density ρR and J are formally independent state variables in diffusion coupled deformation theory, and hence, the mixed derivative can be interpreted as the chemomechanical coupling coefficient,
(126)
The second identity is a version of the Maxwell relation.
While treated as independent variables in the thermomechanical balance equation, material density ρR and “local” volume J are coupled in the constitutive model [Eq. (108)] for the simple elastic solid. They only appear combined in the product ρ = J−1ρR [Eq. (20)]. As a result, the partial derivatives in Eqs. (124)(126) can all be expressed as derivatives to spatial density ρ. Let us first see how this works out for Σ. Changing J derivatives at constant ρR to ρ derivatives, Eq. (126) is transformed to
(127)
where the explicit specification of the variables kept constant in the partial derivative have been suppressed. Substitution of Eq. (113) for the density derivative of the chemical potential gives
(128)
with the chemostriction coefficient κΛ defined in Eq. (115). Setting Λ/ = 0, we obtain
(129)
The local density contribution μ̄ is a function of density only with a derivative determined by the compressibility. This follows from Eqs. (114) and (118),
(130)
and is yet another formulation of the local Gibbs–Duhem relation [Eq. (28)] valid for LDA functionals. Replacing the density derivative of p̄ by the compressibility
(131)
we have
(132)
Substituting in Eq. (129) yields
(133)
This expression for the chemical potential differentiated to spatial density will reappear frequently, which is why we have given it its own symbol η in the last step, where we also have inserted a reminder that it is still a partial derivative. The isochoric deformation is fixed. ρ2η can be interpreted as an effective bulk modulus corrected for the effect of chemostriction. The correction is positive since the elastic energy eΛ can be expected to be larger (or equal) to zero. The solid becomes harder to compress.
Returning to Eq. (126), Σ can now be written in terms of the coefficient η of Eq. (133). Substituting in Eq. (127), we obtain
(134)
Evidently, Σ is always negative. Next, the chemistry modulus Γ of Eq. (125) is subjected to similar manipulation with a similar result. The ρR derivative at constant J can again be converted to a derivative to spatial density ρ, and we can write
(135)
The bulk modulus of Eq. (124) is equally proportional to η. First, the J derivative at constant ρR is exchanged for a ρ derivative,
(136)
Then, substituting Eq. (120) and using Eq. (130), we have
(137)
and hence, we can write the bulk modulus of Eq. (124) as
(138)
Having converted all ρR and J partial derivatives to ρ derivatives, we can now assemble the differential for the chemical potential (still keeping Ci fixed),
(139)
Applying the same procedure to the differential of the volumetric pressure, we find
(140)
Comparing Eqs. (139) and (140), we must conclude that
(141)
This is of the same form as Eq. (130) for the local density pressure p̄ and chemical potential μ̄. The differential increments of pJ and μ again satisfy a local Gibbs–Duhem relation [Eq. (28)], now including a shear elastic contribution. The extension of the Gibbs–Duhem relation to volumetric pressure in a solid is an intentional feature of Eq. (108). It is a consequence of making the shear modulus Λ in Eq. (109) depend on spatial density ρ rather than reference density ρR.
The built in liquid-like behavior is also reflected in the chemomechanical response coefficients KJ, Γ, and Σ. Multiplying with factors of density to obtain compatible physical dimensions, we define a 2 × 2 Hessian,
(142)
Substitution of Eqs. (134), (135), and (138) yields
(143)
We immediately see that K has a zero eigenvalue. The other eigenvalue is
(144)
which is positive definite. These eigenvalues characterize the response to normal variations in ρR and J at constant Ci. Increasing ρR amounts to locally adding particles at fixed volume. Scaling up J corresponds to an increase in local volume at a fixed particle number. An increase in ρR compensated by the same increase in J conserves density ρ [Eq. (20)]. This has no effect on the energy of our model and must therefore give a zero eigenvalue. The complementary combination of increasing ρR and scaling back J leads to higher spatial density and therefore a change in energy. The eigenvalue is finite. A change in J with a value opposite of what was needed for the zero eigenvalue will double the deformed density. This explains factor 2 in Eq. (144).
Continuing with the response to a variation in isochoric deformation, we define the tensor
(145)
Substituting Eq. (113) and using Eq. (A2),
(146)
Corresponding to a first order change in chemical potential,
(147)
Adding to Eq. (139), we for the total differential of the chemical potential find
(148)
A similar derivation for the pressure response gives
(149)
Finally, we investigate the chemomechanical response in the natural equilibrium state. This state is stress- and strain-free,
(150)
The solid nature of the system becomes apparent in the finite response to shear deformation,
(151)
where the “0” superindex has been omitted from the quantities on the r.h.s by the understanding that they are to be evaluated in the natural equilibrium state. In particular, setting = 0 in Eq. (148), for a mechanical action at constant chemical potential carried out with the natural state as the initial state, we have
(152)
Equation (152) is a key result in our chemomechanical analysis. It shows that chemostriction (κΛ ≠ 0) couples density and shear strain under conditions of constant chemical potential. Shear stress can therefore induce changes in density, while the solid remains in chemical equilibrium with a liquid reservoir. This is impossible for a system under hydrostatic pressure for which density can only be changed by a change in chemical potential as is the common behavior for simple liquids.
The mechanochemical response coefficients of Sec. VI B are complemented by a set of elastic moduli formally defined as the Cauchy–Green tensor derivative of the second Piola stress tensor,
(153)
Working this out for Eq. (100) of TRR is a daunting task even for simple isotropic elastic solids. The analysis will therefore be restricted to the small strain expansion relative to the equilibrium state [Eq. (150)]. All these are necessary for the application in Sec. VII B. The model system is a body of solid material of rectangular shape with an axis frame oriented along edges. In this geometry, the left and right Cauchy–Green tensors are identical. As usual, in continuum mechanics, x, y, and z coordinates are labeled by the numerals m = 1, 2, 3. Written in the coordinate index form, the isochoric tensor [Eq. (96)] is
(154)
where λ1, λ2, λ3 are the principle stretches. Because of the requirement of unit determinant for isochoric deformation, the trace
(155)
is not a quadratic function of the principal stretches. Setting λm = 1 + ϵm and expanding in strains ϵm to the second order, we obtain
(156)
The second line is a spectral representation in terms of normalized eigenmodes, an “A” mode with eigenvalue (4/9)5 and two independent “E” modes with eigenvalue (4/9)2. There is no term linear in ϵm. The elastic energy is parabolic in ϵm with a minimum of eΛ = Λ(tr Bi − 3)/2 = 0 at ϵ1 = ϵ2 = ϵ3 = 0.
For the stress tensor, we also need the deviatoric component of Bi as defined in Eq. (104),
(157)
Now, the Taylor expansion does generate a term first order in ϵm,
(158)
Equation (158) can be used to cast the deviatoric Cauchy stress tensor of Eq. (123) in the conventional Lamé representation of linear elasticity,
(159)
The superscript i has been added as a reminder that Eq. (159) only accounts for stress due to isochoric deformation. Both μi and λi are therefore proportional to the parameter Λ introduced in Eq. (109),
(160)
confirming that Λ can indeed be interpreted as the shear modulus G (also indicated by μ). λi on the other hand is not the full Lamé λ. Indeed, evaluating the bulk modulus according to the usual textbook expression,2 
(161)
and substituting the moduli of Eq. (160) gives Ki = 0, which would indicate that the system is not stable. However, we have ignored the isotropic term in the stress tensor, as given in Eq. (105). The corresponding contribution to the bulk modulus is obtained from the small strain expansion of the volumetric pressure pJ. Using Eq. (138), we first write dpJ as a differential of J,
(162)
Next, substituting the small strain approximation to J,
(163)
we can write the full linearized Cauchy stress tensor as
(164)
which implies that the proper Lamé λ is equal to the sum ρ2η − 2Λ/3. Substituting in the standard expression Eq. (161) for the bulk modulus, the 2Λ/3 terms cancel and we are back at Eq. (138) for KJ. Note, however, that for finite chemostriction, the bulk elastic response is still affected by shear strain [see Eq. (133) for η].
Of interest is also the Young modulus E associated with the linear elastic model defined by Eq. (164). Consulting again Ref. 2, we use the relation
(165)
Similarly, choosing a for our application convenient equation for the Poisson ratio, we have
(166)
Evidently, ν < 1/2 and decreases with increasing magnitude Λ > 0 of the isochoric elastic energy. Young’s modulus of Eq. (165) vanishes for Λ = 0. This is the tensile stress response expected for a liquid. However, Poisson’s ratio as predicted by using Eq. (166) remains finite for vanishing Λ tending to ν = 1/2. A Poisson ratio of a half is the correct limiting behavior for an incompressible solid (K → ∞). However, our Λ = 0 system is compressible (finite K) but has lost shear rigidity (E = 0) as is the characteristic for a liquid. Regardless of this fundamental difference, Eq. (166) assigns the same ν = 1/2. We have no explanation for this paradox other than that Eq. (166) may not be meaningful for liquids.
Consider a sharp dividing surface separating a solid from its melt. The solid and liquid are modeled as homogeneous continua. The dividing surface at this level of modeling acts as discontinuity in the constitutive equations. The solid and the liquid are in chemical equilibrium and mechanical equilibrium, exchanging mass and volume. The chemical potential and pressure of the liquid phase are those of the simple liquid of Sec. IV,
(167)
For the solid, we use the isotropic elastic continuum model of Sec. VI. The local density chemical potential and pressure are of the same form as for the liquid,
(168)
The parameterization is different and therefore the equilibrium density at coexistence (ρsρl). In this respect, a simple liquid–solid interface is similar to a liquid–vapor interface. However, unlike Eq. (167), we no longer can set μs=μ̄s and ps=p̄s. The reason is the contribution of the elastic energy term in Eq. (109) controlled by the density dependent shear Lamé Modulus Λ. The mechano-chemical potential of the solid was given in Eq. (113) and repeated in the notation used in this section for later reference,
(169)
Similarly copying Eq. (120), the volumetric pressure in the solid is written as
(170)
where the substript J has now been omitted. Adding the deviatoric stress tensor of Eq. (123) gives the full non-linear Cauchy stress tensor of the solid,
(171)
The main question investigated in this section is how shear strain affects the discontinuity in density when crossing the dividing surface. This is not only a matter of change in the volume occupied by the solid. Chemical equilibrium maintains the same chemical potential across the interface and as a result strain can also induce a transfer of mass between the two phases. This is diffusion coupled deformation in action at a dividing surface. Setting μlμs and substituting Eqs. (167) and (169), we have
(172)
This is the chemical boundary condition. In addition, the system must also satisfy a mechanical boundary condition involving the Cauchy stress tensor at the dividing surface or, more precisely, the traction vector relative to the surface normal n,
(173)
where Tl=p̄l1 is the Cauchy stress tensor on the liquid side of the interface. Substituting Ts of Eq. (171) and rearranging, we have
(174)
The pressure pa added to the liquid reference pressure p̄l represents the normal traction force due to an externally applied load.

The liquid plays the role of a reservoir with fixed chemical potential μ̄l and pressure p̄l. These quantities are not degrees of freedom but parameters. Equations (172) and (174) are therefore a set of coupled chemomechanical equations for the equilibrium reference density ρR and elastic degrees of freedom J and Ci of the solid.

The diffusion coupled deformation theory based on the multiplicative decomposition [Eq. (96)] of the Cauchy–Green tensor is intrinsically non-linear. To avoid thermochemical inconsistencies, the chemo-mechanical surface balance equations [Eqs. (172) and (174)] were therefore derived strictly staying within the non-linear framework. With these equations established, it is now safe to expand in the strain, assuming that it is small. The deviatoric isochoric left Cauchy–Green tensor dev Bi determining in Eq. (174) the shear stress will be replaced by its linear approximation [Eq. (158)]. The elastic energy eΛ appearing in both Eqs. (172) and (174) is proportional to the trace of Bi, which is in the lowest order a quadratic function of strain, as verified in Eq. (156). Linearizing Eq. (172) in the isochoric deformation, we must conclude that to the first order, the local chemical potentials remain equal (μ̄l=μ̄s) for finite pa. Similarly, omitting from Eq. (174), the terms proportional in eΛ reduces the mechanical equilibrium condition to
(175)
where dev Bi is understood to be the linear approximation [Eq. (158)].
Surface tension is not taken into account in this first exploratory study. It can therefore be assumed that in the absence of external loads, the solid will be free of shear stress (we will return to this important point in Sec. VIII). Hence, zero pa sets Bi = 0 and what is left of Eq. (175) balances the local pressures in the two phases (p̄s=p̄l). However, a finite load not only induces shear strain but also perturbs the equilibrium solid density ρs and with it the local pressure. p̄s. As a result, p̄s and p̄l no longer cancel. With the liquid acting as a reservoir, p̄l is constant whatever the value of pa. Therefore, to first approximation, p̄sp̄l=(dps/dρs)Δρs. Changes in density are second order in the strain [Eq. (152)], and we set p̄l=p̄s in Eq. (175), which simplifies it to
(176)
As expected, the induced strain is proportional to the applied stress.
While Eq. (176) may look familiar, it is not the regular response of a stable linear isotropic elastic continuum. The strain is entirely controlled by the deviatoric Cauchy stress tensor [Eq. (159)]. However, this means, as noted in Sec. VI C, that the bulk modulus K = 0. The system is infinitely compressible. This clearly extreme mechanical behavior is, in fact, consistent with classical nucleation theory, as will be further discussed in Sec. VII C. However, some rigidity remains. It shows up, for example, in the simple geometry used in the study of diffusive creep of gels:9,44 A typical experiment consists in placing a bar of gel surrounded by the bath of solvent on a rigid substrate, which constrains stretching in the y and z direction: λ2 = λ3 = 1, or in linear approximation, ϵ1 = ϵ, ϵ2 = ϵ3 = 0 (see Fig. 2). A weight is placed on the upper surface, applying a load in the x direction. This will squeeze out some of the absorbed liquid. In what is best considered as a thought experiment, we exchange the block for a bar of simple solid and the liquid by its melt. In this one-dimensional geometry, Eq. (176) is reduced to
(177)
FIG. 2.

Constrained axial loading of a solid bar (dark gray) immersed in a bath of its melt (light gray). Solid and liquid phases are in chemomechanical equilibrium (double half arrows). The bar is supported by a rigid substrate (black) inhibiting deformation in the two and three direction (stretches λ2, λ3 are unity). A weight is placed on the upper surface, applying a stress force pa compressing the bar in the one direction leading to a strain ϵ < 0. Experiments of this type are standard for the study of the swelling of polymeric gels. (A detailed discussion can be found in the textbook by Doi.44). Here, this setup is transferred in a thought experiment to a solid–liquid interface of a simple one-component system.

FIG. 2.

Constrained axial loading of a solid bar (dark gray) immersed in a bath of its melt (light gray). Solid and liquid phases are in chemomechanical equilibrium (double half arrows). The bar is supported by a rigid substrate (black) inhibiting deformation in the two and three direction (stretches λ2, λ3 are unity). A weight is placed on the upper surface, applying a stress force pa compressing the bar in the one direction leading to a strain ϵ < 0. Experiments of this type are standard for the study of the swelling of polymeric gels. (A detailed discussion can be found in the textbook by Doi.44). Here, this setup is transferred in a thought experiment to a solid–liquid interface of a simple one-component system.

Close modal
Reading the linear approximation to devB11i from Eq. (158), we for the strain ϵa induced by the applied load pa find
(178)
The bar responds to the constrained tensile stress with a modulus
(179)
The elastic constant M is known in the seismology literature as the P-wave modulus, determining the velocity of longitudinal waves under conditions of negligible shear motion. The P-wave modulus of Eq. (179) should be compared to the value obtained for a bar without the additional non-mechanical boundary condition (see Fig. 3). The linear elastic properties are now determined by the Lamé stress tensor of Eq. (164), giving
(180)
Comparing to Eq. (179), we see that chemical equilibrium has eliminated the local density bulk modulus term 1/κρ, consistent with the vanishing bulk modulus under these conditions. Still, the modulus for axial loading remains finite and positive, suggesting that a new mechanical equilibrium can be established under loading.
FIG. 3.

Constrained axial contraction of a rectangular bar decomposed in volumetric compression followed by isochoric shear. The stretches λi, i = 1, 2, 3 for each operation are given below. Indicated are also the corresponding linear elastic moduli: M is the so-called P-wave modulus, K is the bulk modulus, and G is the shear modulus. In the linear approximation, M = K + 4G/3 [cf. Eq. (180)]. For an open system exchanging mass with a liquid reservoir, as realized in the system of Fig. 2, the bulk modulus vanishes (K = 0), eliminating the free energy cost of the volumetric compression (see Sec. VII B). The shear modulus G > 0 remains finite and therefore also M.

FIG. 3.

Constrained axial contraction of a rectangular bar decomposed in volumetric compression followed by isochoric shear. The stretches λi, i = 1, 2, 3 for each operation are given below. Indicated are also the corresponding linear elastic moduli: M is the so-called P-wave modulus, K is the bulk modulus, and G is the shear modulus. In the linear approximation, M = K + 4G/3 [cf. Eq. (180)]. For an open system exchanging mass with a liquid reservoir, as realized in the system of Fig. 2, the bulk modulus vanishes (K = 0), eliminating the free energy cost of the volumetric compression (see Sec. VII B). The shear modulus G > 0 remains finite and therefore also M.

Close modal
Similar to electrostriction, the most noticeable effect of chemostriction is a change in density. Under chemical equilibrium conditions, this change in solid density is described by using Eq. (152). This relation implies that chemostriction forces, being proportional to the elastic energy, are a non-linear effect and can therefore be ignored in linearized equilibrium equations, as in Sec. VII B. Still, for a given first order strain obtained independently of chemostriction, we can compute the second order change in density. Such a non-self-consistent decoupled approach is often applied for a first exploration of non-linear effects. The non-linear quantity of interest is the shear elastic energy eΛ. Inserting the strain of Eq. (178) into Eq. (157), we find
(181)
The elastic energy increases with the square of the perturbation as is the expected behavior in linear response theory. Evaluation Δρs of Eq. (152) for the elastic energy of Eq. (181) gives
(182)
The density always decreases whether the applied stress is compressive (pa > 0) or tensile (pa < 0). For a tensile load, a decrease in density would be the normal response because stretching at fixed lateral dimensions increases the volume of the bar. Under compression, a reduction of solid density can only be achieved by conversion of the solid material to liquid, which is an admissible degree of freedom for an open system. The change in density is accompanied by a change in volumetric pressure, which can be computed using Eq. (140),
(183)
By the substitution of Eq. (181), neglecting higher order terms quadratic in eλ, we obtain
(184)
In Sec. VI B, we argued that expressions such as Eq. (183) relating the differential of volumetric pressure and chemical potential can be regarded as a solid variant of a local Gibbs–Duhem relation. It must be kept in mind, however, that Eq. (184) is specific to our simple constitutive model. In this context, note also that a discontinuity in volumetric pressure is not a violation of mechanical equilibrium. Force balance in a planar geometry is controlled by the normal pressure. Due to elastic rigidity, these two measures of pressure can deviate, as explained in Sec. II C.

Having determined the changes in geometry, we are ready for a closer investigation of the energetics and stability. First of all, there is the issue of the zero bulk modulus. We were led to this conclusion by the linear analysis of Sec. VII B. However, the lack of resistance to volumetric expansion is an inherent property of the original non-linear constitutive model, as formulated in Eqs. (108) and (109). With equal principle stretch in all directions, the isochoric Cauchy–Green tensor Bi is unity. The trace should be 3, as confirmed by evaluating Eq. (155) for λ1 = λ2 = λ3. Similarly, Eq. (157) for the deviatoric stress is the zero tensor. The shear elastic energy is eΛ = 0, and the Cauchy stress tensor becomes Ts=p̄s1. There is no difference with the hydrostatic stress tensor Tl=p̄l1 of the liquid other than the actual density dependence of the LDA pressure p̄.

This is what the volumetric–isochoric factorization of deformation built into the constitutive model of the solid was designed for. Under load and strain-free chemomechanical equilibrium, p̄s=p̄l and μ̄s=μ̄l. The liquid/solid dividing surface can be displaced, expanding or shrinking the volume of the solid, at no free energy cost, provided that the shape is not changed. This picture is consistent with classical nucleation theory. When there is no net decrease in free energy to be gained by enlarging the volume occupied by the growing phase, the nucleus will collapse under the stress exerted by the surface tension. Surface tension is not accounted for in our model, and hence, the volume of the solid body is arbitrary. Clearly, the introduction of surface tension is badly needed to make the diffusion coupled deformation scheme proposed here more realistic. Further discussion of this point is deferred to the concluding remarks in Sec. VIII.

The case for stability under axial loading is not as rigorous. The problem is the non-self-consistent treatment of the chemostriction energy in the linearization scheme applied in Sec. VII B. Consider a change in shape preserving density. The cost of energy is the shear energy eλ, which is quadratic in the strain. This is the only contribution to the energy accounted for in Sec. VII B. However, as argued above, chemostriction perturbs the density. To the second order in the strain, the change in solid density Δρs is proportional to the elastic energy, as is indicated in Eq. (182). The corresponding contribution to the energy density can be estimated as μsΔρs and is of the same order as eΛ. The total shear stress induced free energy density change is therefore given by
(185)
For a stable solid (μs < 0), the deformation process will be exergonic (Δf > 0), assuming as before that the chemostriction coefficient κΛ > 0 (the magnitude Λ of the shear interaction increases with density). Axial loading will still come at an energy cost, but clearly, this will need to be confirmed by a rigorous self-consistent non-linear treatment.

This paper reports on a continuum mechanics investigation of the mechanical response of a solid in equilibrium with its liquid phase at the melting line. The liquid is modeled by the usual free energy density depending on local density only (no gradients). Steric repulsion and cohesion in the solid are represented by a similar local free energy density. Allowing for a different parameterization, this constitutive design makes the solid resemble the liquid as closely as possible, which is appropriate for a solid in equilibrium with its melt. The crucial difference is that the local density functional in the solid is extended with an elastic energy term accounting for shear deformation. The magnitude of the elastic term is controlled by a non-linear Lamé shear modulus varying with density. This elastic coefficient therefore effectively acts as an order parameter. A finite value turns a liquid into a solid without explicit modeling of positional ordering of atoms. Interface energy and surface tension have not been included in the present minimal model.

Chemo-mechanical equilibrium between solid and liquid is imposed by requiring that the chemical potential and normal component of the Cauchy stress tensor are continuous across the interface. The expressions for the chemical potential and stress tensor of the simple continuum solid are derived by applying diffusion coupled deformation theory. This theory was originally developed for the study of diffusion induced stress in multi-component systems, but following a suggestion by Baek and Srinivasa,8 the formalism is applied here to a one-component two phase system. Consistent with classical nucleation theory, we find that the bulk modulus of a solid body immersed in a bath of its melt vanishes. However, the solid keeps its shear rigidity. This is demonstrated by evaluating the elastic modulus for constrained axial loading of a rectangular bar. This example also illustrates how density dependence of the shear modulus gives rise to chemostriction. This somewhat unconventional name has been chosen because of the similarity to electrostriction in a dielectric fluid polarized by an external electric field. Chemostriction is manifested in a decrease in equilibrium density and pressure of the solid under mechanical loading. Similar to electrostriction, the effect is quadratic in the applied force and therefore the same for tensile dilation and compression.

The application is the subject of the final section with the larger part of this paper taken up by a detailed and inevitably technical explanation of the methodology. The reason that this was felt appropriate is that the diffusion coupled solid mechanics and constitutive theory as applied here have a number of non-standard aspects. There is first of all the question of the background against which diffusive motion is measured. In the theory of solvent permeation of gels, the reference is provided by the polymeric network. In Sec. II, we proposed that in a one-component system, the convective flow as defined in continuum mechanics can be used as a reference for migration of mass. Changing over to a Lagrangian picture, this argument leads us to treat material mass density in the reference space as an independent degree of freedom in addition to the deformation. The diffusion coupled deformation theory as presented by Gurtin, Fried, and Anand (GFA)19 is then applied without modification, leading to expressions for the chemical potential and stress tensor in terms of derivatives of the free energy density.

The constitutive model of the solid also makes non-standard assumptions about reference. This concerns the reference for shear deformation. The elastic coefficients in the model, including the non-linear shear modulus, are functions of the spatial (deformed) density. Similar to a liquid, specification of a material space is not needed. Current space is the (constitutive) material space. To implement this idea for the shear motion, the deformation gradient tensor was factorized in a volumetric expansion and an isochoric deformation following a mathematical scheme again taken over from the book of GFA. Shear deformation is effectively defined as displacive motion of covectors tangent to surfaces of constant volume density as quantified by J. However, isochoric deformation, while preserving local volume, can change inhomogeneous density. The example in the present application is Eq. (152), linking the change in density of the solid at constant chemical potential to the energy stored in the shear strain. As a result, shear deformation can modify the shear modulus, creating a self-consistency problem. This complication was evaded in the present application by linearization of the final mechanical equilibrium equations. The coupling to density was shown to be a second order effect and can be ignored for in the determination of the first order shear. However, in future fully non-linear applications (see below), these inconsistencies can no longer be ignored. Hopefully, we can learn from anelasticity theory where coupling of deformation to an evolving structure in material space is a major topic of research. The relatively “accessible” example is thermo-elasticity.19 Thermal expansion leads to another stress-free reference state.

The thermomechanical formalism and non-linear constitutive model developed in this paper are applied in a thought experiment, probing the mechanical response of a block of solids in chemomechanical equilibrium with its liquid phase. We first verify that enforcing equilibrium mass exchange across the dividing surface sets the bulk modulus of the solid body to zero. The application then continues showing how the mechanical response is stabilized by the remaining shear rigidity. This is illustrated by evaluating the elastic modulus for constrained axial loading of a solid bar in the bath of its melt. This analysis is carried out by linearizing the constitutive elastic model of the solid. The resulting strain is then used to estimate the chemostriction due to loading. The chemostriction is a consequence of the density dependence of the shear modulus and is a quadratic effect.

Clearly, our elementary model system is far removed from a phase equilibrium of a realistic solid–liquid system, such as ice cubes in a beaker of water. First of all, what is missing dearly is surface energy and surface stress. Moreover, for the description of the shape of a crystal, the surface energy will also have to be orientation dependent. This level of modeling is required to address questions such as the generalization of the Young–Laplace equation to solid interfaces.50 Furthermore, the migration of realistic solid–liquid phase boundaries can be expected to be hindered by kinetic barriers, such as activated accretion. However, we believe that the key message of the thermomechanical exercise reported in this paper still stands, namely that chemical and mechanical equilibrium must be treated separately at interfaces.

A simple elastic model, extended with surface tension, could already be useful for the investigation of soft capillarity. Similar to externally applied surface loads, surface tension will equally deform the interior of a soft solid. This is the subject of the new and fascinating field of elasto-capillarity.45–47 Solid nucleation will be a more ambitious goal. Again, the question can be asked how the capillary stress is modified by chemical equilibrium at the dividing surface. What will be the effect on the nucleation barriers? Elastocapillarity is not part of the analysis in most microscopic numerical studies of solid nucleation48–51 and freezing of confined nanophases.52 It could be important because at the melting line, solids are soft. The study of nucleation will also force us to reconsider another aspect of sharp interface modeling applied in this paper. The boundary between solid and liquid, while flexible in the deformed system, is fixed in material space, defining the reference volume and shape of the solid body. This issue touches on the fundamental question about the relation between thermodynamic tension and elastic stress.53 

We conclude with a very brief comment on the relation to the DFT of inhomogeneous fluids. Density gradients can generate deviatoric stress. The well-known example is the Korteweg capillary stress tensor.54 This tensor can be derived from the square gradient (Cahn–Hilliard) free energy density by a Lagrangian procedure similar to the one applied in this paper. The free energy density is pulled back to a reference space and differentiated with respect to the deformation gradient, yielding a Piola stress tensor. Pushed forward to current space, the result is identical to the Korteweg stress tensor.55 There is finite stress but no strain (F = 1) as is also consistent with the freedom of choice of reference space. This raises the question that can stress in an inhomogeneous fluid, in principle, generate spontaneous (“eigen”) strain, given a suitably accurate density functional? In other words, can a variational minimization of a pure DFT functional extended with a search for the optimal shear deformation tensor break symmetry finding F1? If it does, is this point where the fluid is no longer a fluid but a solid? This possibility is suggested by the observation that the higher order gradient functionals used in phase field crystal theory31–35 exhibit an elastic response. The same is true for the non-local density functionals used in the DFT of crystals.32,56,57 Unfortunately, evaluation of the stress tensor of a pattern forming (Hohenberg–Swift) density functional33 is very demanding and probably not very practical, at least not without further coarse graining.29 Here, a hybrid strain–density functional could be an option. The strain could play a similar role to the one-electron orbitals in the Kohn–Sham approach to the DFT of the electronic structure.

Similar questions concerning spontaneous strain in formally liquid systems arise in electric double layers modeled by Poisson–Boltzmann theory. The Maxwell stress tensor, generated by the electrostatic potential gradient, also comes, in general, with a deviatoric component. Can highly charged double layers develop shear elasticity and show chemo-mechanical response resembling the solid bar investigated in the work presented here? Another matter is possible coupling between Maxwell and steric stress. This requires a distinction between charge and mass density, which is usually not explicitly made in modified Poisson–Boltzmann theory.58 How are these electromechanical effects reflected in the differential capacitance profile? As mentioned in the Introduction, these questions, which are beginning to be explored,59,60 were the original motivation for the present investigation and are already the subject of ongoing research.23 Furthermore, a continuum mechanics based analysis of stress may also hopefully be useful for the interpretation of the results of surface force balance experiments carried out in electrically active soft matter.61 

Daan Frenkel is acknowledged for his interest and critical comments. No external funding was received for the research reported in this paper.

There are no conflicts of interest to disclose.

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

Listed are a selection of the equations needed in the derivation of Eq. (105) for the Cauchy stress tensor under multiplicative decomposition [Eq. (92)] of the deformation gradient tensor F. Detailed full derivations can be found in Sec. 55 of the book of GFA. Equation numbers in italic refer to the equations in that section. First, we give the C (Cauchy–Green) derivatives of three important functions of C: the determinant of F,
(A1)
of the trace of C,
(A2)
and of the isochoric factor Ci Eq. (96),
(A3)
where I is the fourth-order identity tensor. These identities are used in the derivation of Eq. (100) for the second Piola stress tensor. In the book of GFA, this is Eq. (55.7), which is written out in full as follows:
(A4)
All derivatives are at constant material density ρR, which has been suppressed as an argument of ψ̃R. Substitution of the short-hand notations of Eqs. (98) and (99) yields Eq. (100).
The multiplicative decomposition is given a particularly useful physical meaning by Eq. (103) for the trace of the Cauchy stress tensor Eq. (102) [Eq. (55.9)]. The proof proceeds by showing that the trace of T0 of Eq. (106) vanishes. Following gain from the book of GFA, we now write, suppressing all mentioning of arguments to ψ̃R,
(A5)
confirming that T0 is the deviatoric component of the Cauchy stress tensor [Eq. (102)].
1.
J.-P.
Hansen
and
I. R.
McDonald
,
Theory of Simple Liquids
, 4th ed. (
Academic Press
,
Oxford
,
2013
).
2.
M. V.
Lubarda
and
V. A.
Lubarda
,
Intermediate Solid Mechanics
(
Cambridge University Press
,
Cambridge
,
2020
).
3.
R.
Evans
, “
The nature of the liquid-vapour interface and other topics in the statistical mechanics of nonuniform, classical fluids
,”
Adv. Phys.
28
,
143
200
(
1979
).
4.
J. S.
Rowlinson
and
B.
Widom
,
Molecular Theory of Capillarity
(
Dover Edition
,
New York
,
2002
).
5.
F. C.
Larché
and
J. W.
Cahn
, “
The interaction of composition and stress in crystalline solids
,”
Acta Metall.
33
,
331
357
(
1985
).
6.
M. E.
Gurtin
, “
Generalized Ginzburg-Landau and Cahn-Hilliard equations based on a microforce balance
,”
Physica D
92
,
178
192
(
1996
).
7.
S.
Shi
,
J.
Markmann
, and
J.
Weissmüller
, “
Verifying Larché-Cahn elasticity, a milestone of 20th-century thermodynamics
,”
Proc. Natl. Acad. Sci. U. S. A.
115
,
10914
10919
(
2018
).
8.
S.
Baek
and
A. R.
Srinivasa
, “
Diffusion of a fluid through an elastic solid undergoing large deformation
,”
Int. J. Non-Linear Mech.
39
,
201
218
(
2004
).
9.
W.
Hong
,
X.
Zhao
,
J.
Zhou
, and
Z.
Suo
, “
A theory of coupled diffusion and large deformation in polymeric gels
,”
J. Mech. Phys. Solids
56
,
1779
1793
(
2008
).
10.
F. P.
Duda
,
A. C.
Souza
, and
E.
Fried
, “
A theory for species migration in a finitely strained solid with application to polymer network swelling
,”
J. Mech. Phys. Solids
58
,
515
529
(
2010
).
11.
S. A.
Chester
and
L.
Anand
, “
A coupled theory of fluid permeation and large deformations for elastomeric materials
,”
J. Mech. Phys. Solids
58
,
1879
1906
(
2010
).
12.
L.
Anand
, “
A Cahn-Hilliard-type theory for species diffusion coupled with large elastic-plastic deformations
,”
J. Mech. Phys. Solids
60
,
1983
2002
(
2012
).
13.
S.
Baek
and
T. J.
Pence
, “
Inhomogeneous deformation of elastomer gels in equilibrium under saturated and unsaturated conditions
,”
J. Mech. Phys. Solids
59
,
561
582
(
2011
).
14.
G. J.
Templet
and
D. J.
Steigmann
, “
On the theory of diffusion and swelling in finitely deforming elastomers
,”
Math. Mech. Complex Syst.
1
,
105
128
(
2013
).
15.
F.
Yang
, “
Interaction between diffusion and chemical stresses
,”
Mater. Sci. Eng., A
409
,
153
159
(
2005
).
16.
Z.
Cui
,
F.
Gao
, and
J.
Qu
, “
A finite deformation stress dependent chemical potential and its application to lithium ion batteries
,”
J. Mech. Phys. Solids
60
,
1280
1295
(
2012
).
17.
C. V.
Di Leo
,
E.
Rejovitzky
, and
L.
Anand
, “
A Cahn-Hilliard type phase-field theory for species diffusion coupled with large elastic deformations: Application to phase-separating Li-ion electrode materials
,”
J. Mech. Phys. Solids
70
,
1
29
(
2014
).
18.
F. Q.
Yang
,
Y.
Li
,
B. L.
Zheng
, and
K.
Zhang
, “
Interaction between stress and diffusion in lithium-ion batteries: Analysis of diffusion-induced buckling of nanowires
,” in
Handbook of Mechanics and Materials
, edited by
C. H.
Hsueh
(
Springer Nature Singapore
,
2018
), pp.
1
20
.
19.
M. E.
Gurtin
,
E.
Fried
, and
L.
Anand
,
The Mechanics and Thermodynamics of Continua
(
Cambridge University Press
,
Cambridge
,
2010
).
20.
A.
Kovetz
,
Electromagnetic Theory
(
Oxford University Press
,
Oxford
,
2000
).
21.
M.
Doi
, “
Onsager principle in polymer dynamics
,”
Prog. Polym. Sci.
112
,
101339
(
2021
).
22.
P.
Podio-Guidugli
, “
A virtual power format for thermomechanics
,”
Continuum Mech. Thermodyn.
20
,
479
487
(
2009
).
23.
M.
Sprik
, “
Electric-field-based Poisson-Boltzmann theory: Treating mobile charge as polarization
,”
Phys. Rev. E
103
,
022803
(
2021
).
24.
R. W.
Ogden
,
Non-Linear Elastic Deformations
(
Oxford University Press, Dover Publications
,
1997
).
25.
M.
Sprik
, “
Continuum model of the simple dielectric fluid: Consistency between density based and continuum mechanics methods
,”
Mol. Phys.
119
,
e1887950
(
2021
).
26.
P. C.
Martin
,
O.
Parodi
, and
P. S.
Pershan
, “
Unified hydrodynamic theory for crystals, liquid crystals and normal fluids
,”
Phys. Rev. A
6
,
2401
2420
(
1972
).
27.
P. D.
Fleming
III
and
C.
Cohen
, “
Hydrodynamics of solids
,”
Phys. Rev. B
13
,
500
516
(
1976
).
28.
P. M.
Chaikin
and
T. C.
Lubensky
,
Principles of Condensed Matter Physics
(
Cambridge University Press
,
Cambridge
,
1995
).
29.
A.
Skaugen
,
L.
Angheluta
, and
J.
Viñals
, “
Dislocation dynamics and crystal plasticity in the phase-field crystal model
,”
Phys. Rev. B
97
,
054113
(
2018
).
30.
A.
Acharya
and
J.
Viñals
, “
Field dislocation mechanics and phase field crystal models
,”
Phys. Rev. B
102
,
064109
(
2020
).
31.
K. R.
Elder
and
M.
Grant
, “
Modeling elastic and plastic deformations in nonequilibrium processing using phase field crystals
,”
Phys. Rev. E
70
,
051605
(
2004
).
32.
K. R.
Elder
,
N.
Provatas
,
J.
Berry
,
P.
Stefanovic
, and
M.
Grant
, “
Phase-field crystal modeling and classical density functional theory of freezing
,”
Phys. Rev. B
75
,
064107
(
2007
).
33.
H.
Emmerich
,
H.
Löwen
,
R.
Wittkowski
,
T.
Gruhn
,
G. I.
Tóth
,
G.
Tegze
, and
L.
Gránásy
, “
Phase-field-crystal models for condensed matter dynamics on atomic length and diffusive time scales: An overview
,”
Adv. Phys.
61
,
665
743
(
2012
).
34.
N.
Pisutha-Arnond
,
V. W. L.
Chan
,
K. R.
Elder
, and
K.
Thornton
, “
Calculations of isothermal elastic constants in the phase-field crystal model
,”
Phys. Rev. B
87
,
014103
(
2013
).
35.
M.
Salvalaglio
,
A.
Voigt
, and
K. R.
Elder
, “
Closing the gap between atomic-scale lattice deformations and continuum elasticity
,”
npj Comput. Mater.
5
,
48
(
2019
).
36.
L. A.
Pozhar
,
K. E.
Gubbins
, and
J. K.
Percus
, “
Generalized compressibility equation for inhomogeneous fluids at equilibrium
,”
Phys. Rev. E
48
,
1819
1822
(
1993
).
37.
V.
Romero-Rochín
and
J. K.
Percus
, “
Stress tensor of liquid-vapor states of inhomogeneous fluids
,”
Phys. Rev. E
53
,
5130
5136
(
1996
).
38.
R.
Lovett
and
M.
Baus
, “
Two molecular scale force distributions associated with a planar interface
,”
J. Chem. Phys.
97
,
8596
8605
(
1992
).
39.
R.
Lovett
and
M.
Baus
, “
The thermodynamic forces in an interface
,”
Adv. Chem. Phys.
102
,
1
38
(
1997
).
40.
C.
Walz
and
M.
Fuchs
, “
Displacement field and elastic constants in nonideal crystals
,”
Phys. Rev. B
81
,
134110
(
2010
).
41.
J. M.
Häring
,
C.
Walz
,
G.
Szamel
, and
M.
Fuchs
, “
Coarse-grained density and compressibility of nonideal crystals: General theory and an application to cluster crystals
,”
Phys. Rev. B
92
,
184103
(
2015
).
42.
A.
DiCarlo
and
P.
Podio-Guidugli
, “
From point particles to body points
,”
Math. Eng.
4
,
1
29
(
2022
).
43.
B. D.
Coleman
and
W.
Noll
, “
The thermodynamics of elastic materials with heat conduction and viscosity
,”
Arch. Ration. Mech. Anal.
13
,
167
178
(
1963
).
44.
M.
Doi
,
Soft Matter Physics
(
Oxford University Press
,
Oxford
,
2013
).
45.
J. H.
Weijs
,
B.
Andreotti
, and
J. H.
Snoeijer
, “
Elasto-capillarity at the nanoscale: On the coupling between elasticity and surface energy in soft solids
,”
Soft Matter
9
,
8494
8503
(
2013
).
46.
R. W.
Style
,
A.
Jagota
,
C.-Y.
Hui
, and
E. R.
Dufresne
, “
Elastocapillarity: Surface tension and the mechanics of soft solids
,”
Annu. Rev. Condens. Matter Phys.
8
,
99
118
(
2017
).
47.
R. W.
Style
and
Q.
Xu
, “
The mechanical equilibrium of soft solids with surface elasticity
,”
Soft Matter
14
,
4569
4576
(
2018
).
48.
A.
Cacciuto
,
S.
Auer
, and
D.
Frenkel
, “
Breakdown of classical nucleation theory near isostructural phase transitions
,”
Phys. Rev. Lett.
93
,
166105
(
2004
).
49.
B.
Cheng
and
M.
Ceriotti
, “
Bridging the gap between atomistic and macroscopic models of homogeneous nucleation
,”
J. Chem. Phys.
146
,
034106
(
2017
).
50.
P.
Montero de Hijes
,
K.
Shi
,
E. G.
Noya
,
E. E.
Santiso
,
K. E.
Gubbins
,
E.
Sanz
, and
C.
Vega
, “
The Young-Laplace equation for a solid-liquid interfaces
,”
J. Chem. Phys.
153
,
191102
(
2020
).
51.
J. F.
Lutsko
, “
How crystals form: A theory of nucleation pathways
,”
Sci. Adv.
5
,
eaav7399
(
2019
).
52.
K. E.
Gubbins
,
Y.
Long
, and
M.
Śliwinska-Bartkowiak
, “
Thermodynamics of confined nano-phases
,”
J. Chem. Thermodyn.
74
,
169
183
(
2014
).
53.
E.
Fried
and
M. E.
Gurtin
, “
Coherent solid-state phase transitions with atomic diffusion: A thermomechanical treatment
,”
J. Stat. Phys.
95
,
1361
1427
(
1999
).
54.
D. M.
Anderson
,
G. B.
McFadden
, and
A. A.
Wheeler
, “
Diffuse interface methods in fluid mechanics
,”
Annu. Rev. Fluid Mech.
30
,
139
165
(
1998
).
55.
V. A.
Eremeyev
and
H.
Altenbach
, “
Equilibrium of a second-gradient fluid and an elastic solid with surface stresses
,”
Meccanica
49
,
2635
2643
(
2014
).
56.
M.
Oettel
,
S.
Görig
,
A.
Härtel
,
H.
Löwen
,
M.
Radu
, and
T.
Schilling
, “
Free energies, vacancy concentrations, and density distribution anisotropies in hard-sphere crystals: A combined density functional and simulation study
,”
Phys. Rev. E
82
,
051404
(
2010
).
57.
J. F.
Lutsko
and
C.
Schoonen
, “
Classical density-functional theory applied to the solid state
,”
Phys. Rev. E
102
,
062136
(
2020
).
58.
M. Z.
Bazant
,
B. D.
Storey
, and
A. A.
Kornyshev
, “
Double layer in ionic liquids: Overscreening versus crowding
,”
Phys. Rev. Lett.
106
,
046102
(
2011
).
59.
A. C.
Maggs
and
R.
Podgornik
, “
General theory of asymmetric steric interactions in electrostatic double layers
,”
Soft Matter
12
,
1219
1229
(
2016
).
60.
J. P.
de Souza
and
M. Z.
Bazant
, “
Continuum theory of electrostatic correlations at charged surfaces
,”
J. Phys. Chem. C
124
,
11414
11421
(
2020
).
61.
C. S.
Perez-Martinez
and
S.
Perkin
, “
Surface forces generated by the action of electric fields across liquid films
,”
Soft Matter
15
,
4255
4265
(
2019
).
Published open access through an agreement with University of Cambridge