1 Introduction

After the centre-of-mass energy, perhaps the most important figure of merit of any particle collider is the time integral of the luminosity delivered to the particle or nuclear physics experiments. Generally, integrated luminosity is a measure of the data collected and defines the precision of measurements and the potential for new discoveries.Footnote 1

The product of luminosity \(\mathcal {L}\) and the appropriate cross section is the rate of events of a specific type. During a fill of a circular collider, the phase-space distributions and intensities of the beams evolve in time under the combined influences of several inter-dependent physical processes, e.g. the beam–beam collisions (or luminosity losses), intra-beam scattering (IBS), and synchrotron radiation damping. As a consequence of the intensity losses, \(\mathcal {L}\) eventually decreases with time. To maximise \(\int \mathcal {L}\,dt\), non-collisional losses must be minimised and the beam sizes should be kept small. To understand, predict, and optimise the physics data collected, it is therefore very important to have reliable models of the quantitative time evolution of the beams and \(\mathcal {L}\). In this article, we present two such models: Collider Time Evolution (CTE) and the Multi-Bunch Simulation (MBS). They have been developed for heavy-ion operation at the CERN Large Hadron Collider (LHC) [1], although they are also applicable to other hadron colliders.

There have been numerous studies of the time evolution of beams and collider luminosity in the past, e.g. those in Refs. [2,3,4,5,6,7,8,9,10,11], with similar applications to this article. Most of them rely on the numerical solution of systems of coupled ordinary differential equations (ODEs) that model the time evolution of a few key beam parameters, most commonly the transverse beam emittances, the bunch lengths, and intensities. These models use additional assumptions on the beam distributions, which are typically assumed to be Gaussian in all three degrees of freedom. The beam distributions can also be tracked by solving Fokker–Planck equations for distribution functions, as in Refs. [5, 6]. Past studies of the heavy-ion luminosity at the LHC were carried out in Ref. [12] using ODEs and in Ref. [10], where a particle tracking simulation was used in addition to the ODE method. That study used an early version of the CTE code, which has since been further developed and significantly improved.

Since the LHC is used for the benchmark, we give first a brief introduction to heavy-ion operation of the LHC in Sect. 2. Then, we present the simulation codes, CTE in Sect. 3 and MBS in Sect. 4, and compare their results to data from a large number of fills in the 2018 LHC Pb–Pb run in Sect. 5. Finally, we use the simulations to predict the performance in future LHC heavy-ion runs in Sect. 6.

2 Heavy-ion operation of the LHC

Table 1 Pb beam parameters at the start of collisions in the LHC, as foreseen in the LHC design report (two experiments illuminated) [1], as achieved in 2018 [14, 15], and as envisaged for future runs in Run 3 (presently planned from 2022) and in HL-LHC [13, 16]
Full size table

The LHC is a circular collider, designed to collide protons and nuclei at a beam energy of \(7\,Z\) TeV; so far, up to \(6.5\,Z\) TeV has been achieved.Footnote 2 The LHC uses two counter-rotating beams, called B1 and B2, that collide at four interaction points (IPs), inside the experiments ATLAS [17] at IP1, ALICE [18] at IP2, CMS [19] at IP5, and LHCb [20] at IP8. Typically, the LHC has operated for about one month per year with heavy-ion beams, mainly fully stripped Pb nuclei, in addition to the main physics programme with proton–proton collisions. The initial aim was mainly to provide Pb–Pb collisions to ALICE, which is specialised in heavy-ion physics, but over time all the LHC experiments have joined the heavy-ion collision programme.

Four Pb–Pb runs were carried out so far (in 2010, 2011, 2015, and 2018) [15, 21, 22], resulting in \(\int \mathcal {L}\,dt\) \(\approx 1.5\) nb\(^{-1}\) at ALICE, thus surpassing the initial target of 1 nb\(^{-1}\) for the first 10 years of operation [18] in spite of more experiments than initially foreseen sharing the luminosity and a lower-than-nominal beam energy ( in Run 1 (2010–2013)Footnote 3 and in Run 2 (2015–2018), as opposed to the nominal ). So far, ALICE needed to be levelled at in order not to exceed the event rate limit of the detector. The high-luminosity experiments ATLAS and CMS collected in total about 2.5 nb\(^{-1}\) in Run 1 and Run 2, since no luminosity levelling was needed, or only done at a much higher value [14]. LHCb was the last experiment to join the heavy-ion programme, integrating about 0.25 nb\(^{-1}\) in Run 2.

In addition to Pb–Pb collisions, a completely new mode of operation with proton–nucleus collisions was put in place [23,24,25,26]. Runs with p–Pb, not foreseen at the LHC design stage, were carried out in 2012, 2013, and 2016. About 250 nb\(^{-1}\) was gathered in ATLAS and CMS and 75 nb\(^{-1}\) in ALICE.

The key parameters of the LHC in Pb–Pb mode are shown in Table 1, where we highlight first the design parameters [1] and those achieved in 2018 [15]. Notably, most of the LHC design parameters have been reached or surpassed; in particular, the 2018 luminosity was more than a factor 6 higher than the design value. This was achieved mainly thanks to a factor 3 larger bunch population from the injectors and a larger number of bunches. A new filling pattern with a shorter bunch spacing of 75 ns, compared to the previous 100-ns spacing, was adopted halfway through the 2018 run.

After the very successful first two Runs, the heavy-ion programme is scheduled to continue in the future, typically with short runs of about 1 month per year of Pb–Pb or p–Pb collisions. The projected future parameters [13] are shown in Table 1. The physics goals and luminosity targets of the experiments are specified in Refs. [27, 28]. The LHC performance in these future runs will benefit from upgrades of the LHC injectors [29] and the high-luminosity LHC (HL-LHC) [30]. Thanks to the injector upgrades, an even shorter bunch spacing of 50 ns will be achievable in the future. To safely store the higher beam intensity without an increase of spurious beam dumps and quenches, it is planned to upgrade the LHC collimation system [31,32,33,34]. Upgrades of the ALICE detectors [35] will allow levelling at a significantly higher hadronic event rate of 50 kHz corresponding to a luminosity of . The feasibility of this luminosity also relies on the installation of new collimators [30] to intercept collision products that risk to quench impacted magnets, in particular from bound-free pair production (BFPP) and to a lesser extent from electromagnetic dissociation [12, 36,37,38,39,40]. The luminosity in LHCb will still be limited by the BFPP losses, since no new collimators are foreseen to be installed in IR8.

