Solution (source code)

= Solution
{c}

Let $\epsilon$ here be a constant infinitesimal transformation parameter, unrelated to the <dimensional regularization> parameter in question 1. The matter and Maxwell terms undergo an ordinary infinitesimal <gauge transformation> with $\alpha=\epsilon\phi$, and their variations cancel as before. The <gauge fixing> contribution is
$$
\delta\mathcal L_{\mathrm{gf}}=-\frac\epsilon\xi B\partial^2\phi.
$$
The extra real <scalar field> kinetic term varies by
$$
\delta\left(\frac12\partial_\mu\phi\partial^\mu\phi\right)=\partial_\mu\phi\partial^\mu\delta\phi=-\frac\epsilon\xi(\partial_\mu\phi)(\partial^\mu B).
$$
Adding the two gives a total derivative:
$$
\boxed{\delta\mathcal L=-\frac\epsilon\xi\partial_\mu(B\partial^\mu\phi),\qquad\delta S=0.}
$$
The last equality uses vanishing boundary contributions, for example sufficiently fast decay of the fields. No field equations were used. This <gauge-scalar global symmetry of covariantly fixed QED> leaves the action invariant even though the density need not be invariant pointwise; the scalar's variation supplies exactly the missing boundary-term cancellation.