Systematic investigations on ion dynamics with noises in Paul trap

In the long time limit, the mean value of the ion position x and the speed v must be zero. Following [31, 36, 38], we choose to expand the stochastic average of the variance in terms of Fourier series basis set, i.e.

$$\begin{eqnarray} &&\left(\begin{array}{l} \left\langle x^2\right\rangle \\[2pt] \left\langle v^2\right\rangle \\ \langle x v\rangle \end{array}\right) = \left(\begin{array}{c} \sum_n x_n \mathrm e^{-\mathrm i n \Omega t} \\ \sum_n v_n \mathrm e^{-\mathrm i n \Omega t} \\ \sum_n c_n \mathrm e^{-\mathrm i n \Omega t} \end{array}\right). \end{eqnarray} \tag{ 6 }$$

The Langevin equation can be then transformed to the Fokker–Planck equation [39]. Specifically, equation (2) turns out to be

$$\begin{eqnarray} Q_n^{-} x_{n-1}+Q_n x_n+Q_n^{+} x_{n+1} & = &\delta_{n 0} D / \gamma, \end{eqnarray} \tag{ 7 }$$

in which

$$\begin{eqnarray} Q_n^{-} & = & -q \frac{\Omega^2}{4}+\frac{\mathrm{i} \Omega^3 q\left(n-1\right)}{4\left(2 \gamma-\mathrm{i} n \Omega\right)}, \end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray} Q_n& = & -\frac{n^2 \Omega^2}{2}+\frac{a \Omega^2}{4}-\mathrm{i} \frac{n \gamma \Omega}{2}-\mathrm{i} \frac{\Omega^3 n a}{4\left(2 \gamma-\mathrm{i} n \Omega\right)}, \end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray} Q_n^{+}& = & -q \frac{\Omega^2}{4}+\frac{\mathrm{i} \Omega^3 q\left(n+1\right)}{4\left(2 \gamma-\mathrm{i} n \Omega\right)}. \end{eqnarray} \tag{ 10 }$$

In order to further simplify equation (7), we introduce new variables $S^{\pm}_n$, which satisfy:

$$\begin{eqnarray} S_n^{+} & = & x_{n+1} / x_n,\quad \left(n\unicode{x2A7E} 0\right); \end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray} S_n^{-} & = & x_{n-1}/ x_n, \quad \left(n\unicode{x2A7D} 0\right). \end{eqnarray} \tag{ 12 }$$

Therefore, we can arrive at a recursion relation of $S^{\pm}_n$:

$$\begin{eqnarray} &&S_{n \mp 1}^{\pm} = -\frac{Q_n^{\mp}}{Q_n+Q_n^{\pm} {S_n^{\pm}}}, \end{eqnarray} \tag{ 13 }$$

and, in turn, the variance $\left\langle x^2 \right\rangle$, after the averaging over time, can be shown to be

$$\begin{eqnarray} \overline{\left\langle x^2\right\rangle} \equiv x_0 & = &\frac{4}{\Omega^2\left(a-2 q \operatorname{Re}\left(S_0^{+}\right)\right)-\frac{q \Omega^3}{\gamma} \operatorname{Im}\left(S_0^{+}\right)}\left(\frac{D}{\gamma}\right) \nonumber\\ & = &\frac{4 \, \overline{\left\langle v^2\right\rangle}}{\Omega^2\left(a-2 q \operatorname{Re}\left(S_0^{+}\right)\right)}, \end{eqnarray} \tag{ 14 }$$

where the time-averaged mean kinetic energy is given by

$$\begin{eqnarray} \overline{\left\langle v^2\right\rangle} \equiv v_0 & = &\frac{1}{1-\frac{\Omega}{\gamma} \frac{q \, \operatorname{Im}\left(S_0^{+}\right)}{\left(a-2q \operatorname{Re}\left(S_0^{+}\right)\right)}}\left(\frac{D}{\gamma}\right). \end{eqnarray} \tag{ 15 }$$

Now, the calculations of $\overline{\langle x^2\rangle}$ and $\overline{\langle v^2\rangle}$ can be carried out from equations (14) and (15), which requires us to know the value of $S_0^{+}$ beforehand. It is clear that, for given values of parameters a, q, and γ, $S_0^{+}$ should take a certain specific value since the values of $\overline{\langle x^2\rangle}$ and $\overline{\langle v^2\rangle}$ would otherwise be different for the same set of parameters according to equations (14) and (15). This will in turn lead to different results for the same physical situation of the ion motion. In other words, a fixed value of $S_0^{+}$ should be associated with the trap parameters.