All upgrades of the injectors, collider, and experiments relevant to the heavy-ion programme are planned to be implemented for the machine start-up after long shutdown 2 (LS2), presently foreseen for 2022, making the full heavy-ion “HL-LHC” performance available [13]. A proposal for future running with lighter nuclei has also been put forward [27].

3 Collider Time Evolution (CTE)

3.1 General features

The CTE code, initially presented and benchmarked with data from the Relativistic Heavy-Ion Collider (RHIC) in Ref. [10], and further developed in Refs. [41, 42] and for this article, tracks bunches of macro-particles using a one-turn map. Typically, we simulate \(2\times 10^4\) macro-particles per bunch, which are tracked in a 6D phase space. On every turn, various physical effects are sequentially applied to the bunches and the particle coordinates are updated accordingly.

The physical effects in CTE, which can be turned on individually in any combination, include betatron motion, longitudinal motion, beam–beam collisions with optional luminosity levelling, IBS, radiation damping and quantum excitation, machine aperture, extra generic losses modelled through a non-collisional lifetime, generic emittance blowup modelled through a rise time, and stochastic cooling.Footnote 4

The user can choose to simulate a smaller number of turns than the real machine turns to increase the computational speed. The strengths of the physical processes, e.g. collision probabilities or amplitude changes from IBS or radiation damping, are then scaled up according to the real time of the simulated fill. For the results shown in this article, we typically substitute \(4\times 10^8\) machine turns, corresponding to about 10 h of LHC operation, by \(2\times 10^4\) simulation turns.

The strengths of the physical processes are scaled further to account for the smaller number of macro-particles in the simulated bunches than real particles in the LHC bunches. Mostly, we use \(2\times 10^4\) macro-particles per bunch. This gives a fast execution time of typically a few minutes for a 10-h LHC fill on a standard desktop PC, without notable loss in precision.

To model the transverse and longitudinal motion, the particle coordinates are updated deterministically through a one-turn map based on the machine tune, chromaticity and longitudinal parameters such as RF voltage and momentum compaction [10]. Any particles outside of the stable area of the RF bucket or a user-defined transverse aperture cut are removed. CTE also includes linear betatron coupling in a thin-lens approximation.

In CTE, radiation damping and quantum excitation in each plane u are modelled as a combination of a random excitation and a deterministic decay, with the latter given by the emittance damping time \(T_{\mathrm {rad},u}\). We expand the damping coefficient to first order in \(T_\mathrm {rev}/T_{\mathrm {rad},u}\) with \(T_\mathrm {rev}\) being the revolution time, as described in detail in Ref. [10, 43]. It should be noted that in the LHC, the blowup from IBS is much stronger than the quantum excitation.

3.2 Beam–beam collisions

To model the collisions between two bunches, we calculate a collision probability \(P_1\) for each macro-particle, and a sampled random number determines whether the particle is removed. \(P_1\) can be calculated either through a numerical solution of the full overlap integral of the single particle trajectory with the opposing bunch distribution, or through a faster approximate method, where the opposing transverse bunch distribution is fitted to a Gaussian, although the actual longitudinal distribution is still accounted for. This allows some of the integrals to be evaluated analytically, resulting in a better computational performance. The detailed derivations and resulting equations are given in Ref. [10].

The collision schedule, determined by the bunch filling pattern, can be modelled either through a fast and simplified approach, where the full beam is represented by one macro-bunch and the collision probability at each IP is scaled by the number of real bunches colliding, or through a more detailed, but slower approach where one macro-bunch is tracked per beam and per beam–beam equivalence class, as defined in [44]. For the CTE studies in this article, we rely on the first approach with one macro-bunch per beam and the approximated \(P_1\), since it provides speed and simplicity in the set-up and, as shown in Sect. 5, the results are nevertheless in very good agreement with LHC data.

The effect of emittance blowup from core depletion in the collisions, first discovered in Refs. [10, 45] and further developed in Ref. [46], is automatically accounted for. As in the real machine, the collision probability in CTE is higher for macro-particles in the core, meaning that more particles are removed in the core than in the tails, and hence, the effective emittance grows.

3.3 Intra-beam scattering

The effect of IBS is modelled by giving each macro-particle a random kick per turn. The distribution of the kicks is inferred from the instantaneous growth times in the transverse and longitudinal planes, which are re-calculated on every turn. The growth times are first calculated approximating the bunch distributions as Gaussian, which allows the use of the well-established calculations detailed below. However, in order to account for non-Gaussian longitudinal profiles, the longitudinal kicks are modulated by the local particle density following the method in Refs. [8, 9]. This allows the evolution of highly non-Gaussian longitudinal bunch profiles, as in RHIC [10]. Optionally, the user can specify a mixing factor between the horizontal and vertical growth rates; as this factor is small although not well known for the LHC, we have generally set it to zero.

Several IBS models are available to choose from, which we denote

  • Piwinski smooth [47],

  • Piwinski lattice [47], similar to Piwinski smooth but using the full lattice and optics instead of a smooth approximation,

  • modified Piwinski [48],

  • Bane [49], a high-energy approximation of the well-known,

  • Bjorken–Mtingwa model [50],

  • Nagaitsev [51], a re-formulation of Bjorken–Mtingwa, with some assumptions that allow faster numerical calculations.

  • the user can also provide an external file with pre-calculated IBS growth rates from any external model, sampled on a grid of bunch dimension values, from which CTE makes an online interpolation.

Table 2 Calculated IBS rise times \(T_\mathrm {IBS}\) for the emittances \(\epsilon \) (with rise times for the bunch dimensions being a factor 2 larger) for the Pb beam parameters in Table 1 assumed for future LHC runs
Full size table

