An ion trap is a device that uses a spatially inhomogeneous electric field to store, separate, and study the quantum properties of charged particles.
Here is the simplest equation for the electric potential in which a particle might be held
\[\DeclareMathOperator{\Grad}{grad}
\DeclareMathOperator{\Div}{div} \Phi = \frac{U}{2r_0^2} \left( \alpha x^2 + \beta y^2 + \gamma z^2 \right),\]
where $U$ is the voltage in the setup, $r_0$ is a typical size of the field inhomogeneity, $\alpha$, $\beta$, $\gamma$ are dimensionless constants.
This potential is created in a vacuum by external electrodes, therefore, it should satisfy the Poisson equation $\Delta \Phi = 0$ (we can substitute $\vec{E} = - \Grad \Phi$ in Maxwell's equation $\Div \vec{E} = \rho / \varepsilon_0$ for vacuum, i.e. for $\rho=0$), where $\Delta$ is the Laplacian. The Laplacian of $\Phi$ is defined by:
\[\Delta \Phi = \frac{\partial^2 \Phi}{\partial x^2} + \frac{\partial^2 \Phi}{\partial y^2} +\frac{\partial^2 \Phi}{\partial z^2}.\]
The above expression shows that it's impossible to create such a static potential where the particle is stable in all three directions. This problem is solved by adding an oscillating with radio frequency electric field to the static field:
\[ \Phi = \frac{U}{2r_0^2} \left( \alpha x^2 + \beta y^2 + \gamma z^2 \right) + \frac{V}{2r_0^2} \cos( \omega_\text{rf} t ) \left( \tilde{\alpha} x^2 + \tilde{\beta} y^2 + \tilde{\gamma} z^2 \right),\]
where $V$ is the amplitude of the AC voltage in the setup, $\tilde{\alpha}$, $\tilde{\beta}$, $\tilde{\gamma}$ are dimensionless constants.
It is important, that $\omega_\text{fr}$ the frequency of the radio band, since this frequency range corresponds to wavelengths much larger than the characteristic size of the ion trap ($r_0 \approx 1~\text{mm}$) and therefore, wave effects can be neglected (in particular, this generally allows us to introduce a potential for an electric field).
So, if $\omega_\text{fr} < \omega_\text{fr,max}$ then the coefficients $\tilde{\alpha}$, $\tilde{\beta}$ and $\tilde{\gamma}$ satisfy the same equation, as the coefficients $\alpha$, $\beta$ and $\gamma$. The motion along each axis is independent and can be described by the Mathieu equation:
\[ \frac{d^2 x}{d \xi^2} + (a_x - 2q_x \cos 2\xi)x = 0, \tag{1}\]
where $\xi$ is a dimensionless variable which depends on the time, $a_x$ and $q_x$ are dimensionless constants.
Consider that the only force acting on the ion is the force of interaction with an electric field $\vec{E}$.
Solutions of the Mathieu equation for different constants $a_x$ and $q_x$ can be stable and unstable. You can study this issue yourself with the program Mathieu.py by setting $a_x$ and $q_x$ in the code (lines 7 and 8).
This program implements the simplest numerical integration (Euler's method). During the program cycle, the variable $\xi$ runs through the value from $0$ to $\xi_\text{max}$ in increments of $\Delta\xi$. According to the known values of $x$, $\xi$ the value of $d^2 x / d\xi^2$ is calculated with respect to the equation (1) at each iteration of the loop. Next, the approximation formulas
\[ \Delta \left( \frac{d x}{d \xi} \right) = \frac{d^2 x}{d \xi^2} \Delta \xi, \quad \Delta x = \frac{d x}{d \xi} \Delta \xi,\]
are used to calculate the values $x$ and $d x / d \xi$ for the next step of the program.
To study the stability for a large number of pairs of parameters $a_x$ and $q_x$, we will develop the following criterion for the automatic analysis. We will consider the solution stable if the maximum value of $x_{\text{max},1/2}$ reached in half of the simulation time differs less than $10\%$ from the maximum value of $x_\text{max}$ during the entire simulation time.
This algorithm is implemented in the program Mathieu-stab.py. In 31-37 lines you can define the ranges for $a_x$ and $q_x$ to study the stability. If you don't have programming skills in Python it is highly recommended to run this program ones with high resolution for both axes (it may possibly take several minutes) and manually write down some coordinates of the the stability region boundaries. Our criteria of stability are not perfect, so the graph may have some artifacts. You can manually precisely examine any region of parameters with the Matheiu.py.
Note that replacing $q_x\to-q_x$ does not affect the stability of the solution. This can be shown as follows: since $-q_x \cos (2\xi) = q_x \cos (-\pi +2 \xi)$, the simultaneous replacement of $q_x\to -q_x$, $\xi \to \xi - \pi/2$ does not change the form of the equation. It only changes the initial conditions, so it does not affect the stability of the solution.
The first common implementation of an ion trap is a cylindrical symmetric 3D RF trap. This trap has the geometry shown in the figure and the following selection of constants corresponds to this geometry.
$\alpha=1$ $\beta=1$ $\gamma=-2$ $\tilde{\alpha}=1$ $\tilde{\beta}=1$ $\tilde{\gamma}=-2$
The second common implementation of an ion trap is the linear RF trap, which is used for mass spectrometry. This trap has the geometry shown in the figure and the following set of constants corresponding to it:
$\alpha=1$ $\beta=-1$ $\gamma=0$ $\tilde{\alpha}=1$ $\tilde{\beta}=-1$ $\tilde{\gamma}=0$
