Wormhole embedding diagrams

Question

I am trying to reproduce the embedding diagrams for the evolution of a Schwarzschild wormhole described in this paper.

Following the paper notation, we denote the Kruskal coordinates by $(v,u)$. For a constant Kruskal time $v_0$, the embedding diagram is given parametrically (see Eqs. (12) and (14) in the paper) by

$$ \vec{x}(u,\phi,v_0)=\left(r(u,v_0)\cos\phi,r(u,v_0)\sin\phi,z(u,v_0)\right) $$

where $r(u,v_0)=1+W\left(\frac{u^2-v_0^2}{e}\right)$, $W$ is the Lambert function. The function $z(u,v_0)$ is given by the integral $$ z(u,v_0) = 2\int_0^u\mathrm{d}u \left(\frac{1}{r}e^{-r}-\left(\frac{ue^{-r}}{r}\right)^2\right)^{1/2}$$

With $r=r(u,v_0)$ (see Eqs. (22) and (24)).

I have tried to implement this in Mathematica:

r[u_, v0_] := 1 + ProductLog[(u^2 - v0^2)/E];
zprime[u_, v0_] := 2 Sqrt[(1/r[u, v0] Exp[-r[u, v0]] - (u Exp[-r[u, v0]]/(r[u, v0]))^2)];

zsolution = ParametricNDSolveValue[{z'[u] == zprime[u, v0], z[0] == 0}, 
  z, {u, -1, 1}, {v0}]

wormhole[v0_?NumericQ] := ParametricPlot3D[{r[u, v0] Cos[phi], r[u, v0] Sin[phi], zsolution[v0][u]},
    {u, -1, 1}, {phi, 0, 2 Pi}, 
  SphericalRegion -> True, BoxRatios -> 1, FaceGrids -> None, 
  Mesh -> 10, PlotStyle -> {Opacity[0.8]}, 
  PlotTheme -> {"Classic", "ClassicLights"}, Boxed -> False, 
  Axes -> False, ImageSize -> Automatic, ViewPoint -> Front]

which reproduce the time evolution for $|v_0|<1$, that is, the formation and the closure of the wormhole.

Table[wormhole[v0], {v0, {-0.9999, -0.9, -0.7, 0, 0.7, 0.9, 0.9999}}]

but neither the instant of pinching off $v_0=1$ nor the separation of the two exterior spaces $v_0>1$. The problem is that Mathematica finds a singularity, but I don't know how to solve this problem.

QUESTION: How can I reproduce the embedding diagrams for the times $|v_0|\geq 1$, as done in the paper?

Thanks in advance!

Note: I am using units of $r/2M$, not $r/M$ as used in the paper.