I have a Lorentz model as $$\varepsilon_1 = \frac{\omega_p^2(\omega_0^2-\omega^2)}{(\omega_0^2-\omega^2)^2+\gamma^2\omega^2},$$ $$\varepsilon_2 = \frac{\omega_p^2 \gamma \omega}{(\omega_0^2-\omega^2)^2+\gamma^2\omega^2}.$$
When the frequency $\omega$ near the resonance frequency $\omega_0$, the above functions may be approximated by $$\varepsilon_1 \approx \frac{\omega_p^2(\omega_0-\omega)/2\omega_0}{(\omega_0-\omega)^2+(\gamma/2)^2},$$ $$\varepsilon_2 \approx \frac{\omega_p^2 \gamma /4\omega_0}{(\omega_0-\omega)^2+(\gamma/2)^2}.$$
I want to know how to obtain the above approximations?
Too long for a comment.
Looking at the approximations, they look like Padé approximants built around $\omega=\omega_0$ but they are not.
Considering $\varepsilon_1$, its $[1,2]$ Padé approximant is
$$\varepsilon_1=\frac{2 \,\omega_p^2}{\gamma ^2 \omega_0} ~~\frac{ (\omega_0 -\omega )}{1-\frac{3 }{2 \omega_0 }(\omega_0 -\omega )+ \frac{\left(\gamma ^2+16 \omega_0 ^2\right) }{4 \gamma ^2 \omega_0 ^2}(\omega_0 -\omega )^2}$$ which does not have much to do with what is in the post.
Considering $\varepsilon_2$, its $[0,2]$ Padé approximant is $$\varepsilon_2=\frac{\omega_p^2}{\gamma \omega_0}~~\frac 1{1-\frac 1{ \omega_0} (\omega_0-\omega)+\frac 4{\gamma^2}(\omega_0-\omega)^2}$$
What I suppose is that your writing probably hides to us some imaginary part.