I. The argument using normal families: $\DeclareMathOperator{\re}{Re}$
We let $b$ denote one of the two square roots of $a$. Now we note that for $z\in \mathbb{C}\setminus \{0\}$ we have
$$\re z \gtrless 0 \iff \re \frac{1}{z} \gtrless 0,$$
and hence $$f(z) = \frac{1}{2}\left(z + \frac{a}{z}\right)$$ maps the half-planes $H^+ := \{ z : \re (z/b) > 0\}$ and $H^- := \{ z : \re (z/b) < 0\}$ to themselves, as well as the (extended) line $L := \{ z : \re (z/b) = 0\} \cup \{\infty\}$. So we can consider the families of iterates of $g = f\lvert_{H^+}$, the restriction of $f$ to $H^+$, and $h = f\lvert_{H^-}$, the restriction of $f$ to $H^-$. Let $\mathscr{G} = \bigl\{ g^n : n \in \mathbb{N}\setminus \{0\} \bigr\}$ and $\mathscr{H} = \bigl\{ h^n : n \in \mathbb{N}\setminus \{0\} \bigr\}$. Then $\mathscr{G}$ is a family of holomorphic functions with values in $H^+$, and $\mathscr{H}$ is a family of holomorphic functions with values in $H^-$.
Since half-planes are biholomorphically equivalent to the unit disk, the families $\mathscr{G}$ and $\mathscr{H}$ are normal. [Note: the two families are in fact locally bounded, so the little Montel theorem asserts the normality without the detour via the unit disk, but showing that directly is non-obvious.]
Next we note that $f$ has the two fixed points $b$ and $-b$ in $\mathbb{C}$ (in the Riemann sphere, $\infty$ is a third fixed point), and since
$$f'(z) = \frac{1}{2}\left(1 - \frac{a}{z^2}\right),$$
we have $f'(\pm b) = 0$, so $b$ and $-b$ are attractive fixed points. Hence the sequence of iterates $(f^n)$ converges uniformly to the constant $b$ on some neighbourhood of $b$, and it converges uniformly to the constant $-b$ on some neighbourhood of $-b$. By the identity theorem, it follows that every subsequence $(g^{n_k})$ of $(g^n)$ that converges locally uniformly on $H^+$ converges to the constant $b$, and ditto every subsequence $(h^{n_k})$ of $(h^n)$ that converges locally uniformly on $H^-$ converges to the constant $-b$. By the normality of $\mathscr{G}$ and $\mathscr{H}$, it then follows that $g^n \to b$ locally uniformly on $H^+$ and $h^n \to -b$ locally uniformly on $H^-$.
So we see that the sequence $(z_n)$ converges to $b$ when $z_0 \in H^+$, and it converges to $-b$ when $z_0 \in H^-$. For $z_0\in L$, the sequence either converges to $\infty$ - when $z_0$ is one of countably many points on $L$ such that $z_n = 0$ for some $n\in \mathbb{N}$ - or not at all.
Remark: For this specific example, the argument by normality is not the best, we can obtain a more precise view with simpler means, as illustrated below, or in user21820's answer. However, in geometrically more complicated situations, such an argument can be very powerful while still remaining (comparatively) simple.
II. The argument by the normal form:
As above, we let $b$ denote one of the square roots of $a$. The rational function $f$ from above has two special points, the attractive fixed points $b$ and $-b$. In such a situation, it is often illuminating to move the special points to $0$ and $\infty$ by conjugating the function with a Möbius transformation. The simplest Möbius transformation mapping $b\mapsto 0$ and $-b \mapsto \infty$ is $$S(z) = \frac{z-b}{z+b}.$$
The inverse is $$S^{-1}(w) = b\frac{1+w}{1-w},$$ and we compute
$$(f\circ S^{-1})(w) = \frac{1}{2}\left( b\frac{1+w}{1-w} + b\frac{1-w}{1+w}\right) = \frac{b}{2} \frac{(1+w)^2 + (1-w)^2}{1-w^2} = b\frac{1+w^2}{1-w^2}$$ and finally
$$(S\circ f \circ S^{-1})(w) = \frac{b\frac{1+w^2}{1-w^2}-b}{b\frac{1+w^2}{1-w^2}+b} = \frac{b(1+w^2)-b(1-w^2)}{b(1+w^2)+b(1-w^2)} = \frac{2bw^2}{2b} = w^2.$$
The behaviour of the iterates of $g = S\circ f \circ S^{-1}$ is now very easy to understand. It is immediately obvious that $g^n(w_0) \to 0$ for $\lvert w_0\rvert < 1$ and $g^n(w_0) \to \infty$ for $\lvert w_0\rvert > 1$, and even the behaviour on the unit circle is better understandable. We see that $g^n(w_0)$ stabilises at $1$ if and only if there are $m\in \mathbb{Z}$ and $k\in \mathbb{N}$ such that
$$w_0 = \exp \left(2\pi i \frac{m}{2^k}\right),$$
for all other $w_0$ on the unit circle the sequence is not convergent, it becomes periodic if $w_0 = \exp (2\pi i r)$ with a rational $r$, and does not repeat for other starting values on the unit circle.
Now we can use that $f^n = S^{-1}\circ g^n \circ S$ to obtain the explicit formula
$$f^n(z) = b\frac{1 +\left(\frac{z-b}{z+b}\right)^{2^n}}{1 - \left(\frac{z-b}{z+b}\right)^{2^n}},$$
which is essentially the same as in user21820's answer.