An in-depth comparison of the IBS models is beyond the scope of this paper, but we show as an example in Table 2 the growth times obtained for the future Pb beam parameters at the start of a fill in Table 1. The CTE calculations are compared with off-line calculations using both the MAD-X program [52], which uses a development of the Conte–Martini model [53], and a standalone implementation of Bjorken–Mtingwa. It can be seen that the latter gives very similar results to both Bane and Nagaitsev, as well as modified Piwinski. The growth time in the longitudinal plane is about a factor 2 shorter than in the horizontal plane. Most models with flat, uncoupled beam optics give a very long negative vertical growth time which is typically unimportant for the dynamics and not given. For the sake of comparison, the energy spread has been forced to the value in Table 1 in all computations in Table 2, while the energy spread calculated internally during CTE simulations is typically slightly shorter, since CTE exactly matches the longitudinal phase space based on the input bunch length and the longitudinal Hamiltonian. For all the studies shown in this paper, we have used the Nagaitsev model.

4 Multi-Bunch Simulation (MBS)

The second simulation code, MBS, is based on a very large set of coupled ODEs, which model the evolution of the bunch parameters of every single bunch in both beams. It can therefore model the real collision schedule according to the applied filling scheme. On the other hand, it relies on the assumption of Gaussian bunch distributions, as opposed to CTE where the tracked particles can assume any distribution.

For bunch j, the ODEs for the bunch intensity \(N_j\) and RMS beam dimensions \(\sigma _{u,j}\) in plane u are implemented as [54]:

$$\begin{aligned}&{\dot{N}}_j(t)=- \sum _{i\in \mathrm {IPs} }\left( \sigma \mathcal {L}_{ij} + \frac{N_j (t)}{T_\mathrm {IP,i}}\right) - \frac{N_j (t)}{ T_\mathrm {nc}} \, , \end{aligned}$$
(1)
$$\begin{aligned}&{\dot{\sigma }}_{u,j}(t)=\left( \frac{1}{2 T_{u,\mathrm {IBS}}} - \frac{1}{2 T_{u,\mathrm {rad}}} \right) \sigma _{u,j}(t) \, . \end{aligned}$$
(2)

Equation (2) holds for the planes \(u=x,y,z\), giving in total four ODEs per bunch, which translates into more than 5000 ODEs for the fully filled LHC machine. As shown in Eqs. (1)–(2), MBS includes the effects of luminosity burn-off, IBS through the rise time \(T_{u,\mathrm {IBS}}\), radiation damping through the damping time \(T_{u,\mathrm {rad}}\), and generic extra losses through a lifetime \(T_\mathrm {nc}\) for all bunches and \(T_\mathrm {IP,i}\) for IP-dependent losses. As we define \(T_{u,\mathrm {rad}}\) and \(T_{u,\mathrm {IBS}}\) to refer to changes in the emittances, a factor 2 is included in Eq. (2) to get the corresponding changes for the bunch dimensions. Note that \(\mathcal {L}\) is calculated per IP and bunch and is a function of the parameters of the collision partner for bunch j at IP i, denoted by the ~ symbol and subscript ij as

$$\begin{aligned} \mathcal {L}_{ij}=\mathcal {L}_{ij} (\sigma _{x,j},\sigma _{y,j},\sigma _{z,j},N_j,{\tilde{\sigma }}_{x,ij},{\tilde{\sigma }}_{y,ij},{\tilde{\sigma }}_{z,ij} {\tilde{N}}_{ij}). \end{aligned}$$
(3)

This mechanism effectively couples the evolution of all bunches via Eq. (1). The determination of the collision partners relies on the input of the full filling pattern.

MBS calculates \(T_{u,\mathrm {rad}}\) with the usual formalism including the radiation integrals [55, 56], while \(T_{u,\mathrm {IBS}}\) is calculated with the so-called completely integrated modified Piwinski model (CIMP) [57, 58]. It accounts for the full lattice and can be used for very fast numerical evaluations. Earlier comparisons have shown a very good agreement with Bjorken–Mtingwa [57]. As given in Table 2, it gives for the LHC a similar longitudinal growth as the Nagaitsev model in CTE, while the horizontal growth time is about 10% shorter, meaning a slightly stronger IBS effect. The user can apply ad hoc correction factors to manually alter the IBS strength and a parameter to mix the horizontal and vertical growth rates. This has not been done in this article.

In MBS, Eqs. (1)–(2) are integrated numerically using the explicit Euler algorithm as discussed in Ref. [54]. The time steps \(\Delta t\) can be specified and in this work we use .

Luminosity levelling can be activated as in CTE, in which a scaling factor \(S_i=\text {min}(1,\mathcal {L}_{ti}/\mathcal {L}_{\mathrm {pot},i})\) is applied to \(\mathcal {L}_{ij}\) in Eq. (1) to keep it at a constant, user-defined target \( \mathcal {L}_{ti}\) at IP i. Here, \(\mathcal {L}_{\mathrm {pot},i}\) is the maximum potential luminosity at IP i calculated over all bunches.

Additional, generic losses can be accounted for through the non-collisional lifetime \(T_\mathrm {nc}\). Bunches colliding at IP i can be assigned an extra IP-specific lifetime \(T_\mathrm {IP,i}\), which has been shown to increase the agreement with early data [54].

Several physical effects are neglected in both MBS and CTE. Firstly, beam–gas collisions could, apart from the risk of increased experimental backgrounds [59, 60], lead to both emittance growth and intensity decay [56]. However, with the very good vacuum levels achieved in the LHC, these processes are weak compared to the strong burn-off in the collisions and the emittance growth from IBS. However, beam–gas interactions are likely responsible for a substantial part of the observed non-collisional losses discussed in Sect. 5 and can easily be included in both CTE and MBS as a generic non-collisional lifetime.

