Impact of Non-Gaussian Beam Profiles in the Performance of Hadron Colliders

At the Large Hadron Collider (LHC), the interplay between a series of effects, including intrabeam scattering (IBS), synchrotron radiation, longitudinal beam manipulations, two beam effects (beam-beam, e-cloud) and machine non-linearities, can change the population of the core and tails and lead to non-Gaussian beam distributions, at different periods during the beam cycle. By employing generalised distribution functions, it can be demonstrated that the modified non-Gaussian beam profiles have an impact in the beam emittance evolution by itself and thereby to the collider luminosity. This paper focuses on the estimation of beam distribution modification and the resulting rms beam size due to the combined effect of IBS and synchrotron radiation by employing a Monte-Carlo simulation code which is able to track 3D generalised particle distributions (SIRE). The code is first benchmarked over classical semi-analytical IBS theories and then compared with measurements from the LHC at injection and collision energies, including projections for the HighLuminosity LHC (HL-LHC) high brightness regime. The impact of the distribution shape on the evolution of the bunch characteristics and machine performance is finally addressed.


I. INTRODUCTION
The performance of a high-energy hadron collider such as the LHC is heavily based on the preservation of the injected emittances, under the influence of several degrading mechanisms. In this respect, an emittance evolution model was constructed including the effects of intrabeam scattering (IBS), synchrotron radiation, elastic scattering and luminosity burn-off (while at collision) [1]. This simple model is based on semi-analytical approaches which assume Gaussian beam distributions, in particular for IBS. The bunch characteristics evolution predicted by this model revealed discrepancies, as compared to the measurements, translated to differences in the luminosity predicted by the model as compared to the experimental estimations [1][2][3]. One of the possible reasons for these differences could be attributed to the fact that the bunch profiles appear to be non-Gaussian both at injection and collision energies, i.e. 450 GeV and 6.5 TeV respectively. The aim of this study is to quantify the impact of the beam distribution shape on the emittance and luminosity evolution of hadron colliders. In order to illustrate this, and employing generalised distribution functions, the luminosity of non-Gaussian beams is determined in a closed form. The generalisation of the luminosity estimate for arbitrary distributions does not only permit its comparison to the usual Gaussian beam estimate but also the extension of classical results for the impact of non-Gaussian beam distributions to the luminosity (see Section II). This motivates the investigation of the emittance evolution beyond the classical analytical formulas for modelling IBS, which are based on 3D Gaussian beam assumptions [4]. In this respect, a Monte Carlo multi-particle simulation code for IBS and Radiation Effects (SIRE) [5,6], is employed and compared to LHC data from Run2. A brief description of the code, and its benchmarking with the existing IBS analytical approaches and simulations for lepton rings is presented in Section III. Extending previous benchmarking studies for the LHC with respect to IBS theories [7], a detailed comparison of the Bjorken-Mtingwa (B-M) IBS theoretical model with the SIRE code for both injection and collision energies is presented for the nominal LHC using the Batch Compression Merging and Splitting (BCMS) [8,9] beam and the high luminosity LHC (HL-LHC) [10] beam parameters (Section IV). Finally, measured data from the LHC corresponding to non-Gaussian longitudinal beam profiles are compared with the expectations of the rms emittance evolution given by the theoretical B-M analytical formalism [4] and the SIRE code (Section V).

