Analytical N beam position monitor method

Measurement and correction of focusing errors is of great importance for performance and machine protection of circular accelerators. Furthermore LHC needs to provide equal luminosities to the experiments ATLAS and CMS. High demands are also set on the speed of the optics commissioning, as the foreseen operation with β -leveling on luminosity will require many operational optics. A fast measurement of the β-function around a storage ring is usually done by using the measured phase advance between three consecutive beam position monitors (BPMs). A recent extension of this established technique, called the N-BPM method, was successfully applied for optics measurements at CERN, ALBA, and ESRF. We present here an improved algorithm that uses analytical calculations for both random and systematic errors and takes into account the presence of quadrupole, sextupole, and BPM misalignments, in addition to quadrupolar field errors. This new scheme, called the analytical N-BPM method, is much faster, further improves the measurement accuracy, and is applicable to very pushed beam optics where the existing numerical N-BPM method tends to fail.


I. INTRODUCTION
In recent years the field of optics measurement and correction is growing in interest with the design of pushed optics like the high-luminosity LHC (HL-LHC) upgrade and next generation light sources.A review of the progress in the last years is given in [1].The original method to determine the β function from turn-by-turn phase data makes use of three adjacent beam position monitors (BPMs) [2].The betatron phase of the BPMs is derived from the harmonic analysis.The phase advance between BPMs is then used to calculate β-functions according to the formula [2] where ϕ ij ¼ ϕðs j Þ − ϕðs i Þ are the measured phase advances between BPMs j and i; the superscript (m) denotes model values.The phase measurement is independent of BPM calibration and transverse misalignments.This method is sensitive to the position of the BPMs relative to each other.If the phase advance between two of them is too close to nπ, errors in the phase measurement get strongly enhanced.In order to measure the transversal displacement of the beam as accurately as possible, in the LHC the beam is excited by an AC-dipole [3].

A. Original N-BPM method
To avoid such cases with unsuitable phase advances and to improve statistics, BPMs can be skipped and more combinations can be used and averaged over with appropriate weights.The N-BPM method was developed [4] to implement this feature.It is illustrated in Fig. 1.If we use a range of N BPMs, where the probed BPM is fixed at position s i , there are n ¼ ðN − 2ÞðN − 1Þ=2 combinations of BPMs.The N-BPM method was successfully used in the LHC [4], as well as in the storage rings of ALBA [5] and ESRF [6].In comparison to the original three BPM method, there is a huge gain in precision in the interaction regions (IRs) of colliders, where neighboring BPMs have unsuitable phase advances with respect to each other.Systematic errors from the model and random phase uncertainties are taken into account separately where only statistical uncertainties are calculated analytically.Systematic errors are determined using Monte Carlo simulations.
The β function of the lth combination reads Unfortunately, the more lattice elements there are between BPMs the more sources of errors may lie between them, Eqs. ( 1) and ( 2) being derived under the assumption of having error-free regions between BPMs.Furthermore the different results for the β function at one location are not independent of each other which also has an impact on the quality of the final result.
Especially in the LHC IRs, where the phase advance between the BPMs is close to 0 or π, the three BPM method fails.At the same time a high precision of the measurement of beta beating is needed to control β Ã and to optimize the aperture.
To find the best estimation for the measured β we calculate a weighted average where the weights g l are calculated using a least-squares estimation.The sum runs over all n BPM combinations (details are worked out in Appendix A).The weights are where V denotes the covariance matrix of the β l .The uncertainty of the averaged β is then given by β is a function of the measured phase advances which are subject to measurement uncertainties.We can get the error matrix V from the single variances by Þ is a diagonal matrix consisting of the variances of the phases and T is the Jacobian matrix δϕ ¼ 0 meaning that the derivatives are evaluated with all phase advance errors set to zero.The correlation of statistical errors (coming from BPM noise) is We place here intentionally a superscript ϕ to highlight that the T matrix only includes statistical errors coming from the uncertainty of the phase measurement.δ α β denotes the Kronecker δ.Including systematic errors, the total error matrix is where V stat ¼ T ϕ MT ϕ−1 and the total uncertainty of the averaged β is then given by In the existing numerical N-BPM method only the statistical errors are calculated analytically while V syst and hence σ 2 syst are evaluated from Monte Carlo simulations of lattices with errors.For that the statistics over many simulations is gathered.Even with a highly parallelized code on a multicore machine this procedure takes about 1 hour for a "fast" set of 1000 simulations.Since the computation time scales linearly with the number of simulations, this can take up to 10 hours for 10 4 simulations which were used in the post processing of a measurement.

