In recent years, ion traps have been used to confine small numbers of charged
particles for experimentation. In particular, the properties of confined
crystals and one-component plasmas have been of interest to researchers. Laser
cooling techniques have allowed one-component crystals to be formed in ion
traps. These crystals are interesting in connection with storage rings for
particle accelerators, as well as for potential applications to computer
components. On the theoretical side, one-component plasma crystals have some
unique properties, and the behavior of ions in a dynamic trap can be
classically chaotic, or it can be ordered, as in a crystal. Thus, properties of
the order-chaos transition can be studied in these systems.
The purpose of this thesis is to provide a greater understanding of the
properties of crystals that can be formed in the Paul Trap. The various types
of structures that are seen in different regions of the stability diagram will
be discussed and the approximate borders between these will be given.
One of the most popular ion traps used today is the Paul Trap. It is a dynamic
electric trap formed by a parabolic ring electrode and two parabolic endcaps.
There is an AC and DC voltage applied between the ring and the end caps, to set
up a time-varying electric field in the trap's center. It is well-known that
such a field produces a strong confining force on ions within it, provided that
the frequency and strength are appropriately tuned to the particle species. The
equation of motion for a single particle in the Paul Trap is,
m\frac{d^2}{dt^2} \left( \begin{array}{c} x \\ y \\ z \end{array} \right)\ +
2e\frac{[U_0 + V_0\cos(\omega_0 t)]}{r_0^2 + 2z_0^2}\left( \begin{array}{c} x \\ y \\ -2z \end{array} \right)\ = 0
Here, $U_0$ and $V_0$ are the DC and AC voltages, respectively, $r_0$ is the
radius of the ring, and $z_0$ is the distance from the center of the trap to
one of the endcaps. This equation can be rewritten in the mathematical form of
a Mathieu Equation by substituting
a=\frac{8eU_0}{m\omega_0^2 (r_0^2 + 2z_0^2)}\ \ ,\ \ q=\frac{4eV_0}{m\omega_0^2 (r_0^2 + 2z_0^2)}
And by replacing $t$ with $\tau\equiv \omega_0 t/2$. We now have 3 uncoupled
Mathieu-type equations,
\frac{d^2}{d\tau^2} \left( \begin{array}{c} x \\ y \\ z \end{array} \right)\ +
[a + 2q\cos(2\tau)]\left( \begin{array}{c} x \\ y \\ -2z \end{array} \right)\ = 0
The Mathieu equation has been studied extensively in the literature, and it is
known to have stable, quasi-periodic solutions only in certain regions of the
parameter space. In this case, the regions are bounded not only by the standard
Mathieu instability lines, for $a$ and $q$, but because of the $-2$ coefficient
on the $z$ component, the Paul Trap has a $z$-instability border bounded by the
Mathieu instability lines for $-2a$ and $-2q$. This gives rise to the stability
diagram for the Paul Trap, shown in figure 1.
[insert figure 1 here]
Figure 1: The stability diagram of the Paul Trap.
Most research on the Paul Trap has focused on the fundamental stability region,
(region A in figure 1). This is because at low values of $a$ and $q$, the trap
field can be approximated by a harmonic pseudo-potential. The pseudo-oscillator
frequency in the radial direction is given by
\omega_{xy}\ = \mu_x \frac{\omega_0}{2}
Where
\mu_x\ \equiv\ \sqrt{\frac{1}{2} q^2 + a}
Likewise, the pseudo-oscillator frequency in the $z$-direction is
\omega_z\ = \mu_z \frac{\omega_0}{2}
And
\mu_z\ \equiv\ \sqrt{2(q^2 - a)}
This pseudo-oscillator, it should be noted, is only a description of the
average, or \emph{secular} motion of the particle. There is also a
\emph{micromotion}, which the particle undergoes during a single cycle of the
trap, at the trap frequency $\omega_0$. Ordinarily the micromotion is quite
small compared to the secular motion of the particle, but at low temperatures
such as those found in the crystalline phase, the micromotion can make up most
of the particle's kinetic energy (assuming there are multiple particles in the
trap). One important consequence of the micromotion is that it causes the
system to heat when the particles' micromotions are out of phase. This will be
discussed in more detail in the next chapter.
When more than one particle is added to the trap, the equations of motion are
modified as follows,
\frac{d^2}{d\tau^2} \left( \begin{array}{c} x_i \\ y_i \\ z_i \end{array} \right)\ + [a + 2q\cos(2\tau)]\left( \begin{array}{c} x_i \\ y_i \\ -2z_i \end{array} \right)\ =\ C\sum_{j\not= i}^N q_iq_j\frac{(\vec{x}_i - \vec{x}_j)}{|\vec{x}_i\,- \vec{x}_j|^3}
Where $C=1/m\omega_0^2 \pi \epsilon_0$, $q_i$ represents the charge of
particle $i$, $x_i$ represents the $x$~coordinate of particle $i$ and $N$ is
the total number of particles. The $q_i$ and $q_j$ in the right hand side
of~(8) will generally be the same, though it is possible to include particles
of different species if the trap is set appropriately (or a particle and its
antiparticle, which allow the masses to be the same). If particles with
different masses are included, the $a$ and $q$ for each mass will be different
(as will $C$), thus requiring that the trap be tuned within a region where it
is stable for all the particle species.
In order to study crystalline structures, it is necessary to cool the particles
enough so that they crystallize. This is done in the simulations by adding a
damping factor, $-\gamma\omega_0\vec{v}_i$ to the right hand side of~(8), where
$0\leq \gamma\,< 1/2\pi$ and $\vec{v}_i$ signifies the appropriate component of
the velocity vector of particle $i$. For ions, $\gamma$ could be generated in a
real trap through the use of laser cooling. It is not clear exactly how cooling
could be achieved in a system of electrons, but for the purposes of studying
the theoretical properties of electron crystals, this is not important. The
important point is that the damping required to cool a system sufficiently will
vary depending on the number of particles and the strength of the trap field.
This dependence is due to rf heating, which is stronger for higher density
plasmas if the micromotions of the particles are out of phase. This means that
once a system has crystallized, slowly removing the damping will leave the
crystal intact, provided it is a stable configuration. Heating and crystal
stability will be discussed in more detail in the next chapter.
We now have 4 free parameters in our system, the trap parameters $a$ and $q$,
the number of particles, $N$, and the damping $\gamma$. Due to classical
scaling, the mass and charge of the particles, assuming they are all of the
same species, does not affect the structure of the system, only the scale.
This means that as long as quantum effects are negligible, which will be
discussed later for the case of electrons, we are free to choose whichever
particle species we are interested in, with the resulting structures being
revelant to other choices. Since we are interested in some properties of
electron crystals, we chose to use electrons for all the simulations.
The simulations were done with the use of a standard fourth order Runge-Kutta
integration scheme. The problem has $6N$ degrees of freedom and can be
separated into that many first order differential equations. Below is the
system of equations used in the program.
$x_i(t_{n+1}) = \Delta t \cdot\ \dot{x}_i(t_n)$
$y_i(t_{n+1}) = \Delta t \cdot\ \dot{y}_i(t_n)$
$z_i(t_{n+1}) = \Delta t \cdot\ \dot{z}_i(t_n)$
$\dot{x}_i(t_{n+1}) = \Delta t \cdot\ \left(\frac{e^2}{4m\pi \epsilon_0} \sum_{j\not=i}^N \frac{x_i(t_n) - x_j(t_n)}{\left| \vec{x}_i(t_n) -\vec{x}_j(t_n) \right|^3}\ -\ \frac{x_i(t_n)\,e(U_0 + V_0\sin(\omega_0 t_n))}{mr_0^2}\ - \gamma \dot{x}_i(t_n)\right)$
$\dot{y}_i(t_{n+1}) = \Delta t \cdot\ \left(\frac{e^2}{4m\pi \epsilon_0} \sum_{j\not=i}^N \frac{y_i(t_n) - y_j(t_n)}{\left| \vec{x}_i(t_n) -\vec{x}_j(t_n) \right|^3}\ -\ \frac{y_i(t_n)\,e(U_0 + V_0\sin(\omega_0 t_n))}{mr_0^2}\ - \gamma \dot{y}_i(t_n)\right)$
$\dot{z}_i(t_{n+1}) = \Delta t \cdot\ \left(\frac{e^2}{4m\pi \epsilon_0} \sum_{j\not=i}^N \frac{z_i(t_n) - z_j(t_n)}{\left| \vec{x}_j(t_n) -\vec{x}_j(t_n) \right|^3}\ +\ \frac{2z_i(t_n)\,e(U_0 + V_0\sin(\omega_0 t_n))}{mr_0^2}\ - \gamma \dot{z}_i(t_n)\right)$
Here $t_n$ signifies the value of time at the current time step, $\Delta t$ is
the size of the time step, $\vec{x}_i$ is the location of particle $i$, and
$x_i$ is the $x$ position of particle $i$. In addition, $r_0$ is the radius of
the ring electrode in the trap, with $z_0$, the distance from the trap center
to one of the endcaps, being set to $r_0/\sqrt{2}$. Since the integrator is a
fourth order Runge-Kutta, the equations are not iterated exactly as written
above, but the formulae above are the ones used within the integrator, with the
values of $t_n$ and the positions and velocities appropriately adjusted by the
integrator. Notice that the actual values of the voltages $U_0$ and $V_0$ and
the frequency $\omega_0$ are used, and that the positions and velocites are
measured in SI units. As previously noted, the values of $a$, $q$ and $\gamma$
are the only ones that matter in the structure of the results.
In order to ensure sufficient accuracy in the calculations, 1000 steps were
taken per cycle of the trap frequency. Since printing out results at each of
these steps would lead to massive data files, we decided to print results only
when the trap field was off, that is, in general when we looked at dynamic
simulations, we were strobing the system at the end of every cycle. Doing this
allows the structure of the system to be viewed without the micromotion. In
instances when the micromotion was important, it was necessary to run the trap
for only a few cycles in order to view the micromotion with sufficient
accuracy.
The program that integrates the system needed to be set up in such a way that
comparisons of the results for different values of the free parameters could be
easily compared to theoretical models for the trap. Since the pseudo-potential
approximation given above uses the trap frequency in addition to $a$ and $q$,
it is convenient to set the trap frequency (and size) to a constant value and
modify the voltages to change the trap parameters. The trap size was set to
$r_0=2mm$, with $z_0=r_0/\sqrt{2}$ as stated above. In order to make the
voltages reasonable, we defined the trap frequency to be a constant value such
that $V_0=100V$ when the trap parameters are set to $q=0.2$ and $a=0.0$. This
is the value of the trap parameters at which our first simulations were done
and is nicely placed within the fundamental stability region of the trap. This
gives the value of $\omega_0\approx 6.62\cdot 10^9$. The exact value can be
calculated from the equation for $q$ given above, if needed.
Now that I have described the system being studied and the simulations being
used, the next chapters will discuss the results of simulations and the
theoretical models developed or found in the literature. First crystals will be
discussed in more detail in the context of the Paul Trap, and then Wigner
crystals in particular will be discussed. Finally, some of the more complex
phenomena in the trap will be discussed.
The Paul Trap has several different regions of crystallization, corresponding
to different crystal structures found there. Even for two particles, where it
would be expected that the particles would align with either the $z$-axis or
the $x$-$y$~plane, it has been shown that there are stable configurations with
a set angle to the $z$-axis. As more particles are added to the trap, the
possibilities become more and more complex, but they essentially fit into a
few categories of similar structures, with some transition states near the
borders.
Before beginning a discussion of the particular configurations observed in our
simulations, it is necessary to determine what we mean by a crystalline
configuration, and to discuss the process of creating these crystals.
Ordinarily, a crystal is defined as stable, static configuration of particles
with some periodic structure. Since the Paul Trap has a dynamic trap potential,
the crystals formed int he trap cannot be static, there is always a breathing
motion caused by the varying trap field. Therefore, in order to define a
crystal in the Paul Trap, we fall back on the stability and periodicity of the
configuration. In general, we can call a configuration a crystal if the
particles return to the same positions after some period of time, generally one
cycle of the trap field. Since the particles returning to the same positions
implies that the energy of the configuration remains constants, we can also
introduce a new condition in the Paul Trap, which is that a crystal should not
heat. Therefore, a brief look at rf heating in the Paul Trap is in order.
If a gas of charged particles is placed into a Paul Trap, the system will heat.
The kneading of the gas by the time-varying field causes a larger number of
collisions than a static filed would produce, and allows the particles to gain
energy over time. The heating process in the Paul Trap is not completely
understood. When the cloud is very sparse, the heating rate drops off, as would
be expected, due to the decreasing number of collisions per time interval.
Before this regime is reached, however, the heating rate plotted against radius
levels off for a period, for reasons which are not known, and before this
levelling off period, there is a steady increase in the heating rate from very
small radii until it levels off. The reasons for a lower heating rate in
very dense clouds are not clear at this time either. For our purposes, it is
sufficient to realize that the gas heats, and that for larger numbers of
particles and stronger trap fields the heating is stronger. This means that
more intense damping is required to form crystals when there are more particles
or when the trap field is stronger. It also means that there can be regions in
which stable crystals will not form because any slight phase difference in the
micromotions of the particles will cause run away heating that will melt the
crystal when damping is removed. Very strong damping can flash freeze a gas of
particles, but such a configuration is only a crystal if the damping can be
adiabatically removed without the crystal melting.
Since we are demanding that a configuration remain stable and not heat when the
damping is removed, our first task was to confirm that such a condition could
fullfilled at all for a significant number of particles in the trap. We chose
300 particles for our first investigations, and used the program described in
the previous chapter to attempt to create a crystal. Figure 1 shows an
$x$-$z$~plot of the crystal after 1000 trap cycles, with the trap parameters
set to $a=0.0$ and $q=0.2$, using a damping of $\gamma =10^{-2}$.
[insert figure 1 here]
Figure 1: $x$-$z$~plot of the 300 particle crystal at $a=0.0$ and $q=0.2$
Several initial conditions were run to be sure that this was a representative
configuration. The calculations for this crystal were then extended for an
additional 1000 trap cycles with damping to allow it to completely settle. Now
to test that the crystal did not heat, the damping was removed entirely, and
the crystal was calculated out for another 3000 trap cycles. Removing the
damping all at once is faster in terms of computation time, but is somewhat
rough treatment of the crystal. If the damping was too high initially, which
can again be due to computational time constraints, the particles in the
crystal will have a phase lag that can cause the crystal to melt. In this case
the only effect of the phase lag was that the period of the crystal breathing
was slightly off from the coded trap cycle time. This can be seen in figure 2,
which is a plot of the energy of the crystal with respect to trap cycle. Note
that since the data is printed out only at the end of each cycle, the effects
of the micromotion are not seen in the first 2000 cycles. The periodic effect
seen after the damping is removed is caused by the period of the crystal
breathing being 1/1000 of a cycle longer than the cycle, hence the system is
essentially being strobed at a different point in the breathing cycle at the
end of each trap cycle. Thus the data is a slow motion replay of what is
actually happening over the course of a single cycle. The double period is
caused by the trap compressing the crystal along the radial direction, then
going back to zero field, then compressing the crystal in the $z$~direction and
finally returning to the zero field state at the end of the cycle. The reason
this causes the energy to go down is because the energy plotted is that of the
particle configuration itself, which is the total energy when the trap is at
zero field, but which gives up energy to the trap when it is switched on. That
is, the system of particles gains negative potential energy from the trap field
when it is not zero.
[insert figure 2 here]
Figure 2: Plot of the energy of a crystal versus time (trap cycles). The
damping on the crystal was removed after 2000 cycles. The periodic change is
caused by the phase lag, but the crystal does not heat overall.
>From figure 2, we can see that the crystal does not heat when the damping is
removed, and hence we can feel justified calling it a crystal. We can now
determine if a particular configuration is a crystal by removing the damping
adiabatically and seeing if the system heats. In general, doing so will cause
a crystal-like configuration to start to move a little and then to suddenly
break up into a gas, but if the system heats at all when the damping is reduced
it is likely it will melt eventually.
The structure of the crystal in figure 1 is somewhat interesting. It was
previously known that at larger particle numbers, crystals in the Paul Trap
show an onion layer structure of concentric shells. Specifically, in the
isotropic conditions around $a=q^2/2$, such structures were known to exist. It
was assumed that oblate and prolate versions of the onion layer structure would
be the standard configuration below and above the isotropic line, respectively.
Indeed, if the data for this crystal is examined, a remnant of the onion layers
is observed. The crystal in figure 1, however, appears to have another
structure as well, which is a flat layering in the $z$~direction. Layered
structures may be of interest to researchers studying crystalline beams in
particle accelerators, or perhaps as a method of storing information for
computers. We can see definite separations of 7 layers, one in the center and 3
above and below in the $x$-$z$~projection. The separation between these layers
is also evident from an $r$-$z$~plot of the data, which is shown in figure~3.
[insert figure 3 here]
Figure 3: An $r$-$z$~plot of the crystal in figure~1, showing the separation
of the shells into layers. The onion shell structure is still evident, though
broken up.
For comparison, figure~4 shows an $r$-$z$~plot of a crystal formed closer to
the isotropic conditions, at $a=0.01$, $q=0.2$.
[insert figure 4 here]
Figure 4: An $r$-$z$~plot of $a=0.01$, $q=0.2$ for comparison with figure~3.
This plot shows that the configuration in figure~3 is broken up compared to the
normal onion layer structure.
The distance between the layers in figure~1 is close to the two-particle
separation in the pseudo-oscillator in the $z$-direction. This implies that the
onion shell configuration can be flattened without any real structural change
until the deformation causes the $z$-separation of nearest neighbors to be on
the order of the two-particle separation. It is somewhat strange that the
shells break up into well defined layers despite the fact that all the
particles with $z$~positions near zero form a ring in the outer shell. What
appears to be happening is that the onion shells themselves break up as the
structure becomes more oblate, which can be seen by comparing figure~4 with
figure~5, which is another $r$-$z$~plot, this time of a spherical case at
$a=0.045$, $q=0.3$.
[insert figure~5 here]
Figure 5: an $r$-$z$~plot of the crystal formed at $a=0.045$, $q=0.3$. This
plot shows 4 concentric, unbroken shells.
>From figure~5, we can see that the 300 particle sphere has 4 onion shells, with
the central shell containing only a few particles. On the other hand, figure~4
shows that as the crystal becomes more oblate, the central shells disappear,
feeding the outer shells, so that only 2 well defined shells remain, with some
extra particles in an approximate third shell in the center.
By comparing these two figures with figure~3, we can see that this process
continues, with only two shells now remaining. In addition to losing shells,
when the crystal gets too flat, the shells are forced to separate into disks,
pressing more particles into a particular plane, while leaving a noticeable
break in the vertical direction. Some disks are just the rings that would be
left by taking a horizontal slice through a hollow spheroid, while others are
full disks formed by flattening out a hemispherical section. Because there are
two shells, some of the layers observed in figure~1 have a slice of each of the
shells in them. This causes the layers to appear to be only partially filled,
but it seems likely that the configuration of the layers, as shown in figure~6
will not fill in much more, but rather that new layers will start. This is
due to the fact that the gaps between particles in the interior and exterior of
the layers directly above and below the central one are caused by the gap
between shells.
[insert figure 6 here]
Figure 6: The structure of the layers of figure~1, shown in an
$x$-$y$~projection. Only the central layer and those above the origin are
shown, for clarity. The layers below the origin are essentially mirrors of
those above.
Previous calculations of crystals had focused on the spherical case, assuming
that the structure remained the same, only deformed, until it came near the
two dimensional regime. Our simulations show that there is an intermediate
phase with layered structure, which had been reported in [Blumel's paper in
conference proceedings] but not further developed. Data was gathered on the
region in which these layered crystals appear, with the results shown in
$a$-$q$~space in figure~7, along with the stability lines and the spherical
line. Attempts to gather data for higher values of $q$ ran into a regime in
which crystals could not be formed, that is, extremely high damping produced
configurations which melted when the damping was slowly removed. Some of the
discoveries in this region will be discussed in chapter~5.
[insert figure 7 here]
Figure 7: Parameter space plot of the stability lines and the spherical line
with data points that show runs with and without disk layering.
In 1934, Eugene Wigner wrote a paper on the quantum mechanics of electrons on
the surface of a metal. At the end of the paper, he predicted that at low
enough temperature, the eletrons would enter a purely classical, crystalline
state, which is now known as a Wigner crystal.
Wigner crystals have since been observed on the surface of He and in the
interface between two semiconductors. The properties of these systems have
been studied both experimentally and theoretically, and more recently, computer
simulations have been used. In most cases, the crystals were studied
in the presence of a substrate, like in the original system discussed in
Wigner's paper. The Paul Trap, however, gives us the opportunity to study
two dimensional electron crystals in the absence of a substrate. In this
chapter, I will discuss briefly some of the experimental research that has been
done previously, before moving on to a discussion of theory and my own
numerical results.
The first system in which Wigner crystals were experimentally observed was in
the suface state electrons of a He film. [Insert experimental discussion here]
One of the interesting properties of the two dimensional electron crystals in
the Paul Trap is that by tuning the trap appropriately, the crystals can
undergo a 2D-3D transition. Numerical results of where this transition occurs
were reported in [Schiffer's Paper]. %insert footnote here
The transition was explored for the value of the trap anisotropy
$\alpha\equiv\frac{\omega_z^2}{\omega_{xy}^2}$ as a function of particle
number~$N$. The best fit for the data at large $N$ was found to be
$\alpha\approx c(N-2)^{\beta}$ where $c=1.073$ and $\beta=.52$.
After observing this transition in the MD simulations, I managed to come up
with a reasonably good analytical model for the transition line in the
parameter space of the trap. The model assumes a constant lattice spacing $d$
and looks at the stability of an ion at the exact center of the trap. The
potential for our test charge due to the coulomb interaction is given by
\vspace{10mm}
V_c=\int_{\frac{d}{2}}^R \int_0^{2\pi} (\frac{Q\sigma}{4\pi \epsilon_0})\,\frac{\rho \,d\rho \,d\theta}{\sqrt{\rho^2\ + z^2}}
Where $Q$ is the charge of the test particle and $\sigma\equiv 4e/\pi\,d^2$ is
the charge density of the lattice. The stability condition for the charge is
that it becomes unstable when, in the limit of small $z$, the total force is in
the same direction as $z$. We would like to find the value of $\omega_z$, the
pseudo-oscillator frequency in the $z$-direction, at which this occurs. This
gives us the following condition:
-\frac{dV_c}{dz}\, - m\omega_z^2\,z = 0
Integrating (1) and then taking the derivative with respect to $z$, we obtain
\frac{Q\sigma}{2\epsilon_0}\,(\frac{z}{\sqrt{d^2/4 + z^2}} - \frac{z}{\sqrt{R^2 + z^2}}) = m\omega_z^2\,z
If we divide through by z and then take the limit where $R~\gg~d~\gg~z$, we
have $\omega_z\,=(\frac{4e^2}{\pi\epsilon_0 md^3})^\frac{1}{2}$. Using
$\omega_z\,=\frac{\omega_0}{2}\,\sqrt{2q^2 - 2a}$ and the previous equation,
and substituting $d^3=\frac{2e^2}{m\pi\epsilon_0\omega_0^2\,(a+q^2/2)}$, which
is the equation for the two-particle spacing, we get
(q^2-a)=4(a+q^2/2)
which simplifies to $a=-q^2/5$. Looking at the data, it was clear that this did
not coincide with the observed transition line, so I looked for a way to
improve my approximation. The formula for $d$ was the most likely source of
error, since the lattice constant for 300 particles is not the same as the
two-particle spacing. To get a rough estimate of the real lattice constant, I
took the average nearest neighbor separation of a particular 2D crystal, and
divided by the two-particle spacing for the appropriate trap parameters. I did
this for 5 different crystals and averaged the results, obtaining the answer
$d_{300}\approx\,.685d_2$, where the subscripts denote total particle number.
Since $\frac{d_2^3}{d_{300}^3}\,=3.11$, we can revise (4) to
$(q^2-a)=4*3.11(a+q^2/2)$ and simplify to get the revised transition line
a\approx\,-.388q^2
This line fits the data quite nicely, as shown in figure 1.
In [Dubin's article], %insert footnote here
a more elaborate analytical analysis of the value for which the transition
occurs was presented. The analysis used an approximation of the crystal as a
regular hexagonal lattice in an anisotropic harmonic oscillator potential. By
analyzing the correlation energy $E_{corr}$, the author determined that the
ground state configuration depended only on the lattice type and the
dimensionless parameter $\sigma a_0^2$, where $\sigma$ is the number of charges
per unit area and $a_0$ is a unit of measure on the order of the two-particle
separation for the same trapping potentials, that is $a_0\equiv (\frac{q^2}{2\pi \epsilon_0 m\omega_{z}^2})^{\frac{1}{3}}$.
By approximating the finite system of $N$ particles as a fluid of total charge
$Ne$, we can get number density of the fluid as a function of radius:
\sigma (r)=\frac{3N}{2\pi R^2}\sqrt{1-\frac{r^2}{R^2}}
where $R$ is the radius of the fluid, given by $R=(3\alpha Na_0^3/16)^{\frac{1}{3}}$. Here $\alpha$ is still the trap anisotropy. The planar crystal is no
longer the minimum energy configuration at $\sigma a_0^2=w_1$, where $w_1=10/9$
. Setting $\sigma a_0^2$ equal to this value, and plugging the equation for $R$
into (6), with $r=0$, to get the smallest value of $\alpha$ at which the
planar crystal is unstable, we find the following formula for $\alpha_3$, the
transition value, as a function of $N$:
\alpha_3(N)=(\frac{96N}{\pi^3w_1^3})^{\frac{1}{2}}
Notice that at large $N$ this formula is very close to the numerical results
given before, where $\alpha_3 =1.073(N-2)^{.52}$ was the best fit. Also note
that for $N=300$ we can solve for $a$ as a function of $q$, since
$\alpha\equiv\frac{2(q^2-a)}{(a+q^2/2)}$. Doing so yields $a=-.393q^2$, which
is remarkably close to the approximate value given in equation~(5). The three
lines are plotted in figure~1 along with the results of my MD simulations for
300 particles, which are discussed below.
\psfig{figure=fig1.ps,height=6cm,angle=270}
Figure 1: $x$-axis: $q$, $y$-axis: $a$. Three lines: The line given by
Schiffer, the line given by Dubin and the line given by our model. 2 sets of
points: the points where monolayer crystals are formed and the points where
they are not.
The data points in figure~1 were determined by running the molecular dynamics
program at the indicated trap parameters with 300 particles for 1000 cycles at
$10^{-2}$ damping. At the end of 1000 cycles, the $x$-$z$ projection of the
particle positions was examined. If the total extent of the particles in the
$z$-direction was small (generally less than $10^{-10}$ meters) and the
particles appeared to be randomly distributed in the $z$-direction, then the
$z$~positions were assumed to be due to rounding errors, and the point ($a,q$)
was considered to be a monolayer crystal point. If the particles were not
randomly distributed, or if the $z$~extent was large (on the order of a
micron), then the point ($a,q$) was considered a no-monolayer point. Typical
non-random distributions were when there was a clear line of particles along
the $x$-axis of the $x$-$z$ plot but also quite a few particles sitting off of
this line. Figure~2 is an example of one such case, notice that the total
extent is still quite small, yet this type of distribution cannot be considered
a 2D~crystal.
\psfig{figure=fig2.ps,height=6cm,angle=270}
Figure 2: An $x$-$z$ plot of a non-monolayer crystal with only central
instability ($q=0.25$, $a=-0.025$)
Since each point represents a run which started from random
initial conditions, there is no indication that structures such as the one in
figure~2 are due to the effects of a crystal melting above the transition line.
Rather, it seems likely that figure~2 is the configuration seen very near the
transition line, where the 2D~crystal is unstable only in the center, and thus
the majority of the particles are still confined to a 2D~lattice. If we look
at an $x$-$y$ projection of the structure in figure~2, it is clear that the
hexagonal lattice is still preserved over most of the crystal, as shown in
figure~3.
\psfig{figure=fig3.ps,height=6cm,angle=270}
Figure 3: An $x$-$y$ plot of the structure in figure~2, showing the hexagonal
lattice configuration is still present.
The presence of configurations such as the one shown in figures~2~and~3, while
still in keeping with the stability model used to present the analytical
transition line, nonetheless lead to the question of what the structure of this
transition line may be. The stability analysis predicts a well-defined line
where the 2D~hexagonal lattice ceases to be the minimum energy configuration.
The model used in this analysis, however, is only an approximation of the real
situation in the Paul Trap. At small $a$ and $q$, the approximation of a
harmonic oscillator potential remains mostly valid for the secular motion of
the particles, but the micromotion may actually play a role in the transition.
We know that the micromotion will cause heating in non-crystalline states, so
it is possible that there are regions below the transition line, which cannot
form crystals. Also, it is frequently the case that in non-linear systems, the
border between two phases is fractal.
In light of these possibilities, I endeavored to illuminate the transition line
with greater resolution and wider scope. I was unsuccessful in this attempt,
but I will nonetheless present the method I used and the problems I encountered
for the benefit of future researchers. In order to gain a greater resolution,
it became clear that it would be necessary to write a program, which would be
capable of modifying the trap parameters and determining whether the new point
was a monolayer crystal point or not. The method I used to implement this
was to start with an already formed crystal in one of the two phases, and then
slowly change the trap parameters in such a way as to move towards some point
on the transition line. To begin with, I chose one of the monolayer crystals,
slowly removed the damping, and then used that as the starting condition, from
which I modified the trap parameters slightly every few cycles, and checked to
see if any of the particles had gained an appreciable $z$~displacement. This
method was unviable because the planar configuration is always a local minumum,
so the crystal would remain flat far up into the 3D~region. In order to get
the crystal away from the local minumum, I attempted to perturb the velocities
in the $z$-direction every 10~cycles. This led to a build-up of disorder in the
crystal and meant that I would need damping to allow the system to settle.
Since damping the crystal would require many more cycles each time the trap
parameters were changed, and since it seemed desirable to explore the
transition without damping as a parameter, I abandoned this method in favor of
starting from above the transition line.
By starting from above the transition line and moving slowly downwards, it
seemed that I would be able to determine the structure of the border in a
relatively short amount of time. I expected that the configuration the
particles were in just above the transition line would collapse when the trap
was modified such that the line was crossed. If this had been the case, it
would have been sufficient to check whether or not the furthest-displaced
particle in the $z$-direction was below some critical level, that is, we would
expect to see a sharp drop-off in the maximum $z$-displacement when the
transition line was crossed. There was no sharp transition, however, and so
the simple $z$ cut-off failed to give useful results. I attempted to develop
some other quantitative criterion for the transition based on statistical
analysis, but nothing worked well enough to give a better picture than that of
figure~1. It seems that the most accurate way to explore the transition would
be to start from a 2D~crystal, change the trap parameters, perturb the crystal,
allow the crystal to settle, and repeat. This method would be time-consuming,
but it could be automated and might result in a better resolution for the
border between 2D~crystals and 3D~crystals.
With the 2D/3D~transition better understood, we turned to some more general
questions about these 2D~crystals. Since our interest in these structures is
not as ion lattices, but rather as Wigner crystals, that is, \emph{electron}
crystals, the question of whether or not our classical simulations are at all
valid for electrons becomes important. Ordinarily, in the Paul Trap, as long as
the particles are all of the same species, the only free parameters are the two
trap parameters and the damping. This ceases to be true when quantum effects
become measurable, however, so an analysis of these crystals has to be done to
determine whether or not they are in the classical regime.
To begin with, we write down the Hamiltonian for a particle at the center of
our crystal, where the charge density is greatest and hence, where quantum
effects are likely to be strongest.
\hat{H}\ =\ \frac{p^2}{2m}\ + \frac{1}{2} m\omega_r^2 \rho^2\ + \frac{1}{4\pi \epsilon_0} \sum_{i=1}^N \frac{e^2}{|\vec{r_i}\ - \vec{\rho} |}
If we expand the interaction term to second order around $\rho =0$, we get the
following, denoting the interaction as $V_c$:
V_c \approx\ \frac{1}{4\pi \epsilon_0}\ \left[ \sum_{i=1}^N \frac{e^2}{r_i}\
+ \rho \sum_{i=1}^N \frac{e^2}{r_i^2}\,\cos \theta_i\ + \frac{1}{2} \rho^2
\sum_{i=1}^N \frac{e^2}{r_i^3}\,(3\cos^2 \theta_i\ - 1) \right]
Since $\sum_{i=1}^N \frac{e^2}{r_i}$ is a constant, we can simply redefine the
zero point of our potential to remove it. In addition,
$\sum_{i=1}^N \frac{e^2}{r_i^2}\,\cos \theta_i$ is just the force acting on the
particle in the center from all the other particles in the lattice, which is
zero by definition. Thus, the first two terms of~(9) drop out, and plugging the
remaining term back into~(8) leaves us with a harmonic potential around $z=0$.
\hat{H}\ = \frac{p^2}{2m}\ + \frac{1}{2} m\omega_r^2 \rho^2\ +
\frac{1}{2} \rho^2\ \frac{1}{4\pi \epsilon_0}\ \sum_{i=1}^N \frac{e^2}{r_i^3}\,(3\cos^2 \theta_i\ - 1)
By defining an effective fequency,
\omega_{eff}^2\ \equiv\ \omega_r^2\ + \frac{1}{m}\,\frac{1}{4\pi \epsilon_0}\,\sum_{i=1}^N \frac{e^2}{r_i^3}\,(3\cos^2 \theta_i\ - 1)
We can immediately write down the wavefunctions:
\Psi_n(\rho) = A_n\left[ \frac{1}{\sqrt{2m}} \left(\frac{\hbar}{i}
\frac{d}{d\rho}\ + im\omega_{eff} \rho \right)^n\right] \exp^{\frac{-m\omega_{eff}}{2\hbar}\,x^2}
>From (12), we can see that the oscillator length,
$b\equiv\,\sqrt{\frac{\hbar}{m\omega_{eff}}}$ needs to be smaller than the
lattice spacing near $\rho =0$. To find $b$ we need to know $\omega_{eff}$, so
we again turn to the fluid approximation and approximate the sum in~(11) as an
integral over a disk from one lattice spacing to the full radius. The lattice
spacing defined in terms of the number density given in~(6) is
$d=\sqrt{\frac{2}{\sqrt{3} \sigma}}$. The equation for $\omega_{eff}$ now becomes
\omega_{eff}^2\ = \omega_r^2\ + \frac{e^2\sigma}{4m\pi \epsilon_0} \int_0^{2\pi} \int_d^R \frac{(3\cos \theta\ - 1)}{r^3}\,rdrd\theta
Evaluating the integral and making the appropriate substitutions for $\sigma$,
$R$ and $d$, we arrive at the final (approximated) formula for $\omega_{eff}$,
\omega_{eff}\ =\left(\sqrt{\frac{3\sqrt{3} N}{\pi}}\ - 1\right)^{\frac{1}{2}} \omega_r
With some more calculation using the formula for $b$, we finally arrive at the
formula for the value of $\omega_r$ at which $b$ becomes equal to $d$:
\omega_r\ = \frac{m\pi e^4}{144\sqrt{3} \hbar^3 \epsilon_0^2}\ \frac{\left(\sqrt{\frac{3\sqrt{3} N}{\pi}}\ - 1\right)^{\frac{3}{2}}}{N}
For $N=300$ particles, this formula gives $\omega_r\,=2.7\cdot 10^{16}$ which
quite a long ways out of range. Of course if we were interested in
$b\approx d/10$, or any other factor, it scales as the factor to the sixth
power, hence for 300 particles, the one-tenth case would require
$\omega_r\,=2.7\cdot 10^{10}$, which still out of range, but considerably less
so. Therefore, we can be fairly certain that the classical calculations of
these 2D~crystals are valid.
Now, we briefly turn our attention to the possibility of reaching the quantum
regime with 2D~electron crystals in the Paul Trap. From equation (15), we can
see that there will be no significant overlap of wave functions, assuming the
approximations hold, anywhere in the fundamental stability region for 300
particles. We may, however, be able to reach the quantum regime with larger
crystals. Looking at the equation, we see that at large $N$ it scales as
$N^{\frac{-1}{4}}$. This means that in order to lower the required $\omega_r$
by a factor 10, we need to increase the number of particles by a factor of
10000. The borders of the fundamental stability region and the heating rate of
large numbers of particles will prevent 2D~quantum electron crystals from being
formed.
The subject of cylindrically confined ions has been studied previously, and is
of special interest to particle beam physics. Since particle beams use a
similar method of confinement to that used in the Paul Trap, the properties of
ions strongly confined in the radial direction, but loosely confined in the
vertical one were briefly invesitigated. Since much has been written on both
the infinite and finite cases of this system, the research done in our
simulations was mainly to verify the existence of one-dimensional crystals in
the Paul Trap. This chapter will present some of the research that has been
done previously and comment on the results found by our investigation.
In [Hasse and Schiffer] it was shown that the cylydrically confined, infinite
coulomb lattice exhibits shells which increase in number as the linear density
of the lattice increases. In addition, there are many transitions between
different lattice types as the increasing density changes the minimum energy
configuration. The first of these transitions is the transition from one
dimension to two, followed by the transition from two to three dimensions, then
a series of structural changes between lattice types, which continues even past
the point where a new inner shell is created. The discussion here will focus
only on the transition from one to two dimensions and from two to three. For a
more in depth discussion of the lattice type changes, I refer the reader to the
paper itself.
The paper by Hasse and Schiffer examined the infinite coulomb lattice using a
harmonic restoring force, $F=-qK\vec{r}$ in the radial direction and a periodic
bounding box to allow numerical calculation. The equations of motion of
particles experiencing the restoring force and mutual coulomb repulsion were
integrated numerically, with the contribution from electrons outside the
bounding box included as the periodic image of one of the electrons in the box.
The linear density, $\sigma$, was rescaled using the charge, $q$, and the
Wigner-Seitz radius, $a=(3q/2K)^{\frac{1}{3}}$, to be $\lambda=\sigma a/q$.
This allows the structure of the results to be independent of the size of the
bounding box, the number of particles and the restoring force and charge, since
all these parameters are included in $\lambda$. In our to form a lattice, the
sum of the kinetic energies of all the particles in each of the coordinate
directions was required to be within 20% of some value. If the energy for a
particular coordinate was outside that range, the velocity of that coordinate
was rescaled for each of the particles, so that it would be within the range.
After there were no corrections of this type for some number of time steps, the
target kinetic energy was reduced, and the process was continued to generate
the final crystal. It should be noted that this allowed the lattice
configurations quite a bit of manuverability in finding the minimum energy
state, but it would be unlikely to work in the Paul Trap, because the constant
heating would cause the energy corrections to be too frequent for the target
energy to ever be reduced, hence the damping method used in our simulations.
It was found that the transition from a one-dimensional string to a
two-dimensional zig-zag pattern occured at $\lambda=.709$ and the 2D-3D
transition occured at $\lambda=.964$. To convert these values to a transition
line in the parameter space of the Paul Trap, it is necessary to determine the
value of the trap anisotropy, $\alpha$, as a function of number of particles,
at which a given lambda is reached. This can be done approximately as follows.
First, we assume that the pseudo-potential approximation is valid, which will
hopefully be true for some decent number of particles in the trap. It might
break down, however, if the anisotropy required comes very close to the border
of the stability region. Secondly, we can get an approximate value of sigma by
assuming that the finite string has a constant particle spacing, which we will
call $d$. These assumptions give us the following relations:
a\equiv \left( \frac{3e}{2K} \right)^{\frac{1}{3}} = \left( \frac{3e^2}{2m\omega_r^2} \right)^{\frac{1}{3}}
m\omega_z^2 \frac{N}{2} d = 2e^2 \frac{1}{d^2} \sum_{n=1}^N \frac{1}{n^2}
\sigma \equiv \frac{e}{d} = e \left( \frac{m\omega_z^2 N}{4e^2} \right)^{\frac{1}{3}} \left( \sum_{n=1}^N \frac{1}{n^2} \right)^{-\frac{1}{3}}
\lambda \equiv \frac{\sigma a}{e} = \left( \frac{3N\alpha}{8} \right)^{\frac{1}{3}} \left( \sum_{n=1}^N \frac{1}{n^2} \right)^{-\frac{1}{3}}
\alpha = \lambda^3 \frac{8}{3N} \sum_{n=1}^N \frac{1}{n^2}
And we can use
a = \frac{2-\frac{\alpha}{2}}{2+\alpha} q^2
To determine the line in parameter space. Equation~(5) can lead to very small
values of $\alpha$ for large $N$, but the approximations break down as it nears
the stability border. Values of $\alpha$ that are smaller than $\sim0.015$
give incorrect results for equation~(6) even at $q\approx 0.2$. This limits the
number of particles which can be used to compare the results in our simulations
with those in the paper. Since 300 particles doesn't give a useful comparison,
we chose to use fifty particles, with the fifty particle crystal shown below
for $a=0.0393$ and $q=0.2$.
[insert figure 1 here]
Figure 1: String crystal of 50 particles at $a=0.0393$, $q=0.2$
Numerical results show that the approximation given by (1-6) is reasonably
accurate for low values of $q$ and $N$, however, it breaks down for higher
ones. The results reported in [Schiffer, transitions] give a best-fit line of
$\alpha \approx .395 N^{-1.73}$ for the transition to two dimensions. This line
does not fit with the observed existence of the crystal in figure~1, for which
$\alpha \approx .0236$. The prediction of equation~(5) is more accurate in this
case, though, as noted, the numerical best-fit is more accurate for large $N$.
It seems that the region of one-dimensional crystals in the Paul Trap is quite
small for large number of particles and is very close to the stability border,
which makes it hard to find at larger values of $q$. It is possible, however,
to form such crystals after some hunting for the appropriate trap settings.
The previous chapters has dealt with the crystalline formations present in
various regions of the Paul Trap for relatively low values of $q$. One reason
for this focus is that the pseudo-potential approximation can be used with
reasonably good accuracy for these low values. The other reason for focusing on
low strength trap settings is that the heating of a large number of particles
is slow enough to be damped away in a small number of cycles. The higher
strength regions of the trap are harder to get a theoretical handle on and in
general do not allow for any large crystals, but can nonetheless provide some
interesting phenomena and are worthy of a more in depth look. The results
presented below show a couple of the interesting effects that can be seen in
the stronger regions of the trap.
The first of the effects that we found at stronger trap settings showed up as a
part of a survey of the layered crystal structures described in chapter 2. To
determine the approximate border between the layered crystal region and the
onion shell region, we conducted a number of runs at different trap settings
ranging upward in strength. The runs at the highest values of $q$ failed to
crystallize due to the increased heating, but one of these runs, at $q=.35$,
$a=.01$, showed the following in an $x$-$z$~plot after 1000 cycles at $10^{-2}$
damping:
[insert figure 1 here]
Figure 1: $x$-$z$~plot of the run at $a=0.01$, $q=0.35$. The plot shows two
groups of particles that have separated off from the main group.
The plot shows two groups of particles separated from the main group. The
separation is unexpected in an approximately harmonic trap, but nonetheless is
a stable phenomenon. Different random initial conditions show that the effect
is caused by the trap conditions, and extending the run shows that it is not
transitory. In addition, several points taken nearby show that there is an area
of the trap, which is small but not unreasonably so, in which the effect is
seen.
Our first theory of how the effect could occur was that the potential had two
extra minima in the region of the separated groups, and that the damping was
keeping the system cool enough to create a noticeable separation. Taking a look
at the effect of the micromotion and watching for 100 steps per cycle over 10
cycles, it became apparent that there was some else going on. The particles
move in a pattern which has a period of 4 cycles. If you take the formation in
figure~1 as the starting point, the entire configuration becomes more prolate,
with the separations remaining. After a quarter-cycle, the formation begins to
flatten out again, returning to approximately the original position after a
second quarter-cycle. The particles in the separated groups then proceed to
move through the main group as it spreads outward in the radial direction. With
limited particle exchange with the main group, the particles move all the way
through the middle and reach a configuration similar to the one in figure~1 at
the end of one cycle. The formation then begins to become prolate again, as
expected, returning to look like figure~1 after another half of a cycle. Then,
rather than sliding through the middle again like in the previous cycle, the
particles that are separated and the particles in the middle flatten into three
near disks, with well-defined separation. The system then returns to look like
figure~1 at the end of the second cycle. At this point, the particles which
were originally on the top have moved to the bottom, so the total period is 4
cycles.
It was suggested to us that this could be a hydrodynamic mode of the system
corresponding to some particular ratio of the trap strengths. Since it has an
integer period when measured in trap cycles, it seems likely that it is a
resonance between one of the two pseudo-oscillators and the trap frequency, and
the orientation of the motion implies that it is the $z$-direction oscillator.
The standard pseudo-oscillator potential in not entirely accurate, but it does
give a ratio of the trap frequency to the $z$~frequency of 4.2:1, which is
close to the expected 4:1, but off by a bit more than is comfortable. Looking
at the variation in this ratio over the region in which the effect is observed
shows that the region itself would be much smaller than the uncertainty in the
calculated ratio. Also, it does not seem to exist all along a line pertaining
to an approximate 4:1 ratio, but rather is localized in a particular small
region. This implies that the ratio to the radial trapping potential is also
important, though the exact relationship between these ratios and the region is
not clear at this time.
An aspect of the Paul Trap, which has not been explored in any depth in the
literature is the structures that can be found in higher order stability
regions. The largest of these regions is around $q=1.9$, $a=-1.3$. In the
interest of exploring this region of the trap, we ran simulations at the trap
parameters above, for 2, 3, 6 and 10 particles. The intent was to see if there
was anything we could call a crystal being formed in this trap region. Even for
two particles, a damping of $10^{-2}$ was insufficient to overcome the heating.
At a damping of $10^{-1}$, the 2 particles form a crystal configuration which
is stable if the damping is slowly removed. For more than 2 particles, the
crystal-like configurations produced by strong damping seem to break up if the
damping is lowered, even slowly. Nonetheless, the particles all seem to occupy
a quasi-crystalline state, which appears to simply be unstable without
sufficient damping. The quasi-crystalline state is show below in
figures~2~and~3.
[insert figure 2 here]
Figure 2: A phase-space plot of the $x$-coordinates of 10 particles at $q=1.9$,
$a=-1.3$. Damping had to remain at nearly $10^{-1}$ to keep the crystal stable.
[insert figure 3 here]
Figure 3: A phase-space plot of the $z$-coordinates of the system in figure~2
These phase-space plots are for a single cycle of the trap. The system appears
to have stable periodic orbits in two-particle shells. Why this is the case is
not clear, but the two and six particle configurations are similar. For some
reason the configuration breaks up if the damping is removed even quite slowly.
It is difficult to tell from this preliminary investigation just what is going
on in this region of the trap, but it is clear that there are interesting
dynamics to be studied.
The Paul Trap is a popular method of studying ion dynamics because of its
ability to trap a single species of particle and also because of the ease of
changing its shape. To fully utilize this tool, it is important to understand
what can done with it. It has been the purpose of this thesis, and the research
that lead to it, to develop a better understanding of the types of crystalline
structures which can be found in the trap.
We have used molecular dynamics simulations and theoretical insight to allow us
to split the small $q$ region of the Paul Trap into different regions in which
there are one-dimensional, two-dimensional zig-zag, three-dimensional prolate
onion layer, spherical onion layer, oblate onion layer, oblate layered and
two-dimensional flat crystals. Figure 1 shows the divisions.
[insert figure 1 here]
Figure 1: The regions of the small $q$ section of the Paul Trap stability
diagram. The upper and lower lines are the stability borders. The regions
closest to the upper border are the one-dimensional and zig-zag regions. The
spherical case is just a single line, with the region above corresponding to
the prolate case and the region immediately below to the oblate case. The
next lower region is the layered region, and the lowest region is the 2DEC
region.
In addition to the stability diagram, we have found a couple of interesting
preliminary results in regions for higher values of $q$, and hopefully have
illuminated some possible areas of future research such as the higher order
stability regions and the fluid dynamics of charged plasmas in the Paul Trap.