Description of the luminosity evolution for the CERN LHC including dynamic aperture effects. Part II: application to Run 1 data

In recent years, modelling the evolution of beam losses in circular proton machines starting from the concept of dynamic aperture its time evolution has been the focus of intense research. Results from single-particle non-linear beam dynamics have been used to build simple models that proved to be in good agreement with beam measurements. These results have been generalised, thus opening the possibility to describe also the luminosity evolution in a circular hadron collider. In a companion paper (Part I), the derivation of a scaling law for luminosity, which includes both burn off and pseudo-diffusive effects, has been carried out. In this paper, the proposed models are applied to the analysis of the data collected during the CERN Large Hadron Collider (LHC) Run 1. A data set referring to the proton physics runs for the years 2011 and 2012 has been analysed and the results are proposed and discussed in detail in this paper.


Introduction
The unavoidable non-linear magnetic field errors that plague the dynamics of charged particles of modern superconducting colliders, inducing new and potential harmful effects, require the development of new approaches to perform $ Research supported by the HL-LHC project more powerful analyses and to gain insight in the beam dynamics. It is the case of the critical revision of the concept of dynamic aperture (DA) and of its dependence on time [2,3] for the case of single-particle effects. Such a scaling law was later successfully extended to the case in which weak-strong beam-beam effects are taken into account [4]. More importantly, this scaling law paved the way to describe the time evolution of beam losses in a circular particle accelerator under the influence of non-linear effects [5], verified experimentally using data from CERN accelerators and the Tevatron, which is at the heart of a novel method to measure experimentally the DA in a circular ring [6].
The description of the luminosity evolution in a circular collider profited from this novel framework. The first attempts to derive a new model are reported in [7,8], whereas a more complete and accurate modelling is described in [1].
The proposed approach is put in action by analysing a selection of luminosity data collected at the LHC during Run 1, in particular for the proton physics runs in the years 2011 and 2012, which is the focus on this paper, where the result of these analyses is discussed in detail. Furthermore, it is worth mentioning that the scaling law [5] has been used also in the analysis of dedicated beam-beam experiments performed at the LHC [9,10].
Two points are worth stressing. Firstly, the focus of this paper is a test of the descriptiveness of the novel model, which is probed by checking the quality of the agreement with the LHC data. Note that the issue of predictiveness will be addressed in a different paper. Secondly, the arguments used to build the proposed model are rather general, thus implying that they should not be applicable to LHC only, but to circular colliders in general.
It is recalled that the starting point is the expression of luminosity, which is a key figure-of-merit for colliders that, neglecting the hourglass effect, reads [11] L = γ r f rev k b n 1 n 2 4 π * β * F (θ c , σ z , σ * ), where γ r is the relativistic γ-factor, f rev the revolution frequency, k b the number of colliding bunches, n i the number of particles per bunch in each colliding beam, * is the RMS normalised transverse emittance, and β * is the value of the beta-function at the collision point. The total beam population is defined as N j = k b n j . Different bunches have different collision schedules, meaning that they collide in different interaction points. Hence, the following analysis could have been performed on a bunch-by-bunch basis and considering bunches with the same collision schedule as members of the same class. However, a simplified approach has been applied, namely the total intensity N j has been rescaled by , where k b,ATLAS,CMS represents the number of bunches colliding in the two high-luminosity experiments. The underlining assumption is that the effects generated by the collisions in the two other low-luminosity experiments can be neglected. This is supported by the difference in typical peak luminosities, whose ratio to those of the high-luminosity experiments is The factor F accounts for the reduction in volume overlap between the colliding bunches due to the presence of a crossing angle and is a function of half the crossing angle θ c and the transverse and longitudinal RMS dimensions σ * , σ z , respectively according to [11]: Note that σ * = β * * /(β r γ r ), where β r is the relativistic β-factor. Equation (1) is valid in the case of round beams ( * x = * y = * ) and round optics (β * x = β * y = β * ). For our scope, Eq. (1) will be recast in the following form [11]: in which the dependence on the total intensity of the colliding beams is highlighted and the other quantities are included in the term Ξ .
The plan of the paper is the following: in section 2 a global overview of the LHC data is presented, including a discussion of the time-dependence of some beam parameters, and the choice made for the data analysis. The application of the proposed model is presented in three different sections each dealing with one of the three main observables considered, namely the time evolution of the luminosity over a physics fill (section 3), the integrated luminosity over a physics fill (section 4), and the optimal duration of a physics fill (section 5). Conclusions are drawn in section 6, whereas the detailed discussion of the important features of the numerical models proposed in this paper is presented in Appendix A.
2. LHC data from Run 1

