Simulating a quadrupole

Before actually constructing a DIY-MS, extensive simulations should be done. This is to assure that at least the mathematical side of the project is correctly worked out and no(more) time and money is wasted on something that could not work in the first place. A first step would therefor be the simulation of the quadrupole mass filter.

Luckily, the ion-movement in the quad follows basic physical principles and (like all the good stuff in engineering) can be described by a differential equation, namely Mathieus linear ordinary second order differential equation. Solutions of this ODE describe ion movements in the quad, which can be seen to perform either stable oscillations or quickly "shoot" to infinity.

With calculus programs like Matlab or GNU Octave, these ion movements (for an ideal quadrupole) can be numerically calculated and plotted. I set up a quick Octave script to plot some trajectories for stable and unstable ions, like this stable ion trajectory for a stable ion of m/z 28:

and this unstable ion of m/z 29 under otherwise same conditions:

Interestingly but not surprisingly, even for an unstable ion, either the X-Z or Y-Z plane is still stable. Explanation can be found when looking at the aq stability diagram, provided for example by this publication by Ma & Taylor. Basically, for small a and q either x or y are stable and the stability diagram is the overlapping region of these stable zones.

Playing around with this small simulation additionally showed me how fragile these stable trajectories actually are. Even a RF-frequency offset of as little as 1% caused the otherwise stable ion of m/z 28 in this simulation to collide. This raises some concern about the electronics-side and overall tolerances of this project, wich have to be designed and manufactured to absolute precision as it seems. At least if a useful resolution of any kind ist to be achieved.

Anyway, this simulation is the first step to simulating an ESI-Quad mass spectrometer and a nice first step towards the ultimate goal.

For further reading I suggest the linked article by Ma & Taylor. The Skript for the simulation can be found below. Octave is free to use software and in many cases very similar to Matlab.

Cheers
Marco

# This skript for GNU Octave simulates ion trajectorys through ideal electric quadrupoles
# using an numeric approach to Mathieus ODE. All calculations are performed using SI-units.
# The quadrupole and initial values can be set up through the val-variable. 
# The skript is an early proof-of-concept an may contain major errors. 
# Initial values are taken from S. Taylor "Simulation of ion trajectories through the mass filter
# of a quadrupole mass spectrometer"
 
 clear
 graphics_toolkit ('gnuplot');
 lsode_options("integration method", "stiff"); 
 
 global val e ukg vz tmax phasepi ax qx
 
 val = [ 0028    # 1 m/z  Ion-specific m/z
   0020    # 2 U  Offset voltage
   0123.5    # 3 V  RF voltage
   2*pi*2.02e+6       # 4 omega RF frequenzy f ... omega=2*pi*f
   0.00275   # 5 r0  Quad formfactor: inner radius
   1   # 6 U beschl. Accelerationvoltage for the ion
   10e-9   # 7 dt  Timeinterval for simulation
   0.15   # 8 zmax Quad formfactor: length in z
   0.0003   # 9 x0   Entry offset in x for ion when entering the quad
   0.0000   #10 vx0  Initial velocity in x  
   0.0003   #11 y0   Entry offset in y for ion when entering the quad
   0.0000   #12 vy0  Initial velocity in y 
   pi/4];   #13 phase Ion entry phase
 
 e=1.602176462e-19;     #Electroncharge
 ukg=1.66054e-27;     #1u=1*ukg kg
 
 vz=sqrt(1/val(1)*2*val(6)*e/ukg);   #Velocity in z
 tmax=val(8)/vz;      #Time for flight through
  
 phasepi=val(13);
 
 ax=1/val(1)*e/ukg*4*val(2)/val(5)^2/val(4)^2#ax-Value in stability diagram
 qx=1/val(1)*e/ukg*2*val(3)/val(5)^2/val(4)^2#qx-Value in stability diagram
 ax/qx       #show a/q 

 function xdot=pendx(x,t)    #Function describing movement in x
  global ax qx phasepi
  xdot(1)=x(2);
  xdot(2)=-(ax+2*qx*cos(2*t+phasepi))*x(1);
 endfunction
 
 function xdot=pendy(x,t)    #Function describing movement in y
  global ax qx phasepi
  xdot(1)=x(2);
  xdot(2)=-(-ax-2*qx*cos(2*t+phasepi))*x(1);
 endfunction
 
 solx=lsode("pendx",[val(9) val(10)],t=[0:val(4)*val(7)/2:val(4)*tmax/2]); #Solve x movement
 soly=lsode("pendy",[val(11) val(12)],t=[0:val(4)*val(7)/2:val(4)*tmax/2]); #Solve y movement

 Z=[0:val(8)/(tmax/val(7)):val(8)];   #Calculate a z for every tau
 Points=tmax/val(7)      
 Zres=size(Z)
 Tres=size(t) 
 
 plot(Z,solx(:,1));     #Plot the Graph
 ylim([-val(5) val(5)]);
 xlim([0 val(8)]);
 title ('Ion Trajectory in X-Z-Axis');