Fundamental Plasma Parameters

This tutorial shows how to calculate some fundamental parameters that are typically used to characterize the ionospheric plasma from AMISR data files. The equations used here are sourced from the NRL Plasma Formulary p28-29 [NRL Plasma Formulary, 2013] and Schunk & Nagy, 2009, p45 and p141 [Schunk and Nagy, 2009]. Note that many of the parameters shown in this tutorial are time dependent, but only single profiles are shown to demonstrate more clearely how they vary in altitude.

import h5py
import os
import urllib.request
import numpy as np
import scipy.constants as const
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import matplotlib.colors as colors

Download files that will be used in this tutorial.

# Download the file that we need to run these examples
filename = '../../data/20200207.001_lp_5min-fitcal.h5'

if not os.path.exists(filename):
    url='https://data.amisr.com/database/dbase_site_media/PFISR/Experiments/20200207.001/DataFiles/20200207.001_lp_5min-fitcal.h5'

    print('Downloading data file...')
    urllib.request.urlretrieve(url, filename)

    print('...Done!')

Note on Ions: There are multiple relevant ion species in the ionosphere. Ion mass \(m_i\) can be found by multiplying the IonMass array by the mass of a proton. The Fits Array includes the fractional ion composition, which can be used to calculate the density of indiviual ion species by multipying by the electron density (see Ion Composition for more details). These operations may require some array broadcasting and expand the size of the array to include a species dimension.

with h5py.File(filename, 'r') as h5:
    ne = h5['FittedParams/Ne'][:]
    frac = h5['FittedParams/Fits'][:,:,:,:,0]
    ion_mass = h5['FittedParams/IonMass'][:]

ni = np.expand_dims(ne,axis=-1)*frac[:,:,:,:-1]
mi = ion_mass*const.m_p
print(ne.shape, ni.shape)

for am, m in zip(ion_mass, mi):
    print(f'{am} AMU = {m} km')

mass_labels = {16:r'O$^+$', 32:r'O$_2^+$', 30:'NO$^+$', 28:'N$_2^+$', 14:'N$^+$'}

(188, 11, 74) (188, 11, 74, 5)
16.0 AMU = 2.67619508152e-26 km
32.0 AMU = 5.35239016304e-26 km
30.0 AMU = 5.0178657778499995e-26 km
28.0 AMU = 4.68334139266e-26 km
14.0 AMU = 2.34167069633e-26 km

The ion density array has the standard dimensions for record, beam, and bin, plus an additional one for the five ion species. Many of the general equations shown below include \(e_\alpha\), which is used to indicate absolute value of charge the species. For the ionospheric regimes relevant to ISR measurments, ions normally only have +1 charge states, so \(e_\alpha = e\) and \(e_\alpha\) will be implicitly replaced with just the elementary charge, \(e\), in coding examples.

# Load some parameters from the file that are only needed for plotting
with h5py.File(filename, 'r') as h5:
    beamcodes = h5['BeamCodes'][:]
    bidx = np.argmax(beamcodes[:,2])
    alt = h5['Geomag/Altitude'][:]
    utime = h5['Time/UnixTime'][:,0]
time = utime.astype('datetime64[s]')
tidx = 20

Collision Frequency

The collision frequencies are contained directly within the AMISR files.

with h5py.File(filename, 'r') as h5:
    nu_i = h5['FittedParams/Fits'][:,:,:,:-1,2]
    nu_e = h5['FittedParams/Fits'][:,:,:,-1,2]
print(nu_e.shape, nu_i.shape)

(188, 11, 74) (188, 11, 74, 5)

fig = plt.figure()
ax = fig.add_subplot(111)
for i in range(len(ion_mass)):
    ax.plot(nu_i[tidx,bidx,:,i], alt[bidx], label=mass_labels[int(ion_mass[i])])
ax.plot(nu_e[tidx,bidx,:], alt[bidx], label='e-')
ax.legend()
ax.set_xscale('log')
ax.set_xlabel('Collision Frequency (Hz)')
ax.set_ylabel('Altitude (m)')

Text(0, 0.5, 'Altitude (m)')

Debye Length

The Debye length is the minimum distance over which plasma can exhibit collective behavior.