Other influences on the beam include noise in feedback systems, RF cavities or magnet power supplies, and dynamic aperture losses due to magnet nonlinearities and beam–beam effects. These effects depend on imperfections that are not well known in quantitative detail. As the simulation result for Pb operation is dominated by burn-off and very close to the measurements, as shown in Sect. 5, we conclude that these effects are minor in LHC heavy-ion operation. Other studies suggest that they may be important for protons [61]. Nevertheless, losses due to these additional effects can be subsumed in the generic non-collisional lifetime in both CTE and MBS, and any effect on the emittance can be included in the generic emittance risetime in CTE.

5 Analysis of the 2018 Pb–Pb run

5.1 Simulation set-up

To compare the simulations with data, 30 out of the 46 physics fills in the 2018 LHC Pb–Pb run [15] at a beam energy of were simulated in detail with CTE and MBS. The remaining 16 fills were de-selected due to either missing logged data, very short fill lengths, or non-standard operational procedures interfering with the evolution. The simulations cover the so-called stable-beam period in each fill, which starts at the point in the operational cycle where the beams are brought into collision and left evolving and it continues until the fill is terminated by a beam dump.

The starting conditions for the simulations were taken from the 2018 operational parameters in Table 1, except for the parameters that varied between fills, which were instead extracted from logged data. The starting intensities, emittances, and bunch lengths varied due to fluctuations in the injector performance as well as small variations between injection and stable beams in the LHC. For CTE, using the approach with a single macro-bunch per beam, the average bunch parameters were used, while for MBS, individual values were extracted for each bunch. Furthermore, the filling scheme and hence the number of bunches colliding at each IP were changed several times in order to first have a gradual ramp-up and then deploy the new 75-ns scheme in the second half of the run. The luminosity levelling target was also gradually increased from one fill to the next, in order to carefully probe the limits from losses due to BFPP.

From the starting point of the colliding beams, the simulated beam parameters and luminosity evolved independently and no further input from measurements was used. Since non-colliding bunches in the machine showed on average 100-h lifetime, we applied the same non-collisional lifetime in the simulations but did not include any IP-dependent lifetimes in MBS or any additional emittance blowup. A total burn-off cross section for particle removal was assumed, consisting of the contributions shown in Table 3. As can be seen, it is dominated by the electromagnetic processes (BFPP and EMD), and the hadronic part is a minor contribution.

Because of the very large \(\sigma _{\mathrm {tot}}\), the LHC Pb ion operation is in a strong burn-off regime, where the total number of injected Pb ions, \(N_{1,2}\) in either beam, determines the maximum possible integrated luminosity,

$$\begin{aligned} \sum _{\mathrm {expts}} \int _{0}^{T_\mathrm {f}} \mathcal {L}(t)\, {\mathrm{d}}t \le \frac{\min (N_1,N_2)}{\sigma _{\mathrm {tot}}} \end{aligned}$$
(4)

where equality is approached in the limit of vanishing non-collisional losses and exhaustion of the lesser beam as the fill length \(T_\mathrm {f}\rightarrow \infty \). In typical LHC heavy-ion fills, the ratio of the two sides of this inequality, the luminous efficiency, exceeds 50–60%.

Table 3 Burn-off cross sections for various interactions between colliding Pb–Pb beams at , as in the 2018 LHC operation, for future Pb–Pb operation at , and for future 7 p–Pb collisions
Full size table

During the first half of the 2018 run, there was an error in the local coupling correction around IP2 [67], which significantly reduced the ALICE luminosity but was later corrected. For the affected fills, we used an effective \(\beta ^{*}\)-value of 0.9 m in the simulations. This value, which reproduces the measured luminosity in the simulated fills, is based on the luminosity scans in [68]. It is the \(\beta ^{*}\) that would give the same emittance in the ALICE scans as the value inferred from the scans at ATLAS and CMS.

Fig. 1
figure 1

The measured evolution (solid lines) of key observables (instantaneous luminosity, integrated luminosity, average bunch intensity \(N_b\), and average transverse emittances and bunch lengths) during two typical Pb–Pb fills (7453—left, 7472—right) from the 2018 LHC run, compared to simulation results from CTE (dashed lines) and MBS (dotted lines). The IP5 luminosity is not shown as it is almost identical to the one in IP1

Full size image
Fig. 2
figure 2

The measured evolution (solid lines) of key observables (instantaneous luminosity, integrated luminosity, average bunch intensity \(N_b\), and average transverse emittances and bunch lengths) during two typical Pb–Pb fills (7477—left, 7490—right) from the 2018 LHC run, compared to simulation results from CTE (dashed lines) and MBS (dotted lines). The IP5 luminosity is not shown as it is almost identical to the one in IP1

Full size image

5.2 Results

Results for a few typical fills with different starting conditions are presented in Figs. 1 and 2. For each fill, the simulated and measured online evolution of the instantaneous and integrated luminosity can be seen, as well as the beam intensity, and an excellent agreement is found in these quantities. It should be noted, though, that the end of the luminosity levelling at IP2 occurs slightly later in MBS than in CTE or the logged data. The other simulated fills are very similar in the level of agreement.

A fair agreement is found also in the emittance and bunch length evolutions. This comes in spite of an uncertainty on the measured emittances from the synchrotron light monitor (BSRT), which was never calibrated in detail for Pb beams in the 2018 run, and which sometimes has missing data as, e.g. in fill 7477 in Fig. 2. Nevertheless, using the measured emittance values as initial conditions at the start of each fill results in a very good agreement in the key observables in Fig. 1, strongly suggesting that the real error on the emittance measurement is small.

The horizontal emittance generally grows initially, while the vertical emittance shrinks, since radiation damping dominates in the vertical plane. The simulated horizontal emittance from CTE typically grows more than the measured one, while the simulated vertical emittance shrinks more. This is likely due to the assumptions of zero IBS mixing between the planes in the simulation. An IBS mixing value can be empirically fitted to generally yield a good agreement with the measured emittances as shown in Fig. 3 and without significantly affecting the simulated luminosity. However, as the real value in the machine is unknown we have refrained from doing so in the other results.

Fig. 3
figure 3

