Magnetic Bottle

In this worksheet we investigate charged-particle trajectories in magnetic fields. The magnetic field of a current loop was calculated in the worksheet Currentloop.mws. The magnetic field lines were calculated using a solver for systems of first-order differential equations. Then the current loop results were extended by combining the fields to form a Helmholtz pair, and a solenoid. Here we investigate charged-particle trajectories by integrating Newton's equations with the Lorentz force. We begin with the study of trajectories in a magnetic field caused by a single current loop.

Subsequently we assemble a magnetic field made up of four current loops. Two of the loops are larger and form a Helmholtz coil. Two loops that are smaller are added to generate a bottleneck. This will be a simple model of a magnetic bottle.

In principle, we work in MKSA (SI) units. We do, however use a charge-over-mass ratio of e / m =1, i.e., we do not specify the mass as 1 kg, and the charge unit as 1.6 * 10^(-19) C. For real-life applications this would have to be corrected. For the qualitative analysis of charged-particle trajectories this is not a problem.

>

Charged-particle trajectories in the field of a current loop

We carry out the solution of Newton's equation for a charged particle with the Lorentz force. The usual textbook discussions deal only with homogeneous fields, where the constant magnetic field imposes circular motion in a plane perpendicular to the magnetic field direction. In general this leads to spiralling motion for particles which possess a velocity component aligned with the field.

The motion of charged particles in inhomogeneous magnetic fields is much more complicated. We are particularly interested in a demonstration that charged particles which move in the direction where the magnetic field becomes stronger are slowed down and turn around, i.e., that they prefer to travel in the direction where the magnetic field becomes weaker. This has implications for cosmic charged particles (picking up energy in galactic magnetic fields), and charged particles travelling from the Sun to the Earth, which spiral towards the poles (and turn around where the field becomes too strong), and collide with each other to produce optical airglow emissions (aurora, or northern lights).

When we calculated B _ x we really obtained the transverse component of the magnetic field, as we worked specifically in the x - z plane. For a solution to Newton's equations we should generate three Cartesian components to be calculated anywhere in space (for any x , y , z ). Given how the Lorentz force acts (at right angles to the velocity and to the B field) we cannot expect the motion to be confined to a plane for the case of a non-homogeneous magnetic field.

We define a procedure to calculate the Cartesian field components for the particular case of a current loop defined in Bxa,Bza. We start over and define the functions needed for the magnetic field components.

> restart; with(linalg): with(plots):

Warning, new definition for norm

Warning, new definition for trace

> Qp:=(Q+R)^2+a^2: Qm:=(Q-R)^2+a^2:

> Qr:=2*sqrt(Q*R/Qp);

> Bxex:=i_L*a/(5*10^(6))/(sqrt(Qp)*Q)*((R^2+Q^2+a^2)*EllipticE(Qr)/Qm-EllipticK(Qr));

> Bxa:=unapply(subs(i_L=1000,R=0.005,Bxex),Q,a):

> Bzex:=i_L/(5*10^(6))/sqrt(Qp)*(EllipticE(Qr)*(R^2-Q^2-a^2)/Qm+EllipticK(Qr));

> Bza:=unapply(subs(i_L=1000,R=0.005,Bzex),Q,a):

To avoid the problems associated with taking the limit at x=0,z=0, we add a tiny constant to the cylindrical coordinate r .

> eta:=10^(-Digits):

> Bfield:=proc(x,y,z) local r,Btr,Bx,By,Bz,cphi,sphi; global Bxa,Bza,eta;

> r:=sqrt(x^2+y^2)+eta;

> Bz:=evalf(Bza(z,r)); Btr:=evalf(Bxa(z,r));

> cphi:=evalf(x/r); sphi:=evalf(y/r);

> By:=Btr*sphi; Bx:=Btr*cphi;

> [Bx,By,Bz]; end:

> Bfield(0.002,0.003,0.004);

> Bfield(0.0,0.00,0.004);

