Counting Degrees of Freedom in Field Theories

I have very limited knowledge on this, but I can try to offer a partial answer.

The 4-potential $A^\mu$ has four degrees of freedom (d.o.f.) but two of these are unphysical and can be eliminated exploiting electromagnetism's invariance under gauge transformations $A^\mu \rightarrow {A^\prime}^\mu = A^\mu + \partial^{\,\mu} f$. For instance, we can take $\partial_\mu {A^\prime}^\mu \equiv 0$ as long as $f$ satisfies $\square\,f=-\partial_\mu A^\mu$:

${A^\prime}^\mu = A^\mu + \partial^{\,\mu} f \Longrightarrow \partial_\mu{A^\prime}^\mu \equiv 0 = \partial_\mu A^\mu + \square\,f \Longrightarrow \square\,f = -\partial_\mu A^\mu$

Condition $\partial_\mu {A^\prime}^\mu \equiv 0$, known as Lorenz gauge, gives you one equation relating $A_0, A_1, A_2$ and $A_3$, therefore it removes one out of four d.o.f.

The last unphysical d.o.f. is removed because there is still some gauge freedom to be explored. To ensure $\partial_\mu {A^\prime}^\mu \equiv 0$ all we have to do is use a function $f$ satisfying $\square\,f = -\partial_\mu A^\mu$. It is clear any function $f_0$ satisfying $\square\,f_0\equiv 0$ can be added to $f$ while keeping $\partial_\mu {A^\prime}^\mu \equiv 0$ intact. This last freedom removes one more d.o.f. as we can take ${A^{\prime\prime}}^\mu = {A^\prime}^\mu + \partial^{\,\mu} f_0$ and choose an appropriate solution of $\square\,f_0\equiv 0$.

The last paragraph may be best understoond with an example. In free space, using the Lorenz gauge $\partial_\mu {A^\prime}^\mu \equiv 0$, the 4-potential satisfies $\square\, {A^\prime}^\mu = 0$. Consider the solution for a plane wave propagating along the $z$ axis: ${A^\prime}^\mu = \mathcal{A}\, \varepsilon^\mu \cos(kz-\omega t)$, where $\mathcal{A}$ is the wave's amplitude, $\omega=|k|$ (I'm using natural units where $c\equiv 1$), and $\varepsilon^\mu$ is the polarization vector. You can check gauge condition $\partial_\mu {A^\prime}^\mu \equiv 0$ is satisfied for the polarization vector $\varepsilon^\mu$ chosen as $\varepsilon^\mu_{(x)} = (0,1,0,0)$ or $\varepsilon^\mu_{(y)} = (0,0,1,0)$, or some combination of these. These correspond to the two physical d.o.f. associated with horizontal and vertical polarizations.

For this example, the two unphysical d.o.f. are related to polarization vectors $\varepsilon^\mu_{(t)} = (1,0,0,0)$ and $\varepsilon^\mu_{(z)} = (0,0,0,1)$, corresponding to "temporal" and longitudinal polarizations, respectively. Although ${A^\prime}^\mu$ with $\varepsilon^\mu = \varepsilon^\mu_{(t)}$ or $\varepsilon^\mu = \varepsilon^\mu_{(z)}$ satisfies the equation of motion $\square\,{A^\prime}^\mu=0$, it won't satisfy the Lorenz gauge $\partial_\mu {A^\prime}^\mu \equiv 0$. Nevertheless, if the polarization vector is taken as the combination $\varepsilon^\mu = c_1 \varepsilon^\mu_{(t)} + c_2 \varepsilon^\mu_{(z)}$, for constants $c_1$ and $c_2$, Lorenz gauge is satisfied once we impose $c_1 = c_2$. Because of that, we are left with only one unphysical d.o.f. Finally, choosing $f_0$ (from previous paragraphs) as $\mathcal{A}\,\frac{c_1}{\omega}\sin(kz-\omega t)$, you can check the solution ${A^{\prime\prime}}^\mu = \mathcal{A}\,c_1(\varepsilon^\mu_{(t)}+\varepsilon^\mu_{(z)})\cos(kz-\omega t) + \partial^{\,\mu} f_o$ identically vanishes, meaning the last unphysical d.o.f. doesn't propagate indeed.

This is what is meant by the degrees of freedom of the electromagnetic field: plane wave solutions have two possible polarizations (or combinations of such), both are spatial and transverse to the propagation direction. In the second quantized theory, this means photons' spin one ($s=1$) has only two possible orientations relative to its motion: parallel ($m_s=+1$, positive helicity) or antiparallel ($m_s=-1$, negative helicity). Notice this has a very different meaning compared to the statement "field theories have infinite d.o.f.", which is related to describing an infinite number of points in spacetime (as opposed to a finite number in particle mechanics).

You really should split your question. I will answer the part where you do not understand how counting of degrees of freedom work.

Basically we count the number of propagating (physical) degrees of freedom per point of spacetime. Of course, the total number of degrees of freedom is infinite because spacetime is continuous and has an infinite number of points, but to ask for the number of degrees of freedom per spacetime point is a reasonable demand to make. Bear in mind that we only care about physical degrees of freedom by which we mean those that can be properly normalized.

