Description of the luminosity evolution for the CERN LHC including dynamic aperture effects. Part I: the model

In recent years, modelling the evolution of beam losses in circular proton machines starting from the evolution of the dynamic aperture 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 this paper, the focus is on the derivation of scaling laws for luminosity, which include both burn off and additional pseudo-diffusive effects. It is worthwhile stressing that time-dependence of some beam parameters can be taken into account in the proposed framework. The proposed models are applied to the analysis of a subset of the data collected during the CERN Large Hadron Collider (LHC) Run~1 in a companion paper (Part II).


Introduction
Since the advent of the generation of superconducting colliders, the unavoidable non-linear magnetic field errors have plagued the dynamics of charged particles inducing new and potential harmful effects. This required the development of new approaches to perform more powerful analyses and to gain insight in the beam dynamics. Parenthetically, in a number of cases, rather than developing new tools, accelerator physics imported existing tools from other domains of science, such as celestial mechanics, dealing with similar problems. Among the several studies that flourished in the nineties, that on the concept of dynamic aperture (DA), i.e. the region in phase space where bounded motions occur, was particularly fruitful. The literature on this topic is very large, and most of it beyond the scope of the results presented in this paper. Nevertheless, it is worth mentioning the work done on the scaling law of the DA as a function of time [2,3] for the case of single-particle beam dynamics. Indeed, such a scaling law was later successfully extended to the case in which weak-strong beam-beam effects are added to the beam dynamics [4]. More importantly, such a scaling law was proposed to describe the time evolution of beam losses in a circular particle accelerator under the influence of non-linear effects [5], and the proposed model was verified experimentally, using data from CERN accelerators and the Tevatron. Note that such a scaling law for beam intensity as a function of time is at the heart of a novel method to measure experimentally the DA in a circular ring [6].
The model developed represents a bridge between the concept of DA, which is rather abstract, and the beam losses observed in a particle accelerator. Clearly, the next step was to extend the model to describe the luminosity evolution in a circular collider. The first attempts are reported in [7,8]. However, in those papers the DA scaling law was used without disentangling the contribution of burn off. Although the results were rather encouraging, to recover the correct physical meaning of the model parameters it was necessary to include as many known effects as possible.
This limitation is removed in the model discussed in this paper. In fact, the proposed scaling law is combined with the well-know intensity decay from particles' burn off so that a coherent description of the physical process is provided.
Moreover, it is worth stressing that the proposed model can be generalised so to consider a time-dependence for some of the beam parameters describing the luminosity evolution, such as emittance.
The plan of the paper is the following: in section 2 the general model of luminosity evolution is presented, starting from the case with burn off and timeindependent beam parameters (section 2.1). In section 3 the proposed model is presented, which includes also pseudo-diffusive contributions from non-linear effects. This model is then shown in action in section 4, where these concepts are applied to the problem of the analysis of the distribution of the integrated luminosity whereas in section 5 the analysis of the optimal physics fill duration is dealt with. The implications of using different models, namely plain burn off or burn off together with pseudo-diffusive phenomena are also addressed.
Conclusions are drawn in section 6 while additional detail is presented in the Appendices. Additional refinements to the burn-off model aimed at including time-dependent emittance variations are discussed in Appendix A, while computations for the case of plain burn off in conjunction with explicit time dependence of some of the beam parameters are reported in Appendix B. The detail of the derivation of the proposed model is worked out in Appendix C.

