I'm going to assume that the ellipse has the equation
$$\frac{x^2}{a^2} + \frac{y^2}{b^2} = 1$$
since that's the more standard assumption. Yours has $a$ and $b$ reversed and I'm not sure if you meant it that way or if that was a typo on your part (I think it was a typo since your expression for the eccentricity matches the standard one). I'll also assume $a \ge b$. The case $b \ge a$ is handled similarly.
Now imagine that you rotate the ellipse around its long axis, that is, around the $x$ axis, and focus your attention on the strip that results from rotating a small arc of the ellipse, located at $(x,y)$ before the rotation. That elliptical arc has a length $ds$ that depends on its location, so it's a function of $(x,y)$. In fact, it's not hard to show that that arc-length is given by
$$ds = dx\,\sqrt{1 + (\frac{dy}{dx})^2}$$
Anyway, the strip resulting from the rotation of that little elliptical arc has a circular shape and, therefore, an area approximately equal to
$$dA = 2\pi\,\mbox{radius} \times ds$$
and you can see from the figure that the radius is just $y$, so
$$dA = 2\pi\,y\,ds$$
Now, given the equation at the top, we find
$$\frac{2x\,dx}{a^2} + \frac{2y\,dy}{b^2} = 0$$
so
$$\frac{dy}{dx} = -\frac{b^2}{a^2}\,\frac{x}{y}$$
and
$$ds = \frac{1}{a^2}\,\frac{dx}{y}\,\sqrt{b^4x^2 + a^4y^2}$$
As promised, $ds$ depends on $(x,y)$. Putting all of the above together, we find
$$dA = dx\,\frac{2\pi}{a^2}\,\sqrt{b^4x^2 + a^4y^2}$$
The area of the entire surface of revolution is then twice the integral of the above expression, from $x=0$ to $x=a$. Twice because we're integrating over only half the ellipse:
$$A = \frac{4\pi}{a^2}\int_{x\,=\,0}^{x\,=\,a} \sqrt{b^4x^2 + a^4y^2}\,dx$$
We still need to eliminate $y$, but that's easy. From the equation at the top, we find
$$y^2 = b^2 - \frac{b^2}{a^2}\,x^2$$
and then:
$$A = 4\pi\,\frac{b}{a}\int_{x\,=\,0}^{x\,=\,a} \sqrt{a^2 - (\frac{a^2 - b^2}{a^2})\,x^2}\,dx$$
The quantity
$$\frac{a^2 - b^2}{a^2}$$
is none other than the ellipse's eccentricity $\varepsilon$. So, finally, we have
$$A = 4\pi\,\frac{b}{a}\int_{x\,=\,0}^{x\,=\,a} \sqrt{a^2 - \varepsilon^2x^2}\,dx$$
Now use the parametrisation $x = a\,\sin\theta$ (Why $\sin$ instead of $\cos$? Because it makes the math easier down below. Shouldn't it be $\cos$, though? Not necessarily. Note that $x$ is now a dummy integration variable and we can choose any substitution we want) to get
$$A = 4\pi\,ab\,\int_{\theta\,=\,0}^{\theta\,=\,\pi/2} \sqrt{1 - \varepsilon^2\sin^2\theta}\,\cos\theta\,d\theta$$
Next set $\sin\phi = \varepsilon\sin\theta$ so $\cos\phi\,d\phi = \varepsilon\cos\theta\,d\theta$ and
$$A = 4\pi\,\frac{ab}{\varepsilon}\,\int\cos^2\phi\,d\phi$$
(I omitted the integration limits but will get back to them below)
To integrate $\cos^2\phi$, we can use the fact that $\cos(2\phi) = \cos^2\phi - \sin^2\phi = 2\cos^2\phi - 1$. Thus,
$$\cos^2\phi = \frac{1 + \cos(2\phi)}{2}$$
and
$$\int\cos^2\phi\,d\phi = \int\frac{1 + \cos(2\phi)}{2}\,d\phi =
\frac{\phi}{2} + \frac{\sin(2\phi)}{4}
$$
Now back to the integration limits. Note that $\theta = 0$ implies $\sin\phi = 0$, thus $\phi = 0$, and $\theta = \pi/2$ implies $\sin\phi = \varepsilon$, that is, $\phi = \arcsin(\varepsilon)$. Note also that $0 \le \varepsilon \le 1$ since $a \ge b$.
So then we get
$$A = 2\pi\,\frac{ab}{\varepsilon}\,(\phi + \frac{\sin(2\phi)}{2})\,\big|_{0}^{\arcsin(\varepsilon)} =
2\pi\,\frac{ab}{\varepsilon}\,\big[\,\arcsin(\varepsilon) + \frac{\sin(2\arcsin(\varepsilon))}{2}\,\big]
$$
Then, using $\sin(2\phi) = 2\sin\phi\cos\phi$, we find
$$A = 2\pi\,\frac{ab}{\varepsilon}\,\big(\,\arcsin(\varepsilon) + \varepsilon\sqrt{1-\varepsilon^2}\,\big)
$$
Finally, using the definition of the eccentricity, we get
$$A =
2\pi\,\frac{ab}{\varepsilon}\,\big(\,\arcsin(\varepsilon) + \varepsilon\,\frac{b}{a}\,\big) =
2\pi\,b^2\,\big(1 + \frac{a}{b}\,\frac{\arcsin(\varepsilon)}{\varepsilon} \,\big)
$$
which is the expression you wanted to prove.
