Hamiltonian formalism of the massive vector field I am currently working through a problem concerning the massive vector field. Amongst other things I have already calculated the equations of motion from the Lagrangian density $$\mathcal{L} = - \frac{1}{4} F_{\mu\nu} F^{\mu\nu} + \frac{1}{2} m^2 A^\mu A_\mu,$$ where $F_{\mu\nu} = \partial_\mu A_\nu - \partial_\nu A_\mu$, which are
\begin{align}
\partial_\mu F^{\mu\nu} + m^2 A^\nu = 0. 
\end{align}
Here the sign convention is $(+,-,-,-)$.
The problem then leads me through some calculations to end up with a Hamiltonian. Basically one defines the canonical momentum and from the equations of motion it follows that $A^0 = \frac{1}{m^2} \partial_i \Pi_i $ (where from here on summation convention is used for repeated indices regardless of their position). Basically this means that $A^0$ is not a dynamical variable and can be eliminated in terms of $\Pi_i$. By using this and the fact that $\Pi_i (x) = \partial_0 A^i (x) + \partial_i A^0 (x)$, one can find the following Hamiltonian:
\begin{align}
H = \int d^3 \vec{x}\;  \mathcal{H} = \int d^3 \vec{x}\; \left(\frac{1}{2} \Pi_i \Pi_i + \frac{1}{2m^2} \partial _ i \Pi_i \partial _j\Pi_j + \frac{1}{2} \partial_i A^j (\partial_i A^j - \partial_j A^i ) + \frac{m^2}{2} A^i A^i \right).
\end{align}
Long story short: I am now supposed to calculate the Hamiltonian equations of motion from this and show that they lead to the same ones that I got from the Lagrangian.
Now it is not clear to me what form the Hamiltonian equations of motion should have here. The way they are written on Wikipedia (https://en.wikipedia.org/wiki/Hamiltonian_field_theory) with only the time derivatives on the left hand side won't lead to the same equations of motion, right?
EDIT: Thanks to the answer by GRrocks I think I got it now.
\begin{align}
-\partial_0 \Pi^k & = - \partial_0 \left(\partial_0 A^k + \partial_k A^0  \right) =\frac{\delta \mathcal{H}}{\delta A^k} =  \\
&= m^2 A^k - \frac{1}{2} \partial_i \partial_i A^k - \frac{1}{2}\partial_i \partial_i A^k + \frac{1}{2} \partial_j\partial_k A^j - \frac{1}{2} \partial_j\partial_k A^j = \\
&= m^2 A^k - \partial_i \partial_i A^k + \partial_j\partial_k A^j 
\end{align}
and so
\begin{align}
\partial_0 \partial_0 A^k - \partial_i \partial_i A^k + m^2 A^k + \partial_0 \partial_k A^0 + \partial_i \partial_k A^i = 0
\end{align}
which is indeed equal to the Lagrangian equations of motion. My question now is what the equations $\partial_0 A ^i = \frac{\delta \mathcal{H}}{\delta \Pi_i} $ are if I already get the full Lagrangian equations of motion with $-\partial_0 \Pi^k  =\frac{\delta \mathcal{H}}{\delta A^k} $. What am I missing?
 A: Hint: Hamilton's equations of motion here are exactly the same as they are in classical mechanics, with ordinary derivatives replaced by functional derivatives.
This is because in general, the Hamiltonian(not the hamiltonian density) $H(t)=H[\psi(\cdot,t),\dot{\psi}(\cdot,t)]$ is a functional of the fields and conjugate momenta at a given time slice, and at that time slice, the fields and momenta obey the poisson bracket(read:commutator) relation familiar from classical mechanics(where the $H=H(q,p)$ is just a function). These coordinates $q,p$ are promoted to fields in QFT, and thus  derivatives wrt them become functional derivatives.
So, just take the functional derivatives of the hamiltonian you've written down, and put them in the functional derivative version of classical equations($\partial H/\partial p=\dot{q}$ etc)
A: There is no need to eliminate the $A_0$ field$^1$. The long story short is that the Hamiltonian Lagrangian density$^2$
$$\begin{array}{ccc} {\cal L}_H~=~\vec{\Pi}\cdot\dot{\vec{A}} - {\cal H}&\stackrel{\vec{\Pi}}{\longrightarrow} & {\cal L}~=~-\frac{1}{4}F_{\mu\nu}F^{\mu\nu}-\frac{1}{2}m^2 A_{\mu}A^{\mu} \cr\cr
\downarrow A_0& &\downarrow A_0\cr\cr
{\cal L}^R_H~=~\vec{\Pi}\cdot\dot{\vec{A}} - {\cal H}^R&\stackrel{\vec{\Pi}}{\longrightarrow} & {\cal L}^R~=~\frac{1}{2}\dot{A}_i\left(\delta^{ij}+\frac{\partial^i \partial^j}{m^2-\nabla^2} \right)\dot{A}_j-\frac{1}{2}\vec{B}^2-\frac{1}{2}m^2\vec{A}^2
\end{array} \tag{1} $$
for the real Proca theory reduces to its Lagrangian counterpart (up to total derivative terms) when one integrates out/eliminates the momenta $\vec{\Pi}$. Therefore the Hamiltonian & Lagrangian EOMs must agree, even after elimination of the $A_0$ field. Moreover the diagram (1) commutes since the order of elimination shouldn't matter.
In eq. (1) the Hamiltonian density is
$$\begin{align} {\cal H}~=~&\frac{1}{2}\vec{\Pi}^2 +\frac{1}{2}\vec{B}^2+\frac{1}{2}m^2 A_{\mu}A^{\mu}-A_0 \vec{\nabla}\cdot\vec{\Pi} \cr\cr
 &\downarrow A_0\cr\cr 
{\cal H}^R~=~&\frac{1}{2}\Pi^i\left(\delta_{ij}-\frac{\partial_i \partial_j}{m^2} \right)\Pi^j +\frac{1}{2}\vec{B}^2+\frac{1}{2}m^2 \vec{A}^2,\end{align}\tag{2}$$
and the magnetic field is
$$ B_i~=~\frac{1}{2}\epsilon_{ijk}F_{jk}, \qquad \vec{B}^2~=~\frac{1}{2}F_{ij}F_{ij} . \tag{3}$$
--
$^1$ If one integrates out/eliminates $A_0$, one would no longer be able to derive its EOM
$$A_0~\approx~-\frac{1}{m^2}\vec{\nabla}\cdot\vec{\Pi} .\tag{4}$$
$^2$ NB. This answer uses the opposite sign convention $(-,+,+,+)$ so that the position of spatial indices doesn't matter.
