Any ideas how to evaluate this integral over the sphere?
$$ \int_{S^2} \frac{dx \wedge dy}{z} e^{1 - (x^2 + y^2)} = \underbrace{\cdots }_{ (0,0,1)} + \underbrace{\cdots }_{ (0,0,-1)}$$
I considered using Green's Theorem or Stokes' Theorem. The sphere could be read as the boundary of the ball $B = \{ x^2 + y^2 + z^2 < 1\}$ but what could the differential form be? Perhaps we could find the divergence $\nabla \cdot F$.
Another idea could be to exploit the rotation invariance:
\begin{eqnarray*} x&=& \;\,x \cos \theta + y \sin \theta\\ y &=& -x \sin \theta + y \cos \theta \end{eqnarray*}
and the integrate over the remaining interval $[-1,1]$.
Let's check that our 2-form is defined at the poles of the hemispheres. If $z \approx 0$ let's say $z = \varepsilon \ll 1$ then $(x,y) = \sqrt{1 - \varepsilon^2} (\cos \theta, \sin \theta)$. Then using the polar coordinaes form of area $dx \wedge dy = r \, dr \wedge d\theta$:
\begin{align} \frac{dx \wedge dy}{z} & = \frac{dx \wedge dy}{\sqrt{1 - (x^2 + y^2)}} \\[10pt] & \approx \frac{d \big( (1 - \frac{1}{2}\varepsilon^2)(\cos \theta) \big)\wedge d \big( (1 - \frac{1}{2}\varepsilon^2)(\sin \theta)\big)}{\varepsilon} = \frac{\varepsilon (-2\varepsilon) \,d\varepsilon \wedge d\theta}{\varepsilon} \end{align}
And it looks like the epsilon's cancel out. A more careful check should use the Mean value theorem.
Indeed you can apply Stokes's Theorem. As @md2perpe observed, $\dfrac{dx\wedge dy}z$ is the area $2$-form $\sigma$ on the unit sphere, so it can be expressed instead as $$\sigma = x\,dy\wedge dz + y\,dz\wedge dx + z\,dx\wedge dy.$$ (This is the usual form of the area $2$-form built out of the unit normal vector field.) By the way, you can check this directly: Since $x^2+y^2+z^2=1$, on $S^2$ we have $x\,dx + y\,dy + z\,dz = 0$, so \begin{align*} z\sigma &= x\,dy\wedge z\,dz + z\,dz\wedge y\,dx + z^2\,dx\wedge dy \\ &= (x\,dy - y\,dx)\wedge (-x\,dx-y\,dy) + z^2\,dx\wedge dy \\ &= (x^2+y^2+z^2)\,dx\wedge dy = dx\wedge dy, \end{align*} as desired.
Next simplify by writing $e^{1-(x^2+y^2)} = e^{z^2}$. So, by Stokes's Theorem, $$\int_{S^2} e^{z^2}\sigma = \int_D e^{z^2}(3+2z^2)dx\wedge dy\wedge dz,$$ where $D$ is the unit ball $\{x^2+y^2+z^2\le 1\}$. We're going to end up pretty much the same place. Doing the $dx\,dy\,dz$ integral in that order, we end up with $$\int_{-1}^1 e^{z^2}(3+2z^2)\big(\pi(1-z^2)\big)\,dz.$$ With the usual integration by parts, you'll end up with some multiple of $\int_0^1 e^{t^2}\,dt$, as @md2perpe observed directly.
P.S. What do those $\dots$ on the right side of your equation mean?