II. MOTIVATION-IMPACT OF NON-GAUSSIAN DISTRIBUTIONS ON LUMINOSITY
The performance of a collider is determined by the luminosity which, for two beams colliding head-on, is given by [11]: with N 1,2 representing the number of particles per bunch for each beam, N b the total number of colliding bunches, f rev the revolution frequency and ρ the beam density distribution functions for each plane and beam. Based on Eq. (1), assuming Gaussian beams that collide head-on, the luminosity is expressed as [11]: In order to achieve high luminosity, high intensity bunches and small beam sizes are required. The horizontal and vertical beam sizes of two colliding Gaussian bunches are given by: where (σ 1x , σ 1y ) and (σ 2x , σ 2y ) are the transverse rms beam sizes of beam 1 and beam 2, respectively. Based on the transverse and longitudinal bunch profile measurements, it has been observed that the particle distributions in the LHC, both at collision and injection energies, appear to have shapes that differ from the ones of a normal distribution [12][13][14]. At the LHC injection energy, the emittance evolution is dominated by the IBS effect, both in the horizontal and in the longitudinal plane, while no effect is expected in the vertical plane [15] where dispersion is minor. From Run 2 data, it is observed that in many cases the transverse bunch profiles appear to be non-Gaussian during the whole injection plateau [12]. At the LHC collision energy, the IBS effect is weaker, while synchrotron radiation damping becomes more pronounced. The bunch profiles at collisions appear to have non-Gaussian tails, as well. In fact, during the energy ramp, the bunches that are blown up longitudinally in order to avoid instabilities due to the loss of Landau damping [16], arrive at the start of collisions with a clearly non-Gaussian shape [14]. By assuming that a particle distribution is Gaussian when this is not the case, not only the rms beam size may be underestimated or overestimated, but also its impact on performance parameters, such as the luminosity. Therefore, it is important to use appropriate fitting functions (or some type of interpolation algorithm) on the beam profile in order to properly address this discrepancy. A generalized Gaussian function, called the q-Gaussian [17], can be employed for fitting more accurately bunch profiles with shapes that differ from the ones of a normal distribution (see Appendix A for the properties of this distribution function). The parameter q describes the weight of the tails as compared to the core, ranging from light tailed ones for q < 1 (including the square distribution for q → −∞) and extending to a heavy tailed ones for q > 1, passing through the Gaussian distributions in the limit of q → 1. This distribution is actually a stationary solution of a generalised Fokker-Plank equation which can cover a full spectrum of statistical behaviors of dynamical systems, from sub to super-diffusion Levy-type processes [18]. In view of quantifying the impact of non-Gaussian distributions, the luminosity is estimated through Eq. (1) by using the specified probability density functions. Assuming that the two beams are identical and that they collide head-on, the luminosity for q-Gaussian distribution functions in the transverse plane is given by: for σ qG x,y being the rms beam sizes (see Appendix A) in the transverse plane, for both beams. The factors I qG x,y which depend on the parameter q in the respective planes and the details of the calculation are presented in Appendix B, together with the validation of the luminosity estimation for q-Gaussian distributions (shown in Fig. 22). By comparing this equation to the standard luminosity formula for Gaussian beams with identical rms sizes, the significance of the tail population contribution on luminosity can be established and parameterised through q. This is illustrated in Figure 1 FIG. 2: The q-Gaussian density distribution function for a light tailed (blue, q < 1), a Gaussian (black, q = 1) and a heavy tailed (red, q > 1) bunch profile, having identical q-Gaussian rms beam sizes.
where the luminosity variation normalised to the corresponding one for Gaussian beams (L G ) is parameterized with the parameter q of the q-Gaussian fitting function, characterising the weight of the tails in the transverse plane, for fixed q-Gaussian rms beam sizes in all planes, assuming head on collisions (i.e. no dependence of the luminosity on the longitudinal beam size, see Appendix B). The bunch profiles corresponding to a light tailed (q < 1), a Gaussian (q = 1) and a heavy tailed (q > 1) distribution, having identical beam sizes, are plotted in Fig. 2. As q (and β qG ) vary for fixed q-Gaussian rms beam sizes (based on Eq. (A5) in Appendix A), the luminosity varies as well with respect to the one estimated for purely Gaussian beams. Practically, if the tails of a distribution differ by 10% compared to the ones of a Gaussian distribution (i.e. q = 0.9 or q = 1.1), the luminosity value can be overestimated or underestimated by 5%. It is then clear that, for two beams colliding head-on, the shape of the transverse distributions has a significant impact to the estimated luminosity, in particular for the LHC experiments which target a ∼ 2 % accuracy in their estimates [19]. The impact of non-Gaussian distributions on the luminosity was firstly discussed by Hereward [20]. In FIG. 3: In order to estimate the divergence of the luminosity for non-Gaussian distributions from the one for Gaussian densities, the Lw < w 2 > 1/2 (see Eq. 5) is plotted versus the weight of the tails q, for light tailed (left) and the heavy tailed (right) distributions, i.e. q < 1 and q > 1, respectively. The points corresponding to a parabola and a Gaussian distribution, firstly presented in [20,21], are plotted. particular, the luminosity integrals were calculated for several distributions as the rectangular and the parabolic which correspond to a q-Gaussian with q → −∞ and q = 0, respectively. Assuming that the two beams are identical and that they collide head-on, these distributions were used as examples to estimate the discrepancy from the luminosity for Gaussian densities. This discrepancy is calculated after solving the L w that is the integral over the density functions in one plane, being identical for both beams, as: with w = x, y and for σ w being the transverse rms beam size, since these solutions correspond to the transverse plane.
It was found that [20,21]: L w < w 2 > 1/2 = 0.2887 , for a rectangular distribution 0.2683 , for a parabolic distribution 1 2 √ π = 0.2821 , for a Gaussian distribution (6) In fact, this approach already identified the existence of a minimum for a light tailed parabolic distribution, which becomes obvious by employing the q-Gaussian, as observed in Fig. 3, where L w < w 2 > 1/2 is plotted versus q for q < 1, i.e. light tails (left), and q > 1, i.e. heavy tails (right). The results for q-Gaussian distributions (grey curves) are in perfect agreement with the case studies discussed in [20]. This is also true for a rectangular distribution which corresponds to a q-Gaussian with q → −∞ and is beyond the range of the left plot of Fig. 3. For heavy tailed distributions, there is no upper limit for the constant of Eq. (6), as already inferred by Hereward [20]. In Fig. 3 (right) the case of a heavy tailed distribution with q = 1.65 is denoted by a red square. The extreme case of q → 5/3 corresponds to a q-Gaussian whose rms size goes to infinity (i.e. Levy distributions, see Appendix A).
The sensitivity of the luminosity on the distribution as generalised by employing the q-Gaussian function justifies the need of carefully studying the evolution of distributions in hadron colliders. For the LHC, a luminosity model was constructed [1]. The evolution of the emittance includes the effects of IBS, synchrotron radiation, elastic scattering, betatron coupling, noise and burn-off. Although this model has a relative agreement with respect to the measured luminosities by the experiments (i.e. ATLAS [22] and CMS [23]), there is still some room for improvement [3]. Indeed, the emittance evolution for IBS was based on the module of MAD-X [24] following the B-M theory which assumes Gaussian beam distributions. The extension of this to non-Gaussian distributions as observed in the LHC may shed light to the origin of the remaining discrepancy between the model and the measurements.

