Frequency-selective harmonic retrieval for Schottky mass spectrometry

Nuclear mass measurements by means of Schottky mass spectrometry critically rely on an accurate determination of revolution frequencies of the circulating ions in a storage ring. Such a harmonic retrieval problem is conventionally tackled via the periodogram of the Schottky data, where the ion peaks are identiﬁed and their spectral locations are obtained by ﬁttings. However, the discrete frequency grid of the periodogram has unfortunately hampered a ﬁne resolution of two closely spaced harmonics. We thereby propose a method based on the state space representation in the frequency domain to overcome this limit. Moreover, its frequency-selective merit has allowed the method to focus only on a narrow band and thus greatly reduced the computational cost while still retaining superb accuracy. With the real Schottky data from an isochronous-Schottky beam time at the experimental cooler-storage ring in Lanzhou, the accuracy of the retrieved harmonics is demonstrated to be around 1 ppm, as limited by the anisochronism effect of the ion optics.


I. INTRODUCTION
Schottky mass spectrometry (SMS), along with isochronous mass spectrometry, is a precision ring-based spectroscopic technique for direct nuclear mass measurements [1,2].The merits of large ring acceptance and small beam emittance due to the electron cooling have allowed SMS to precisely map a sizable area of nuclides on the nuclear chart during one beam time (see, e.g., Refs.[3,4]).To date, about 300 mass values of nuclei in ground and isomeric states have been obtained by SMS [1].Recently, SMS was successfully performed in combination with the isochronously tuned ion optics of the experimental cooler-storage ring (CSRe) at the Institute of Modern Physics (IMP) [5].Because the electron cooling is no longer required, this technical advancement has potentially pushed the lower lifetime limit of measurable exotic nuclei downwards into the millisecond regime.The measurement precision of the isochronous SMS will be comparable to the cooling SMS after employing a position-resolving Schottky resonator to correct for the hampering anisochronism effect [6,7].
Either way, an accurate determination of revolution frequencies of the stored ions has always been a central focus in SMS data analysis since they reflect the mass relationships between different ion species [8].It is well known that the power spectral density (PSD) of an ion's Schottky signal in the frequency domain manifests an infinite number of spikes equally distanced by its revolution frequency, where each spike is called a harmonic and indexed by its harmonic number [9].Hence, the revolution frequency can be obtained from any harmonic by dividing its spectral location by the corresponding harmonic number.Within this context, the present work addresses the harmonic retrieval problem, which aims to accurately determine the frequency of one specific harmonic.
The widest adopted approach to the harmonic retrieval problem is to directly identify the spectral peak on the PSD spectrum and extract its location by Gaussian fitting (see, e.g., Refs.[3,4]).Usually, the PSD spectrum is estimated by the periodogram, which is the modular square of the windowed discrete Fourier transform (DFT) of the Schottky data [10].By varying the window function, one can find a balance between the spectral leakage and the peak width, depending on whether the particular purpose is detecting a weak peak or resolving adjacent peaks [11].In practice, it is also necessary to average several similar PSD spectra to further reduce the noise fluctuation, which is beneficial in certain cases with a poor signal-to-noise (SN) ratio, such as single ion detection [12].
While this nonparametric method is quite intuitive, it critically depends on the PSD's estimate, which changes substantially across different window functions.What is even worse, the method considers only the DFT's modulus and discards its argument, which may lead to an unsuccessful discrimination of two closely spaced harmonics (see Sec. III for details).As a complementary approach, a few ab initio methods, such as the multiple signal classification (MUSIC) [13] and the estimation of signal parameters via rotational invariance techniques (ESPRIT) [14], have been proposed to overcome the aforementioned limitations.They are all based on the state space representation of the data: a mathematical technique that was first developed by R. Kalman [15] and has become a canonical method in control engineering ever since [16].Those methods can yield adequate results by exploiting the data's internal structure but meanwhile suffer from the heavy computational burden when calculating the data's autocorrelation in the time domain, as well as from the intraharmonic interference if the data contain too many harmonics.
It was then shown independently by two separate groups that a state-space-based method could likewise be developed in the frequency domain [17,18].As a result, the computational cost can be greatly reduced by means of the fast Fourier transform (FFT) algorithm, and the intraharmonic interference can be well controlled due to the harmonics' frequency-isolating property.Having been inspired by those frequency-domain methods, we present in this work an improved harmonic retrieval solution with a favorable frequency-selective merit as an alternative approach to the determination of revolution frequencies in SMS data analysis.

