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$$