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.

You can write the (scaled) interaction part of the action as:
$$S_I \equiv \int_{\mathbb R^d}d^d \mathbf r \ \phi(\mathbf r)\int_{\mathbb R^d}d^d \mathbf r' \ K(\mathbf r-\mathbf r') \ \phi(\mathbf r')$$
Let's take the inner integral over $\mathbf r'$ first (I will call it $\mathcal I$ to make things easier). Expanding $\phi(\mathbf r')$ around $\mathbf r$ gives :
$$\mathcal I \equiv\int_{\mathbb R^d}d^d \mathbf r' \ K(\mathbf r-\mathbf r') \ \phi(\mathbf r') \approx \int_{\mathbb R^d}d^d \mathbf r' \ K(\mathbf r-\mathbf r') \ \bigg(\phi(\mathbf r)+ \sum_{i=1}^d (x_i'-x_i)\partial_i \phi(\mathbf r) \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ +\frac 12\sum_{i=1}^d\sum_{j=1}^d (x_i'-x_i)(x_j'-x_j)\partial_i \partial_j \phi(\mathbf r) \bigg)$$
Now take the integral inside to get:
$$\mathcal I\approx \phi(\mathbf r) \int_{\mathbb R^d}d^d \mathbf r' \ K(\mathbf r-\mathbf r') \ + \sum_{i=1}^d \partial_i \phi(\mathbf r)\int_{\mathbb R^d}d^d \mathbf r' \ (x_i'-x_i) K(\mathbf r-\mathbf r') +\frac 12\sum_{i=1}^d\sum_{j=1}^d \partial_i \partial_j \phi(\mathbf r) \times\int_{\mathbb R^d}d^d \mathbf r'(x_i'-x_i)(x_j'-x_j)K(\mathbf r-\mathbf r') \bigg)$$
Now assuming that the coupling is homogenous, $K(\mathbf r-\mathbf r')\equiv K(\mathbf r'-\mathbf r)$. With that in mind, and also changing variables $\mathbf R \equiv \mathbf r'-\mathbf r$, we get:

$$\mathcal I \approx \phi(\mathbf r) \int_{\mathbb R^d}d^d \mathbf R \ K(\mathbf R) \ + \sum_{i=1}^d \partial_i \phi(\mathbf r)\int_{\mathbb R^d}d^d \mathbf R \ R_i K(\mathbf R) +\frac 12\sum_{i=1}^d\sum_{j=1}^d \partial_i \partial_j \phi(\mathbf r) \ \ \ \ \times\int_{\mathbb R^d}d^d \mathbf R \ R_iR_jK(\mathbf R) \bigg)$$
You can relate each of the integrals over $\mathbf R$ to the Fourier transform of $K(\mathbf R)$ defined as $\tilde K(\mathbf q) \equiv \int_{\mathbb R^d} d^d \mathbf R \ K(\mathbf R)\exp(-i \mathbf{q} . \mathbf R)$:

- First integral:
$$\int_{\mathbb R^d}d^d \mathbf R \ K(\mathbf R) = \int_{\mathbb R^d}d^d \mathbf R \ K(\mathbf R) \ e^{-i \mathbf q. \mathbf R} |_{\mathbf q =0} = \tilde K(\mathbf 0)$$
- Second integral:
$$\int_{\mathbb R^d}d^d \mathbf R \ R_i K(\mathbf R) =0$$
Because of the integrand being odd as you mentioned.

- Third integral:

For this one we first note that as you mentioned the integral is zero for all different $i,j$. For $i=j$, first note that:
$$\frac{\partial^2}{\partial q_i^2} \int_{\mathbb R^d}d^d \mathbf R \ K(\mathbf R) \ e^{-i \mathbf q. \mathbf R} = \int_{\mathbb R^d}d^d \mathbf R \ (-i)(-i) R_i R_i K(\mathbf R) \ e^{-i \mathbf q. \mathbf R} \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \\ = - \int_{\mathbb R^d}d^d \mathbf R \ R_i^2 K(\mathbf R) \ e^{-i \mathbf q. \mathbf R} $$
Which implies:
$$\int_{\mathbb R^d}d^d \mathbf R \ R_i^2 K(\mathbf R)=-\frac{\partial^2}{\partial q_i^2} \int_{\mathbb R^d}d^d \mathbf R \ K(\mathbf R) \ e^{-i \mathbf q. \mathbf R}|_{\mathbf q=0}=-\frac{\partial^2}{\partial q_i^2}\tilde K(\mathbf q) |_{\mathbf q=0} $$
Now if you assume that the coupling is also isotropic, i.e. $\exists \mathcal K : K(\mathbf R) \equiv \mathcal K(|\mathbf R|)$, the Fourier transform of $K$ will become a single variable function, meaning that the third integral is just $-\tilde K''(0)$.

In summary, $\mathcal I$ is:
$$\mathcal I \approx \phi(\mathbf r) \tilde K(0) \ - \frac 12\sum_{i=1}^d \partial_i^2 \phi(\mathbf r) \tilde K''(0)$$
Thus, the interaction term in the action is:
$$S_I = \int_{\mathbb R^d}d^d \mathbf r \ \phi(\mathbf r) \mathcal I = \int_{\mathbb R^d}d^d \mathbf r \ \phi(\mathbf r)\bigg(\phi(\mathbf r) \tilde K(0) \ - \frac 12\sum_{i=1}^d\partial_i^2 \phi(\mathbf r) \tilde K''(0)\bigg)$$

$$=\tilde K(0)\int_{\mathbb R^d}d^d \mathbf r \ \phi^2(\mathbf r) - \frac {\tilde K''(0)}2\sum_{i=1}^d \int_{\mathbb R^d}d^d \mathbf r \ \phi(\mathbf r) \ \partial_i^2 \phi(\mathbf r) $$
Integrating by parts in the second term results in (boundary terms vanish because $\phi(\mathbf r) \to 0$ as $|\mathbf r| \to \infty$ so that the integrals converge):
$$S_I =\tilde K(0)\int_{\mathbb R^d}d^d \mathbf r \ \phi^2(\mathbf r) + \frac {\tilde K''(0)}2\sum_{i=1}^d\int_{\mathbb R^d}d^d \mathbf r \ \partial_i \phi(\mathbf r) \ \partial_i \phi(\mathbf r)$$
$$=\tilde K(0)\int_{\mathbb R^d}d^d \mathbf r \ \phi^2(\mathbf r) + \frac {\tilde K''(0)}2 \int_{\mathbb R^d}d^d \mathbf r \ \sum_{i=1}^d (\partial_i \phi(\mathbf r))^2$$

$$=\int_{\mathbb R^d}d^d \mathbf r \ \bigg( \tilde K(0) \phi^2(\mathbf r) + \frac {\tilde K''(0)}2 \ \big( \partial \phi(\mathbf r) \big)^2 \bigg)$$
Plugging this in the full action finally gives:
$$S[\phi] =\int_{\mathbb R^d}d^d \mathbf r \ \bigg( \frac {\tilde K''(0)}8 \ \big( \partial \phi(\mathbf r) \big)^2 + \left(\frac {\tilde K(0)}{4}- \frac 12\right) \phi^2(\mathbf r) + \frac 1{12} \phi^4(\mathbf r) \bigg) $$
Notice that the coefficient of the quadratic term can change sign with temperature (through $\tilde K$), which is a sign of a *phase transition*.

## Best Answer

In physics, infinite number of degrees of freedom means the state or configuration of a system cannot be given completely by finite number of variables, but requires infinite number of variables. These do not need to correspond to any physical space axes.

Often the infinite number of variables is due to working with field. Field is function of position, and since there is infinite number of positions, the field has infinite number of degrees of freedom.