II. METHOD
Owing to the periodic revolutions of the stored ions, the digitized Schottky signal y n of length N can be modeled as a superposition of J sinusoids: where a j = α j e iθ j and ω j ∈ [−π, π ) are the complex amplitude and the angular frequency of the jth sinusoid, respectively.The physical revolution frequency f r is related to the angular frequency ω by where f s is the sampling rate of the Schottky data, f l is the center frequency set by the local oscillator, and h is the harmonic number.It is reasonable to require that all {ω j } J j=1 are distinct.The harmonic retrieval is, in fact, to estimate every ω j from {y n } N−1 n=0 .

A. State space representation
In light of the state space representation, the Schottky data in Eq. ( 1) can be rewritten in a recursive form: where x n is the state variable with the initial condition and The state space representation of Eq. ( 1) is actually not unique.For instance, with a nonsingular matrix T ∈ C J×J , we construct The trio (b , x 0 , A ) is also a valid representation.It can be verified that A's range and eigenvalues are invariant under this transformation [19].In other words, they are the intrinsic properties of the Schottky data.Indeed, the dimension of A's range, denoted by rank(A), is just the number of sinusoids, and the arguments of A's eigenvalues are the corresponding angular frequencies.
Hereinafter, we will switch from the time domain to the frequency domain by means of DFT: where z N = e i2π/N is the Nth primitive root of unity.Note that the set { 2πk N } N−1 k=0 defines the DFT grid.By virtue of DFT's linear property, Eq. ( 3) is transformed to which merely says that the frequency spectrum ỹk is a superposition of subspectra xk , of which each element is one sinusoidal spectrum at the kth grid point.Also, by virtue of DFT's shift property, Eq. ( 4) is transformed to where an extra term comprising s = (I − A N )x 0 accounts for DFT's implicit periodic extension of {x n } N−1 n=0 .Equations ( 6) and ( 7) together constitute a state space representation of the Schottky data in the frequency domain.

B. Principal matrix equation
By recursively substituting Eq. ( 6) into Eq.(7), one can establish a set of equations for the geometric progression {z lk N ỹk } L l=1 in terms of xk and s, expressed as where C ∈ C L×J is a Vandermonde matrix and S ∈ C L×L encloses all the s-related terms.Its explicit expression is suppressed as it is of no interest for the following derivation.While Eq. ( 8) is valid for any grid point, it is advantageous, however, to select only a small subset, denoted by {k m } M m=1 and tagged as the central band, and to leave the rest as the peripheral band.This is exactly the frequency-selective merit of the proposed method in this work, which can further reduce the computational cost by neglecting the rather wide peripheral band.Note that the tagged grid points do not necessarily have to be contiguous.The proposed method applies also for a few broken intervals as a whole central band.
Once the central band is defined, Eq. ( 8) is promoted to include all the relevant {k m } M m=1 , resulting in where and the Vandermonde matrix To eliminate the unknown matrix S in Eq. ( 9), one can use null projection by right multiplying Z's kernel [19]: It can be verified that ZP = O, P H = P, P 2 = P.
Equation ( 9) is thus reduced to As the whole DFT grid is intentionally divided into two bands, the matrices defined on the grid are correspondingly partitioned as where the subscripts c and p denote the central and peripheral subspaces, respectively.Note that rank , which is the number of sinusoids in the central band, and X p = O as X is natively defined in the central subspace.
Consequently, Eq. ( 13) can be rewritten in a simpler form: which is actually irreducible if and only if rank(X c P) = J c , whose necessary condition is rank(P) J c , i.e., M J c + L [19].Since Eq. ( 14) plays a key role in the proposed method, it is especially named the principal matrix equation.