Example of simulated emittances from CTE for fill 7472, shown in Fig. 1, but with using an empirically fitted mixing value between the IBS growth rates in the horizontal and vertical planes. We used \(T_{x,\mathrm {IBS,eff}}\) =0.65\(T_{x,\mathrm {IBS}}\) + 0.35 \(T_{y,\mathrm {IBS}}\) with the analogous mixing in the vertical plane

Full size image
Table 4 Calculated emittance rise times and damping times for IBS, radiation damping, and core depletion, as well as the total effective rise time over all processes i, calculated as \(1/T_\mathrm {eff}=\sum _i 1/T_i\)
Full size table

In MBS, the emittances stay smaller than in CTE, in spite of a slightly stronger IBS growth (Table 2). This is due to the core depletion, as can be understood from Table 4, where we compare the strengths of the main effects. For the horizontal plane in a typical 75-ns fill with 733 bunches, IBS is largely counteracted by radiation damping at the start of the fill, with a residual combined effective risetime \(T_{x,\mathrm {eff}}\) \(\approx \)17–18 h. However, if the core depletion rise time of about 28 h is also included, the effective total growth time goes down to about 11 h. Since core depletion is included by construction in the CTE collision routine, the effective emittance grows much more than in MBS which does not include this effect. On the other hand, the bunch length grows significantly more in MBS than in CTE, with CTE showing a better agreement with the measurements. This is driven mainly by stronger longitudinal IBS due to the smaller transverse emittance.

Fig. 4
figure 4

Distribution of the ratio of simulated (with CTE—left, MBS—right) to measured integrated luminosity for ATLAS, ALICE and LHCb per fill in the 2018 Pb–Pb run

Full size image

In practical terms, the most important quantity for the simulation benchmark is the integrated luminosity at the end of each fill. The ratios of \(\int \mathcal {L}\,dt\) between simulations and measurements in the 30 simulated fills are shown in Fig. 4 and the mean and standard deviations per experiment in Table 5. For the fills analysed and for both codes, the average ratio is very close to 1 in both IP1 and IP2, with a small standard deviation of a few percent, which we consider an excellent agreement.

At IP8, a the spread is much larger, with a standard deviation of 0.12. The detailed reason for this is not fully clear, although it is likely related to very large uncertainties in the LHCb luminosity calibration [69]. MBS consistently overestimates \(\int \mathcal {L}\,dt\) at IP8. In CTE, it is underestimated for a couple of fills in the early part of the run where very few bunches collided at IP8, while the luminosity in later fills with more colliding bunches is overestimated. This is possibly due to larger differences between individual bunches for the schemes with few collisions, that can only be seen with the MBS approach. The two codes are roughly consistent for all fills with a significant number of collisions at IP8.

In spite of the small discrepancies at IP8, we consider the overall agreement very good. The fact that both codes largely agree, in spite of being based on different underlying principles, strengthens our confidence in the results and shows that they can be used for reliable predictions of the luminosity in future machine configurations.

Table 5 Mean and RMS value of the ratio between simulated and measured integrated luminosity over the 30 analyzed fills in the LHC 2018 Pb–Pb run
Full size table

6 Projected performance in future runs

In this section, we estimate the heavy-ion performance in future LHC runs, based on simulations of typical fills using CTE and MBS, which are extrapolated to a typical one-month run. We discuss Pb–Pb operation in Sect. 6.1 and p–Pb in Sect. 6.3.

6.1 Pb–Pb performance in the HL-LHC baseline configuration

We assume the future running scenario detailed in Ref. [13], with the key parameters shown in Table 1. Some of the most notable differences to the 2018 configuration are the increased number of bunches and collisions, the higher beam energy, and the higher luminosity levelling target at IP2 of . We assume the same levelling at IP1 and IP5 for simplicity, although the experiments there could potentially accept a higher luminosity. At IP8, we assume levelling at to avoid quenches from BFPP. We use \(\sigma _{\mathrm {tot}}\) = 515 b at , with details given in Table 3. As given in Table 4, the processes influencing the emittances are stronger in the future configuration, with an effective initial risetime of about 8 h if core depletion is included.

Fig. 5
figure 5

The simulated HL-LHC Pb–Pb performance from CTE in terms of instantaneous luminosity (top) and integrated luminosity (middle) during a typical fill for the considered filling schemes from Ref. [13], shown together with the evolution of the beam intensity and normalised emittances (bottom). Only B1 is shown, but B2 is fully symmetric

Full size image

We show the simulated luminosity, intensity, and emittance evolutions for a typical future Pb–Pb fill in Fig. 5, with various proposed 50-ns filling schemes from Ref. [13]. The filling scheme names are defined with the customary LHC convention indicating the total number of bunches (b) followed by the number of collisions at ATLAS/CMS (IP1/5), ALICE (IP2), and LHCb (IP8). All new schemes, motivated by a recent request by LHCb for significant Pb–Pb data taking [27], use 1240 bunches but with different distributions of collisions between IP1, IP2, and IP5 on the one hand and IP8 on the other hand. Note that IP8 is slightly displaced from the symmetry point of the LHC ring and therefore requires displaced 50-ns bunch trains to have collisions at the expense of fewer collisions at the other IPs.

In addition to the 50-ns schemes with 1240 bunches, we simulate also the 75-ns filling scheme from 2018, which is a backup in case of production problems of the 50-ns beams. For this scheme, we assume the 2018 starting conditions in Table 1, except that we use a energy.

Fig. 6
figure 6

The time-averaged luminosity as a function of the turnaround time, calculated with Eq. (5) and the CTE simulations in Fig. 5 for the filling scheme

Full size image

As shown in Fig. 5, the intensities and emittances are similar between different 50-ns filling schemes, while different initial parameters are assumed for the 75-ns scheme. The levelling time at IP1 and IP5 is up to about 1 h. IP2 has a slightly longer levelling time because of its smaller crossing angle. At IP8, levelling is only needed for the three filling schemes with most collisions there. The IP8 luminosity is lower than at the other IPs because of its lower levelling target, fewer colliding bunches, the larger \(\beta ^{*}\), and larger crossing angle. The integrated luminosity improves significantly at IP8 for filling schemes that provide it with more collisions, at the price of a moderate reduction for the other experiments. As expected, the 75-ns scheme performs significantly worse than the 50-ns schemes.

