Lorentz Dispersion model: why complex? The differential equation behind the Lorentz Dispersion model is this one:
$$ \ddot x+\omega_S\dot x+\omega_0^2 x=\frac{qE_{0}}{m}\exp(i\omega_F t) $$
where $x(t)$ is the displacement of the electron from its equilibrium position.
Using the Fourier transform method the solution is easily found:
$$x(t)=\frac{qE_{0}}{m}\frac{\exp(i\omega_F t)}{\omega_0^2-\omega_F^2-i\omega_S\omega_F}$$
That is, of course, complex.
Since $$D(t)=\varepsilon_0E(t)+P(t)=\varepsilon_0E(t)+Nqx(t)=\varepsilon E(t) \tag{1}$$ this leads to the known Lorentz Model:
$$ \frac{\varepsilon}{\varepsilon_0}=1+\frac{Nq^2}{m\varepsilon_0}\frac{1}{\omega_0^2-\omega_F^2-i\omega_S\omega_F}$$
Since this is a complex quantity, it has an imaginary part that can be used to justify the fact that materials absorb energy from electromagnetic waves.
Indeed: $k=\omega\sqrt{\varepsilon\mu}=\beta +i\frac{\alpha}{2}$. So
$$E_0 e^{ikx}=E_0 e^{i\beta x}e^{-\alpha x/2} \tag{2}$$
But every single textbook says that
$E_0\exp(i\omega_F t)$ is just a short way to write $E(t)=\mathbb{Re}[E_0\exp(i\omega_F t)]=E_0\cos(\omega_F t)$. Writing again the equation:
$$ \ddot x+\omega_S\dot x+\omega_0^2 x=\frac{qE_{0}}{m}\cos(\omega_F t) $$
This time standard resolution methods lead to a real solution
$$ x(t)=Ae^{-{\omega_S\over 2}t}\cos(\omega_0 t+\phi)+B\cos(\omega_F t-\delta) $$
That will, of course, lead to a real dielectric function $\varepsilon(\omega)$. So the explanation of the absorption is now missing. How do I explain that?
 A: If you phrase the equation of motion as explicitly real valued,
$$ \ddot x+\omega_S\dot x+\omega_0^2 x=F(t)=\frac{qE_{0}}{m}\cos(\omega_F t), $$
then indeed you can write down the solution as
$$ x(t)=Ae^{-{\omega_S\over 2}t}\cos(\omega_0 t+\phi)+B\cos(\omega_F t-\delta) $$
but because you are only interested in the steady-state behaviour you can reduce this to the form
$$ x(t)=B\cos(\omega_F t-\delta).$$
This is completely correct and entirely equivalent to the complex-$x(t)$ phasor formulation. In particular, this solution does encode the absorption: it is hardwired into the phase term $\boldsymbol \delta$.
This is required by the correspondence with the phasor model, since it is easy to derive from its formulation that $\delta$ will vanish if and only if the damping $\omega_S$ is eliminated. However, it is probably more illuminating to tackle this directly, which you do by computing the instantaneous and average power delivered to the charge by the external forcing. This comes through the velocity of the particle,
$$
\dot x(t) = -\omega_F B\sin(\omega_F t-\delta),
$$
and it is essentially of the form
$$
P(t) = \dot x(t) \cdot F(t) = -\frac{\omega_F qE_0}{m}B \sin(\omega_F t-\delta)\cos(\omega_F t).
$$
However, in this form you get an awkward product, with a convoluted interplay of the field with the dephasing, so it's much cleaner to split things up a bit by cracking open the trigonometrics as
\begin{align}
P(t) &= -\frac{\omega_F qE_0}{m}\left[B \cos(\delta)\sin(\omega_F t)\cos(\omega_F t)-B\sin(\delta) \cos^2(\omega_F t)\right].
\end{align}
This looks a bit messy, but you have now reduced the instantaneous power delivered to the charge into two very different contributions:


*

*One is the first term, with coefficient $B \cos(\delta)$ and varying as $\sin(\omega_F t)\cos(\omega_F t)$, which takes both positive and negative values (i.e. power is delivered to, then taken from, the charge, alternatingly) and averages exactly to zero. This is known in circuit parlance as 'reactive power', and it reflects a back-and-forth exchange of energy between the field and the matter (much like capacitors and inductors bounce energy back and forth in an RLC circuit) without any net transfer of energy and hence without any absorption.

*The second term, on the other hand, varies as $\cos^2(\omega_F t)$ and it is always positive, so it will average out to $\frac12$ over each cycle. This directly models the absorption, and not coincidentally it's proportional to $B\sin(\delta)$, which vanishes at $\delta=0$.
Alternatively, you could take a more macroscopic-focused approach, and try and see things through the electric polarization $P_\mathrm{dip}(t)$ generated by this motion of the charges, and here the absorption again shows up as intrinsically tied to the phase shift between the electric field and the electric polarization: waves can propagate without losses only so far as both vectors oscillate completely in sync; if the polarization lags behind the field then you can show that the amplitude will diminish.
With that in mind, you can indeed take the real-valued perspective and say that your forcing as $\cos(\omega_Ft)$ introduces a response in the form
$$
F(t) = F_0 \cos(\omega_Ft) \ \rightarrow \ x(t) = B\cos(\omega_F t-\delta).
$$
However, if you then simply define the "permittivity" $\varepsilon(\omega_F)$ as simply the amplitude ratio $B/F_0$, then yes, you lose information - but only because you're explicitly discarding it when you try to interpret the mathematics. If you were doing things correctly, you'd need to describe the response function as an ordered pair $(B/F_0,\delta)$ that included both the amplitude and the phase, but if you did this you'd discover that it is simply an awkward mapping of the complex number $\frac{B}{F_0}e^{i\delta}$ which directly corresponds with the complex permittivity from the phasor method. And, as usual, at that point the wise thing to do is to just give up and embrace the complex-valuedness of those quantities (while keeping mind that the physical part is still just the real part of the quantities involved.
A: 
That will, of course, lead to a real dielectric function $\varepsilon(\omega)$. So the explanation of the absorption is now missing. How do I explain that?