C. Angular frequency estimation
It can be shown that A has the same relationship with C, which is already given in Eq. ( 8), as A c does with C c , except for changing b to b c .Furthermore, this relationship can be rewritten in a computationally more favorable form: where C c is all but the first row of C c and C c is all but the last row of C c .A merit of Eq. ( 15) is that it is invariant under the change in representation with a nonsingular matrix T c ∈ C J c ×J c .This can be concluded from the transforming expression In addition, such a change in representation also preserves C c 's range.Therefore, we will not specifically distinguish C c from C c as they both denote the same invariant subspace.
According to the principal matrix equation ( 14), C c can be exactly recovered from the singular value decomposition (SVD) of ZYP by comprising the left singular vectors corresponding to all J c nonzero singular values [19].
Once C c is obtained, Eq. ( 15) can be employed to compute A c .However, C c is, in general, a nonsquare matrix, entailing that its inverse is ill defined.Consequently, A c can, at best, be formally solved by means of Moore-Penrose inverse C + c , provided that L J c + 1 [19]: Afterwards, the unknown angular frequencies in the central band can be deduced from A c 's eigenvalues.

D. Noise corruption
To project our mathematical model in Eq. ( 1) to the physical world, an inevitable noise term w n must be added to the right side.This additive noise is assumed to be wide-sense stationary, with its first and second moments given as where E{•} denotes the ensemble average and the second moment is understood in a cyclic sense such that r n is N periodic.Note that such an assumption about the noise is quite general but fairly applicable in most, if not all, practical situations.By means of DFT of Eqs. ( 17) and (18), the noise's first and second moments in the frequency domain read where . According to the Wiener-Khinchin theorem, {ρ k } N−1 k=0 is just the noise's PSD spectrum, which is real and nonnegative [20].
After considering the noise term, it can be shown that the principal matrix equation ( 14) is updated to where In light of the widesense stationarity of the noise, it is advantageous to take the quadratic form of Eq. ( 21), followed by ensemble average: where Eq. ( 19) has been referred to.Furthermore, by virtue of Eq. ( 20), it can be shown that the noise term in Eq. ( 22) has a simple, closed form: ) where the real number p m is the mth diagonal element of the Hermitian matrix P.
As for the data term on the left side of Eq. ( 22), it is unfeasible to handle the ensemble average in practice.We 053310-3 will then use the time average • , which is an average over a number of similar data segments acquired at different times, to approximate the ensemble average.As a result, the sinusoidal term on the right side of Eq. ( 22) is estimated as Similar to using Eq. ( 14), Eq. ( 24) can instead be employed to obtain C c , except that SVD ought to be replaced by the eigenvalue decomposition (EVD) owing to the Hermitian structure of C c X c PX H c C H c .Since the noise contributes nonlocally to the whole PSD spectrum of the Schottky data, the right side of Eq. ( 24) is usually of full rank.Hence, all the eigenvalues are nonzero.It is reasonable, however, to take those eigenvectors corresponding to the J c greatest eigenvalues, which are, meanwhile, significantly greater than the remaining ones, and form a matrix to estimate C c .The goodness of this estimate certainly depends on the SN ratio and the time-averaging number, just like the conventional periodogram does for the PSD estimation.

