Hamiltonian Method

Part 2: Point Dipole

M. E. Mauel
Columbia University

Off[General::spell1,General::spell];

Introduction

The normalized Hamiltonian is

h = (pz^2 + pr^2)/2 + v

Simple Magnetic Mirror

b[r_,z_] = Sqrt[(r^2) + 4 (z^2)]/((((r^2) + (z^2)))^2)

vDipole = ( 1/r - r/(z^2 + r^2)^(3/2) )^2/2

Plot[Evaluate[vDipole /. z -> 0],  {r, 0.001, 4.0}, PlotRange -> {0, 0.1}]

ContourPlot[vDipole, {r, 0.001, 2.0},  {z, -1, 1},ContourShading -> False, Contours -> {0.001, 0.01, 0.02, 0.03, 0.04}, AspectRatio -> Automatic, PlotPoints -> 40]

hMirror = h /. v -> vDipole

Hamilton's Equations of Motion

Radial Force

rForce[r_,z_]=-D[hMirror,r]// Simplify

Axial Force

zForce[r_,z_]=-D[hMirror,z]// Simplify

Full Equations of Motion

prEq = D[pr[t],t] == rForce[r[t],z[t]];
pzEq = D[pz[t],t] == zForce[r[t],z[t]];

The equations of motion for {r, z} are easy:

rEq = D[r[t],t] == pr[t];
zEq = D[z[t],t] == pz[t];

ϕEq = D[ϕ[t],t] == 1/r[t]^2 - 1/(r[t]^2 + z[t]^2)^(3/2)

The list of unknown functions are:

p = {r[t], ϕ[t], z[t], pr[t], pz[t]};

The initial conditions are:

initial = {r[0] == r0, ϕ[0] == ϕ0, z[0] == z0, pr[0] == pr0, pz[0] == pz0}

The full equations of motions (12 equations in all!) are:

eqns = {rEq,ϕEq, zEq, prEq, pzEq} ~ Join~ initial

Sanple Solution (Low Energy)

initial = { ϕ0 -> 0.0,z0 -> 0.0, r0 -> 1.0, pz0 -> 0.01, pr0 -> 0.05} // N

energy = ((pz0^2) + (pr0^2))/2 /. initial

wellPlot = ContourPlot[-vDipole, {r, 0.001, 2.0},  {z, -1, 1},ContourShading -> True, Contours -> {-energy}, AspectRatio -> Automatic, PlotPoints -> 50];

tmax = 400; sol1 = NDSolve[Evaluate[eqns /.initial],p, {t, 0, tmax}]

Viewing the orbit in the {r, z} plane

p1 = ParametricPlot[Evaluate[{r[t],z[t]} /. sol1], {t, 0, tmax}, PlotPoints-> 500, Frame-> True, AspectRatio -> Automatic]

Show[wellPlot,p1];

Viewing the orbit in the {x, y} plane

pxy1 = ParametricPlot[Evaluate[{r[t] Cos[ϕ[t]], r[t] Sin[ϕ[t]]} /. sol1], {t, 0, tmax}, PlotPoints-> 500, Frame-> True, AspectRatio -> Automatic]

p1 = ParametricPlot[Evaluate[{r[t] Cos[ϕ[t]], r[t] Sin[ϕ[t]]} /. sol1], {t, 0, 6.0}, PlotPoints-> 500, Frame-> True, AspectRatio -> Automatic]

Viewing the orbit in three-dimensions

p2 = ParametricPlot3D[Evaluate[{r[t] Cos[ϕ[t]], r[t] Sin[ϕ[t]], z[t]} /. sol1], {t, 0, tmax},PlotPoints-> 500]

Kinetic Energy

dϕ[t_] = D[ϕ[t]/. First[sol1],t]

Plot[Evaluate[(((pr[t])^2) + ((pz[t])^2) + (((r[t] dϕ[t]))^2))/2 /. sol1], {t, 0, tmax}, PlotPoints->500]

(Approximate) Magnetic Moment

Plot[Evaluate[(((pr[t])^2) + (((r[t] dϕ[t]))^2))/b[r[t],z[t]] /. initial/. sol1], {t, 0, tmax}, PlotPoints->500]

Fourier Analysis of Trajectory

Plot[Evaluate[r[t] /. sol1], {t, 0, tmax}, PlotPoints->500]

Plot[Evaluate[r[t] /. sol1], {t, 0, tmax/10}, PlotPoints->500]

Plot[Evaluate[z[t] /. sol1], {t, 0, tmax}, PlotPoints-> 100]

Sanple Solution (Higher Energy)

initial = { ϕ0 -> 0.0,z0 -> 0.0, r0 -> 1.0, pz0 -> 3 *0.01, pr0 -> 3 *0.05} // N

