Solution (source code)

= Solution

Choose zero energy in the unstrained state. The required <strain energy density> is
$$
\boxed{W(e)=\frac\lambda2(\operatorname{tr}e)^2+\mu e:e.}
$$
Its <variation> is $\delta W=(\lambda\operatorname{tr}(e)I+2\mu e):\delta e$, giving the prescribed <stress> as the <tensor> <gradient>. This variational notation handles the symmetric off-diagonal components without a spurious factor of two.

Decompose $e=e^{\rm dev}+(\operatorname{tr}e)I/3$ into <orthogonal> deviatoric and spherical parts. Then
$$
\boxed{W=\mu\,e^{\rm dev}:e^{\rm dev}+\frac{3\lambda+2\mu}{6}(\operatorname{tr}e)^2.}
$$
The two squares vanish together only for $e=0$. Thus $\mu>0$ and $3\lambda+2\mu>0$ are sufficient for <positive definiteness>. They are necessary by testing a nonzero <trace>-free <strain> and a pure <dilation> separately. Equivalently both <shear modulus> $\mu$ and <bulk modulus> $K=\lambda+2\mu/3$ must be positive.