According to this answer, $$z_\infty:=\lim_n z_n=I:=\int_0^\infty F(s)G(s)\,ds,$$ where $$F(s):=\prod_{k=1}^\infty\frac{1}{\sqrt{1+2s/k^3}},\quad G(s):=\sum_{k=1}^\infty\frac k{k^3+2s}.$$
Mathematica can express $F$ and $G$ in terms of functions built-in in Mathematica (and these expressions should be rather straightforward to verify), and then the Mathematica command NIntegrate numerically evaluates $z_\infty=I$ as $\approx1.99218$ -- close to $2$, but not $2$; see the image of the corresponding Mathematica notebook below.
Using the facts that (i) $F$ and $G$ are positive, decreasing, and convex, and hence $FG$ is so, and that (ii) Mathematica can find the values of all its built-in functions with any degree of accuracy, it is rather straightforward to show that, in fact, $I<2$.
