Worksheet 9: Motion of Charged Particles in Electromagnetic Fields

In this worksheet we continue exploring electromagnetism using Mathematica and Fortran. We will now study the motion of charged particles in crossed electric and magnetic fields. This will necessitate solving ordinary differential equations (ODE's). This will be our last worksheet for this semester. We hope you enjoyed learning more about Mathematica and being introduced to Fortran. We have by no means covered either one exhaustively, and hopefully you will continue learning more about computational physics throughout your education and career. You will have three class periods for this worksheet.

Problem Formulation

Consider the motion of a particle with charge q with initial velocity $% \vec{v}_{0}$ in an electromagnetic field with electric field $\vec{E}(\vec{r}% ,t)$ and magnetic field $\vec{B}(\vec{r},t)$. In general, these can be any function of coordinates and time, but in this worksheet we will restrict our attention to the case when the direction of the fields is fixed, and their magnitudes oscillate harmonically with a frequency $\omega$ and amplitudes $% \delta E$ and $\delta B$ around E₀ and B₀. If we let the angle between $\vec{E}$ and $\vec{B}$ be $\theta$ and the phase displacement between their oscillation $\varphi$, then it is quite general to assume that:

 

$$\vec{B} B_{0}\left( 1+\delta B\cdot \cos \omega t\right) \cdot \hat{z}$$ = (1)

$$\vec{E} E_{0}\left[ 1+\delta E\cdot \cos \left( \omega t+\varphi \right) \right] \cdot \left( \hat{x}\sin \theta +\hat{z}\cos \theta \right)$$ =

In the worksheet we will consider only specific simpler cases of this general case. The classical motion of the particle in such a field is determined by Newton's equation of motion, where the force is given by:

$$\vec{a}=\frac{\vec{F}}{m}=\frac{q}{m}(\vec{E}+\vec{v}\times \vec{B})$$ (2)

As usual, we will be interested mostly in qualitative understanding of the motion, so set $\frac{q}{m}=1$. We can now write this as an ODE for $\vec{r}=% \vec{r}(t)$ with a given initial velocity $\vec{v}(t=0)=\vec{v}_{0}$ and starting position $\vec{r}(t=0)=\vec{r}_{0}$:

 

$$\frac{d^{2}\vec{r}(t)}{dt^{2}} \frac{q}{m}\left[ \vec{E}+\frac{d\vec{r}(t)% }{dt}\times \vec{B}\right]$$ = (3)

$$\vec{r}(0) \vec{r}_{0}\text{, }\frac{d\vec{r}(0)}{dt}=\vec{v}_{0}$$ =

More formally, equation is what is known as a second order initial value problem (IVP), and there are many numerical techniques aimed at solving this type of extremely common equation. Mathematica 's function NDSolve is equipped with a suite of such techniques, and there many Fortran libraries that solve problems of the type as well.

Numerical Solution of IVP's: Euler Method

First Order IVP's

When approaching a problem numerically, it is usually wise to write in the most general (yet specific) way, so that the same numerical algorithm can be used for a variety of problems. A typical second order IVP has the form:

 

$$\frac{d^{2}\vec{y}}{dx^{2}} \vec{f}(x,\vec{y},\frac{d\vec{y}}{dx})$$ = (4)

$$\vec{y}(x_{o}) \vec{y}_{0}\text{, }\frac{d\vec{y}(x_{o})}{dx}=\vec{y}% _{0}^{\prime }$$ =

Look at equation and identify all the variables in the above equation before proceeding. Most numerical methods are tailored for the first order IVP,

 

$$\frac{d\vec{y}}{dx} \vec{f}(x,\vec{y})$$ = (5)

$$\vec{y}(x_{o}) \vec{y}_{0}$$ =

since any problem of order n, such as (n=2), can be cast as a first order problem by introducing auxiliary variables for the n-1 lower derivatives. For example, eq. can be written in the form of eq. by introducing a new variable for the first derivative $\vec{y}% ^{\prime }$ (velocity):

 

$$\frac{d\vec{y}^{\prime }}{dx} \vec{f}(x,\vec{y},\vec{y}^{\prime })\text{, }\vec{y}^{\prime }(x_{o})=\vec{y}_{0}^{\prime }$$ = (6)

$$\frac{d\vec{y}}{dx} \vec{y}^{\prime }\text{, }\vec{y}(x_{o})=\vec{y}_{0}% \text{ }$$ =

Although this is formally system of first order IVP's, it is in fact equivalent to a formal first order IVP for the variable that appends the first derivative with the variable itself, $\vec{Y}=[\vec{y}^{\prime }\vdots \vec{y}]$,

 

$$\frac{d\vec{Y}}{dx} \vec{F}(x,\vec{Y})\text{, }\vec{Y}(x_{o})=[\vec{y}% _{0}^{\prime }\vdots \vec{y}_{0}]$$ = (7)

$$\text{where }\vec{F}(x,\vec{Y}) \left[ \vec{f}(x,\vec{y},\vec{y}^{\prime })\vdots \vec{y}^{\prime }\right]$$ =

so that any method that can solve can also solve and thus also . Before proceeding, rewrite equation as a first order IVP (either in the form or ) and show your result to your teacher or TA.

Euler's Method

From here on we will only discuss problem , having in mind the previous discussion. The simplest numerical method for solving such an equation is Euler's method, which steps through the independent variable (time) in small time steps dx from the current value x, using a first-order Taylor expansion to approximate the dependent variable y at the next time step x+dx:

$$\vec{y}(x+dx)\approx \vec{y}(x)+\frac{d\vec{y}(x)}{dx}=\vec{y}(x)+\vec{f}% \left[ x,\vec{y}(x)\right]$$

It is clear that once we have an equation for $\vec{f}$ we can start from the initial point x₀ and step through time until any desired point to obtain an approximate solution for the IVP. Euler method is very slow and a very small step size dx is needed to achieve good results, but it is prototypical of the commonly used methods, such as Runge-Kutta methods (see Numerical Recipies for a description, or the file /classes/phy201/RK4.f90).