The dielectric function $\varepsilon(\omega)$ that has non-zero imaginary part was introduced to relate complex phasors $\tilde{P}(t)$ and $\tilde{E}(t)$ in linear media. The complex $\varepsilon$ makes sense only when working with such complex quantities. It gives magnitude and angle of $\tilde{P}(t)$ with respect to $\tilde E(t)$ according to the formula
$$
\tilde{P}(t)=(\epsilon(\omega) - \epsilon_0) \tilde{E}(t).
$$
If the electric field and polarization are not represented by such complex phasors, but by real functions of time instead, there is no such function $\varepsilon(\omega)$ that would allow us to relate these functions by
$$
P(t)=(\epsilon(\omega) - \epsilon_0) E(t).
$$
This is because for harmonic functions $P(t)$ and $E(t)$ these two are not usually in phase, but there is certain phase lag between the two oscillations. It is this phase difference that indicates there is absorption or emission taking place in the medium. This can be verified by taking into account this phase delay when solving Maxwell's equations - it can be shown it implies that amplitude of electric field in the medium decreases along the direction of propagation of the EM wave.
Complex representation just expresses the phase relation of the two functions differently and allows easy algebraic manipulations. But it is not necessary; all results can be obtained via real quantities, if we wish. This may take a little more time as it may be necessary to recall and employ formulae for sin(x+y) or cos(x+y).
A: There are a few unrelated issues. First, you've assumed that if a function is real, then its Fourier transform is also real. But the definition of the Fourier transform is
$$\widetilde{f}(\omega) = \int f(t) e^{-i\omega t} dt$$
which explicitly contains complex numbers; $\widetilde{f}(\omega)$ will generically be complex even if $f(t)$ is real. So your final sentence isn't right.
You might protest that the Fourier transform can be written in terms of purely real quantities. For example, we can decompose $f(t)$ in terms of sines and cosines. But this is clunky because now the Fourier transform of $f(t)$ is two functions, one for sine and one for cosine. So the complex exponential is a very useful convention.
The real issue is that you've got the reasoning backwards. We don't say the electron absorbs energy because $\epsilon(\omega)$ has an imaginary part. It absorbs energy simply because it does a damped oscillation, as you can see without doing any Fourier transforms at all! Then you can show that, if you use the common Fourier transform convention above, then the fact that $f(t)$ exponentially decays translates to $\widetilde{f}(\omega)$ having an imaginary part.
A: By using the Fourier transform, one is implicitly restricting to the solution(s) that can be Fourier transformed. So the transient part of the general solution, i.e., $A e^{-\omega_S t / 2} \cos(\omega_{0}t + \phi)$ with $A$ and $\phi$ being arbitrary constants, is automatically discarded.
Discarding this term is perfectly fine as long as we are interested in the steady-state behavior.

Let me add more because the point OP intended to raise is different from what is addressed on the above.
The imaginary part of the proportionality constant between $x(t)$ and $E(t)$ is associated with damping only in the context of a time-harmonic source $E(t) = E_0 e^{i\omega_{F} t}$. With this kind of $E(t)$, $d/dt$ is equivalent to multiplying $i\omega_{F}$, and, therefore, the relation between $E(t)$ and $x(t)$ is linear. Also, a first-order derivative with a real coefficient in the equation of motion, which corresponds to damping, gives rise to an imaginary part in the proportionality constant.
The complex exponential $e^{i\omega_{F} t}$ is perfectly fine as a source function, because we can construct any $E(t)$ as a superposition of complex exponentials via the Fourier transform. [To be precise, the integral of $|E(t)|$ over the real numbers should be finite for the Fourier transform to be well-defined, but this condition is satisfied by any reasonable driving field.] The restriction that $E(t)$ is a real function simply leads to the relation between the transformed quantities, i.e., $\tilde{E}(\omega)^{\ast} = \tilde{E}(-\omega)$.
If $E(t) \propto \cos(\omega_{F} t)$, the resulting solutions simply doesn't have the nice properties described above. In this case, in the first place, the relation between $E(t)$ and $x(t)$ is not linear because $x(t)$ is proportional to $\cos(\omega_{F}  t - \delta)$. Of course, one can define the proportionality constant between the amplitudes of these cosine functions, but this constant being real has nothing to do with lack of damping. As stated above, such a connection can be made only if we are considering the propertionality constant arising from the harmonic source $E(t) = E_0 e^{i\omega_{F} t}$.