General considerations
The starting point is the expression of luminosity, which is a key figure-ofmerit for colliders and, neglecting the hourglass effect, reads 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 and the fact that not all bunches are colliding in the high-luminosity experimental points is taken into account by introducing a scale factor.
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 the half crossing angle θ c and the transverse and longitudinal RMS dimensions σ * , σ z , respectively according to: 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: in which the dependence on the total intensity of the colliding beams is highlighted and the other quantities are included in the term Ξ .
Under normal conditions, i.e. excluding any levelling gymnastics or dynamicbeta effects, only the emittances and the bunch intensities can change over time. Therefore, Eq. (1) is more correctly interpreted as peak luminosity at the beginning of the fill, as in general L is a function of time. When the burn off is the only relevant mechanism for a time-variation of the beam parameters, it is possible to estimate the time evolution of the luminosity, which turns out to be derived from the following equation where σ int represents the cross section for the interaction of charged particles and the two colliding beams have been assumed to be of equal intensity. The value used is 73.5 mb for 3.5 TeV and 76 mb for 4 TeV [9,10] for protons, representing the total inelastic cross-section. Here, n c stands for the number of collision points. The solution is given by where N i stands for the initial intensity. Equation (5) implies that the luminosity evolves as In the most general case, where both beams can have different intensities, the intensity evolution is described by the following equations For reasons that will become clear later, a different time variable is considered, τ being an adimensional variable representing the number of turns, where a shift of the origin of τ with respect to t has been introduced. In the following, the derivative with respect to τ is indicated by˙, while indicates the derivative with respect to t.
The solution of Eq. (7), indicated as N bo 1,2 (τ ) to highlight that it only includes the burn off contribution, can be obtained by re-writing: and from which one finds Equation (11) has two solutions depending on the value of ξ. If ξ = 0 then where N i = N i,1 = N i,2 stands for the initial beam intensity.
Otherwise, if ξ = 0 then where ξ = N i,1 − N i,2 and N r = Ni,2 Ni,1 . Note that the result does not depend on the sign of ξ, although, from a computational point of view a negative exponent is preferred, as τ can grow to very large values. It is possible to enforce this by This gives the same result as in Eq. (13) using N i,max and N i,min instead of N i,1 and N i,2 . Note also that Eq. (12) can be easily recovered from Eq. (13) by expanding the exponential and carefully setting As an example, in Fig. 1 the behaviour of the intensity and luminosity evolution is shown for various cases with equal or different intensities for the two beams. The luminosity evolution is normalised to the value corresponding to the solution with equal intensities. The different behaviour is clearly visible. The case with equal intensities turns out to be the best scenario as far as the luminosity is concerned.
To compare the results, the constraint N i,1 N i,2 = N 2 has been imposed.
The differences are clearly visible. Firstly, the time behaviour is exponentiallike for the case of different intensities, thus generating shorter fills. Secondly, in case of unequal intensity one beam is completely burnt while the other is not, which makes the use of the protons in the collider inefficient. Furthermore, the evolution of the relative luminosity, i.e. the luminosity normalised to the corresponding value of the case with equal beam intensities, shows that the most favourable situation is the one featuring equal beam intensities, as any imbalance in beam intensity worsen the collider performance. This is due to the fact that the burn off is the same for both beams, hence the maximum fill length is given by the time needed to burn the weaker beam.
Whenever additional time dependence in the luminosity evolution needs to be taken into account, and assuming that a model for such a time-dependence is available [12,13], then the solutions (12) and (13) can be extended to take into account these effects. The detail is presented in Appendix A and Appendix B.

