Motion of Trapped Protons in the Earth’s Magnetic Field

The Earth has a substantial magnetic field that extends outward for many thousands of kilometers. It is thought that this magnetic field is generated by circulating currents within a spinning metal liquid core — the dynamo effect. This magnetic field closely resembles a dipole field; however, there is a tilt between the spin axis of the Earth and that of the magnetic dipole. There are also other asymmetries in the Earth’s magnetic field that require the use of software models to describe fully. These models are constructed using data obtained from orbiting spacecraft.

Having an accurate model of the Earth’s magnetic field is very important in studies of the Earth’s deep interior, its crust, and its ionosphere and magnetosphere. Therefore, a standard model is maintained by the International Association of Geomagnetism and Aeronomy (IAGA). This model is referred to as the International Geomagnetic Reference Field (IGRF). The IGRF model is updated every few years and is currently on its 12th generation (Ref. 1).

The influence of the solar wind means that the full extended shape of the Earth’s magnetosphere is very different from a dipole. However, within a few Earth radii, simpler models of the Earth’s magnetic field are sufficient.

Charged particles in the Earth’s magnetic field travel in helical paths around the field lines. The angle between the direction of the magnetic field and a particle's trajectory is referred to as the pitch angle. As particles move along the magnetic field lines, they approach a pole in the magnetic field. Near the poles, the particles get closer to earth, causing the magnetic field to increase in magnitude and thus causing the pitch angle to increase. If the pitch angle reaches 90 degrees while the particle is still outside of the atmosphere (>100 km), the particle motion reverses direction back along the field line. The point at which this occurs is referred to as a mirror, or bounce, point. The particle similarly bounces off the other mirror point in the other hemisphere.

Particles are considered trapped if they are not lost to the Earth’s atmosphere during this motion. Particles can also be liberated due to scattering by electromagnetic waves.

The latitude of a mirror point (λm) is related to the particle’s pitch angle at the equatorial plane, or its equatorial pitch angle (αeq). As the equatorial pitch angle increases, the mirror point latitude decreases. Any particles that have equatorial pitch angles that give mirror points outside the atmosphere are said to be outside the atmospheric loss cone and are trapped.

Particles that bounce from mirror point to mirror point also exhibit a drift motion around the Earth, switching from field line to field line. This is due to the fact that the magnetic field magnitude is increasing as the particle moves toward the Earth, so that its gyration is not circular but has in fact a smaller radius on the side closer to the Earth.

As electrons and protons drift in opposite directions around the Earth, an electric current (ring current) is set up. The magnitude of the ring current increases during solar storms, and its effect can be measured on the ground as a weakening of the measured magnetic field. The disturbance storm time index, Dst, is one such measure of this ring current and is used to assess the severity of magnetic storms from the Sun.

These motions of trapped charge particles lead to “belts” of energetic charged particles in the near Earth environment and are referred to as the Van Allen radiation belts. They extend from about 1000 to 60,000 km (~0.15 to ~10 Earth radii) above the Earth’s surface and therefore pose a real threat to space-based microelectronics.

Model Definition

Mathematically, the IGRF model consists of the Gauss coefficients, which define a spherical harmonic expansion of the magnetic scalar potential:

where r is the radial distance from the Earth’s center, L is the maximum degree of the expansion, Φ is East longitude, θ is colatitude (polar angle), Re is the Earth’s radius, glm and hlm are Gauss coefficients, and Plmcos θ are the Schmidt normalized associated Legendre functions of degree l and order m.

The model uses a simple sphere of radius Re to represent the Earth within a larger spherical simulation domain of radius 5Re where this particle trajectories are computed. The geometry is shown in Figure 1.

The Magnetic Force feature in the Charged Particle Tracing interface and the Particle Tracing for Fluid Flow interface includes a built-in option to compute magnetic fields using the data from the IGRF model. To access this data, in a 3D model select from the list in the settings window for the Magnetic Force feature.

Figure 2 shows the IGRF magnetic field lines. The difference from a simple dipole is not really evident in this figure, but asymmetries exist in both latitude and longitude.

Figure 1: Geometry of the simulation domain, which extends from radius Re to 5Re.

Figure 2: IGRF magnetic field lines.

The model uses the Charged Particle Tracing interface only, together with a Time Dependent study.

Results and Discussion

A single 10 MeV proton is released from the equatorial plane at a distance of 2Re from the center of the Earth. It is released with an equatorial pitch angle of 30 degrees. The three components of its motion (gyration, bounce, and drift) are clearly visible in Figure 3. The timescales for the drift motion are much longer than that of the bounce motion, which in turn is much longer than the gyration period.