III. INTRABEAM SCATTERING THEORIES, OBSERVATIONS AND SIMULATIONS
One of the statistical processes causing a spreading of particles in phase space or a continuous increase of beam emittance is the small angle multiple Coulomb scattering, called Intrabeam scattering (IBS), which plays an important role in e + /e − damping rings, high intensity/low energy light sources [25] and high intensity hadron [26] and ion [27] circular machines. The IBS theory for accelerators was firstly introduced by Piwinski [28] and extended by Martini [29], establishing a formulation called the standard Piwinski method. Later, Bjorken and Mtingwa (BM) [4] used a different approach to describe the effect, taking into account the strong focusing effect. The Modified Piwinski method [30] that includes the strong focusing effect, was developed by Bane. Some approximations of these theories are the high energy one by Bane [30] and the completely integrated modified Piwinski [31]. A different approach developed by Lebedev for hadron beams is based on a Boltzmann type integro-differential equation and includes betatron coupling [32].
The analytical models that describe the IBS effect [4,33] assume Gaussian beam distributions. The stationary solution of the Fokker-Planck equation is a particle distribution that is Gaussian in the phase space. However, taking into account the effects of IBS, radiation damping and quantum excitation but also other diffusive mechanisms [34], there is no evidence that the distribution remains Gaussian. Therefore, it is important to develop analytical formulas and simulation tools that calculate the interplay between these effects for any distribution.
For the Relativistic Heavy Ion Collider (RHIC), the IBS growth rates were calculated and benchmarked with experimental data using the distribution function evolution (based on the Fokker-Planck equations), extending the usual approach of employing a conventional Gaussian-like distribution [35]. In this respect, IBS growth rates were calculated for a bi-gaussian distribution, which was interesting for studying the possibility of using electron cooling in RHIC [36]. Later, a model which is suitable for IBS calculations for arbitrary distribution functions and its comparison to experimental data was presented in [37]. The IBS effect was also studied for high-brightness electron linac beams which appear to be non-Gaussian, especially in the longitudinal plane [38]. For low-emittance high-intensity electron storage rings, the interplay between intrabeam scattering and wake-field forces is discussed in [39]. In particular, the calculation of the IBS growth rates and the estimation of the emittance growth is discussed in detail in [40], where the importance of knowing the formation of the distribution tails is underlined, referring also to the use of Monte-Carlo methods.
In order to simulate the impact of a distribution shape on the emittance evolution when considering IBS and radiation effects the SIRE (Software for IBS and Radiation Effects) [5,6] code was developed by Vivoli and Martini at CERN. A similar Monte Carlo approach (IBStrack [41]) was implemented also in the collective effects simulation tool CMAD [42,43]. Both algorithms have as their basis MOCAC (MOnte CArlo Code), a Monte-Carlo code developed by Zenkevich and collaborators [44,45], which calculates the IBS effect for arbitrary distributions, by representing the beam through a large number of macro-particles occupying points in the 6-dimensional phase space. Being an extension of MOCAC, SIRE was developed to simulate the evolution of the beam particle distributions, taking into account the effects of IBS, synchrotron radiation and quantum excitation. For the evolution of the LHC bunch parameters and the shape of the bunch profiles presented in this paper, the simulations are performed using SIRE. As the physics and implementation details of the code are extensively explained in [6], we only briefly summaries some of its key features.

A. Software for IBS and Radiations Effects (SIRE)
For the IBS calculations, SIRE [5,6] uses the classical Rutherford cross section, which is closer to the Piwinski formalism [33]. It uses as input the optics functions at different locations of the lattice in order to determine the trajectories of the particles in phase space. Instead of using the 6 coordinates for position and momentum, the two Courant-Snyder and longitudinal invariants and the 3 phases (betatron and synchrotron) are used. For a linear ring, the 3 invariants are conserved between points around the lattice and can only be changed by the effects of IBS, synchrotron radiation and quantum excitation, while the phases are chosen randomly at each given point of the lattice. The time steps for which the IBS and radiation effects are called should be specified such that they are larger than the revolution time and smaller than the damping/growth times. Dividing the total time by the time steps shows how frequently the IBS, synchrotron radiation and quantum excitation routines are called. The quantum excitation is implemented by adding to the 6 coordinates of each macro-particle a Gaussian random noise component.
Depending on the elapsed time, the synchrotron radiation damping acts on the invariants of the macro-particles as an exponential decrement. The routine introduced for this reason is called after the calculation of the IBS effect at each iteration. Using small iteration time steps dt (which are much smaller than the damping times and for which the emittances change adiabatically), the evolution of the transverse emittance and energy spread due to the effects of IBS and synchrotron radiation can be obtained by solving the coupled differential equations: with x,y0 and σ p0 being the zero current (without the effect of IBS) equilibrium transverse emittances and energy spread, respectively. The τ x,y , τ p are the synchrotron radiation damping times and T x,y , T p the IBS growth times. The algorithm SIRE uses to calculate IBS is similar to that implemented in MOCAC, where the beam is represented by a large number of macro-particles occupying points in the 6-dimensional phase space. The default distribution defined in SIRE by using a random number generator, is the Gaussian and is given in action angle variables. In order to introduce a different distribution, the proper random deviates should be generated or the action angle variables of all macroparticles for the desired distribution should be provided. After specifying the total beam population and the number of macro-particles, the initial distribution can be tracked. The particle distribution in all planes can be saved as often as requested during the simulation time. Currently, the output file gives the evolution of the emittance in all planes for the specified time steps.
After providing the beam distribution and the optics along a lattice, the beam is geometrically divided according to the specified number of cells for each plane. The macro-particles are assigned to each cell according to their geometrical position. For each lattice point defined in the optics file, the 3 phases of each macro-particle are randomly chosen and the position and momentum of the macro-particles are calculated. Based on the classical Rutherford cross section, intra-beam collisions between pairs of macro-particles are calculated in each cell. The momentum of particles is changed due to scattering. The number of macro-particles and cells, i.e. the number of collisions each macro-particle experiences, is chosen so as to give accurate results for a reasonable computational time. The scattering angles for each collision are determined. In order to get the mean value of the emittance and momentum deviation changes for all particles, we have to integrate over the phase space volume of betatron coordinates, momentum deviations and azimuthal positions of the interacting particles. The beam distribution is then updated based on the new invariants of the macro-particles. For a specified number of time steps which practically shows how frequently the IBS and synchrotron radiation routines are called, the beam distribution is updated and the rms beam emittances are recomputed, giving finally as output the emittance evolution in time. The simulation proceeds to the next lattice point and continues until the end time is reached.
A lattice compression technique named "lattice recurrences" has been implemented to speed up the calculations [5]. Since the increase of the invariants due to IBS is linear to the first order in the traveling time along an element, the IBS kicks with optics functions differing less than a specified precision value are considered equivalent. For the corresponding group of elements, the IBS effect is evaluated only once, resulting to a reduction of the computational time.