\[ \lambda_D = \sqrt{\frac{\epsilon_0 k_B T_e}{n_e e^2}} \]

with h5py.File(filename, 'r') as h5:
    ne = h5['FittedParams/Ne'][:]
    te = h5['FittedParams/Fits'][:,:,:,-1,1]

lam_D = np.sqrt((const.epsilon_0*const.k*te)/(ne*const.e**2))
print(lam_D.shape)

(188, 11, 74)

/tmp/ipykernel_2643/3888609315.py:5: RuntimeWarning: invalid value encountered in sqrt
  lam_D = np.sqrt((const.epsilon_0*const.k*te)/(ne*const.e**2))

fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(lam_D[20,bidx,:], alt[bidx])
# ax.legend()
ax.set_xscale('log')
ax.set_xlabel('Debye Length (m)')
ax.set_ylabel('Altitude (m)')

Text(0, 0.5, 'Altitude (m)')

Thermal Velocity

\[ v_{T\alpha} = \sqrt{\frac{k_B T_\alpha}{m_\alpha}} \]

with h5py.File(filename, 'r') as h5:
    ion_mass = h5['FittedParams/IonMass'][:]
    te = h5['FittedParams/Fits'][:,:,:,-1,1]
    ti = h5['FittedParams/Fits'][:,:,:,:-1,1]

v_te = np.sqrt(const.k*te/const.m_e)
v_ti = np.sqrt(const.k*ti/(ion_mass*const.m_p))
print(v_te.shape, v_ti.shape)

(188, 11, 74) (188, 11, 74, 5)

/tmp/ipykernel_2643/3741107914.py:6: RuntimeWarning: invalid value encountered in sqrt
  v_te = np.sqrt(const.k*te/const.m_e)
/tmp/ipykernel_2643/3741107914.py:7: RuntimeWarning: invalid value encountered in sqrt
  v_ti = np.sqrt(const.k*ti/(ion_mass*const.m_p))

fig = plt.figure()
ax = fig.add_subplot(111)
for i in range(len(ion_mass)):
    ax.plot(v_ti[tidx,bidx,:,i], alt[bidx], label=mass_labels[int(ion_mass[i])])
ax.plot(v_te[tidx,bidx,:], alt[bidx], label='e-')
ax.legend()
ax.set_xscale('log')
ax.set_xlabel('Thermal Velocity (m/s)')
ax.set_ylabel('Altitude (m)')

Text(0, 0.5, 'Altitude (m)')

Gyrofrequency

The frequency at which charged particles gyrate around the magnetic field line. This is also called the cyclotron frequency.

\[ \omega_{c\alpha} = \frac{e_\alpha B}{m_\alpha} \]

Becasue there are no time-dependent parameters in the gyrofrequency, it only has a single value for each beam, bin, and species in the array (no record dimension).

with h5py.File(filename, 'r') as h5:
    B = h5['Geomag/Babs'][:]
    ion_mass = h5['FittedParams/IonMass'][:]

omega_ce = const.e*B/(const.m_e)
omega_ci = const.e*np.expand_dims(B,axis=-1)/(ion_mass*const.m_p)
print(omega_ce.shape, omega_ci.shape)

(11, 74) (11, 74, 5)

fig = plt.figure()
ax = fig.add_subplot()
for i in range(len(ion_mass)):
    ax.plot(omega_ci[bidx,:,i], alt[bidx], label=mass_labels[int(ion_mass[i])])
ax.plot(omega_ce[bidx], alt[bidx], label='e-')
ax.legend()
ax.set_xscale('log')
ax.set_xlabel('Gyrofrequency (Hz)')
ax.set_ylabel('Altitude (m)')

Text(0, 0.5, 'Altitude (m)')

Gyroradius

The gyroradius is the radius at which charged particles gyrate around the magnetic field.

\[ r_\alpha = \frac{v_{T\alpha}}{\omega_{c\alpha}} = \frac{\sqrt{k_B T_\alpha m_\alpha}}{eB} \]

with h5py.File(filename, 'r') as h5:
    B = h5['Geomag/Babs'][:]
    ion_mass = h5['FittedParams/IonMass'][:]
    te = h5['FittedParams/Fits'][:,:,:,-1,1]
    ti = h5['FittedParams/Fits'][:,:,:,:-1,1]