energy = ((pz0^2) + (pr0^2))/2 /. initial

wellPlot = ContourPlot[-vDipole, {r, 0.001, 2.0},  {z, -1, 1},ContourShading -> True, Contours -> {-energy}, AspectRatio -> Automatic, PlotPoints -> 50];

tmax = 400; sol1 = NDSolve[Evaluate[eqns /.initial],p, {t, 0, tmax}, MaxSteps -> 20000]

Viewing the orbit in the {r, z} plane

p1 = ParametricPlot[Evaluate[{r[t],z[t]} /. sol1], {t, 0, tmax}, PlotPoints-> 500, Frame-> True, AspectRatio -> Automatic];

Show[wellPlot,p1]

Viewing the orbit in the {x, y} plane

pxy1 = ParametricPlot[Evaluate[{r[t] Cos[ϕ[t]], r[t] Sin[ϕ[t]]} /. sol1], {t, 0, tmax}, PlotPoints-> 500, Frame-> True, AspectRatio -> Automatic]

Viewing the orbit in three-dimensions

p2 = ParametricPlot3D[Evaluate[{r[t] Cos[ϕ[t]], r[t] Sin[ϕ[t]], z[t]} /. sol1], {t, 0, tmax},PlotPoints-> 500]

Kinetic Energy

dϕ[t_] = D[ϕ[t]/. First[sol1],t]

Plot[Evaluate[(((pr[t])^2) + ((pz[t])^2) + (((r[t] dϕ[t]))^2))/2 /. sol1], {t, 0, tmax}, PlotPoints->500]

(Approximate) Magnetic Moment

Plot[Evaluate[(((pr[t])^2) + (((r[t] dϕ[t]))^2))/b[r[t],z[t]] /. initial/. sol1], {t, 0, tmax}, PlotPoints->500]

Fourier Analysis of Trajectory

Plot[Evaluate[r[t] /. sol1], {t, 0, tmax}, PlotPoints->500]

Plot[Evaluate[r[t] /. sol1], {t, 0, tmax/10}, PlotPoints->500]

Plot[Evaluate[z[t] /. sol1], {t, 0, tmax}, PlotPoints-> 100]

Sanple Solution (Even Higher Energy)

initial = { ϕ0 -> 0.0,z0 -> 0.0, r0 -> 1.0, pz0 -> 4 *0.01, pr0 -> 4 *0.05} // N

energy = ((pz0^2) + (pr0^2))/2 /. initial

wellPlot = ContourPlot[-vDipole, {r, 0.001, 2.0},  {z, -1, 1},ContourShading -> True, Contours -> {-energy}, AspectRatio -> Automatic, PlotPoints -> 50];

tmax = 250; sol1 = NDSolve[Evaluate[eqns /.initial],p, {t, 0, tmax}, MaxSteps -> 20000]

Viewing the orbit in the {r, z} plane

p1 = ParametricPlot[Evaluate[{r[t],z[t]} /. sol1], {t, 0, tmax}, PlotPoints-> 500, Frame-> True, AspectRatio -> Automatic];

Show[wellPlot,p1]

Viewing the orbit in the {x, y} plane

pxy1 = ParametricPlot[Evaluate[{r[t] Cos[ϕ[t]], r[t] Sin[ϕ[t]]} /. sol1], {t, 0, tmax}, PlotPoints-> 500, Frame-> True, AspectRatio -> Automatic]

Viewing the orbit in three-dimensions

p2 = ParametricPlot3D[Evaluate[{r[t] Cos[ϕ[t]], r[t] Sin[ϕ[t]], z[t]} /. sol1], {t, 0, tmax},PlotPoints-> 500]

Kinetic Energy

dϕ[t_] = D[ϕ[t]/. First[sol1],t]

Plot[Evaluate[(((pr[t])^2) + ((pz[t])^2) + (((r[t] dϕ[t]))^2))/2 /. sol1], {t, 0, tmax}, PlotPoints->500]

(Approximate) Magnetic Moment

Plot[Evaluate[(((pr[t])^2) + (((r[t] dϕ[t]))^2))/b[r[t],z[t]] /. initial/. sol1], {t, 0, tmax}, PlotPoints->500]

Fourier Analysis of Trajectory

Plot[Evaluate[r[t] /. sol1], {t, 0, tmax}, PlotPoints->500]

Plot[Evaluate[r[t] /. sol1], {t, 0, tmax/10}, PlotPoints->500]

Plot[Evaluate[z[t] /. sol1], {t, 0, tmax}, PlotPoints-> 100]

Summary

The computer allows computation of charged particle motion even when the small gyroradius approximation becomes invalid.

As the gyroradius approaches the gradient scale length of the magnetic field, then the magnetic moment, μ is no longer constant.