Do primary first class constraints change the electric field in the Hamiltonian form of Maxwell's theory? In my understanding of Dirac's theory of constrained Hamiltonians, the primary (and also the secondary) first class constraints are generators of canonical transformations that do not change the physical state: the electric field is part of the physical state so it has zero response to a primary first class constraint. However, a paper http://arxiv.org/abs/1310.2756 recently appeared which says that the primary first class constraints change the physical state. The paper gives a direct calculation which I'll reproduce below.
Using the notation in Dirac's Lectures on Quantum Mechanics, the p's are $B^{\mu}$ and the q's are the electromagnetic potentials $A_{\mu}$. The primary first class constraints are $B^{0}\approx 0$. So, the generator of the primary first class constraints is,
$$
G=\int d^{3}x \xi(x) B^{0}(x)
$$ 
The response of the electromagnetic field is given by the PB,
$$
\frac{dA_{\mu}}{d\epsilon}=[A_{\mu},G]=\delta^{0}_{\mu}\xi(x)\ .
$$
The paper defines the electric field as,
$$
E_{r}=A_{r,0}-A_{0,r}
$$
and denies any relation between $E_{r}$ and the canonical momenta $B^{r}$ until the dynamical equation $\dot{q}=[q,H]$ has been used. The paper gets the response of the electric field to the primary first class constraint as,
$$
\frac{dE_{r}}{d\epsilon}=\frac{\partial}{\partial t}\frac{dA_{r}}{d\epsilon}-\frac{\partial}{\partial x^{r}}\frac{dA_{0}}{d\epsilon}=-\xi_{,r}
$$
and this is troubling me because the response should be zero. 
I thought I understood constrained Hamiltonians but now I'm not sure, please help. 
 A: The problem lies in what we learn about good old constrained dynamics from traditional Dirac approach is not complete and is somehow inconsistent, and the above is one example of this. This was the message of Pitts' paper mentioned in the question above, who reviewed a bunch of previous work on this very matter. I will mention couple of references from that paper which should clarify the problem and solve it.
What is wrong is that the generator of gauge transformation in a certain constrained system is not an isolated first class constraint, but a certain linear combination of first class constraints. This is what the mentioned Pitts' paper tries to emphasize. Somehow, people worked by inertia ever since Dirac's lecture notes, and never got to notice that a single first class constraint does not generate a gauge transformation, but a ''bad change'' (as Pitts calls it) as you have stated in your question. Then came Castellani in 1982, and in his paper "Symmetries in constrained systems" in Annals Phys. 143, p. 357 (1982) formulated a generator of gauge symmetries as a well-defined linear combination of first class constraints. This paper is very insightful and I recommend it as a starting point when one starts with constrained systems. There he derives an algorithm which determines the form of a gauge generator, and then he derives gauge generators for a simple toy model, for General Relativity, and for Yang-Mills theories (of which Electromagnetism is a special case, so the results also apply to the question above). All of them are a linear combination of first class constraints.
There is also a nice discussion and possibly a very detailed answer to the exact question posted above, in paper by Pons, as well as in Sundermeyer's book "Symmetries in Fundamental Physics".
The bottom line is: isolated first class constraints (primary or secondary or tertiary...) do not, in general, generate gauge transformations. But they are each a part of a gauge generator, which is defined as a linear combination of these constraints.
A: The fundamental problem here is that many people, and also Pitts in his paper, are not careful about what theory they are currently talking about. "Quantization of Gauge Systems" by Henneaux and Teitelboim is actually careful about this, and their chapter 3 shows the correct resolution of this problem, even though Pitts cites it as an example for those who do not recognize the problem.
The claim "first-class constraints generate gauge transformations" is in the context of the extended action
$$ S_E[q^i, p_i, \lambda^i] = \int \left(p_i \dot{q}^i - H - \lambda^i \gamma_i\right)\mathrm{d}t,$$
where the $\gamma_i$ are the first-class constraints, $\lambda^i$ the Lagrange multipliers enforcing the constraints and for simplicity we assume there are no second-class constraints (as is the case in the example of electromagnetism). In this formulation, solutions to the equations of motion are tuples $(q(t), p(t), \lambda(t))$, and symmetries act on all these dynamical variables. The extended action is invariant under the infinitesimal local transformation $\delta_{\epsilon} F = \epsilon^i \{F,\gamma_i\}$ for $F(q,p)$ any function of the $q$ and $p$ only if one additionally lets the Lagrange multipliers transform as
$$ \delta_{\epsilon}\lambda^i = \dot{\epsilon}^i + {C^i}_{jk}\lambda^j\epsilon^k,$$
where $\{\gamma_i,\gamma_j\} = {C^k}_{ij}\gamma_k$. This set of symmetries reduces to the symmetries of the non-extended (canonical) action
$$ S_C[q^i,p_i,\bar{\lambda}^i] = \int \left(p_i \dot{q}^i - H - \lambda^{i'} \gamma_{i'}\right),$$
where the $i'$ index now only runs over the primary constraints, only after imposing the gauge conditions $\lambda^j = 0$ for all $j$ where $\gamma_j$ is not primary. It is the symmetries of the canonical action, not of the extended action, that directly translate to symmetries of the original Lagrangian action. The residual gauge symmetries of this action are those that preserve the conditions $\lambda^j = 0$ for the non-primaries, and are in general generated by a specific subset of combinations of the first-class constraints that other sources call the gauge generator(s).
For the concrete example Pitts is complaining about, the gauge transformations of the extended action of free electrodynamics are (see also H/T's chapter 19):
\begin{align} 
\delta A_0 & = \epsilon^1 & \delta A_i & = \partial_i\epsilon^2 \\
\delta \lambda^1 & = \dot{\epsilon}^1&  \delta\lambda^2 & = \dot{\epsilon}^2 - \epsilon^1
\end{align}
for two arbitrary functions of spacetime $\epsilon^1, \epsilon^2$. The preservation of the gauge condition $\lambda^2 = 0$ imposes $\epsilon^1 = \dot{\epsilon^2}$, so we are left with the residual gauge symmetry
$$ \delta A_\mu = \partial_\mu \epsilon^2 \quad \delta \lambda^1 = \ddot{\epsilon}^2,$$
which now has the familiar form for a gauge transformation of the 4-potential $A_\mu$ of Lagrangian electrodynamics.
The thing that seems to trip people up is that in the extended formalism, the quantity $E^i = F^{0i} = \partial^0 A^i - \partial^i A^0$, which they would like to identify as the electric field, is gauge-variant under transformations with $\epsilon^1 \neq \dot{\epsilon}^2$. How can the extended theory be "faithful" to our real, physical system if it turns gauge-invariant quantities into gauge-variant ones?
To see how, let us inspect the change of $E^i$ under an arbitrary transformation:
\begin{align} 
\delta E_i & = \partial_0 \delta A_i - \partial_i \delta A_0 \\
& = \partial_i (\dot{\epsilon}^2 - \epsilon^1)
\end{align}
We note that this is precisely the spatial derivative of the transformation behaviour of $\lambda^2$, so this means that $E^i - \partial^i\lambda^2$ is a gauge-invariant observable, that turns into the residual gauge-invariant $E^i$ upon use of the gauge condition $\lambda^2 = 0$. So the extended theory does contain a gauge-invariant observable that is the electric field, it is just that its expression contains the auxiliary variable $\lambda^2$ which we have to eliminate to see whether or not this theory is equivalent to the Lagrangian theory that doesn't know about the $\lambda^i$.