Figure 3: Particle trajectory of a single 10 MeV proton in Earth’s magnetic field. Particle gyration, bounce and drift motion are clearly visible.

Due to the large changes in the magnitude of the magnetic field in different parts of the simulation domain, care has to be taken to ensure an appropriate time step is used for the particles. An good indication of how well the time steps are resolving the gyration of the particles is to plot the particle energy as a function of time and observe little or no change. Figure 4 plots the relative change in the particle kinetic energy as a function of time. This relative error is of the order 10−3 and so you can be certain you are resolving most of the particle motion.

Figure 4: Relative error in particle kinetic energy.

There are many aspects of the particle motion that can be investigated using this model, but the scope of this tutorial is limited to investigating how the mirror point latitude is affected by the equatorial pitch angle of the released protons.

As mentioned above, as the equatorial pitch angle of a particle decreases, the mirror point latitude is expected to increase. At the limits, a particle with an equatorial pitch angle of 90 degrees would remain in the equatorial plane whereas one with a pitch angle of zero degrees would travel directly along the field line without bouncing.

Figure 6 shows the mirror point latitude versus equatorial pitch angle for 15 particles with equatorial pitch angles ranging from 5 degrees to 80 degrees.

Figure 5: Trajectories of several trapped protons. The color expression corresponds to the equatorial pitch angle of each proton.

Figure 6: Mirror point latitude (λm) against particle equatorial pitch angle (Ea).

References

1. International Geomagnetic Reference Field website, http://www.ngdc.noaa.gov/IAGA/vmod/igrf.html.

Particle_Tracing_Module/Charged_Particle_Tracing/trapped_protons

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select .

Click .

Click

In the tree, select .

Click

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Value	Description
Re	6371.2[km]	6.3712E6 m	Radius of the Earth
E0	10[MeV]	1.6022E-12 J	Initial particle energy
alpha	30[deg]	0.5236 rad	Equatorial pitch angle

Geometry 1

Add a sphere with radius Re within a larger simulation domain of radius 5*Re.

In the window, under click .

In the window for , locate the section.

From the list, choose .

Sphere 1 (sph1)

In the toolbar, click

In the window for , locate the section.

In the text field, type Re.

Sphere 2 (sph2)

In the toolbar, click

In the window for , locate the section.

In the text field, type 5*Re.

Difference 1 (dif1)

In the toolbar, click

and choose .

Select the object only.

In the window for , locate the section.

Find the subsection. Select the

toggle button.

Click the

button in the toolbar.

Select the object only.

Click

Click the

button in the toolbar.

Click the

button in the toolbar.

Definitions

Hide the outermost boundary to facilitate setup of the Charged Particle Tracing interface.

Hide for Geometry 1

In the window, expand the node.

Right-click and choose .

In the window for , locate the section.

From the list, choose .

On the object , select Boundaries 1–4, 9, 10, 13, and 16 only.

Ball 1

Add a selection containing all boundaries on the Earth’s surface.

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Locate the section. In the text field, type Re.

Charged Particle Tracing (cpt)

In the window, under click .

In the window for , locate the section.

From the list, choose .

In this example, the particles are relativistic protons.

Select the check box.

Particle Properties 1

In the window, under click .

In the window for , locate the section.

From the list, choose .

Magnetic Force 1

Exert a magnetic force on the particles using the Earth’s magnetic field.

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. From the B list, choose .

Release from Grid 1

Set up the feature to release a single 10 MeV proton at a distance of 2 Earth radii along the x-axis. Set its velocity x-component to zero, and its y- and z- components using the equatorial pitch angle.

In the toolbar, click

and choose .

In the window for , locate the section.

In the qx,0 text field, type 2*Re.

Locate the section. From the list, choose .

In the E0 text field, type E0.

Specify the L0 vector as


0	x
sin(alpha)	y
cos(alpha)	z

Auxiliary Dependent Variable 1

Define two auxiliary dependent variables for the mirror point latitude and equatorial pitch angle. These auxiliary dependent variables will be disabled in Study 1 but will be needed in Study 2.

In the toolbar, click

and choose .

In the window for , locate the section.

In the text field, type Lm.

Auxiliary Dependent Variable 2

In the toolbar, click

and choose .

In the window for , locate the section.

In the text field, type Ea.

Velocity Reinitialization 1

The feature can be used to detect when a particle bounces, and set the auxiliary dependent variable Lm to the instantaneous value of the latitude at that point.

In the toolbar, click