Therefore, one can use the recursion equation (13) to evaluate $S_0^{+}$ with a wide range of parameters by using the numerical method. To do this, we choose an integer $n_{\max}$ as the times of the backward recursion and a random number $S_{n_{\max}}^{+}$ as the starting value of the recursion. By using equation (13), we can calculate the value of any value of $S_n^{+}$ for which $0 \unicode{x2A7D} n \lt n_{\max}$, including the particular $S_0^{+}$ that we aim at. It is worth noting that it is unnecessary to calculate the cases for n < 0 because the case of $n = -k$ is exactly equivalent due to $x_k = x_{-k}^*$ if one considers equation (13) for n = k. Hence we replace $S_n^{+}$ with S_n in the following text.

To validate this method, one should check whether it is possible to obtain the same value of S₀ when $n_{\max}$ and $S_{n_{\max}}$ are varied. In the recursion method, S₀ is the recursion ending value, while $n_{\max}$ and $S_{n_{\max}}$ is the times and the starting value for the recursion, respectively. The results are shown in figure 1, in which we have set $\Gamma = \frac{\gamma}{\Omega}$ as a dimensionless parameter for the damping coefficient. As can be seen, by varying $n_{\max}$ and $S_{n_{\max}}$, it is indeed possible to obtain the same value of S₀ when $n_{\max} \gt 3$ (clearly, values of $n_{\max}$ that are too small cannot be used). In practice, an intermediate value of $n_{\max} = 10$ can be chosen, since larger values are usually not necessary and take longer time to compute. Meanwhile, we find that the choice of $S_{n_{\max}}$ is not important, thus setting $S_{n_{\max}} = 1$ will be sufficient. It seems that different values of $S_{n_{\max}}$ corresponds to the same S₀, but we can find the recursion for n from smaller to bigger is rather unstable and only some particular values of S₀ can make the whole series converge.

Zoom In Zoom Out Reset image size

The above method can produce a converged value of S₀ under various parameters a, q, and Γ, but one has to make sure this method can indeed give us the right result. To check whether the S₀ calculated by this method is correct, for a given S₀, we can further use equations (14) and (15) to calculate $\overline{\langle x^2\rangle}$ and $\overline{\langle v^2\rangle}$ at a wide range of parameters and compare our results with those in literature.

In figure 2, we show the average of the kinetic energy as a function of q for various values of a. We observe curved upward-opening lines, which indicate a rapid divergence of the kinetic energy that can even result in negative values for certain parameters. Consequently, it can be deduced that the ion motion is unstable in during these circumstances. It is easy to find the parameter boundary which corresponds to a finite value by the numerical method. Within the boundaries, the ion motion can be stable and the Mathieu's equation has positive and finite eigenvalues, as shown in figure 3, which agrees with the results of the Mathieu's equation calculated by the conventional methods, such as numerically solving differential equations or using series expansions to solve the eigenvalue problem [40]. When γ increases, the range of the stable area also increases, which implies that the ion motion can be stable more easily. This is consistent with the physical picture. It should be noted that our proposed method is much faster and significantly more accurate than conventional methods [37].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The most convincing and direct way to verify our recursion method is to perform numerical simulations by the molecular dynamics and compare the results with those from equations (14) and (15). To do this, we first convert the equation of motion of the ion into a homogeneous system of first-order differential equations. We introduce the dimensionless parameters:

$$\begin{eqnarray} T & = &\frac{1}{2}\Omega t, \end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray} \Gamma & = & \frac{\gamma}{\Omega}, \end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray} D^{^{\prime}}& = &\frac{D}{l_0^2 \Omega}, \end{eqnarray} \tag{ 18 }$$

with l₀ being a constant with the dimension of length. There may be a multiplying constant in the noise, but we ignore it and only consider the relative value. We use the difference of the Wiener process to mathematically describe the white noise . With W_T representing the Wiener process [41], we have

$$\begin{eqnarray} \mathrm d v+2 \Gamma v \mathrm d T+\left(a-2 q \cos 2 T\right) x &\mathrm d T = \mathrm d W_T, \end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray} &\mathrm d x = v \mathrm d T. \end{eqnarray} \tag{ 20 }$$

