Connection between the vanishing of the conjugate momentum $\pi^0$ and non-existence of propagator for the free EM field If we try to quantize the free electromagnetic field without a gauge fixing term added to the Lagrangian density $\mathscr{L}=-\frac{1}{4}F_{\mu\nu}F^{\mu\nu}$, (i)zeroth component of the conjugate momentum density $\pi^0=$ vanishes, and (ii)also the propagator doesn't exist where the reason for its non-existence is usually attributed to the operator $(g^{\lambda\mu}\square-\partial^\mu\partial^\lambda)$ being non-invertible (singular). 
My question is whether (i)the vanishing of $\pi^0$ and (ii) the non-existence of the propagator is related? 
I think they are related because the problem of non-existence of the propagator doesn't arise in the case of a scalar field or Dirac field where $\pi^0\neq 0$. Moreover, the fixation of gauge solves both the problems at one shot. But I'm not sure where is the relation connection these two problems.
 A: Gauge theories become constrained Hamiltonian systems when passing from the Lagrangian $L(q,\dot{q},t)$ to the Hamiltonian $H(q,p,t)$ where $p = \frac{\partial L}{\partial \dot{q}}$. Generically, you get a constrained Hamiltonian system whenever the matrix/operator with components
$$ \frac{\partial^2 L}{\partial \dot{q}^i\partial\dot{q}^j}$$
is singular, i.e. non-invertible, i.e $\det(\frac{\partial^2 L}{\partial \dot{q}\partial\dot{q}}) = 0$. As you correctly observe, this is already the case for a massive vector field.
So, let us look at the Lagrangian of a generic vector field:
$$ L(A,\dot{A}) = \int (\frac{1}{4}F^{\mu\nu}F_{\mu\nu} + \mu^2 A^\mu A_\mu) \mathrm{d}^3x$$
Regardless of $\mu$, the canonical momenta are
$$ \pi^\mu = F^{\mu 0}$$
so we have the primary constraint1
$$\pi^0 \approx 0\tag{1}$$
always. The canonical Hamiltonian reads
$$ H = \int (\frac{1}{2}\pi^i\pi_i - \frac{1}{4} F^{ij}F_{ij} - A_0 \partial_i \pi^i -\mu^2 A^\mu A_\mu)\mathrm{d}^3 x$$
and, for consistency of the constraint, we incur a secondary constraint
$$ \dot{\pi}^0 = \{\pi^0,H\} = \partial_i\pi^i + \mu^2 A_0 \approx 0 \tag{2}$$
using that the only non-zero Poisson bracket of $\pi_0$ is $\{A_0,\pi^0\} = 1$. The nature of the theory is now very different depending on whether or not $\mu^2 = 0$.
If $\mu^2\neq 0$, then $(2)$ does not, effectively, impose a restriction on $\partial_i \pi^i$. Solving the two constraints means just putting $A_0 = -\frac{1}{\mu^2}\partial_i\pi^i$ and $\pi^0 = 0$, meaning we reduce the phase space dimension by two (effectively forgetting there was a pair of coordinates $(A_0, \pi^0)$) and are in a constraint-free theory. On this reduced phase space, canonical quantization may now proceed as usual, and we do not have gauge degrees of freedom left. In particular, canonical quantization delivers a propagator, since the operator is invertible on the degrees of freedom that are left.
If $\mu^2 = 0$, then $(2)$ is just Gauss' law $\nabla \cdot \vec E = 0$ since $\pi^i = F^{i0} = E^i$. Although it is in principle possible to locally solve this constraint and again pass to a reduced phase space (with lower dimension than before), this is not actually feasible or desirable in practice. Hence, canonical quantization as for unconstrained systems is not possible, and we do not obtain a propagator if we try to calculate one, since there are gauge degrees of freedom left.

1$\approx$ denotes weak equalities which hold on the constraint surface, but are not identically zero throughout the entire phase space.
A: This is indeed the same problem, from different viewpoints. The deep source of the problem is that the number of $A_\mu$ fields is higher than the actual physical degrees of freedom. Here are a few remarks:


*

*The operator you are trying to invert does not have an inverse because it
is just a projection operator(you can check this by direct computation, the
square of the operator is the unity operator), it projects away unphysical 
degrees of freedom. 

*It is intuitive that the propagator should not exist without gauge fixing, since there
are a number of different $A_\mu$ configurations corresponsing to the same physical state. Once you (completely) fix the gauge though, the propagator should exist, that is the time evolution of the free field should be unique.

*$\Pi^0$ does not exist because not every component of the $A_\mu(x)$ vectorpotential
corresponds to a physical degree of freedom. Nevertheless the Euler-Lagrange equations
of $A_0$ have to be satisfied, but these do not containt derivative of $A_0$, they are not dynamical equations, but constraints. 

*Once you fix a gauge you can solve for this constraint, and use the actual degrees of freedom for quantization. For example, in axial gauge $A_3=0$, this gives a new constraint. The dynamical degrees of freedom will be $A_1,A_2$ and $\Pi_1,\Pi_2$. The 
others you have to express by these by solving the contraints.
This is not just a problem with canonical quantization, if you try functional integral quantization, you run into the same problem, and  you have to fix the quage there as well.
