Homework 3

due 10/03/2018 at 11:59pm

Problem 1

  1. Write a version of binary_Coulomb_collisions that can handle any arbitrary number of particles.
  2. Modify plot_2D function so that you can print the results.

Problem 2

  1. Write a function that shows the trajectory of many electrons inside the lattice of fixed ions and use the right boundary conditions to generate plasma oscillations
  2. Measure the oscillation frequency and compare it to analytical results.

Solutions

Let's setup the main parameters and objects for this homework.

In [1]:
import math
from math import sqrt as sqrt
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline 
In [2]:
class charged_particle:
    m=9.10938356e-31 #charge 
    q=-1.60217662e-19 #mass
    x=np.zeros(3) #location
    v=np.zeros(3) #speed

Problem 1

We use here an array of object to pass all the particles to the function. Much faster and more general than the function binary_Coulomb_collision in the text of Ch03. Note that we added a cut-off radius r_cutoff to limit large acceleration between time steps. This allows to increase the time steps between iterations but it also increase the error on the trajectory.

In [7]:
def Coulomb_collisions(charged_particles,time_frame,n_steps,r_cutoff):
    N=len(charged_particles)
    loc=np.zeros((3,n_steps,N))
    dt=(time_frame[1]-time_frame[0])/n_steps
    for i in range(n_steps):
        for j in range(N):
            E=np.zeros(3)
            for k in range(N):
                R=charged_particles[j].x-charged_particles[k].x
                r=np.linalg.norm(R)
                if (r>r_cutoff):
                    E+=charged_particles[k].q/r**2*R/r
            charged_particles[j].v+=1/(4*np.pi*8.85e-12)*E*charged_particles[j].q/charged_particles[j].m*dt
        for j in range(N):
            charged_particles[j].x+=charged_particles[j].v*dt # we compute all the velocities then we compute the locations
            loc[:,i,j]=charged_particles[j].x
    return loc

That's plot_2D from the text.

In [8]:
def plot_2D(r,n_interpolated=1,plot_size=3,plot_axes='axes'):
    n = len(r[0,:,0])
    k = n_interpolated
    fig, ax = plt.subplots(1, 1, figsize=(plot_size, plot_size))
    # Now, we draw our points with a gradient of colors.
    for i in range(len(r[0,0,:])):
        x_interpolated = np.interp(np.arange(n * k), np.arange(n) * k, r[0,:,i])
        y_interpolated = np.interp(np.arange(n * k), np.arange(n) * k, r[1,:,i])
        ax.scatter(x_interpolated, y_interpolated, c=range(n*k), linewidths=0,
                   marker='o', s=2*plot_size, cmap=plt.cm.jet,) 
    #compute autoscale parameters
    xc=(r[0,:,:].max()+r[0,:,:].min())/2.
    x_low=xc-(r[0,:,:].max()-r[0,:,:].min())/2.*1.1
    x_high=xc+(r[0,:,:].max()-r[0,:,:].min())/2.*1.1
    yc=(r[1,:,:].max()+r[1,:,:].min())/2.
    y_low=yc-(r[1,:,:].max()-r[1,:,:].min())/2.*1.1
    y_high=yc+(r[1,:,:].max()-r[1,:,:].min())/2.*1.1
    #set autoscale parameters
    ax.set_xlim(x_low,x_high)
    ax.set_ylim(y_low,y_high)
    if (plot_axes!="axes"):
        ax.set_axis_off()
    else:
        plt.xlabel("x (m)")
        plt.ylabel("y (m)")

We now set the initial conditions

In [24]:
N=25
cp=np.empty(N,dtype=object)
for i in range (N):
    cp[i]=charged_particle()
    cp[i].x=np.array([(np.random.random_sample()-0.5)*30e-3,(np.random.random_sample()-0.5)*30e-3,0])
    cp[i].v=np.array([100*(np.random.random_sample()-0.5),100*(np.random.random_sample()-0.5),0])
time_frame=[0,4e-5]
n_steps=200

And then we run the code and plot the results.

In [25]:
r=Coulomb_collisions(cp,time_frame,n_steps,1e-10)
plot_2D(r,n_interpolated=1,plot_size=8)

It worked!

Problem 2

This is a difficult problem to adjust to get the right oscillations. The particles have to be far apart so the electrostatic forces are weak enough to allow for large time steps. We need an acceptable cut-off radius when computing the electrostatic force so the particles do not acquire large velocities between two consecutive time steps. Note that this is really to allow for fast computations rather than enabling oscillations to occur. Here we use $n_0=10^{13}m^{-3}$ which gives an oscillation period of 40ns, a time scale that has worked for us in previous problems. We also try to keep the size of the lattice at a scale that we have also used successfully. We used a length of $100 \mu m$. We take a radius cut-off on the order of $1\mu m$, a thousand smaller than the scale length of our lattice.

In [27]:
N=3
length=1e-4
n0=N**2/length**3
omega=sqrt(n0*1.6e-19**2/(8.85e-12*9e-31))
f=omega/(2*3.1416)
T=1/f
print('Density',n0,'m^-3 Frequency',f,'Hz Period',T,'s')
Density 8999999999999.998 m^-3 Frequency 27068704.021020308 Hz Period 3.694303204259229e-08 s

Now, starting small, we slowly increase the velocity so that oscillations start to appear. It is important to realize that these oscillations correspond to motion between two consecutive ions, not motion around one single ion. Check the difference between the velocity where oscillations between two ions occur ($v_x=6km/s$) and revolutions around a single ion ($v_x=4km/s$), or complete escape ($v_x=8 km/s$). It is important to note that complete escape would not happen if our lattice was repeating to infinity in all directions.

With these parameters we get particles that start at a given position $x$ and get back top that position in 60ns, which is close to the plasma frequency. We do not expect to find the exact period due to the limitation in lattice size.

In [28]:
delta_x=length/10
delta_y=length/10
lattice=np.linspace(-length,length,N)
lattice=np.linspace(-length,length,N)
cp=np.empty(0,dtype=object)
l=0
for i in range (N):
    for j in range (N):
        cp=np.append(cp,charged_particle())
        cp[l].x=np.array([lattice[i],lattice[j],0])
        cp[l].m=1.67e-27
        cp[l].q*=-1
        cp[l].v=np.array([0.,0,0])
        l+=1
for i in range (N):
    for j in range (N):
        cp=np.append(cp,charged_particle())
        cp[l].x=np.array([lattice[i]+delta_x,lattice[j]+delta_y,0])
        cp[l].v=np.array([6000.,0,0])
        l+=1
time_frame=[0,6e-8]
n_steps=200
r=Coulomb_collisions(cp,time_frame,n_steps,3e-6)
plot_2D(r,n_interpolated=3,plot_size=8)