As usual, the optimal fill length (time spent in collision), \(T_\mathrm {f}\), can be calculated to maximise the average luminosity \(\mathcal {L}_\mathrm {avg}\)

$$\begin{aligned} \mathcal {L}_\mathrm {avg}(T_\mathrm {f})\,=\, \frac{\int _{0}^{T_\mathrm {f}} \mathcal {L}(t)\, {\mathrm{d}}t}{T_\mathrm {f}+T_\mathrm {ta}}, \end{aligned}$$
(5)

where \(\mathcal {L}(t)\) is the instantaneous luminosity given by the simulation and \(T_\mathrm {ta}\) is the turnaround time, i.e. the time between the dump and the start of the collisions in the next fill. In reality, the turn-around time is the sum of a minimum irreducible value and a random value whose distribution reflects the general operational efficiency but cannot generally be predicted before the previous fill is dumped. However, for the sake of simplicity and consistency with the corresponding treatment for protons [30, 70], we assume a typical value of \(T_\mathrm {ta}\,=\,200\) min, based on the detailed time estimates in Ref. [16].

As an example, \(\mathcal {L}_\mathrm {avg}(T_\mathrm {f})\) is shown in Fig. 6 for the scheme. It turns out that the optimum fill time, \(T_\mathrm {f,opt}\), does not differ much between experiments. For the 50-ns schemes, it is around 3.8–3.9 h at IP1, 4.0–4.1 h at IP2, and 3.6–4.5 h at IP8. The spread in \(T_\mathrm {f,opt}\) is largest at IP8, but the curve of \(\mathcal {L}_\mathrm {avg}\) is also very flat (Fig. 6). Therefore, for Pb–Pb, we can adopt the value of \(T_\mathrm {f,opt}\) calculated for IP2 for each filling scheme.

The integrated luminosity \(\mathcal {L}_\mathrm {tot}\) in one Pb–Pb run is then estimated as

$$\begin{aligned} \mathcal {L}_\mathrm {tot}\, = \, \mathcal {L}_\mathrm {avg}(T_\mathrm {f,opt}) \times T_\mathrm {run} \times \eta , \end{aligned}$$
(6)

where \(T_\mathrm {run}\) is the total time allocated to the physics run. Typically, one month per year is allocated for heavy-ion operation. Allowing the first week for commissioning, we assume \(T_\mathrm {run}=24\) days.

The factor \(\eta \) in Eq. (5) is the operational efficiency, which should account for downtime and unavailability of the machine, premature fill aborts, occasional longer \(T_\mathrm {ta}\), and, most importantly, the build-up of performance to the ideal during the few weeks of the run. Conventionally, and conservatively, we take \(\eta \,=\,0.5\), as for HL-LHC proton operation [70]. In the 2018 heavy-ion run, a higher \(\eta \) was achieved, when the machine availability was exceptionally 85% after the initial commissioning [14]. Since \(\eta \) takes account of the build-up of luminosity that occurs during these short runs, it is less than the machine availability. Note also that \(\eta \,=\,0.5\) has been typical during some proton runs.

The calculated \(\mathcal {L}_\mathrm {tot}\) per one-month run using these parameters and assumptions is shown in Table 6 for the different filling schemes in Ref. [13]. Results without brackets are based on CTE, but corresponding MBS simulations have also been performed and these results are shown in square brackets. The results based on the two simulation codes, using fundamentally different models, agree within 5% except for the filling schemes with very few collisions at IP8. This strengthens our confidence in the results.

For the 50-ns schemes, \(\mathcal {L}_\mathrm {tot}\) varies in the range 2.2–2.6 nb\(^{-1}\) at IP1 and IP5, and in the range of 2.4–2.8 nb\(^{-1}\) at IP2, which has a smaller net crossing angle. With the backup 75-ns scheme (last row in Table 6), the loss in \(\mathcal {L}_\mathrm {tot}\) at these experiments is about 20–30% per run.

At LHCb, \(\mathcal {L}_\mathrm {tot}\) depends strongly on the filling scheme. With the highest number of collisions considered for 50-ns beams, about 0.5 nb\(^{-1}\) can be collected per run. This is significantly higher than the 0.35 nb\(^{-1}\) for the 75-ns backup scheme . At the same time, \(\mathcal {L}_\mathrm {tot}\) at the other experiments is also significantly higher with 50 ns. Therefore, we conclude that a 50-ns scheme should always be preferred if available, provided the projected beam quality can be achieved.

Table 6 Integrated luminosity (given in nb\(^{-1}\)) during a one-month Pb–Pb run at each experiment for the considered filling schemes from Ref. [13], assuming an operational efficiency of \(\eta \) = 0.5 h and 24 days available for physics operation in Eqs. (5)–(6)
Full size table

6.2 Sensitivity study and performance enhancements

So far, we have considered the baseline scenario, but using the simulations and Eqs. (5)–(6) we have also investigated the sensitivity of \(\mathcal {L}_\mathrm {tot}\) to changes in the baseline parameters and assumptions in the simulations. One representative filling scheme was selected, , and the inputs were varied in CTE. First, we have checked the influence of the IBS model. Figure 7 shows the evolution of the horizontal emittance for each of the built-in IBS models. Minor variations are observed in the emittance. The luminosity evolution remains very close to Fig. 5 in all cases due to the strong burn-off regime and is therefore not shown. The calculated \(\mathcal {L}_\mathrm {tot}\) per one-month run changes by up to 1% between these runs, where the largest difference is seen for the Piwinski smooth model, which is expected to have a poorer accuracy as it does not account for the full lattice. Further studies show that the IBS mixing factor between the transverse planes has an even smaller impact on \(\mathcal {L}_\mathrm {tot}\).

Fig. 7
figure 7

The evolution of the horizontal emittance during a typical fill as simulated with CTE using different IBS models. The beam parameters for future runs in Table 1 and the filling scheme were assumed

Full size image
Fig. 8
figure 8

