Looking at Matter from Two Scales


I'm currently taking the Plasma Physics MOOC from EPFL. One of the interesting aspects of the plasma state is that it is, over large scales, electrically neutral but on shorter scales indexed by the Debye length, locally non-neutral.

While discussing this over lunch running with a colleague, it occurred to me that this is something that can be nicely illustrated by a simulation using + and - charges. So this is what this notebook is about.

Simulating a medium with + and - particles

To simulate a medium, we just need a distribution of particles. Our distribution has a fundamental property: it has a scale by which the structure is measured. In our case, we will use a lattice of $(x,y)$ points that are all 1 unit apart. To make things look more realistic, we will add some gaussian noise on the order of .2 to each coordinate we draw.

In [1]:
import numpy as np
In [2]:
x_lattice_coords = np.arange(0, 100, 1.)
y_lattice_coords = np.arange(0, 80, 1.)

x, y = np.meshgrid(x_lattice_coords, y_lattice_coords)

We also define our center, for further use.

In [51]:
center = (50, 40)

Let's visualize our points:

In [4]:
%matplotlib inline
import matplotlib.pyplot as plt
In [5]:
In [6]:
In [7]:

Now, let's add a bit of noise to the lattice:

In [8]:
x += np.random.normal(size=x.shape)
y += np.random.normal(size=y.shape)

Our lattice is complete. Let's plot each point in it:

In [9]:
plt.plot(x.ravel(), y.ravel(), 'o')

Finally, let's assign some charges to our lattice. We do this by drawing random numbers between 0 and 1. Numbers above 0.5 get a + charge, while the others get a - charge.

In [10]:
charges = np.where(np.random.random(size=x.shape) > 0.5, 1, -1)
In [11]:
plt.scatter(x.ravel(), y.ravel(), c=charges.ravel())

As you can see, there's a lot of blue and red, meaning we've randomized the charges in a good manner.

The point of representing matter as a collection of particles over a square domain is to see what happens when we look into the matter from up close and what happens when we look at it from a distance. So let's zoom in!

Zooming into the matter

We will write a function that allows us to look at the previous plot in various ways.

In [12]:
from ipywidgets import interact
In [13]:
def zoom_in(box_size=(0.1, 10)):
    plt.scatter(x.ravel(), y.ravel(), c=charges.ravel())
    plt.xlim(center[0] - box_size / 2,
             center[0] + box_size / 2)
    plt.ylim(center[1] - box_size / 2,
         center[1] + box_size / 2)

It turns out we don't see the same thing when we look from up close at the charges and when we look from afar. We can do the following visualization to illustrate this point:

In [14]:
plt.figure(figsize=(10, 10))
for index, box_size in enumerate([3, 10, 25, 100]):
    plt.subplot(2, 2, index + 1)
    plt.title('zoom box size = {} units'.format(box_size))

In other words, the "average picture" or the "summed picture" of our particles depends on the scale we look at it.

We can animate this zooming effect easily using code by Jake Vanderplas. To embed animations inside the notebook, we have to define the following functions:

In [15]:
from tempfile import NamedTemporaryFile
import base64 

 {0}" type="video/mp4">
 Your browser does not support the video tag.

def anim_to_html(anim):
    if not hasattr(anim, '_encoded_video'):
        f = NamedTemporaryFile(suffix='.mp4', delete=False)
        anim.save(f.name, fps=20, extra_args=['-vcodec', 'libx264', '-pix_fmt', 'yuv420p'])
        video = open(f.name, "rb").read()
        anim._encoded_video = base64.b64encode(video).decode('utf-8')
    return VIDEO_TAG.format(anim._encoded_video)
In [16]:
from IPython.display import HTML

def display_animation(anim):
    return HTML(anim_to_html(anim))
In [17]:
from matplotlib import animation
In [18]:
fig, ax = plt.subplots()

def init():
    ax.scatter(x.ravel(), y.ravel(), c=charges.ravel())
# animation function.  This is called sequentially
def animate(i):
    box_size = np.linspace(0.1, 50, num=FRAMES)[i]
    ax.set_xlim(center[0] - box_size / 2,
             center[0] + box_size / 2)
    ax.set_ylim(center[1] - box_size / 2,
         center[1] + box_size / 2)
# call the animator.
anim = animation.FuncAnimation(fig, animate, init_func=init,
                               frames=FRAMES, interval=100)

# call our new function to display the animation