This problem pertains to viscous as opposed to inviscid flow and the Bernoulli equation is not relevant.
It is also assumed that the flow is steady and unidirectional with vanishing radial and axial velocity components, $u_r = u_z = 0$, and the azimuthal velocity component $u_\varphi(\sigma)$ is a function only of the radial coordinate $\sigma$.
In this case, the $\varphi$- component of the Navier-Stokes equations reduces to
$$\frac{1}{\sigma} \frac{\partial}{\partial \sigma}\left(\sigma \frac{\partial u_\varphi}{\partial \sigma} \right) - \frac{u_\varphi}{\sigma^2} = 0,$$
with general solution $a\sigma + b \sigma^{-1}$. The constants $a$ and $b$ are found by applying the boundary conditions $u_\varphi(R) = 0$ and $u_\varphi(\kappa R) = \Omega_i\kappa R$ yielding
$$ a = - \frac{\Omega_i \kappa^2}{1 - \kappa^2}, \,\,\, b = \frac{\Omega_i \kappa^2R^2}{1 - \kappa^2}, \,\,\,u_\varphi = - \frac{\Omega_i \kappa^2}{1 - \kappa^2}\left(\sigma - \frac{R^2}{\sigma} \right)$$
(You obtained the correct functional form for $u_\phi$ but the coefficient is incorrect -- which you will see by checking the boundary conditions.)
Introducing the dimensionless variable $\zeta = \sigma/R$ we obtain
$$u_\varphi = - \frac{ \kappa^2R\Omega_i}{1 - \kappa^2}\left(\zeta - \frac{1}{\zeta} \right)$$
Free Surface
By symmetry the pressure field $p(\zeta,z)$ is independent of $\varphi$ and depends only of the radial coordinate $\zeta$ and the axial coordinate $z$.
The radial component of the Navier-Stokes equations reduces to a balance between the radial pressure gradient and the centripetal acceleration, given with $\zeta = \sigma/R$ by
$$\frac{\partial p}{\partial \zeta} = \frac{\rho u_\varphi^2}{\zeta}$$
Substituting for $u_\phi$ we get,
$$\frac{\partial p}{\partial \zeta} = \rho\frac{ (\kappa^2R\Omega_i)^2}{(1 - \kappa^2)^2}\left(\zeta - \frac{2}{\zeta} + \frac{1}{\zeta^3} \right),$$
and after integrating both sides,
$$p(\zeta,z) = \frac{\rho}{2}\frac{ (\kappa^2R\Omega_i)^2}{(1 - \kappa^2)^2}\left(\zeta^2 - 4 \log \zeta - \frac{1}{\zeta^2} \right) + C(z)$$
Here the integration constant $C(z)$ can depend on $z$.
The axial component of the Navier-Stokes equations reduces to $\frac{\partial p}{\partial z} = \rho g$ which implies that $C(z) = \rho g z + C_0$, and
$$p(\zeta,z) = \frac{\rho}{2}\frac{ (\kappa^2R\Omega_i)^2}{(1 - \kappa^2)^2}\left(\zeta^2 - 4 \log \zeta - \frac{1}{\zeta^2} \right) + \rho g z + C_0$$
At the outer cylinder we have $\zeta = 1$ and free surface height $Z_R$. Consequently,
$$p(1,z_R) = \rho gz_R + C_0 = p_A ,$$
where $p_A$ is atmospheric pressure. At a radial position $\zeta$ where the free surface height is $z_\zeta$ we also have atmospheric pressure. Hence,
$$p_A = p(\zeta,z_\zeta) = \frac{\rho}{2}\frac{ (\kappa^2R\Omega_i)^2}{(1 - \kappa^2)^2}\left(\zeta^2 - 4 \log \zeta - \frac{1}{\zeta^2} \right) + \rho g z + \underbrace{C_0}_{= p_A - \rho g z_R}$$
Eliminating $p_A$ from both sides and rearranging we get
$$z_R - z_\zeta = \frac{1}{2g}\frac{ (\kappa^2R\Omega_i)^2}{(1 - \kappa^2)^2}\left(\zeta^2 - 4 \log \zeta - \frac{1}{\zeta^2} \right) $$
