Can we always express the EM-Field Hamiltonian as (possibly time dependent) pair of annihilation and creation operators? The EM-Field Hamiltonian is, in principle, a functional (with a chosen operator ordering) that is defined on operator-fields $\hat{A}(x)$ and $\partial_\mu \hat{A}$. If you carry out the calculations and use definitions of $\hat{B}$ and $\hat{E}$, you'll arrive at:
$$
\hat{H} = \int d^3x \frac{\epsilon_0}{2} (\hat{\vec{E}}^2 + c^2 \hat{\vec{B}}^2)
$$
I'll take this as the Definition of the Hamiltonian in future calculations. 
For the free field, the Ansatz $\hat{\vec{A}} = \vec{e}(\hat{a_{\vec{k}}}e^{i(\vec{k}\vec{x} - \omega_{\vec{k}}t)} + \hat{a_\vec{k}}^\dagger e^{-i(\vec{k}\vec{x} - \omega_{\vec{k}}t)})$ satisfies this wave equation. Using linearity, one can superimpose all the solutions, plug them into the definition of the hamiltonian, and arrives at:
$$
\hat{H} = \sum_{\vec{k}, \vec{\lambda}} \hat{a}_{\vec{k}, \lambda}^{\dagger}\hat{a}_{\vec{k}, \lambda} \omega_{\vec{k}} \hbar
$$ 
Now my question is: Can I also use this Hamiltonian in an interacting theory? (For example an EM-Field coupled to an atom). I'm asking because the wave equations that the Heisenberg operators do change. Superimposing the creation and annihilation operators, as shown above, is no longer a solution to the field equation, so I can't express $\vec{E}$ and $\vec{B}$ no longer, just using $\hat{a}$ and $\hat{a}^{\dagger}$?
How can I still motivate $$
\hat{H} = \sum_{\vec{k}, \vec{\lambda}} \hat{a}_{\vec{k}, \lambda}^{\dagger}\hat{a}_{\vec{k}, \lambda} \omega_{\vec{k}} \hbar
$$ 
to be the right Hamiltonian?
EDIT: It is clear to me that in case of Interaction, there will be an additional Interaction Term, for example, something like 
$\hat{\vec{x}} \hat{\vec{E}} \frac{e}{\hbar}$. I'm clear of that fact. I however want to know if one can always expand quantities like $\hat{\vec{E}}$ in Terms of creation and annihilation operators. 
For example: The full Hamiltonian for an electron interacting with EM-Field would be (assuming dipole approximation):
$$
\hat{\vec{p}} \frac{1}{2m} + V(\hat{\vec{x}}) + \hat{\vec{x}} \hat{\vec{E}}(\vec{x}_{Atom}) \frac{e}{\hbar} +\int d^3x \frac{\epsilon_0}{2} (\hat{\vec{E}}^2 + c^2 \hat{\vec{B}}^2)$$
I want to know if this can in general be expressed using creation and annihilation operators instead of $\hat{\vec{E}}$ and $\hat{\vec{B}}$, for example like:
$$
\hat{\vec{p}} \frac{1}{2m} + V(\hat{\vec{x}}) + \hat{\vec{x}} \hat{\vec{E}}(\vec{x}_{Atom}) \frac{e}{\hbar} +  \sum_{\vec{k}, \vec{\lambda}} \hat{a}_{\vec{k}, \lambda}^{\dagger}\hat{a}_{\vec{k}, \lambda} \omega_{\vec{k}} \hbar$$
 A: In a Fock space, every operator can be expressed in terms of creation and annihilation operators, though usually there appear nonquadratic terms. This is the subject of second quantization, and can be made mathematically rigorous.
However, in an interacting theory, the underlying Hilbert space is not a Fock space (Haag's theorem: the interaction picture does not exist), hence a description in terms of creation and annihilation operators (which act on Fock spaces only) is no longer meaningful. The attempt to do so leads to the well-known infinities. Renormalization destroys the Fock space and with it the expression of the Hamiltonian and other generators of the Poincare group in terms of creation and annihilation operators.
Note: One can always construct the interaction picture in quantum mechanics with a finite number of classical degrees of freedom. But in relativistic quantum field theory the corresponding would-be unitary transformation does not exist because it diverges violently when the cutoff of a regularized version is removed, and unlike for the S-matrix elements, renormalization does not rescue the situation. This is called Haag's theorem. A comprehensive modern treatment is given in the PhD thesis  by Lutz Klaczynski (2016). See also 
[this PhysicsOverflow discussion] 
(https://www.physicsoverflow.org/22400).
A: OP is asking whether an expression of the form
$$
{H} = \sum_{\vec{k}, \vec{\lambda}} {a}_{\vec{k}, \lambda}^{\dagger}{a}_{\vec{k}, \lambda} \omega_{\vec{k}} 
$$
holds in an interacting QFT. The answer is, in general,
No.
For one thing, such a Hamiltonian would not, in general, be time-independent, inasmuch as
$$
\frac{\mathrm d}{\mathrm dt}a_{\vec k,\lambda}(t)\neq 0
$$
For example, for a Klein-Gordon field,
$$
\frac{\mathrm d}{\mathrm dt}a_{\vec k}(t)\sim\int \mathrm dx\ (\partial^2+m^2)\phi
$$
while for a Dirac field,
$$
\frac{\mathrm d}{\mathrm dt}a_{\vec k,\lambda}(t)\sim\int\mathrm dx\ (\not\!\partial+im)\psi
$$
neither of which vanishes in interacting theories. These formulae are very important in scattering theory and can be found, for example, in Srednicki's book on QFT (§§ 5, 41).
This proves that, in general, the interacting Hamiltonian cannot be expressed as a quadratic function of the creation and annihilation operators. The former is time-independent while the latter are not.

It bears mentioning, however, that if you regard all the objects as interacting picture operators (rather than Heisenberg), then an expansion of $H$ in terms of creation/annihilation operators is valid, but it has terms of higher order in these objects (because non-free theories include creation and annihilation phenomena, where several particles scatter into other particles). This expansion is valid, in fact, for any operator, in any theory. The proof of this claim can be found in Weinberg's book on QFT, §4.2.

Towards the edit
OP asks whether one can always expand any field in terms of creation and annihilation operators. If the field is free, or an interacting field in the interaction picture, the answer is obviously yes. If the field is interacting and in the Heisenberg picture, the answer is still yes, as I discuss in this PSE post. But there is a very important difference: in the latter case, these creation and annihilation operators will not be time-independent, in which case their interpretation as operators that create and destroy particles loses its meaning -- such an interpretation is no longer possible/consistent. This is coherent the fact that  in interacting theories there is no clear notion of particles unless we consider asymptotic (that is, free) fields.
In any case, it is very important to remark that the expression in the OP
$$
 \frac{\vec{p}^2}{2m} + V(\vec{x}) + \vec{x}\cdot \vec{E}(\vec{x}) {e} +\frac12\int (\vec{E}^2 + \vec{B}^2)\mathrm d^3\vec x
$$ 
is an effective, low-energy approximation to QED, which is only valid in the non-relativistic regime. As such, it doesn't contemplate the possibility of having dynamical photons -- there is no back-reaction from the atom to the photons. The former feel the effect of the latter but not the other way around: the photons remain oblivious to the presence of the electron and are therefore essentially free (as the gauge field for QED is non-abelian). In this sense, the expansion of $\vec E,\vec B$ in terms of creation and annihilation operators is perfectly justified: these operators are free to all practical purposes.