E. Remarks on numerical computation
Although the Vandermonde matrix Z defined in Eq. ( 11) has a well-structured analytic expression, it is notorious for its ill-conditioned behavior in numerical computation, owing to the exponential inflation of the condition number [21].Hence, special attention must be paid to the numerical stability of the proposal method.
For instance, it is disfavored to directly employ Eq. ( 12) to compute the projection matrix P, as the matrix inversion is computationally very costly and, unfortunately, rather inaccurate.Instead, matrix factorization such as SVD or, more efficiently, QR factorization should be employed to obtain P. In short, let Z H be factorized as where Q and Q together form a unitary matrix and R is an upper triangular matrix.It can be shown by substituting Eq. ( 25) into Eq.( 12) that Note that the effective rank(Z) can sometimes be less than its theoretical value because of the limited machine precision.Therefore, a rank revealing QR factorization is recommended to obtain Q [22].
Owing to the same consideration of the numerical stability, it is advised to separate the multiplication of Z and Z H from EVD when estimating C c by using Eq.(24).That is, let U c ∈ C M×J c be a part of the eigenvectors of YPY H − D, which corresponds to the J c most prominent eigenvalues.The product ZU c also estimates C c .
Once the estimate of C c is obtained, Eq. ( 15) will be used to solve for A c .However, the estimation uncertainty propagates to both sides of Eq. (15).Hence, the Moore-Penrose inverse used in Eq. ( 16), which is, in fact, the solution in the leastsquares sense, should be replaced by the total-least-squares solution to account for this change [19].In short, let V c ∈ C 2J c ×2J c be all the right singular vectors of the compound In the end, the algorithm of the proposal method for the frequency-selective harmonic retrieval is summarized as follows.
(1) Start with the Schottky data y n .
(2) Transform y n to the frequency counterpart ỹk by means of FFT.
(3) Select the central band in which J c harmonics are to be retrieved, and construct Y by virtue of Eq. ( 10).
(5) Find Z's kernel Q by virtue of Eq. ( 25), and construct P by virtue of Eq. ( 26). ( 6) With the a priori noise's PSD spectrum ρ k , construct D by virtue of Eq. ( 23).(7) Compute YPY H − D, followed by EVD to find U c , which is the eigenvectors corresponding to the J c most prominent eigenvalues.
(8) Estimate C c by ZU c , and solve for A c by virtue of Eq. ( 15) in the total-least-squares sense.
(9) End with the retrieved angular frequencies, which are the arguments of A c 's eigenvalues.
Note that the proposed algorithm depends on two free parameters, L and M, which should fulfill the relationship M J c + L 2J c + 1. Numerical experiments have suggested that L = J c + 1 and 20 M 30 can usually yield the most accurate results.As a demonstration, a PYTHON implementation of the algorithm can be found in Ref. [23].

III. EXAMPLE
In this section, the superiority of the proposed statespace-based method over the conventional periodogram-based method will be presented with synthetic data for the following formula: where f is fixed to be 0.2 while δ f can be either 0.002 or 0.004, the complex noise w n is white Gaussian, with a variance of 0.1, and the amplitude α is the same for the two harmonics and is chosen to be 0.0356 such that the SN ratio within the interval [0.18, 0.22] is −5 dB.This interval will later be used as the central band for the proposed method.
In the conventional approach to the harmonic retrieval, the PSD spectrum of the data is estimated by its periodogram, which is merely the modular square of the data's windowed DFT.The window function can trade the spectral leakage for the frequency resolution, where the latter can best be achieved without any windows or, equivalently, with an implicit boxcar window [11].The finest frequency resolution is, unfortunately, limited by the granularity of the DFT grid, which is 1/N, with N being the total number of grid points.If the frequencies of two equally strong harmonics differ by less than 1/N, they will be indistinguishable on the periodogram.Such a phenomenon, as demonstrated in Fig. 1, is akin to the Rayleigh criterion of resolving two optical point sources from their diffraction pattern [24].The solid lines in Fig. 1 show the periodogram of synthetic data of length N = 500, after time averaging 20 similar spectra and zooming in on the vicinity of the two harmonics.As the separation δ f decreases from 0.004 to 0.002, the two distinct peaks degenerate into one peak, prohibiting a clear resolution.
As part of employing the proposed method, the parameter L is selected to be 3 because J c = 2, and the predefined central band shown by the shaded regions in Fig. 1 corresponds to M = 21.Moreover, the data length and the time-averaging number are set to be the same as those used for the periodogram.The retrieved harmonics with an accuracy around 10 −4 by means of the proposed method are indicated by the dashed lines in Fig. 1 and are also listed in Table I.Notably, the eigenvalues of YPY H − D are plotted in a decreasing order in the insets of Fig. 1, which can help estimate the parameter J c in case it cannot be determined in advance.It is evident that the first two eigenvalues significantly surpass the remaining ones regardless of whether the Rayleigh limit of the separation δ f is reached, which indeed validates J c = 2 in this example.In this sense, the proposed method can achieve a superresolution, exceeding the Rayleigh limit of the periodogram-based method.TABLE I. Harmonic retrieval result of two harmonics defined in Eq. ( 27) when they are separated by δ f = 0.004 and δ f = 0.002, respectively.

