Your expression is correct, but the result is not gauge invariant. Noether's theorem can only be applied to global symmetry, and cannot be applied to local invariance (aka gauge redundancy).
The physical symmetry group acts on a physical state in the Hilbert space by transforming it into other physical states. However, physical states should form a trivial representation of the gauge group.
In what follows, the following conventions will be used:
$$\left\{g_{\mu\nu}\right\}=\mathrm{diag}(1,-1,-1,-1)$$
$$\left\{\partial_{\mu}\right\}=(\partial_{t},\nabla)\quad\left\{A_{\mu}\right\}=(\phi,\mathbf{A})$$
$$\left\{F^{\mu\nu}\right\}=\begin{pmatrix}
0 & -E^{1} & -E^{2} & -E^{3} \\
E^{1} & 0 & -B^{3} & B^{2} \\
E^{2} & B^{3} & 0 & -B^{1} \\
E^{3} & -B^{2} & B^{1} & 0
\end{pmatrix}$$
- Now take the $U(1)$-gauge theory as an example:
\begin{align}
\mathcal{L}&=-\frac{1}{4}F_{\mu\nu}F^{\mu\nu}+\bar{\Psi}(iD\!\!\!\!/-m)\Psi \tag{1.0} \\
&=\frac{1}{2}(\|\mathbf{E}\|^{2}-\|\mathbf{B}\|^{2})+\bar{\Psi}(iD\!\!\!\!/-m)\Psi,
\end{align}
where $$F^{0i}=(\frac{\partial\mathbf{A}}{\partial t}+\nabla\phi)_{i}=-E^{i},\quad F^{ij}=-\sum_{k=1}^{3}\epsilon^{ijk}B^{k},$$
and $D_{\mu}=\partial_{\mu}-iA_{\mu}$ is the covariant derivative.
The equations of motion of the above spinor-QED are
$$\partial_{\mu}F^{\mu\nu}=j^{\nu} \tag{1.1}$$
$$(iD\!\!\!\!/-m)\Psi=0 \tag{1.2}$$
where the spinor current $$j^{\mu}=-\bar{\Psi}\gamma^{\mu}\Psi \tag{1.3}$$
is conserved on-shell, because $\partial_{\mu}j^{\mu}=\partial_{\mu}\partial_{\nu}F^{\mu\nu}=0$.
Global $U(1)$-Symmetry
In fact, the current $j^{\mu}$ is also the Noether current associated with the global $U(1)$ symmetry $$\Psi(x)\rightarrow e^{i\Lambda}\Psi(x)\quad\mathrm{or}\quad\delta\Psi(x)=i\Lambda\Psi(x)$$
$$\bar{\Psi}(x)\rightarrow\bar{\Psi}(x)e^{-i\Lambda}\quad\mathrm{or}\quad\delta\bar{\Psi}(x)=-i\bar{\Psi}(x)\Lambda$$
$$A^{\mu}(x)\rightarrow A^{\mu}(x)$$
where $\Lambda\in\mathbb{R}$ is an arbitrary constant. This can be easily verified. Using Noether's (1st) theorem, the following current $$J^{\mu}=\frac{\partial\mathcal{L}}{\partial(\partial_{\mu}\Psi)}\delta\Psi+\delta\bar{\Psi}\frac{\partial\mathcal{L}}{\partial(\partial_{\mu}\bar{\Psi})}=i\Lambda\frac{\partial\mathcal{L}}{\partial(\partial_{\mu}\Psi)}\Psi-i\Lambda\bar{\Psi}\frac{\partial\mathcal{L}}{\partial(\partial_{\mu}\bar{\Psi})}=-\Lambda\bar{\Psi}\gamma^{\mu}\Psi,$$
is conserved on-shell. Up to a constant factor, the above Noether current equals to the current (1.3) from the equations of motion. It's conservation law is a further constraint on the space of solutions to the equations of motion.
$U(1)$-Gauge Invariance
The action (1.0) as well as EOM (1.1) and (1.2) are also invariant under the following $U(1)$-gauge transformations:
$$\Psi(x)\rightarrow e^{i\Lambda(x)}\Psi(x)\quad\mathrm{or}\quad\delta\Psi(x)=i\Lambda(x)\Psi(x)$$
$$\bar{\Psi}(x)\rightarrow\bar{\Psi}(x)e^{-i\Lambda(x)}\quad\mathrm{or}\quad\delta\bar{\Psi}(x)=-i\bar{\Psi}(x)\Lambda(x)$$
$$A^{\mu}(x)\rightarrow A^{\mu}(x)+\partial^{\mu}\Lambda(x)\quad\mathrm{or}\quad \delta A^{\mu}(x)=\partial^{\mu}\Lambda(x)$$
As explained previously, the Noether's theorem only applies to global symmetries. However, there's no obstacles to plug the above expressions of variations into the expression of the Noether current. One obtains
\begin{align}
K^{\mu}(x)&=\frac{\partial\mathcal{L}}{\partial(\partial_{\mu}\Psi)}\delta\Psi+\delta\bar{\Psi}\frac{\partial\mathcal{L}}{\partial(\partial_{\mu}\bar{\Psi})}+\frac{\partial\mathcal{L}}{\partial(\partial_{\mu}A_{v})}\delta A_{v} \\
&=-\Lambda(x)\bar{\Psi}\gamma^{\mu}\Psi-F^{\mu\nu}\partial_{\nu}\Lambda(x) \\
&=\Lambda(x)j^{\mu}(x)-F^{\mu\nu}\partial_{\nu}\Lambda(x) \\
&=-\partial_{\nu}[F^{\mu\nu}\Lambda(x)],
\end{align}
where in the last line the equations of motion (1.1) has been used. Clearly, the current $K^{\mu}$ is conserved, but it also depends on the choice of gauge due to the appearance of the function $\Lambda(x)$. As is explained in the beginning, gauge redundancies are not physical symmetries. Physical observables should remain unchanged under different choices of gauge. Thus, the Noether current $K^{\mu}$ is not a physical observable. Another way to look at it is to convert the charge $Q=\int d^{3}xK^{0}$ into a surface integral at infinity. But since $\Lambda(x)$ in general does not always vanish at infinity, one obtains an infinite family of redundant charges related by gauge transformations, and so they serve no practical purposes.
On the other hand, the Noether current $J^{\mu}$ is invariant under the above gauge transformations, thus corresponds to conserved physical charges. It is also clear from the expression of $K^{\mu}$ that the global $U(1)$ symmetry corresponds to the "global part" of the $U(1)$-gauge group that leaves the gauge field $A^{\mu}$ invariant.
Hamiltonian Formalism
In the Hamiltonian formalism, there should be a $3+1$ foliation of spacetime, and all fields are evaluated at a time slice.
In the fermionic sector, since classical fermions are anti-commuting Grassmann number-valued fields, one encounters an ambiguity when defining the canonical momentum of $\Psi(x)$. Specifically, one has the following two definitions:
$$\Pi_{R}=\mathcal{L}\frac{\overset{\leftarrow}{\delta}}{\delta\dot{\Psi}}\quad\mathrm{or}\quad\Pi_{L}=\frac{\overset{\rightarrow}{\delta}}{\delta\dot{\Psi}}\mathcal{L},$$
which are related via $\Pi_{R}=-\Pi_{L}$ because of the anti-commuting relation $$\left\{\Psi,\frac{\overset{\rightarrow}{\delta}}{\delta\Psi^{\dagger}}\right\}=0.$$
Both the left-wing and right-wing definitions serve equal definitions of the Hamiltonian density because the Lagrangian density should remain the same in both conventions as $$\mathcal{L}=\Pi_{R}\dot{\Psi}-\mathcal{H}=\dot{\Psi}\Pi_{L}-\mathcal{H}.$$
In order to maintain consistent, in the following only the right-wing definition is used, and $\Pi_{R}\equiv\Pi$. In this convention, the Poisson bracket between any two functionals is defined as $$\left\{F,G\right\}_{PB}=\int d^{3}x\left\{F\left(\frac{\overset{\leftarrow}{\delta}}{\delta\Psi(\mathbf{x})}\frac{\overset{\rightarrow}{\delta}}{\delta\Pi(\mathbf{x})}+\frac{\overset{\leftarrow}{\delta}}{\delta\Pi(\mathbf{x})}\frac{\overset{\rightarrow}{\delta}}{\delta\Psi(\mathbf{x})}\right)G\right\},$$
and so the fundamental Poisson bracket in the fermionic sector reads $$\left\{\Psi(\mathbf{x}),\Pi(\mathbf{y})\right\}_{PB}=\delta(\mathbf{x}-\mathbf{y}).$$
With the above expression of the fermionic Poisson bracket, one expects that the time evolution of any functional $F$ on the phase space should be $$\frac{dF}{dt}=\frac{\partial F}{\partial t}+\left\{F,H\right\}_{PB},$$
where $H$ is the Hamiltonian. From the fermionic sector in Lagrangian (1.0), one can easily check that the canonical momentum is $$\Pi=\mathcal{L}\frac{\overset{\leftarrow}{\delta}}{\delta\dot{\Psi}}=i\bar{\Psi}\gamma^{0}.$$
So the Hamiltonian of the fermionic sector is
\begin{align}
H_{f}&=\int d^{3}x\left(\Pi\dot{\Psi}-\mathcal{L}_{f}\right) \\
&=\int d^{3}x\left(-\Pi\gamma^{0}(\pmb{\gamma}\cdot\nabla+im)\Psi-A_{\mu}\bar{\Psi}\gamma^{\mu}\Psi\right) \\
&=\int d^{3}x\left(-\Pi\gamma^{0}(\pmb{\gamma}\cdot\nabla+im)\Psi+A^{0}\rho+\mathbf{A}\cdot\mathbf{j}\right),
\end{align}
where $\mathcal{L}_{f}=\bar{\Psi}(iD\!\!\!\!/-m)\Psi$ is the Lagrangian density of the fermionic sector, and
$$\rho=-\bar{\Psi}\gamma^{0}\Psi=i\Pi\Psi$$
$$\mathbf{j}=-\bar{\Psi}\pmb{\gamma}\Psi=i\Pi\gamma^{0}\pmb{\gamma}\Psi.$$
In the bosonic sector, one observes that the time derivative of the gauge potential $\phi=A^{0}$ does not appear in the Lagrangian. Thus, as a second order formalism the Lagrangian of the electromagnetic field is singular. The canonical momentum of $A^{i}$ is $$\Pi^{i}=\frac{\partial\mathcal{L}}{\partial\dot{A}_{i}}=-\dot{A}^{i}+\partial^{i}A^{0}=E^{i}.$$
For any functionals defined on this "naive" phase space of photons, the Poisson bracket is defined as $$\left\{F,G\right\}_{PB}=\sum_{i=1}^{3}\int d^{3}x\left\{\frac{\delta F}{\delta A^{i}(\mathbf{x})}\frac{\delta G}{\delta\Pi^{i}(\mathbf{x})}-\frac{\delta F}{\delta\Pi^{i}(\mathbf{x})}\frac{\delta G}{\delta A^{i}(\mathbf{x})}\right\},$$
and so the fundamntal Poisson bracket in the bosonic sector reads $$\left\{A_{i}(\mathbf{x}),\Pi_{j}(\mathbf{y})\right\}_{PB}=\delta_{ij}\delta(\mathbf{x}-\mathbf{y}). \tag{1.4}$$
So the Hamiltonian of the pure electromagnetic field is
\begin{align}
H_{\mathrm{em}}&=\int d^{3}x\left(\Pi^{i}\dot{A}_{i}-\mathcal{L}_{\mathrm{em}}\right) \\
&=\int d^{3}x\left(\Pi_{i}(\partial^{0}A^{i}-\partial^{i}A^{0})+\Pi_{i}\partial^{i}A^{0}-\mathcal{L}_{\mathrm{em}}\right) \\
&=\int d^{3}x\left(-\Pi_{i}E^{i}+\partial_{i}(\Pi^{i}A^{0})-A^{0}\partial_{i}\Pi^{i}-\mathcal{L}_{\mathrm{em}}\right) \\
&=\int d^{3}x\left(\frac{1}{2}\|\mathbf{E}\|^{2}+\frac{1}{2}\|\mathbf{B}\|^{2}-A^{0}\nabla\cdot\mathbf{E}\right) \\
&=\int d^{3}x\left(\frac{1}{2}\|\mathbf{\Pi}\|^{2}+\frac{1}{2}\|\nabla\times\mathbf{A}\|^{2}-A^{0}\nabla\cdot\mathbf{\Pi}\right).
\end{align}
Thus, the total Hamiltonian is given by
$$H=H_{\mathrm{em}}+H_{f} \tag{1.5}$$
$$=\int d^{3}x\left(\frac{1}{2}\|\mathbf{\Pi}\|^{2}+\frac{1}{2}\|\nabla\times\mathbf{A}\|^{2}-A^{0}(\nabla\cdot\mathbf{\Pi}-\rho)-\Pi\gamma^{0}(\pmb{\gamma}\cdot\nabla+im)\Psi+\mathbf{A}\cdot\mathbf{j}\right).$$
From the above expression, one finds that the Gauss law $\nabla\cdot\mathbf{\Pi}-\rho=0$ appears in the total Hamiltonian density as a constraint, and $A^{0}$ plays the role of a Lagrange multipier, which is arbitrary. Apparently, the total Hamiltonian is gauge invariant. One can easily check that the Hamilton's equations
\begin{align}
\dot{\Psi}&=\left\{\Psi,H\right\}_{PB}=\int d^{3}x\left(\frac{\overset{\rightarrow}{\delta}}{\delta\Pi}H\right) \\
&=-\gamma^{0}(\pmb{\gamma}\cdot\nabla+im)\Psi+iA^{0}\Psi+i\gamma^{0}\mathbf{A}\cdot\pmb{\gamma}\Psi \\
\dot{\mathbf{A}}&=\left\{\mathbf{A},H\right\}_{PB}=-\mathbf{\Pi}-\nabla A^{0} \\
\dot{\mathbf{\Pi}}&=\left\{\mathbf{\Pi},H\right\}_{PB}=\nabla(\nabla\cdot\mathbf{A})-\nabla^{2}\mathbf{A}-\mathbf{j}=\nabla\times\mathbf{B}-\mathbf{j}
\end{align}
reproduce the Dirac equation $(iD\!\!\!\!/-m)\Psi=0$ and the Maxwell equations.
Dirac Constraints
In the above discussion, one obtains the full Hamiltonian (1.5) with two first class contraints $$\Pi^{0}=0\quad\mathrm{and}\quad\nabla\cdot\mathbf{\Pi}-\rho=0. \tag{1.6}$$
According to Dirac's algorithm, first class constraints are generators of gauge redundancies. Take the Gauss constraint as an example, one can define the following smeared generator $$\mathcal{G}[\lambda]=\int d^{3}x\left(\nabla\cdot\mathbf{\Pi}(\mathbf{x})-\rho(\mathbf{x})\right)\lambda(\mathbf{x}),$$
where $\lambda(\mathbf{x})$ is an infinitesimal gauge function. Then, it's easy to check
\begin{align}
&\delta_{\mathcal{G}}A_{i}(\mathbf{y})=\left\{A_{i}(\mathbf{y}),\mathcal{G}\right\}_{PB}=\partial_{i}\lambda(\mathbf{y}), \\
&\delta_{\mathcal{G}}\Pi_{i}(\mathbf{y})=\left\{\Pi_{i}(\mathbf{y}),\mathcal{G}\right\}_{PB}=0 \\
&\delta_{\mathcal{G}}\Psi(\mathbf{y})=\left\{\Psi(\mathbf{y}),\mathcal{G}\right\}_{PB}=i\lambda(\mathbf{y})\Psi(\mathbf{y}) \\
&\delta_{\mathcal{G}}\bar{\Psi}(\mathbf{y})=\left\{\bar{\Psi}(\mathbf{y}),\mathcal{G}\right\}_{PB}=-i\lambda(\mathbf{y})\bar{\Psi}(\mathbf{y})
\end{align}
from which one recognizes that it is indeed just an infinitesimal time-independent gauge transformation. Thus, on the "naive" phase space, the gauge transformation generated by $\mathcal{G}$ shifts $\mathbf{A}$ field by total derivatives, and leaves $\mathbf{\Pi}$ field unchanged, as is illustrated in the picture below.
Also notice that physical observables do not change under gauge transformations. i.e $$\left\{Q(t),\mathcal{G}[\lambda]\right\}_{PB}=0,$$
where $Q(t)\equiv\int d^{3}x\rho(\mathbf{x})\Lambda$ is the electric charge, and
\begin{align}
\left\{H,\mathcal{G}[\lambda]\right\}_{PB}&=\int d^{3}x\left[\left(\nabla^{2}\mathbf{A}-\nabla(\nabla\cdot\mathbf{A})+\mathbf{j}\right)\cdot\nabla\lambda-\mathbf{j}\cdot\nabla\lambda\right] \\
&=-\int d^{3}x\left(\nabla\times(\nabla\times\mathbf{A})\right)\cdot\nabla\lambda \\
&=-\int d^{3}x\left(\nabla\times\mathbf{B}\right)\cdot\nabla\lambda \\
&=-\int d^{3}x\left[\nabla\cdot(\mathbf{B}\times\nabla\lambda)+\mathbf{B}\cdot(\nabla\times\nabla\lambda)\right] \\
&=0
\end{align}
where in the first line the last term $-\mathbf{j}\cdot\nabla\lambda$ comes from the Poisson bracket in the fermionic sector, and in the last line one uses the fact that $\mathrm{curl}\circ\mathrm{grad}=0$. This implies that the shift in the phase space is a "do-nothing transformation" on physical observables. As is mentioned previously, the existence of such gauge redundancies means the "naive" phase space on which the Poisson bracket (1.4) is defined is too large. The goal is to find gauge fixing conditions to reduce all non-physical degrees of freedom so that all constraints become second class constraints. Since one has $\Pi^{0}=0$ and $A^{0}$ at the moment is arbitrary, it makes sense to impose $A^{0}=0$. Furthermore, considering the identity $$\nabla^{2}\frac{1}{4\pi\|\mathbf{x}-\mathbf{y}\|}=-\delta(\mathbf{x}-\mathbf{y}),$$
one finds that by performing a time-independent gauge transformation $$\mathbf{A}(\mathbf{x})\rightarrow\mathbf{A}(\mathbf{x})+\nabla_{x}\int d^{3}y\frac{1}{4\pi\|\mathbf{x}-\mathbf{y}\|}\nabla_{y}\cdot\mathbf{A}(\mathbf{y})$$
one can always set $\nabla\cdot\mathbf{A}=0$.
In summary, one obtains the following four constraints
\begin{align}
&\varphi_{1}\equiv A^{0}=0 \tag{1.7a} \\
&\varphi_{2}\equiv\Pi_{0}=0 \tag{1.7b} \\
&\varphi_{3}\equiv\nabla\cdot\mathbf{A}=0 \tag{1.7c} \\
&\varphi_{4}\equiv\nabla\cdot\mathbf{\Pi}-\rho=0 \tag{1.7d},
\end{align}
which is known as the coulomb gauge. They are now second class constraints because the matrix of their Poisson brackets is non-singular, i.e $$\mathbf{M}=\begin{pmatrix}
\left\{\varphi_{3}(\mathbf{x}),\varphi_{3}(\mathbf{y})\right\}_{PB} & \left\{\varphi_{3}(\mathbf{x}),\varphi_{4}(\mathbf{y})\right\}_{PB} \\
\left\{\varphi_{4}(\mathbf{x}),\varphi_{3}(\mathbf{y})\right\}_{PB} & \left\{\varphi_{4}(\mathbf{x}),\varphi_{4}(\mathbf{y})\right\}_{PB}
\end{pmatrix}=\begin{pmatrix}
0 & \nabla^{2}_{x} \\
-\nabla^{2}_{x} & 0
\end{pmatrix}\delta(\mathbf{x}-\mathbf{y}).$$
One should notice that the above second class constraints are incompatible with the Poisson bracket (1.4) because $$\left\{\nabla\cdot\mathbf{A}(\mathbf{x}),\Pi_{i}(\mathbf{y})\right\}_{PB}=\partial_{i}\delta(\mathbf{x}-\mathbf{y})\neq 0.$$
The remedy is to introduce the Dirac bracket defined as $$\left\{F,G\right\}_{DB}\equiv\left\{F,G\right\}_{PB}-\sum_{ij}\int d^{3}z\int d^{3}w\left\{F,\varphi_{i}(\mathbf{z})\right\}_{PB}(\mathbf{M}^{-1})_{ij}(\mathbf{z},\mathbf{w})\left\{\varphi_{j}(\mathbf{w}),G\right\}_{PB}.$$
Then, the fundamental equal-time canonical relation of the electromagnetic field should be replaced by
\begin{align}
\left\{A_{i}(\mathbf{x}),\Pi_{j}(\mathbf{y})\right\}_{DB}&=\delta_{ij}\delta(\mathbf{x}-\mathbf{y})-\frac{\partial^{2}}{\partial x^{i}\partial x^{j}}\frac{1}{4\pi\|\mathbf{x}-\mathbf{y}\|} \\
&=\left(\delta_{ij}-\frac{\partial_{i}\partial_{j}}{\nabla^{2}}\right)\delta(\mathbf{x}-\mathbf{y}). \tag{1.8}
\end{align}
Now the Hamilton's equations on the reduced phase space are
\begin{align}
&\dot{\mathbf{A}}=\left\{\mathbf{A},H\right\}_{DB}=-\mathbf{\Pi} \\
&\dot{\mathbf{\Pi}}=\left\{\mathbf{\Pi},H\right\}_{DB}=\nabla\times\mathbf{B}-\mathbf{j}.
\end{align}
One can easily show that on this reduced phase space, the redundant degrees of freedom are completely removed. i.e
\begin{align}
&\left\{\Pi_{i}(\mathbf{x}),\mathcal{G}[\lambda]\right\}_{DB}=0 \\
&\left\{A_{i}(\mathbf{x}),\mathcal{G}[\lambda]\right\}_{DB} \\
&=-\int d^{3}y\left[\lambda(\mathbf{y})\partial^{j}\left(\delta_{ij}-\frac{\partial_{i}\partial_{j}}{\Delta}\right)\delta(\mathbf{x}-\mathbf{y})\right] \\
&=-\int d^{3}y\lambda(\mathbf{y})\left[\partial_{i}\delta(\mathbf{x}-\mathbf{y})-\partial_{i}\delta(\mathbf{x}-\mathbf{y})\right] \\
&=0, \\
&\left\{\Psi(\mathbb{x}),\mathcal{G}[\lambda]\right\}_{DB}=0.
\end{align}
To sum up, on the reduced phase space where the Dirac bracket (1.8) is defined, one has two conserved charges: $$Q=\int d^{3}x\rho(\mathbf{x})\Lambda\quad\mathrm{and}\quad\mathcal{G}[\lambda]=\int d^{3}x\left(\nabla\cdot\mathbf{\Pi}(\mathbf{x})-\rho(\mathbf{x})\right)\lambda(\mathbf{x}).$$
The charge $Q$ generates the global $U(1)$-symmetry acting on $\Psi$, i.e $$\delta_{Q}\Psi(\mathbf{x})=\left\{\Psi(\mathbf{x}),Q\right\}_{PB}=i\Lambda\Psi(\mathbf{x}),$$
and charge $\mathcal{G}[\lambda]$ generates a "do-nothing" transformation.
Canonical Quantization
In quantum mechanics, a pure state means a unit element in the Hilbert space, which can be identified as a ray in the Hilbert space. In other words, the set of pure states in a quantum system can be identified as the projective Hilbert space: $$\mathbb{P}\mathscr{H}\equiv(\mathscr{H}\backslash\left\{0\right\})/\sim,$$
where the equivalence class $\sim$ means for any two vectors $|a\rangle$ and $|b\rangle\in\mathscr{H}$, $|a\rangle\sim|b\rangle$ iff there exists a non-zero complex number $\lambda\in\mathbb{C}^{\ast}$ such that $|a\rangle=\lambda|b\rangle$. The motivation behind the above construction is that in experiments, one can only measure the probability distribution, and the transition probability between two states $|a\rangle$ and $|b\rangle$ is $$\mathcal{T}(|a\rangle,|b\rangle)\equiv\frac{|\langle a|b\rangle|^{2}}{||a\rangle|^{2}||b\rangle|^{2}},$$
from which one finds that scaling by a non-zero complex number does not change the result of transition probability.
For this reason, in quantum mechanics, the symmetry of the given quantum system is defined as a projective automorphism such that $$\mathcal{T}(T\cdot|a\rangle,T\cdot|b\rangle)=\mathcal{T}(|a\rangle,|b\rangle)$$
It was Eugene Wigner who proved the following theorem:
Theorem: Every projective automorphism $T$ on $\mathbb{P}\mathscr{H}$ arises from either a unitary or a anti-unitary operator $U$ on $\mathscr{H}$, and $U$ is determined by $T$ up to a complex phase.
The consequence of this theorem is that if a system has a classical symmetry group, say $\mathcal{O}$, then for each $L\in\mathcal{O}$, there is a unitary (or anti-unitary) operator $\rho(L):\mathscr{H}\rightarrow\mathscr{H}$, such that two consecutive transformations should satisfy $$\rho(L_{2})\rho(L_{1})=\omega(L_{2},L_{2})\rho(L_{2}L_{1}),$$
where $L_{1}$ and $L_{2}\in\mathcal{O}$, and $\omega(L_{2},L_{1})$ is a complex phase. Such a function $\rho$ is known as a projective representation. In mathematical language, what the above theorem is showing is the following short exact sequence: $$1\xrightarrow{}U(1)\overset{\mathrm{diag}}{\hookrightarrow}U(\mathscr{H})\overset{\hat{P}}{\twoheadrightarrow}PU(\mathscr{H})\xrightarrow{}1$$
where $U(\mathscr{H})$ is the group of unitary operators on $\mathscr{H}$ and $PU(\mathscr{H})$ is the projectie unitary group of $\mathscr{H}$.
The essence of the above discussion is that in quantum mechanics the symmetry group one usually refers to is in fact the central extension of the actual symmetry group by $U(1)$! The importance of this fact is that if a transformation acting on a quantum state by giving it a global complex phase, it is deemed as trivial.
According to superselection rules in quantum mechanics, the Hilbert space of theory (1.0) should be a direct sum of the photon sector and the spinor sector, i.e $$\mathscr{H}=\mathscr{H}_{\mathrm{em}}\oplus\mathscr{H}_{f}$$
In the spinor sector, as is shown previously, there's a global $U(1)$-symmetry in the classical theory. One should check whether this symmetry acts trivially on quantum states.
First of all, in the fermionic sector, one should replace the classical anti-commuting Poisson bracket by the quantum anti-commutator $$[\Psi(\mathbf{x}),\Pi(\mathbf{y}]_{+}=-i\hbar\delta(\mathbf{x}-\mathbf{y}),$$
and expand the Dirac field in terms of creating and annihilation operators $$\Psi(t,\mathbf{x})=\int\frac{d^{3}k}{(2\pi)^{3}}\frac{m}{k_{0}}\sum_{a=1}^{2}\left[\mathbf{b}_{a}(k)u^{a}(k)e^{-ikx}+\mathbf{d}^{\dagger}_{a}(k)v^{a}(k)e^{ikx}\right],$$
where $\mathbf{b}_{a}$ is an annihilation operator, and $\mathbf{d}^{\dagger}_{a}$ is a creation operator.
Global $U(1)$ Symmetry
Using the above expansion, one finds that the anti-commutation relation is equivalent to $$\left[\mathbf{b}_{a}(k),\mathbf{b}^{\dagger}_{a^{\prime}}(k^{\prime})\right]_{+}=\left[\mathbf{d}_{a}(k),\mathbf{d}^{\dagger}_{a^{\prime}}(k^{\prime})\right]_{+}=(2\pi)^{3}\frac{k_{0}}{m}\delta(\mathbf{k}-\mathbf{k}^{\prime})\delta_{aa^{\prime}}.$$
Thus, the global $U(1)$ charge $Q$ can be written as
\begin{align}
\mathbf{Q}(t)&=\int d^{3}x:\rho(t,\mathbf{x}):=\int\frac{d^{3}k}{(2\pi)^{3}}\frac{m}{k_{0}}\sum_{a=1}^{2}\left[\mathbf{b}^{\dagger}_{a}(k)\mathbf{b}_{a}(k)-\mathbf{d}^{\dagger}_{a}(k)\mathbf{d}_{a}(k)\right] \\
&=\mathbf{N}_{b}-\mathbf{N}_{d}.
\end{align}
By Wightman axioms, one assumes that there is an invariant vacuum state $|0\rangle$ so that $$\mathbf{Q}|0\rangle=|0\rangle.$$
Next, consider the one-particle state of an electron: $|\mathbf{k}\rangle$. Since $Q$ consists of two parts, one expects that this one-particle state is doubly degenerate. In other words, one expects two linearly independent states $$|\uparrow,a,\mathbf{k}\rangle=\mathbf{b}^{\dagger}_{a}(k)|0\rangle\quad\mathrm{and}\quad|\downarrow,a,\mathbf{k}\rangle=\mathbf{d}^{\dagger}_{a}(k)|0\rangle.$$
According to superselection rules, the above two states do not talk to each other. Indeed, one can check that $$\mathbf{Q}|\uparrow,a,\mathbf{k}\rangle=+|\uparrow,a,\mathbf{k}\rangle\quad\mathrm{and}\quad\mathbf{Q}|\downarrow,a,\mathbf{k}\rangle=-|\downarrow,a,\mathbf{k}\rangle.$$
Now consider the exponential map $$U(\Lambda)=e^{i\Lambda\mathbf{Q}}$$
of the charge operator $\mathbf{Q}$. One finds $$e^{-i\Lambda\mathbf{Q}}\Psi(t,\mathbf{x})e^{i\Lambda\mathbf{Q}}=e^{i\Lambda}\Psi(t,\mathbf{x}).$$
Thus, $U(\Lambda)$ represents a finite global $U(1)$-transformation generated by $\mathbf{Q}$. It's easy to check that $$U(\Lambda)|\uparrow,a,\mathbf{k}\rangle=e^{i\Lambda}|\uparrow,a,\mathbf{k}\rangle\quad\mathrm{and}\quad U(\Lambda)|\downarrow,a,\mathbf{k}\rangle=e^{-i\Lambda}|\downarrow,a,\mathbf{k}\rangle,$$
from which one observes that the action of this global $U(1)$ symmetry on quantum states cannot be factored out as a global complex phase. It acts on the one particle states of an electron in the following way: $$U(\lambda)|a,\mathbf{k}\rangle=\begin{pmatrix}
e^{i\Lambda} & 0 \\
0 & e^{-i\Lambda}
\end{pmatrix}|a,\mathbf{k}\rangle,$$
which is non-trivial.
$U(1)$-Gauge Invariance
Before fixing the gauge, one can verify that on the naive phase space where the naive canonical commutation relation $$\left[A_{i}(\mathbf{x}),\Pi_{j}(\mathbf{y})\right]=-i\hbar\delta_{ij}\delta(\mathbf{x}-\mathbf{y}) \tag{2.0}$$
is defined, the Gauss constraint generates local $U(1)$-gauge transformations on canonical operators.
One can define the following $U(1)$-valued operator $$U[\lambda]=\exp\left(i\mathcal{G}\right)=\exp\left(i\int d^{3}x\left(\nabla\cdot\mathbf{\Pi}(\mathbf{x})-\rho(\mathbf{x})\right)\lambda(\mathbf{x})\right).$$
Then, using the naive canonical commutator, it follows that
\begin{align}
[A_{i}(\mathbf{x}),\mathcal{G}[\lambda]]=\partial_{i}\lambda(\mathbf{x}),\quad &\mathrm{or}\quad U[\lambda]A_{i}(\mathbf{x})U[\lambda]^{\dagger}=A_{i}(\mathbf{x})+\partial_{i}\lambda(\mathbf{x}), \\
[\Pi_{i}(\mathbf{x}),\mathcal{G}[\lambda]]=0,\quad &\mathrm{or}\quad U[\lambda]\Pi_{i}(\mathbf{x})U[\lambda]^{\dagger}=0, \\
[\Psi(\mathbf{x}),\mathcal{G}[\lambda]]=i\lambda(\mathbf{x})\Psi(\mathbf{x}),\quad &\mathrm{or}\quad U[\lambda]\Psi(\mathbf{x})U[\lambda]^{\dagger}=e^{i\lambda(\mathbf{x})}\Psi(\mathbf{x}),
\end{align}
which are indeed finite gauge transformations. Also, one can check that using the naive canonical commutator (2.0), one has
\begin{align}
[\mathbf{Q}(t),\mathcal{G}[\lambda]]=&0,\quad\mathrm{or}\quad U[\lambda]\mathbf{Q}(t)U[\lambda]^{\dagger}=\mathbf{Q}(t), \\
[H,\mathcal{G}[\lambda]]=&0,\quad\mathrm{or}\quad U[\lambda]HU[\lambda]^{\dagger}=H,
\end{align}
which indicate that the $U(1)$-valued operator $U[\lambda]$ is a do-nothing action on quantum observables.
Now, consider an eigenstate $|\psi\rangle$ of the set of quantum observables $\left\{H,\mathbf{P},\mathbf{Q},\cdots\right\}$. For convenience, an arbitrary element of the set of quantum observables is denoted by $\mathbf{Q}_{m}$, whose eigenvalue is denoted as $q_{m}$. Then, it follows that $$U[\lambda]\mathbf{Q}_{m}U[\lambda]^{\dagger}|\psi\rangle=\mathbf{Q}_{m}|\psi\rangle.$$
Since the operator $U[\lambda]$ is $U(1)$-valued, one obtains $$\mathbf{Q}_{m}U[\lambda]^{-1}|\psi\rangle=U[\lambda]^{-1}\mathbf{Q}_{m}|\psi\rangle=q_{m}U[\lambda]^{-1}|\psi\rangle.$$
This implies that if $|\psi\rangle$ is an eigenstate, then $U[\lambda]|\psi\rangle$ is also an eigenstate of the same set of quantum observables. In other words, one has $$U[\lambda]|\psi\rangle=c|\psi\rangle,$$
where $c\in\mathbb{C}^{\ast}$ is a global constant complex factor. Or, one can say $|\psi\rangle\sim U[\lambda]|\psi\rangle$. i.e they belong to the same ray in the projective Hilbert space. Without losing generality, one can safely assume that $c=1$, and the quantum mechanical Gauss law becomes $$U[\lambda]|\psi\rangle=|\psi\rangle.$$
Then, one should immediately realize that the above equation is equivalent to the identity $$\left(\nabla\cdot\Pi(\mathbf{x})-\rho(\mathbf{x})\right)|\psi\rangle=0, \tag{2.1d}$$
where $|\psi\rangle$ is an arbitrary physical state. That is to say, the Gauss constraint (1.7d) is promoted as a week condition (2.1d).
To get rid of this gauge redundancy, one follows the same procedure as is done in the classical field theory, by imposing the other constraints:
\begin{align}
&\varphi_{1}\equiv A^{0}=0 \tag{2.1a} \\
&\varphi_{2}\equiv\Pi_{0}=0 \tag{2.1b} \\
&\varphi_{3}\equiv\nabla\cdot\mathbf{A}=0 \tag{2.1c}
\end{align}
The above three extra quantum mechanical constraints are imposed as strong operator constraints. Now consider a Fourier expansion of the gauge field $$\mathbf{A}(t,\mathbf{x})=\int\frac{d^{3}p}{(2\pi)^{2}2p_{0}}\mathbf{A}(\mathbf{p},t)e^{i\mathbf{p}\cdot\mathbf{x}}.$$
The constraint (2.1c) becomes the transversality condition: $$\mathbf{p}\cdot\mathbf{A}(\mathbf{p},t)=0.$$
This implies that a photon has two independent polarizations. One can introduce unit vectors $\pmb{\epsilon}_{1}(\mathbf{p})$ and $\pmb{\epsilon}_{2}(\mathbf{p})$, such that $$\pmb{\epsilon}_{1}\cdot\pmb{\epsilon}_{2}=\pmb{\epsilon}_{1}\cdot\mathbf{p}=\pmb{\epsilon}_{2}\cdot\mathbf{p}=0,\quad\mathrm{and}\quad(\pmb{\epsilon}_{1})^{2}=(\pmb{\epsilon}_{2})^{2}=1.$$
Then, the canonical quantization of the gauge field takes the form $$\mathbf{A}(x)=\int\frac{d^{3}p}{(2\pi)^{2}2p_{0}}\sum_{\alpha=1}^{2}\pmb{\epsilon}_{\alpha}\left[\mathbf{a}_{\alpha}(\mathbf{p})e^{-ipx}+\mathbf{a}^{\dagger}_{\alpha}(\mathbf{p})e^{ipx}\right],$$
where $p^{2}=0$, and $p_{0}=\|\mathbf{p}\|$, and the creation and annihilation operator satisfy commutation relations: $$[\mathbf{a}_{\alpha}(\mathbf{p}),\mathbf{a}^{\dagger}_{\beta}(\mathbf{k})]=2p_{0}(2\pi)^{3}\delta(\mathbf{p}-\mathbf{k})\delta_{\alpha\beta}.$$
Then, after imposing the quantum mechanical constraints (2.1), one fixes all gauge redundancies. This is the coulomb gauge in QED. It is straightforwad to verify that the passage towards canonical quantization in the coulomb gauge is achieved by replacing the Dirac bracket (1.8) by the commutation relation $$\left[A_{i}(\mathbf{x}),\Pi_{j}(\mathbf{y})\right]=-i\hbar\left(\delta_{ij}-\frac{\partial_{i}\partial_{j}}{\Delta}\right)\delta(\mathbf{x}-\mathbf{y}), \tag{2.2}$$
which is compatible with all the quantum mechanical constraints in (2.1). It is straightforward to check that in the reduced phase space where the canonical commutation relation (2.2) is defined, all gauge redundancies are removed. i.e
\begin{align}
[\Pi_{i}(\mathbf{x}),\mathcal{G}[\lambda]]=0,\quad &\mathrm{or}\quad U[\lambda]\Pi_{i}(\mathbf{x})U[\lambda]^{\dagger}=\Pi_{i}(\mathbf{x}), \\
[A_{i}(\mathbf{x}),\mathcal{G}[\lambda]]=0,\quad &\mathrm{or}\quad U[\lambda]A_{i}(\mathbf{x})U[\lambda]^{\dagger}=A_{i}(\mathbf{x}), \\
[\Psi(\mathbf{x}),\mathcal{G}[\lambda]]=0,\quad &\mathrm{or}\quad U[\lambda]\Psi(\mathbf{x})U[\lambda]^{\dagger}=\Psi(\mathbf{x}).
\end{align}