B. The logarithmic Coulomb factor
The IBS growth times have a complicated dependence on the beam properties, due to the coupling of the three planes through dispersion. Some of these properties are the bunch charge and energy, the beam optics and the equilibrium rms horizontal, vertical emittances and the energy spread. The IBS growth times depend also on a logarithmic Coulomb factor which is used to include the contribution of events having a very large and very small impact parameter. The typical way of computing a log factor overemphasises the importance of the very small impact parameter events, for which the tails of the steady-state bunch distributions are non-Gaussian. In the high energy approximation by Bane [25], in order to describe the size of the core of the bunch, the Coulomb log factor is calculated as was first proposed by Raubenheimer [46], i.e. based on a boundary between the contribution to the core and the tails. In B-M, Bane and CIMP methods, the Coulomb factor is defined as the ratio of the maximum r max to the minimum r min impact parameter in the collision of two particles in the bunch, that is (log) ≡ ln(r max /r min ). For typical flat beams, the r max is taken to be equal to the vertical beam size σ y , while r min is taken to be r min = r 0 β x /(γ 2 x ), with r 0 being the classical particle radius. Then, the Coulomb factor can be written as: The formalism by Piwinski always seems to underestimate the IBS effect with respect to the other theoretical models. What diversifies Piwinski's method, is the different definition of the Coulomb factor. In that approach, the maximum impact parameter which is typically taken as the vertical beam size appears. In the high energy limit, with d being the maximum impact parameter, the Coulomb (log) for Piwinski can be written as [47]: where a = σx γ βx x . Comparing the (log) factors of Eq. (8) and (9), we find that condition d = 4σ y in order for the Piwinski approach agreeing with the other models. SIRE uses the "binary collision map" algorithm, conceived by Zenkevich, which allows to reduce the effects of the continuous time dynamical IBS system to a discrete time map in momentum space. For the binary collision events, the maximum impact parameter is taken as the beam height.

C. Benchmarking of IBS theoretical model with SIRE
The performance of hadron machines is limited by the IBS effect causing emittance growth. For lepton machines such as future linear collider Damping Rings, new generation light sources and B-factories, the IBS effect can also be predominant. It is thus important to study the IBS theories in the presence of synchrotron radiation and quantum excitation and benchmark the existing theoretical models and tracking codes with experimental data. In this way, the codes limitations can be identified so that to apply the necessary improvements in order to get better predictions for a machine's operation.
A benchmarking of the IBS theoretical models with Monte-Carlo codes is presented in [48] for lepton rings. The comparison between different theoretical models and SIRE is discussed for the Compact Linear Collider (CLIC) damping ring (DR) [49], having ultra-low emittances which are strongly dominated by IBS, in the presence of synchrotron radiation and quantum excitation. Results of this comparison are presented in Fig. 4, for one turn of the DR lattice. Due to the fact that in SIRE the generation of the distribution is based on a random number generator, the tracking simulations were performed several times, resulting in the one standard deviation error-bars (plotted in green). The classical formalism of Piwinski (red colored) and SIRE are in perfect agreement, as was expected, since SIRE uses the classical Rutherford cross section which is closer to the Piwinski formalism. The Bjorken-Mtingwa (black colored) and Bane (magenta colored) formalisms overestimate the effect compared to Piwinski method mainly to a mismatch of the Coulomb factor used in the different approaches (see Eq. (9)). The IBS theoretical models have been studied in detail and benchmarked with experimental data also for hadron beams over the years [26,27]. A comparison of the LHC data with simulations performed with SIRE is discussed in [7,13]. In this paper, the SIRE simulations, as well as the benchmarking with the B-M formalism and experimental data are discussed in detail for the LHC and the High Luminosity LHC (HL-LHC) [10].

IV. SIMULATIONS PERFORMED WITH SIRE FOR THE LHC
In order to understand the evolution of the bunch characteristics, based on the bunch profile observations, it is important to study the interplay between IBS and radiation effects (synchrotron radiation and quantum excitation) during the full LHC cycle. This is done using SIRE for two cases which are important for the current and future machine performance; the nominal BCMS [8,9] and the HL-LHC [10] parameters. Despite some blow-up in the LHC during the ramp, it is observed that the BCMS beam gives an increase in peak luminosity of around 20%. The HL-LHC is the major LHC upgrade aiming to increase integrated luminosity by at least a factor of 10 compared to the nominal LHC design value (from 300 to 3000 fb −1 ). In order to achieve that, the bunch population needs to be increased and the transverse beam size at the collision points has to be lowered.
Apart from the IBS and synchrotron radiation which are the dominant effects for the emittance evolution in the LHC, a combination of other diffusion mechanisms, like the beam-beam effect, electron-cloud, noise (due to the power converters, the transverse damper, the crab cavities, etc.), and other non-linearities cause emittance growth and/or particle losses [50]. Despite the fact that these mechanisms are not included in SIRE, it is possible to add empirically (i.e. based on observations) their contribution. Practically, there is the option of adding or complementing the variation of the bunch parameters in time. In fact, the simulation studies presented in this paper for the LHC are focused on the 3σ range of the particle distributions and therefore, mechanisms which concern the far tail regime are not taken into account as they are more important for the lifetime of the beam then on distribution evolution .