B. Extension of the N-BPM method
In this paper we introduce a new method that calculates also the systematic errors analytically.On the same computer the analytical N-BPM method takes only 30 seconds to compute the β function for more than 500 BPMs whereas the N-BPM method takes about one hour.It provides a fully analytical calculation of the uncertainties while the precision of the original method depends on the number of simulations.Any source of uncertainties can be taken into account if its analytical expression is known.The method does not depend on the stability of the lattice whereas the Monte Carlo simulations fail if, for some combinations of errors, no closed orbit can be found.This is a limiting factor when the N-BPM method is used for pushed optics with very low β Ã like the HL-LHC.
Equation (2) will be modified in Sec.II for taking into account the presence of quadrupolar field errors as done in [7], as well as transverse misalignments of sextupoles and longitudinal displacement of quadrupoles and BPMs.In FIG. 1. TOP: the 3 BPM method takes only adjacent BPMs.The blue BPM is the one whose β function is being determined.BOTTOM: The N-BPM method allows to skip BPMs in order to use a bigger amount of data to increase statistics and to avoid unsuitable phase advances.
Sec. III we then incorporate these results in the development of a fully analytical N-BPM method which does not require the splitting of the statistical phase uncertainties and the systematic errors.This was a limitation of the N-BPM method described in [4].We eventually compare the methods through simulations for the current LHC lattice design as well as the Achromatic Telescopic Squeeze (ATS) lattice proposed for the HL-LHC upgrade.

II. β FUNCTION FORMULA WITH IMPERFECTIONS
Equation (2) assumes that no error is present in the region between the involved BPMs.A new formula has been developed in [7] that takes quadrupolar errors into account, The sum runs over all elements w between BPMs i and j. δK w;1 is the integrated quadrupolar field error of element w.
The definition of K n follows the MAD [8] convention.Note that Eq. ( 12) is only defined for the case s 1 < s 2 < s 3 i.e. the β function is calculated at the position of the first BPM.
Since we want to use as many BPM combinations as possible we will derive below a form of Eq. ( 11) without this constraint.In [7] the quality of Eq. ( 11) has been assessed for the ESRF storage ring by simulating a lattice with errors and comparing the results of Eq. ( 11) to the beta beating of the simulated lattice.In this paper the same study is repeated for the LHC and its future upgrade the HL-LHC.Table I summarizes those uncertainty estimates for the LHC elements that are relevant for this study.In Fig. 2 a comparison between Eqs. ( 1) and ( 11) is shown in the horizontal plane.The plot shows a histogram of the deviation of the formulas (1) and (11) from the real (i.e., simulated) β function values.The histogram contains the values of all the individual BPMs for one simulated lattice.The top plot shows the case with only quadrupolar field errors whereas in the bottom plot also misalignment errors were included in the lattice.The introduction of errors not taken into account in Eq. ( 12) deteriorates the result.However Eq. ( 11) still yields a clearly better estimate.
Moreover Eq. ( 11) may be easily modified for taking into account a more realistic set of errors-quadrupolar field errors as well as sextupole transverse misalignments, quadrupole longitudinal misalignments, and BPM longitudinal displacements.
While the magnet misalignment errors can be approximated as effective quadrupolar field errors and integrated in Eq. ( 12) the BPM misalignments need a different approach as shown in the next paragraph.[9,10].
2. Top: Difference between the real (simulated) horizontal β-functions and the ones calculated by Eq. ( 1) in red and ( 11) in blue, respectively.Data are from MADX simulations of a lattice with quadrupolar errors only.Bottom: The same quantities are evaluated for the case with additional magnet misalignments.rms denotes the root mean square of the deviations and p-t-p denotes peak-to-peak deviation.

A. Effect of transverse sextupole misalignments
The magnetic field of a sextupole displaced horizontally by Δx reads This induces a quadrupolar field error whose strength δK 1 is This term can be used in Eq. ( 12) to include sextupole offsets in hij .

B. Effect of longitudinal quadrupole misalignments
The effect of a longitudinal displacement δs of a quadrupole magnet can be approximated by leaving the magnet at its original position and introducing two thin magnets at its edges to mimic the displacement, as shown in Fig. 3.In the direction of the displacement there is an additive element with integrated field strength δK 1 ¼ k 1 δs (k 1 being the nonintegrated quadrupole strength), whereas an error −δK 1 is placed at the opposite end.

