SULI Computational Workshop¶

Welcome to the computational workshop at the 2022 Introduction to Fusion Energy and Plasma Physics Course! The goal of this workshop is to give an introduction to PlasmaPy: an open source community-developed Python package for plasma research and education.

Today we'll start by going through astropy.units, which is used quite heavily in PlasmaPy, SunPy, and many other packages in the scientific pythoniverse. We'll then take a tour of plasmapy.particles and end with plasmapy.formulary. We will be filling out the notebook interactively during this session.

If you get stuck at any point, please feel free to ask in the Zoom chat, on Discord, or in PlasmaPy's Element chat. You can also refer to the the completed notebook.

To execute a cell in a Jupyter notebook, press Shift + Enter.

Let's start with some preliminaries¶

Let's execute a few cells to get us ready for the tutorial. If you end up having to restart the notebook, please run the cells in this section again.

If you are on Google Colab, run the next cell with Shift + Enter to install PlasmaPy. When the warning comes up, click Run anyway. This may take ∼30 seconds.

You should not need to restart the runtime at this point, and some warning messages might show up.

Run the next cell with Shift + Enter to perform some initial imports and settings.

If the above cell doesn't work, try running a new cell with the following command to install PlasmaPy:

%pip install plasmapy

Then run the above cell again.

If you are running this on Binder, remember that the notebook will time out after 10 minutes of inactivity.

Astropy units¶

In scientific computing, we often represent physical quantities as numbers.

Representing a physical quantity as a number has risks. We might unknowingly perform operations with different units, like time_in_seconds + time_in_hours. We might even accidentally perform operations with physically incompatible units, like length + time, without catching our mistake. We can avoid these problems (and potentially $327.6M errors) by using a units package.

This notebook introduces astropy.units with an emphasis on the functionality needed to work with plasmapy.particles and plasmapy.formulary. We typically import this subpackage as u.

Unit essentials¶

We can create a physical quantity by multiplying or dividing a number or array with a unit.

This operation creates a Quantity: a number, sequence, or array that has been assigned a physical unit.

We can also create an object by using the Quantity class itself.

We can create Quantity objects with compound units.

We can even create Quantity objects that are explicitly dimensionless.

We can also create a Quantity based off of a NumPy array or a list.

Unit operations¶

Operations between Quantity objects handle unit conversions automatically. We can add Quantity objects together as long as their units have the same physical type.

Units get handled automatically during operations like multiplication, division, and exponentiation.

Attempting an operation between physically incompatible units gives us an error, which we can use to find bugs in our code.

Quantity objects behave very similarly to NumPy arrays because Quantity is a subclass of numpy.ndarray.

Most frequently encountered NumPy and SciPy functions can be used with Quantity objects. However, Quantity objects lose their units with some operations.

Unit conversions¶

The to() method allows us to convert a Quantity to different units of the same physical type. This method accepts strings that represent a unit (including compound units) or a unit object.

The si and cgs attributes convert the Quantity to SI or CGS units, respectively.

Detaching units and values¶

The value attribute of a Quantity provides the number (as a NumPy scalar) or NumPy array without the unit.

The unit attribute of a Quantity provides the unit without the value.

Equivalencies¶

Plasma scientists often use the electron-volt (eV) as a unit of temperature. This is a shortcut for describing the thermal energy per particle, or more accurately the temperature multiplied by the Boltzmann constant, $k_B$. Because an electron-volt is a unit of energy rather than temperature, we cannot directly convert electron-volts to kelvin.

To handle non-standard unit conversions, astropy.units allows the use of equivalencies. The conversion from eV to K can be done by using the temperature_energy() equivalency.

Radians are treated dimensionlessly when the dimensionless_angles() equivalency is in effect. Note that this equivalency does not account for the multiplicative factor of $2π$ that is used when converting between frequency and angular frequency.

Physical constants¶

We can use astropy.constants to access the most commonly needed physical constants.

A Constant behaves very similarly to a Quantity. For example, we can use the Boltzmann constant to mimic the behavior of u.temperature_energy().

Electromagnetic constants often need the unit system to be specified. Code within PlasmaPy uses SI units.

Physical types¶

A physical type corresponds to physical quantities with dimensionally compatible units. Astropy has functionality that represents different physical types. These physical type objects can be accessed using either the physical_type attribute of a unit or get_physical_type().

These physical type objects can be used for dimensional analysis.

Particles¶

The plasmapy.particles subpackage contains functions to access basic particle data, and classes to represent particles.

Particle properties¶

There are several functions that provide information about different particles that might be present in a plasma. The input of these functions is a representation of a particle, such as a string for the atomic symbol or the element name.

