Why should field operators satisfy the classical equations of motion? To quantize a scalar field theory with the action:
$$S=\int \mathrm d^Dx\mathscr{L}(\phi,\partial_\mu\phi)=\int \mathrm dx^0L(\phi,\partial_0\phi)$$
we promote $\phi(\vec{x})$ and $\pi=\frac{\delta L}{\delta(\partial_0\phi)}=\pi(\vec{x})$ to be field operator at fixed time $x^0$ (with the canonical commutation relation $[\phi(\vec{x}),\pi(\vec{y})]=\delta^{D-1}(\vec{x}-\vec{y})$), then we use the classical equation of motion to find $\phi$ and $\pi$ at all time.
What principles lie behind this? Why can we assume that the field operators satisfy the classical equation of motion?

By example if we have a classical theory with $S=\int \,\mathrm dtL(\phi,\partial_t\phi)$, momentum $\pi=\frac{\delta L}{\delta(\partial_t\phi)}$ hamiltonian $H=\pi\phi-L$, the classical equation of motion is:
$$\frac{\mathrm d\phi}{\mathrm dt}=\frac{\delta L}{\delta(\partial_t\phi)}$$
When we do the canonical quantization, we promote $[\phi,\pi]=\delta(\vec{x}-\vec{y})$ and the equation of motion for $\phi$ is $$\frac{\mathrm d\phi}{\mathrm dt}=i[H,\phi]$$
Why should $i[H,\phi]$ has the same form as $\frac{\delta L}{\delta(\partial_t\phi)}$? This fact comes from where?
 A: Suppose first that you deal with free field theories.
You start from canonical commutation or anticommutation relations for coordinate $q(x, t)$ and momentum $p(x, t)$:
$$
\tag 0 [q_{i}(\mathbf x , t), p_{j}(\mathbf y , t)]_{\pm} = i\delta (\mathbf x - \mathbf y)\delta_{ij}, \quad [q_{i}(\mathbf x , t),q_{j}(\mathbf y , t)]_{\pm} = [p_{i}(\mathbf x , t), p_{j}(\mathbf y , t)]_{\pm} = 0
$$
This is in some sense quantum analog of Poisson brackets. 
Such relations reault in important variational identities (since here operators are coordinate dependent, we need to use variational calculus instead of differential one) for bosonic operator $F \equiv F(q, p)$:
$$
\frac{\delta F[q(t),p(t)]}{\delta q_{n}(\mathbf x , t)} \equiv i[p_{n}(\mathbf x , t),F[q(t), p(t)]], \quad \frac{\delta F[q(t),p(t)]}{\delta p_{n}(\mathbf x , t)} \equiv i[F[q(t), p(t)], q_{n}(\mathbf x , t)]
$$
So that for arbitrary $c$-valued functional $F[q(t), p(t)]$ its variation is
$$
\tag 1 \delta F [q(t), p(t)] = 
$$
$$
\sum_{n}\int d^{d}\mathbf r \left( \delta q_{n}(\mathbf r , t)i[p_{n}(\mathbf r , t),F[q(t), p(t)]] + \delta p_{n}(\mathbf r , t)i[F[q(t), p(t)], q_{n}(\mathbf r , t)]\right)
$$
From the mathematical sense of $\hat{H}_{0}$ (it is time translation operator) follows that
$$
\tag 2 q_{n}(\mathbf r, t) = e^{iH_{0}t}q_{n}(\mathbf r , 0)e^{-iH_{0}t}, \quad p_{n}(\mathbf r, t)e^{iH_{0}t}p_{n}(\mathbf r , 0)e^{-iH_{0}t}
$$
From $(1)$ and $(2)$ follows canonical Heisenberg EOM:
$$
\tag 2\frac{dF[q(t), p(t)]}{dt} = i[F, H_{0}]
$$
Now, let's assume case of interacting theory. We may interact $Q(\mathbf x , t), P(\mathbf x , t)$ from free field coordinates and momentum by using Heisenberg picture:
$$
\tag 3 Q(\mathbf r , t) = e^{iHt}q(\mathbf r , 0)e^{-iHt}, \quad P(\mathbf r , t) = e^{iHt}p(\mathbf r , 0)e^{-iHt},
$$
where $H$ is the full hamiltonian. Since this is similarity transformation (and it commutes with $H$), then the full hamiltonian is a functional of $Q, P$ as well as it was functional of $q, p$:
$$
H[Q, P] = e^{iHt}H[q, p]e^{-iHt} = H[q, p]
$$ 
Next, due to the same reason $Q, P$ obeys the same canonical relations $(0)$ and EOM $(2)$ as free theory $q, p$, but with $H$ instead of $H_{0}$.
A: How to see it in canonical quantization: All operators $\mathcal{O}$ in a quantum theory fulfill the Heisenberg equations of motion
$$ \frac{\mathrm{d}}{\mathrm{d}t}\mathcal{O}(t) = \mathrm{i}[H,\mathcal{O}(t)]$$
where $H$ is the Hamiltonian density and which is exactly the quantum version of the classical Hamiltonian equations of motion. So the fields, as they are operators, indeed must obey the classical equations of motion.
How to see it in path integral quantization: Write $\phi'(x) = \phi(x) + \epsilon\delta(y-x)$ and observe the path integral measure is invariant under this. Expand the integrand to first order in $\epsilon$, and deduce
$$ \int\left(\frac{\delta S}{\delta\phi(x)} + J(x)\right)\mathrm{e}^{\mathrm{i}S[\phi]+J\phi}\mathcal{D}\phi = 0$$
which is known as the Schwinger-Dyson equation. Setting $J=0$ gives $\delta S/\delta\phi = 0$ (inside the path integral, which is the PI version of "as an operator equation"), which is exactly the classical equation of motion.