A. Reduced lattice
As mentioned earlier, one of the inputs required by SIRE are the optical functions along the ring. As the LHC is a very long accelerator of about 27 km, with a very large number of elements in the sequence (more than 11000), SIRE requires an extremely long computational time to track the distribution for all the elements along the ring. Aiming to reduce the computational time, a study was first performed in order to identify the optimal minimum number of critical IBS kicks around the lattice, without affecting the overall effect. The IBS growth rates were firstly calculated for the full optics of the LHC, using the IBS module of the Methodical Accelerator Design code (MAD-X) [24] which is based on the Bjorken-Mtingwa formalism. Figure 5 shows the IBS growth rates in the longitudinal (green), the horizontal (blue) and the vertical (magenta) plane. Taking into account the strong IBS kicks along the ring, various lattices with a reduced number of elements, including the case of the smooth lattice approximation, were tested. Then, using the IBS module of MAD-X, the emittance evolution was calculated for several sets of beam parameters to assure that the choice of the elements is valid for a wide range of regimes, for which the IBS impact may be weaker or stronger. Finally, the optimal lattice chosen consists of only 92 elements whose positions are denoted by red dashed lines in Fig. 5. Figure 6 shows the emittance (left) and the bunch length (right) growth during 30 min at injection energy, for the nominal BCMS beams, with initial emittance and 4σ bunch length that are respectively x0 = 1.5 µmrad and σ s0 = 1 ns, having a bunch population of 1.2 × 10 11 protons. The black solid line refers to the case of the full lattice, while the red dashed one to the reduced lattice with the 92 elements. The magenta dotted line corresponds to the case of the smooth lattice approximation for which a lattice with a unique element, having the optics that represent in the best possible way the mean optics of the full lattice, is considered. The agreement of the full and the reduced lattice is very good in all planes. On the other hand, by using the smooth lattice approximation the IBS effect is slightly underestimated, in particular, in the longitudinal plane. In the next, since the results for the reduced and the full lattice agree also in SIRE, the reduced lattice is used as an input for the simulation code. After choosing the optimal number of cells and macro-particles, the computational time in the case of the reduced lattice is almost 20 times shorter than the one of the full LHC lattice.

B. Convergence studies
For a specified set of input beam parameters, various scans should be performed for different combinations of number of macro-particles and cells in order to find the optimal values which provide a fast tracking and at the same time, guarantee that the scattering process leads to accurate results. In these terms, in order to avoid having a very small number of macro-particles per cell, the total number of cells is calculated based on the optimal minimum number of macro-particles per cell. For n x , n y , n z being the number of cells in the horizontal, vertical and longitudinal plane, respectively, it is assumed that in the transverse plane there is a correlation between the number of cells ratio and the beam sizes ratio, meaning that n x /n y =σ x /σ y . Therefore, for n mp being the total number of macro-particles, the number of macro-particles per cell is: A scanning of the total number of cells is performed for an example set of beam parameters to be used as an input for tracking. Based on Eq. (10), by keeping the total number of macro-particles constant, the different combinations of cell numbers determines the number of macro-particles per cell. Figure 7 (left) shows the dependence of the emittance variation (ratio of final versus initial value) in the horizontal (blue) and longitudinal (green) plane on the number of macro-particles per cell, for a specified time duration. The value of the number of macro-particles per cell after which the variation of the emittances in both planes remains constant is chosen as the optimal minimum value. After specifying this value, a scanning is performed for a fixed number of macro-particles, in order to choose the number of cells to be used, firstly in the longitudinal plane. Then, the number of cells in the horizontal plane can be calculated using Eq. (10), when knowing the ratio of the beam sizes in the transverse plane 1 . In Fig. 7, the dependence of the emittances variation is plotted versus the number of cells in the longitudinal (middle) and horizontal (right) plane, for a specific time duration. It can be noticed that the variation of the emittances remains constant after a certain number of cells in the longitudinal and horizontal plane that is, for the example set of beam parameters, 350 and 13 cells, respectively. Finally, the number of cells in the vertical plane can be calculated by n y = n x /(σ x /σ y ).
C. Benchmarking of SIRE with the B-M IBS theoretical model SIRE has the advantage to accept any type of distribution as an input. If requested, it also gives as output the distribution at any stage of the tracking. In order to benchmark the code with the analytical formulation of B-M for the LHC, a Gaussian distribution was tracked for two sets of bunch parameters which are summarized in Table I for both the injection (450 GeV) and collision energy (6.5 TeV). The first case corresponds to the nominal BCMS [8,9] LHC beams, having a significantly lower transverse beam size with respect to the nominal production scheme. The second case corresponds to the HL-LHC [10,51] parameters, for which the bunch population is very high. The input optics functions used for tracking correspond to the ones of the aforementioned reduced lattice.

