Wormhole embedding diagrams

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.

• Just to spell out the failure, NDSolveValue[ {Derivative[1][z][u] == 2 Sqrt[-(( E^(-2 - 2 ProductLog[(u^2 - v0^2)/E]) u^2)/(1 + ProductLog[(u^2 - v0^2)/E])^2) + E^(-1 - ProductLog[(u^2 - v0^2)/E])/( 1 + ProductLog[(u^2 - v0^2)/E])], z[0] == 0}, z, {u, -1, 1} ] fails when v0 is greater than or equal to 1. Feb 19, 2018 at 16:39
• From the documentation: ProductLog[z] has a branch cut discontinuity in the complex z plane running -infinity from to -1/E. Mar 7, 2018 at 10:13