The problem of (non-perurbative) quantization of interacting fields is formidable. Not only it contains an infinite number of degrees of freedom, it is nonlinear due to the interaction, thus its classical solutions do not consist of a Hilbert space, in addition, it suffers from the ultraviolet and infrared problems encountered already in perturbation theory.
In principle, to perform a quantization, one needs a space of "modes”: classical solutions of the equations of motion to expand the field with respect to which. In the free case, the modes can be chosen as wave packets. It is true that in many cases in the interacting field theory, there is a one to one correspondence of the interacting solutions to free asymptotic solutions, and this is the principle according to which QED is quantized. But in QED we know these correspondences only in perturbation theory, where the interaction is assumed to be small. In the case of strong interactions one needs to know these maps explicitly to come up with quantitative answers.
Thus in practice, one needs to resort to approximations in treating interacting electromagnetic field quantizations in the presence of charges or matter in general. There are several known approximation schemes such as:
1) Quantization of the electromagnetic field in macroscopically modeled materials such as dielectrics. For linear dielectrics, here, we still have mode expansions but not wave packets in general. In addition the gauge fixing conditions will allow for longitudinal modes in addition to the transversal ones.
2) Quantization of the electromagnetic field in the presence of charges, such as interaction with atoms in theories of radiation emission. Here the electromagnetic field is expanded in a wave packet basis, but non-perturbative approximations of the nonlinear interacting equations are applied such as the rotating wave approximation