1) As a warm-up exercise, try to write a simple code to integrate the location and velocity of a non-relativistic electron feeling the Lorentz force in constant electric and magnetic fields. The particle can move in two dimensions. The electric field is in the plane of motion, while the magnetic field is out of the plane. For simplicity, you can start with a Runge-Kutta method (we will learn that it's not that great for this problem), and consider two cases: E smaller than B by magnitude and E greater than B (cgs units). Plot the orbit of the particle, and check quantities that should be conserved (if any). You can use Python, IDL, Matlab, Mathematica, or anything else that can provide quick plots to check your results.

2) Boris pusher

Particle in cell method requires accurate integration of relativistic particle motion. The time integration method of choice is usually leapfrog, but proper centering of quantities is needed. Consider a discretization of the nonrelativistic Lorentz force equation:

Boris' method for solving this implicit equation has three steps, involving defining new variables:

(this is half of the electric push)

(this is a rotation)

(this is the second half electric push)

First, prove that the second step above is the rotation of a vector. Second, prove that the new variables are related as:

This expression can be used to evolve particle trajectories. Implement this pusher in a compiled language of your choice (or Matlab/Python/IDL, etc), and test it on drift motion in constant electric and magnetic fields. Consider two-dimensional motion in x-y plane, with Bz and Ey, and try several ratios of Ey/Bz including smaller and larger than 1 (cgs units). For small Ey/Bz the expressions above should suffice, but otherwise, the Boris pusher has to be rewritten in the relativistic form, returning the Lorenz factors.

For pure magnetic rotation (Ey=0), check what happens as the timestep delta t increases. Is the integrator stable? What happens to the larmor radius of the partilce?

3) Using the FDTD scheme described in lecture, write a 1D solver for hyperbolic Maxwell equations in vacuum (dB/dt = - curl E, and dE/dt = curl B). Integrate the propagation of a transverse electromagnetic wave in a periodic box, and check whether the scheme is dispersive.

For the rest of the problems, you will need to install one of two codes, XOOPIC, or Tristan-MP. If you only have marginal interest in kinetic simulations, then you should only play with this web code for two-stream instability (problem 5 below): Two-Stream instability online.
You can also get similar functionality from an iOS application "es1" which will run electrostatic PIC on your phone or ipad! (It's available on App Store. I don't know if there is Android version too).

XOOPIC stands for X-windows Object Oriented Particle In Cell code. It's easy to run and easy to set up on Unix, but may have some graphics issues with Tk libary on the Mac (we'll keep you updated on this). If this is an issue, the easiest way is to install it on a remote Linux machine that you can access through VNC or X-windows tunnel. The packages are here: XGRAFIX-2.70.5 ,
Xoopic-2.70.2 .
Once configured and compiled, it will make an executable xoopic, which can be run with an input file from the directory input (or files below).

The user manual for XOOPIC sister comercial code OOPIC-PRO is here, which also describes the format of the input files for XOOPIC is here. A more full-fledged package is TRISTAN-MP, which is developed at Princeton and now is available for download and testing: from Github. It is recommended for serious enthusiasts. See Github page for intallation instructions. If you are using XOOPIC you can start with 4. For Tristan, it's easier to start with problem 5 and modify input files and user files in the distribution, then returning to prob 4.

4) Plasma oscillation project.

In this problem you will simulate the plasma oscillation in electron-ion plasma with stationary background ions. Start XOOPIC woth the input file plasmaosc.inp . In XOOPIC, start it with xoopic -i pathtofile/plasmaosc.inp , where "pathtofile" is where you put the downloaded file. There are several ways of initializing the plasma oscillation. One is to initialize the electron density with a spatially dependent profile, but this can give a rather noisy result, and requires initial Poisson equation solve. Another is to give electrons an initial kick. We follow the latter method, and give electrons a sinusoidal initial velocity. The sinusoid is initialized in slices (there is no easy way to create a spatially dependent velocity via the input file without changing the code). The number of wavelengths in the box is set by parameter "wl". In Tristan-MP, set up an oscillating velocity distribution in init_particle_distribution_user.

Read the plasmaosc.inp file and make sure you understand its structure. You will be modifying these input files later. Note that we initialized only electrons. Where are the ions?

Hit Run, and plot energy diagnostics ("TE-beige,..." diagnostic in the diagnostics window, TE stands for total energy), and the phase space("ux vs x for all species").
Measure the frequency of plasma oscillation from TE diagnostic, and compare it with the analytic estimate (plasma frequency is defined in the input file, albeit in MKS...). How well is energy conserved? This is a hard question to answer, because, unfortunately, the latest version of XOOPIC has introduced an underflow bug in the energy diagnostics. As you can see, the kinetic energy seems to turn to zero abruptly in part of the cycle -- this is a diagnostics problem. Hopefully, this is easier to answer with Tristan-MP.

Try different numerical parameters: number of particles per cell (numElectronsPerCell), the spatial resolution (ds, which is now set to be 1/20th of the skin depth), and timestep. See if you can make the timestep times omegap larger than 2. Does the code explode?

A simple extension of this problem is to study the plasma oscillation with a background magnetic field (called a hybrid oscillation). This can be initialized by adding B03init (in tesla) to the Control section of the input file.

5) Electron two-stream instability

Now we study two electron streams, one moving left, the other right. You can modify the input file, leaving only two load statements, one for each stream, or you can looking at how this is done in the two_Stream_ee_em.inp file in the oopic distribution. You can also use this file twostreamhw.inp. Give electrons velocity of 1e8m/s, and look at the phase space of the electrons. The linear theory of the instability is described in Birdsall and Langdon's book, excerpt of which is provided here BLTwoStream.pdf .