Luminosity evolution including pseudo-diffusive effects
An efficient modelling of the luminosity evolution in a real collider can be obtained either by means of numerical tracking or by means of analytical or semianalytical models, grasping the main features of the luminosity evolution. The first approach has been used, e.g. in Ref. [11], where the luminosity evolution is studied by means of direct numerical simulations taking into account all relevant physical processes. The second approach, has been used in Refs. [12,13], where phenomenological fit models were proposed and applied with success to the characterisation of luminosity evolution in the Tevatron machine. The functional form of the proposed models was suggested by considerations on scaling laws of key quantities such as emittances. Another example is given in Ref. [14], where a refined semi-analytical model including several phenomena affecting emittance variation has been derived in view of optimising the performance of future circular colliders. However, none of the models studied included the effect of non-linear motion and this is at the heart of the approach proposed in Ref. [7]. The basis for such a model is the evolution of the dynamic aperture (DA) with time in a hadron collider. The analysis of single-particle tracking results showed that the time evolution of the DA follows a simple law [2,3], whose justification is not entirely phenomenological. Recently, this approach has been successfully applied to the analysis of intensity evolution in hadron machines [5]. So far, however, the results were obtained in the case of singleparticle simulations or whenever the conditions in a particle accelerator were not under the influence of any collective effect. To extend the proposed scaling law to luminosity evolution, it is necessary to show that it is valid also in the presence of beam-beam effects. This is the case at least for the results of numerical simulations in the weak-strong regime, as discussed in Ref. [4], thus opening the possibility to justify the proposed interpretation.
The proposed approach is a refinement of what is presented in Ref. [7] and assumes that all possible pseudo-diffusive effects can be modelled by a scaling of the intensity with time as where The parameters D ∞ , b, κ are normally fitted to the experimental data and the variable τ represents the turn number and satisfies τ ∈ [1, +∞[. It is worthwhile stressing some properties of the parameters as highlighted in Refs. [3,5], where two regimes were identified: • in 4D systems the three quantities D ∞ , b, κ are all positive [2,3]. This corresponds to having a stable region in phase space for an arbitrarily long time and hence being in the condition of applicability of the Kolmogorov-Arnold-Moser (KAM) and Nekhoroshev theorems.
• in 4D systems with tune modulation or 6D models, it is possible that there is no stable region even for a finite number of turns [3]. This corresponds to one of the two cases: The first case represents a situation with global chaoticity [3]. In particular, the fact that tracking data could be fitted with negative κ has been already observed several years ago (see discussion in Ref. [3]). The latter is compatible with a situation in which the stable KAM area shrinks to zero and the escape to infinity is governed by a Nekhoroshev-like behaviour.
Clearly • The fit parameters in Eq. (15) might depend on the beam intensity. However, if one assumes that the overall intensity variation over one physics fill is not too large, it is then reasonable to consider that the pseudo-diffusive effects are, to a good extent, almost constant.
Then, under these assumptions, it is justified to describe the intensity evolution The termsD i represent the intensity-independent pseudo-diffusive effects. The reason for a shift of origin of τ with respect to t is now clear: it is needed so that t = 0 corresponds to τ = 1 and then, D(1) = +∞ according to Eq. (15), which is the expected situation.
Equation (7) becomes and that the explicit solution has been assumed to be of the form (14) [4,5,7].
Therefore, one obtains The detail of the derivation can be found in Appendix C, but under the assumptions that the initial beam intensities are the same as well as the terms D j , then an explicit expression at the lowest order in ε (see Eq. (10)) can be given for both intensity and luminosity, namely and where L i is the initial value of the luminosity, given by L i = Ξ N 2 i . In case ε = ε(τ ) is time-dependent, an analogous derivation can be applied as shown in Appendix B and Appendix C. The result is very similar to the time-independent case, namely and where µ(τ ) is defined in Eq. (B.5). Note that the time-independent and the timedependent solutions are related by two substitutions, namely (τ − 1) ε → µ(τ ) and D 2 (τ ) → D 2 (τ ) (τ )/ (1). The second substitution deserves some comments. The scaling law (14) implies that D(τ ) is measured in terms of beam sigma, under the assumption that sigma does not change with time. In case of emittance growth, this assumption is violated and the proposed substitution simply takes into account the time-dependence of sigma (which scales with respect to the emittance as σ ∼ √ ). Note that this approach is possible because it assumes a complete independence of the pseudo-diffusive and emittance growth effects. This is fully justified by the assumption that the pseudo-diffusive effects are acting on the tails, while emittance growth affects the beam core. Note that these replacements have a general validity and will be true also for what follows.

Integrated luminosity over a physics fill
The models analysed in the previous sections can be used to derive some useful scaling laws for the integrated luminosity as a function of the length of the physics fill. Indeed, assuming the simple case of equal intensities for both beams, it is possible to obtain for the burn off part which is obtained by means of Eq. (6) and the very definition of ε. Note the so that one can normalise the integrated luminosity aŝ It is straightforward to express the integrated luminosity as a function of τ , starting from Eqs. (3) and (12), and integrating over [1, τ ], namely .
HenceL bo int (t) and L bo int (τ ) are related by a scale factor 1/f rev : Note that because L bo norm (t) and L bo norm (τ ) are now equal: Furthermore, by using the normalised turn variablē L bo norm can be recast in the following form Hence, L bo norm (τ ) has a very simple scaling law in terms ofτ . This allows comparing experimental data from physics runs with different beam parameters, such as β * , crossing angle, bunch intensity, and number of bunches (see [1]).
This treatment can be generalised to the case when the initial intensities are not the same, i.e. using the notation of the previous section, when ξ = 0. In this case, the final result reads In this case the scaling law then depends on two parameters, i.e. a re-scaled turn numberτ and N r .
As mentioned in the previous section, it is easy to extend this derivation to To include pseudo-diffusive effects it is enough to apply the computations made before to the general solution of the intensity-evolution equation, based on the sum of components N bo 1,2 (τ ) and N pd 1,2 (τ ), hence giving where L bo norm stands for the burn off component of the luminosity evolution derived in Eqs. (27) or (34). It is worth recalling that the terms N bo 1,2 and N pd 1,2 depend on the small parameter ε. It is therefore possible to find an expansion of L pd in terms of ε and to truncate it at an appropriate order, which, in our case, will be the first one: The generalisation to different intensities is straightforward: In case ε is time-dependent, it should be included in the integrand.

Optimal physics fill duration
A relevant quantity for evaluating the performance of a circular collider is the optimal fill length. A detailed discussion about the possible optimisation of this quantity has been presented in Ref. [14], but there, non-linear effects were not taken into account. Indeed, the models developed in the previous section allow providing estimates of this quantity including also the impact of non-linearities. In the following, the relation is considered, 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. According to the convention set in previous sections the time is always expressed in terms of turn number. With these assumptions n represents the number of physics fills in one year of operation and the total integrated luminosity is given by where only the burn off has been taken into account. The optimal fill length τ fill can be obtained by setting to zero the derivative of L bo tot,norm with respect to τ and it reads where the exponential terms in Eq. (34) have been approximated considering that the arguments are small and the power series has been truncated to first order in ε. Clearly, the expression for ξ = 0 in Eq. (41) tends to that for ξ = 0. Equations (41) provide the scaling laws of the optimal fill length when the burn off is the mechanism for intensity reduction. Indeed, τ bo fill scales with intensity and the turnaround time, but only as the square root. As an example, in Fig. 2 the optimal fill length is shown as a function of τ ta for different configurations with equal and unequal beam intensities. Also in this case N i,1 N i,2 = N 2 has been set, to allow for a direct comparison of the various configurations of the constraint.
It is also possible to estimate how L bo tot,norm (τ bo fill ) scales with respect to, e.g. τ ta , thus obtaining The increase in L bo tot,norm generated by a reduction of τ ta scales only as the square root of τ ta , whereas L bo tot,norm decreases with increasing τ ta as the inverse of the turnaround time. All this indicates how easily the performance of a collider can be spoilt. The value of L bo tot,norm (τ bo fill ) depends on N i , which represents either the common beam intensity or the highest one. The next step is to evaluate how the term L pd changes the conclusions concerning the optimal fill time. In this case the total normalised luminosity reads where L bo tot,norm (τ ) is the term for which the optimal fill time is τ bo fill . The value τ fill that maximises L tot,norm (τ ) is obtained by setting to zero the derivative of L tot,norm (τ ). An approximate expression can be obtained by developing the first-order derivatives of L bo tot,norm (τ ) and L pd tot,norm (τ ) taking into account that, by definition,L bo tot,norm (τ bo fill ) = 0. . (45)

Conclusions
In this paper, a model for luminosity evolution based on the scaling law of dynamic aperture as a function of time, in addition to particle burn off, has been presented and discussed in detail.
Estimates of the luminosity evolution have been worked out, together with prediction for the scaling law of the integrated luminosity as a function of the fill length and of the optimal fill length. While most physical parameters have been considered as being time-independent, generalisation to cases in which these parameters vary with time is provided in the appendices. Non-negligible differences with respect to simpler models based on burn off only have been clearly found.
Application of the proposed approach to collider's data is discussed in a companion paper [1], where a subset of the LHC Run 1 data is analysed and discussed in detail.

Appendix A. Refinement of the burn off model
The model that includes only the burn off can be improved to take into account also intensity-independent emittance variation, i.e. either growth or damping. These phenomena include beam-residual gas scattering and radiation damping. Here, the following form for the variation of the normalised emit- where the term −λ 1 * describes a damping of the beam emittance, and λ 2 a growth. The emittance evolution is given by the following expression * (τ ) = where * 1 = * (1). The combined effect of damping and growth makes it possible to define an equilibrium emittance * eq given by * eq = The relative variation of the normalised emittance is readily computed to be * (τ ) − * Under the assumptions mentioned at the beginning of this section, it is 1 It is worth stressing that throughout the paper the concept of emittance always refers to the normalised emittance. This has an effect on the values of the standard expressions for the coefficients λ i .
possible to give the expressions of λ i as [25,26] where r p , m p c 2 , E, ρ x are the classical proton radius, the proton mass, energy, and bending radius, respectively. Typical values for the LHC at 4 TeV are λ 1 ≈ 7.52 × 10 −10 . Furthermore, the components of λ 2 can be expressed as [25,26]  with where β x , α x , γ x are the Twiss parameters [27] and D x , D x the horizontal dispersion and its derivative, respectively. For typical LHC parameters corresponding to the runs at lower energies in 2011 and 2012, λ 2,rad ≈ 1.38 × 10 −23 m, hence completely negligible.
The second component of λ 2 is expressed by where N a , R, T, P stand for the Avogadro number, the gas constant, the temperature in Kelvin, and the pressure in Pascal, respectively. The properties of the residual gas are taken into account via A and Z, the atomic and the charge number, respectively. Finally, the impact of the ring's optics is taken into account by the term β representing the average beta-function, while p is beam momentum.
In Table A.1 the parameters corresponding to the various gas species in the LHC for nominal conditions [28] are listed, together with the corresponding values of λ 2,scat . It is now possible to verify that the condition * eq / * 1 << 1, required to neglect the contribution of λ 2 with respect to that of λ 1 in the emittance evolution, is always satisfied at the LHC with nominal parameters from Table A.1. Moreover, it is also easy to check that even the impact of λ 1 , assuming Run 1 beam parameters, can be correctly neglected. Nonetheless, for the sake of completeness, the complete discussion of the considered emittance growth model is reported in Appendix B.

Appendix B. Evolution of beam intensity and luminosity with burn off and additional time dependency
Let us consider the homogeneous part of Eq. (18) assuming that the luminosity is changing not only due to particles' burn off, but also for other timedependent effects, such as emittance growth. In this case it is possible to recast the equation in the form Once more, the solutions are of two types, as they can satisfy 2) can be solved using standard methods, giving Note that this solution has the same structure as the time-independent case in Eq. (12) after making the substitution 3) can also be solved with standard methods and one obtains which again holds for ξ = 0 independently of the sign of ξ.
where N r has been defined before, i.e. N r = N i,min /N i,max . Once more, the possibility of evaluating in closed form the integral of ε(τ ) is the condition to find a closed form expression for µ(τ ) and hence for the integrated luminosity.
Normalisation factors for the integrated luminosity can be found as so that a very simple expression for L bo norm (τ ) can be found, namely It is readily seen that the expressions for the normalised integrated luminosity feature rather simple scaling laws with respect to µ(τ ), very much like the timeindependent case.
The possibility of finding an expression for the optimal physics fill length can also be explored. To this aim the total luminosity is given by where the symbols have the same meaning as in the main body of the paper.
The optimal fill length can be computed by looking for the solutions of the equation: The derivative of the burn off contribution to the integrated luminosity can be computed to be (B.14) Therefore, the optimal fill time τ bo fill (taking only the burn off into account) is the solution of the following equations According to Eqs. (B.4) and (B.7) the intensity evolution reads for the case of equal intensity and whenever the initial intensities are different. It is worth noting that the timedependence of the beam intensity becomes much more involved when emittance evolution is included in the model.
The normalisation derived in this Appendix can be applied and by using Eq. (B.11) one obtains can be done only numerically. As a last remark, it is worth stressing that the split of the luminosity evolution into the sum of two parts according to Eq. (36) remains unchanged.

Appendix C. Derivation of the proposed model
Equations (18) can be cast in the form of a Riccati equation [15], although this does not provide any useful information about the properties of its solutions.
Therefore, a different approach has been applied, i.e. constructing the solution as the sum of two components, namely which, applied to Eq. (18) gives By replacing the developments in power series of N bo j (τ ), N pd j (τ )  into Eq. (C.2) and using the Cauchy product of series, one obtains and The solution of the system (C.5) is trivial and is given by (C.7) Equation (C.6) can be solved by re-casting it into form (9) thus giving  The functions N pd j (τ ) depend on the parameters ξ m , which can generate a difference between N pd 1 (τ ) and N pd 2 (τ ) for τ = 1, given that N pd 1 (1) = 0. As the N pd j (τ ) represent the components related with the pseudo-diffusive effects, it is reasonable to set ξ m = 0, m ≥ 1 as any intensity imbalance between the two beams is taken care of by the functions N bo j (τ ) and N pd 0,j (τ ). If additional time dependencies are present, i.e. ε = ε(τ ), extra care should be taken as a plain power expansion of N bo j (τ ), N pd j (τ ) in ε is not appropriate. Instead, we can write ε as: where ε 1 = ε(1) its initial value. Note that ε 1 1 ,ε(τ ) ∼ 1 , (C.11) implying we can expand N bo j (τ ), N pd j (τ ) in function of the small (and constant) parameter ε 1 : N pd m,j (τ ) ε m 1 , j = 1, 2 . (C.12) As before, we replace these developments into Eq. (C.2) (now with ε → ε(τ )).
This leaves Eq. (C.5) unchanged, so the O(ε 0 1 ) solution is still valid as given by Eq. (C.7). For the higher order terms we can use the same steps as before, in order to finally get the full solution with time-dependent ε(τ ):