> NE:=diff(x(t),t)=p_x(t),diff(p_x(t),t)=p_y(t)*Bfield(x(t),y(t),z(t))[3]-p_z(t)*Bfield(x(t),y(t),z(t))[2],diff(y(t),t)=p_y(t),diff(p_y(t),t)=p_z(t)*Bfield(x(t),y(t),z(t))[1]-p_x(t)*Bfield(x(t),y(t),z(t))[3],diff(z(t),t)=p_z(t),diff(p_z(t),t)=p_x(t)*Bfield(x(t),y(t),z(t))[2]-p_y(t)*Bfield(x(t),y(t),z(t))[1]:

We start a trajectory in the x - z plane with no momentum in the x -direction, and non-vanishing components in the y - and z -directions.

> IC:=x(0)=0.02,p_x(0)=0,y(0)=0,p_y(0)=0.001,z(0)=0,p_z(0)=0.001;

> sol:=dsolve({NE,IC},{x(t),p_x(t),y(t),p_y(t),z(t),p_z(t)},numeric,output=listprocedure):

Error, (in DEtools/convertsys) unable to convert to an explicit first-order system

This problem is solved in Maple6.

> sol(1);

For faster calculations and to be able to work in Maple5 we introduce our own Runge-Kutta procedure used before for the fieldlines.

> RKn:=proc(f,ic::list,N,x1,h) local xn,yn,ynew,m1,m2,m3,m4,ii,nn;

> xn:=op(1,ic):

> yn:=[]:

> for nn from 1 to N do:

> yn:=[op(yn),op(1+nn,ic)]: od:

> for ii from 1 to round((x1-xn)/h) do:

> m1:=[]:

> for nn from 1 to N do:

> m1:=[op(m1),evalf(f(xn,yn)[nn])]; od:

> m2:=[]:

> for nn from 1 to N do:

> m2:=[op(m2),evalf(f(xn+h/2,evalm(yn+h/2*m1))[nn])]; od:

> m3:=[]:

> for nn from 1 to N do:

> m3:=[op(m3),evalf(f(xn+h/2,evalm(yn+h/2*m2))[nn])]; od:

> xn:=xn+h;

> m4:=[]:

> for nn from 1 to N do:

> m4:=[op(m4),evalf(f(xn,evalm(yn+h*m3))[nn])]; od:

> ynew:=[]:

> for nn from 1 to N do:

> ynew:=[op(ynew),evalf(yn[nn]+h/6*(m1[nn]+2*(m2[nn]+m3[nn])+m4[nn]))]; od: yn:=ynew;

> od:

> [xn,op(ynew)];

> end:

We define a RHS procedure for the problem.

> fn:=proc(t,yn) local x,y,z,p_x,p_y,p_z,BF; global Bfield;

> x:=yn[1]: p_x:=yn[2]: y:=yn[3]: p_y:=yn[4]: z:=yn[5]: p_z:=yn[6]: BF:=Bfield(x,y,z);

> [p_x,p_y*BF[3]-p_z*BF[2],p_y,p_z*BF[1]-p_x*BF[3],p_z,p_x*BF[2]-p_y*BF[1]];

> end:

The argument list contains the external RHS, a list that starts with the initial time (independent variable), and the initial values, the number of dependent variables, the final value of the time (independent variable), and the time-step (step-size of the independent variable).

>

It seems that we can use a large step. (at least intially).

We prepare for a graph of the current loop which causes the magnetic field and the trajectory:

> PL4:=spacecurve([0.005*cos(f),0.005*sin(f),0],f=0..2*Pi,color=blue):

We pick initial conditions for a charged particle located above the current loop without initial momentum in the z -direction, and a small amount of motion in the x - y plane. The reader is encouraged to play with the initial conditions patiently in order to understand the outcomes.

> Rt:=[]: dt:=0.5: x0:=0.002: px0:=-0.0: y0:=eta: py0:=-0.00005: z0:=0.005: pz0:=0.0: for i from 1 to 1000 do:

