Micah P. Dombrowski / Jul 21 2018

Project 4b: Wave propagation in magnetized plasmas along B -- Faraday Rotation

Imported and adapted from the UCLA Plasma Simulation Group's JupyterPIC Notebooks 

Before anything else, we need to change to the stored notebook's directory, which contains input files for the PIC simulator.

cd /notebooks/faraday-rotation

In this project you will study the propagation of an EM wave into a magnetized plasma while moving along the magnetic field. In such a plasma R and L waves exist. As was discussed in class the R and L waves are circularly polarized. So if you launch a linearly polarized wave in to the plasma you have to first decompose it into right and left handed circularly polarized waves to see what happens when it enters the magnetized plasma. If the frequency of the wave is larger than the  ωR\omega_R (which is larger than  ωL\omega_L for the highest density then both R and L waves will exit and it will propagate through. We will consider this case for this project (this situation could also arise if the frequency of the EM wave was below  ωR\omega_R but below  ωce\omega_{ce} and above  ωL\omega_L).

However, if the frequency of the EM wave is lower then  ωR\omega_R at the highest density but higher than  ωL\omega_L and  ωce\omega_{ce} then the R wave will be reflected and the L wave will propage through. And if the plasma then ends a circularly polarized wave will exit the plasma. It will be right handed circularly polarized.

If the frequency of the EM wave is lower than  ωL\omega_L at the highest density but lower than  ωce\omega_{ce} then the L wave will be reflected and the R wave will propagate through. And if the plasma then ends a circularly polarized wave will exit the plasma. It will be left handed circularly polarized.

All of these cases can be studied in this project if ones changes the input files. However, for this project we will examine Faraday rotation so the EM wave will have a frequency above both  ωR\omega_R and  ωL\omega_L.

1.
Simulations with a Particle-in-Cell Code

In this project you simulate plasmas that have a uniform applied magnetic field in the  x1^\hat{x_1} (or  x^\hat{x}) direction. The value of B corresponds to  ωce/ω0\omega_{ce}/\omega_0=.5 where the incident EM wave has a frequency of  ω0\omega_0. Note that this is different normalization than before. You can think of this as where we are sending the same EM wave into plasmas with different plasma densities but with the same magnetic field.

Each plasma electron is initialized with positions (only in  xx or what we call  x1x_1) such that the density is non-uniform. The ions are initialized at the same positions but they have an infinite mass. Each electron is also initialized with velocities ( v1v_1,  v2v_2,  v3v_3)=(.05c,.05c,.05c) or momentum ( mv1mv_1,  mv2mv_2,  mv3mv_3) from a Maxwellian in each direction. The particles then begin to move in the self-consistent fields that their current and charge density produce. The density will start at 0, rise rapidly to a value of  n0=.14n_0=.14 or  n0=.42n_0=.42 (make sure you understand the units) at  x1x_1=300 and then fall back to 0 at  x1x_1=480.

  • The length of the simulation box is 1000  c/ω0c/\omega_0 
  • The simulation will run for a time 1100  1/ω01/\omega_0.
  • The simulation uses 50,000 particles.

An EM wave that is linearly polarized in the  E3E_3 direction is launched from the vacuum. You will be looking at plots of the electric field in the two transverse directions,  (E2,E3)(E_2,E_3).

# Here is the density profile,
# you can change the length of the plasma or the external magnetic field 
# to achieve the perfect change of polarization:

# Here we plot the density profile:    
from scipy.optimize import fsolve 
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

npeak=0.14
xstart = 300
xend = 480.0
xmax =1000.0
xarray=np.arange(0,xmax,10.0)
nx=xarray.shape[0]
den_array=np.zeros(xarray.shape[0])

for ik in range(1,nx):
    if (xarray[ik] < xstart):
        den_array[ik]=0
    elif (xarray[ik] < xend):
        den_array[ik]=npeak
    else:
        den_array[ik] = 0.0
        
plt.plot(xarray, den_array,'b',label='Plasma Density Profile')

plt.xlabel('Z $[c/\omega_{pe}]$'); plt.ylabel('Density $[n_c]$')
plt.title('Plasma Density Profile'); plt.legend()
plt.xlim([0,xarray[nx-1]]); plt.ylim([0,1.25*npeak])
plt.legend(loc=0)
plt.grid(b=True, which='major', axis='both')
plt.show()
import osiris
%matplotlib inline

2.
Case 1: 90 degree rotation

# case 1:  90 degree rotation:
dirname = 'faraday-rotation'
osiris.runosiris(rundir=dirname,inputfile='faraday-rotation.txt')
dirname = 'faraday-rotation'
osiris.plot_tx(rundir=dirname, b0_mag= 0.3, plot_or=2, show_theory=False,cmap='RdBu') 
osiris.plot_tx(rundir=dirname, b0_mag= 0.3, plot_or=3, show_theory=False,cmap='RdBu')

For this case,  n0=ωp02/ω02=.14n_0=\omega_{p0}^2/\omega_0^2=.14. Here we plot  E3E_3 and  E2E_2.

  • Can you explain the information shown in these two plots?

Next we make plots of  E2(t)E_2 (t) and  E3(t)E_3 (t) at positions in the vacuum both before ( x1=40x_1=40) and after ( x1=500x_1=500) the plasma.

from h5_utilities import *

e2=read_hdf(dirname+'/e2.h5')
e3=read_hdf(dirname+'/e3.h5')

nt = e2.data.shape[0]
tmax=e2.axes[1].axis_max
dt=tmax/(nt)
time_axis=np.arange(0,tmax,dt)

f,(ax1,ax2,ax3,ax4)=plt.subplots(4,sharex=True,sharey=True)
ax1.plot(time_axis,e2.data[:,200],label='e2, incident')
ax1.legend()
ax2.plot(time_axis,e2.data[:,2500],label='e2, transmitted')
ax2.legend()
ax3.plot(time_axis,e3.data[:,200],label='e3, incident')
ax3.legend()
ax4.plot(time_axis,e3.data[:,2500],label='e3, transmitted')
ax4.legend()
plt.show
  • Can you explain what you see in these plots?

3.
Below we look at the polarization of the EM wave as it moves across the plasma

Next we make a graph where we plot dots for the values of the electric field in the  E2E_2 vs.  E3E_3 plane. Each dot corresponds to a different time. The colors correspond to different values of  x1x_1. Blue is for an  x1x_1 in the vacuum before the EM wave enters the plasma, red is for an  x1x_1 in the plasma, and green is for an  x1x_1 in the vacuum after the wave exits the plasma.

  • the blue line is the polarization as a function of time at  x1=40x_1=40 
  • the red line is the polarization as a function of time at  x1=420x_1=420 
  • the green line is the polarization as a function of time at  x1=560x_1=560 
dirname = 'faraday-rotation'

from h5_utilities import *

e2=read_hdf(dirname+'/e2.h5')
e3=read_hdf(dirname+'/e3.h5')

plt.scatter(e2.data[:,200],e3.data[:,200],c='b',label='in vacuum',s=2)
plt.scatter(e2.data[:,2100],e3.data[:,2100],c='r',label='in plasma',s=2)
plt.scatter(e2.data[:,2800],e3.data[:,2800],c='g',label='in vacuum (out)',s=2)
plt.legend()
plt.xlabel('$e_2$',fontsize=18)
plt.ylabel('$e_3$',fontsize=18)
plt.show()
  • Can you explain this graph?

4.
Case 2: 270 degree rotation --> this is achieved by changing the plasma density here

We next do a simulation in which everything is the same except the plasma density is increased to  n0=ωp02/ω02=.42n_0=\omega_{p0}^2/\omega_0^2=.42.

# case 1:  90 degree rotation:
dirname = 'faraday-rotation-270'
osiris.runosiris(rundir=dirname,inputfile='faraday-rotation-270.txt')
dirname = 'faraday-rotation-270'
osiris.plot_tx(rundir=dirname, b0_mag= 0.3, plot_or=2, show_theory=False,cmap='RdBu') 
osiris.plot_tx(rundir=dirname, b0_mag= 0.3, plot_or=3, show_theory=False,cmap='RdBu')

Here we plot  E3E_3 and  E2E_2.

  • Can you explain the differences in these plots as compared to the lower density case?

Next we make plots of  E2(t)E_2 (t) and  E3(t)E_3 (t) at positions in the vacuum both before ( x1=40x_1=40) and after ( x1=500x_1=500) the plasma.

from h5_utilities import *

e2=read_hdf(dirname+'/e2.h5')
e3=read_hdf(dirname+'/e3.h5')

nt = e2.data.shape[0]
tmax=e2.axes[1].axis_max
dt=tmax/(nt)
time_axis=np.arange(0,tmax,dt)

f,(ax1,ax2,ax3,ax4)=plt.subplots(4,sharex=True,sharey=True)
ax1.plot(time_axis,e2.data[:,200],label='e2, incident')
ax1.legend()
ax2.plot(time_axis,e2.data[:,2500],label='e2, transmitted')
ax2.legend()
ax3.plot(time_axis,e3.data[:,200],label='e3, incident')
ax3.legend()
ax4.plot(time_axis,e3.data[:,2500],label='e3, transmitted')
ax4.legend()
plt.show
  • Can you explain the differences in these plots from the low density case?

5.
Below we look at the polarization of the EM wave as it moves across the plasma

Next we make a graph where we plot dots for the values of the electric field in the  E2E_2 vs.  E3E_3 plane. Each dot corresponds to a different time. The colors correspond to different values of  x1x_1. Blue is for an  x1x_1 in the vacuum before the EM wave enters the plasma, red is for an  x1x_1 in the plasma, and green is for an  x1x_1 in the vacuum after the wave exits the plasma.

  • the blue line is the polarization as a function of time at  x1=40x_1=40 
  • the red line is the polarization as a function of time at  x1=420x_1=420 
  • the green line is the polarization as a function of time at  x1=560x_1=560 
dirname = 'faraday-rotation-270'

from h5_utilities import *

e2=read_hdf(dirname+'/e2.h5')
e3=read_hdf(dirname+'/e3.h5')

plt.scatter(e2.data[:,200],e3.data[:,200],c='b',label='in vacuum',s=2)
plt.scatter(e2.data[:,2100],e3.data[:,2100],c='r',label='in plasma',s=2)
plt.scatter(e2.data[:,2800],e3.data[:,2800],c='g',label='in vacuum (out)',s=2)
plt.legend()
plt.xlabel('$e_2$',fontsize=18)
plt.ylabel('$e_3$',fontsize=18)
plt.show()
  • Can you explain the differences in these plots as compared to the lower density case?