We can provide a number that represents the atomic number.

We can also provide standard symbols or the names of particles.

The symbols for many particles can even be used directly, such as for an alpha particle. To create an "α" in a Jupyter notebook, type \alpha and press tab.

We can represent isotopes with the atomic symbol followed by a hyphen and the mass number. Let's use half_life to return the half-life of a radioactive particle in seconds as a Quantity.

We typically represent an ion in a string by putting together the atomic symbol or isotope symbol, a space, the charge number, and the sign of the charge.

Functions in plasmapy.particles are quite flexible in terms of string inputs representing particles. An input is particle-like if it can be used to represent a physical particle.

Most of these functions take additional arguments, with Z representing the charge number of an ion and mass_numb representing the mass number of an isotope. These arguments are often keyword-only to avoid ambiguity.

Particle objects¶

Up until now, we have been using functions that accept representations of particles and then return particle properties. With the Particle class, we can create objects that represent physical particles.

Particle properties can be accessed via attributes of the Particle class.

Antiparticles¶

We can get antiparticles of fundamental particles by using the antiparticle attribute of a Particle.

We can also use the tilde (~) operator on a Particle to get its antiparticle.

Ionization and recombination¶

The recombine() and ionize() methods of a Particle representing an ion or neutral atom will return a different Particle with fewer or more electrons.

When provided with a number, these methods tell how many bound electrons to add or remove.

If the inplace keyword is set to True, then the Particle will be replaced with the new particle.

Custom particles¶

Sometimes we want to use a particle with custom properties. For example, we might want to represent an average ion in a multi-species plasma. For that we can use CustomParticle.

Many of the attributes of CustomParticle are the same as in Particle.

If we do not include one of the physical quantities, it gets set to numpy.nan (not a number) in the appropriate units.

CustomParticle objects are not yet able to be used by many of the functions in plasmapy.formulary, but are expected to become compatible with them in PlasmaPy 0.9.0.

Particle lists¶

The ParticleList class is a container for Particle and CustomParticle objects.

By using a ParticleList, we can access the properties of multiple particles at once.

We can also create a ParticleList by adding Particle and/or CustomParticle objects together.

We can also get an average particle.

ParticleList objects are also expected to become compatible with many functions in plasmapy.formulary by PlasmaPy's 0.9.0 release.

Particle categorization¶

The categories attribute of a Particle provides a set of the categories that the Particle belongs to.

The is_category() method lets us determine if a Particle belongs to one or more categories.

If we need to be more specific, we can use the require keyword for categories that a Particle must belong to, the exclude keyword for categories that the Particle cannot belong to, and the any_of keyword for categories of which a Particle needs to belong to at least one.

The valid_categories attribute of is_category() for any Particle gives a set containing all valid categories.

In the next release of PlasmaPy, there will be an is_category() method on ParticleList too.

Creating functions with particle inputs¶

It's often the case that we'll want to write functions that depend on particles. For this need, we can use @particle_input. This decorator processes arguments that have an annotation of Particle. Let's write a function f to try this out.

Let's pass this function a string that represents a proton.

Here, "proton" was transformed by @particle_input into a Particle before being passed into the function.

Nuclear reactions¶

We can use plasmapy.particles to calculate the energy of a nuclear reaction using the > operator.

If the nuclear reaction is invalid, then an exception is raised that states the reason why.

Ionization state data structures¶

The ionization state distribution for an element refers to the fractions of that element at each ionic level. For example, the charge state of helium might be 10% He$^{0+}$, 70% He$^{1+}$, and 20% He$^{2+}$. Let's store this information in an IonizationState object, and specify the number density of the element too.

The ionization state distribution is stored in the ionic_fractions attribute.

Because we provided the number density of the element as a whole, we can get back the number density of each ionic level.

We can also get the electron density under the assumption of quasineutrality.

We can use the summarize() method to get information about the ionization state.

Now let's look at some actual ionization state data for filament material in an interplanetary coronal mass ejection (CME) observed by the Advanced Composition Explorer (ACE) near 1 AU, averaged over about an hour. The data were extracted from Figure 4 in Gilbert et al. (2012). These data are noteworthy because there is information from very low charge states to very high charge states.

Now let's use this information as an input for IonizationStateCollection: a data structure for the ionization states of multiple elements.

We can summarize this information too, but let's specify the minimum ionic fraction to print.

Now let's visualize the ionization state for each of these elements, remembering that these values were averaged over about an hour.

Question: Does this plasma look like it all came from the same temperature?

PlasmaPy formulary¶