You correctly state that photons can be off-shell but they are only those involved in internal processes. External photons are always on-shell. Moreover, gauge invariance is a physical property. External fields which you measure in your laboratory should be independent of your chosen gauge. In other words, the S-matrix should be gauge-invariant. On the other hand, there is nothing that stops me from having gauge-broken internal processes if ultimately I can make the S-matrix gauge-invariant. Therefore, the word "physical" should almost always give you a picture of external on-shell gauge-invariant quantities.

So yes, gauge redundancy kills one degree of freedom, and when we are talking about propagating physical degrees of freedom, one more is killed on-shell. You have to understand how that happens. It is not that every time you see an equation of motion, a degree of freedom is killed. Killing of degrees of freedom requires an elaborate process of imposing constraints on the equation of motion known as gauge-fixing. And this has to be done on a case by case basis.

For example, consider the four equations of motion (separated into temporal and spatial sets) for the massless photon $A^\mu = (\phi, \vec A)$ describing four on-shell degrees of freedom as follows.

\begin{align*} -\Delta \phi + \partial_t \vec\nabla\cdot\vec A = 0\,,\\ \square \vec A - \vec\nabla(\partial_t\phi-\vec\nabla\cdot\vec A) = 0\,.\\ \end{align*}

Since these equations exhibit a gauge symmetry $A_\mu \to A'_\mu := A_\mu + \partial_\mu \alpha_1(x)$, we can try to fix the gauge by choosing $\alpha_1$ such that, for instance, it is a solution of $\square \alpha_1 = -\vec\nabla\cdot\vec A$, giving us

\begin{align*} \Delta \phi' = 0\,, \\ \square \vec A' - \vec\nabla\partial_t\phi' = 0\,. \\ \\ \vec\nabla\cdot\vec{A}'=0\,.\\ \end{align*}

We have selected a divergence-free field, the so-called Coulomb gauge. Under this choice, the electric potential becomes non-propagating, that is there are no kinetic terms in the Lagrangian for it (observe that $\Delta \phi' = 0$ does not have any time derivatives).

In momentum space, this gauge condition reads $\vec p \cdot \vec \epsilon = 0$ where $ \vec \epsilon$ is the polarisation vector (Fourier transform of the magnetic potential). There are three solutions to this constraint. Choosing a frame in which $p^\mu = (E,0,0,E)$, we find that the three polarisation vectors are

$$ \epsilon^\mu_1 = (0,1,0,0), \qquad \epsilon_2^\mu=(0,0,1,0), \qquad \epsilon_t^\mu = (1,0,0,0) $$

The third polarisation is time-like and therefore cannot be normalized. It is unphysical, and we have to get rid of it. Luckily, the gauge symmetry is not exhausted. There are more available choices of gauge transformations which preserve the Coulomb gauge $\vec p \cdot \vec \epsilon = 0$. For example, we could go from $A'_\mu \to A_\mu:= A'_\mu + \partial_\mu \alpha_2(x)$ such that $\Delta \alpha_2 = 0,\ \partial_t \alpha_2 = - \phi'$ which preserves the divergence and sets $\phi = 0$.

Note that this time we have to make sure that this gauge transformation happens on-shell, namely that $\Delta \phi = 0$, otherwise this gauge-fixing will be inconsistent because $\Delta \alpha_2 = 0 \Rightarrow$ $0 = \Delta \partial_t\alpha_2 = - \Delta\phi' \ne 0$ off-shell. In other words, requiring $\phi = 0$, or equivalently $\epsilon^0 = 0$, in order to get rid of unphysical degrees of freedom requires us to be on-shell.

To summarize, we made an off-shell gauge choice $\vec p \cdot \vec \epsilon = 0$, an on-shell gauge choice $\epsilon^0 = 0$ and our equation of motion became $p^2 = 0$. Having exhausted our gauge choices, we find only two physical polarization modes or degrees of freedom.

Now, you understand that merely having an equation of motion does not eat up a degree of freedom. To find the correct number of degrees of freedom, keep on making gauge choices (producing independent constraint equations), some off-shell and some on-shell, until you exhaust your gauge freedom. Then check how many degrees of freedom you are left with. If you notice any unphysical guy showing up, most likely you haven't used up all your gauge freedom and you still have enough flex to shoot this guy dead. Then, count all that you are left with. That's your answer.

The only counting of DOF one has is made in the Hamiltonian analysis, i.e. one starts off with a valid Lagrangian (density) and computes the Hamiltonian (density). If constraints are found (such as the Yang-Mills fields for an SU(N) gauge group), then the general formula simply returns the net number of fields (which can be taken to be a parametrization of the reduced phase space under the Dirac bracket).

Example: The U(1) gauge field [classical electromagnetism]. The number of Lagrangian fields (DOF) = $4 \times \infty$. The Hamiltonian analysis returns 2 constraints (one primary, one secondary, both 1st class), therefore the true number of DOF of the U(1) gauge field is $ 2\times\infty$.