Finite difference implicit schema for wave equation 1D not unconditionally stable? The wave equation 1D with constant density is defined as:
\begin{equation}
\frac{\partial^2 U}{\partial t ^2} = V^2 \frac{\partial^2 U}{\partial x ^2}
\label{eqa}
\end{equation}
And the implicit difference schema:
\begin{equation}
U_j^{n+1}  =  \left( \frac{\Delta t  V_j^n}{\Delta s} \right) ^2 \left(  U_{j+1}^{n+1} - 2 U_j^{n+1} + U_{j-1}^{n+1}  \right) + 2 U_j^n - U_j^{n-1}
\label{eq:b}
\end{equation}
Where $\Delta x = \Delta s $.
Von Newman stability analysis.
Since $ u(x,t) = e^{i(wt+kx)} $ is a solution for the discrete schema above can be written as:
\begin{eqnarray}
u(x_j,t_n) &=& e^{i(wt_n+kx_j)} \nonumber \\
&=& e^{iwn\Delta t} e^{ikj\Delta s} \nonumber \\
&=& \epsilon^n e^{ikj\Delta s}
\label{eq:c}
\end{eqnarray}
Where $\epsilon = e^{iw \Delta t} $ and $ i = \sqrt{-1} $ . 
To maintain stability we should make sure that $ \epsilon \leq 1 $ not growing exponentially with increasing time steps. So applying the previous equation in the implicit differences schema we can analyse the growth.
\begin{multline}
 U_j^{n+1}  =  \left( \frac{\Delta t  V_j^n}{\Delta s} \right) ^2 \left(  U_{j+1}^{n+1} - 2 U_j^{n+1} + U_{j-1}^{n+1}  \right) + 2 U_j^n - U_j^{n-1} \\
\epsilon^{n+1} e^{ikj\Delta s}  =  r^2 \left(  \epsilon^{n+1} e^{ik(j+1)\Delta s}  - 2 \epsilon^{n+1} e^{ikj\Delta s} + \epsilon^{n+1} e^{ik(j-1)\Delta s}  \right) \\ 
+ 2 \epsilon^n e^{ikj\Delta s}  - \epsilon^{n-1} e^{ikj\Delta s} \\
 1 = r^2 (e^{ik\Delta s} -2 + e^{-ik\Delta s}) + 2 \epsilon^{-1} - \epsilon^{-2} \\
\epsilon^{-2} - 2 \epsilon^{-1} -  r^2 (e^{ik\Delta s} -2 + e^{-ik\Delta s}) + 1 = 0  \\
\epsilon^{-2} - 2 \epsilon^{-1} - 4 r^2 \sin^2 \left( k \Delta s / 2 \right) + 1 = 0   
 \label{eq:d}
\end{multline}
Where $ r = \left( \frac{\Delta t  V_j^n}{\Delta s} \right) $
At the previous equation we can substitute $ \phi = \epsilon^{-1} $ turning it to a second degree bellow:
\begin{eqnarray}
\phi^2 - 2 \phi -4 r^2 \sin^2 \left(k\Delta s/2 \right) +1 = 0  \nonumber \\
\phi^2 - 2 \phi + c = 0  \nonumber
\end{eqnarray}
With:
    $$ c = 1 -4 r^2 \sin^2 \left(k\Delta s/2 \right) $$
Has roots $ \phi^{'} $ and $ \phi^{''} $ as bellow. Replacing $c$ also.
\begin{eqnarray}
 &=& \frac{2 \pm \sqrt{4 - 4c}}{2} = 1 \pm \sqrt{1-c} \nonumber \\
 &=& 1 \pm \sqrt{4 r^2 \sin^2 \left(k\Delta s/2 \right) } \nonumber \\
 &=& 1 \pm 2 r \sin \left(k\Delta s/2 \right)  \nonumber \\
 \phi^{'} &=& 1 + 2 r \sin \left(k\Delta s/2 \right)   \\
 \phi^{''} &=& 1 - 2 r \sin \left(k\Delta s/2 \right)   
\end{eqnarray}
Going back to $\epsilon$ we get :
\begin{eqnarray}
 \epsilon^{'} &=& \frac{1}{1 + 2 r \sin \left(k\Delta s/2 \right)}   \\
\epsilon^{''} &=& \frac{1}{1 - 2 r \sin \left(k\Delta s/2 \right) }   
\end{eqnarray}
But with $r \sin \left(k\Delta s/2 \right) \to 0$ we get that $\epsilon^{''} \to \infty $
Should it not be unconditionally stable? What is wrong in here?



Update: After awesome answer by Jitse, fixing the signal we get:
\begin{eqnarray}
 \phi^{'} &=& 1 + i 2 r \sin \left(k\Delta s/2 \right)   \\
 \phi^{''} &=& 1 - i 2 r \sin \left(k\Delta s/2 \right)   
\end{eqnarray}
Analyzing the modulus, since $  \phi^{'} $ and $  \phi^{''} $ are conjugate pairs of the same complex number the modulus is the same for both.
We get:
$$ \| \phi^{'} \| = \| \phi^{''} \| = \sqrt{1+ 4 r^2 \sin ^2 \left(k\Delta s/2 \right)}$$
So since $ \| \phi \| \geq 1 \ \ \forall \ r, \ k,\ \Delta s $ than 
$$ \| \epsilon \| = \frac{1}{ \sqrt{1 + 4 r ^2 \sin ^2 \left(k\Delta s/2 \right)}} $$
$$ \| \epsilon \| \leq 1 \ \ \forall \ r, \ k,\ \Delta s$$
Awesome, unconditionally stable!!

Your method for Neumann stability analysis is correct, and indeed I would expect the scheme to be unconditionally stable.
I do not understand the step where you introduce the sine. I think that
$$ e^{ix} - 2 + e^{-ix} = -4\sin^2(\tfrac12x). $$
The minus sign seems to have been lost in your computation. 