Set frequency
Retrieved frequency Deviation

IV. APPLICATION
Although the proposed method surpasses the conventional method in certain cases, the latter is still complementary as it reveals some a priori information for the former.For instance, the estimated PSD spectrum by the periodogram can coarsely indicate the spectral location of the target ions to help define the central band for the proposed method.Moreover, the noise's PSD spectrum, which is an important input to the proposed method, can be estimated by the periodogram, either of particular background data with no ion signals or of normal measurement data where the background can be extracted in situ by SVD-based data smoothing [25].In this section, the applicability of the proposed method in SMS data analysis will be demonstrated with real data from one beam time.
In an isochronous SMS experiment conducted in the heavy-ion research facility in Lanzhou (HIRFL) [26], the primary beam 58 Ni 19+ was accelerated by the synchrotron CSRm to an energy of 393.165MeV/u, followed by bombarding onto a 15-mm Be target.The projectile fragments along the isospin line T z = 1 were selected and purified by the fragment separator RIBLL2, and injected into the storage ring CSRe.The ring's lattice was set to an isochronous mode on the T z = 1 nuclides with the transition point γ t = 1.313.The Schottky signals of the stored ions were detected by the Schottky resonator installed in CSRe's one straight section [27].During the data acquisition, the local oscillator was tuned at f l = 243.2MHz to match the resonant frequency of the detector, and a rather wide frequency span of 2 MHz was set to accommodate all the potential T z = 1 nuclides.
Figure 2 shows the periodogram estimate of the PSD spectrum of the Schottky data after 5-s time averaging.The Kaiser window with a parameter of 20 is employed to compute the periodogram [28].Clearly visible are many sharp ion peaks, in addition to the rather broad resonant background.A few selected peaks, which belong to some high-yield T z = 1 nuclides, are highlighted by the shaded regions in Fig. 2, with close-ups in the insets.Other peaks are their transverse bands or belong to contaminating ions that do not fulfill the isochronous condition.
The zoomed spans are intentionally kept the same for all the peaks and will also be used as the central bands for the proposed method to retrieve their corresponding harmonics.Such a span translates to the parameter M = 23.In addition, we will naturally set the number of harmonics J c = 1 for each central band.However, as already hinted by the shape and the width of the 58 Ni 28+ peak, J c is instead set as 2. This is due to the significant space-charge effect in the quite intense 58 Ni 28+ beam (around 2 × 10 6 ions) generating a parasitic peak adjunct to the main peak [29].The noise's PSD is estimated by fitting another background spectrum with the Lorentzian function superposed on a possible linear trend.
The harmonics retrieved by means of the proposed method are indicated by the dashed lines in the insets of Fig. 2 and, more explicitly, listed in Table II.In total, six species of bare ions with T z = 1 are identified with an in-house program [30]: 58 Ni 28+ , 56 Co 27+ , 54 Fe 26+ , 52 Mn 25+ , 50 Cr 24+ , and 48 V 23+ .Owing to their adjacent positions on the nuclear chart and hence their similar mass-to-charge ratios m/q, they all have comparable revolution frequencies f r .As a result, their spectral peaks on the PSD spectrum manifest a structured pattern.Specifically, the peaks on the left shoulder of the resonant background all have the harmonic number h = 161, and h = 162 for the peaks on the right shoulder.
Figure 2 shows that 54 Fe 26+ and 52 Mn 25+ are both present at the 161st and the 162nd harmonics.For those two nuclides, their revolution frequencies can be obtained in two ways, either by calculating the spectral distance between their neighboring harmonics or by virtue of Eq. (2).Hence, the self-consistency of the retrieved spectral locations can be tested as follows.First, the spectral distances between the two harmonics are 1.506 359 and 1.505 005 MHz for 54 Fe 26+ and 52 Mn 25+ , respectively.By presuming these are the actual f r , h can correspondingly be calculated to be 161.003and 162.003 same for 54 Fe 26+ and 52 Mn 25+ .Since these numbers are very close to integers, the harmonic retrieval result is indeed self-consistent.Furthermore, the accuracy of the result can be tested by relating to the physical reality.First, each nuclide's f r is obtained by virtue of Eq. ( 2) with the boldfaced values in Table II.Since those nuclides were all under the isochronous condition, their respective f r is independent of the orbital length in the first order approximation [8].Therefore, all the ions can be assumed to circulate on the same orbit and experience the same magnetic field.It can then be shown that f r is determined by m/q with just two parameters, a and b [31]: Figure 3 shows the f r -m/q plot for the six identified nuclides, where the nuclear masses are taken from the latest atomic mass evaluation AME2016 with a proper deduction of the orbital electrons' contribution [32].It is remarkable that the simple model in Eq. ( 28) can well fit the data points by using the harmonic retrieval result as a key input.The relative residual is only about 1 ppm, or even less, for all the nuclides.This excellent goodness of fit has supportively confirmed the validity of the proposed method.Note that the obtained accuracy aligns well with a previous machine study on the isochronous correction at CSRe, as characterized by the relative uncertainty of the revolution time of an isochronous nucleus, reported to be 1.34 × 10 −6 in Ref. [33].
Finally, the applicability of the proposed method can be tested against weak peaks in the periodogram shown in Fig. 2, where a zoomed region around −750 kHz is shown in Fig. 4. For the single peak on the left and the double peak on the right 053310-6 FIG. 3. Top: revolution frequencies of the six nuclides as a function of their mass-to-charge ratios.The solid line is the fitting result of the data points by virtue of the model in Eq. (28).Bottom: fitting residuals relative to the corresponding fitted values. of this region, separate parameter sets are fed to the proposed method to retrieve their respective harmonics.Specifically, M = 21, J c = 1 are set for the single peak, and M = 25, J c = 2 are set for the double peak, while the corresponding central bands are highlighted by the shaded regions in Fig. 4. The harmonic retrieval results, as indicated by the dashed lines in Fig. 4, have proved with their good alignments with the corresponding peaks that the proposed method is also applicable in poor SN ratio cases.

