Plot dipole magnetic field lines in 3D space

Question

I would like to reproduce the neutron star (NS) dipole magnetic field in the 3D space. I have also a constraint, the line should start from the inner edge of an accretion disk and go toward the pole of the NS magnetic axis. In the attached figure, I have realized the simple case where the magnetic axis coincides with the rotation axis, which is orthogonal to the plane, where the disk lies. I would like to extend this configuration to the case where the magnetic axis is inclined by an angle beta with respect to the rotation axis. In this particular case, I have to draw the magnetic field lines that are closer to the magnetic axis. Therefore some of them will be connected with the North pole, while the other with the South pole. I report below my code in Mathematica:

RS = 6; (*NEUTRON STAR RADIUS*)
RM = 40; (*MAGNETOSPHERIC RADIUS=INNER EDGE OF THE DISK*)

ROUT = RM + 10; (*OUTER EDGE OF THE DISK*)
\[Beta] = 
  30*Degree; (* INCLINATION OF THE MAGNETIC AXIS*)

i = 60*Degree; (*INCLINATION OF THE OBSERVER*)

NS = SphericalPlot3D[RS, {x, 0, 2*Pi}, {y, 0, Pi}, 
   ColorFunction -> Black, PlotStyle -> {Gray, Opacity[1]}, 
   Mesh -> None, PlotPoints -> 50]; (*NEUTRON STAR*)

DISK = ParametricPlot3D[{r*Sin[Pi/2]*Cos[\[Phi]], 
    r*Sin[Pi/2]*Sin[\[Phi]], r*Cos[Pi/2]}, {r, RM, ROUT}, {\[Phi], 0, 
    2*Pi}, Mesh -> None, PlotPoints -> 50]; (*DISK*)

raxis = Graphics3D[{Red, Arrowheads[0.02], 
    Arrow[Tube[{{0, 0, -10}, {0, 0, 10}}, 
      0.02]]}] ; (*ROTATION AXIS*)
maxis = 
 Graphics3D[{Blue, Arrowheads[0.02], 
   Arrow[Tube[{{10*Sin[Pi + \[Beta]], 0, 
       10*Cos[Pi + \[Beta]]}, {10*Sin[\[Beta]], 0, 10*Cos[\[Beta]]}}, 
     0.02]]}] ; (*MAGNETIC AXIS AXIS*)

OBSERVER = 
 Graphics3D[{PointSize[Large], 
   Point[{0, ROUT*Sin[i], ROUT*Cos[i]}]}];
lineofsight = 
  Graphics3D[{Dashed, 
    Line[{{0, 0, 0}, {0, ROUT*Sin[i], 
       ROUT*Cos[i]}}]}]; (*OBSERVER AT INFINITY*)

NN = 100; (*NUMBER OF MAGNETIC FIELD LINES*)

f[j_] := -Pi/2 + Pi*j/NN; (*PHI COORDINATE*)

PN[x_, j_] := (RM*Sin[x]^2)*{Sin[x]*Cos[f[j]], Sin[x]*Sin[f[j]], 
   Cos[x]}; (*POINTS OF THE MAGNETIC FIELD LINES*)

ld = Table[
  ParametricPlot3D[PN[x, j], {x, 0, Pi/2}], {j, 0, NN, 
   1}]; (*MAGNETIC FIELD LINES DISK-NORTH POLE*)
lind = 
 Table[ParametricPlot3D[(-PN[x, j]), {x, 0, Pi/2}], {j, 0, NN, 
   1}]; (*MAGNETIC FIELD LINES DISK-SOUTH POLE*)

Show[NS, DISK, raxis, maxis, OBSERVER, lineofsight, ld, lind, 
 PlotRange -> All, Axes -> False, Boxed -> False, ImageSize -> 800]

I have some questions: 1) Is it correct what I have plotted? 2) How do I extend this plot for an inclined magnetic field? Thank you very much in advance for your great help.

Answer 1

1) I think it's correct.

2) Did I understand you right that you want to get it this way:

Does it look better?

You will get the solution soon but I want to solve it on my own first. And I want to combine it in a readable way. I did a lot of steps to reach it.

I'm really sorry. I figured out that the amplitude factor 'Sin[[Phi]]^2' has to be adapted by the inclination angle like I did here:

    li = Graphics[Table[Line[{{0, 0}, 1.2 {Cos[\[Beta]], Sin[\[Beta]]}}], 
    {\[Beta], 0 \[Degree], 60 \[Degree], 15 \[Degree]}]];
    pp = ParametricPlot[{{Sin[\[Phi]], Cos[\[Phi]]}, 
    Table[Sin[\[Phi]]^2 {Sin[\[Phi]], Cos[\[Phi]]}/
    Cos[\[Beta]]^2, {\[Beta], 0, 60 \[Degree], 
    15 \[Degree]}]}, {\[Phi], 0, \[Pi]/2}, 
    PlotLabel ->
    Style["Dipole field vs. inclination angle", Larger, Black]];