We can use the fourth-order Runge–Kutta (R–K) method to integrate the differential equations. In practice, we compute the ion motion 1000 times with different random forces and add them up for the statistical average of x² and v². The implementation of the random force relies on the Python's random number selection. The force satisfies a normal distribution and applies to each ion independently (see appendix B). The results of $\langle x^2 \rangle$ are shown in figure 4.

Zoom In Zoom Out Reset image size

To verify our recursion method in the last subsection, we need to consider the limit of the long time and include both the statistical averaging and the time averaging. We need to choose a time range large enough for different parameters. The results are shown in figure 5, from which we can conclude that the results from equations (14) and (15) agree completely with those calculated by the numerical simulations. The small deviations from the numerical simulations are attributed to the fact that we cannot simulate an infinite time for the time averaging. These comparisons show that our recursion method to calculate S₀ is undoubtedly correct. It is noteworthy that $\langle x^2\rangle$ starts to increase when the damping parameter Γ is larger than a specific value around 0.6. This can be attributed to the fluctuation–dissipation theorem, which indicates that the fluctuation is proportional to the dissipation. It should be noted that the application of the WKB method below will give similar results in the case of a large Γ.

Zoom In Zoom Out Reset image size

The above discussions mainly focus on S₀ and the average of the variance. Actually, we can also verify equations (13) for other terms in the Fourier series of S_n , and the time evolution of x² can be also examined. When one calculates S₀ by the recursion method, at the same time we obtain S_n (obviously, for a given n, we should take a large enough $n_\textrm{{max}}$ to ensure a converged S_n ). Meanwhile, in the long time limit, $\langle x^2\rangle$ has the following Fourier decomposition,

$$\begin{equation} \langle x^2 \rangle = \sum_n \left( X_n \cos{2 n T}+X_n^{^{\prime}} \sin{2 n T} \right), \end{equation} \tag{ 21 }$$

where $X_n = x_n+x_n^* = 2x_{n}^{r}$ ($x_{n}^{r}$ is the real part of x_n ) due to the fact that $\langle x^2 \rangle$ is a real number (cf equation (6)). As an example, we verify the factors of the cosine terms. According to equations (11) and (12), we have

$$\begin{equation} X_n = 2x_{n}^{r} = 2{\mathrm{Re}}\left(x_0S_0S_1{\ldots}S_n\right). \end{equation} \tag{ 22 }$$

Following equation (21), one can make a discrete Fourier transform to $\langle x^2 \rangle$ in the long time limit:

$$\begin{align} X_n& = \frac{2}{T_0}\int_{t}^{t+T_0}\langle x^2 \rangle \cos\left(2 n T\right)\mathrm d T, \nonumber\\ & = \sum_{k = 1}^{T_0/\Delta T}\langle x^2 \rangle_T \cos\left(2 n \left(t+k\Delta T\right)\right) \Delta T. \end{align} \tag{ 23 }$$

As x₀ has been proven to be accurate and reliable, we can compare $2x_n/x_0$ using the recursion method and $X_n/x_0$ using the numerical method. The results are listed in table 1, which corresponds to the parameters given in the caption. Note that the errors arising from the fact that the discrete Fourier transform itself has a relative error of about $(\Delta t)^2\sim10^{-4}$ since we have chosen $\Delta t = 0.01$. As can be seen from table 1, the difference between the recursion method and the numerical simulation is around 10⁻⁴. One does not need to further verify the case for a bigger $n~(n\gt4)$ since $X_n/x_0$ will be smaller than this error limit.

Table 1. Fourier series of different orders for trap parameters $a = -0.2$ and q = 0.8, where two cases of the damping coefficient $\Gamma = 0.01$ and 0.1 have been considered. For each case, the first row is for the theoretical method and the second row for the numerical method a total time T = 1000 and a stepsize $\Delta T = 0.01$. The data are directly shown in the table with the same number of valid digits, regardless of the error.

Γ	n	0	1	2	3	4
	${2x_{n}^{r}}/{x_0}$	1.0000	−3.9000 $\times 10^{-1}$	5.6683 $\times 10^{-2}$	−4.1249 $\times 10^{-3}$	1.7897 $\times 10^{-4}$
0.01	${X_n}/{x_0}$	1.0000	−3.9003 $\times 10^{-1}$	5.6694 $\times 10^{-2}$	−4.0472 $\times 10^{-3}$	2.6688 $\times 10^{-4}$
	Δ	0	3 $\times 10^{-5}$	1 $\times 10^{-5}$	8 $\times 10^{-5}$	9 $\times 10^{-5}$
	${2x_{n}^{r}}/{x_0}$	1.0000	−3.8372 $\times 10^{-1}$	5.4872 $\times 10^{-2}$	−3.9330 $\times 10^{-3}$	1.6782 $\times 10^{-4}$
