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.