which gives a result of:

But this doesn't work if your are not perpendicular to the magnetic axis. Therefore I figured out that the amplitude should be:

    ((Subscript[\[ScriptCapitalR], di]*
    Sin[\[Alpha]]^2)/(Cos[\[Beta]] Cos[\[Phi]]^2 + Sin[\[Phi]]^2)^2)

Unfortunately there is still an error for great inclination angles of [Beta]. It even looks better like you can see here:

Maybe you will continue a bit more and see what I have done so far. Here it is:

Subscript[\[ScriptCapitalR], ns] = 
  6; (*Radius of neutron star.*)
\[Beta] = 
 30 \[Degree]; (*Inclinition of the magnetic axis.*) 
\[Gamma] = 
 60 \[Degree]; (*Inclination of the observer.*)
Subscript[\
\[ScriptCapitalR], di] = 30; (*Inner disk radius \[LongEqual] \
magnetospheric radius.*)
Subscript[\[ScriptCapitalR], do] = 
 Subscript[\[ScriptCapitalR], di] + 10; (*Outer disk radius.*)
rmx = 
 RotationMatrix[\[Beta], {0, 1, 
   0}]; (*Rotation matrix for arrow & field.*)
nst = 
 Graphics3D[{Gray, 
   Sphere[{0, 0, 0}, Subscript[\[ScriptCapitalR], 
    ns]]}]; (*Neutron star.*)
axs = 
 Graphics3D[{Red, Line[{{0, 0, 0}, {25, 0, 0}}], Green, 
   Line[{{0, 0, 0}, {0, 25, 0}}], Blue, 
   Line[{{0, 0, 0}, {0, 0, 
      25}}]}]; (*Koordiante center & direction.*)
arr = 
 Arrow[Tube[{{0, 0, -15}, {0, 0, 15}}, 
   0.2]]; (*Arrow for axis.*)
rax = 
 Graphics3D[{Orange, arr}]; (*Rotation axis.*)
max = 
 Graphics3D[{Blue, 
   GeometricTransformation[arr, rmx]}]; (*Magnetic axis.*)
pos = {0, 
  Subscript[\[ScriptCapitalR], do]*Sin[\[Gamma]], 
  Subscript[\[ScriptCapitalR], do]*
   Cos[\[Gamma]]}; (*Observer position.*)
obs = 
 Graphics3D[{PointSize[Large], Point[pos]}]; (*Observer.*)
los = 
 Graphics3D[{Dashed, 
   Line[{{0, 0, 0}, 
     pos}]}]; (*Line of sight.*)
nof = 40; (*Number of field lines.*)

fld[\[Alpha]_, \[Phi]_] := ((
   Subscript[\[ScriptCapitalR], di]*
    Sin[\[Alpha]]^2)/(Cos[\[Beta]] Cos[\[Phi]]^2 + 
     Sin[\[Phi]]^2)^2)*{Sin[\[Alpha]]*Cos[\[Phi]], 
   Sin[\[Alpha]]*Sin[\[Phi]], Cos[\[Alpha]]};
fdl = Table[
  fld[\[Alpha], \[Phi]], {\[Phi], -\[Pi]/2, \[Pi]/2, \[Pi]/
   nof}, {\[Alpha], 0, \[Pi]/2, \[Pi]/50}]; (*Field data +.*)
fdr = 
 fdl.ScalingMatrix[{-1, 1, -1}]; (*Field data -.*)
ldb = {Line[fdl], 
  Line[fdr]}; (*Line data both \[PlusMinus].*)
gld = 
 Graphics3D[{Blue, 
   GeometricTransformation[ldb, 
    rmx]}]; (*Field lines as graphic.*)
disk = 
 ParametricPlot3D[{r*Cos[\[Phi]], r*Sin[\[Phi]], 0}, {r, 
   Subscript[\[ScriptCapitalR], di], Subscript[\[ScriptCapitalR], 
   do]}, {\[Phi], 0, 2 \[Pi]}, Mesh -> None, PlotPoints -> 50];
plt = Show[gld, axs, nst, disk, rax, max, los, obs, PlotRange -> All, 
  Axes -> False, Boxed -> False, ImageSize -> 500]

Answer 2

I'm sorry it took a while to figure it out. The amplitude of the field had to be adapted to the magnetosphere radius.

