You only need to choose a branch cut so that all of the branch points lie on the cut. Remember that a branch point is a point where the function is discontinuous when traversing an arbitrarily small circuit about the point. By choosing a branch cut like this you're essentially acknowledging this behavior and choosing specifically where this discontinuity will occur.
The points $z = \pm 1$ are clearly branch points, and they are the only finite branch points. The point at infinity may also be a branch point. In this case it isn't (shown below), but if it were then you would need to choose your cut so that it contains all three branch points ($1$, $-1$, and $\infty$). You could do this by choosing it to be the real interval $(-\infty,1]$, but there are infinitely many other choices you could make too.
To see that the point at infinity is not a branch point, rewrite the function as
$$
f(z) = z\,\left(\frac{1}{z}-1\right)^{3/5}\left(\frac{1}{z}+1\right)^{2/5}.
$$
Let $|z| > 1$. By doing this, we can ensure that neither of the two quantities $a = \frac{1}{z}-1$ or $b = \frac{1}{z}+1$ will make a circuit about the origin, regardless of what $z$ is doing. They are both restricted to unit disks centered at $-1$ and $1$, respectively. Indeed, $|a+1|<1$ and $|b-1|<1$.
Now, with this condition on $|z|$, if $z$ makes a circuit about the origin (equivalently, a circuit about infinity), the change in the argument of $f$ is precisely $2\pi$, and the function remains unchanged. We have thus continuously traversed an arbitrarily small circuit about infinity, and hence infinity is not a branch point of $f$.
As a result, we only need to choose a branch cut which includes $1$ and $-1$. The real inverval $[-1,1]$ is the simplest choice.
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
In case this link Defining a holomorphic branch of $\log(z^2-1)$ does not address your specific question:
If you simply cut along the negative real axis, the resulting domain will still include the point $z=+1$ at which no branch of $\log(z^2-1)$ can be defined. Any cut must eliminate bot $z=1$ and $z=-1$.
But there is another way to approach obtaining a single valued holomorphic version of $\log f(z)$, where for simplicity we will take $f(z)$ to be holomorphic except at a finite number of zeros and / or poles. We can require a domain $\Omega$ in which every path $\gamma \subset \Omega$ from $z_0$ to $z$ gives the same value, $$ \log z = \log(z_0) + \int_\gamma \frac{f'(z)}{f(z)} ~ dz. $$ This can only arise if all paths joining $z_0$ to $z$ are homotopic in $\Omega$, that is each can de deformed continuously into the other without encountering poles of $f'(z)/f(z)$.
That will be the case if, as above, the domain excludes all the zeros of $f(z)$ and the paths $\gamma$ cannot circulate the poles of $f'(z)/f(z)$. There may be exceptional cases where the sum of their residues is zero, and then there may be more flexibility, but we will momentarily exclude this case. Then, to avoid circulation, a sufficient condition is that the domain be simply connected.
In our case $f(z) = z^2-1$ and the residue sum for $f'(z)/f(z)$ is non-zero, so our cut must prevent a circuit to include $z =\pm 1$. That still leaves many choices. One is the domain that is cut along the real axis from $z=-\infty$ to $z=1$, so that $\Omega = \mathbb C \setminus (-\infty, 1]$. Another, used in the link, is to take $\Omega = \mathbb C \setminus \Big( (-\infty, -] \cup [1,\infty) \Big).$ Both give different holomorphic versions of the $\log$ function.
Any cut can be made to work if it meets these criteria: choose the domain $\Omega$, take a point $z_0 \in \Omega$ and its corresponding value $\log z_0$ and define, $$ \log(z) = \log(z_0) + \int_{z_0}^z \frac{2w}{w^2-1}~dw = \log(z_0) + \int_{z_0}^z \frac{1}{w+1} + \frac{1}{w-1}~dw. $$ This will give a single valued holomorphic branch. The value is independent of the path by Cauchy's theorem.
Coming back to an example where the residues do cancel, consider $\displaystyle \log\frac{z-1}{z+1}$.
Here $f'(z)/f(z) = 2/(z^2-1)$, the residues are $\pm1$ and the corresponding cut for $\log z$ need only exclude jointly the points $z = \pm 1$. A commonly used domain in this case is $\Omega = \mathbb C \setminus [-1,1]$, although you could again use $\Omega = \mathbb C \setminus \Big( (-\infty,-1] \cup [1,\infty) \Big)$.
Addition: I should have mentioned that for $f'(z)/f(z)$ every zero of $f$ has residue 1 and every pole residue -1, so the sum of the residues is simply the number of zeros less the number of poles.