0.1	${X_n}/{x_0}$	1.0000	−3.8389 $\times 10^{-1}$	5.4880 $\times 10^{-2}$	−3.8444 $\times 10^{-3}$	2.6091 $\times 10^{-4}$
	Δ	0	17 $\times 10^{-5}$	1 $\times 10^{-5}$	9 $\times 10^{-5}$	10 $\times 10^{-5}$

As was recently proposed by Conangla et al [32], one can apply the WKB method to calculate the variance of the position x(t). The advantage of the WKB method is that one can examine the ion motion at any timescale, not restricted to the case of the long time limit. Furthermore, the analytical expressions can tell us how the ion motion will change in response to variations in the trap and noise parameters. In particular, we consider the overdamped situation where the dc term is neglected, i.e. one takes a = 0 in equation (2).

To simplify the subsequent derivation, we denote $T = \Omega t$, $\Gamma = \frac{\gamma}{\Omega}$, $D^{\prime} = \sigma^2$, and $\lambda_1 = -\frac{q^2}{8\Gamma^3}$. Following [32], at time T, one can obtain an approximate solution to equation (2):

$$\begin{equation} x\left(T\right) \approx \frac{\sigma}{\Gamma} \int_0^T \mathrm e^{\lambda_1\left(T-s\right)} \eta_s \mathrm d s = \frac{\sigma}{\Gamma} \int_0^T \mathrm e^{\lambda_1\left(T-s\right)} \mathrm d W_s. \end{equation} \tag{ 24 }$$

According to the Itô isometry [42], the variance of the approximate solution equals to

$$\begin{equation} E\left[x\left(T\right)^2\right] = \left(\frac{\sigma}{\Gamma}\right)^2 \int_0^T \mathrm e^{2 \lambda_1 \left(T-s\right)} \mathrm d s = \left(\frac{\sigma}{\Gamma}\right)^2 \frac{1-\mathrm e^{2 \lambda_1 T}}{2|\lambda_1|}, \end{equation} \tag{ 25 }$$

where E stands for the ensemble average, with the same meaning of the symbol '$\langle \rangle$'.

Now we can consider two limiting cases. Firstly, in the long time limit, one has

$$\begin{equation} E\left[x\left(T\right)^2\right] < \frac{\sigma^2}{2 \Gamma^2\left|\lambda_1\right|}\approx\frac{4\sigma^2 \Gamma}{m q^2} = Y_0, \end{equation} \tag{ 26 }$$

which indicates that the variance is a constant Y₀ independent on T. Secondly, in the limit of $T \rightarrow 0$, one gets

$$\begin{equation} E\left[x\left(T\right)^2\right] \approx \frac{\sigma^2}{\Gamma^2} T. \end{equation} \tag{ 27 }$$

According to equation (26), one immediately arrives at the conclusion that $E\left[x(T)^2\right]$ has an upper limit Y₀ which is proportional to the damping parameter Γ. This fact can be verified by the recursion method discussed in the previous subsections. In doing this, for the overdamped case where a = 0, we calculate $E\left[x(T)^2\right]$ at different Γ and for various values of q. The results are presented in figure 6, which clearly shows that, above a certain value of Γ, $E\left[x(T)^2\right]$ increases as Γ grows, which are consistent with our assumptions (overdamped). In particular, for all the curves, $E\left[x(T)^2\right]$ indeed shows a linear relationship with Γ when it is sufficiently large. At this point, we have further confirmed the correctness of our recursion method. In addition, we mention that the WKB method can be also used in the case of the colored noise, as ill be shown below.

Zoom In Zoom Out Reset image size

The colored noise is no longer Markovian. The type of noise chosen in the current work is a classical one with a Gaussian correlation function, i.e.

$$\begin{equation} E\left[X_t X_u\right] \approx \frac{\sigma^2}{2 |\lambda^{\prime}|} \mathrm e^{-\lambda^{\prime}|u-t|}, \end{equation} \tag{ 28 }$$

where $\lambda^{\prime}$ represents the rate of the mean reversion and σ is the noise intensity. In this case, the colored noise is equivalent to the OU-process with the stochastic differential equation (SDE) given by [34]

