# Refined betatron tune measurements by mixing beam position data

Panos Zisopoulos, Yannis Papaphilippou, Jacques Laskar

► To cite this version:

Panos Zisopoulos, Yannis Papaphilippou, Jacques Laskar. Refined betatron tune measurements by mixing beam position data. Physical Review Accelerators and Beams, 2019, 22 (7), pp.071002. 10.1103/PhysRevAccelBeams.22.071002. [https://arxiv.org/abs/1907.07100]

HAL Id: hal-02123355

https://hal.science/hal-02123355v1

Submitted on 28 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Refined betatron tune measurements by mixing beam position data

P. Zisopoulos* and Y. Papaphilippou

CERN, CH-1211 Geneva, Switzerland

J. Laskar

IMCCE, Observatoire de Paris, 77 Avenue Denfert-Rochereau, 75014 Paris, France

<figure><img src="image_2.png" /></figure>

(Received 29 March 2019; published 8 July 2019)

The measurement of the betatron tunes in a circular accelerator is of paramount importance due to their impact on beam dynamics. The resolution of these measurements, when using turn-by-turn (TBT) data from beam position monitors, is greatly limited by the available number of turns in the signal. Because of decoherence from finite chromaticity and/or amplitude detuning, the transverse betatron oscillations appear to be damped in the TBT signal. On the other hand, an adequate number of samples is needed, if precise and accurate tune measurements are desired. In this paper, a method is presented that allows for very precise tune measurements within a very small number of turns. The theoretical foundation of this method is presented with results from numerical and tracking simulations and experimental TBT data which are recorded at electron and proton circular accelerators.

DOI: 10.1103/PhysRevAccelBeams.22.071002

# I. INTRODUCTION

Betatron tune measurements [1] are used as a reliable diagnostic of transverse beam dynamics. The measurement of the working point of a circular accelerator is an essential procedure in order to optimize performance and reduce particle losses. In such an accelerator, each beam position monitor (BPM) records the betatron oscillations of the centroid in the transverse plane for many turns. In order to measure the tunes, the beam needs to perform coherent betatron oscillations after transverse excitation from its closed orbit. The turn-by-turn (TBT) transverse oscillations of the centroid are recorded at each BPM, yielding a discrete signal which can be analyzed with algorithms that perform Fourier analysis.

The figure of merit in tune measurements is the resolution in the frequency space, i.e., the smallest difference of adjacent harmonics that can be identified in the discrete Fourier transform spectrum of a BPM signal. In the noise-free regime, the resolution is defined from the error $\epsilon(N)$ in the betatron tune estimation from $N$ turns, which follows a power law of the form [2]

$$ \epsilon(N) = |Q(N) - Q_o| \propto \frac{1}{N^l}, \quad (1) $$

where $Q(N)$ is the estimated tune within a number of turns $N$, $Q_o$ is the actual betatron tune, and $l$ is an exponent which determines the speed of convergence of the tune measurement to the actual tune value.

As is suggested from Eq. (1), a large number of turns is vital for precise betatron tune measurements. However, in experimental and simulated TBT data, the useful number of turns is greatly limited from decoherence due to finite chromaticity and/or amplitude detuning [3,4]. Because of this mechanism, the expectation value of the transverse beam position, as recorded at the BPMs, is strongly damped with respect to the number of turns. Consequently, a choice of a large number of turns is followed from a substantial reduction of the signal-to-noise ratio in the TBT signal. In addition, decoherence induces a significant dependence of the betatron tunes on the number of turns. Apart from nonlinearities, a large number of turns is not always available, e.g., in the case of accelerator commissioning or due to particular beam dynamics measurements, e.g., placing the working point of the accelerator close to a resonance that may lead to strong particle losses.

For the case of a simple fast Fourier transform (FFT) algorithm [5], the exponent $l$ in Eq. (1) is $l = 1$, which is inadequate for a precise determination of the betatron tunes in a circular accelerator. To overcome such limitations, refined frequency analysis (RFA) methods have been developed [2,6,7] which offer a substantially improved resolution in tune estimations, with respect to a plain FFT. The numerical analysis of fundamental

<footnote>*Also at High Energy Physics Department, Uppsala University, Uppsala, Sweden.</footnote>

<footnote>panagiotis.zisopoulos@cern.ch</footnote>

<footnote>Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.</footnote>
frequencies (NAFF) [8] algorithm is an RFA method, which has been widely used with success in beam dynamics measurements [9–12]. This algorithm offers an enhanced decrease of the error in Eq. (1) after applying a Hann window [13] of the order of $p$ on the data. The resolution has been analytically calculated from Laskar to have the asymptotic approximation for $N \to \infty$ [14]:

$$ \epsilon(N) \propto \frac{1}{N^{2p+2}}, \qquad (2) $$

which offers a substantial improvement in the tune estimation, with respect to a simple FFT, for a given number of turns $N$. Nevertheless, a further enhancement of the resolution would be beneficial for the operation of modern circular accelerators and for analyzing TBT data from tracking simulations. The reason for this is that the usage of a very small number of turns ($N \le 50$) results in a weaker influence of decoherence on the tune determination. In addition, the identification of transient effects, such as ripples from power converters or fast-pulsing magnets, could be possible.

In this paper, it is shown that by analyzing the trajectory of the beam from $M$ BPMs and for $N$ turns, while employing the NAFF algorithm, the resolution of the tune measurements is greatly improved, several orders of magnitude compared to single-BPM analysis. The reconstruction can be performed by vectorizing an $N \times M$ array that contains TBT data for $N$ turns and from $M$ BPMs, in a BPM-by-BPM manner. This transformation results in the simultaneous increase of the number of samples ($MN$ from $N$) and of the sampling rate, since after the vectorization, the sampling rate is $M$ samples per turn, instead of the initial one sample per turn.

This technique is referred to as the mixed BPM method and, when it is employed together with the NAFF algorithm, the TBT error in the tune estimation scales as

$$ \epsilon(N) \propto \frac{1}{M^{2p+1}N^{2p+2}}, \qquad (3) $$

where the gain in resolution with respect to the single-BPM analysis Eq (2) is pronounced from the factor $M^{2p+1}$. The complete form and derivation of Eq. (3) is presented in the Appendix.

The mixed BPM method consists of flattening the $M \times N$ array in a BPM-by-BPM manner, i.e., transforming it into a vector of $1 \times MN$ dimension. However, the transformation introduces two systematic errors: (i) an error in the sampling period of the mixed BPM signal due to the fact that the BPMs are not strictly equidistant, i.e., they are not distributed homogeneously around the ring, and (ii) an error due to the BPM-by-BPM modulation of the beta function. Fortunately, the previous errors are periodic in nature, since they are repeated for every revolution of the beam. This characteristic

allows for the aforementioned analysis of the trajectory of the beam in a BPM-by-BPM manner.

An interesting application of evaluating the tunes from multiple BPMs with a different approach through the continuous Fourier transform [2] can be found in Ref. [15]. The mixed BPM method, combined with the NAFF algorithm, has been introduced and tested at proton and electron rings [16,17] for precise tune measurements with a very small number of turns. Here, the theoretical foundations of this method are presented, along with results from tracking simulations and experimental measurements. The paper is organized as follows.

In Sec. II, the methodology of the mixed BPM scheme is presented, together with analytical expressions for the transformed Fourier spectra and the improved betatron tune resolution. In Sec. III, a numerical simulation is deployed in order to visualize qualitatively and quantitatively the theoretical results of the method. In Sec. IV, tracking simulations are performed and the method is applied for tune measurements with the CERN Proton Synchrotron (PS) ideal lattice, under the influence of linear and nonlinear dynamics. In Sec. V, some experimental results are shown from the PS and the CERN Proton Synchrotron Booster (PSB) and from the Karlsruhe Research Accelerator (KARA) electron light source, where the efficiency of the mixed BPM method is highlighted.

## II. METHODOLOGY

### A. Nonuniform periodic sampling

The transformation of the TBT data from $M$ BPMs and for $N$ turns results in a vector of samples, where $N$ groups each with $M$ samples are created with a periodicity of one turn. However, for the samples of the aforementioned vector, the NAFF algorithm, as with any of the usual Fourier analysis methods, assumes a constant sampling period; i.e., the samples are produced from equidistant BPMs. In general, the sampling intervals $\tau_k$, which are defined as the time interval that the beam needs to transport from BPM $k-1$ to BPM $k$, are not equal, and they can be expressed as

$$ \tau_k = t_k - t_{k-1}, \qquad (4) $$

where $t_k$ and $t_{k-1}$ are the time instants that the beam passes from BPMs $k$ and $k-1$, respectively. The boundary condition $t_0$ can be chosen to be $t_0 = 0$. The periodicity conditions, due to the circular geometry of the ring, suggest that

$$ \tau_k = \tau_{k+M}. \qquad (5) $$

According to Eq. (4), the time instances $t_k$ can be expressed for any $k$ as

$$ t_k = \sum_{n=1}^{k} \tau_n, \qquad (6) $$
and, when the centroid of the beam passes from the BPM k to the BPM k + M, it has circulated the ring once, i.e.,

$$t_{k+M} - t_k = T_o, \qquad (7)$$

where $T_o$ is the revolution period of the beam. Combining Eqs. (6) and (7) yields

$$\begin{gathered} \sum_{n=k+1}^{k+M} \tau_n = T_o, \\ \frac{1}{M} \sum_{n=k+1}^{k+M} \tau_n = \frac{T_o}{M}, \\ \langle \tau_k \rangle = \frac{T_o}{M}. \end{gathered} \qquad (8)$$

which means that, independently of the nonuniform positioning of the BPMs in the ring, the *average* sampling period is bounded. Because of the periodicity condition Eq. (5), the expectation value $\langle \tau_k \rangle$ is constant as $k \to \infty$. The previous expression can be used to construct a relationship for the time instances $t_k$. By introducing an error $\delta_k$, which quantifies the deviation of the BPM k from a longitudinally uniform position, the time instances for $k > 0$ are simply

$$\begin{aligned} t_k &= k \langle \tau_k \rangle + \delta_k \\ &= k \frac{T_o}{M} + \delta_k, \end{aligned} \qquad (9)$$

where the errors $\delta_k$ are considered to be *random independent variables*, given in units of seconds. In the case of single-BPM analysis for $M = 1$, $\delta_1 = 0$. The positions of the BPMs must not overlap; thus, the error values are bounded in

$$-\frac{1}{2M} < \delta_k < \frac{1}{2M}. \qquad (10)$$

As a consequence, the expectation value of the $\delta_k$ errors as $k \to \infty$ is

$$\langle \delta_k \rangle = 0. \qquad (11)$$

Because of the periodicity of the ring, the errors are one-turn periodic, i.e., $\delta_k = \delta_{k+M}$, which implies that the variance of $\delta_k$, $\sigma_{\delta_k}^2$, is also bounded, for all k. The previous conditions highlight the fact the nonuniform sampling of the beam around the ring is *stationary*, which allows the reconstruction of the trajectory of the beam simultaneously from all the BPMs. Therefore, the acquired signal, after the transformation of the TBT data with the mixed BPM method, can be considered as a band-limited signal with a nonuniform but *recurrent* sampling scheme [18].

<figure><img src="image_3.png" /><figcaption>FIG. 1. A hypothetical ring with eight BPMs at longitudinal positions which are marked with red circles. When the mixed BPM method is employed, a sampling error $\delta_k$ is introduced, due to the deviation of the BPM positions from hypothetical locations that divide the circumference of the ring in exactly eight equal parts, marked with blue circles. BPM 1 is set as the reference point.</figcaption></figure>

The nonuniformity of the sampling process can be visualized in the fictitious ring of Fig. 1, where eight BPMs (blue circles) are placed at locations where they divide the circumference of the hypothetical ring in eight equal parts. The injection point is assumed to be at the position of BPM 1 (black marker), and the beam follows the clockwise direction. As the centroid of the beam rotates around the ring during one turn, the eight fictitious and longitudinally symmetric BPMs (blue circles) sample the transverse coordinates with a constant sampling period of $T_o/M$. The actual BPMs (red circles) are situated in longitudinal positions that deviate from the symmetric positions by an offset $\delta_k$ (red arc). This error is assumed to be positive for a real BPM which is downstream from the symmetric position and negative for a BPM which is upstream from the symmetric position.

## B. Frequency spectra

The TBT data from multiple BPMs can be represented with an array. Let the array A contain the spatial and temporal histories of the centroid of the beam:

$$A = \begin{bmatrix} z_{11} & \cdots & z_{1M} \\ \vdots & \ddots & \vdots \\ z_{N1} & \cdots & z_{NM} \end{bmatrix}, \qquad (12)$$

where z can be either x or y transverse planes, N is the number of turns, and M is the number of BPMs. Each column represents the signal from one BPM, with a
<figure><img src="image_4.png" /><figcaption>FIG. 2. Synthesis of a one-dimensional signal of 84 000 samples (top right), from an initial two-dimensional signal of N = 2000 and M = 42 (top left). The first turn from five BPMs are shown (bottom left). The signals from each BPM are represented with different colors. The transformation of the five turns results in 210 samples (bottom right) which increases the sampling rate. Each period is now populated by 42 samples of the original signal, and the periods are shown between the red dashed lines.</figcaption></figure>

sampling period of $\tau_s = 1$ turn per sample. The traditional frequency analysis performs the tune measurements for each column of the array in Eq. (12) and then provides an average tune estimation from the M individual observations. Because of the correlations between the data in A, an increase of the sampling rate of the system is possible by vectorizing the array (12) as

$$ \tilde{A} = [z_{11}z_{12}\dots z_{NM-1}z_{NM}]. \quad (13) $$

An example of the transformation is shown in Fig. 2, where 2000 turns from 42 BPMs (top left) are transformed to a single vector of 84 000 samples (top right). Indeed, the information of the first five turns (bottom left) is transformed to a signal of 210 samples.

The vector $\tilde{A}$ contains NM samples of the TBT data, and the new sampling period becomes $\tilde{\tau}_s = \frac{1}{M}$, i.e., 1 turn per M samples.

Although the sampling process of TBT BPM data is usually described in terms of time, the same procedure can be equally described in space. The mixed BPM signal of Eq. (13) can be acquired by sampling the pseudoharmonic oscillation

$$ z(s) = \mathrm{Re} \left[ \sqrt{\epsilon \beta(m)} e^{i[2\pi\Psi(s) + \phi_0]} \right] \quad (14) $$

along the ring, where s is the longitudinal variable, $\phi_0$ and $\epsilon$ are constants of the motion, $\beta(s)$ is the beta function, and

$$ \Psi(s) = \int_{0}^{s} \frac{ds}{\beta(s)} \quad (15) $$

is the cumulative phase function with the constraint

$$ \Psi(C) = Q_z, \quad (16) $$

where C is the circumference of the ring and $Q_z$ is the transverse betatron tune. The methodology consists of sampling Eq. (14) with M BPMs, which are distributed at different longitudinal positions {$s_1, s_2, \dots, s_M$}. Since the new signal consists of discrete samples, the continuous variable s is dropped in favor of a discrete variable m. The generating function Eq. (14) becomes

$$ z(m) = \mathrm{Re} \left[ \sqrt{\epsilon \beta(m)} e^{i[2\pi\Psi(m) + \phi_0]} \right] \quad (17) $$

and the constraint $\Psi(M+1) = Q_z$.

The integral $\Psi(M+1)$ can be divided into M equal parts, and each part will advance the phase by a constant value of

$$ \Psi_0 = \frac{Q_z}{M}. \quad (18) $$

By introducing the aforementioned sampling errors $\delta_m$ in units of $2\pi$, for each m sample, the cumulative phase $\Psi(m)$ is

$$ \begin{aligned} \Psi(m) &= m\Psi_0 + \delta_m \\ &= m\frac{Q_z}{M} + \delta_m. \end{aligned} \quad (19) $$

With Eq. (19), Eq. (17) is written as

$$ z(m) = \mathrm{Re} \left[ \sqrt{\epsilon \beta(m)} e^{i\{2\pi[m(Qz/M) + \delta_m] + \phi_0\}} \right]. \quad (20) $$

The periodicity of the beta function and of the error $\delta_m$ suggests that $\beta(m+M) = \beta(m)$ and $\delta_{m+M} = \delta_m$; i.e., the signal Eq. (20) is modulated in amplitude and phase with a period of M.

In order to highlight the modulation of the signal, Eq. (20) is expressed as

$$ \begin{aligned} z(m) &= \mathrm{Re}[(\sqrt{\epsilon\beta(m)}e^{i(2\pi\delta_m + \phi_0)})e^{i2\pi m(Qz/M)}] \\ &= \mathrm{Re}[H(m)e^{i2\pi m(Qz/M)}]. \end{aligned} \quad (21) $$

The function $H(m)$ is M periodic, and it can be expanded in discrete Fourier series [19]

$$ H(m) = \sum_{k=-M/2}^{M/2} C_k e^{ik(2\pi/M)m}, \quad (22) $$

with the weight functions for each harmonic defined as

$$ C_k = \frac{1}{M} \sum_{-M/2}^{M/2} H(m) e^{-ik(2\pi/M)m}. \quad (23) $$
Substitution of Eq. (22) in Eq. (21) results in

$$z(m) = \mathrm{Re} \left[ \sum_{k=-M/2}^{M/2} C_k e^{i2\pi[(k+Qz)/M]m} \right]. \quad (24)$$

Inspecting Eq. (24), the following observations can be made.

(i) Since the signal is modulated periodically in amplitude from the optics and in phase from the sampling error, sidebands appear in the spectrum. The carrier frequency is the betatron tune which is found at the $k = 0$ harmonic. The distance from the rest of the harmonics in the frequency spectra is $1/M$, which is the transformed sampling frequency. After transformation to the original frequency space by multiplying the frequencies with $M$, this distance becomes 1. Therefore, the fractional betatron tune can always be recovered, since the uncertainty between the harmonics is exactly one integer tune unit.

(ii) The bandwidth of the signal in Eq. (22) is bounded by the Nyquist frequency. In the case of $M$ BPMs and according to Shannon's sampling theorem [20], the following relationship holds for the betatron tune:

$$Q_z \le \frac{M}{2}, \quad (25)$$

from where it can be deduced that, if the number of BPMs $M$ is at least twice the tune, then the integer part of the tune can be also recovered. Moreover, since the bandwidth of the signal has now become $M$ times larger, the usual discrepancy in the modulo of the fractional tune from single-BPM analysis is no longer present. This implies that tunes with fractional parts above 0.5 can be recovered.

### C. Frequency resolution

The frequency resolution, i.e., the TBT error in the estimation of the betatron tunes, is correlated to the total observation time of a signal. For the mixed BPM method, the derived relationship of the time instances $t_k$, in the case of nonequidistant $M$ BPMs with a sampling error $\delta_k$, is (see Sec. II A)

$$t_k = k \frac{T_o}{M} + \delta_k, \quad (26)$$

with $T_o$ the revolution period. If the total number of samples is $m = MN$, the error $\delta_m = \delta_M$, since the last sample is sampled from the last BPM $M$. Moreover, Eq. (26) can be rewritten for $k = m$ as

$$t_m = m \frac{T_o}{M} g(m), \quad (27)$$

where $\tilde{\delta}_M$ is the sampling error for the BPM $M$, normalized to the revolution period, and the function $g(m)$ is defined as

$$g(m) = \left( 1 + \frac{M \tilde{\delta}_M}{m} \right). \quad (28)$$

The expression in Eq. (27) corresponds to the total observation time of the mixed BPM signal for $m = MN$ samples. This total time is now used for the cases of FFT and NAFF algorithms in order to estimate the error in the frequency analysis of each algorithm.

### 1. FFT

In the case of a simple FFT, the TBT error of the betatron frequency estimation for $m$ number of samples is

$$\Delta\nu(m) = |\nu(m) - \nu_o| = \frac{1}{t_m}, \quad (29)$$

with $\nu(m)$ the time-dependent frequency estimation and $\nu_o$ the true frequency of the TBT data. The error $\epsilon(m)$ in the betatron tune $Q_o$ estimation within $m$ samples is defined as

$$\epsilon(m) = \Delta\nu(m)T_o = \frac{T_o}{t_m}, \quad (30)$$

and substitution of Eq. (27) in the previous expression, with $m = MN$, yields

$$\epsilon(N) = \frac{1}{N + \tilde{\delta}_M}. \quad (31)$$

Since the error $|\tilde{\delta}_M| < \frac{1}{2M} \ll N$, Eq. (31) can be expanded around $\tilde{\delta}_M \approx 0$ to give

$$\epsilon(N) = \frac{1}{N} - \frac{\tilde{\delta}_M}{N^2} + \mathcal{O}(\tilde{\delta}_M^2). \quad (32)$$

Note that this result is consistent with Eq. (1), for $l = 1$ and for single-BPM analysis, i.e., $\tilde{\delta}_M = 0$. Clearly, the mixed BPM method when used with the FFT does not result in any gain in convergence with respect to the single-BPM analysis, $M = 1$. Moreover, the sampling error has been introduced in the estimation of the tunes, although it converges faster to zero than $1/N$. It should be noted that the mixed BPM method still allows the estimation of the integer part of the tune with a simple FFT, provided that the condition in Eq. (25) is fulfilled.

### 2. NAFF

In the case of the NAFF algorithm, the signal under study is treated as *quasiperiodic* [21]. The resolution of the tune measurements scales as $t_m^{-(2p+2)}$ for the total observation time $t_m$ in Eq. (26). For the limiting case of $\tilde{\delta}_M \to 0$ and for $m = MN$ number of samples and by using the theory behind NAFF [14], the analytical relationship of the error in the estimation of the tunes $\epsilon(N)$ is found to be (see the Appendix)
$$ \epsilon(N) = \frac{\overline{C_L}}{M^{2p+1}} \left( \frac{1}{N^{2p+2}} - (2p+2) \frac{\tilde{\delta}_M}{N^{2p+3}} \right), \quad (33) $$

where the factor $\overline{C_L}$ depends on the number of samples $m$, the error $\tilde{\delta}_M$, the order of the Hann window $p$, and the betatron frequencies and amplitudes of the signal under study. In fact, in the case of rational frequencies and/or the presence of multiple harmonics in the signal, the factor $\overline{C_L}$ can diverge rapidly. It is important to mention that the gain in resolution by a factor of $M^{2p+1}$ is followed by a blowup of the factor $\overline{C_L}$, as it can be confirmed from the full derivation of Eq. (33) in the Appendix. Moreover, due to the dependence of $\overline{C_L}$ on the order of the Hann window, the blowup is expected to increase for an increasing $p$.

The transformation of the signal with the mixed BPM method creates a spectrum with multiple harmonics around the main frequency line, as has been shown in Sec. II B. This behavior can also interfere with the convergence of the betatron tune determination for a small number of turns $N$, due to the dependence of the $\overline{C_L}$ on the number of harmonics. As a result, collections of BPMs with the least variation of the optics and the sampling period are expected to produce more precise betatron tune measurements for a very low number of turns $N$.

The contribution of the sampling error $\tilde{\delta}_M$ is negligible, since it is very small by construction and it rapidly approaches zero. Thus, for the mixed BPM method, the two previous effects that can make the factor $\overline{C_L}$ diverge can potentially reduce the improvement in the estimation error of the betatron tunes.

### III. NUMERICAL SIMULATIONS

Numerical simulations are performed with pYNAFF [22], in order to qualitatively investigate the theoretical derivations of the mixed BPM method. Since, in a real machine, the TBT data from a BPM resemble pseudoharmonic oscillations, the numerically simulated signal is chosen to be a superposition of four harmonic terms:

$$ z(m) = \sum_{k=1}^{4} e^{i2\pi Q_k [(m/M) + \delta_m]}, \quad (34) $$

where $m$ is the index of each sample with the constraint $1 \le m \le MN$, $M$ is the number of BPMs, $N$ is the number of turns, $Q_k$ is the tune of the harmonic $k$ in $2\pi$ units, and $\delta_m$ is the sampling error of sample $m$. The complex signal in Eq. (34) contains information on the positions (real part) and the momenta (imaginary part) of the oscillations. In these numerical simulations, the real part of Eq. (34) is used, so as to simulate the types of signals that are acquired normally in an actual accelerator. The projection of the simulated TBT oscillations on the complex plane is presented in Fig. 3. The values of the frequencies are

<figure><img src="image_5.png" /><figcaption>FIG. 3. The projection of the signal used in the numerical simulations, on the complex plane. The oscillations consist of the superposition of four frequencies, with a frequency-to-frequency shift of $\Delta Q = 0.05$.</figcaption></figure>

chosen to be $Q_k = \frac{1}{2\pi} + (k-1)\Delta Q$, and the frequency separation between the harmonics is $\Delta Q = 0.05$.

The goal of the simulation is to use the mixed BPM method with NAFF, in order to measure the frequencies of the signal in Eq. (34). The uncertainty $\epsilon_k$ in the estimation of the $Q_k$ harmonic is defined as

$$ \epsilon_k(N, M) = |Q_{k_o} - Q_k(N, M)|, \quad (35) $$

where $Q_{k_o}$ are the actual tunes of the numerical signal Eq. (34) and $Q_k(N, M)$ is the frequency estimated over $N$ turns and $M$ BPMs. In this analysis, only the results of the $k=1$ harmonic will be presented, i.e., of the main frequency. The $N$ variable is constrained to low values ($N \le 50$) in order to highlight the contribution of the parameter $M$ in Eq. (33). Mixing the BPM data together results in the increase of the sampling rate from one sample per turn to $M$ samples per turn. Therefore, the frequencies $Q_k$ are transformed to $Q_k/M$ which will be referred to as the reduced frequencies.

For some $M$, the $Q_k/M$ ratio could yield a figure equal or almost equal to a rational number, which would result in the inability to reconstruct the original frequency with the convergence of Eq. (33). For example, for $M=20$, $Q_1/20 \approx 1/2$, which is an even resonance. Indeed, this behavior is shown in Fig. 4, where the mixed BPM method is applied and the error in the tune estimation $\epsilon_1$ is measured for an increasing number of BPMs $M$ and for three cases of $N$. Although the general trend shows a decreasing error with respect to $M$ for all cases of $N$, the convergence curves are contaminated due to the appearance of resonance peaks at specific numbers of BPMs $M$. The even resonances are indicated with full lines, while the odd resonances with dashed lines. The comparison of the trends of the curves suggest a decrease of the error for an increasing $N$ as expected; however, the gain in convergence would be more evident in a “nonresonant” case.

In order to bypass these constraints and to demonstrate the convergence of the tune estimation with respect to the
<figure><img src="image_6.png" /><figcaption>FIG. 4. Appearance of resonances in the mixed BPM analysis of the numerical simulations, due to a rational $Q/M$. The measurement error $\epsilon$ is shown in the logarithmic scale, for three cases of $N$, against the number of BPMs $M$. The even resonances are indicated with thick lines, while the odd resonances with dashed lines.</figcaption></figure>

number of BPMs $M$, a varying frequency is introduced, where, for each case of $M$, the $Q_k$ values in Eq. (34) are multiplied with $M$, in order to keep the $Q/M$ ratio constant. In this way, the generated TBT data do not lock on to the aforementioned resonances and no systematic errors are introduced in the analysis.

### A. Tune convergence for increasing $M$

During the simulations, the sampling rate is kept constant and $\tilde{\delta}_m = 0$; i.e., the $M$ BPMs are homogeneously distributed around the fictitious ring, and the optics functions (betatron amplitude and phase advance) are equal at each BPM position. The order of the Hann window is set to be $p = 1$. The results of the mixed BPM measurements are shown in Fig. 5, where the error $\epsilon$ defined in Eq. (35) is plotted in the logarithmic scale with respect to the number of BPMs. The same values of $N$ as in Fig. 4 are considered. Obviously, the dependence of the error $\epsilon$ on the number of

<figure><img src="image_7.png" /><figcaption>FIG. 5. Convergence of the mixed BPM method in the numerical simulations, for three cases of turns, with respect to the number of turns. After introducing the varying frequency scheme, the resonances disappear.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIG. 6. Color map of the tune estimation error with the mixed BPM method for a numerical simulation with PyNAFF. Accurate values of the tune, at the order of $10^{-5}$, can be achieved with either a small number of turns ($N \le 10$) or a larger number of BPMs ($M \ge 50$). The same conclusions can be made for $N \ge 40$ and $M \le 30$.</figcaption></figure>

BPMs $M$ follows a power law, and a comparison with Fig. 4 confirms the absence of the resonance lines and the smooth convergence of the error $\epsilon$. For $N = 10$ and $M = 50$, the error is at the order of $10^{-5}$, while for $N = 30$ and $N = 50$, the error is around $10^{-7}$ and $10^{-8}$, respectively. These low errors are expected due to the use of a smooth quasiperiodic signal, Eq. (34), the absence of additive noise in the signals, and the constant sampling rate and optics. A fit of the convergence $\epsilon$ with a model of the form

$$y(M) = c_1 + c_2 \log_{10}(M), \quad (36)$$

where $y(M) = \log_{10}(\epsilon_1)$, $c_1$ is a constant term, and $c_2$ is the exponent of $M$ in the error $\epsilon$, yields $c_2 = -3$, confirming the theoretical dependence of the convergence $\epsilon$ on $M$ for $p = 1$, as is shown in Eq. (33) for $\tilde{\delta}_m = 0$.

The role of the number of BPMs $M$ and turns $N$ in the error of the frequency estimation with NAFF is explored, by scanning over a range of $(M,N)$ values. The resulting surface can be inspected in Fig. 6, where the color bar represents the error $\epsilon$ in the logarithmic scale and up to 50 turns are taken into account. The distribution of the convergence on the $(M,N)$ surface appears to take hyperbolic shapes as a result of the power law that it obeys. At first glance, it is obvious that by increasing the number of BPMs, for a constant number of turns, the achieved error is gradually decreased. For example, below 10 turns, the error reaches $10^{-4}$ for $M = 50$, and it decreases further at around $10^{-7}$ for $M = 100$. The same order of magnitude for the error can be achieved for $M = 50$ and $N = 40$ turns.

### B. Results for a varying window order

The analytical relationship of the tune estimation error in Eq. (33) suggests a strong dependence on the Hann window order $p$. In frequency analysis, a window is always used in order to reduce the impact of the finite sampling rate and
<figure><img src="image_9.png" /><figcaption>FIG. 7. The coefficients estimated from the fit of the tune error $\epsilon_1$ to the surface of Eq. (37), along with their uncertainty. The $c_1$ coefficient is marked with blue circles, $c_2$ with blue squares, and $c_3$ with red diamonds. The error bars represent the $1\sigma$ standard error of the fit. The theoretical predictions for the coefficients are shown with a blue dashed and thick line for $c_1$ and $c_2$, respectively. The theoretical estimation of $C_L$ is shown with a red line.</figcaption></figure>

finite duration of the signal, i.e., to compensate the leakage effect. For the case of a Hann window, the efficiency of the compensation depends strongly on the order of the window. Although a higher-order window leads to a faster damping of the error $\epsilon$, studies have been performed [23] on the effect of different orders, and the conclusions suggest that in the presence of noise in the BPM signal, e.g., from electronics, larger orders of the window lead to a greater loss of precision. Analytical estimations of the impact of white noise in tune estimations can be found in Ref. [15].

For the numerical simulations presented in this section, frequency analysis is performed again for different values of $(M,N)$ and for different values of the Hann window order $p$. The estimated convergence $\epsilon$ is fitted on the surface:

$$y(M, N) = c_1 \log_{10}(M) + c_2 \log_{10}(N) + c_3, \quad (37)$$

where $y(M, N) = \log_{10}(\epsilon_1)$, $c_1 = 2p + 1$, $c_2 = 2p + 2$, and $c_3 = \log_{10}(C_L)$. The estimated values of the fit coefficients are presented in Fig. 7, where each coefficient is shown for different orders of $p$, along with the theoretical expectations. Indeed, it is evident that the numerical results agree with the theoretical findings. The reconstruction of the TBT signal in Eq. (34) with $M$ BPMs results in a reduction of the error by a factor of $M^{2p+1}$, as is predicted by Eq. (33). At the same time, the $C_L$ coefficient is increasing with respect to the window order $p$, as is also expected from theory.

## IV. TRACKING SIMULATIONS

The efficiency of the mixed BPM scheme is also tested with simulations using MADX-PTC [24] and an optics model of the PS [25]. The parameters of the simulation can be found in Table I. No collective effects are considered in

<table><caption>TABLE I. Parameters of the MADX-PTC tracking simulations with the PS model.</caption><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>Energy</td><td>2.3 [GeV]</td></tr><tr><td>M</td><td>42 BPMs</td></tr><tr><td>Q<sub>x</sub>, Q<sub>y</sub></td><td>6.24, 6.27 [2π]</td></tr><tr><td>Q'<sub>x</sub>, Q'<sub>y</sub></td><td>-5.78, -7.66 [2π]</td></tr><tr><td>Emittance ε<sub>x</sub>, ε<sub>y</sub></td><td>1.0, 0.8 [μm rad]</td></tr><tr><td>Beam size σ<sub>x</sub>, σ<sub>y</sub></td><td>4.6, 3.0 [mm]</td></tr><tr><td>Dimensionality</td><td>4D</td></tr><tr><td>Distribution</td><td>Gaussian</td></tr></tbody></table>

the simulations. The TBT data are recorded at 42 BPMs arranged around the PS lattice. The goal is to use the mixed BPM method for tune measurements with a very small number of turns $N \le 50$, make comparisons with the traditional single-BPM measurement, and explore the influence of the sampling rate error and of the optics variation on the resolution of the tune measurement.

The beam is initially excited, horizontally and vertically, with deflections that correspond to initial amplitudes of $2\sigma$ to $12\sigma$ with a step of $2\sigma$. The centroid oscillations are measured from the average orbit of the particles. The evolution of the centroid of the beam is shown in Fig. 8 with respect to the number of turns $N$, for all the BPMs, and for the different initial conditions. The damping of the oscillations is faster for larger kicks, due to amplitude detuning coming from the nonlinear magnetic elements. The beam exhibits maximum horizontal and vertical amplitude-dependent tune shifts of $(\Delta Q_x, \Delta Q_y) = (1.5 \times 10^{-3}, 2.5 \times 10^{-3})$ for excitations of $2\sigma-12\sigma$.

### A. Sampling of different BPM configurations

The PS ring has a mean radius of $R = 100$ m and consists of ten superperiods, each made of ten combined function magnets. The optics used in these simulations corresponds to the so-called *bare machine* optics, where the low-energy quadrupoles which are used to correct the betatron tunes are not activated and there are no skew elements in the lattice. Comparisons can be made for the efficiency of the mixed BPM method by using different BPM configurations, i.e., groups of BPMs with different levels of variation of the optics and of the longitudinal distance $\Delta s$ between the BPMs. As a consequence, the configurations also have a different sample size, i.e., a different number of BPMs $M$. The beta functions and phase advances for both planes, at the location of the $M = 42$ BPMs (configuration A), are presented in Fig. 9, where the periodicity of the optics, for every four BPMs, is visible. This periodicity is, however, broken in the range of BPMs from 23 to 32. Since the periodicity of the optics is a factor that can potentially interfere with the efficiency of the mixed BPM method, two BPMs which are located at
<figure><img src="image_10.png" /><figcaption>FIG. 8. The turn-by-turn response of the centroid of the beam, for different kick strengths. The top row depicts the horizontal oscillation and the bottom row the vertical. The left column shows the TBT data for an initial amplitude of 4σ, the middle column for 8σ, and the right column for 12σ.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIG. 9. Beta function (top) and phase advance (bottom) with respect to the index of the 42 BPMs of the PS. The horizontal optics are shown in blue and the vertical in orange.</figcaption></figure>

positions that perturb the periodicity are removed, yielding a new configuration (configuration B) with $M = 40$.

Moreover, two interesting collections can be further sampled from configuration B: a total of $M = 30$ BPMs that have almost equal optics, i.e., almost the same values for the beta function and phase advance (configuration C), and a total of $M = 10$ BPMs that have almost equal optics and that are located at equal distances $\Delta s$ from each other (configuration D).

A convenient metric of estimating the variation of the optics and of the longitudinal distance $\Delta s$ for every BPM collection is the standard deviation, which is presented in Table II. From all the configurations, the one which exhibits the least BPM-by-BPM variation in the optics and the longitudinal position is expected to produce the best precision and accuracy in the mixed BPM method.

The Fourier spectra of the horizontal betatron function in one turn are shown in Fig. 10, for all the aforementioned BPM collections. The average value of the beta function has been subtracted from the samples, in order to suppress low-frequency modes. For the 42 BPMs (configuration A), two main peaks arise in the spectrum signifying the almost periodic lattice of the PS. Reducing the number of BPMs to 40 (configuration B) results in a cleaner spectrum with a very well-defined oscillation frequency which corresponds to the tenth harmonic.

The spectra from configuration C are shown in the bottom left figure, where only high-frequency components are present. These components appear due to the fact that the beta functions are almost, but not entirely, equal at the positions of these 30 BPMs. Finally, for the case of configuration D, the spectrum is flat, signifying equal values for the optics at each BPM of this collection. The vertical beta function exhibits similar periodicities.

Because of the one-turn modulations of the optics and of the error in sampling the sampling rate between the BPMs,

<table><caption>TABLE II. The standard deviation of the optics ($\sigma_{\beta_{x,y}}$ for the horizontal and vertical beta functions and $\sigma_{\mu_{x,y}}$ for the horizontal and vertical phases advance) and of the longitudinal distances between the BPMs $\sigma_{\Delta s}$, for the BPM configurations A–D.</caption><thead><tr><th>Configuration</th><th><i>M</i></th><th>σ<sub>β<sub>x</sub></sub><br>[m]</th><th>σ<sub>β<sub>y</sub></sub><br>[m]</th><th>σ<sub>μ<sub>x</sub></sub><br>[2π]</th><th>σ<sub>μ<sub>y</sub></sub><br>[2π]</th><th>σ<sub>Δs</sub><br>[m]</th></tr></thead><tbody><tr><td>A</td><td>42</td><td>4.64</td><td>4.64</td><td>0.044</td><td>0.042</td><td>4.30</td></tr><tr><td>B</td><td>40</td><td>4.53</td><td>4.49</td><td>0.030</td><td>0.030</td><td>4.27</td></tr><tr><td>C</td><td>30</td><td>0.24</td><td>0.40</td><td>0.11</td><td>0.11</td><td>7.19</td></tr><tr><td>D</td><td>10</td><td>0.052</td><td>0.013</td><td>0.000 61</td><td>0.0012</td><td>0.45</td></tr></tbody></table>
<figure><img src="image_12.png" /><figcaption>FIG. 10. Fourier spectra of the horizontal beta functions signifying the periodicities of the optics for each BPM configuration. The name of each configuration is shown in the legend.</figcaption></figure>

sidebands are expected to appear in the Fourier spectra of the mixed BPM signal. For the current chosen BPM configurations, the Fourier spectra for 24 turns are shown in Fig. 11. In the top plot, the appearance of sidebands is shown around the main peak of betatron tune, with an exponentially decreasing amplitude.

By using configuration B, only one sideband appears, and the amplitude of the main peak increases due to the absence of any other spectral components, leading to a much cleaner signal. Furthermore, by using configuration C, the spectral quantity is almost identical. In the case of configuration D, the signal exhibits only one frequency. This component is the betatron frequency; however, due to aliasing which arises from the small number of BPMs, the integer part of the tune is found to be 3 instead of 6. Nevertheless, an important result is that the fractional part of the tune can always be determined, regardless of the number of the BPMs.

<figure><img src="image_13.png" /><figcaption>FIG. 11. The Fourier spectra of the mixed BPM signal for the different configurations of BPMs that are used in the analysis. The name of each configuration is shown in the legend.</figcaption></figure>

## B. Tune measurements

### 1. Precision

A straightforward method for estimating the precision is to measure the difference in betatron tune between consecutive TBT measurements. The results can be visualized in Fig. 12, where the TBT convergence of the measured horizontal (top row) and vertical (bottom row) tune values is shown for the first 50 turns. The results from two different sets of initial conditions are shown: the $4\sigma$ excitation (left column) and the $12\sigma$ excitation (right column). The different colors correspond to three BPM configurations. The case for a single BPM is shown in blue, the case for configuration A, with all the available 42 BPMs, is shown in green, and the case for configuration D, with ten BPMs which are equidistant and have similar optics values, is shown in orange.

For the horizontal plane, the NAFF algorithm estimates the tunes correctly from the first six turns for all BPM configurations and for both initial excitations. For the $4\sigma$ excitation, the single BPM precision is around $10^{-2}$ after six turns, while the mixed BPM method exhibits a precision of $10^{-3}$ for the 42 BPMs of configuration A and even lower for the ten BPMs of configuration D. At 50 turns, the precision is below $10^{-7}$, whereas for the single-BPM case it is at around $10^{-4}$. Regarding the $12\sigma$ excitation, the mixed BPM method is found to be more precise than the single-BPM analysis. However, for all configurations, the convergence of the measured betatron tunes is heavily modulated by the strong decoherence.

The same conclusions can be also drawn for the vertical tunes. Indeed, for the first 50 turns, mixing the BPM data together results in improved precision. For the $4\sigma$ case, results from configuration D exhibit the highest precision,

<figure><img src="image_14.png" /><figcaption>FIG. 12. The precision of the tune measurements for the horizontal (top row) and vertical (bottom row) planes. The results for the single-BPM analysis are shown in blue circles, the mixed BPM scheme with configuration A in green squares, and with configuration D in orange. The initial excitation of the beam corresponds to $4\sigma$ (left column) and $12\sigma$ (right column).</figcaption></figure>
while for the 12σ case, both mixed BPM collections present the same level of precision, almost 2 orders of magnitude better than the single-BPM case.

Similar improvements in precision of the transverse tunes estimation are also found for the rest of the cases of initial excitations.

## 2. Accuracy

A common problem for estimating the accuracy of betatron tune measurements is that the real value of the tune is not known a priori. In principle, the reference value of the tune could be extracted from the transfer matrices of the particle tracking program; however, this would introduce many systematic errors due to decoherence.

An intermediate solution to this problem could be to estimate the betatron tunes of a single particle. These values could then be used as a reference for accuracy measurements. The choice of the particle can be done in the following way: The particle is initially assigned to have a vanishing transverse amplitude; i.e., it is situated at the point $(x, y) = (0, 0)$, where x and y are the horizontal and vertical phase space coordinates, respectively. At $N = 0$, a finite transverse impulse is given to the particle so as to generate betatron oscillations. Then the particle is tracked for many turns (around 1000), and the betatron tunes are estimated for each transverse plane. These values can then be used as the reference tunes. The reason for placing the particle at $(x, y) = (0, 0)$ is that, in this way, nonlinear amplitude detuning is minimized. Furthermore, using a large number of turns minimizes the uncertainty in the estimation of the tunes due to Fourier analysis, which is not affected by decoherence since a single particle is used.

By referring to the single-particle tune estimations as $Q_0$, the TBT accuracy $\epsilon(N)$ can be evaluated as

$$ \epsilon(N) = |Q(N) - Q_0|, \qquad (38) $$

where $Q(N)$ is the TBT tune estimations of the centroid of the beam. The results for the 4σ case are plotted in Fig. 13. In the top plot, the error in the horizontal tune estimation is shown with respect to the number of turns, where the improvement of accuracy for the mixed BPM cases is obvious. For the case of configuration A (42 BPMs), the accuracy is at the order of $10^{-3}$ already after 12 turns, one order of magnitude better than the single-BPM case, while for configuration D (ten BPMs) the accuracy is improved by 2 orders of magnitude. The effect of sampling rate error and modulation of the optics for configuration A is reduced at the very first turns, as predicted from Eq. (33). In this case, there is no obvious gain between configuration C (30 BPMs) and configuration D (ten BPMs). In fact, there is only a marginal increase in accuracy after the very first turns for the latter configuration.

In the bottom picture, the same estimations of the vertical tunes are shown, with similar improvement in the accuracy when the mixed BPM scheme is used. The fact that the

<figure><img src="image_15.png" /><figcaption>FIG. 13. Comparison of the accuracy for the tune measurements with single-BPM (blue curve) and mixed BPM configurations for configuration A (green curve), configuration B (red curve), configuration C (purple curve), and configuration D (orange curve). The excitation corresponds to $4\sigma$ for the horizontal (top) and vertical (bottom) tunes.</figcaption></figure>

single-BPM case for the vertical tune measurements exhibits faster convergence than the horizontal ones can be attributed to the optics. Indeed, by inspecting the working point, the horizontal tune is very close to the fourth-order resonance, which can be excited by strong sextupoles or octupoles.

Similar observations can be made for the case of a larger initial excitation. The accuracy plots are shown for the case of $12\sigma$ in Fig. 14, for the same configurations as in the previous case of initial excitation. In this case, the impact of decoherence is larger; however, the horizontal tunes (top plot) can be resolved at $N = 6$ turns at the order of $10^{-4}$, 2 orders of magnitude better than the single-BPM case, when configuration B is used. The rest of the BPM configurations perform better than the single-BPM case. For the vertical tunes (bottom plot), the mixed BPM method exhibits more

<figure><img src="image_16.png" /><figcaption>FIG. 14. Comparison of the accuracy for the tune measurements with single-BPM (blue curve) and mixed BPM configurations for configuration A (green curve), configuration B (red curve), configuration C (purple curve), and configuration D (orange curve). The excitation corresponds to $12\sigma$ for the horizontal (top) and vertical (bottom) tunes.</figcaption></figure>
accurate results as well. The estimation of the vertical tune when using a single BPM cannot converge before $N = 50$ turns, whereas all the BPM configurations for the mixed BPM method converge at around $N = 25$ turns. The loss of accuracy in the vertical tune determination from a single BPM after $N = 6$ turns is indeed attributed to the inability in this case to converge to the expected value. The increased effect of decoherence is expected to contribute to this effect as well. As a matter of fact, due to decoherence, a loss in accuracy is observed for configuration C at $N = 6$ turns.

The interpretation of the previous results is that the number of BPMs and the periodicity of the optics are important factors in the performance of the mixed BPM method, especially for the very first turns ($N < 20$). The use of configuration B (40 BPMs), which exhibits a well-defined periodicity in the optics with respect to configuration A (42 BPMs), allows for the most accurate tune estimations at six turns and for all initial excitations. The level of the optics variation is also found to be important; for example, configuration D (ten BPMs) performs better than configuration A (42 nonequidistant BPMs) in all cases of excitation. However, configuration B outperforms configuration D, especially for the estimation of the vertical tune, due to the larger number of BPMs (40 versus ten). For all BPM collections, the contribution of the sampling rate error is found to be less important than the optics, as is also expected from Eq. (33). From the same equation comes a fast TBT reduction of all the error contributions (optics and sampling rate variation). Indeed, the present results exhibit a similar behavior for $N > 20$, for all mixed BPM collections, reaching an accuracy level 2–3 orders of magnitude better than the single-BPM results.

### C. Equidistant BPMs

The periodic variation of the sampling period can result in the loss of convergence for a very small number of turns. In order to analyze this effect, a new collection of ten BPMs which are not equidistant with each other is sampled from configuration A (42 BPMs). This collection is intentionally chosen to have a large sampling rate error in order to perform comparisons with configuration D (ten BPMs which are equidistant). The BPMs that belong to configuration D are referred to as *regular*, and the BPMs with large sampling rate error as *irregular*. The normalized error $\tilde{\delta}_k$ in the sampling instance of the BPM $k$ can be estimated by using information of the longitudinal position of the BPMs, $s_k$. Indeed, from Eq. (9), the normalized error can be rewritten as

$$ \tilde{\delta}_k = \frac{t_k}{T_o} - \frac{k}{M}, \qquad (39) $$

where $t_k$ is the time that the beam needs to reach BPM $k$, $T_o$ is the revolution period, and $M$ is the number of BPMs in one turn. The following expressions can be used to express $\tilde{\delta}_k$ in terms of the azimuthal position $s_k$ of each BPM:

<figure><img src="image_17.png" /><figcaption>FIG. 15. The normalized error $\tilde{\delta}_k$ (top) and the horizontal beta function $\beta_x$ (bottom) for ten irregular BPMs (blue circles) and ten regular BPMs (orange squares). The average value of $\tilde{\delta}_k$ has been extracted.</figcaption></figure>

$$ C = \beta c T_o, \qquad (40) $$

$$ s_k = \beta c t_k, \qquad (41) $$

where $c$ is the speed of light, $\beta$ is the Lorentz factor, and $C$ is the circumference of the ring. Combining the previous expressions with Eq. (39) results in

$$ \tilde{\delta}_k = \frac{s_k}{C} - \frac{k}{M}. \qquad (42) $$

The normalized error $\tilde{\delta}_k$ for the aforementioned BPM collections is shown in Fig. 15 (top), together with the model horizontal beta function (bottom) for the irregular (orange line) and the regular (blue line) BPM collections. The oscillation of the error is negligible for the regular case, whereas for the irregular collection the oscillation of the error is obvious. In addition, the horizontal beta function of the regular case seems almost constant, while the irregular case exhibits larger deviations from a constant value. For the irregular case, both the root-mean-square (rms) values of the error $\tilde{\delta}_k$ and of the horizontal beta function $\beta_x$ are almost 50 times larger than the respective rms values of the regular case.

The difference in the tune measurements of the two collections, relative to the tune measurements of the regular BPMs, are shown in Fig. 16, for the horizontal plane (blue) and the vertical plane (green) of the $4\sigma$ excitation. The trend of the curves exhibit a reduction of the relative difference, and, at $N = 30$ turns, it is around $10^{-6}$ for the horizontal and one order of magnitude larger for the vertical. The relative error converges at around $N = 50$ turns. As a result, precise tune measurements are possible even for a collection of BPMs with large asymmetry in the longitudinal position.

## V. EXPERIMENTAL APPLICATION OF TBT

Over the past years, various measurements in proton and electron storage rings have been undertaken, where the
<figure><img src="image_18.png" /></figure>

FIG. 16. The absolute difference between the tune measurements by mixing ten regular and ten irregular BPMs, normalized to the tunes for the ten mixed regular BPMs, with respect to the first 50 turns. The horizontal tune is shown in blue and the vertical in green. The vertical axis is in the logarithmic scale.

mixed BPM method has been used with success for precise tune determination over a small number of turns [16,17, 26,27]. Because of the unavoidable noise of the experimental TBT, collections of BPMs that exhibit small or no variation in the optics and/or the sampling rate do not offer any substantial improvement in the convergence. This observation agrees with the detailed measurements that were performed at the European Synchrotron Radiation Facility by using the mixed BPM method [16]. Some distinct cases of experimental measurements are presented in this section.

### A. PSB

At the PSB, TBT data are analyzed from beams that are excited by a single kick, horizontally and vertically, at 160 MeV. Recently, the PSB has been undergoing many upgrades in view of the LHC injectors upgrade (LIU) project [28]. The ability for precise tune measurements is very important for such a low-energy machine, where collective effects and instabilities can drive the betatron tunes on resonances. The interesting feature of the PSB lattice is that it has a 16-fold symmetry and 16 BPMs; i.e., one BPM is located at each superperiod, and the optics exhibit no modulations. In addition, the BPMs are equidistant, with an almost equal longitudinal distance of about $\Delta s = 10$ m. As a result, precise measurements could be performed with the mixed BPM method, in a few turns. The results of horizontal and vertical tune measurements, from the aforementioned TBT data, are presented in Figs. 17 and 18. Remarkably, the betatron tunes can be extracted with only ten turns, and the convergence at 20 turns is about $10^{-5}$. Noise had been filtered from the experimental TBT by using singular value decomposition (SVD) [29], which helps as well in the fast and precise measurement of the tunes.

### B. KARA

The KARA ring (previously named ANKA) is a very flexible electron light source, operating at 2.5 GeV. The ring

<figure><img src="image_19.png" /></figure>

FIG. 17. Tune measurements with respect to the number of turns *N* at the PSB ring 1 by using horizontal (top) and vertical (bottom) data.

<figure><img src="image_20.png" /></figure>

FIG. 18. Convergence of the tune measurements at the PSB ring 1 for the horizontal (top) and vertical (bottom) planes.

has a fourfold symmetry, and it is equipped with 39 BPMs. Recently, the prototype of the superconducting wiggler for the Compact Linear Collider (CLIC) damping rings [30,31] was installed at KARA in order to perform tests. In parallel, optics measurements have been performed with the wiggler at the maximum field (2.9 T). The measurements of the horizontal tune with the mixed BPM method are shown in Fig. 19, along with the statistical uncertainty. The method proved to be very efficient also in this case, since the horizontal tunes can be evaluated at around 20 turns, with a convergence of below $10^{-3}$. This measurement is important in order to quantify any quadrupolar effects of the wiggler on the horizontal tune. Indeed, the results show that a slight horizontal tune shift is present when the field of the wiggler changes from 0 to 2.9 T. The tune shift is at the order of about $\Delta Q = 4 \times 10^{-3}$.

### C. PS

The PS is one of the most indispensable parts of the LHC injector complex. At injection, the PS receives bunches from the PSB at a kinetic energy of 1.4 GeV, and then it
<figure><img src="image_21.png" /><figcaption>FIG. 19. Horizontal betatron tune measurements at the KARA light source under the operation of the CLIC superconducting wiggler. The case for the field of the wiggler at 2.9 T is shown in red, while the zero field case is shown in blue. The error bars are one standard deviation from the mean value of the measurements.</figcaption></figure>

accelerates them up to 14–26 GeV. During injection, TBT data are gathered from 43 BPMs around the ring for the mixed BPM method. In order to investigate the injection process, which involves a very fast injection bump [32] that lasts for about 500 turns, a scanning window of 40 turns is applied to the data. In this way, any variations of the *instantaneous* betatron tunes can be measured. The gating functionality of the BPM system at the PS is used, in order to perform the analysis bunch by bunch. A preanalysis is also performed to characterize the mixed BPM method at the PS. With only 40 turns, the convergence of the tunes is of the order of $10^{-4}$, which makes the method very precise. It should be noted that similar rates of convergence cannot be reached with the traditional PS frequency analysis tools for such a short time window.

The results of the scanning window analysis can be seen in Fig. 20, where the estimations of the horizontal (thick lines) and vertical (dashed lines) tunes for four bunches (bunch 1 in magenta, bunch 2 in red, bunch 3 in green, and bunch 4 in blue) coming from the PSB are shown. The observation of the periodic modulation of the tunes is evident. The depth of the modulation is such as to have a maximum tune shift at the order of $10^{-2}$ for both planes. Because of the simultaneous presence of the effect in both planes and the large horizontal amplitude of the beam in the injection area, the effect is attributed to feed-down effects. More specifically, persistent sextupolar fields are created on the surface of the vacuum chamber of the injection bumpers and during injection; they modulate the quadrupolar component of the machine, resulting in the visible TBT tune shift. This observation is significant for the LIU project [33], because such large tune shifts can result in particle losses and reduction of the machine's efficiency as an injector. This effect could have not been evaluated without the possibility to estimate the tune and its evolution in a small number of turns, as provided by the mixed BPM method.

<figure><img src="image_22.png" /><figcaption>FIG. 20. Instantaneous betatron tune measurements with the mixed BPM method, during the injection process at the PS. The estimation of the horizontal tunes is shown in thick lines and of the vertical tunes in dashed lines. The analysis is performed for four bunches (bunch 1 in magenta, bunch 2 in red, bunch 3 in green, and bunch 4 in blue) by using a sliding window of 40 turns.</figcaption></figure>

## VI. CONCLUSIONS

The efficiency of the mixed BPM scheme, when combined with NAFF, has been demonstrated for numerical, tracking, and experimental data. For all cases, the method has been proven to be extremely efficient for precise tune measurements with a very small number of turns (below 50).

Because of amplitude and phase modulations, which are caused from the periodicity of the optics and of the sampling error at the locations of the BPMs, the estimation of the betatron tunes can be affected, but results show that, as the number of turns increases, this contribution is reduced.

The mixed BPM method, which combines the TBT data from *M* BPMs simultaneously, cannot be used with a simple FFT due to the linear dependence of the frequency resolution to the number of samples. When a refined Fourier analysis method is used, e.g., NAFF with a Hann window of the order of *p*, an improvement of $M^{-(2p+1)}$ is derived in the analytical expressions. Because of the periodic change of the sampling period at each BPM, an error is introduced in the tune estimations which is negligible, and it is rapidly reduced.

Numerical simulations confirm the efficiency of the mixed BPM method for precise measurements of the tune *Q* within a small number of turns. The change of the sampling rate in the mixed BPM method results in the reduced tunes *Q*/*M*, and, when this fraction is close to a rational number, the error in the tune measurement is significantly increased. By implementing a varying frequency scheme, the expected dependence $M^{-(2p+1)}$ is recovered in these simulations.

The mixed BPM technique is also used for tracking data from the lattice of the PS machine. Comparisons with single-BPM analysis validate the improvement in accuracy and precision, even for BPMs with asymmetries in the optics and/or the longitudinal position. Depending on
magnitude of the modulation of the optics and the sampling rate of the M BPMs, a decrease of the precision in the tune measurements is observed for a very small number of turns; however, this discrepancy is reduced quite fast, as the number of turns increases.

In the case of experimental data, maximization of the number of BPMs is often preferred, rather than using fewer BPMs with symmetries in the optics and/or to the longitudinal position. The reason for this preference is the existence of noise in the TBT data, which reduces significantly the efficiency of both the single and mixed BPM methods. However, studies show that, by keeping the order of the Hann window as low as possible, the impact of noise in the tune estimations of the NAFF algorithm is reduced. Although noisy experimental TBT data can be sometimes impossible to be Fourier analyzed, the application of appropriate filtering algorithms (e.g., SVD) can reduce the impact of noise on frequency analysis. In such cases, the mixed BPM method has always proved to be very precise in the estimation of the betatron tunes.

Useful applications of the mixed BPM method have been demonstrated with data obtained from the PSB and PS proton rings, where the method has been proven capable of performing precise tune measurements with a very small number of turns. As a matter of fact, a large periodic tune shift has been discovered recently at the PS, during the injection process, which could not have been observed with the single-BPM analysis. The scheme of the mixed BPMs is also applied to data from the KARA light source, where betatron tune measurements were used for modeling the beam's response during operation of the CLIC superconducting wiggler.

The mixed BPM method has the potential to be a very useful tool, because it offers a substantial improvement in tune estimations for a very small number of turns. In addition, future studies could investigate improvements on the method, such as the interpolation of nonuniformly sampled signals [34], in order to minimize the contribution of the various periodic modulations at the position of the BPMs.

Finally, it should be noted that, in cases of tracking simulations, the mixed BPM method can be employed for many applications that require small number of turns, e.g., for the construction of fast frequency maps or sliding windows in the TBT data for measuring the instantaneous tunes. In cases of the absence of noise from the data, schemes with BPMs that are equidistant and/or periodic in the optics should be preferred, because they can provide even faster convergence after the first few turns.

## ACKNOWLEDGMENTS

We acknowledge the PS and PSB OP teams and the KARA operators for being extremely helpful with the acquisition of the data. This work has been supported by the CLIC project and the low emittance rings (LER) network of EUCARD2 and RULE of ARIES.

# APPENDIX: FREQUENCY RESOLUTION OF THE MIXED BPM METHOD WITH NAFF

The main success of NAFF is the ability to provide an approximation of a Kolmogorov-Arnold-Moser [35] quasiperiodic function $f(t)$, which can be a solution of a Hamiltonian system. The function $f(t)$ can be approximated as

$$f(t) \approx \sum_{k=1}^{n} \alpha_k e^{i\omega_k t} \quad (A1)$$

with $n$ the number of terms and $\alpha_k$ and $\omega_k$ the complex amplitude and real angular frequency of the $k$th harmonic, respectively.

Physically, the function $f(t)$ can be a trajectory of a nonlinear quasiperiodic system, similar to the oscillation of the beam's centroid, along the lattice of an accelerator. The NAFF algorithm can recover the fundamental eigenfunctions of $f(t)$ in a very short time. For instance, the error in the estimation of the main frequency $\nu_1 = \frac{\omega_1}{2\pi}$ of Eq. (A1) has been mathematically proven by Laskar to be [36]

$$|\nu(T) - \nu_1| = \frac{C_L}{T^{2p+2}} + O\left(\frac{1}{T^{2p+2}}\right), \quad (A2)$$

for $T \to \infty$, where $T$ is the total observation time, $\nu(T)$ is the time-dependent frequency estimation of NAFF for the main harmonic, $p$ is the order of the Hann window, and the constant $C_L$ given by

$$C_L = c_0 \sum_{\vec{k}=(1,0,\dots,0)}^{n} \frac{\mathrm{Re}\{\alpha_k\}}{\Omega_k^{2p+1}} \cos(\Omega_k T), \quad (A3)$$

in units of $s^{-(2p+1)}$, with $\Omega_k = \langle \vec{k}, \vec{\nu} \rangle - \nu_1$, $\vec{k} = \{k_1, k_2, \dots, k_n\}$ the basis vector and $\vec{\nu} = \{\nu_1, \nu_2, \dots, \nu_n\}$ the frequency vector. Their inner product $\langle \vec{k}, \vec{\nu} \rangle = k_1\nu_1 + k_2\nu_2 + \dots + k_n\nu_n$ forms the resonance space of the quasiperiodic solution. The symbol $\vec{k} - (1, 0, \dots, 0)$ means that, from the summation of the harmonics, the first-order resonance $\{k_1, k_2, \dots, k_n\} = \{1, 0, \dots, 0\}$ is excluded. The term $\mathrm{Re}\{\alpha_k\}$ the real part of the complex Fourier amplitude $\alpha_k$ of Eq. (A1), in units of amplitude, and the dimensionless constant

$$c_0 = \frac{(-1)^{p+1} \pi^{2p} (p!)^2}{|\alpha_1| |\phi''(0)|}, \quad (A4)$$

with $|\alpha_1|$ the absolute value of the complex amplitude of the first harmonic for $k=1$, and the constant

$$\phi''(0) = -\frac{2}{\pi^2} \left( \frac{\pi^2}{6} - \sum_{k=1}^{p} \frac{1}{k^2} \right). \quad (A5)$$
In the original derivation of Eqs. (A2) and (A3) in Ref. [36], the normalization to the amplitude of the first harmonic $|\alpha_1|$ is implied. Note that, in the case of rational frequencies $\vec{\nu}$, the quasiperiodic orbits cannot be defined and $\Omega_k \to 0$, which makes Eq. (A3) diverge rapidly. In addition, the existence of multiple harmonics in the signal Eq. (A1) can reduce the rate of convergence to the actual frequencies as well. On the contrary, if the signal $f(t)$ is composed of only one periodic term, i.e., only one harmonic, the error $|\nu(T) - \nu_1|$ vanishes, since the summation term of Eq. (A3) goes to zero.

In the practical case of frequency analysis of a quasi-periodic signal, where $m$ samples are gathered with a uniform sampling period $\tau_s$, the total observation time is $T = m\tau_s$. As a result, Eq. (A2) is modified as

$$|\nu(m) - \nu| = \frac{C_L}{\tau_s m^{2p+2}}, \qquad (A6)$$

and Eq. (A3) becomes dimensionless, with the form

$$C_L = c_0 \sum_{k-(1,0...0)}^{n} \frac{\mathrm{Re}\{\alpha_k\}}{(\tau_s \Omega_k)^{2p+1}} \cos(\Omega_k m \tau_s). \qquad (A7)$$

Next, the mixed BPM scheme is explored with the NAFF algorithm. Under this transformation, in the case of $M$ BPMs that record $N$ turns of the betatron oscillation of the beam, the total observation time is (see Sec. II A)

$$T = m \frac{T_o}{M} + \delta_M, \qquad (A8)$$

for $m$ number of samples. The term $T_o$ is the revolution period of the beam, and $\delta_M$ is the deviation of the position of the last BPM from a fictitious position, which would be symmetric for all the $M$ BPMs. The previous expression can be used for estimating the error in the betatron tune measurement with NAFF.

The error in the betatron tune estimation $\epsilon(m)$, with respect to the error in the frequencies $\Delta\nu(m) = |\nu(m) - \nu|$, for a total of $m$ samples is defined as

$$\epsilon(m) = \Delta\nu(m)T_o. \qquad (A9)$$

The total time in Eq. (A8) can be used to define a varying sampling period. Since the total observation time is $T = m\tau_s$ and from Eq. (A8)

$$T = m \frac{T_o}{M} g(m), \qquad (A10)$$

where $g(m) = 1 + \frac{M\tilde{\delta}_M}{m}$ and $\tilde{\delta}_M$ is the error of the last BPM normalized to the revolution period, the sampling period $\tau_s$ can be expressed as

$$\tau_s = \frac{T_o}{M} g(m). \qquad (A11)$$

Substitution of the previous expression to Eq. (A6) and combining the result with Eq. (A9), the error $\epsilon(m)$ in the betatron tune estimation for the mixed BPM method is written as

$$\epsilon(m) = \frac{M}{m^{2p+2} g(m)^{2p+2}} \overline{C_L}(m), \qquad (A12)$$

where the $\overline{C_L}(m)$ factor is

$$\overline{C_L}(m) = c_0 \sum_{k-(1,0...0)}^{n} \frac{1}{\Omega_k^{2p+1}} \mathrm{Re}(\alpha_k) \cos[mg(m)\overline{\Omega_k}] \qquad (A13)$$

and the term $\overline{\Omega_k}$ is

$$\overline{\Omega_k} = \frac{T_o \Omega_k}{M}. \qquad (A14)$$

In the case of $m = MN$ samples, the error $g(MN)$ is

$$g(MN) = 1 + \frac{\tilde{\delta}_M}{N}, \qquad (A15)$$

and Eq. (A12) becomes

$$\epsilon(MN) = \frac{\overline{C_L}(MN)}{M^{2p+1} N^{2p+2}} \left(1 + \frac{\tilde{\delta}_M}{N}\right)^{-(2p+2)}. \qquad (A16)$$

Since the error term $\tilde{\delta}_M \ll N$, the previous expression is expanded around $\tilde{\delta}_M \approx 0$, which yields

$$\epsilon(MN) = \frac{\overline{C_L}(MN)}{M^{2p+1}} \left( \frac{1}{N^{2p+2}} - (2p+2) \frac{\tilde{\delta}_M}{N^{2p+3}} \right). \qquad (A17)$$

The expression in Eq. (A17) testifies that the convergence is improved by a factor of $M^{2p+1}$ with the mixed BPM method. The contribution of the small error $\tilde{\delta}_M$ converges to zero rapidly enough, so as to be negligible in the frequency analysis with NAFF and the mixed BPM method.

[1] M. Serio, Transverse betatron tune measurements, in *Frontiers of Particle Beams; Observation, Diagnosis and Correction*, edited by M. Month and S. Turner (Springer, Berlin, 1989), pp. 65–93.

[2] R. Bartolini, M. Giovannozzi, W. Scandale, A. Bazzani, and E. Todesco, Precise measurement of the betatron tune, Part. Accel. **55**, 1 (1996).

[3] R. E. Meller, A. W. Chao, J. M. Peterson, S. G. Peggs, and M. Furman, Decoherence of kicked beams, Superconducting Super Collider Laboratory, Report No. SSC-N-360, 1987.

[4] S. Y. Lee, Decoherence of the kicked beams II, Indiana University, Report No. SSCL-N-749, 1991.
[5] J. W. Cooley and J. W. Tukey, An algorithm for the machine calculation of complex Fourier series, *Math. Comput.* **19**, 297 (1965).

[6] E. Asséo, J. Bengtsson, and M. Chanel, Absolute and High Precision Measurements of Particle Beam Parameters at CERN Antiproton Storage Ring LEAR Using Spectral Analysis with Correction Algorithms, in *Proceedings of the 4th European Signal Processing Conference, Saint Martin d'Hères, France* (1988), pp. 1317-1320.

[7] M. Gasior and J. L. Gonzalez, AIP Conf. Proc. **732**, 276 (2004).

[8] J. Laskar, C. Froeschle, and A. Celletti, The measure of chaos by the numerical analysis of the fundamental frequencies. Application to the standard mapping, *Physica (Amsterdam)* **56D**, 253 (1992).

[9] J. Laskar, Frequency map analysis and particle accelerators, *Proceedings of the 2003 Particle Accelerator Conference, Portland, OR, Conf. Proc. C030512*, 378 (2003).

[10] C. Skokos, L. Farvacque, Y. Papaphilippou, E. Plouviez, J. L. Revol, A. Ropert, and J. Laskar, Experimental frequency maps for the ESRF storage ring, in *Proceedings of the 9th European Particle Accelerator Conference, Lucerne, 2004* (EPS-AG, Lucerne, 2004).

[11] Y. Papaphilippou, Detecting chaos in particle accelerators through the frequency map analysis method, *Chaos* **24**, 024412 (2014).

[12] P. Zisopoulos, F. Antoniou, Y. Papaphilippou, A. Streun, and V. Ziemann, Frequency maps analysis of tracking and experimental data for the SLS storage ring, in *Proceedings of the 5th International Particle Accelerator Conference (IPAC 2014): Dresden, Germany* (JACoW, Geneva, Switzerland, 2014), THPRO076, https://doi.org/10.18429/JACoW-IPAC2014-THPRO076.

[13] F. J. Harris, On the use of windows for harmonic analysis with the discrete fourier transform, Proc. IEEE **66**, 51 (1978).

[14] J. Laskar, Frequency analysis for multi-dimensional systems. Global dynamics and diffusion, *Physica (Amsterdam)* **67D**, 257 (1993).

[15] Y. Alexahin, E. Gianfelice-Wendt, and W. Marsh (Proton Driver Collaboration), Tune evaluation from phased BPM turn-by-turn data, Conf. Proc. C100523, MOPE084 (2010).

[16] C. Skokos, J. Laskar, and Y. Papaphilippou, Precise tune measurements from multiple beam position monitors, Conf. Proc. C070625, 3913 (2007).

[17] P. Zisopoulos, F. Antoniou, E. Hertle, A.-S. Müller, and Y. Papaphilippou, Transverse tunes determination from mixed BPM data, in *Proceedings of the 6th International Particle Accelerator Conference (IPAC 2015): Richmond, Virginia* (JACoW, Geneva, Switzerland, 2015), TUPJE042, https://doi.org/10.18429/JACoW-IPAC2015-TUPJE042.

[18] J. Yen, On nonuniform sampling of bandwidth-limited signals, *IRE Transactions on Circuit Theory* **3**, 251 (1956).

[19] A. V. Oppenheim and R. W. Schafer, *Discrete-Time Signal Processing*, 3rd ed. (Pearson, Upper Saddle River, NJ, 2010).

[20] C. E. Shannon, Communication in the presence of noise, Proc. IRE **37**, 10 (1949).

[21] J. Moser, On the theory of quasiperiodic motions, *SIAM Rev.* **8**, 145 (1966).

[22] F. Asvesta, N. Karastathis, and P. Zisopoulos, PyNAFF: A (C)Python module that implements NAFF, https://github.com/nkarast/PyNAFF.

[23] S. Kostoglou, N. Karastathis, Y. Papaphilippou, D. Pellegrini, and P. Zisopoulos, Development of computational tools for noise studies in the LHC, in *Proceedings of the 8th International Particle Accelerator Conference (IPAC 2017): Copenhagen, Denmark* (JACoW, Geneva, Switzerland, 2017), THPAB044, https://doi.org/10.18429/JACoW-IPAC2017-THPAB044.

[24] F. Schmidt, E. Forest, and E. McIntosh, Introduction to the polymorphic tracking code: Fibre bundles, polymorphic taylor types and exact tracking (unpublished).

[25] A. Blas et al., The PS complex as proton pre-injector for the LHC: Design and implementation report, 2000.

[26] P. Zisopoulos and M. Gąsior, M. Serluca, and G. Sterbini, Fast bunch by bunch tune measurements at the CERN PS, in *Proceedings of the International Particle Accelerator Conference (IPAC 2017): Copenhagen, Denmark* (JACoW, Geneva, Switzerland, 2017), MOPAB122, https://doi.org/10.18429/JACoW-IPAC2017-MOPAB122.

[27] M. Serluca et al., The (7,7) optics at CERN PS, CERN Proc. **2**, 43 (2017).

[28] K. Hanke et al., Status and plans for the Upgrade of the CERN PS booster, in *Proceedings of the 6th International Particle Accelerator Conference (IPAC 2015): Richmond, Virginia* (JACoW, Geneva, Switzerland, 2015), THPF090, https://doi.org/10.18429/JACoW-IPAC2015-THPF090.

[29] C.-x. Wang, Model independent analysis of beam centroid dynamics in accelerators, Ph.D. thesis, Stanford University, Physics Department, 1999.

[30] Y. Papaphilippou et al., Conceptual design of the CLIC damping rings, in *Proceedings of the 3rd International Particle Accelerator Conference, New Orleans, LA, 2012* (IEEE, Piscataway, NJ, 2012), Vol. C1205201, pp. 1368-1370.

[31] A. Bernhard et al., A CLIC damping wiggler prototype at ANKA: Commissioning and preparations for a beam dynamics experimental program, in *Proceedings of the International Particle Accelerator Conference (IPAC 2016): Busan, Korea* (JACoW, Geneva, Switzerland, 2016), WEPMW002, http://dx.doi.org/10.18429/JACoW-IPAC2016-WEPMW002.

[32] H. Serluca, W. Bartmann, V. Forte, M. Fraser, and G. Sterbini, Injection septa position and angle optimisation in view of the 2 GeV liu upgrade of the CERN PS, CERN Proc. **2**, 51 (2017).

[33] M. Meddahi et al., LHC injectors upgrade (LIU) project at CERN, in *Proceedings of the 6th International Particle Accelerator Conference (IPAC 2015): Richmond, Virginia* (JACoW, Geneva, Switzerland, 2015), THPF093, https://doi.org/10.18429/JACoW-IPAC2015-THPF093.

[34] S. Maymon and A. Oppenheim, Sinc interpolation of nonuniform samples, *IEEE Trans. Signal Process.* **59**, 4745 (2011).

[35] H. S. Dumas, *The KAM Story: A Friendly Introduction to the Content, History, and Significance of Classical Kolmogorov-Arnold-Moser Theory* (World Scientific, Hackensack, NJ, 2014).

[36] J. Laskar, Frequency map analysis and quasiperiodic decompositions, arXiv:math/0305364.