I am attempting to re-create a Mathieu stability diagram like the one shown in a paper by Leary and Schmidt [1]:

Mathieu stability diagram for a quadrupole mass filter

I expected that I could use MathieuC to generate this graph by assuming that instability occurs when the function returns a complex number, so I try

Quiet@DiscretePlot[Re@a /. FindRoot[MathieuC[a, q, 1], {a, 0.2}], {q, 0, 1, 0.1}]

which is wrong for two reasons - it produces a plot not even close to the desired output and the output is dependent on what value I use for z (in the example above, 1).

For those interested, the impetus behind this problem is to explore how Mathematica can be used to simulate the behavior of a quadrupole mass spectrometer.


  1. Leary, J. J.; Schmidt, R. L. Quadrupole Mass Spectrometers: An Intuitive Look at the Math. J. Chem. Educ. 1996, 73 (12), 1142. DOI: 10.1021/ed073p1142.
  • 1
    $\begingroup$ I did this for my Ph.D. and would be happy to send you the code. Shall we talk about it in the chat? I think the question is TL otherwise. Your key reference is V. I. Baranov, J. Am. Soc. Mass Spectrom. 14, 818 (2003). Here's what it's supposed to look like. The contours are indicative of the oscillation amplitude. My icon is a picture of the QMF and one of its trajectories, by the way. $\endgroup$ Commented Apr 25, 2014 at 5:26
  • $\begingroup$ @OleksandrR. I'd be very interested in chatting about this further; let me read through the paper you suggested first. My intentions here are primarily pedagogical and my goal is to expose undergrads to MSD design in an interactive manner. $\endgroup$ Commented Apr 25, 2014 at 13:31
  • $\begingroup$ The paper I suggested is quite technical and not really a good pedagogical introduction to the QMF. There actually is not a lot of useful literature about it, except for the original papers: W. Paul und H. S. Steinwedel, Z. Naturforschung 8a, 448 (1953); W. Paul und M. Raether, Z. Physik 140, 262 (1955). These are both in German but are better than anything I could find in English. The definitive reference (in English) is Quadrupole Mass Spectrometry and Its Applications by Peter Dawson, but this book is expensive and I wasn't able to find any copies of it at the time. $\endgroup$ Commented Apr 25, 2014 at 14:25
  • $\begingroup$ can someone provide some (physical, mathematical) interpretations on the three graphs illustrating the three stable regions? It could be also in terms of the stability of the Mathieu diagrams, you know. i am profoundly interested in this. Thank you very much! $\endgroup$
    – April
    Commented Jun 27, 2017 at 14:35
  • $\begingroup$ I am interested in how one can interpret the three graphs illustrating the three regions - first/conventional, second and third/intermediate region. What are these closer views on these three regions are illustrating? Thank you very much! $\endgroup$
    – user51078
    Commented Jul 12, 2017 at 15:00

1 Answer 1


I think that this question is too localized as it concerns the physics of a specific scientific instrument. Nonetheless, it is upvoted, so here I provide an answer for the benefit of the voters. I would still be happy to discuss this in the chat.

The mathematics of the quadrupole mass filter is more complicated than you might think. Basically, your assumption is not correct; the stability condition is actually

$$-C(a,q,\xi) \dot{S}(a,q,\xi) \in \mathbb{R}$$

which in fact does not depend on the formal parameter $\xi$ (which can be thought of as a sort of dimensionless time variable defining the phase of the oscillating electric field). The oscillation amplitude is proportional to the reciprocal of this quantity.

Additionally, the QMF itself has two planes of symmetry, $a=0$ and $q=0$, and we require stability in all four quadrants. However, the solutions of the Mathieu equation are symmetric about $a=0$ anyway, and we can get away in practice with placing just one mirror plane along the diagonal. So, with

$$P(a,q)=-C(a,q,\xi=0) \dot{S}(a,q,\xi=0)$$

we can plot the (symmetrized) stability diagram simply from the contours of $P(a,q) P(-a,-q)$ as follows:

p[a_, q_] := -MathieuC[a, q, 0] MathieuSPrime[a, q, 0]
 p[a, q] p[-a, -q], {q, 0.05, 0.95}, {a, 0.00, 0.25}, 
 MaxRecursion -> 3, RegionFunction -> Function[{x, y, f}, f > 0], 
 ColorFunction -> (ColorData["DarkRainbow"][1 - #] &),
 AspectRatio -> 1/GoldenRatio

This takes a while to run and produces a large output, which looks as follows:

plot of the stability parameter in the first Mathieu stability region

Although quite aesthetically pleasing and perhaps useful as a visual guide to the variation of the oscillation amplitude in $(a,q)$ space, this is not a very practical way to plot the diagram, and perhaps more importantly, it does not really provide any mathematical insight.

A better approach is to plot our diagram as the domain over which the elliptic sine and cosine functions

$$ce_r(\xi,q)=C(a_r,q,\xi) \\ se_r(\xi,q)=S(b_r,q,\xi)$$

are periodic. This is true for characteristic values $a_r(r,q)=$ MathieuCharacteristicA[r, q] and $b_r(r,q)=$ MathieuCharacteristicB[r, q]. $r$ is an index that numbers the elliptic trigonometric functions according to the number of roots they possess in the interval $0\le\xi\le\pi$. The elliptic sine is an odd function, so there is no $se_0$; the first periodic elliptic sine is thus $se_1$.

In other words,

 { MathieuCharacteristicA[0, q],  MathieuCharacteristicB[1, q] (* upright *),
  -MathieuCharacteristicA[0, q], -MathieuCharacteristicB[1, q] (* reflected*)}, 
 {q, 0, 1}, PlotRange -> {All, {0.0, 0.3}},
 PlotStyle -> {
   Directive[Thick, Blue], Directive[Thick, Red],
   Directive[Thick, Dashed, Blue], Directive[Thick, Dashed, Red]
 Filling -> Table[{2 n + 1 -> {{2 n + 2}, Directive[Opacity[1/2], Purple]}}, {n, 0, 1}]

which looks like this:

plot of the boundaries of the first Mathieu stability region

Its apex is at:

FindRoot[MathieuCharacteristicA[0, q] + MathieuCharacteristicB[1, q], {q, 0.6}]
(* -> {q -> 0.7059960697133709`} *)

{a -> MathieuCharacteristicB[1, q]} /. %
(* -> {a -> 0.23699399311247354`} *)

What was the point of all that mathematical discussion just to plot this, you might ask? Well, actually, there are many stable regions located apart from another in $(a,q)$ space, any of which might be used in a practical instrument. Infinitely many, in fact. Here are the three for which instruments have been constructed so far:

   {Tooltip[ MathieuCharacteristicA[r,     q],  Subscript[ce, r    ]], 
    Tooltip[ MathieuCharacteristicB[r + 1, q],  Subscript[se, r + 1]],
    Tooltip[-MathieuCharacteristicA[r,     q], -Subscript[ce, r    ]], 
    Tooltip[-MathieuCharacteristicB[r + 1, q], -Subscript[se, r + 1]]},
   {r, 0, 1}
  ], {q, 0, 8}, Evaluated -> True,
 PlotRange -> {All, {0, 4}},
 PlotStyle -> {
   Directive[Thick, Blue], Directive[Thick, Red],
   Directive[Thick, Dashed, Blue], Directive[Thick, Dashed, Red]
 Filling -> Table[{2 n + 1 -> {{2 n + 2}, Directive[Opacity[1/2], Purple]}}, {n, 0, 3}]

plot of the boundaries of the first three Mathieu stability regions

Here's a closer look at them.

The first/conventional region:

first Mathieu stability region

The second/"r.f.-only" region:

second Mathieu stability region

The third/"intermediate" region:

third Mathieu stability region

We can also plot the ion trajectories anywhere in the $(a,q)$ plane if we wish, using the expressions given in V. I. Baranov, J. Am. Soc. Mass Spectrom. 14, 818 (2003). (Beware! This paper contains some mistakes.) But the expressions are complicated and that is beyond the scope of this post.

  • 1
    $\begingroup$ Brilliant and descriptive. While I agree that the question is TL, one cannot argue the ease with which relatively simply Mathematica functions can be used to describe a difficult system. $\endgroup$ Commented Apr 25, 2014 at 14:36
  • 1
    $\begingroup$ @bobthechemist yes, the special function implementation in Mathematica is really excellent. These functions are far from easy to evaluate and the fact that Mathematica can do it reliably and to high precision is a testament to a lot of hard work that went into programming them. $\endgroup$ Commented Apr 25, 2014 at 14:41
  • 1
    $\begingroup$ That's an excellent explanation. Separation in elliptic coordinates is already quite a bit trickier than polar coords... $\endgroup$
    – Jens
    Commented Apr 25, 2014 at 16:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.