The plasmapy.formulary subpackage contains a broad variety of formulas needed by plasma scientists across disciplines, in particular to calculate plasma parameters.

Plasma beta in the solar atmosphere¶

Plasma beta ($β$) is one of the most fundamental plasma parameters. $β$ is the ratio of the thermal plasma pressure to the magnetic pressure:

How a plasma behaves depends strongly on $β$. When $β ≫ 1$, the magnetic field is not strong enough to significantly alter the dynamics, so the plasma motion is more gas-like. When $β ≪ 1$, magnetic tension and magnetic pressure dominate the dynamics.

Let's use plasmapy.formulary.beta to compare $β$ in different parts of the solar atmosphere.

Solar corona¶

Let's start by defining some plasma parameters for an active region in the solar corona.

When we use these parameters in beta, we find that $β$ is quite small which implies that the corona is magnetically dominated.

Solar chromosphere¶

Next let's calculate $β$ for the chromosphere. Bogod et al. (2015) found that the quiet chromosphere ranges from ∼40–200 G. We can get the temperature and number density of hydrogen from model C7 of Avrett & Loeser (2007) for 1 Mm above the photosphere.

When $B$ is small, plasma $β$ is not too far from 1, which means that both magnetic and plasma pressure gradient forces are important when the chromospheric magnetic field is relatively weak. When near the higher range of $B$, $β$ is small so that the magnetic forces are more important than plasma pressure gradient forces.

Quiet solar photosphere¶

Let's specify some characteristic plasma parameters for the solar photosphere.

When we calculate β for the photosphere, we find that it is an order of magnitude larger than 1, so plasma pressure gradient forces are more important than magnetic forces.

Sunspot photosphere¶

The magnetic field in the solar photosphere is strongest in sunspots, so we would expect β to be lowest there. Let's estimate some plasma parameters for a sunspot.

When we calculate β, we find that both pressure gradient and magnetic forces will be important.

Plasma parameters in Earth's magnetosphere¶

The Magnetospheric Multiscale Mission (MMS) is a constellation of four identical spacecraft. The goal of MMS is to investigate the small-scale physics of magnetic reconnection in Earth's magnetosphere. In order to do this, the spacecraft need to orbit in a tight configuration. But how tight does the tetrahedron have to be? Let's use plasmapy.formulary to find out.

Physics background¶

Magnetic reconnection is the fundamental plasma process that converts stored magnetic energy into kinetic energy, thermal energy, and particle acceleration. Reconnection powers solar flares and is a key component of geomagnetic storms in Earth's magnetosphere. Reconnection can also degrade confinement in fusion devices such as tokamaks.

In the classic Sweet-Parker model, reconnection occurs when oppositely directed magnetic fields are pressed towards each other in a plasma in an elongated current sheet where they can resistively diffuse. The reconnection rate is slow because of the bottleneck associated with conservation of mass.

The inertial length for a particle species is the characteristic length scale for getting accelerated or decelerated by forces in a plasma.

When the reconnection layer thickness is shorter than the ion inertial length, $d_i ≡ c/ω_{pi}$, collisionless effects and the Hall effect enable reconnection to be fast (Zweibel & Yamada 2009). The inner electron diffusion region has a thickness of about the electron inertial length, $d_i≡c/ω_{pe}$.

Our goal: calculate $d_i$ and $d_e$ to get an idea of how far the MMS spacecraft should be separated from each other to investigate reconnection.

Length scales¶

Let's choose some characteristic plasma parameters for the magnetosphere.

Let's calculate the ion inertial length, $d_i$. On length scales shorter than $d_i$, the Hall effect becomes important as the ions and electrons decouple from each other.

The ion diffusion regions should therefore be a few hundred kilometers thick. Let's calculate the electron inertial length next.

The electron diffusion region should therefore have a characteristic length scale of a few kilometers, which is significantly smaller than the ion diffusion region.

We can also calculate the gyroradii for different particles

The four MMS spacecraft have separations of ten to hundreds of kilometers, and thus are well-positioned to investigate Hall physics during reconnection in the magnetosphere.

Frequencies¶

We can also calculate some of the fundamental frequencies associated with magnetospheric plasma.

1D Maxwellian distribution function¶

Let's try out Maxwellian_1D(), also from plasmapy.formulary. As a first example, let's get the probability density of finding an electron with a speed of 500 km/s if we have a plasma at a temperature of $3×10^4$ K:

Let's now plot the 1D Maxwellian distribution function for protons and alpha particles.

Thermal Speed¶