Subscript[\[ScriptCapitalR],ns] = 6; (*Radius of neutron star.*)
\[Mu] = 25 \[Degree]; (*Inclination of the magnetic axis.*) 
\[Omega] = 60 \[Degree]; (*Inclination of the observer.*)
Subscript[\[ScriptCapitalR], di] = 30; (*Inner disk radius \[LongEqual] magnetospheric radius.*)
Subscript[\[ScriptCapitalR], do] = Subscript[\[ScriptCapitalR], di] + 10; (*Outer disk radius.*)
Subscript[\[ScriptCapitalR], dm] = Mean[{Subscript[\[ScriptCapitalR], do], Subscript[\[ScriptCapitalR],di]}]; (*Mean disk radius.*)
Subscript[\[ScriptCapitalR], dd] = (Subscript[\[ScriptCapitalR], do]-Subscript[\[ScriptCapitalR], di])/2; (*Delta disk radius.*)
(*Angle between the field line and the magnetosphere.*)
calc\[Upsilon][\[Beta]_, \[Gamma]_] := 2 ArcTan[1/2 Cot[\[Beta]/2] Sec[\[Gamma]] (-1 + Tan[\[Beta]/2]^2 + Cos[\[Gamma]] Tan[\[Beta]/2] Sqrt[Cot[\[Beta]/2]^2 Sec[\[Gamma]]^2 (1 - 2 Tan[\[Beta]/2]^2 + 4 Cos[\[Gamma]]^2 Tan[\[Beta]/2]^2 + Tan[\[Beta]/2]^4)])];
Subscript[\[ScriptCapitalA], f][\[Gamma]_] := 1/Cos[calc\[Upsilon][\[Mu], \Gamma]]]^2; (*Field amplitude correction.*)
rmx = RotationMatrix[\[Mu], {0, 1, 0}]; (*Rotation matrix for arrow & field.*)
nst = Graphics3D[{Gray, Sphere[{0, 0, 0}, Subscript[\[ScriptCapitalR], ns]]}]; (*Neutron star.*)
axs = Graphics3D[{Red, Line[{{0, 0, 0}, {25, 0, 0}}], Green, Line[{{0, 0, 0}, {0, 25, 0}}], Blue, Line[{{0, 0, 0}, {0, 0, 25}}]}]; (*Coordinate center & direction for orientation.*)
arr = Arrow[Tube[{{0, 0, -15}, {0, 0, 15}}, 0.2]]; (*Arrow for axis.*)
rax = Graphics3D[{Orange, arr}]; (*Rotation axis.*)
max = Graphics3D[{Blue, GeometricTransformation[arr, rmx]}]; (*Magnetic axis.*)
pos = {0, Subscript[\[ScriptCapitalR], do]*Sin[\[Omega]], Subscript[\[ScriptCapitalR], do]*Cos[\[Omega]]}; (*Observer position.*)
obs = Graphics3D[{PointSize[Large], Point[pos]}]; (*Observer.*)
los = Graphics3D[{Dashed, Line[{{0, 0, 0}, pos}]}]; (*Line of sight.*)
nof = 40; (*Number of field lines.*)
fld[\[Alpha]_, \[Phi]_] := Subscript[\[ScriptCapitalA], f][\[Phi]]*Subscript[\[ScriptCapitalR], di]*Sin[\[Alpha]]^2*{Sin[\[Alpha]]*Cos[\[Phi]], Sin[\[Alpha]]*Sin[\[Phi]], Cos[\[Alpha]]};
fdl = Table[fld[\[Alpha], \[Phi]], {\[Phi], -\[Pi]/2.002, \[Pi]/2.002, (0.999 \[Pi])/nof}, {\[Alpha], 0, \[Pi]/2, \[Pi]/50}]; (*Field data +.*)
fdr = fdl.ScalingMatrix[{-1, 1, -1}]; (*Field data -.*)
ldb = {Line[fdl], Line[fdr]}; (*Line data both \[PlusMinus].*)
gld = Graphics3D[{Blue, GeometricTransformation[ldb, rmx]}]; (*Field lines as graphic.*)
disk = ParametricPlot3D[{(35 + 5 Cos[\[Phi]]) Sin[\[Upsilon]], (35 + 5 Cos[\[Phi]]) Cos[\[Upsilon]], 1 Sin[\[Phi]]}, {\[Upsilon], 0, 2 Pi}, {\[Phi], 0, 2 Pi}, PlotStyle -> {Orange}, Mesh -> None];
plt = Show[gld, axs, nst, disk, rax, max, los, obs, PlotRange -> 1.2 {{-Subscript[\[ScriptCapitalR], do], Subscript[\[ScriptCapitalR], do]}, {-Subscript[\[ScriptCapitalR], do], Subscript[\[ScriptCapitalR], do]}, {-Subscript[\[ScriptCapitalR], do], Subscript[\[ScriptCapitalR], do]}}, Axes -> False, Boxed -> False, ImageSize -> 560]

Finally the field lines fit the magnetosphere equator independent of incline.