General considerations
The models derived in the companion paper [1] will be applied to the analysis of the LHC performance data collected during Run 1. Detailed information on this topic can be found in Refs. [15][16][17][18], while in Ref. [19] a preliminary analysis was made, without focusing on models to describe the luminosity and its time evolution. Here, the focus will be on the proton physics run and the data analysed can be found at [20]. As an example, the evolution of some key parameters is shown in Fig. 1 as a function of the fill number, which is an incremental integer number representing in a unique way the physics fill. The peak luminosity L and the total beam intensities are shown in the upper row of Fig. 1, while * 1 is reported in the middle row. The latter can be derived from the knowledge of L and the beam parameters entering in Eqs. (1) and (2) according to * with the assumption that the emittances of the two beams are equal as well as the value of β * . In the lower row the evolution of β * and k b is also shown. Data for the 2011 and 2012 runs are shown in the left and right columns, respectively.
The total beam intensity has been increasing throughout the 2011 run, levelling out in 2012. Correspondingly, the peak luminosity has increased also because of the reduction of * 1 and β * and the increase of the number of bunches k b . A step decrease in * 1 is clearly seen in the 2011 data and corresponds to the progressive reduction of the controlled transverse emittance blow-up applied at the SPS [21]. Since that change, the emittance is fairly constant even during the 2012 run. It is worth mentioning that the beam brightness is defined by the LHC injectors' complex and is basically constant, which implies a linear relationship between intensity and transverse emittance. The data shown in Fig. 1 are also used in the following analysis of the luminosity evolution. Among the full data set available from [20] a selection has been considered including only the fills that resulted in successful physics runs, the so-called stable beams, of a total duration exceeding 10 3 s and featuring N i,1,2 > 10 13 p. Such a filtering allows removing data corresponding to beam commissioning stages or low luminosity fills, which would not be representative of the typical LHC performance. Additionally we only select those fills that have a number of bunches k b > 1300 to exclude ion runs.
The analysis of this data set showed that the difference in beam intensity at the beginning of a physics fill is rather small, at the level of few percent [19].
Hence, the simplifying assumption N i,1 = N i,2 is fully justified and is used in the following. Using L i = Ξ N i,1 N i,2 we can calculate ε from the initial luminosity where n c = 2 because the vast majority of protons are burnt in the two high-luminosity interaction points (see the previous comment on the relative luminosities of the various LHC experiments), and the total inelastic cross-section for proton-proton collisions is σ int is 73.5 mb for 3.5 TeV and 76 mb for 4 TeV [13,14] for protons, representing the total inelastic cross-section.