> t0:=(i-1)*dt; sol:=RKn(fn,[t0,x0,px0,y0,py0,z0,pz0],6,t0+dt,dt); x0:=sol[2]: px0:=sol[3]: y0:=sol[4]: py0:=sol[5]: z0:=sol[6]: pz0:=sol[7]: Rt:=[op(Rt),[x0,y0,z0]]: od:

> PL5:=pointplot3d(Rt,axes=boxed,color=green,labels=[x,y,z]):

> display(PL4,PL5);

> Rt(nops(Rt));

An interesting feature of Loretnz-force trajectories is that they are not symmetric under time reversal. Under time reversal velocities are reversed, and one ought to be able to run through the trajectory from the end to the beginning. This is usually a check for the accuracy of the integration scheme. To accomplish this with the Lorentz force, we have to change the sign of the charge.

> q_p:=-1;

> fn:=proc(t,yn) local x,y,z,p_x,p_y,p_z,BF; global Bfield,q_p;

> x:=yn[1]: p_x:=yn[2]: y:=yn[3]: p_y:=yn[4]: z:=yn[5]: p_z:=yn[6]: BF:=Bfield(x,y,z);

> [p_x,q_p*(p_y*BF[3]-p_z*BF[2]),p_y,q_p*(p_z*BF[1]-p_x*BF[3]),p_z,q_p*(p_x*BF[2]-p_y*BF[1])];

> end:

> Rt:=[]: dt:=0.5: px0:=-px0: py0:=-py0: pz0:=-pz0: for i from 1 to 1000 do:

> t0:=(i-1)*dt; sol:=RKn(fn,[t0,x0,px0,y0,py0,z0,pz0],6,t0+dt,dt); x0:=sol[2]: px0:=sol[3]: y0:=sol[4]: py0:=sol[5]: z0:=sol[6]: pz0:=sol[7]: Rt:=[op(Rt),[x0,y0,z0]]: od:

> PL5:=pointplot3d(Rt,axes=boxed,color=green,labels=[x,y,z]):

> display(PL4,PL5);

> Rt[nops(Rt)];

Our algorithm for the reverse propagation indeed took us back to the beginning of the trajectory.

Let us propagate further to see where the particle goes. We accumulate the trajectory points into the existing list Rt.

> for i from 1 to 1000 do:

> t0:=(i-1)*dt; sol:=RKn(fn,[t0,x0,px0,y0,py0,z0,pz0],6,t0+dt,dt); x0:=sol[2]: px0:=sol[3]: y0:=sol[4]: py0:=sol[5]: z0:=sol[6]: pz0:=sol[7]: Rt:=[op(Rt),[x0,y0,z0]]: od:

> PL5:=pointplot3d(Rt,axes=boxed,color=green,labels=[x,y,z]):

> display(PL4,PL5);

> Rt[nops(Rt)];

The charged particle spiralled out again, but in a different direction.

We observe qualitatively some spiralling of the charged-particle orbits, and preferred motion in the z direction towards weaker fields. This can be used for trapping of charged particles (magnetic bottle, confinement in fusion reactors), which is discussed further below.

A number of open questions remains:

1) we have simply followed the mathematical convention for the traversal of the current loop. What happens to the magnetic field components if the current direction (direction in which the curve integral is carried out) is reversed?

2) we have considered positively charged particles; what happens if the sign of the Lorentz force is changed?

3) most importantly: we have been sloppy as far as the units are concerned, i.e., we are just looking qualitatively at the trajectories. For real-life applications we would need to understand better what our units (m=1, e=1) imply, i.e., how do we convert to SI units?

>

>

Simple magnetic bottle

We now proceed with the model of a magnetic bottle based on four current loops. First we restart and define the analytic expressions for the magnetic field of a single loop.

> restart; with(linalg): with(plots):

Warning, new definition for norm

Warning, new definition for trace

We pick for the Helmholtz coil radius:

> R_Hh:=0.005;

> Qp:=(Q+R)^2+a^2: Qm:=(Q-R)^2+a^2:

> Qr:=2*sqrt(Q*R/Qp);