omega_ce = const.e*B/const.m_e
v_te = np.sqrt(const.k*te/const.m_e)
re = v_te/omega_ce
omega_ci = const.e*np.expand_dims(B,axis=-1)/(ion_mass*const.m_p)
v_ti = np.sqrt(const.k*ti/(ion_mass*const.m_p))
ri = v_ti/omega_ci
print(re.shape, ri.shape)

(188, 11, 74) (188, 11, 74, 5)

/tmp/ipykernel_2643/3104341878.py:8: RuntimeWarning: invalid value encountered in sqrt
  v_te = np.sqrt(const.k*te/const.m_e)
/tmp/ipykernel_2643/3104341878.py:11: RuntimeWarning: invalid value encountered in sqrt
  v_ti = np.sqrt(const.k*ti/(ion_mass*const.m_p))

fig = plt.figure()
ax = fig.add_subplot(111)
for i in range(len(ion_mass)):
    ax.plot(ri[tidx,bidx,:,i], alt[bidx], label=mass_labels[int(ion_mass[i])])
ax.plot(re[tidx,bidx,:], alt[bidx], label='e-')
ax.legend()
ax.set_xscale('log')
ax.set_xlabel('Gyroradius (m)')
ax.set_ylabel('Altitude (m)')

Text(0, 0.5, 'Altitude (m)')

Plasma Frequency

This describes the response of particles to time-varying electric fields.

\[ \omega_{p\alpha} = \sqrt{\frac{n_\alpha e_\alpha^2}{\epsilon_0 m_\alpha}} \]

with h5py.File(filename, 'r') as h5:
    ne = h5['FittedParams/Ne'][:]
    frac = h5['FittedParams/Fits'][:,:,:,:,0]
    ion_mass = h5['FittedParams/IonMass'][:]
    
ni = np.expand_dims(ne,axis=-1)*frac[:,:,:,:-1]
omega_pe = np.sqrt((ne*const.e**2)/(const.epsilon_0*const.m_e))
omega_pi = np.sqrt((ni*const.e**2)/(const.epsilon_0*ion_mass*const.m_p))
print(omega_pe.shape, omega_pi.shape)

(188, 11, 74) (188, 11, 74, 5)

/tmp/ipykernel_2643/3299216418.py:7: RuntimeWarning: invalid value encountered in sqrt
  omega_pe = np.sqrt((ne*const.e**2)/(const.epsilon_0*const.m_e))
/tmp/ipykernel_2643/3299216418.py:8: RuntimeWarning: invalid value encountered in sqrt
  omega_pi = np.sqrt((ni*const.e**2)/(const.epsilon_0*ion_mass*const.m_p))

fig = plt.figure()
ax = fig.add_subplot(111)
for i in range(len(ion_mass)):
    ax.plot(omega_pi[tidx,bidx,:,i], alt[bidx], label=mass_labels[int(ion_mass[i])])
ax.plot(omega_pe[tidx,bidx,:], alt[bidx], label='e-')
ax.legend()
ax.set_xscale('log')
ax.set_xlabel('Plasma Frequency (Hz)')
ax.set_ylabel('Altitude (m)')

Text(0, 0.5, 'Altitude (m)')

Inertial Length

\[ d_\alpha = \frac{c}{\omega_{p\alpha}} \]

with h5py.File(filename, 'r') as h5:
    ne = h5['FittedParams/Ne'][:]
    frac = h5['FittedParams/Fits'][:,:,:,:,0]
    ion_mass = h5['FittedParams/IonMass'][:]

ni = np.expand_dims(ne,axis=-1)*frac[:,:,:,:-1]
omega_pe = np.sqrt((ne*const.e**2)/(const.epsilon_0*const.m_e))
omega_pi = np.sqrt((ni*const.e**2)/(const.epsilon_0*ion_mass*const.m_p))
di = const.c/omega_pi
de = const.c/omega_pe
print(de.shape, di.shape)

(188, 11, 74) (188, 11, 74, 5)

