Solution (source code)

= Solution

Write $D=\partial_\mu A^\mu$. The additional <gauge fixing> variation is $\delta(-D^2/2)=-D\partial_\mu\delta A^\mu$. <Integration by parts> gives the <Euler-Lagrange field equation> $\partial_\mu F^{\mu\nu}+\partial^\nu D=0$, hence $\Box A^\nu=0$. This equation follows from the gauge-fixed density without imposing $D=0$ as a separate operator identity.

Direct differentiation of the density as printed gives the <canonical momentum> components conjugate to the lower-index fields:
$$
\boxed{\pi^\nu=\frac{\partial\mathcal L}{\partial\dot A_\nu}=-F^{0\nu}-\eta^{0\nu}D,\qquad \pi^0=-D,\quad\pi^i=\dot A_i-\partial_iA_0.}
$$
There is a boundary-term convention to reconcile with part (c). The <Feynman-gauge Maxwell kinetic density after a boundary-term subtraction> is
$$
\mathcal L'=-\frac12\partial_\mu A_\nu\partial^\mu A^\nu,
\qquad
\mathcal L=\mathcal L'+\partial_\mu K^\mu,
\qquad
K^\mu=\frac12\bigl(A_\nu\partial^\nu A^\mu-A^\mu D\bigr).
$$
It has the same <Euler-Lagrange field equations>, but its <canonical momentum> components are $\pi'^\nu=-\dot A^\nu$. The mode expansion supplied in part (c) is the expansion of these latter momenta. Thus it is valid after the stated boundary-term subtraction; it is not the direct derivative of the original density. For example, a time-independent spatially varying $A_0$ with $A_i=0$ gives $\pi^i=-\partial_iA_0$ but $\pi'^i=0$.

The <boundary-term shift of canonical field momenta> is a <canonical transformation>. If spatial boundary terms vanish, $B=\int d^3x\,K^0=\int d^3x\,A_0\partial_iA_i$, and $\pi^\nu=\pi'^\nu+\delta B/\delta A_\nu$. This explains why the two conventions have the same dynamics while their <momentum> formulas differ.