> Bxex:=i_L*a/(5*10^(6))/(sqrt(Qp)*Q)*((R^2+Q^2+a^2)*EllipticE(Qr)/Qm-EllipticK(Qr));

> Bxa:=unapply(subs(i_L=1000,R=R_Hh,Bxex),Q,a):

> Bzex:=i_L/(5*10^(6))/sqrt(Qp)*(EllipticE(Qr)*(R^2-Q^2-a^2)/Qm+EllipticK(Qr));

> Bza:=unapply(subs(i_L=1000,R=R_Hh,Bzex),Q,a):

We define a second function for the loops with smaller radius, which are introduced to confine the charged particles. The idea is to create a higher density of magnetic field lines inside the Helmholtz coils. For the confining field radius we choose:

> R_cf:=0.002;

We choose twice the current in the confining coils to obtain a large increase in the magnetic field at the bottlenecks.

> Bxb:=unapply(subs(i_L=2000,R=R_cf,Bxex),Q,a):

> Bzb:=unapply(subs(i_L=2000,R=R_cf,Bzex),Q,a):

To avoid the problems associated with taking the limit at x =0, z =0, we add a tiny constant to the cylindrical coordinate r .

> eta:=10^(-Digits):

> Bfield:=proc(x,y,z) local r,Btr,Bx,By,Bz,cphi,sphi; global Bxa,Bza,eta;

> r:=sqrt(x^2+y^2)+eta;

> Bz:= evalf(Bza(z+0.5*R_Hh,r)+Bza(z-0.5*R_Hh,r)+Bzb(z+0.5*R_Hh,r)+Bzb(z-0.5*R_Hh,r)); Btr:=evalf(Bxa(z+0.5*R_Hh,r)+Bxa(z-0.5*R_Hh,r)+Bxb(z+0.5*R_Hh,r)+Bxb(z-0.5*R_Hh,r));

> cphi:=evalf(x/r); sphi:=evalf(y/r);

> By:=Btr*sphi; Bx:=Btr*cphi;

> [Bx,By,Bz]; end:

> Bfield(0.002,0.003,0.004);

> Bfield(0.0,0.00,0.004);

For fast calculations and to be able to work in Maple5 we introduce our own Runge-Kutta procedure used before for the fieldlines rather than relying on dsolve[numeric].

> RKn:=proc(f,ic::list,N,x1,h) local xn,yn,ynew,m1,m2,m3,m4,ii,nn;

> xn:=op(1,ic):

> yn:=[]:

> for nn from 1 to N do:

> yn:=[op(yn),op(1+nn,ic)]: od:

> for ii from 1 to round((x1-xn)/h) do:

> m1:=[]:

> for nn from 1 to N do:

> m1:=[op(m1),evalf(f(xn,yn)[nn])]; od:

> m2:=[]:

> for nn from 1 to N do:

> m2:=[op(m2),evalf(f(xn+h/2,evalm(yn+h/2*m1))[nn])]; od:

> m3:=[]:

> for nn from 1 to N do:

> m3:=[op(m3),evalf(f(xn+h/2,evalm(yn+h/2*m2))[nn])]; od:

> xn:=xn+h;

> m4:=[]:

> for nn from 1 to N do:

> m4:=[op(m4),evalf(f(xn,evalm(yn+h*m3))[nn])]; od:

> ynew:=[]:

> for nn from 1 to N do:

> ynew:=[op(ynew),evalf(yn[nn]+h/6*(m1[nn]+2*(m2[nn]+m3[nn])+m4[nn]))]; od: yn:=ynew;

> od:

> [xn,op(ynew)];

> end:

We define a RHS procedure for the problem.

> fn:=proc(t,yn) local x,y,z,p_x,p_y,p_z,BF; global Bfield;

> x:=yn[1]: p_x:=yn[2]: y:=yn[3]: p_y:=yn[4]: z:=yn[5]: p_z:=yn[6]: BF:=Bfield(x,y,z);