$$\begin{equation} \mathrm d X_t = -\lambda^{^{\prime}} X_t \mathrm dt+{\sigma} \mathrm d W_t, \end{equation} \tag{ 29 }$$

in which W_t stands for the Wiener process and thus $\mathrm dW_t$ means a white noise.

We can now consider the influence of using different types of noises. In the overdamped case, if we ignore the difference of the constant term F(t) in equation (2), the function is the same as that discussed in the previous section. Thus, we can continue to use the WKB method, with appropriate changes. Following [32], we can get an approximate solution to x(T) in the integral form:

$$\begin{equation} x\left(T\right) = \frac{1}{\Gamma} \int_0^T \exp \left[\lambda_1\left(T-s\right)\right] X_s \mathrm d s \\. \end{equation} \tag{ 30 }$$

To carry out the above integral, we make a slight change to equation (29), i.e.

$$\begin{equation} X_s \mathrm ds = \frac{1}{\lambda^{^{\prime}}}\left(-\mathrm d X_s+{\sigma} \mathrm d W_s\right), \end{equation} \tag{ 31 }$$

which can be inserted into equation (30), yielding

$$\begin{equation} \begin{split} x\left(T\right) & = \frac{1}{\Gamma \lambda^{^{\prime}}} \int_0^T \exp \left[\lambda_1\left(T-s\right)\right]\left(\sigma \mathrm d W_s-\mathrm d X_s\right), \\ & = \frac{1}{\Gamma \lambda^{^{\prime}}}\left\{\sigma \int_0^T \exp \left[\lambda_1\left(T-s\right)\right] \mathrm d W_s-\left.\exp \left[\lambda_1\left(T-s\right)\right] X_s\right|_0 ^T\right. \left.-\lambda_1 \int_0^T \exp \left[\lambda_1\left(T-s\right)\right] X_s \mathrm d s\right\}, \\ & = \frac{1}{\Gamma \lambda^{^{\prime}}}\left\{\sigma \int_0^T \exp \left[\lambda_1\left(T-s\right)\right] \mathrm d W_s-X_T-\lambda_1 \Gamma x\left(T\right)\right\},\\ \end{split} \end{equation} \tag{ 32 }$$

where we have assumed $X_0 = 0$ because the starting of the OU-process is zero. Since now we have an equation having x(T) on both sides, we can make transposition and solve the dynamics of the ion with parameters and the noise function as follows:

$$\begin{equation} x\left(T\right) = \frac{\lambda^{^{\prime}}}{\lambda^{^{\prime}}+\lambda_1} \cdot \frac{1}{\Gamma \lambda^{^{\prime}}}\left\{\sigma \int_0^T \exp \left[\lambda_1\left(T-s\right)\right] \mathrm d W_T+X_T\right\}. \end{equation} \tag{ 33 }$$

Then, we can calculate the variance $E\left[x(T)^2\right]$ by the statistical averaging:

$$\begin{align} E\left[x\left(T\right)^2\right] & = \frac{1}{\left(\lambda^{^{\prime}}+\lambda_1\right)^2} \frac{1}{\Gamma^2} E\left\{\left\{\sigma \int_0^T \exp \left[\lambda_1\left(T-s\right)\right] \mathrm d W_s-X_T\right\}^2\right\}, \nonumber\\ & = \frac{1}{\left(\lambda^{^{\prime}}+\lambda_1\right)^2} \frac{1}{\Gamma^2} \left\{E\left\{\sigma^2\left\{\int_0^T \exp \left[\lambda_1\left(T-s\right)\right] \mathrm d W_s\right\}^2\right\}\right. \nonumber\\ & \quad \left.-E\left\{2 \sigma \int_0^T \exp \left[\lambda_1\left(T-s\right)\right] \mathrm d W_s \cdot X_T\right\}+E\left[X_T^2\right] \right\},\nonumber\\ & = E_1\left(T\right)+E_2\left(T\right)+E_3\left(T\right). \end{align} \tag{ 34 }$$

in which there are three terms. As has been shown in the previous section, the first term can be calculated by the Itô isometry [42] to be

$$\begin{equation} E_1\left(T\right) = \frac{1}{\left(\lambda^{^{\prime}}+\lambda_1\right)^2} \frac{1}{\Gamma^2}\left(\sigma^2 \frac{1-\mathrm e^{2 \lambda_1 T}}{2 \mid \lambda_{1} |}\right). \end{equation} \tag{ 35 }$$

