The fact that using matrices for fractional linear transformations creates a group action (turns matrix multiplication into function composition) may seem like divine inspiration if projective spaces are not mentioned to motivate the situation.
Given a vector space $V$, there is an associated projective space $\mathbb{P}^1(V)$ which collects together all of the "lines," i.e. one-dimensional subspaces. When $V=\mathbb{R}^n$ or $\mathbb{C}^n$ we write $\mathbb{RP}^{n-1}$ or $\mathbb{CP}^{n-1}$. In the case that $n=2$, we have a complex projective line $\mathbb{CP}^1$.
We may put an equivalence relation on $\mathbb{C}^2\setminus0$, saying $v\sim w$ whenever $v$ and $w$ generate the same complex one-dimensional subspace (i.e. $\mathbb{C}v=\mathbb{C}w$). Since $(x,y)\sim\lambda(x,y)$ for any $\lambda\ne0$, we may scale any vector to $(1,z)$ (if its $x$ coordinate was nonzero) or $(0,1)$ (if $x=0$). Thus, one way of interpreting $\mathbb{CP}^1$ is the Riemann sphere $\widehat{\mathbb{C}}=\mathbb{C}\cup\{\infty\}$. Topologically, this is indeed a sphere, and is the one-point compactification of the complex plane $\mathbb{C}$. To speak in a more uniform manner, we can interpret $(1,\infty)$ to be the special point $(0,1)$.
Since $\mathrm{GL}_2(\mathbb{C})$'s action on $\mathbb{C}^2\setminus 0$ commutes with scalar multiplication, and the one-dimensional subspaces (equivalence relations of $\sim$) are precisely the orbits of $\mathbb{C}^\times$'s action, there is an induced action of $\mathrm{GL}_2(\mathbb{C})$ on $\mathbb{CP}^1$. Indeed, if $(x,y)\sim(1,z)$ we have
$$\begin{pmatrix} a & b \\ c & d \end{pmatrix}\cdot (x,y)=(ax+by,cx+dy)\sim\left(\frac{az+b}{cz+d},1\right) $$
Thus, using $\mathbb{C}\cup\{\infty\}$ as a set of representatives for the elements of $\mathbb{CP}^1$, the induced action of $\mathrm{GL}_2(\mathbb{C})$ is by linear fractional transformations.
Stereographic projection creates a diffeomorphism $\widehat{\mathbb{C}}\simeq S^2\subset\mathbb{R}^3$. Explicitly, think of $\mathbb{C}\times\mathbb{R}$ as three-space $\mathbb{R}^3$. The line through the north pole $(0,1)$ and a complex number $(z,0)$ is parametrized by $t(z,0)+(1-t)(0,1)=(tz,1-t)$. In order for a point on the line to be on the unit sphere centered at the origin, we require $t^2|z|^2+(1-t)^2=1$ which after some algebraic manipulation implies that $t=2/(|z|^2+1)$, and therefore the point is
$$ z\mapsto \left(\frac{2z}{|z|^2+1},\frac{|z|^2-1}{|z|^2+1}\right).$$
This creates a map $\mathbb{C}\to S^2$. We can extend this to $\widehat{\mathbb{C}}$ if we can define where to send $\infty$, which can be done by taking the limit $z\to\infty$, yielding the north pole $(0,1)$.
Thus, we can transport the action of $\mathrm{GL}_2(\mathbb{C})$ over to $S^2$, but not every such matrix will act in a way that comes from rotating $S^2$ within $\mathbb{R}^3$. In fact, it will do so precisely when it acts by an isometry on the sphere $S^2$. (Rotations restrict to isometries of course, and it is a standard exercise to prove isometries come from orthogonal transformations, but $\mathrm{GL}_2(\mathbb{C})\to\mathrm{O}(3)$ must land in the connected component $\mathrm{SO}(3)$ since $\mathrm{GL}_2(\mathbb{C})$ is connected.)
Over at Qiaochu's question here I provide an elementary explanation of why matrices from $\mathrm{SU}(2)$ do indeed induce rotations on the space $\mathbb{R}^3$ via stereographic projection. I also mention how the three stabilizer subgroups of $\mathrm{SU}(2)$ associated to the special points $0,1,i$ correspond to one-parameter subgroups of coordinate axis rotations in $\mathbb{R}^3$ (and also to Pauli matrices), which implies the map $\mathrm{SU}(2)\to \mathrm{SO}(3)$ is surjective, and also mention the kernel is $\{\pm I\}$.
To sum up, have a map $F:\mathrm{GL}_2(\mathbb{C})\to\mathrm{Diff}(S^2)$ which restricts to a map $\mathrm{SU}(2)\to\mathrm{SO}(3)$ (note that $\mathrm{SO}(3)$ can be identified with a subgroup of $\mathrm{Diff}(S^2)$ due to the earlier remark about isometries of $S^2$ corresponding to transformations of $\mathbb{R}^3$).
Say we have two solutions $\alpha,\beta$ so to the short exact sequence
$$ 1\to K\to G\to H\to 1. $$
That is, $\alpha,\beta:G\to H$ are onto Lie group homomorphisms with kernel $K$ (a Lie subgroup). Then they are related by an automorphism $\gamma$ of $H$. In particular, consider $\overline{\alpha},\overline{\beta}:G/K\to H$ (afforded by the first isomorphism theorem) and define $\gamma=\overline{\beta}\circ\overline{\alpha}^{-1}$, so that $\beta=\gamma\circ\alpha$.
In the case of $H=\mathrm{SO}(3)$, every automorphism $\gamma$ is inner, i.e. is conjugation $\gamma_g(h)=ghg^{-1}$ for some rotation $g\in\mathrm{SO}(3)$. To see this, note that any automorphism $\gamma$ of $\mathrm{SO}(3)$ induces one of its lie algebra $\mathfrak{so}(3)\cong(\mathbb{R}^3,\times)$ which must be by a rotation $g$ since $\mathrm{SO}(3)$ is the symmetry group of the three-dimensional cross product. Note that $\mathfrak{so}(3)\cong(\mathbb{R}^3,\times)$ is not just an isomorphism of lie algebras but also $\mathrm{SO}(3)$-representations, so $\gamma$ and $g$ agree acting on $\mathfrak{so}(3)$ (where the adjoint action is used). Then since the matrix exponential $\exp:\mathfrak{so}(3)\to\mathrm{SO}(3)$ is surjective (as $\mathrm{SO}(3)$ is compact, but it can be proved as an elementary exercise too), their actions can be lifted to $\gamma$ and conjugateion-by-$g$ on $\mathrm{SO}(3)$, which must agree too.
(I'm sure there's a better way than this, but this is what comes to mind.)
By the way, other matrices in $G=\mathrm{GL}_2(\mathbb{C})$ can be interpreted as acting on $S^2\subset\mathbb{R}^3$ too, but we require some more freedom. In particular, we have an Iwasawa decomposition $G=KAN$ with subgroups $K=\mathrm{SU}(2)$, $A$ the positive real diagonal matrices, and $N$ the unitriangular complex matrices. We've already seen the interpretation of $\mathrm{SU}(2)$ acting, $A$ acts on $\mathbb{C}$ via scalings / dilations / homotheties, which can be interpreted as the restriction of corresponding scalings of $\mathbb{R}^3$, and finally $N$ acts on $\mathbb{C}$ by translations, which can be interpreted as restrictions of planar translations acting on $\mathbb{R}^3$ (remember we're identifying $\mathbb{R}^3$ and $\mathbb{C}\times\mathbb{R}$).
This combo of rotations, translations and scalings is the full menu of transformations seen in the Möbius Transformations Revealed video, in case anybody wondered.
Best Answer
I became curious about what that picture means. Here is my speculation: given a Mobius transformation as a matrix $ G \in SL(2,\mathbb{C}) $, it corresponds (by logarithm) to an element $ g \in \mathfrak{sl}_2 $, which then gives a vector field on $ \mathbb{R}^2 $ (by $ \frac{d}{dt} e^{g t} \big|_{t = 0}) $
Here are details for an easy example. Let $$ g= \begin{pmatrix} 0 & w \\ w & 0 \end{pmatrix} $$ where $ w \in \mathbb{C} $. He is in $ \mathfrak{sl}_2 $ because his trace is zero. Exponentiating: $$ e^{ t g} = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} + t \begin{pmatrix} 0 & w \\ w & 0 \end{pmatrix} + O(t^2) $$ This corresponds to the Mobius transformation: \begin{align} z &= \frac{z + t w }{z t w + 1} + O(t^2)\\ &= z + t \; w (1 - z^2) + O(t^2) \end{align} The $ O(t) $ term is the vector field. He has zeros at $ z = \pm 1 $, which means $ G = \exp(g) $ will have fixed points at $ z = \pm 1 $. With different choices of $ w $, $ g $ upon exponentation gives all Mobius transformations with those fixed points.
Writing $ w = a + b i $ and $ z = x + i y $, and taking the real and imaginary parts to get the $ \hat{x} $ and $ \hat{y} $ components: $$ w ( 1 - z^2 ) \sim \left( 2 b x y + a( 1 - x^2 + y^2) \right) \hat{x} + \left( - 2 a x y + b(1 - x^2 + y^2) \right) \hat{y} $$
If I set $ b = 0 $, then I get pictures like the one you posted. With $ a = 0 $, I get pictures that look like the B-field lines of two wires.
You can check that the divergence and curl of are nonzero, so it won't be an E field or B field on the nose.
However, you can also check that after rescaling the vector field by a function, then it is the electric field of two point charges. So the stream lines are indeed the same!
EDIT: The electric field of one point charge (in 2d!) can be written as: $$ \frac{( x - x_0) \hat{x} + (y - y_0) \hat{y} }{(x-x_0)^2+(y-y_0)^2} $$ Putting one charge at $ (1,0) $ and another charge at $ (-1,0) $, I get: $$ \frac{( 2 - 2 x^2 + 2 y^2) \hat{x} + (4 x y ) \hat{y} }{x^4 + 2 x^2 (-1 + y^2) + (1 + y^2)^2}$$ So you see, this is the same as the vector field from above with $ a = 0, b = 2 $ multiplied by the scalar function $ 1 / \; ( x^4 + 2 x^2 (-1 + y^2) + (1 + y^2)^2 ) $.