I tried to calculate the time-like geodesics of the Morris-Thorne Wormhole $ [1]$ , $ [2]$ , $ [3]$ for redshift function $ \Phi(r) = 0$ and $ b(r) = \sqrt{r_{0} r}$ . But I don’t know for sure if all the plots that I made are correct.

So, the metric is:

$ $ ds^2 = -e^{2\Phi(r)}dt^2+\Biggr\{ 1-\frac{b(r)}{r}\Biggr\}^{-1}dr^2+r^2[d\theta^2+sin^2(\theta)d\phi^2]\tag{1}$ $

And the plots are:

The various spheres are the 2-spheres due to spherical symmetry of the metric, and the central sphere is the sphere of the radius $ r_{0}$ , therefore is the throat.

My doubt is:

Am I seeing a particle orbiting the throat of a Morris-Thorne Wormhole? In other words is it correct to say that the particle approches from infinity, orbits the throat and returns to the same universe?

$ $ * * * $ $

$ [1]$ LOBO.F.S.N; *Exotic solutions in General Relativity*,arxiv:0710.4474v1, 2007

$ [2]$ LOBO.F.S.N; *Wormholes, Warp Drives and Energy Conditions*,Springer, Vol.189, 2017

$ [3]$ VISSER.M; *Lorentzian Wormholes*, Springer, AIP Press, 1995

**Appendix: ***Mathematica* code:

`n = Length[coords]; a = 0; r0; \[CapitalPhi][r] = 0; B[r] = (Sqrt[r*r0]); tt = -(Exp[ 2*\[CapitalPhi][ r] ]); rr = (1)/(1 - ((B[r])/(r))); \[Theta]\[Theta] = r^2; \[CurlyPhi]\[CurlyPhi] = r^2*(Sin[\[Theta]]*Sin[\[Theta]]); t\[CurlyPhi] = 0; metric = {{tt, 0, 0, 0}, {0, rr, 0, 0}, {0, 0, \[Theta]\[Theta], 0}, {0, 0, 0, \[CurlyPhi]\[CurlyPhi]}}; metric // MatrixForm inversemetric = Simplify[Inverse[metric]]; inversemetric // MatrixForm christoffel := christoffel = Simplify[Table[ 1/2*Sum[inversemetric[[i, s]]*(D[metric[[s, j]], coords[[k]]] + D[metric[[s, k]], coords[[j]]] \[Minus] D[metric[[j, k]], coords[[s]]]), {s, 1, n}], {i, 1, n}, {j, 1, n}, {k, 1, n}]] listchristoffel := Table[If[UnsameQ[christoffel[[i, j, k]], 0], {ToString[\[CapitalGamma][i, j, k]], christoffel[[i, j, k]]}], {i, 1, n}, {j, 1, n}, {k, 1, j}] TableForm[ Partition[DeleteCases[Flatten[listchristoffel], Null], 2], TableSpacing -> {2, 2}] geodesic := geodesic = Simplify[Table[\[Minus]Sum[ christoffel[[i, j, k]] coords[[j]]' coords[[k]]', {j, 1, n}, {k, 1, n}], {i, 1, n}]] listgeodesic := Table[{"d/d\[Tau]" ToString[coords[[i]]'], " =", geodesic[[i]]}, {i, 1, n}] TableForm[listgeodesic, TableSpacing -> {2}] max\[Tau] = 750; ivs = {0, 0, 0.088}; ics = {0, 6.5, \[Pi]/2, 0}; r0 = 1; computeSoln[max\[Tau]i_, ivsi_, icsi_] := Block[{ivs, ics, i, \[Chi], tmp, soln}, ics = icsi; ivs = Join[{\[Chi]}, ivsi]; op1 = Table[coords[[i]] -> ics[[i]], {i, 0, n}]; tm = metric /. op1; tmp = ivs.(tm.ivs); \[Chi]slv = Solve[tmp == uinvar, \[Chi]]; ivs[[1]] = Last[\[Chi] /. \[Chi]slv]; op = {Derivative[1][t] -> Derivative[1][t][\[Tau]], Derivative[1][r] -> Derivative[1][r][\[Tau]], Derivative[1][\[Theta]] -> Derivative[1][\[Theta]][\[Tau]], Derivative[1][\[CurlyPhi]] -> Derivative[1][\[CurlyPhi]][\[Tau]], t -> t[\[Tau]], r -> r[\[Tau]], \[Theta] -> \[Theta][\[Tau]], \[CurlyPhi] -> \ \[CurlyPhi][\[Tau]]}; deq = Table[ coords[[i]]''[\[Tau]] == Simplify[geodesic[[i]] /. op], {i, 1, n}]; deq = Join[deq, Table[coords[[i]]'[0] == ivs[[i]], {i, 1, n}], Table[coords[[i]][0] == ics[[i]], {i, 1, n}]]; soln = NDSolve[deq, coords, {\[Tau], 0, max\[Tau]i}]; soln] uinvar = \[Minus]1; sphslnToCartsln[soln_] := Block[{xs, ys, zs}, xs = r[\[Tau]] Sin[\[Theta][\[Tau]]] Cos[\[CurlyPhi][\[Tau]]] /. soln; ys = r[\[Tau]] Sin[\[Theta][\[Tau]]] Sin[\[CurlyPhi][\[Tau]]] /. soln; zs = r[\[Tau]] Cos[\[Theta][\[Tau]]] /. soln; {xs, ys, zs}] udotu[solni_, \[Tau]val_] := Block[{x\[Alpha], u\[Alpha]}, x\[Alpha] = Table[coords[[i]][\[Tau]] /. solni, {i, 1, n}] // Flatten; u\[Alpha] = D[x\[Alpha], \[Tau]]; x\[Alpha] = x\[Alpha] /. \[Tau] -> \[Tau]val; u\[Alpha] = u\[Alpha] /. \[Tau] -> \[Tau]val; u\[Alpha].((metric /. Table[coords[[i]] -> x\[Alpha][[i]], {i, 1, n}]).u\[Alpha])] coordlist[\[Tau]in_] := Table[ToString[coords[[i]]] <> " = " <> {ToString[coords[[i]][\[Tau]in] /. soln // First]}, {i, 1, n}] soln = computeSoln[max\[Tau], ivs, ics]; xyzsoln = sphslnToCartsln[soln]; Join[{"Final Coordinates:"}, coordlist[max\[Tau]]] // TableForm Join[{{"", "", "", "u.u values"}}, Table[{"\[Tau]=", ToString[i], "->", udotu[soln, i]}, {i, 0, max\[Tau], max\[Tau]/5}]] // TableForm {Plot[Evaluate[ Table[coords[[i]][\[Tau]] /. soln, {i, 1, n}]], {\[Tau], 0, max\[Tau]}, AxesLabel -> {"\[Tau]", "Coordinate"}, PlotLegends -> {"t", "r", "\[Theta]", "\[CurlyPhi]"}, PlotRange -> {0, 30}], Show[ParametricPlot[ Evaluate[{xyzsoln[[1]], xyzsoln[[2]]} // Flatten], {\[Tau], 0, max\[Tau]}, AspectRatio -> 1, PlotStyle -> Gray], Graphics[{Red, Circle[{0, 0}, 2]}]]} Shape FUNCTION; (1)/(((r)/((Sqrt[r*r0]))) - 1) Integrate[(1)/(((r)/((Sqrt[r*r0]))) - 1), r] max\[Tau] = 2000; ivs = {\[Minus]0.08, .035, .0359}; ics = {0, 2, Pi/4, 0.2}; xyzsoln = sphslnToCartsln[computeSoln[max\[Tau], ivs, ics]]; angle = ParametricPlot3D[ Evaluate[Re[xyzsoln] // Flatten], {\[Tau], 0, max\[Tau]}, AxesLabel -> {x, y, z}, DisplayFunction -> Identity]; sphhoriz = SphericalPlot3D[{1, 2, 3, 4}, {\[Theta], 0, Pi}, {\[Phi], 0, 2*Pi}, PlotRange -> {{-10, 10}, {-10, 10}, {-30, 30}}, BoxRatios -> {1, 1, 1}, PlotTheme -> "Classic", BoxRatios -> 2, Mesh -> None, PlotStyle -> Opacity[0.2]]; sphhoriz2 = SphericalPlot3D[{1, 2, 3, 4}, {\[Theta], 0, Pi}, {\[Phi], 0, (3/2)*Pi}, PlotRange -> {{-10, 10}, {-10, 10}, {-30, 30}}, BoxRatios -> {1, 1, 1}, PlotTheme -> "Classic", BoxRatios -> 2, Mesh -> None, PlotStyle -> Opacity[0.3]]; Show[angle, sphhoriz, DisplayFunction -> $ DisplayFunction, PlotRange -> {{-10, 10}, {-10, 10}, {-10, 10}}] Show[angle, sphhoriz2 , DisplayFunction -> $ DisplayFunction, PlotRange -> {{-10, 10}, {-10, 10}, {-10, 10}}] `