Observed time-dependence of beam parameters and other assumptions for data analysis
The analyses presented in Ref. [1] included the situation when some beam parameters are changing during the fill, as it can be the case for the rms bunch length σ z or any of the two transverse normalised emittances * x,y . Equations (1) and (2) show that while σ z has an impact on F , only, the transverse normalised emittances * x,y affect both F and the peak luminosity. The measured data revealed that the variation of σ z over a typical physics fill does not exceed ≈ 7 %. Such a small variation is understandable by considering that at the collision energy of 3.5 TeV or 4 TeV, which are the values for the 2011 and 2012 runs, respectively, the damping generated by synchrotron radiation is not too strong. Therefore, the time-dependence of σ z can be safely neglected in the analyses presented in the following sections.
The time-dependence of ε needs to be assessed to decide the approach to be applied to the data analyses and the outcome of these investigations, based on the selected fills of the 2011 and 2012 runs, is shown in Fig. 2.
The data have been fitted using an exponential function and the result is given by ∆ ε(t) = 34.69 e −0.1358 t − 35.39 (6) where t is expressed in hours and ∆ε in percent. The fit quality is given by (8) in section 3 for the definition of this quantity). For the majority of fills ∆ ε does not exceed ≈ 30 % and it has been decided to perform the numerical analyses assuming a time-independent ε. Note that the fit models presented in section 3 have been cross-checked against models in which ε had been assumed to be time-dependent and the differences have been found small, thus confirming that the assumption made is appropriate. It is worth mentioning that it is planned to apply the most general model described in Ref. [1] to the description of data for the LHC Run 2 [22]. Another aspect to consider is whether some of the parameters entering in the proposed models should be different for the two beams. In fact, in Ref. [1] the descriptive models could include beam-specific values for both the initial beam intensity and the pseudo-diffusive effects. A close inspection of the Run 1 data [19] reveals that for a typical physics fill the quantity 2 |N i,1 − N i,2 |/(N i,1 + N i,2 ) does not exceed ≈ 10 %. Hence, in the analysis reported in the following sections, the two initial beam intensities have always been assumed equal and the corresponding model for the burn-off part has been used [1]. Given that a similar estimate holds also for the intensities at the end of the physics fills, the pseudo-diffusive effects have been assumed to be the same for both beams.

Time evolution of luminosity over a fill
The first step in the analysis of the LHC Run 1 data is the fit of the pseudodiffusive component of the luminosity evolution based on the expression, which was derived in [1], given by where L i is the initial value of the luminosity, given by L i = Ξ N 2 i , and N i the corresponding initial value of the beam intensity. . Also shown is R 2 adj , the so-called adjusted coefficient of determination given by where N is the sample size, ν the number of fit parameters, Σ 2 the sum of residues squared defined in Eq. (A.1), and is the total variance of the data withȳ the average over all y i . Note that R 2 adj compares the fit under consideration to the simplest fit, i.e. a constant line through the mean. When R 2 adj 1 (or possibly even negative), the fit is of poor quality as the mean of the data provides a better fit than the proposed model. A good fit has R 2 adj → 1, indicating that the residues are small compared to the data variance.
If we look at the results in Fig. 3, we notice that all fits are of particular good quality, as all except one have R 2 adj > 90 %, while for all fits from 2011 this is even R 2 adj > 99 %. There is a clear distinction between the results for 2011 and those for 2012, both in spread, but also in behaviour, as the yearly average value of D ∞ is negative for 2011 whereas it is positive for 2012. The larger spread in the fit parameters for the 2012 run might be generated by the transverse instabilities that plagued that run [23]. It is interesting to stress that from Fig. 3 no clear systematic difference between the fitted models based on the ATLAS or CMS data is found as the variations change on a fill-by-fill basis. The differences observed are not always within the error bars, which might be underestimated. Moreover, some differences are to be expected as β * waist position correction was not fully mastered, the first being at the level of 5 % [24]. A further refinement has been applied during Run 2 lowering the uncertainty on β * at the level of about 1 % [25].
Following the discussion in Appendix A, it is useful to fit the data to a slightly adapted model, which has a reduced set of parameters. To this end, we selected three different configurations: one where we fix κ = 2 (according to the Nekhoroshev estimate [29]) and fit b and D ∞ ; one where we fix D ∞ = 0 (as it is approximately the average of Run 1) and fit b and κ; and one where we fix both κ = 2 and D ∞ = 0 and fit only b, thus leaving only one model parameter. The resulting weighted average parameter values are listed in Table 1

