I am trying to compute the full set of second class constraints that specify (1+1)-dimensional $U(1)$ scalar QED after complete gauge fixing. (Specifically, I would like to use the Coulomb $\partial_1 A_1=0$ combined with temporal gauge $A_0=0$). Since this is an Abelian gauge theory I suppose that this is possible. The way I would like the model to be described is using periodic boundary conditions on the fields. (With a universe of spatial circumference $L$.) However, the Dirac algorithm does not seem to terminate and the secondary constraints keep piling up:
The Lagrangian that I am using is $$L = \int_{-L/2}^{L/2} dx \, \mathcal L = \int_{-L/2}^{L/2} dx \,\left( -\frac 1 4 F_{\mu\nu}F^{\mu\nu} + |D_\mu\phi|^2 - V(\phi) \right),$$ with metric convention $\eta=\text{diag}(1,-1)$, with $F_{\mu\nu} = \partial_\mu A_\nu - \partial_\nu A_\mu$ and where $D_\mu = \partial_\mu-igA_\mu$. Hence, $|D_\mu\phi|^2 = {D_\mu}^* D^\mu = (\partial_\mu+igA_\mu)\phi^*(\partial^\mu-igA^\mu)\phi$). For completeness, the equations of motion (and the definition of the conserved current $J^\mu$, i.e. $\partial_\mu J^\mu = 0$) are: $$\partial_\mu F^{\mu\nu} = J^\nu \equiv ig((D^\nu\phi)^*\phi - \phi^* D^\nu \phi)$$ $$D_\mu D^\mu \phi = - \frac{\partial V}{\partial \phi} \quad \text{(and its conjugate.)}$$
From here on out I will be using lowered indices for the most part (and I will be leaving out the arguments).
The momenta including our primary constraint appear to be (Note, $\approx$ is Dirac's weak equality, stating that this must hold on the constraint surface) $$\Pi_0 = \frac{\delta \mathcal L}{\delta (\partial_0 A_0)} \approx 0,$$ $$E_1 = \frac{\delta \mathcal L}{\delta (\partial_0 A_1)} = F_{01} = \partial_0 A_1 - \partial_1 A_0,$$ $$\pi^* = \frac{\delta \mathcal L}{\delta (\partial_0 \phi)} = (D_0\phi)^*,$$ $$\pi = \frac{\delta \mathcal L}{\delta (\partial_0 \phi^*)} = D_0\phi.$$
This means our primary Hamiltonian is (after doing the Legendre transform and crossing some terms out) $$H_p = \int_{-L/2}^{L/2} dx \, \left( \frac 1 2 E_1^2 + |\pi|^2 + |D_1\phi|^2 + V(\phi) + A_0 [J_0-\partial_1 E_1] \right),$$ where $J_0 = ig(\pi^* \phi - \pi \phi^*)$. (I will be assuming equal time commutation relations between the fields and their momentas from this point onwards, i.e. $\{\phi(x), \pi(y)\} = \delta(x-y)$. Here $\{\cdot,\cdot\}$ is the Poisson bracket.)
Let us now perform the (naive) Dirac algorithm (I won't treat $A_0$ as a Lagrange multiplier, Gauss's law will then be our secondary constraint). We begin with our primary constraint $\Gamma_1 \equiv \Pi_0 \approx 0$. Preservation of this constraint in time gives us our secondary constraint: $\Gamma_2 \equiv \dot \Gamma_1 = \{\Pi_0, H_p\} = J_0-\partial_1 E_1 \approx 0$, which is Gauss's law for this model. There are not further constraints after this point since $\dot \Gamma_2 = \{\Gamma_2, H_p\} = 0$ (strongly zero, since every term crosses out). This primary and secondary constraint together form a set of two first class constraints ($\{\Gamma_1, \Gamma_2\} = 0$).
Now we add by hand our first canonical gauge fixing condition as another constraint, the goal of adding these further constraints is ultimately for the set of them all to become second class (that each constraint has a non-vanishing Dirac bracket with at least one other constraint). This first condition will be to the Coulomb gauge condition $\Gamma_3 \equiv \partial_1 A_1 = 0$. Preservation in time of this constraint gives us $\Gamma_4 \equiv \dot \Gamma_3 = \{\Gamma_3, H_p\} = \partial_1 E_1 - \partial_1^2 A_0 \approx 0$. Usually (such as in free Maxwell in (3+1) dimensions), this leads to the conclusion that $A_0 = 0$.
This appear to be where it goes wrong: In this (1+1)-dimensional model this constraint (as well as Gauss's law) can be solved, it does not appear that $A_0 = 0$ is at all equivalent (So how can we then use the temporal gauge condition $A_0=0$?). More importantly, this is normally the last constraint at which the algorithm terminates, leaving us with four second class constraints. If we, however, simply take this constraint at face value we have to continue the algorithm and compute its preservation in time, which does not seem to vanish: $\Gamma_5 \equiv \dot \Gamma_4 = \{\Gamma_4, H_p\} = ig \partial_1 J_1 \approx 0$. In other words, this constraint is not a combination of the prior ones, so we need to add it. We can keep going like this for a long time adding more and more constraints, for instance $\Gamma_6 \sim O(g^3)$. Something appears to be going wrong, but what is it?
I have a vague feeling this might have something to do with the fact that the zero component the conserved current $J_0$ contains a covariant derivative $D_0$ and hence the field $A_0$, which would not have happened if the scalar would be replaced by, say, a fermion.