The evolution of the instantaneous luminosity at IP2 when varying the bunch intensity (left), the transverse emittance (middle), and the non-collisional lifetime \(T_\mathrm {nc}\) during a typical fill as simulated with CTE. The beam parameters for future runs in Table 1 and the filling scheme were assumed

Full size image
Fig. 9
figure 9

The change in \(\mathcal {L}_\mathrm {tot}\) during a full one-month run as a function of the change in bunch intensity (left) and transverse emittance (right) as calculated from CTE simulations. The beam parameters for future runs in Table 1 and the filling scheme were assumed

Full size image

We have also studied the influence of the non-collisional lifetime \(T_\mathrm {nc}\), the bunch intensity N, and the transverse emittance \(\epsilon _{xy}\). So far, we used \(T_\mathrm {nc}\,=\,100\) h from the fit of the non-colliding bunches in the 2018 run. The source of these losses is not fully understood, and a scan was performed over \(T_\mathrm {nc}\) between 20 h and 220 h. The bunch intensity and emittance were varied within a \(\pm \,20\)% interval of the design values in Table 1.

The simulated \(\mathcal {L}\) at IP2 is shown in Fig. 8, and the resulting change in \(\mathcal {L}_\mathrm {tot}\) during a one-month run is shown in Fig. 9. Since the beam lifetime from the luminosity burn-off evolves in the range of 5–8 h, the negative influence of a lower \(T_\mathrm {nc}\) becomes clearly visible only for very low values, and for \(T_\mathrm {nc}>100\) h, the gains of longer lifetimes are very small. The estimated \(\mathcal {L}_\mathrm {tot}\) changes by less than 2% for all tested values of \(T_\mathrm {nc}\) except \(T_\mathrm {nc}\,=\,40\) h, giving a 4% reduction, and \(T_\mathrm {nc}\,=\,20\) h, which reduces \(\mathcal {L}_\mathrm {tot}\) by 10%. We consider that such low values of \(T_\mathrm {nc}\) are very unlikely and would only occur in case of serious problems with the machine.

On the other hand, Figs. 8–9 show that N strongly affects both \(\mathcal {L}\) and \(\mathcal {L}_\mathrm {tot}\), although \(\mathcal {L}_\mathrm {tot}\) is closer to a linear N dependence than the \(N^2\) dependence of the instantaneous \(\mathcal {L}\). A 10% change in N results in a 10–14% change in \(\mathcal {L}_\mathrm {tot}\). As expected, the dependence on \(\epsilon _{xy}\) is weaker, and the variations of \(\mathcal {L}_\mathrm {tot}\) stay within about 5% except at IP8, where the effect is slightly larger.

Finally, as shown in Eq. (6), \(\mathcal {L}_\mathrm {tot}\) varies linearly with the operational efficiency \(\eta \), which has a rather large uncertainty. Gains in \(\eta \) could have a very significant impact.

We have also investigated changes in the machine configuration whose feasibility is yet to be demonstrated. Nevertheless, these studies can indicate paths to higher luminosity. First, \(\beta ^{*}\) might be squeezed to smaller values. The lower limit, still to be quantified in detail, is given by the protection of the available aperture, as shown for proton operation in Refs. [71, 72].

An optimistic extrapolation of aperture measurements with protons [73] leads us to study with CTE a configuration with \(\beta ^{*}= 0.4\,\text{ m }\) at IP1, IP2, and IP5, and \(\beta ^{*}= 0.5\,\text{ m }\) at IP8, while keeping the baseline crossing angles. This relies also on the assumption that there is some margin to decrease the normalised beam–beam separation, which in the baseline is assumed to be larger than in the 2018 proton operation. The CTE results are illustrated in Fig. 10, which shows \(\mathcal {L}_\mathrm {tot}\) per one-month run as a function of the number of bunches colliding per IP, together with baseline performance in Table 6. In the configuration with smaller \(\beta ^{*}\), we observe a gain in \(\mathcal {L}_\mathrm {tot}\) of 6–9% at IP1, IP2, and IP5, while the gain at IP8 is up to a factor 2, depending on filling scheme. The gain is much smaller with more collisions at IP8, since the levelling time approaches \(T_\mathrm {f,opt}\).

Fig. 10
figure 10

Calculated integrated luminosity during a one-month Pb–Pb run at each experiment, as a function of the number of colliding bunches, comparing the baseline configuration (filled markers) with a reduced-\(\beta ^*\) configuration (open markers). The luminosity has been calculated using CTE simulations and Eqs. (5)–(6), assuming an operational efficiency of \(\eta \) = 0.5 and 24 days available for physics operation

Full size image

We have also investigated other performance enhancements for comparison. If the IP1/5 net half-crossing angle could be reduced to as at IP2, a gain of 6–7% at IP1/5 and a loss of 2–4% at IP2 are predicted. The crossing angle has a larger effect at IP8, where a net half angle of would bring a gain in \(\mathcal {L}_\mathrm {tot}\) of around 13%. It could also be envisaged not to use luminosity levelling at IP1/5, since the detectors are not limited and there is still some margin to the BFPP quench limit. However, using the above assumptions, the simulated gain in \(\mathcal {L}_\mathrm {tot}\) is limited at about 1–2%, with a similar loss at IP2/8. In this case, \(T_\mathrm {f,opt}\) is significantly different between IP1/5 and IP2. Shorter fills with a \(T_\mathrm {f,opt}\) calculated for IP1/5 could be compensated with a filling scheme that redistributes some collisions to IP2. Such a scheme could potentially provide gains for all experiments, but a detailed study is left as future work.

All these paths to increased \(\mathcal {L}_\mathrm {tot}\) still need to be demonstrated through further studies.

6.3 p–Pb performance in the HL-LHC baseline configuration

Fig. 11
figure 11

The simulated HL-LHC p–Pb performance from MBS in terms of instantaneous luminosity (first row) and integrated luminosity (second row) during a typical fill for the considered filling schemes from Ref. [13], shown together with the evolution of the intensity and normalised emittances of the Pb ions in B1 (third row) and protons in B2 (fourth row)

Full size image