Remembering that $E[X_t] = 0$ is considered in this work, we can thus drop the second term $E_2(T)$. For the third term, we know that the variance of the OU-process is given by

$$\begin{equation} E\left[X_T^2\right] = \frac{\sigma^2}{2\lambda^{^{\prime}}}\left(1-\mathrm e^{-2\lambda^{\prime} T}\right), \end{equation} \tag{ 36 }$$

from which we consider the case when the noise becomes stable. In this case, since $\lambda^{\prime} \rightarrow +\infty$, one thus has $E[X_T^2] = \frac{\sigma^2}{2\lambda^{^{\prime}}}$. After substituting these results into equation (34), we finally arrive at

$$\begin{equation} E\left[x\left(T\right)^2\right] = \frac{1}{\left(\lambda^{^{\prime}}+\lambda_1\right)^2} \frac{1}{\Gamma^2}\left(\sigma^2 \frac{1-e^{2 \lambda_1 T}}{2 \mid \lambda_{1} |}+\frac{\sigma^2}{2 \lambda^{^{\prime}}}\right). \end{equation} \tag{ 37 }$$

Now, let us first consider the case of the long time limit, which also means $\left|\lambda_1\right| T \rightarrow +\infty$. Similarly, under these conditions, it is easy to show that $E\left[x(T)^2\right]$ can also reach a finite limit, i.e.

$$\begin{equation} \begin{split} E\left[x\left(T\right)^2\right]& = \frac{1}{\left(\lambda^{^{\prime}}+\lambda_1\right)^2} \frac{\sigma^2}{2 \Gamma^2} \frac{\lambda^{^{\prime}}+\lambda_1}{| \lambda_1 \mid \lambda^{^{\prime}}}, \\ & = \frac{1}{\left(\lambda^{^{\prime}}+\lambda_1\right) \lambda^{^{\prime}}} \frac{\sigma^2}{2 \Gamma^2\left|\lambda_1\right|}. \\ \end{split} \end{equation} \tag{ 38 }$$

We further assume $\lambda^{^{\prime}}T \gg |\lambda_1|T$, which means that the noise becomes stable much faster than the ion motion reaches the above limit. Thus, we finally get the following results:

$$\begin{equation} E\left[x\left(T\right)^2\right] = \frac{\sigma^2}{2 \Gamma^2\left|\lambda_1\right|} \cdot \frac{1}{\lambda^{^{\prime} 2}}\left[1-\frac{\lambda_1}{\lambda^{\prime}}+o\left(\frac{\lambda_1}{\lambda^{\prime}}\right)\right]. \end{equation} \tag{ 39 }$$

We approximate the analytical solution in order to show the relationship between the variance and the parameters distinctly, which can be verified by numerical simulations. In principle, we can still use the R–K method to integrate the motion of the ion, but we first need to generate the colored noise. We calculate the value of the colored noise by solving the SDE for the OU-process according to equation (29). Considering that the function itself contains the white noise, one judges that the error from the algorithm of the numerical solution to the SDE has little influence on our results. Hence, we choose to use the forward Euler method:

$$\begin{equation} X_{T+\Delta T} = X_{T}-\lambda^{^{\prime}}X_T \Delta T+\sigma \eta\left(T\right)\sqrt{\Delta T}, \end{equation} \tag{ 40 }$$

where $\eta(T)$ is the white noise. For the convenience of calculations, we assume that $\eta(T)$ is independent of time and subject to a standard normal distribution. The stepsize $\Delta T$ is chosen to be the same as that we use to calculate the motion of the ion.

In figure 7, we demonstrate the value of the colored noise under two different values of rate of the mean reversion $\lambda^{\prime}$. One finds that for a larger $\lambda^{\prime}$, as shown in figure 7(a), the volume $\lambda^{\prime} T$ quickly rises, which means that the noise becomes to be stable so fast that we cannot observe any apparent change of the process over time. The behavior is very similar to the case of the white noise. On the contrary, when $\lambda^{\prime} T$ becomes small, as shown in figure 7(b), it resembles the Wiener process, as expected.

Zoom In Zoom Out Reset image size

Now, we consider the case where $\lambda^{\prime}$ is much larger than $\frac{1}{T}$. We first calculate the noise X_T and supply it to the R–K integration algorithm to the following equations:

$$\begin{eqnarray} \mathrm d v+2 \Gamma v \mathrm d T-2q \cos\left( 2 T \right)x &\mathrm d T = 2X_T \mathrm dT, \end{eqnarray} \tag{ 41 }$$

$$\begin{eqnarray} &\mathrm d x = v \mathrm d T. \end{eqnarray} \tag{ 42 }$$

In figure 8(a), we show our result from the numerical calculation, which appears similar to the case of the white noise. This is not surprising, as can be seen from equation (39), one has a factor of $\frac{\sigma^2}{2 \Gamma^2\left|\lambda_1\right|}$, which is the same as the case of the white noise. The second factor $\frac{1}{\lambda^{^{\prime} 2}}$ only contributes a multiple constant. The last factor $\left[1-\frac{\lambda_1}{\lambda}+o\left(\frac{\lambda_1}{\lambda}\right)\right]$ approaches to 1 and thus does not largely change our result when $\lambda^{^{\prime}} \gg |\lambda_1|$. Overall, the result should resemble the case of the white noise and may only increase slightly due to the last term since $\lambda_1\lt0$.

Zoom In Zoom Out Reset image size

In figure 8(b), we compare the results under various values of $\lambda^{\prime}$. As can be seen, one observes that the average of the variance of the ion position gradually becomes smaller when one increases the value of $\lambda^{\prime}$. To verify the effect of the factor $\frac{1}{\lambda^{^{\prime} 2}}$ in equation (39) for the case of $\lambda^{^{\prime}}\gg |\lambda_1|$, we calculate $\langle x^2 \rangle \times \lambda^{^{\prime} 2}$ as a function of time. As shown in figure 8(c), we find that the four curves tend to overlap with each other as $T\rightarrow \infty$, which confirms that our results for $\frac{1}{\lambda^{^{\prime} 2}}$ are reliable.

Then, we move to examine the motion of the ion under different damping strengths by changing the damping coefficient Γ. We have chosen four different values of Γ around 6.0 but with small differences, because the computation time for a converged result increases rapidly as Γ becomes large due to the fact that the volume is $\lambda_1 T$ with $\lambda_1 = -\frac{2q^2}{\Gamma^3}$. The results are shown in figure 8(d), from which we find that the four curves converge to similar values when $T\rightarrow \infty$. Again, please note that the time needed for the convergence increases with the increase of Γ. These observations are similar to those made for the case of the white noise.

Finally, we consider the case of $T\rightarrow 0$, and it is unsuitable to use equation (37) to make approximations because equation (36) may result in entirely different outcomes. Instead, we are able to simulate this case by choosing a smaller $\lambda^{\prime}$, and we find that the noise appears to be distinct, as shown in figure 7(b). Then, we calculate the variance of the position with $\lambda^{\prime} = 0.01$, and the result is shown in figure 9. Significantly, the behavior of the outcome is markedly different from those previously discussed, as different tendencies are observed in three distinct regions. Specifically, when $T\rightarrow \infty$, the curve indeed gradually converges, because the noise becomes stable so that the previous approximation remains valid. In the case of $T\rightarrow 0$, the variance $\langle x^2 \rangle$ increases nonlinearly, while it is approximately proportional to the time when T is in an intermediate range.

Zoom In Zoom Out Reset image size

Now, we will discuss the limit case of $\lambda^{\prime} \rightarrow 0$ for the OU-process, in which case figure 7(b) is found to be very similar to the Wiener process. Hence, we can make approximation to equation (26) and the noise becomes a Wiener process. In particular, we denote it as W_t [41]:

$$\begin{equation} \mathrm d X_t = {\sigma} \mathrm d W_t. \end{equation} \tag{ 43 }$$

We have used the inverse of the Wiener process to mathematically describe the white noise. Therefore, the Wiener process has a higher order of continuity than the white noise, which will simulate the case where the noise is disturbed and keeps being accumulated.

As the Wiener process is the limit case of the OU-process, we can find the ion motion using the method in the preceding subsection. Specifically, we use the WKB method and make similar procedures as those to equation (30). By integrating by parts, we have

$$\begin{equation} \begin{split} x\left(T\right) & = \frac{\sigma}{\Gamma} \int_0^T \mathrm e^{\lambda_1\left(T-s\right)} W_s \mathrm d s, \\ & = \frac{\sigma}{\Gamma \lambda_1} \int_0^T\left(\mathrm e^{\lambda_1\left(T-s\right)}-1\right) \mathrm d W_s. \end{split} \end{equation} \tag{ 44 }$$