At the LHC injection energy (450 GeV)
The evolution of the horizontal emittance (left), the vertical emittance (middle) and energy spread (right) after 1 h at injection energy (450 GeV), where the IBS effect is dominant, are presented in Fig. 8 for the nominal BCMS case and in Fig. 9 for the HL-LHC parameters. The red and the blue lines correspond to the analytical calculations of the MAD-X [52] IBS routine (based on the B-M formalism) and to the SIRE results, respectively. Due to the fact that in SIRE the generation of the distribution is based on a random number generator, the tracking simulations were performed several times, resulting in the two standard deviation spread that are plotted in light blue. Table II summarizes the IBS growth of the transverse emittances and energy spread, for the nominal BCMS and HL-LHC parameters, as computed by the SIRE code and the B-M analytical formalism in MAD-X.   In the horizontal and longitudinal plane the IBS effect is dominant, while in the vertical plane, it is minor. Even though the SIRE simulation algorithm and the B-M analytical formalism make use of different approaches to calculate the IBS effect (SIRE uses the classical Rutherford cross section which is closer to the Piwinski formalism), they seem to agree very well during the 1 h time at injection energy. In the longitudinal plane, there is a small difference observed for longer time-spans. Such differences can be explained by the fact that SIRE reshapes the beam distributions, in a self-consistent way, after each collisional process, while the B-M IBS formalism assumes Gaussian beam distributions throughout the calculation.
The variation of the initially Gaussian particle distributions within 1 h at injection energy is shown in logarithmic scale in Fig. 10 and Fig. 11 for the nominal BCMS and the HL-LHC case, respectively. The initial and final (after 1 h) distributions in the horizontal (left), vertical (middle) and longitudinal (right) plane, are denoted by blue and red stars, respectively. They are fitted with the Gaussian (dashed line) and the q-Gaussian (solid line) functions. The fitting results of the initial and final distributions are presented in Table III for the nominal BCMS case and in Table IV for the HL-LHC case.
As was expected from the results shown in Fig. 8-9 concerning the IBS growth, the horizontal and longitudinal rms beam sizes get larger as time evolves, while the vertical one does not change. The vertical distributions remain Gaussian since q ≈ 1. For both the nominal and the HL-LHC case, the q parameter of the horizontal and longitudinal distributions is decreased. This can be explained by the fact that, due to IBS, the core of the distributions is blown up in such a way that it covers up the initially Gaussian tails of the input distributions, which remain less affected. In the longitudinal plane the decrease in q is more significant for the HL-LHC case. This indicates that the stronger IBS is, the more the core is blown up. Since for a light tailed distribution (q < 1) the Gaussian fit overestimates the rms value, the resulted beam sizes are slightly larger than in the case of the q-Gaussian fit. Comparing the root mean square error (RMSE) values of the two fitting functions for the final non-Gaussian bunch profiles shows that the q-Gaussian fit is better, in particular for the horizontal plane.  Since at collision energy IBS becomes weaker and synchrotron radiation starts playing an important role, it is the interplay between these effects that determines the evolution of the bunch characteristics. In this respect, for the benchmarking of the B-M IBS theoretical model with SIRE at collision energy, apart from the IBS, the radiation effects (synchrotron radiation and quantum excitation) are also taken into account. It should be mentioned that for the results presented in the following plots the intensity is assumed to be constant (which is actually true for the few non-colliding bunches, during physics fills). Figure 12 shows the horizontal emittance (left), the vertical emittance (middle) and energy spread (right) evolution after 10 h at collision energy for the nominal BCMS case, while Fig. 13 shows the evolutions for the HL-LHC parameters. The red and the blue lines correspond to the analytical calculations of the MAD-X [52] IBS routine (based on the B-M formalism) and to the SIRE results, respectively. The two standard deviation error-bars for the simulation results are plotted in light blue. Table V   After a few hours at collisions, the B-M analytical formalism and the simulations start differentiating. In order to understand whether these differences are explained by the fact that SIRE reshapes the beam distributions after each collisional process and the B-M IBS formalism assumes always Gaussian beam distributions, the bunch parameters given by SIRE at 5 h are used as input for the IBS and synchrotron radiation calculations in MAD-X (Gaussian bunches). The red dotted lines in Fig. 12 and Fig. 13 represent the results of these tests. Even if giving as input to MAD-X exactly the same bunch parameters as in SIRE, there is clear divergence of the MAD-X results (red dotted lines) with SIRE right after the 5 h at collisions. This divergence is much larger than the one observed during the first hours at collisions. After 5 h at collision energy, the beam in SIRE has been reshaped enough so that IBS and radiation processes act differently as compared to Gaussian MAD-X distributions. Consequently, the differences observed between the B-M analytical formalism and the simulations are expected because MAD-X assumes always Gaussian distribution, in contrast to SIRE that takes into account the variation of the bunch shape throughout the calculation.
Due to the fact that the IBS effect is minor in the vertical plane, the strong synchrotron radiation damping mechanism leads to a clear reduction of the vertical emittance. However, the variation of the horizontal emittance and energy spread is determined by the interplay of IBS growth with synchrotron radiation damping. For the nominal BCMS parameters, these variations are very small after 10 h at collision energy (Table V). For the HL-LHC case, having the same initial horizontal emittance but the double bunch population compared to the nominal BCMS parameters (Table I), the IBS effect prevails over synchrotron radiation in the horizontal plane after almost 3 h (Fig. 13 (left)). As can be seen in Fig. 13 (right) this in not the case for the longitudinal plane, for which the initial bunch length of 1.2 ns compared to the 1 ns in the nominal case, renders IBS weaker than synchrotron radiation and, results in the decrease of the energy spread.  The evolution of the initially Gaussian (in all planes) particle distributions within 10 h at collision energy is shown in logarithmic scale in Fig. 14 and Fig. 15 for the nominal BCMS and the HL-LHC case, respectively. The initial and final (after 10 h) distributions in the horizontal (left), vertical (middle) and longitudinal (right) plane, are denoted by blue and red stars, respectively. They are fitted with the Gaussian (dashed line) and the q-Gaussian (solid line) functions. The fitting results of the initial and final distributions are presented in Table VI for the nominal BCMS VII: Initial and final (after 10 h) fitting results for the horizontal, vertical and longitudinal bunch profiles shown in Fig. 15, for the HL-LHC parameters case at collision energy (7 TeV case and in Table VII for the HL-LHC case. The RMSE values of the two fitting functions show that when the final bunch profiles are strongly non-Gaussian, the q-Gaussian fitting results should be considered. In this respect, the evolution of the particle distributions in all planes for the nominal BCMS and HL-LHC cases is discussed based on the q-Gaussian results. The horizontal beam sizes do not change after 10 h at collision energy because the blow up caused by IBS is balanced out by the synchrotron radiation damping. However, there is a transformation of the horizontal distributions' shape for which the tails become less populated (q < 1). In the longitudinal plane both the beam size and the q parameter are reduced, meaning that synchrotron radiation prevails over IBS and the core is blown up due to IBS giving underpopulated tails. In the vertical plane, the dominant synchrotron radiation damping results in a smaller beam size without changing much the formation of the tails, so the distribution remains Gaussian.

V. BUNCH PROFILE MEASUREMENTS IN THE LHC
The transverse diagnostic instruments for measuring the bunch profiles in the LHC are the betatron matching monitor [53], the Beam Gas Ionization (BGI) monitor [54], the Beam Gas Vertex (BGV) monitor, the Beam Wire Scanners (WS) [55] and the Beam Synchrotron Light Monitor (BSRT) [56]. Compatibly with high intensity and high energy operation, the BSRT is the only instrument offering non-invasive, continuous and bunch-by-bunch measurements of the LHC beams. The BSRT is calibrated with respect to the WS during dedicated low beam intensity runs. 2 The LHC is equipped with two synchrotron radiation monitors (one per beam) used to characterize the transverse and longitudinal beam distributions. Due to the significant diffraction patterns at the tails of the BSRT transverse distributions there is no clear picture about the shape of the tails and so, they are often assumed to be Gaussian.
A parameter that is generally used to measure the longitudinal emittance in circular accelerators is the bunch length. The bunch length is given by the projection of the distribution function on the phase axis, which is known as the bunch profile or line density. It is operationally measured by the LHC Beam Quality Monitor (BQM) [57] which uses a wall current monitor pick-up (WCM) [58] to acquire the longitudinal profiles. Additionally, the longitudinal synchrotron radiation monitor (BSRL) [59], which uses the same synchrotron light source as the BSRT, continuously measures the longitudinal distribution of charges in the beams. The scopes connected to the WCM pick-ups can acquire longitudinal bunch profiles of both beams during a full LHC cycle.
A. Comparison between experimental data, the SIRE and the B-M analytical formalism, at the LHC collision energy (6.5 TeV) The longitudinal bunch manipulations performed during the ramp to avoid instabilities due to the loss of Landau damping [16], produce bunches that arrive at collision energy with a clearly non-Gaussian longitudinal shape. In addition, the transfer functions of the pickups and cables were measured and are used for deconvolution [60], resulting in some cases in tails which are asymmetric or have ripples. By assuming that these profiles are Gaussian may lead in underestimating or overestimating the actual bunch length. For the studies presented in this paper, these profiles are fitted using the q-Gaussian function. An example showing the evolution of the q parameter for the longitudinal profile of a bunch train during 11.5 h at collisions (6.5 TeV) in the LHC is presented in Fig. 16. It is clear that with such q parameter values, corresponding to non-Gaussian tails, the rms beam size cannot be accurately estimated by using the Gaussian function. The increase of the q parameter means that the longitudinal distributions with the underpopulated tails (q < 1) at the start of collisions, become more Gaussian (q → 1) as time evolves. This is a general statement that can be made for the longitudinal distribution observed at the collision energy of the LHC. The evolution of the longitudinal particle distribution of a single bunch that is picked out of the train of bunches is shown in Fig. 17 for the time period of 11.5 h. The initial bunch profile (plotted in blue) is fitted with the Gaussian and the q-Gaussian functions that give different rms beam sizes because of the dependence of the standard deviation on the q parameter (Eq. (A5)). The fitting results are used to generate a Gaussian and a q-Gaussian distribution to be tracked in SIRE in order to compare the experimental observations with the results of the code. In Fig. 18, the initial (at the start of collisions) and the final (after 11.5 h) longitudinal bunch profiles, as observed in the LHC (left) and as calculated by SIRE (right) for an initially q-Gaussian simulated profile, are denoted by blue and red stars, respectively. They are plotted in logarithmic scale and they are fitted with the Gaussian (dashed line) and the q-Gaussian (solid line) functions. The reduction of the bunch population with time due to burn-off and the extra (on top of IBS) transverse emittance blow-up observed in the machine, are taken into account for the simulation. The transverse distributions are assumed to be Gaussian, since at collisions the shape of their tails is not clear due to diffraction. The fitting results are presented in Table VIII. Even if there seems to be no change at the tails of the simulated distribution, in reality the profiles become more Gaussian. Within 11.5 h at collisions, the rms beam size of the measured bunch profile and of the corresponding tracked distribution is reduced by 21% and by 19%, respectively. This shows a very good agreement between the experimental data and the simulations performed with SIRE. Figure 19 shows in logarithmic scale the initial (blue stars) and the final (red stars) horizontal (left) and vertical (right) bunch profiles as calculated by SIRE, fitted with the Gaussian (dashed line) and the q-Gaussian (solid line) functions. As can be seen in Table IX, the simulations showed no change in the transverse beam sizes and that is because the extra (on top of IBS) transverse emittance blow-up is included. The effect of IBS together with the extra blow-up assumed, widens the core of the horizontal bunch in such a way that the q parameter is decreased by around 10% within these 11.5 h. Since IBS is negligible in the vertical plane, the fact that the vertical bunch profile remains Gaussian indicates that the interplay between the synchrotron radiation damping and the extra blow-up do not change the tails of the distribution.
The 4σ-bunch length evolution when assuming Gaussian (left) and q-Gaussian (right) initial distributions is shown in Figure 20. The blue line corresponds to the SIRE calculations and the red line to the results given by the IBS module of MAD-X [24] which is based on the analytical formulation of B-M and always assumes Gaussian distributions. The bunch length evolution, together with the two standard deviation error-bars, when fitting the data with the Gaussian and the q-Gaussian functions is represented by a black and a grey line, respectively. The bunch length values differ for the two distribution functions used due to the fact that, for a light tailed distribution the rms value is overestimated by fitting a Gaussian. When assuming a Gaussian distribution, the bunch length evolution calculated by the code gets closer to the measured data. For the q-Gaussian case the agreement between data and simulations is excellent. This is a remarkable result, taking into account the fact that no assumptions are being made in the simulations apart from identical initial conditions with respect to the experimental ones. In agreement with the results presented in the previous section, the divergence between the SIRE and the MAD-X for longer time-spans is something to be expected since the distribution shape in SIRE is updated, while in MAD-X it remains Gaussian. At collisions, the divergence between the luminosity model [1] and the measured luminosity by the experiments becomes more pronounced as time evolves [3]. Actually, the predicted luminosity by the model is always larger compared to the measured one by the end of collisions. As calculated by SIRE, the weight of the horizontal bunch profile tails is decreased in time (see Table IX) and as explained in Section II, for lighter tails and a constant beam size, the luminosity is expected to become lower. It is then clear that by taking into account the luminosity change due to the variation of the transverse distribution tails, the model predictions can be significantly improved.

VI. SUMMARY AND OUTLOOK
In the LHC, the interplay between a series of effects can lead to distributions with non-Gaussian tails. Since the rms value of a distribution can be underestimated or overestimated by using a simple Gaussian function, the use of appropriate fitting functions to accurately estimate the beam size and the behavior of the tails is necessary. The impact of non-Gaussian distribution shapes on the estimated luminosity is discussed. One of the next steps is to improve the luminosity model, that is currently based on Gaussian distributions, by taking into account the actual shape of the bunch profiles. In this way, it is possible to get more accurate luminosity predictions. Already, for the operational scenario of the High Luminosity LHC upgrade [51], a non-Gaussian bunch length estimation is being considered.
The way IBS and radiation effects act depends on the shape of the bunch profiles. Aiming to quantify the impact of the distribution's shape on the emittance evolution, a multi-particle tracking code called SIRE, is used. The benchmarking of the B-M analytical formalism with SIRE showed a very good agreement for the first couple of hours at the injection (450 GeV) and collision (6.5 TeV) energies of the LHC, even if they make use of different approaches to calculate the IBS effect. The differences observed for longer time-spans are expected, since in SIRE the particle distributions are updated, while MAD-X always assumes Gaussian distributions. The results obtained from the simulations encourage the idea of using the code for tracking distributions coming from experimental data, in order to study the impact of the distributions shape on the evolution of the bunch characteristics. After the comparison with experimental data, the fact that SIRE takes into account the change of the particle distribution showed that it is a very useful tool for estimating the actual bunch parameters evolution in the machine. The contribution of effects such as betatron coupling, noise and electron-cloud, to the emittance growth are planned to be included in the simulation code, in order to complement the existing semi-analytical emittance evolution model.

Appendix A: The q-Gaussian distribution function
The q-Gaussian [17] which is used to describe more accurately bunch profiles with tails that differ from the ones of a normal distribution, has a probability density function given by: The q-exponential function is given by: The parameter q describes the weight of the tails, in the sense that the larger its value, the heavier the tails become, as presented in Fig. 21. In the limit of q → 1, the normal distribution is obtained. The distribution is characterized as "light" tailed when q < 1 and as "heavy" tailed when q > 1. The normalization factor C q differs for specific ranges of the q parameter and it is written as: , for − ∞ < q < 1 √ π , for q = 1 , for 1 < q < 3 .

(A3)
The parameter β qG is a real positive number. As the normal distribution, the q-Gaussian is an even function taking its maximum at x = 0, where For a certain q value, the higher is the value of β qG , the larger is the maximum of the probability density function, as can be observed in Fig. 21. The standard deviation which also differs for specific ranges of the q parameter, is given by: σ qG = 1 β qG (5−3q) , for q < 5/3 ∞ , for 5/3 ≤ q < 2 undefined , for 2 ≤ q < 3 . (A5) In the heavy tail regimes, the distribution is equivalent to the Student's t-distribution with a direct mapping between q and the degrees of freedom ν (Eq. (A6)). Statistically the q-Gaussian is a scaled re-parametrization of the Student's t-distribution [61] for which the parameter ν is constrained to be a positive integer related to the sample size. The advantage of the q-Gaussian function is that, by introducing the parameters q and β qG , a generalization of the Student's t-distribution to negative and non-integer degrees of freedom is possible, where: q = ν + 3 ν + 1 with β qG = 1 3 − q . (A6)