Molecular Dynamics Simulation

The method of molecular dynamics simulations (MD) is one of the routine tools in the theoretical study of a many-body system. This computational method gives direct numerical solutions of the equations of motion of particles and calculates their time dependent behaviors, which allow us to check the models, offer insights to the experimentalists, assist in the interpretation of new results, predict the properties of the dynamics and other complicated phenomena that cannot be found out in other ways. At the same time, they reveal hidden details behind the experimental measurements. It acts as a bridge between microscopic length and time scales and the macroscopic world of the laboratory, and between theory and experiment.

Equations of motion

Molecular dynamics simulations are based on solving the classical equations of motion, which for trapped ions in a linear Paul trap may be written as  

where mi and ri are mass and position vector of the ith particle, respectively.

is the total force acting on the ith particle, which may depend on the positions and velocities of all other particles, and the time t. There are several contributions to this force.


The first term Fitrap is the force from the trapping potential and can be written as


For the rf electric field it can be weell approximated as a time averaged harmonic pseudopotential, The motion of the particle is separated into slow and fast parts (called macromotion and micromotion, respectively).

The second term FiCoulomb is due to the Coulomb repulsive force between each pair of particle in the trap.

The stochastic force Fistochastic results from the random interactions of particles with surroundings, such as collisions with residual gas, the scattering light, electric field noise, and so forth. In simulation all the factors are modelled in a simple way by velocity “kicks” applied to each particle k of species i and at each integration step.

The last term Filaser acts only on laser cooled particles. This force can be simplified to be a constant light pressure force independent of the ion velocity plus a linear viscous damping force FL=-bdx/dt with friction coefficients in the range b=(1.2-8)*10-22 kg/s.

Ion number determination:

In the simulation, the ion number can be adjusted arbitrarily, and the parameters of external electric fields can be modified near the experimental values conveniently, so we see different shape and structure of the laser cooled ion ensemble in real-time, and select the best set to fit the CCD image of experiment. Since the ion number of laser cooled ions is determined, the counting coefficient of a channeltron ion detector in the experiment is calibrated. Now the ion number of sympathetically cooled complex molecular ions can also be determined.

Fig. 1 is an example: Through the comparison of simulation and experiment we know there are 27 barium ions and 3 isotopes. While it is easy to count the number of barium ions in the case, it is impossible to count them in Coulomb crystal consisting of several hundred or thousand ions, there you rely on the simulations..



Fig. 1. Compare the simulation to the experiment to determine the number of ions,
through the comparison we know there are 27 138Ba+ and 3 invisible isotopes

Temperature determination:

Another important feature of the simulations is the determination of the final temperatures of different ion ensembles. Therefore, we simulate the pure laser cooled ion ensemble with different temperatures, and compare them to the CCD image to know which temperature the experiment is close to. Fig. 2 shows a simulated pure barium ion ensemble (300 ions) with different temperatures. Second, we keep the parameters set in the first step such as cooling and heating rate of the barium ions. Then we add the sympathetically cooled complex molecular ions, which are hotter and will heat the barium ion ensemble. By adjusting the heating rate of the complex molecular ions we can get the barium ion ensemble images with different final temperatures. Finally, we compare the simulated images with the CCD image of barium ion ensemble (with the complex molecular ions in the trap simultaneously) and get the right heating rate for the complex molecular ions. Since we know all the parameters, the final temperature of the complex molecular ion ensemble is determined. Fig. 3 is a good example.






Fig. 2 Simulated ion ensemble with different temperatures



CCD image of a pure barium ion crystal, through this we find the heating rate of barium ions



CCD image of barium ion ensemble with AF350 ions simultaneously, through this we find the heating rate for the complex molecular ions.

Fig. 3. Simulation of a three-species Coulomb crystal containing 830 laser-cooled 138Ba+ ions (blue) at a temperature of 25 mK, 420 sympathetically cooled ions of barium isotopes (red) and 200 protonated AlexaFluor350 molecules (AF350, green) at 120 mK