C. Effect of BPM misalignments
An error, δs i , in the longitudinal position of a BPM affects the evaluation of Eqs.(11) and (12) which rely on the model values of β and ϕ at the nominal position of the BPM.To determine the effect we start with the definition of the phase advance We can approximate the phase error and the resulting β shift at the position up to first order in δs i .By α i we denote the α function at the position of element i, defined as We have to rederive an equation similar to Eq. ( 11) by taking into account the considerations of the preceding sections.The steps of the derivation are elaborated in Appendix B. The final formula up to first order reads: where δK w;1 now includes quadrupolar-like errors coming from sextupole misalignments and quadrupole longitudinal misalignments, as described in the previous sections.
Having defined the set I as so that an element with index w ∈ I lies between elements i and j, Eqs. ( 19) and (20) hold for every combination i, j, k of the BPMs.By doing so, we do not need to distinguish the three cases where the probed BPM is in the middle, left or right.All these considerations can be put into Eq.( 19) and used to get a more accurate β function.To verify the validity of Eq. ( 19), its horizontal β-functions are compared to the ones simulated by MADX along with the ones inferred from Eq. (1), this time including sextupole radial offsets and BPMs longitudinal shifts.
The result is shown in Fig. 4. The accuracy is now as good as the one of Eq. ( 11) when only quadrupolar field errors were introduced in the lattice, which in turn is much greater than the old formula, Eq. (1).FIG. 3. The top sketch shows the displaced quadrupole (solid gray) relative to the original position (dashed).In the bottom sketch one can see the quadrupole at its original position with thin magnets on both ends.

III. THE ANALYTICAL N-BPM METHOD A. Calculation of the correlation matrix
The Jacobian T of Eq. ( 7) can be split into blocks for the uncertainties of phase T ϕ , quadrupole field T K and BPM misalignment T s .T ϕ is the same of Eq. ( 8).For the quadrupolar field errors we get with : The contribution from the BPM misalignment is calculated analogously: In Eqs. ( 23) and ( 25) the minus and plus signs refer to the horizontal and vertical case, respectively.

B. Removal of bad BPM combinations
Since a phase advance ϕ ij ≈ nπ results in an enhancement of phase measurement errors and in the extreme case numerically unstable values, a filtering was introduced.Instead of keeping a constant number of combinations as in [4] we set a threshold for bad phase advances.A phase advance Δϕ is considered bad if Δϕ ∈ ½nπ − δ; nπ þ δ for n ∈ N and a given threshold δ.If any of the four phase advances ϕ ij l ; ϕ ik l ; ϕ m ij l ; ϕ m ik l in Eq. ( 2) is bad, the corresponding BPM combination is disregarded in the calculation of the weighted mean.This allows us to still use several combinations but skipping those which are numerically unstable.The current value for the threshold is δ ¼ 2π × 10 −2 .The use of fewer combinations results in a lower computation time.
To test the analytical N-BPM method and compare it to the original 3-BPM and the Monte Carlo N-BPM method, a large set of LHC lattices with β Ã ¼ 40 cm with randomly distributed errors is generated and a measurement is simulated by tracking a single particle via polymorphic FIG. 4. Accuracy of the horizontal β-function evaluated via Eq.( 1) and ( 19) with the effect of magnets and BPM misalignments taken into account.The accuracy of Eq. ( 19) is similar to the one of Eq. ( 11) with quadrupolar field errors only (Fig. 2).FIG. 5. Comparison of the 3-BPM method, the original N-BPM method and the analytical N-BPM method (denoted as A.N-BPM) for the nominal LHC lattice at collision, with β Ã ¼ 40 cm.Bottom: histogram of the difference to the real β-function in percent.Top: the average of the error bars and the accuracy spread (width of a standard distribution fit to the distribution in the bottom plot) in percent.The analytical N-BPM method has the best accuracy both in the arcs and in the IRs.Data have been cleaned of outliers.
tracking code (PTC) [11].The random errors are created from Table I and a Gaussian noise of σ x ¼ 0.1 mm is applied to the BPM signal.No singular value decomposition cleaning is applied since it would clean the artificial noise too efficiently [5].The excitation amplitude is 0.8 mm at a β function of about 120 m.The tracked particle positions are then analyzed by the three methods (3-BPM, N-BPM, and analytical N-BPM) and the respective deviation from the real horizontal β function is shown in the bottom plot of Fig. 5.The analytical N-BPM method includes the filtering of phase advances.Especially in the IR, where neighbouring BPMs have often unsuitable phase advances, the N-BPM and analytical N-BPM method yield more accurate values.
The top diagram of Fig. 5 shows that the error of the 3-BPM method is considerable larger, whereas the N-BPM and analytical N-BPM method are very accurate with similar accuracies in the IR and arcs.

