Derive the energy.
Find the driving force.

Why an eigenstrain produces elastic energy.

Cut an inclusion-shaped region from a stress-free matrix and let it transform freely by \(\varepsilon^0\). It no longer fits its original hole. Restoring it to the cavity requires elastic loading; welding it to the matrix enforces coherence; releasing the auxiliary tractions lets the joined body relax to mechanical equilibrium. The residual elastic strain is therefore \(\varepsilon^{\mathrm{el}}=\varepsilon-\varepsilon^0\), not the transformation strain itself.

Eshelby solves this construction for an ellipsoidal inclusion. Mura’s eigenstrain formalism generalizes the source of internal stress. Khachaturyan replaces a single idealized inclusion by shape functions or continuous fields describing arbitrary numbers, shapes and variants, and evaluates their collective relaxation in Fourier space.

Conventions before calculation.

The derivation assumes linear elasticity, homogeneous stiffness \(C_{ijkl}\), periodic fields and zero applied mean stress. Let \(\boldsymbol c=(c_1,\ldots,c_{N_c})\) collect the composition fields and \(\boldsymbol\eta=(\eta_1,\ldots,\eta_{N_\eta})\) collect the structural order parameters. We expand their general eigenstrain map in \(P\) tensorial contributions:

Each \(\theta_p\) may depend on several composition and order-parameter fields. The derivation below acts on the \(P\) scalar interpolation fields and therefore does not require one particular phase-field parameterization. The mean strain relaxes to \(\langle\varepsilon^0_{ij}\rangle\). The elastic energy is therefore determined by the fluctuating interpolation fields

Our Fourier-transform convention is

The derivation is written for spatial dimension \(d\), so the same notation covers two- and three-dimensional implementations. A prime on a Fourier integral, \(\int^{\prime}\), excludes the zero mode \(\mathbf k=\mathbf 0\).

References for the derivation: J. D. Eshelby’s 1957 ellipsoidal-inclusion solution; T. Mura’s eigenstrain formulation in Micromechanics of Defects in Solids; and A. G. Khachaturyan, Theory of Structural Transformations in Solids, Chapter 7, especially §7.2, “Strain Energy of Multiphase Alloys,” and §7.3, “Strain-Induced Interactions between Coherent New Phase Inclusions.” A detailed multivariant development appears in Saswata Bhattacharyya, Evolution of Multivariant Microstructures with Anisotropic Misfit: A Phase Field Study, Appendices A and B. The homogeneous-stiffness derivation on this page is distinct from the later inhomogeneous spectral iteration developed in our CICP paper.

From mechanical equilibrium to \(B_{pq}\).

Step 1: Mean and periodic strain.

The mean terms cancel because the mean strain can relax. This cancellation would change under imposed mean strain or other loading conditions.

Step 2: Solve the equilibrium equation mode by mode.

For k ≠ 0, define n = k/|k|, σ⁽ᵖ⁾ = C:ε⁽ᵖ⁾ and the acoustic tensor Q. Let a⁽ᵖ⁾ = σ⁽ᵖ⁾n. Mechanical equilibrium eliminates displacement through the inverse acoustic tensor Ω = Q⁻¹.

With the Fourier convention above, differentiating displacement multiplies its transform by +ik. The zero mode is handled by the mean-strain condition, not by inverting Q at k = 0.

Step 3: Substitute the relaxed strain into the energy.

The first term is the unrelaxed eigenstrain contribution. The second subtracts the energy released by a compatible displacement field. For a stable stiffness tensor, the relaxed quadratic form is nonnegative. Bₚq = Bqₚ and Bₚq(n) = Bₚq(−n).

A full derivation obtains this expression using Parseval’s theorem and cancellation of the mean terms.

The variation of both factors.

Varying the quadratic energy produces two terms: one from θ̃ₚ and one from θ̃*q. Real fields obey θ̃*(k) = θ̃(−k). Combining this reality condition with the evenness and symmetry of B makes the two half-contributions equal.

The factor ½ in the energy therefore disappears from its derivative. Here i labels the selected field; p is summed. We write the symmetry explicitly as Bₚᵢ = Bᵢₚ to avoid confusing the free field index with a summation index.

Derivatives for each composition and order-parameter field.

The kernel first gives \(g_p=\delta F_e/\delta\theta_p\). The constitutive eigenstrain map then supplies the chain rule for every underlying field:

For the specialization \(\theta_1=c\) and \(\theta_{b+1}=\eta_b^2\), these reduce to \(g_1/N_v\) and \(2\eta_b g_{b+1}/N_v\). Nᵥ is optional when the free-energy functional is already per unit volume. Evaluate the derivatives of \(\theta_p\) and multiply by \(g_p\) in real space; multiplying their Fourier transforms pointwise would incorrectly replace a product by a convolution.

The free-energy derivative and the evolution law.

Each composition derivative enters a conserved evolution law, while each structural order-parameter derivative enters a nonconserved evolution law. With cross mobilities allowed,

The expressions above provide only the elastic parts of these derivatives. Chemical and gradient-energy contributions must also be included.

The same calculation in both languages.

Assume \(P\) real interpolation fields and a precomputed real, even, symmetric \(P\times P\) kernel with \(B(\mathbf k=0)=0\). These snippets compute \(g_p=\delta F_e/\delta\theta_p\); they do not choose the constitutive map \(\theta_p(\boldsymbol c,\boldsymbol\eta)\), construct the kernel or advance time. NumPy uses axes (field, x, y); Julia uses (x, y, field).

# theta.shape = (P, nx, ny); B.shape = (P, P, nx, ny)
# Build theta from every composition and order-parameter field.
h = np.fft.fft2(theta, axes=(-2, -1))
g_hat = np.einsum('ipxy,pxy->ixy', B, h)
g = np.fft.ifft2(g_hat, axes=(-2, -1)).real
# Apply the Jacobian of theta(c, eta, ...) in real space:
# dFe_dc[a]   = sum_p g[p] * dtheta_p_dc[a] / Nv
# dFe_deta[b] = sum_p g[p] * dtheta_p_deta[b] / Nv
# Elastic energy, with NumPy's default FFT normalization:
N = theta.shape[-2] * theta.shape[-1]
Fe = 0.5 * dx * dy / N * np.sum((h.conj() * g_hat).real)

# theta: (nx,ny,P); B: (nx,ny,P,P)
# Build theta from every composition and order-parameter field.
h = fft(theta, (1, 2))
g_hat = zeros(ComplexF64, size(h))
P = size(theta, 3)
for i in 1:P, p in 1:P
    g_hat[:,:,i] .+= B[:,:,i,p] .* h[:,:,p]
end
g = real.(ifft(g_hat, (1, 2)))
# Apply the Jacobian of theta(c, eta, ...) pointwise in real space.
N = size(theta,1) * size(theta,2)
Fe = 0.5 * dx * dy / N * sum(real.(conj.(h) .* g_hat))

My numerical check of the energy derivative.

For a small perturbation hᵢ(r), the central difference [Fₑ(θ + ah) − Fₑ(θ − ah)]/(2a) should approach ∫ ∑ᵢ gᵢhᵢ dV. To check the η derivative, perturb η and rebuild θᵢ = ηᵢ² before evaluating each energy. This catches a missing factor of two, an FFT normalization error, or an incorrect field-index contraction.

Run the single-field elastic example ↗

The worked spinodal example has one composition field and a scalar kernel. The general formulas allow any number of coupled composition and order-parameter fields; the four-field model is only one specialization. The single-field example isolates the same numerical idea.