Integrated luminosity over a physics fill
The second step consists of establishing the model for the integrated luminosity delivered in a single fill for physics. In Ref. [1] it has been shown that, under the assumption of considering burn-off phenomena only, the integrated luminosity over a fill can be expressed as L bo norm (τ ) =τ 1 +τ (10) with an appropriate rescaling and by using a normalised time given byτ = ε N i (τ − 1). Whenever pseudo-diffusive effects are taken into account, then one can assume the following form for the luminosity evolution where L bo norm stands for the burn off component of the luminosity evolution [1] and L pd can be expressed at first order in the small parameter ε as where N i is the initial intensity.
The plot of L int (the integrated luminosity over a physics fill) is shown in The large spread observed for the 2011 data is due to the change of beam and ring parameters, i.e. transverse emittance, intensity, and β * occurred during the year 2011 (see Fig. 1), whereas the situation in terms of beam parameters has been much more stable throughout the 2012 run. The impact of the proposed normalisation of both L int and τ according to Eq. (10) is also shown in Fig. 4 (lower). The spread of the data points is almost completely removed and a sort of universal curve is appearing, with a similar shape for both 2011 and 2012 data. The normalised luminosity allows for an easy recognition of outlying data points, which have been removed (in total 12 data points) when fitting the data to the evolution model.
The luminosity data are given as a function of time instead of number of turns, hence the least computationally expensive way to obtain the pseudodiffusive contribution to luminosity from these data is to rearrange Eq. (11) into and use Eq. (5) to evaluate ε for every data point. When the data, calculated from Eq. (13), has been fitted to the model (12), it is recast into the normalised luminosity using Eq. (11).
As a first investigation, the pseudo-diffusion model has been fitted to the complete Run 1 dataset. This is shown on the left side of Fig. 5, and the values of the fit parameters including the associated errors are reported in Table 2.
From Fig. 5 it is clear that including the pseudo-diffusive effects is very efficient to recover a nice agreement between model predictions and measured data, as the discrepancy between the burn off-only model and the measurements is rather large ifτ > 0.1. Furthermore, it is important to note that the fit parameters of the pseudo-diffusive component correspond, given that all parameters are positive, to a situation in which a stable region exists in phase space for an arbitrarily long time. Note also that κ is extremely close to 2, the estimate given in the proof of the Nekhoroshev theorem.
The pseudo-diffusive effect on a yearly basis is shown in the right plot of  A comparison of the fit parameters for the corresponding cases reported in Tables 1 and 2 shows that the values are compatible, within the errors, for the case of three-parameter fit, while the compatibility degrades as the number of fit parameters is reduced, the case with fixed κ being more compatible between the non-integrated and integrated luminosity models, than that with D ∞ . This confirms once more that fixing κ is the best option among those with reduced fit parameters.