For p–Pb operation, we base our running scenario on the detailed considerations in Refs. [13, 54]. We assume the same filling patterns as for Pb–Pb. This is an approximation, since the injection sequence for protons is different to that for Pb. While proton injection is expected to be rather flexible we account for the possibility that it cannot perfectly match the Pb filling pattern, by subtracting 5% of the total integrated luminosity as in Ref. [54].

In the simulations, B1 is assumed to contain Pb and B2 protons, but interchanging the beams does not influence the outcome at this level. We assume the same future Pb beam parameters as for Pb–Pb (Table 1) and that the proton beam has intensities chosen to be at \(3\!\times \! 10^{10}\) protons per bunch and normalised emittances around .

We include in MBS a phenomenological 100-h non-collisional lifetime for the Pb beam as in Sect. 6.1. In addition, as explained in [54] and based on fits to 2016 p–Pb data, an IP-dependent lifetime of 38 h is included per collision for bunches colliding in IP1 and IP5, 48 h in IP2 and 317 h in IP8. For the proton beam, a 5842-h non-collisional lifetime is used for all bunches, as observed in 2016 [54]. It is assumed that IP2 is levelled at \(\mathcal {L}\,=\,\) , a factor 5 higher than in 2016, thanks to the ALICE upgrade [35, 74]. As in 2016, IP1, IP5, and IP8 are not levelled. The assumed burn-off cross sections for p–Pb are shown in Table 3. Because of their strong dependence on Z, the electromagnetic processes are almost negligible in comparison with Pb–Pb operation.

Table 7 Integrated luminosity (in nb\(^{-1}\)) during a one-month p–Pb run at each experiment for the considered filling schemes from Ref. [13], using Eqs. (5)–(6)
Full size table

The simulated time evolution of the luminosity and beam parameters from MBS are shown in Fig. 11. The proton intensity remains almost unchanged as the luminosity losses are small in comparison to the total intensity. The Pb beam loses about one order of magnitude of its intensity in 10 h.

The IP1/5 luminosity starts from a peak of \(\mathcal {L}\,=\,\) , about twice what was achieved in 2016, and then decays rapidly. IP2 needs to be levelled significantly longer than with Pb–Pb, with typical levelling times of 6–7 h. The optimal fill time, calculated to maximise \(\mathcal {L}_\mathrm {avg}\) in Eq. (5) and assuming \(T_\mathrm {ta}\,=\,200\) min, is therefore rather different between IP2 (typically around 8 h) and IP1/5 (typically around 4.8 h). The fill time used for the calculation of \(\mathcal {L}_\mathrm {tot}\) for p–Pb is therefore estimated as the geometric mean of the two optima, typically .

The calculated \(\mathcal {L}_\mathrm {tot}\) per one-month run is given in Table 7 with the same assumptions as for Pb–Pb. Results from CTE and MBS agree to within 5% at IP1/5 and IP2. Slightly larger discrepancies of about 15% are found for the schemes with few collisions at IP8.

At IP1 and IP5, \(\mathcal {L}_\mathrm {tot}\) is in the range of 530–690 nb\(^{-1}\) per run for the 50-ns schemes, which is significantly higher than the 310 nb\(^{-1}\) obtained at IP2, mainly because only IP2 is levelled. The loss in \(\mathcal {L}_\mathrm {tot}\) with the 75-ns backup scenario is about 20–40% at IP1/5, but only about around 10% for IP2, thanks to the levelling.

At IP8, the maximum simulated \(\mathcal {L}_\mathrm {tot}\) in a one-month run over the filling schemes is about 150 nb\(^{-1}\). As for Pb–Pb, we could envisage improvements such as a smaller \(\beta ^{*}\) and smaller crossing angles if a higher \(\mathcal {L}_\mathrm {tot}\) will be needed.

The above estimates are for p–Pb operation as envisaged until the end of Run 4. Studies with MBS of the performance in collisions between protons and lighter nuclei, as proposed for after Run 4 [27], are presented in Ref. [75].

7 Conclusions

We have presented two independent simulation codes for beam and luminosity evolution in heavy-ion colliders, CTE and MBS, which are based on different underlying principles. CTE tracks bunches of macro-particles, subject to various physical effects, while MBS numerically solves a system of ordinary differential equations. MBS has the advantage of individually simulating every single bunch in the machine, while CTE can instead track non-Gaussian distributions and has more physical effects included, such as emittance blowup from collisions and machine aperture.

A thorough benchmark with data from 30 physics fills in the 2018 LHC Pb–Pb run shows an excellent agreement between both codes and measurements of the luminosity and intensity evolution. Given the starting conditions and fill duration, the integrated luminosity can typically be predicted within a few percent at the main experiments. Minor discrepancies are observed in emittances and bunch length, where the latter could be corrected in CTE by introducing a suitable coupling between the horizontal and vertical emittances.

The codes have then been used to predict the performance in future LHC heavy-ion runs with updated beam parameters, filling schemes, and operational assumptions. We estimate an integrated luminosity per one-month Pb–Pb run of about 2.2–2.6 nb\(^{-1}\) in ATLAS (IP1) and CMS (IP5), 2.4–2.8 nb\(^{-1}\) in ALICE (IP2) and up to about 0.5 nb\(^{-1}\) in LHCb (IP8), assuming a 50% operational efficiency and 24 days of physics operation after the initial commissioning.

With p–Pb, the integrated luminosity per one-month run is estimated to lie in the range of 530–690 nb\(^{-1}\) at ATLAS and CMS, at about 310 nb\(^{-1}\) at ALICE, and up to about 150 nb\(^{-1}\) at LHCb.

A sensitivity study shows that the results are rather robust against uncertainties in the IBS model and changes in the non-collisional lifetime. On the other hand, they are rather sensitive to any change in bunch intensity and operational efficiency and, to a lesser extent, to variations of the initial emittances. The present projections might thus be exceeded, e.g. if the operational efficiency or injected bunch intensity is higher. Other ways to gain integrated luminosity could include smaller \(\beta ^{*}\), smaller crossing angles, or modified filling schemes.