> [p_x,p_y*BF[3]-p_z*BF[2],p_y,p_z*BF[1]-p_x*BF[3],p_z,p_x*BF[2]-p_y*BF[1]];

> end:

The argument list contains the external RHS, a list that starts with the initial time (independent variable), and the initial values, the number of dependent variables, the final value of the time (independent variable), and the time-step (step-size of the independent variable).

It seems that we can use a large step dt=0.5. The reader should check occasionally whether a reduction in the time step results in reproduced trajectories. This is the price to be paid for using a non-adaptive algorithm.

We prepare for a graph of the current loop which causes the magnetic field and the trajectory:

> PL1:=spacecurve([R_Hh*cos(f),R_Hh*sin(f),0.5*R_Hh],f=0..2*Pi,color=blue):

> PL2:=spacecurve([R_Hh*cos(f),R_Hh*sin(f),-0.5*R_Hh],f=0..2*Pi,color=blue):

> PL3:=spacecurve([R_cf*cos(f),R_cf*sin(f),0.5*R_Hh],f=0..2*Pi,color=red):

> PL4:=spacecurve([R_cf*cos(f),R_cf*sin(f),-0.5*R_Hh],f=0..2*Pi,color=red):

We pick initial conditions for a charged particle located in between the current loops with some initial momentum in the z -direction, and a small amount of motion in the x - y plane. We start the trajectory in the middle of the bottle (z=0) so that we are in a relatively weaker magnetic field region with our initial conditions.

The reader is encouraged to play with the initial conditions patiently in order to understand the outcomes.

> p_long:=-0.00004; p_trans:=0.00002;

> Rt:=[]: dt:=0.50: x0:=0.0005: px0:=-0.0: y0:=eta: py0:=p_trans: z0:=0.0: pz0:=p_long: for i from 1 to 400 do:

> t0:=(i-1)*dt; sol:=RKn(fn,[t0,x0,px0,y0,py0,z0,pz0],6,t0+dt,dt); x0:=sol[2]: px0:=sol[3]: y0:=sol[4]: py0:=sol[5]: z0:=sol[6]: pz0:=sol[7]: Rt:=[op(Rt),[x0,y0,z0]]: od:

> PL5:=pointplot3d(Rt,axes=boxed,color=green,labels=[x,y,z]):

> display(PL1,PL2,PL3,PL4,PL5);

p_long=-0.00004; p_trans=0.00002 work. Too much longitudinal momentum: violates the condition below. Too much transverse momentum gets the charged particle out of the strong-field region sideways.

J.D. Jackson in Classical Electrodynamics derives an expression for the trapping condition. The ratio of the longitudinal to transverse initial velocity of the particle has to be bounded from above by an expression that involves the ratio of magnetic field strengths at maximum as compared to the middle region.

v_L(0)/v_T(0) < sqrt(B_max/B-1)

The above example shows that the trajectories can be complicated. We can check whether the confinement result is consistent with the above formula. The reader should change the initial conditions to see that not all of them lead to particle trapping.

> abs(p_long/p_trans);

> B0:=norm(Bfield(0.,0.,0.));

> B_m:=norm(Bfield(eta,eta,0.0024));

> sqrt(B_m/B0-1);

We can make the particle escape by exceeding this ratio substantially:

> p_long:=-0.00008; p_trans:=0.00002;

> Rt:=[]: dt:=0.50: x0:=0.0000: px0:=-0.0: y0:=eta: py0:=p_trans: z0:=0.0: pz0:=p_long: for i from 1 to 400 do:

> t0:=(i-1)*dt; sol:=RKn(fn,[t0,x0,px0,y0,py0,z0,pz0],6,t0+dt,dt); x0:=sol[2]: px0:=sol[3]: y0:=sol[4]: py0:=sol[5]: z0:=sol[6]: pz0:=sol[7]: Rt:=[op(Rt),[x0,y0,z0]]: od:

> PL5:=pointplot3d(Rt,axes=boxed,color=green,labels=[x,y,z]):

> display(PL1,PL2,PL3,PL4,PL5);

>

>

>