The thermal_speed() function can be used to calculate the thermal velocity for a Maxwellian velocity distribution. There are three common definitions of the thermal velocity, which can be selected using the method keyword, which are defined for a 3D velocity distribution as

• "most_probable"
  

$v_{th} = \sqrt{\frac{2 k_B T}{m}}$

• "rms"
  

$v_{th} = \sqrt{\frac{3 k_B T}{m}}$

• "mean_magnitude"
  

$v_{th} = \sqrt{\frac{8 k_B T}{m\pi}}$

The differences between these velocities can be seen by plotting them on a 3D Maxwellian speed distribution:

Emission of Thermal Bremsstrahlung by a Maxwellian Plasma¶

Bremsstrahlung is electromagnetic radiation produced by the deceleration of a charged particle when deflected by another charged particle.

The thermal_bremsstrahlung() function calculates the bremsstrahlung spectrum emitted by the collision of electrons and ions in a thermal (Maxwellian) plasma. This function calculates this quantity in the Rayleigh-Jeans limit where $\hbar\omega \ll k_B T_e$.

Let's start by setting the plasma density, temperature, and ion species.

Let's create an array of frequencies over which to calculate the bremsstrahlung spectrum and convert these frequencies to photon energies for the purpose of plotting the results.

Let's calculate the spectrum and then plot it.

Coulomb logarithms¶

Coulomb collisions are collisions between two charged particles where the interaction is governed solely by the electric fields from the two particles. Coulomb collisions usually result in small deflections of particle trajectories. The deflection angle depends on the impact parameter of the collision, which is the perpendicular distance between the particle's trajectory and the particle it is colliding with. High impact parameter collisions (which result in small deflection angles) occur much more frequently than low impact parameter collisions (which result in large deflection angles).

The minimum and maximum impact parameters ($b_\min$ and $b_\max$, respectively) represent the range of distances of closest approach. While a typical Coulomb collision results in only a slight change in trajectory, the effects of these collisions are cumulative, and it is necessary to integrate over the range of impact parameters to account for the effects of Coulomb collisions throughout the plasma. The Coulomb logarithm accounts for the range in impact parameters for the different collisions, and is given by

But what should we use for the impact parameters?

Usually $b_\max$ is given by the Debye length, $λ_D$, which can be calculated with Debye_length(). On length scales $≳ λ_D$, electric fields from individual particles get cancelled out due to screening effects. Consequently, Coulomb collisions with impact parameters $≳ λ_D$ will rarely occur.

The inner impact parameter $b_\min$ requires more nuance. One possibility would be to set $b_\min$ to be the impact parameter corresponding to a 90° deflection angle, $ρ_⟂$, which can be calculated with impact_parameter_perp(). Alternatively, $b_\min$ could be set to be the de Broglie wavelength, $λ_{dB}$, calculated using the reduced mass, $μ$, of the particles. Typically,

The impact_parameter() function in plasmapy.formulary simultaneously calculates both $b_\min$ and $b_\max$. Let's estimate $b_\min$ and $b_\max$ for proton-electron collisions in the solar corona.

When we can calculate the Coulomb logarithm, we find that it is ∼20 (a common value for astrophysical plasma).

Our next goals are visualize the Coulomb logarithm and impact parameters for different physical conditions. Let's start by creating a function to plot $\ln{Λ}$ against temperature.

Next we will create a function to plot $b_\min$ and $b_\max$ as functions of temperature, along with the $ρ_⟂$ and $λ_{dB}$.

Let's now plot the Coulomb logarithm and minimum/maximum impact parameters over a range of temperatures in a dense plasma.

To investigate what is happening at the sudden change in slope, let's plot the impact parameters against temperature.

At low temperatures, $b_\min$ is $ρ_⟂$. At higher temperatures, $b_\min$ is $λ_{dB}$. What happens if we look at lower temperatures?

The Coulomb logarithm becomes negative! Let's look at the impact parameters again to understand what's happening.

This unphysical situation occurs because $b_\min > b_\max$ at low temperatures.

So how should we handle this? Fortunately, PlasmaPy's implementation of Coulomb_logarithm includes the methods described by Gericke, Murillo, and Schlanges (2002) for dense, strongly-coupled plasmas. For most cases, we recommend using method="hls_full_interp" in which $b_\min$ is interpolated between $λ_{dB}$ and $ρ_⟂$, and $b_\max$ is interpolated between $λ_D$ and the ion sphere radius:

The Coulomb logarithm approaches zero as the temperature decreases, and does not become unphysically negative. This is the expected behavior, although there is increased uncertainty for Coulomb logarithms $≲4$.

And that's it for PlasmaPy! Thank for you staying (or at least scrolling) to the end!