C. HL-LHC
The ATS optics [12] is the baseline choice for the HL-LHC and our optics measurement tools have to be prepared for the challenges imposed by such an optics.In Fig. 6 the three methods are compared in the same way as in Fig. 5.The excitation amplitude was 0.8 mm at a β function of 127 m.For the current HL-LHC collision optics (β Ã ¼ 15 cm) the performance of N-BPM and analytical N-BPM method is again better than the 3-BPM method, especially in the IRs.All three methods are, however, about a factor two more inaccurate than for the β Ã ¼ 40 cm optics of LHC, in agreement with Fig. 7 of [4].
In the post-processing of the data taken during the LHC Machine Development measurement (MD) [13] for testing the ATS principle with a β Ã ¼ 10 cm optics, the analytical N-BPM method was used for the first time with filtering of bad phase advances.Figure 7 demonstrates that the analytical N-BPM method deals well with the ATS MD optics.Monte Carlo simulations were not possible for this optics.
Figure 8 shows the precision of the final results for the β Ã ¼ 10 cm optics of both beams compared to the simulations of the HL-LHC lattice with β Ã ¼ 15 cm.To ease the comparison the error bars are shown in the bottom plot.They are slightly larger for the real measurement than those in simulations.We believe that the use of lower beam excitation to ensure machine protection is behind these larger error bars.Figure 7 shows also that the 3-BPM method has many outliers and error bars up to several kilometers caused by bad phase advances.Large error bars have been excluded for the mean shown in Fig. 8.
The Monte Carlo simulations failed for low β Ã optics and so we were not able to use the original N-BPM method during ATS MDs.This is another advantage of the FIG.comparison of Fig. 5 between the three methods for the HL-LHC β Ã ¼ 15 cm ATS optics.The analytical N-BPM method yields clearly better results, both in the IRs and in the arcs.Compared to the β Ã ¼ 40 cm optics shown in Fig. 5, the β function reconstruction is less accurate: This was also demonstrated for the β Ã ¼ 20 cm optics in [4].FIG. 7. Horizontal β beating of beam 1 at β Ã ¼ 10 cm during the ATS MD 2016.Top: 3-BPM method.Bottom: analytical N-BPM method.The 3-BPM method suffers from bad phase advances and has many outliers.The regions of high β beating at around 8000 m and 23000 m lie in the telescopic arcs.
analytical N-BPM method that it is able to evaluate the systematic errors independently of the success of particle tracking.

IV. CONCLUSION
A new method for the measurement of β and α functions has been developed based on the existing N-BPM method.A fully analytical calculation of the covariance matrix provides a faster and more accurate measurement of β and α functions.The analytical N-BPM method also avoids the complications from failing to find closed optics that occur in the Monte Carlo simulations needed by the existing N-BPM method.This stability with respect to the choice of optics model makes it more suitable for low β Ã optics.Simulations show that, together with a filtering of BPM combinations according to the phase advances, the method is optimal for the HL-LHC upgrade.

ACKNOWLEDGMENTS
Many thanks go to the LHC Optics Measurement and Correction (OMC) team, especially T. Persson for fruitful discussions and J. Coello de Portugal for introducing one of the authors into the Python programming language.

APPENDIX A: LEAST-SQUARES ESTIMATION
We search for the weights g i of Eqs. ( 3)-( 4) through a least-squares optimization.We can rewrite the weighted mean, ⃗ β, as system of linear equations where ϵ i is the error of measurement β i , i.e. the difference to the weighted mean.ĝi are the weights and n is the number of measured β values.We seek a set of weights g i for which the squared errors in Eq. (A1) are minimal.Since the different β i are correlated and have a covariance matrix Cov½ ⃗ β ≠ diagðσ 2 1 ; …; σ 2 n Þ we have to apply the theory of generalized least-squares estimation [14].
The covariance matrix of a set of random variables ω i is defined as where E½ ⃗ω is the expected value of the random variables ⃗ω.Equation (A2) can be expressed component-wise:

Error propagation
We have a set of unperturbed phase measurements fϕ 1 ;…; ϕ n g and nonobservable parameters fK 1;1 ; …; K 1;m ; s 1 ; …; s ν ; x 1 ; …x μ g with corresponding errors Δϕ α ; ΔK 1;β ; Δs γ ; Δx κ .We collect all parameters into a vector Using a Taylor expansion we can propagate the errors in these parameters where the Einstein summation convention is used and the derivatives are with respect to Ω i : ∂ i β l ≡ ∂β l ∂Ω i .In the following we omit the argument ð ⃗ Ω 0 Þ for the β function and its derivatives at the unperturbed position ⃗ Ω 0 .For the calculation of the covariance matrix, we truncate Eq. (A5) to second order and derive the two summands of Eq. (A2) separately.
We assume that there are no systematic errors in the measurement instruments and so E½ΔΩ i ¼ 0 and the middle term vanishes.We can conclude where terms proportional to the square of the variance go into the remainder OðΔΩ 3 Þ.
To calculate the term E½β l β m more steps are needed: Finally, subtracting Eq. (A7) from Eq. (A8), we obtain In the last step we made again the assumption In matrix notation Eq. (A9) reads where If the parameters ⃗ Ω are uncorrelated, the covariance matrix of ⃗ ΔΩ simplifies to and Eq.(A9) reads with the variance being defined by 2. The generalized least-squares estimator After computing the covariance matrix of β, we can calculate the generalized least-squares estimator by solving the translated system where the translation X is such that X T X ¼ V ¼ Cov½ ⃗ β.The least-squares estimation searches for a minimum of We solve Eq. (A19) in the index notation: WEGSCHEIDER, LANGNER, TOMÁS, and FRANCHI PHYS.REV.ACCEL.BEAMS 20, 111002 (2017) 111002-8 where we renamed the summation index in the numerator to highlight a solution of the system: The uncertainty of the weighted average is then APPENDIX B: β FUNCTION FROM PHASE WITH QUADRUPOLE AND MISALIGNMENT ERRORS In this Appendix we briefly present the derivation of Eqs. ( 19) and (20).It is basically a repetition of Appendix A of [7] without the constraint i < s j < s k and with BPM longitudinal misalignments from Eqs. ( 16) and ( 17) included.The derivation is based on the resonance driving terms (RDTs) formalism of [15] and references therein.As in [7] it is assumed that coupling RDTs are negligible, tune lines from nonlinear RDTs are small and Hamiltonian octupolar-like terms can be neglected.The important RDTs and Hamiltonian terms for the derivations are for the horizontal plane.It is understood that in the above equations the β function and phase are the horizontal ones.
For the vertical plane the RDT f 0200;j and the Hamiltonian term h 0011 shall be used while both terms have the opposite sign.For the sake of brevity the vertical plane will not be regarded in this derivation.The sums run over all W quadrupole errors (including misalignments that act as quadrupolar field errors).The indices jklm of the RDTs f jklm;j will be suppressed in the following discussion.Under the assumptions of Eq. ( 16), the RDT f j ≡ f 2000;j behaves as follows Subtracting two Hamiltonian terms i and j yields We can use the ratio of the components m 11 and m 12 of the transport matrix Mðs i ; s j Þ (see [2], chapter 2) to calculate the β function: From [7] the β and α functions and the phase advance are calculated as: with h ij ≡ h 1100;ij from Eqs. (B3) and (B4).Introducing the BPM misalignments as described in Eqs. ( 16), ( 17) changes the equations above to: up to first order in δs i .f j can be expressed in terms of f i , as developed in [16][17][18] Thus, we can eliminate the f i dependence in the second line of Eq. (B7) Before putting everything into Eq.(B5), we have to expand cot We can still simplify 2h ij − 4R½ ĥij e 2iϕ m j a bit: Combining Eqs.(B7), (B10) and the result for ḡij into (B5) yields for the horizontal plane.We note here that for the vertical plane ḡij would have the opposite sign and β and ϕ would of course be the vertical β function and phase.For a fixed probed BPM i we can exchange j → k in Eq. (B14).We can solve for β i by taking the difference of the resulting equations:

FIG. 8 .
FIG. 8.A comparison of the error bars of the β Ã ¼ 10 cm optics of the October 2016 MD.Top: histogram of the error bars of beam 1. Center: histogram of the error bars of beam 2. Bottom: mean of the size the error bars.The mean of the analytical N-BPM method is a factor 4 more accurate than the 3 BPM method.The third set of values shows the mean of the error bars of the simulations as shown in Fig. 6 but for the whole ring.

TABLE I .
Estimates for the LHC gradient and misalignment errors.MQ are focusing and defocusing quadrupoles, MS are sextupoles.The values of the magnet errors are derived from magnetic measurements of