V. CONCLUSIONS
We have presented a frequency-domain approach to the harmonic retrieval problem based on the state space representation.This method exploits the sinusoidal structure of the FIG. 4. Zoomed portion of the periodogram in Fig. 2 showing weak peaks.Note that the ordinate scale changes from logarithmic to linear.The shaded regions define central bands for the proposed method to retrieve the corresponding harmonics, while the dashed lines indicate the retrieval results.harmonic and casts only a general assumption of wide-sense stationarity on the noise.Being complementary to the conventional periodogram-based approach, the proposed method is fully window independent and takes both DFT's modulus and argument into account.Its frequency-selective merit can greatly reduce the computational cost while still being able to yield accurate results.Moreover, its superresolving capability makes the proposed method in particular superior in certain tricky situations, such as when a few harmonics are closely clustered, where the periodogram-based method may fail to distinguish them.Although the method itself has been developed within the context of SMS, it can likely find applications also in a broader field of Fourier-transform mass spectrometry and even beyond, such as in ion-cyclotron-resonance spectrometry [34] or in the direction-of-arrival estimation [35].

053310- 4 FIG. 1 .
FIG. 1. Demonstration of the proposed method's superresolution of two closely spaced harmonics separated by δ f = 0.004 (top) and δ f = 0.002 (bottom).The solid line shows the periodogram, and the dashed lines indicate the retrieved harmonics by means of the proposed method.The shaded region tags the central band used by the proposed method.Shown in the inset are the eigenvalues of YPY H − D, ordered decreasingly.

053310- 5 FIG. 2 .
FIG. 2. Schottky data periodogram of length N = 2000 after 5-s time averaging.The shaded regions indicate the central bands, each of which is used by the proposed method to retrieve the harmonic contained therein, whereas the boxed region around −750 kHz is blown up in Fig. 4. Shown in the insets are the zoomed spectra of identified ions, labeled in the top right corner.The dashed lines in all insets indicate the retrieved harmonics.

TABLE II .
Harmonic retrieval result of the identified ion peaks in Fig.2.Note that in the case with two harmonics corresponding to the same ion, only the value in boldface is employed to compute its revolution frequency.