TL;DR: Yes, helical motion around a ringworld is possible. However, it is far from uniform at larger distances (≥ 0.04 AU).
Summary of results:
For a toroidal ringworld with mass
$M_R = 3 M_\text{star} = 3M_\odot$, central radius $ a = 1 \text{ AU}$ and
inner radius $b = 10^{-4} a \simeq 15000\text{ km}$
(hence density $\rho \sim 8800\text{ kg/m}^3$, compare with
$\rho_{\text{Fe}}\sim 7800\text{ kg/m}^3$), the relation between the mean
separation from the central ring of the ringworld and the time
period is given by the following graph:
$\lambda a$ is the distance from the central ring. The equations are fit using the average (square) value of $\lambda$ over 1 year. The median (pentagon) and launch (circle) values of $\lambda$ are also shown. The error bars indicate the minimum and maximum value of $\lambda$. Black → no central star, orange → $M_{\text{star}} = M_{\odot}$.
The trajectory is indeed helical.
We can look at the projected cross-section below. The center of the ring is at $x=-1$ (scaled by $a$, not shown), the central ring of the torus is at $x=0$. It is clear that the "inner distances" are bigger than the "outer distances" as one would naively expect. As earlier, black → no central star, orange → $M_{\text{star}} = M_{\odot}$. The first graph is translucent so you can see both cross-sections by zooming in.
Physics:
For simplicity, consider a point particle; if the satellite is too large
(how large?) there would be complications due to the Roche limit etc. Replace the ringworld (torus for volume calculation) with its central ring for all other calculations.
Let the particle be launched from $(x,y,z)=((1-\lambda) a, 0,0)$.
We only consider the regime $0.01 \leq \lambda \leq 0.1$. The lower limit
prevents the particle from recognizing that the ringworld has been approximated by a ring and the upper limit prevents its orbit from being perturbed
substantially by the star.
Suppose initially the radial velocity is zero.
The tangential velocity for orbit around the star (initially $v_y$)
should roughly be $\sqrt{G M_{\text{star}}/((1-\lambda) a)}$.
The orbital velocity around the ringworld (initially $v_z$) should roughly be
$\sqrt{g((1-\lambda) a,0,0) \lambda a}$ where $g(x,y,z)$ is the magnitude of the net gravitational field (a.k.a. acceleration due to gravity) as a function of position.
We're kind of stuck without a number for $g((1-\lambda) a,0,0)$.
The electric potential for a ring of charge $Q$ is given as (Ref. 1):
$$
V(r,\phi,z) = \frac{1}{4\pi\epsilon_0}\frac{Q}{2\pi}\int_0^{2\pi}\frac{d\phi'}{\sqrt{r^2-2r a\cos(\phi-\phi')+a^2+z^2}}
$$
It is easy to obtain the potential for a graviational ring by substituting $\epsilon_0 \rightarrow -1/(4\pi G)$ and $Q\rightarrow M$ in the equation. One can take the gradient (with a $-$ sign) and find the field numerically.
Or one could go a few steps further and actually calculate everything.
Using the equation shown earlier for one value of $M_R$ (orange in graph), one can extrapolate to other values of $M_R$ using $1/T\propto\sqrt{M_R}$ (as long as $M_R$ is not much smaller than $M_S$) as $$T = 2\pi R/v \sim 2\pi R/\sqrt{g R} \sim 1/\sqrt{M_R}$$
where $R = \lambda a$.
Implementation (Mathematica):
(Everything is in SI units unless mentioned otherwise.)
First we set up the constants. The elliptic integral for $V$ earlier is somewhat nasty and takes a while to simplify, so I simplified it once and replaced the definition with the output of the simplification.
G = 6.674 10^-11;
EarthMass = 5.9722 10^24;
SolarMass = 333000 EarthMass;
RingMass = 3 SolarMass;
AU = 1.508 10^11;
a = 1 AU;
b = 10^-4 a;
day = 24*3600 // N;
year = 365.25 day;
\[Rho] = RingMass/((2 \[Pi] a) (\[Pi] b^2)) (* roughly 8800, Fe \[Rule] 7800 *)
(* Math *)
VRing[r_, \[Phi]_, z_,
MR_] = -G MR/(2 \[Pi]) ((2 Sqrt[(
a^2 + r^2 + z^2 - 2 a r Cos[\[Phi]])/((a - r)^2 +
z^2)] (EllipticF[\[Pi] - \[Phi]/2, -((
4 a r)/((a - r)^2 + z^2))] +
EllipticF[\[Phi]/2, -((4 a r)/((a - r)^2 + z^2))]))/(Sqrt[
a^2 + r^2 + z^2 - 2 a r Cos[\[Phi]]]));
VRingxyz[x_, y_, z_, MR_] =
TransformedField["Polar" -> "Cartesian",
VRing[r, \[Phi], z, MR], {r, \[Phi]} -> {x, y}];
Vtot[x_, y_, z_, MR_, MS_] = -G MS/Norm[{x, y, z}] +
VRingxyz[x, y, z, MR];
gRing[x_, y_, z_, MR_] = -Grad[VRingxyz[x, y, z, MR], {x, y, z}];
gtot[x_, y_, z_, MR_,
MS_] = -Grad[Vtot[x, y, z, MR, MS], {x, y, z}] /.
Abs[p_] Abs'[p_] -> p;
gtotmag[x_, y_, z_, MR_, MS_] = Norm[gtot[x, y, z, MR, MS]];
Let's do a quick sanity check and see if the gravitational field is as expected.
imgWidth = 2160;
plotAndExport[fname_,
plot_] := (Export[NotebookDirectory[] <> fname,
Rasterize[plot, ImageSize -> imgWidth]]; plot);
fieldPlotXLim = 1.5/Sqrt[2]; fieldPlotYLim = fieldPlotXLim;
splot = plotAndExport["field.png", #] &@
Show[StreamPlot[
Chop@(gtot[x1 a, y1 a, 0, RingMass, SolarMass][[1 ;; 2]])
, {x1, -fieldPlotXLim, fieldPlotXLim}, {y1, -fieldPlotYLim,
fieldPlotYLim}
, BaseStyle -> {FontSize -> 24}]];
Looks alright. The first graph shows the "flow" of the field (the arrows sizes are not correct). The magnitude of the field along the $x$-axis is also shown.
Now we implement the solvers for the particle trajectory.
(* Trajectory solvers with initial conditions *)
xi[\[Lambda]_] := a (1 - \[Lambda]);
yi[\[Lambda]_] := 0.;
zi[\[Lambda]_] := 0.;
vxi[\[Lambda]_] := 0.;
vyi[\[Lambda]_] := Sqrt[G SolarMass/Abs[xi[\[Lambda]]]];
vzi[\[Lambda]_, MR_] :=
Sqrt[Abs[a - xi[\[Lambda]]] Norm@
gRing[xi[\[Lambda]], yi[\[Lambda]], zi[\[Lambda]], MR]];
ringSol[\[Lambda]_, MR_, time_] := NDSolve[
Flatten@{Thread[{xs''[t], ys''[t], zs''[t]} ==
gRing[xs[t], ys[t], zs[t], MR]],
xs'[0] == vxi[\[Lambda]], ys'[0] == vyi[\[Lambda]],
zs'[0] == vzi[\[Lambda], MR],
xs[0] == xi[\[Lambda]], ys[0] == yi[\[Lambda]],
zs[0] == zi[\[Lambda]]},
{xs, ys, zs}, {t, 0, time}];
xiFull[\[Lambda]_] := xi[\[Lambda]];
yiFull[\[Lambda]_] := yi[\[Lambda]];
ziFull[\[Lambda]_] := zi[\[Lambda]];
vxiFull[\[Lambda]_] := vxi[\[Lambda]];
vyiFull[\[Lambda]_, MS_] := Sqrt[G MS/Abs[xi[\[Lambda]]]];
vziFull[\[Lambda]_, MR_, MS_] :=
Sqrt[Abs[a - xi[\[Lambda]]] Norm@
gtot[xi[\[Lambda]], yi[\[Lambda]], zi[\[Lambda]], MR, MS]];
fullSol[\[Lambda]_, MR_, MS_, time_, \[Epsilon]_] := NDSolve[
Flatten@{Thread[{xs''[t], ys''[t], zs''[t]} ==
gtot[xs[t], ys[t], zs[t], MR, MS]]
, xs'[0] == vxiFull[\[Lambda]], ys'[0] == vyiFull[\[Lambda], MS],
zs'[0] == (1 + \[Epsilon]) vziFull[\[Lambda], MR, MS]
, xs[0] == xiFull[\[Lambda]], ys[0] == yiFull[\[Lambda]],
zs[0] == ziFull[\[Lambda]]}
, {xs, ys, zs}, {t, 0, time}
];
appendVelocities[solution_] :=
Append[solution, {vx -> xs', vy -> ys', vz -> zs'} /. solution]
We will need a bunch of functions to analyze the time period.
(* Examining the period T of rotation about the ring *)
(* findPeriod and reconstruct copied from \
https://mathematica.stackexchange.com/a/38221/9332 *)
findPeriod[data_, threshold_] :=
Module[{fs, s1, s = {}, i, a0f, af, pf, pos, fr, frpos, fdata,
fdatac, n, per}, n = Length[data];
fs = Fourier[data];
s1 = Drop[fs, -Floor[Length[fs]/2]];
For[i = 1, i < Length[s1], i++,
If[Abs[fs][[i + 1]] > threshold, AppendTo[s, i + 1]]];
a0f = Abs[fs[[1]]]/Sqrt[n];
af = 2/Sqrt[n] Abs[fs][[s]];
pf = Arg[fs][[s]];
{a0f, Transpose[{s, af, pf}]}]
reconstruct[data_, fp_] := Module[{n}, n = Length[data];
Show[ListLinePlot[data, PlotStyle -> Black],
Plot[fp[[1]] +
Sum[fp[[2, j, 2]] Cos[
2 Pi (fp[[2, j, 1]] - 1)/n t - fp[[2, j, 3]]], {j, 1,
Length[fp[[2]]]}], {t, 0, n}, PlotStyle -> Red]]];
getOrbitPeriod[solution_, totalTime_, timeStep_] := Module[{data},
data =
Flatten@Table[
zs[t timeStep] /. solution, {t, 0, totalTime/timeStep}];
(* Not strictly correct as there are many frequencies but good \
enough for first approximation *)
totalTime/(timeStep Sort[
findPeriod[data, 10^8][[2]], #1[[2]] > #2[[2]] &][[1, 1]])];
(* The period T is observed to be linear in \[Lambda] *)
\
\[Lambda]TFit[\[Lambda]list_, Tlist_] :=
LinearModelFit[
Transpose@{\[Lambda]list, Tlist}, \[Lambda], \[Lambda]];
setGraphFontSize = BaseStyle -> {FontSize -> 12};
graphLineWidth = 0.003;
graphMarkerLineWidth = 0.005;
graphMarkerSize = 6;
opacity = 0.5;
polygonMarker[color_, n_] :=
Graphics[{EdgeForm[{Thickness -> graphMarkerLineWidth, color}],
FaceForm[None], Polygon[CirclePoints@n]},
ImageSize -> graphMarkerSize];
coloredListPlot[x_, y_, color_, PM_] :=
ListPlot[Transpose@{x, y}, PlotStyle -> color, PlotMarkers -> PM];
Needs["ErrorBarPlots`"]
\[Lambda]TFitGraph[{\[Lambda]list_, min\[Lambda]_, max\[Lambda]_,
mean\[Lambda]_, median\[Lambda]_}, Tlist_, color_] :=
Module[{model = \[Lambda]TFit[mean\[Lambda], Tlist]},
Show[
Plot[Normal[model], {\[Lambda], 0.01, Max[mean\[Lambda]]}
, PlotStyle -> {color, Dashed, Thickness -> graphLineWidth},
AxesLabel -> {"\[Lambda]", "T (days)"}
, PlotLegends -> SwatchLegend[{color}, {Normal[model]}]
, Evaluate@setGraphFontSize,
PlotRange -> {{0, Automatic}, {0, Automatic}}]
, ErrorListPlot[
(({{#1, #4}, ErrorBar[{#2 - #1, #3 - #1}, {0, 0}]} &) @@ # &) /@
Transpose@{mean\[Lambda], min\[Lambda], max\[Lambda], Tlist}
, PlotStyle -> {color, Thickness -> graphLineWidth},
PlotMarkers -> polygonMarker[color, 4]]
, coloredListPlot[\[Lambda]list, Tlist,
color, {Automatic, graphMarkerSize}]
, coloredListPlot[median\[Lambda], Tlist, color,
polygonMarker[color, 5]]
]]
Finally, we actually run the solvers and see the data.
(* Actually run simulations *)
ringSolutionTime = year;
ring\[Lambda]list = Range[0.01, 0.1, 0.01];
AbsoluteTiming[
ringSolutions =
Flatten@appendVelocities@ringSol[#, RingMass, ringSolutionTime] & /@
ring\[Lambda]list
][[1]]
ringPeriods =
getOrbitPeriod[#, ringSolutionTime, day] & /@ ringSolutions;
{ringMaxDist, ringMinDist, ringMeanDist, ringMedianDist} =
Transpose[distCalc[#, ringSolutionTime, day/24] & /@ ringSolutions];
TableForm@{ring\[Lambda]list, ringMaxDist, ringMinDist,
ringMeanDist, ringMedianDist}
fullSolutionTime = year;
full\[Lambda]list = ring\[Lambda]list + 0.005;
AbsoluteTiming[
fullSolutions =
Flatten@appendVelocities@
fullSol[#, RingMass, SolarMass, fullSolutionTime, 0] & /@
full\[Lambda]list
][[1]]
fullPeriods = getOrbitPeriod[#, fullSolutionTime, day] & /@ fullSolutions;
{fullMaxDist, fullMinDist, fullMeanDist, fullMedianDist} =
Transpose[distCalc[#, fullSolutionTime, day/24] & /@ fullSolutions];
TableForm@{full\[Lambda]list, fullMaxDist, fullMinDist,
fullMeanDist, fullMedianDist}
Making the $T$ vs $\lambda$ plot and seeing the trajectory (graphs in summary).
plotAndExport["Tvl.png", #] &@
Show[
\[Lambda]TFitGraph[{ring\[Lambda]list, ringMinDist, ringMaxDist,
ringMeanDist, ringMedianDist}, ringPeriods, Black]
, \[Lambda]TFitGraph[{full\[Lambda]list, fullMinDist, fullMaxDist,
fullMeanDist, fullMedianDist}, fullPeriods, Orange]
, PlotRange -> {{0, Automatic}, {0, Automatic}}
]
plotAndExport["traj.png", #] &@
Show[GraphicsGrid[{{
trajectory[full\[Lambda]list[[1]], fullSolutions[[1]],
fullSolutionTime/7]
, trajectory[full\[Lambda]list[[1]], fullSolutions[[1]],
fullSolutionTime]
}}]]
References:
- http://physics.oregonstate.edu/portfolios/Activities/EMActivities/ElectricPotentialRing/RingVSolutions070701.pdf