Again, we use the Itô isometry [42] to calculate the variance of the position of the ion:

$$\begin{equation} \begin{aligned} E\left[x\left(T\right)^2\right] & = \left(\frac{\sigma}{\Gamma \lambda_1}\right)^2 \int_0^T\left(\mathrm e^{\lambda_1\left(T-s\right)}-1\right)^2 \mathrm d W_s, \\ & = \left(\frac{\sigma}{\Gamma \lambda_1}\right)^2\left(T+2 \frac{1-\mathrm e^{\lambda_1 T}}{\lambda_1}-\frac{1-\mathrm e^{2 \lambda_1 T}}{2 \lambda_1}\right). \end{aligned} \end{equation} \tag{ 45 }$$

Apparently, by taking different limit of T, one can arrive at the following results:

$$\begin{equation} \begin{aligned} & E\left[x\left(T\right)^2\right] \approx \frac{1}{3}\left(\frac{\sigma}{\Gamma}\right)^2 T^3,\,\, T \rightarrow 0; \\ & E\left[x\left(T\right)^2\right] \approx\left(\frac{\sigma}{\Gamma \lambda_1}\right)^2 T, \,\, T \rightarrow \infty. \end{aligned} \end{equation} \tag{ 46 }$$

For the numerical simulation, we can calculate the Wiener process as we have done in the previous subsection. According to the SDE of the Wiener process, we can use the forward Euler method and arrive at

$$\begin{equation} X_{T+\Delta T} = X_{T}+\sigma \eta\left(T\right)\sqrt{\Delta T}. \end{equation} \tag{ 47 }$$

The result is shown in figure 10, which is very similar to that in figure 7(b). Then, using the noise, we can calculate the motion of the ion by the R–K method for the following equations:

$$\begin{eqnarray} \mathrm d v+2 \Gamma v \mathrm d T-2q \cos\left( 2 T \right)x &\mathrm d T = 2W_T \mathrm dT, \end{eqnarray} \tag{ 48 }$$

$$\begin{eqnarray} &\mathrm d x = v \mathrm d T. \end{eqnarray} \tag{ 49 }$$

Zoom In Zoom Out Reset image size

The result of $E\left[x(T)^2\right]$ is shown in figure 11, from which we find that the curve is partially similar to the front part of figure 9. It is easy to see that $E\left[x(T)^2\right]$ is proportional to the time when T is sufficiently large, but it endlessly increases over time and has no limit.

Zoom In Zoom Out Reset image size

For $T \rightarrow 0$, $E\left[x(T)^2\right]$ increases very slowly and the curve appears to be flat. In order to examine the exponent over T when $T \rightarrow 0$, we calculate $\ln{\langle x^2 \rangle}$ as a function of $\ln{T}$. As shown in the insert of figure 12, the natural logarithm of $\langle x^2 \rangle$ is proportional to $\ln{T}$ when T is larger than 1 ($\ln{T}\gt0$). This confirms that $\langle x^2 \rangle$ has a power law relation with T. To be more visually clear, we can calculate the slope of the line. As can be seen from the main figure of figure 12, it is evident that the value of the slope is only slightly larger than 3, which is apparently due to the effect of higher power terms.

Zoom In Zoom Out Reset image size

With the above discussions in mind, we can come back to the results of the OU-process as shown in figure 9. Obviously, the first and the second part of figure 9 agree with the case of the Wiener process. The first part requires $\lvert \lambda_1 \rvert T\ll 1$ and $\lambda^{\prime} T \ll 1$. The variance of x increases slowly and follows the power law. The second part requires $\lvert \lambda_1 \rvert T\gg 1$ and $\lambda^{\prime} T \ll 1$ when the noise has not been stable and the variance of x is proportional to the time. Finally, the third part requires $\lambda^{\prime} T \gg 1$ when the noise has been stable and the result is the same as the white noise and the colored noise of the OU-process when $\lambda^{\prime} \gg \lvert \lambda_1 \rvert$. If $\lambda^{\prime} \sim \lvert \lambda_1 \rvert$, the second part may not be distinguishable. The first part will be hard to be identified if $\lambda^{\prime} T \gg 1$.

In summary, we have investigated the motion of a single trapped ion subject to a colored noise environment by using an analytical approximation technique. Our theoretical results are consistent with numerical simulations, which indicates the accuracy of our approximation. The approach presented here can be extended to study other systems influenced by the colored noise.