This work introduces novel unconditionally stable operator splitting methods for solving the time dependent nonlinear Poisson-Boltzmann (NPB) equation for the electrostatic analysis of solvated biomolecules. In a pseudo-transient continuation solution of the NPB equation, the nonlinear term is analytically integrated, so that the difficulties in direct treatment of the strong nonlinearity can be bypassed. However, in a pseudo-time NPB computation, the use of large time increments is necessary to reach the steady state efficiently. The existing alternating direction implicit (ADI) methods for the transient NPB equation are known to be conditionally stable, although being fully implicit. To overcome this difficulty, we propose several new operator splitting schemes, in both multiplicative and additive styles, including locally one-dimensional (LOD) schemes and additive operator splitting (AOS) schemes. The proposed schemes become much more stable than the ADI methods, and some of them are indeed unconditionally stable in dealing with solvated proteins with source singularities and non-smooth solutions. By using finite differences in space and implicit integrations in time, the numerical orders of the proposed schemes are found to be one in both space and time. Nevertheless, the precision in calculating the electrostatic free energy is low, unless a small time increment is used. Further accuracy improvements are thus considered, through constructing a Richardson extrapolation procedure and a tailored recovery scheme in treating the vacuum case. After acceleration, the optimized LOD method can produce a reliable energy estimate by integrating for a small and fixed number of time steps. Since one only needs to solve a tridiagonal matrix in each one dimensional subsystem, the overall computation is very efficient. The unconditionally stable LOD method scales linearly with respect to the number of atoms in the protein studies, and is over 20 times faster than the conditionally stable ADI methods.