Optimal physics fill duration
As discussed in Ref. [1], the optimisation of the performance of a circular collider can be performed by maximising the yearly integrated luminosity given by L bo tot,norm (τ ) = where τ ta is the turnaround time, i.e. the time between the end of a physics fill and the beginning of the next one, τ is the fill length that should be optimised, and T is the total time for physics over one year and only the burn off has been taken into account is considered. The optimal fill length τ fill can be obtained by setting to zero the derivative of L bo tot,norm . Of course, the term L pd changes the conclusions concerning the optimal fill time and an approximate expression reads [1] τ fill ≈ τ bo fill −L pd tot,norm (τ bo fill ) L bo tot,norm (τ bo fill ) +L pd tot,norm (τ bo fill ) (15) and it is possible to use Eq. (15) to estimate the optimal physics fill duration for our models, as shown in Fig. 7. Note that there is a clear difference between the case where only burn off is taken into account, and that where pseudodiffusive effects are included. The relation between the various models used to derive τ fill is similar to that for the fit of L pd , namely, for 2011 all four cases behave the same, while for 2012 only the case where two parameters are fixed (κ = 2 and D ∞ = 0) is different from the others.
Then, we carried out the comparison of the estimate for the optimal fill length with the actual fill duration during Run 1, focusing on the situation for the year 2012. The distribution of the actual values of τ fill as a function of τ ta is shown in Fig. 8. The data considered includes all fills for high-luminosity physics, excluding all other cases, e.g. special runs, commissioning periods.
Note, also, that τ ta for a given physics fill is computed as the time between the end of the previous physics fill and the beginning of that under consideration.
All these data have been divided into two groups: a class in which the end of the fill for physics is controlled by operation, so-called programmed dump, and a class in which the end of the fill for physics is triggered by the machine protection system, so-called protection dump. The relevance of this classification is that the first class allows for performing an optimisation of the fill length, while the latter does not. It is also worth mentioning that even in the case of the first class, there are some situations in which a beam dump is indeed triggered by operation, but in view of preventing a protection dump. It is, e.g. the case when the cryogenic conditions are going to be lost in a short while and the operator dumps the beams before a genuine protection dump occurs. This special subset of the first class explains the cases of programmed dumps with rather short τ fill . Figure 8 shows a number of interesting features: τ ta is always larger than ∼ 2 h, the minimum turnaround time based on the performance of the LHC hardware; τ ta can reach rather high values, which indicate a fault that occurred in between fills for physics; for the class of programmed dump the length of the physics fill is clearly different from the corresponding optimal fill length, either because the physics fill is too short or because it is too long; for the class of protection dump some fill are also too long with respect to the optimal duration. While the proposed approach is very useful to qualify whether a single fill is optimal with respect to its duration, the overall optimisation of the yearly performance, in terms of integrated luminosity, of a collider is a much more complex task.

Conclusions
The models proposed in the companion paper [1] have been benchmarked against the data from the LHC Run 1, with special emphasis on the years 2011 and 2012. The inaccuracy to reproduce the LHC data using only burn off has been confirmed by the analysis made, while the proposed models showed a remarkable power in reproducing and describing the observed behaviours of luminosity as a function of time and of integrated luminosity.
A detailed discussion of the potential numerical issues related with the proposed fitting models is presented in Appendix A. It is also shown that the difficulties can be efficiently resolved by reducing the number of fit parameters, by using the estimates provided by the proof of Nekhoroshev theorem for some of them.
Given the encouraging results of the analyses reported in this paper, the data from Run 2 will be considered next. In fact, the higher beam energy that characterises the proton physics in Run 2 opens a new domain in terms of beam behaviour, such as strong longitudinal emittance damping due to synchrotron radiation as well as a burn-off dominated regime.
Finally, it is worth stressing that while in this paper the approach has been to fit the model parameters to measured data, in future the analysis can be shifted to using numerical simulations to provide the input about the dynamic aperture evolution, which is needed in the proposed models, to verify the agreement with observations from the LHC. In this way the descriptiveness of the proposed approach might turn into predictiveness, which could be used to assess the performance for future colliders and in particular the luminosity upgrade of the LHC.      This essential observation leads to the conclusion that decreasing the number of fit parameters might be advisable. A natural option is to fix κ to the value provided by the estimate in the proof of Nekhoroshev theorem [29], i.e. κ = (d + 1)/2, where d is the number of degrees of freedom of the system under consideration. Other possibilities will be considered and discussed in the next sections.
Another important aspect to consider is the estimate of the error associated with the fit parameters. Such an estimate is closely linked with that of the error associated to the luminosity measurement and crucially depends upon it.
Given that it is not so easy to provide a reliable estimate for the total error on luminosity, the so-called Bootstrap Method [27,28] has been applied. Here, the population of residues of the fit is used as a distribution to generate 10 000 realisations with replacements, that are again fitted to the model. The set of realisations of each parameter will then be distributed around the original results of the fit procedure. Therefore, the parameter's distribution can provide an estimate of the error associated with the model parameter by means of the bias-corrected accelerated (BCa) Bootstrap confidence interval [28] at 1σ. For the studies that exhibit degeneracy in the parameter space, it is often the case that the fitted parameter value lies outside the BCa interval, even at 90 %. In such case the BCa interval calculation is less meaningful and we estimate the error by the standard deviation of the distribution of the fit parameters over the realisations instead.