/tmp/ipykernel_2643/3274770970.py:7: RuntimeWarning: invalid value encountered in sqrt
  omega_pe = np.sqrt((ne*const.e**2)/(const.epsilon_0*const.m_e))
/tmp/ipykernel_2643/3274770970.py:8: RuntimeWarning: invalid value encountered in sqrt
  omega_pi = np.sqrt((ni*const.e**2)/(const.epsilon_0*ion_mass*const.m_p))
/tmp/ipykernel_2643/3274770970.py:9: RuntimeWarning: divide by zero encountered in divide
  di = const.c/omega_pi

fig = plt.figure()
ax = fig.add_subplot(111)
for i in range(len(ion_mass)):
    ax.plot(di[tidx,bidx,:,i], alt[bidx], label=mass_labels[int(ion_mass[i])])
ax.plot(de[tidx,bidx,:], alt[bidx], label='e-')
ax.legend()
ax.set_xscale('log')
ax.set_xlabel('Inertial Length (m)')
ax.set_ylabel('Altitude (m)')

Text(0, 0.5, 'Altitude (m)')

Hall Conductivity

\[ \sigma_H = \frac{n_e e^2}{m_e}\frac{\omega_{ce}}{\nu_e^2+\omega_{ce}^2} - \sum_i \frac{n_i e^2}{m_i}\frac{\omega_{ci}}{\nu_i^2+\omega_{ci}^2} \]

with h5py.File(filename, 'r') as h5:
    B = h5['Geomag/Babs'][:]
    ion_mass = h5['FittedParams/IonMass'][:]
    frac = h5['FittedParams/Fits'][:,:,:,:,0]
    nu_i = h5['FittedParams/Fits'][:,:,:,:-1,2]
    nu_e = h5['FittedParams/Fits'][:,:,:,-1,2]
    ne = h5['FittedParams/Ne'][:]

ni = np.expand_dims(ne,axis=-1)*frac[:,:,:,:-1]
omega_ce = const.e*B/(const.m_e)
omega_ci = const.e*np.expand_dims(B,axis=-1)/(ion_mass*const.m_p)
sigma_H = (ne*const.e**2)/const.m_e*omega_ce/(nu_e**2+omega_ce**2) - np.sum((ni*const.e**2)/(ion_mass*const.m_p)*omega_ci/(nu_i**2+omega_ci**2), axis=-1)
print(sigma_H.shape)

(188, 11, 74)

fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(sigma_H[tidx,bidx,:], alt[bidx])
# ax.legend()
ax.set_xscale('log')
ax.set_xlabel('Hall Conductivity (S/m)')
ax.set_ylabel('Altitude (m)')

Text(0, 0.5, 'Altitude (m)')

Pedersen Conductivity

\[ \sigma_P = \frac{n_e e^2}{m_e}\frac{\nu_e}{\nu_e^2+\omega_{ce}^2} + \sum_i \frac{n_i e^2}{m_i}\frac{\nu_i}{\nu_i^2+\omega_{ci}^2} \]

with h5py.File(filename, 'r') as h5:
    B = h5['Geomag/Babs'][:]
    ion_mass = h5['FittedParams/IonMass'][:]
    frac = h5['FittedParams/Fits'][:,:,:,:,0]
    nu_i = h5['FittedParams/Fits'][:,:,:,:-1,2]
    nu_e = h5['FittedParams/Fits'][:,:,:,-1,2]
    ne = h5['FittedParams/Ne'][:]

ni = np.expand_dims(ne,axis=-1)*frac[:,:,:,:-1]
omega_ce = const.e*B/(const.m_e)
omega_ci = const.e*np.expand_dims(B,axis=-1)/(ion_mass*const.m_p)
sigma_P = (ne*const.e**2)/const.m_e*nu_e/(nu_e**2+omega_ce**2) + np.sum((ni*const.e**2)/(ion_mass*const.m_p)*nu_i/(nu_i**2+omega_ci**2), axis=-1)
print(sigma_P.shape)

(188, 11, 74)

fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(sigma_P[tidx,bidx,:], alt[bidx])
# ax.legend()
ax.set_xscale('log')
ax.set_xlabel('Pedersen Conductivity (S/m)')
ax.set_ylabel('Altitude (m)')

Text(0, 0.5, 'Altitude (m)')