Check that the growth rate and the most unstable wavelength correspond to the theory. Then follow the evolution into the nonlinear regime. How does the instability saturate? Look at the phase space and the energy evolution with time. How do electron phase-space holes (vortices) form?

What happens if the beams have larger initial temperature?

Feel free to play with the prepackaged input files in the oopic distribution. Note that the code can do many boundary conditions and r-z cylindrical geometry, very useful for modeling laboratory plasmas.

Next we will try to simulate astrophysical applications, such as electromagnetic streaming instabilities and collisionless shocks.

In the following exercises, you will have to write your own input files. You can start by modifying the twostream input file above. Since I will not be there on Thursday or Friday, below I also provide the versions of the input files that I wrote. If you can avoid it, don't use them. This way you will learn how to use the code. If you are stuck, however, you can use my files to get going with the physics part.

6) Weibel instability in electron beams

Begin with an extension of twostreamhw.inp to the case of two counter-propagating electron beams in the immobile ion background. The streaming direction should be along y, not x, to avoid a bug wtth periodic BCs in XOOPIC. The beams fill all volume, and I recommend to start with a square box of 128x128 cells, with 10 cells/skin depth and at least 10 particles per cell. Watch the development of the instability in space and in magnetic energy, and look at the transition from electrostatic to electromagnetic modes as the transverse size of the box changes (transverse to the direction of motion). Vary the velocity from .1 of the speed of light to very close to speed of light (2.98e8 m/s). Good variables to plot include Ub (B^2 energy), Bz, total energy plot vs time, uy-y momentum space. Description of the instability and its growth rates can be found in Ch.2 and Appendix A of Medvedev & Loeb 1999, ApJ, 526, 697. Check that the growth rate makes sense. The growth rate should be comparable to omega_p; however, for relativistic flow omega_p includes the gamma factor of the flow in the denominaror (the electron mass becomes multiplied by gamma). What happens in the nonlinear stage and what saturates the instability? Same questions can be aswered with TRISTAN by using user_weibel.F90 and input.weibel files. You can stream along x direction.

7) Weibel instatbility in pair flows

Modify 6) to include the mobile positive charges (call them positrons). To do this you need to declare a new species, and load it in the domain. Have them start with the same loading as electrons, e.g. counterstreaming in the y direction. How did the instability change? What determines which filaments merge together? (think about currents they carry). Are filaments neutral or charged? Why? To clarify the last points, you may want to individually plot the up- and down- moving positrons and electrons. To do this, declare and load 4 species. E.g., electronsup, electronsdn, positronsup, positronsdn. It's the same simulation as above, but you just labeled particles according to their initial motion. What patterns do you see in charge separation? Can you explain this in the spirit of Medvedev & Loeb picture? Now change the mass of the positrons and make them 64 times the mass of the electrons. Describe the stages of the instability. For Tristan user file user_weibel.F90 can be used.

8) Weibel instability in the transverse plane

Have two neutral plasmas consisting of electrons and positrons counterstream in the direction normal to the plane of the simulation. Label one stream as plasmaelectrons and plasmaions, and the other as beamelectrons and beamions. Start with all streams of the same density, and send them at half a speed of light at each other. Note how the instability develops. What does magnetic energy look like? (in projection of the surface plot Ub) What is the direction of the generated B field? (you can plot vector B). What is the field topology in the nonlinear stage (both before and well after the peak of the field)? What is the energy of the B field at the peak compared to the initial kinetic energy of the system? Now, keep the density of the "plasma" particles the same, but decrease the density of the "beam" particles (both electrons and positrons), going down to 1/10th of the original (note that smaller than that you will start running out of macroparticles in the beam. You can circumvent this by changing the weight of the beam macroparticles, np2c variable). Describe any qualitative changes in the topology and evolution of the field. How does charge separation work here? Do you see any clumping of charge from "beam" or "plasma" particles? What is carrying the currents?

9) Collisionless shocks.

Modify the two-stream input file to have a conducting reflecting boundary on the left wall, and start with a stream of electrons and positrons, initially filling the whole domain, all moving to the left. This is a way to set up a collisionless shock which has no initial B field. You may need to enlarge the x dimension so there is enough time for the interaction between the incoming and reflected beams. You can also modify the right boundary condition to have a beam emitter which pumps in more particles with time. Can you form a shock? How long do you need to wait? What is mediating the shock? The best way to identify a shock is to look at the density of particles, and ask when the compression exceeds 2 (a factor of 2 you get just due to reflection off the wall). As you change the velocity of the flow from nonrelativistic to relativistic, does the shock mechanism change? Next, you can try adding external magnetic field to see how shock changes in the magnetized regime. You may see strange artifacts in the flow before the shock, and the density behind the shock maybe noisy, making it difficult to recognize what's happening. Another sign of shock formation is thermalization of the flow in the downstream, so you no longer see individual streaming beams in the momentum space plot, x-ux. This exercise is not guaranteed to work, as OOPIC is not well tuned for the relativistic shock regime. But you will get a sense of the physics involved.

You can see the OOPIC manual for reference on the boundary condition syntax.

My input files (never mind the comments in the headers -- I did not update them):
weibelpairs.inp
weibeltrans.inp
shock.inp Tristan is a better fit for this exercise. Use user_shock.F90 and input.shock.