Stream functions and divergence? 
We see that the existence of a stream function guarantees that the vector field has zero 
  divergence or, equivalently, is source free. The converse is also true on simply connected 
  regions of $\mathbb{R}^2$

Why does the region have to be simply connected on $\mathbb{R}^2$ for the converse to be true?
If the vector field has zero divergence,
$$\frac{\partial{f}}{\partial{x}} + \frac{\partial{g}}{\partial{y}} = 0$$
$$\frac{\partial{f}}{\partial{x}} = -\frac{\partial{g}}{\partial{y}}$$
Integrating either side with respect to the appropriate variable would give you the stream line function. Doesn't this show that if a vector field has zero divergence it has a stream function no matter what?
 A: The problem with non-simply-connected domains is that the result of integration may depend on how the path of integration goes around the holes. A better known manifestation of this phenomenon is the absence of a potential function for some irrotational fields.
The standard example for your situation is the radial flow with velocity 
$$\vec V(x,y) = \frac{x}{x^2+y^2}\vec\imath+ \frac{y}{x^2+y^2}\vec\jmath$$
more compactly written as $\vec V=\vec r/|\vec r|^2$.
This is a divergence free flow in $\mathbb R^2\setminus \{(0,0)\}$. (Note that the flux through a circle centered at the origin is $2\pi$ regardless of its radius.)  Yet, there is no streamline function for this flow, for such a function would have a gradient pointing in the tangential direction, which leads to a contradiction when we come back to the original position after following the gradient. 
A: This being said, there is a result in non simply connected domain. If
$$ \textrm{div}D=0,$$ and on the boundary $\Gamma_i$ of each connected component of $\Omega$ there holds $$\int_{\Gamma_i} D\cdot n \textrm{d}\sigma =0,$$ then there exists a function $H$ such that $$D=(-\partial_{x_2} H,\partial_{x_1} H).$$
See Lemma I.1 in Brezis, Bethuel & Hélein's book. The proof is that you can solve the Neumann problem on each hole, and then supplement $D$ with the gradient of the solutions, obtaining a divergence free problem on the corresponding simply connected set.
