The "branch-cuts" along the imaginary axis seem to be a software artifact. The value of $\sqrt{z^2-1}$ is (apparently) computed by first computing $z^2-1$, and then taking the square root, choosing the square root with non-negative real part [and probably the one with non-negative imaginary part when the square root is purely imaginary]. That loses the information where the value whose square root shall be taken comes from, and thus leads to an even function (since $z^2-1$ is even), instead of a continuous branch, which would be an odd function (with the branch-cut $[-1,1]$). To avoid that, you need to remember where $z$ lies, and for $z$ in the left [or in the right, if you want the other branch] half-plane multiply the computed value with $-1$ [also for $z$ in the upper or lower half of the imaginary axis].
Then you should get a nice smooth plot, with a branch-cut on the interval $[-1,1]$.
For branch cuts I thought we had to identify where the argument of the square root is zero, which gives us the branch points and then we are free to take lines from those branch points to infinity in any way we like.
That describes the situation for $\sqrt{z}$ and similar, but generally, there are a few other considerations. If you have a rational function $f$, and consider $\sqrt{f(z)}$, the potential branch points are the zeros and poles of $f$.
But zeros and poles of even multiplicity are not actual branch points (in a neighbourhood of such a point you have $f(z) = (z-z_0)^{2k}\cdot h(z)$ with a holomorphic $h$ that does not vanish at $z_0$, so $\sqrt{h(z)}$ is holomorphic in a neighbourhood of $z_0$ and $\sqrt{f(z)} = (z-z_0)^k\sqrt{h(z)}$ is a well-defined meromorphic function in a neighbourhood of $z_0$), so only the zeros and poles of odd multiplicity are actual branch-points. [The situation is similar for $k$-th roots, zeros and poles whose multiplicity is a multiple of $k$ aren't branch-points, zeros and poles of other multiplicities are; essential singularities of $f$ may or may not be branch-points of $\sqrt[k]{f(z)}$, consider $\exp \frac{1}{z}$ and $z\exp \frac{1}{z}$.] Then you need to connect each branch-point to at least one other branch-point with a branch-cut to obtain a domain on which a holomorphic branch of $\sqrt[k]{f(z)}$ is defined.
In the situation here, $\infty$ is a pole of multiplicity $2$ of $z^2-1$, so we have only the two branch-points $1$ and $-1$ of $\sqrt{z^2-1}$. Any curve connecting the two branch-points gives a region on which a holomorphic branch of $\sqrt{z^2-1}$ is defined, whether that curve passes through $\infty$ or not. The most common choices of branch-cuts are the interval $[-1,1]$ and the other arc of the real axis (plus $\infty$) connecting the two points, $\bigl(\mathbb{R}\setminus (-1,1)\bigr) \cup \{\infty\}$ [the latter branch-cut makes $\sqrt{z^2-1}$ an even function on its domain].
The image shows the branch cuts in red for the first identity on the left and the second on the right.
I don't think splitting the radicand into factors is a good idea here. For $\sqrt{z+1}$ and $\sqrt{z-1}$ (resp. $\sqrt{1+z}$ and $\sqrt{1-z}$), you need a branch-cut connecting $\mp 1$ to $\infty$, but for $\sqrt{z^2-1}$, $\infty$ is in no way a special point. The branch-cuts of a product can [often] be chosen quite differently from the branch-cuts of the factors.
Why is there a branch cut extending the imaginary axis in the second case?
There shouldn't be. It's most likely an artifact of how the values are computed by the software.
How would I tailor the inbuilt python function to choose my own branch cuts?
Write a wrapper that multiplies the result of the inbuilt square root function with $-1$ for $z$ in the appropriate parts of the plane, according to your chosen branch-cuts.
Consider the function
$$
h(z) = \exp\left\{ \int_{i}^{z}\frac{w}{1+w^2}dw \right\}
= \exp\left\{\frac{1}{2}\int_{i}^{z}\frac{1}{w+i}+\frac{1}{w-i}\,dw\right\},
$$
where the integral from $i$ to $z$ is taken along any simple path from $i$ to $z$ in the open region $\Omega$ obtained by subtracting from $C$ the closed segment from $-i$ to $i$. It doesn't matter which path is chosen from $i$ to $z$ because the difference of exponent integrals along two different paths $\gamma_1$, $\gamma_2$ will be a closed path integral with the same winding number around $i$ as around $-i$, resulting in a difference of integrals that is an integer multiple of $2\pi i$. Hence $h(z)$ does not depend on the path chosen from $i$ to $z$, so long as it does not cross the segment from $-i$ to $i$. So $h$ is holomorphic in $\mathbb{C}\setminus [-i,i]$, and
$$
h'(z) = \exp\left\{\int_{i}^{z}\frac{w}{1+w^2}dw\right\}\frac{z}{1+z^2}=h(z)\frac{z}{1+z^2} \\
\implies \left(\frac{h(z)^2}{1+z^2}\right)'=0.
$$
After multiplying by an appropriate constant, the resulting function $g(z)=Ch(z)$ will satisfy $g(z)^2=1+z^2$ everywhere in $\mathbb{C}\setminus[-i,i]$.
So even thought there isn't a logarithm for $1+z^2$ in $\mathbb{C}\setminus[-i,i]$, everything works out because there is a "logarithm modulo integer multiples of $2\pi i$." I'll leave you to make sense of that phrase in whatever way you want.
There's no problem defining an actual holomorphic logarithm in $\mathbb{C}\setminus\{ (\infty,-i]\cup[i,\infty)\}$ because this is a simply connected region where $w/(1+w^2)$ is holomorphic, which gives an antiderivative of $w/(1+w^2)$ in this region, and hence a holomorphic square root function $\sqrt{1+w^2}$ in this region.
Best Answer
It may be useful to consider $w=\sqrt{1-z^2}+iz$. Isolating the square root and squaring, you may note that the $z^2$ terms cancel, and you end up with $$z=\frac1{2i}\Bigl(w-\frac1w\Bigr),$$ and then subsituting this into the equation $\sqrt{1-z^2}=w-iz$ you also find $$\sqrt{1-z^2}=\frac12\Bigl(w+\frac1w\Bigr).$$ Thus both $z$ and $\sqrt{1-z^2}$ are single-valued functions of $w$, so it should be easier to analyze the given function as $\log w$.
There are two values of $w$ for each value of $z$: Replacing $w$ by $-1/w$ leaves $z$ unchanged, and flips the sign of $\sqrt{1-z^2}$.
Note that imaginary $w$ gives real $z$: Your proposed branch cuts in the $z$ plane go along the imaginary axis in the $w$ plane, from $\pm i$ to infinity in opposite directions, but also (if you pick the other branch of the square root) from $\pm i$ to $0$. Thus your branch cuts divides the $w$ plane in two halves along the imaginary axis, and you end up with not having to pick further branch cuts for the logarithm.
(In my first edition of this answer I got a little confused because I was thinking of getting the full Riemann surface for the given function. I hope I managed to fix this before confusing anybody else too much. The analysis I give here is perhaps better for understanding the Riemann surface; it could well be overkill for the branch cut question. Oh well …)