## Free charges in conductors

There’s a great discussion in Griffiths book on electromagnetism (typical book used for junior/senior course for physics majors) about what happens to free charges inside a conductor. He discusses how they will always find their way to the surface, distributing themselves to cancel the field inside the conductor. If they don’t do that, there will be a residual field to push more charges around. This happens until the field is totally cancelled. In a footnote he mentions how this surface effect only happens in 3D. In lower dimensions, like charges on a plate, say, the charges don’t necessarily go to the boundary.

One thing I really like about his discussion is how he tries to get students to think about other ways the charges could distribute themselves to minimize their energy. He points out that, in one sense, it’s strange that the charges all go to the surface. They congregate there next to each other, even though they don’t like to be close to each other. Why don’t a few float to the interior to relieve some of that stress?

I’m teaching this course right now and, as I was preparing my thoughts for this section, I thought it might be fun to show students how we could model this. I’m writing this post for me, for my students, and for anyone else who might be interested in this approach.

### The model

What I decided to do was to model a finite number of charges dumped into a region of space where they were free to move, while interacting with each other via the Coulomb force. When they get to the boundary of the conductor, they feel a force that tries to keep them inside.

As usual, instead of directly calculating all the forces, I decided to approach this through a Lagrangian formalism. I calculate the potential of the system as a function of the locations of all the particles. Then, for each component of each particle, I can determine the equation of motion via the Euler-Lagrange equation:

$\frac{\partial L}{\partial x_i}-\frac{d}{dt}\frac{\partial L}{\partial x_i'}=0$

Where L is the difference between the kinetic and potential energy written in terms of the particle positions ($x_i$) and velocities ($x_i'$). Mathematica is perfectly happy to analytically figure out all of those equations of motion for me and to numerically integrate all of them to let me make pretty pictures (see below). The big trick, of course, is making sure you can correctly calculate the potential energy (typically, the kinetic energy is pretty easy: $\sum \frac{1}{2}m_i x_i'^2$).

The potential I ended up using has three parts:

1. The particles: For every pair of particles there’s a contribution that looks like this: $1/\left|\vec{r}_i-\vec{r}_2\right|$
1. Mathematica’s Subsets command helps with finding all the pairs
2. The wall: I’m using various polar functions ($r=f(\theta)$) to make the walls of the conductors. Following what I did in this post, I made the potential rise as a steep ramp as the position got further outside of the conductor. This forces the particle back into the conductor where this potential term goes back to zero.
3. External field: Sometimes I turn on an external electric field by simply supplying an extra term to the potential energy in the form of some constant times the y-component of the particle. This models a constant field in the -y direction.

Once I have the potential, Mathematica is ready to go.

By the way, one other thing I did was to add in some viscous friction so that the charges would settle down into their lowest energy state. Without this, the charges just continue to oscillate wildly. My students and I used this trick when trying to find the equilibrium state for a complex system of masses and springs. Basically we just set the initial conditions randomly and waited for the particles to settle down. Works like a charm.

### Spherical conductor

The first conductor I tried was a sphere. I dumped in 10 free charges and let them go nuts.

spherical conductor with 10 charges, no external field

It’s hard to see without being able to twist the image around, but, trust me, they end up on the surface.

Here’s what it looks like with an external field (directed to out of the face closest to you):

Spherical conductor with an external field

You can see how all the charges favor the side closest to you, and, once again, they all end up on the surface.

### 2D circular conductor

I wanted to see what would happen if I did the same thing in two dimensions. Here’s the no field case:

circular conductor with 10 charges. No external field.

Now here’s something cool! You can see that the lowest energy state for this case does not have all the charges going to the edge. One of the charges stays in the middle. Just like what Griffiths suggested in his footnote.

Here’s just the last frame for the case when there’s a downward directed external field:

circular conductor with external field (last frame only)

Once again you can see that there’s a charge in left in the interior while all of them are being forced downward, at least a little.

### Getting away from circles and spheres

I wanted to see if the same things would happen if the conductor shapes weren’t so simple. Here’s the result of a wavy sphere:

Non-spherical conductor, no external field

Yep, all on the surface.

Here’s the last frame of the case of a non-circular conductor in the 2D realm:

non-circular conductor in 2D with external field

I like this one because I did 20 charges and you can see that three charges stayed in the interior. Cool stuff, for sure.

I just wanted to finish with a some comments on how using Mathematica to do this was a great experience. First of all, everything you’ve seen was done with 21 lines of code. Some of those are function definitions, some are setting initial conditions, and some are for plotting. If I tried, I could probably get it down to under 10 lines of code (though one line would be quite long).

Next, keeping track of the variables was really easy. The basic format was r[i][j] where i ran from 1-3 to denote x, y, and z, and j counted the particles. That format made keeping track of things great and Mathematica was perfectly happy to take analytical derivatives with respect to each of them (and their time derivatives!).

The most complicated run I did was 20 particles for a 3D case. That meant integrating 60 coupled, second order differential equation. One command, NDSolve, took care of that for me. I didn’t have to convert the second order equations into first order ones, which is a nice bother to skip.

Plotting, and animating, was also very straightforward. Being able to easily adjust the opacity of the sphere, for example, was quite useful.

### Useful?

I hope this is useful to my students this coming week as we talk about this. I hope it’s useful to me in future years. Finally, I hope it’s useful to some of you!

Associate professor of physics at Hamline.
This entry was posted in mathematica, physics. Bookmark the permalink.

### 22 Responses to Free charges in conductors

1. bwfrank says:

I really enjoyed this. Is it possible in certain circumstances for charges to settle in a local-energy minimum, instead of finding a global minimum?

• Andy "SuperFly" Rundquist says:

That does happen sometimes, but I usually fiddle with the strength of the friction term to make sure most of the parameter space is explored in the initial stages.

2. bwfrank says:

This post has sent me down the rabbit hole. I was thinking more about 2D vs 3D, and why the center is an unstable equilibrium in 3D and stable in 2D. I think this is similar to how uniformly charged plate has constant force, where uniformly charged line has 1/z force. Like if our 3D conductor was a box, and our 2D conductor was a square. Then we introduce a new charge near but not exactly in middle. The new charge in 3D is approximately in equilibrium because force from parallel plates doesn’t depend on distance (ignore corner effects). That is, until charges on surface begin to run away, lessening the force on one side leading to run away motion toward the side it was initially off-center. In 2D, however the charge is not in equilibrium to begin with, rather it is driven back toward middle because 1/z force. The rate of run away charges can’t overcompensate for the 1/z dependence of force, and so it oscillates around middle (and settles with dissipation). Does this make any sense?

• Andy "SuperFly" Rundquist says:

You do have to be careful with Ehrenfest’s theorem (no stable equilibrium in electrostatics) when talking about this. However, the differences in the various dimensions are really interesting, and not just for this problem. Even the basics of how the field falls off, as you mention, is pretty cool. There isn’t a dimension where the field falls off in some crazy way, it’s always r^(some integer), which is nice and easy to deal with.

• bwfrank says:

Regarding stability, we are thinking about the situation as really like quasi-static, right? We ignore magnetic forces (despite moving charges), time-delayed electric forces (due to speed of light), and radiating E&M fields (despite acceleration). Merely move all charges and re-calculate coulombs forces. Not static, but not full dynamics either. For qualitative analysis, I’m certainly not considering every particle individually. So instead of calculating coulombs law, I’d calculate 1/z^n force and a new charge density due to run away charges. How to think about run away charges qualitatively is difficult. But anyway, the force is changing with proximity to surface (except 3D), but the effective charge on the surface changes as charges run away. It is a matter of who wins out, distance effect of run away charge effects. So, the question I’m now thinking about is: what fractional dimension “n” denotes the boundary between always run to the boundary vs. always run to the middle?

• Andy "SuperFly" Rundquist says:

whoops, sorry about that. I meant Earnshaw’s theorem, not Ehrenfest. That’s what I get for being too lazy to go check my book.

• Andy "SuperFly" Rundquist says:

First of all, I realized that Ehrenfest isn’t necessarily the right way to think about it. That talks about building a surface with a particular potential on the walls, not charges on a conductor.

You’re right that we’re talking about quasi-static situations (no mag field etc). I’ve had some fun over the last few days asking my wife and sons what they think should happen. My wife is pretty certain that things should go to the surface but my 11-yr-old was convinced that some charges should go into the interior. When I showed him the simulations, he was really interested.

• bwfrank says:

My brain won’t go away from this. I guess I’m still interested in what makes a really convincing argument about this. Up until this morning, I would have made some arguments which I think would have basically been the same whether or not it was 2D or 3D. I would have said something along the lines of the “If there were a non-zero E-field, then…” But, as you point out, any argument like this (one that doesn’t implicitly or explicitly involve attention to the number of dimension) is an argument for a possibility, but not a certainty.

The two dimensional situation certainly doesn’t have a non-zero E-field everywhere inside, and it doesn’t tend toward all charges going to the boundary. There isn’t even a stable static equilibrium where the inside charges are; rather the whole system exhibits a stable “quasi-static” equilibrium. In other words, the configuration as a whole is stable, even if the potential energies for inside charges are at unstable maximums.

Are there other systems like this? where objects reside at unstable equilibrium points, but remain stable because rest of the system responds quickly to perturbations in ways that turn that unstable equilibrium point into a “dynamic stability”?

Sorry to hog your comment space… you’ve just got my brain going.

3. Eli says:

Excellent. I’ve always tried talking through this argument with my students, but the Mathematica method is far better. Any interest in sharing your code?

• Andy "SuperFly" Rundquist says:

I sent Eli my code without comments. I’ll try to put comments in later and post a public version.

4. BlackGriffen says:

Hm, I’m not really convinced. It seems that the viscous friction is somehow permitting it to settle into an unstable equilibrium. My guess is that what you’re seeing is a problem of rounding error and/or insufficient runtime. See, a real viscous force should never reduce the velocities of the particles to zero. This is both in the mathematical asymptotic approach to zero sense and in the physical sense in which the damping force will eventually leave the object in a state of brownian diffusion. In any computer model, however, once the velocity drops below a certain value the particle will stop moving due to rounding error ( roughly v * delta(t) + position = position after rounding). If you’re using double precision floating point, for instance, 1.0 + 2^(-53) = 1.0.

There’s a really easy test to see whether I’m right. Simply put a low velocity cutoff to the viscous force. Then tune the cutoff to try to prove there there’s a critical cutoff, vcrit, below which the charges never stop moving and never stray far from the positions shown and above which they always fall into a moving equilibrium you would expect from theory. If I’m right then no such vcrit exists given sufficient runtime that will depend on how low the vcutoff is. A reasonable value would be, say, Tmax = 3 * D / vcutoff, so that an object traveling at vcutoff has sufficient time to travel/diffuse across the structure (the 3 is arbitrary).

• BlackGriffen says:

Oh, never mind. Here’s the deal – the force confining the charges to move in 2d while you’re using the 3d form of the fields is what’s doing the job. In other words you’re only moving halfway to 2d. If you go all the way, then your potential fields look like logarithms instead of 1/r and you’ll find that the 2d case behaves just like the 3d one. That’s what was throwing me off – I didn’t read that you were using a 3d potential. Never mind my previous post.

Put another way, you’re confusing the dynamics of 3d charges confined to the surface of an infinitely thin plate with real 2d dynamics. In this case it’s the forces that are holding the charge to the surface of the conducting plate that maintain the equilibrium.

• Andy "SuperFly" Rundquist says:

That’s exactly right. I’m using the point-charge form of Coulomb’s law but considering situations where the conductor is a 2D situation (which can certainly happen). The 1D case has a lot of history to it as well, but I haven’t modeled it here. Maybe soon I’ll tackle that.

I assume what you mean by going all the way is to consider several infinite lines of charge viewed end on. Then, of course, you’re right, the form of the Coulomb force changes and you get a different result. What I was doing here was to lay some groundwork to help my students understand Griffiths footnote. I like this was of thinking about it, too, and maybe I’ll bring it up in class today.

As for your other note about the oscillations, certainly the charges are still moving in the last frame of the simulation. I just stop the integration when I can tell that the equilibrium points have been identified.

• BlackGriffen says:

From the way you describe it I actually don’t like Griffith’s discussion. What is a 2d conductor, after all, but a surface? One thing you could mention is that you add more and more charges you’d see them move from the center and begin to occupy intermediate radii. Eventually the average charge density would look like a one of the problems in Jackson. Get them to glance at the “hard stuff.” ;)

• Andy "SuperFly" Rundquist says:

In class today it was interesting that the students wanted to focus on the 3D results. The 2D stuff seemed mildly interesting but they recognized the odd constraints required. We did talk about your version of 2D as well.

The 2D case sounds a little little like the formation of Wigner molecules. For example, if you have a square potential and put 4 electrons in, they repel each other so that one electron sits in each corner. If you add a 5th electron though, it energetically prefers to sit in the middle of the square. I wouldn’t be surprised if there were some similar commensurabiity relations in the finite simulations you did.

• Andy "SuperFly" Rundquist says:

Thanks, Rachel. I think they are the same, since the potential interaction is the same. My students asked in class If I would do the 3D case with a cubical conductor and 9 particles. With 8 they assume they’ll all go to the corners so they’re curious about 9. I’m off to work on that this afternoon!

6. eipi10 says:

Just curious if you’ve gotten around to creating a commented version for posting. I’d be interested in a copy of the source file, whether commented or not.

7. eipi10 says:

Thanks!

8. I like you idea. I did something similar in Visual Python back in 2008 (but having the computer directly calculate all of the forces). Unfortunately, there was no way to display this online, until now. There is now a branch of javascript called glowscript which has allowed me to reproduce these same visuals in a browser.