First of all: this is question lies in the context of Surfaces and Embbedings on differential geometry. More precisely in the context of Kruskal coordinates and how to plot a 3D dynamical Schwarzschild Wormhole $[1]$, $[2]$, $[3]$. The main context is given on reference $[1]$.
PART I) Some Physics and Mathematics
I.1) A "nothing" on Gravity and Black Holes
Black Holes are solutions of Einstein Field Equations. More preciselly, we have then the master (tensor) equation:
$$\textbf{Ric} - \frac{1}{2}R\textbf{g} = \frac{8\pi G}{c^{4}} \textbf{T} \tag{1}$$
The solutions of equation $(1)$ are given by the metric tensors:
$$\textbf{g} = g_{\mu\nu} \textbf{dx}^{\mu} \otimes \textbf{dx}^{\nu} \tag{2}$$
In a coordinate (a.k.a. where we do calculations) chart, we can express the scalar quantity called line element which are preciselly the scalar field when the metric tensor acts on tangent vector fields on a manifold $\mathcal{M}$:
$$\textbf{g}(d\textbf{r},d\textbf{r}) = g_{\mu\nu}dx^{\mu}dx^{\nu} := ds^2\tag{3}$$
Therefore, given a metric tensor on a manifold we caracterize the geometry, in a point-wise fashion, of the Manifold.
An important fact on tensors is something called general covariance: the physics is independent of coordinates; the nature phenomena do not care about coordinates; the equations of physics must to be independent of coordinates; the nature fenomena can be well described by using cartesian coordinate system, and equally well described by any curvilinear coordinate system. The mathematical object that captures this fact are tensors.
Therefore the the abstract object (metric tensor) is independent of coordinate transformations:
$$ \textbf{g'} = \textbf{g} \tag{4}$$
But in a level of coordinate charts, their coordinates do change in a form such as:
$$ g'_{\mu\nu} = \frac{\partial x^{\alpha}}{\partial x'^{\mu}}\frac{\partial x^{\beta}}{\partial x^{\nu'}}g_{\alpha \beta}\tag{5}$$
A Schwarzschild Black Hole have the following metric tensor:
$$ ds^{2} = -\Bigg(1- \frac{2GM}{c^{2}r}\Bigg)c^{2}dt^{2} + \frac{1}{1-\frac{2GM}{c^{2}r}}dr^{2} - r^{2}d\theta^{2} + r^{2}sin^{2}(\theta)d\phi^{2}\tag{6}$$
We can readly see the coordinate chart here: the spherical coordinates extended to spacetime physics;
$$\mathrm{Coordinates} \implies (t,r,\theta,\phi) \tag{7}$$
Therefore we are dealling with a four-dimensional spacetime $\mathcal{M}$.
I.2) Embbebing Procedure of Schwarschild Wormhole
The basic context here is: we need to represent that Four-dimensional spacetime $(6)$ in a Three dimensional Euclidean Space. Then, we must to embbed a surface in a 3D space.
We then "supress" two dimensions, i.e, the fix a moment of coordinate time $t = T \in \mathbb{R}$ and a coordinate angle $\theta = \frac{\pi}{2}$ (in fact this angle gives us the geometry of the "equatiorial plane"). Thefore the differentials are zero: $dt^{2} = dT^{2} = 0 $, $d\theta^{2} = d\frac{\pi}{2}^{2} = 0 $ and then, the metric becomes:
$$ ds^{2} = \frac{1}{1-\frac{2GM}{c^{2}r}}dr^{2} + r^{2}d\phi^{2}\tag{8}$$
The next step is then to "compare" this $2D$ metric tensor with the Cilindrical metric tensor:
$$ ds^{2} = dz^{2} + dr^{2} + r^{2}d\phi^{2} \equiv \Bigg[ 1+\Bigg(\frac{dz}{dr}\Bigg)^{2} \Bigg]dr^{2} + r^{2}d\phi^{2}\tag{9}$$
Then we say that:
$$\Bigg[ 1+\Bigg(\frac{dz}{dr}\Bigg)^{2} \Bigg]dr^{2} = \frac{1}{1-\frac{2GM}{c^{2}r}}dr^{2} \tag{10}$$
Which implies, finally:
$$ \frac{dz}{dr} = \frac{1}{\sqrt{1-\frac{2GM}{c^{2}r}}} - 1 \tag{11}$$
We can integrate this function and discover the very function of interrest here:
$$ z(r) = \pm\int^{r}_{r_{0}} \frac{1}{\sqrt{1-\frac{2GM}{c^{2}r}}} - 1 dr \tag{11}$$
Integrating equation $(11)$ and using the function $\mathrm{RevolutionPlot3D}$ we plot the wormhole-like surface
PART II) My Doubt and What kind of a Mathematica Program I Want to Write
II.1) The Kruskal Numerical Embbeding Integral
There are a lot to cover, but the physical fact is: this wormhole throat isn't someting "constant in time" I mean, in Schwarzshild black holes, the wormhole closes very after the black hole is formed, then using other coordinates, the Kruskal-Szekeres coordinates, we can represent the "dynamical wormhole" in the "passage of time", i.e., fo each value of "kruskal time" we have a different surface!
The coordinates change then for the Kruskal-Szekeres coordinates:
$$\mathrm{Coordinates'} \implies (v,u,\theta,\phi) \tag{12}$$
Where $v$ are called the "Kruskal Time". It's important to mention that this coodinate isn't the time of our clocks, in fact event $t$ is not the time of our clocks, but rather they carries a temporal signature.
In this coordinate chart the Embbeding integral (such as $(11)$ ) becomes:
$$ z(u) = \pm\int_{0}^{u} \frac{dz}{du}du = \pm\int_{r(0)}^{r(u)} \frac{dz}{du}\frac{du}{dr}dr = \pm\int_{r(0)}^{r(u)} \sqrt{g_{rr} - 1} dr =$$
$$ = \pm\int_{r(0)}^{r(u)} \sqrt{\frac{2\Big( e^{\frac{r}{2}} -v_{0}^2 \Big)}{re^{\frac{r}{2}} - 2\Big( e^{\frac{r}{2}} -v_{0}^2 \Big)}} dr \tag{12}$$
So, $v_{0}$ are the constant Kruskal times, therefome for each $v_{0}$ we can have a integral on $u$ ! Also the functions $r(u)$ are in fact:
$$ r(u,v) =: 2m \Biggr[1+\mathrm{LambertW}\Bigg(\frac{u^2-v^2}{e}\Bigg) \Biggr] \tag{13}$$
The usage, in this problem, of the function $(13)$ is more preciselly (for $m=1$ ),
$$ r(u) \equiv r(u,v_{0}) =: 2 \Biggr[1+\mathrm{LambertW}\Bigg(\frac{u^2-v_{0}^2}{e}\Bigg) \Biggr] \tag{14}$$
Therefore, the integral $(12)$ becomes finally :
$$ z(u) = \pm\int_{r(0) = 2 \Biggr[1+\mathrm{LambertW}\Bigg(\frac{-v_{0}^2}{e}\Bigg) \Biggr] }^{r(u)=2 \Biggr[1+\mathrm{LambertW}\Bigg(\frac{u^2-v_{0}^2}{e}\Bigg) \Biggr]} \sqrt{\frac{2\Big( e^{\frac{r}{2}} -v_{0}^2 \Big)}{re^{\frac{r}{2}} - 2\Big( e^{\frac{r}{2}} -v_{0}^2 \Big)}} dr \tag{15}$$
So we NEED to calculate the integral $(15)$, and this is a numerical task. In short we calculate (in some way) the integral $(15)$, and then we plot the very 3D surface:
$$ \textbf{X} = [r(u)cos(\phi), r(u)sin(\phi), z(u)] \tag{16}$$
II.2) My NIntegrate Calculation's ListPlot
So I've done a numerical Integration using NIntegrate but I simply can't use it as a $z[u]$ function! Take a look:
Manipulate[
ListPlot[Table[
NIntegrate[
Sqrt[((2*(Exp[(r)/(2)] - ((v0)^2)) ))/(r*
Exp[(r)/(2)] - ((2*(Exp[(r)/(2)] - ((v0)^2)) )))], {r, \
(2*(1)*(1 + (
LambertW[((( 0^2 ) - (v0^2))/(Exp[
1]))] ))), (2*(1)*(1 + (
LambertW[((( u^2 ) - (v0^2))/(Exp[1]))] )))}], {u, 0, k,
1}], Joined -> True, InterpolationOrder -> 1, Mesh -> All], {k, 1,
25}, {v0, -0.99999, 0.99999}]
This code in fact calculates the integral numerically, but I can't use it in the $\mathrm{ParametricPlot3D}$.
PART III) Some Already Done Embbedings and Suggestions
III.1) My older code based on reference $[4]$
In fact an user $[4]$ have "solved" the problem. I adapted his code a little bit:
r[u_, v0_] := 1 + LambertW[(u^2 - v0^2)/E];
z[u_, v0_] :=
2 Sqrt[(1/
r[u, v0] Exp[-r[u, v0]] - (u Exp[-r[u, v0]]/(r[u, v0]))^2)];
Z = ParametricNDSolveValue[{z'[u] == z[u, v0], z[0] == 0},
z, {u, -1, 1}, {v0}]
Manipulate[
ParametricPlot3D[{r[u, v0] Cos[phi], r[u, v0] Sin[phi],
Z[v0][u]}, {u, -1, 1}, {phi, 0, 2 Pi}, PlotTheme -> "Classic",
PlotRange -> {{-1.5, 1.5}, {-1.5, 1.5}, {-10, 10}},
BoxRatios -> {2, 2, 2}], {v0, -0.99999999999999, 0.99999999999999}]
But he is solving a differential equation, I want to use the integral calculation!
III.2) Hamilton's Dynamical Embbeding GIF
This guy Hamilton have acheived the supreme plot that I want to do, in fact I want to reproduce his dynamical plot, i.e., this question of mine and all my doubts exists just to reproduce his plot. You can take a look on the dynamical wormhole throat in the following:
III.3) My Professor's adivice
My professor told me a interresting possible solution:
Maybe you can try to solve the numerical integral for various values of $u$, and get several points of the type $(u, z)$, and then you can for example interpolate these points to get a smooth function of $z = z (u)$, and use this function in the $\mathrm{ParametricPlot3D}$. Try and see if it works.
Seems to be right, but I simply do not know how to translate get several points of the type $(u, z)$, and then you can for example interpolate these points to get a smooth function of $z = z (u)$ into Mathematica code.
PART IV) My Doubt
So, my numerical integration on section II.2 seems to be correct, therefore it's just a matter of fact on how can I rotate this curve. But the way I wrote the whole thing I cannot use the $\mathrm{RevolutionPlot3D}$, so my doubt is:
How Can I use $\mathrm{RevolutionPlot3D}$ in my code II.2?
$$ * * * $$
$[1]$ Embeddings and time evolution of the Schwarzschild wormhole - https://arxiv.org/abs/1107.4871
$[2]$ Gravitation - Misner, Thorne, Wheeler Pages 612 to 615
$[3]$ Wormholes in spacetime and their use for interstellar travel: A tool for teaching general relativity - https://aapt.scitation.org/doi/10.1119/1.15620
