Impact of Beam-Beam Effects on Absolute Luminosity Calibrations at the CERN Large Hadron Collider

At the Large Hadron Collider (LHC), absolute luminosity calibrations obtained by the van der Meer (vdM) method are affected by the mutual electromagnetic interaction of the two beams. The colliding bunches experience relative orbit shifts, as well as optical distortions akin to the dynamic-$\beta$ effect, that both depend on the transverse beam separation and must therefore be corrected for when deriving the absolute luminosity scale. In the vdM regime, the beam-beam parameter is small enough that the orbit shift can be calculated analytically. The dynamic-$\beta$ corrections to the luminometer calibrations, however, had until the end of Run 2 been estimated in the linear approximation only. In this report, the influence of beam-beam effects on the vdM-based luminosity scale is quantified, together with the associated systematic uncertainties, by means of simulations that fully take into account the non-linearity of the beam-beam force, as well as the resulting non-Gaussian distortions of the transverse beam distributions. Two independent multiparticle simulations, one limited to the weak-strong approximation and one that models strong-strong effects in a self-consistent manner, are found in excellent agreement; both predict a percent-level shift of the absolute pp-luminosity values with respect to those assumed until recently in the physics publications of the LHC experiments. These results also provide guidance regarding further studies aimed at reducing the beam-beam-related systematic uncertainty on beam-beam corrections to absolute luminosity calibrations by the van der Meer method.


Introduction
The determination of the absolute scale of the luminosity delivered to the ALICE, ATLAS, CMS and LHCb experiments at the LHC relies almost entirely [1] on the van der Meer (vdM) method [2,3]. Luminosity calibrations, which must be performed under specially tailored beam conditions, use beam-separation scans, also known as vdM scans, to relate the collision rate measured by a given luminometer to the absolute luminosity inferred from directly measured beam parameters. The proportionality between the measured rate and the absolute luminosity is expressed in terms of a "visible" cross-section denoted σ vis , which is specific to the luminometer and the counting method considered. The precision requirements, which are driven by the LHC physics program, call for vdM-calibration uncertainties in the 1% (0.5%) range at the LHC (HL-LHC). The beam parameters used to determine the absolute luminosity during vdM scans are the particle populations and transverse sizes associated with each colliding-bunch pair. During a beam-separation scan, the mutual electromagnetic interaction between two opposing bunches induces separation-dependent orbit shifts, as well as variations in the transverse size and shape of each bunch, thereby distorting the luminosity-scan curves from which the beam-overlap integrals, or equivalently the convolved transverse beam sizes, are extracted. Depending on the beam conditions under which the scans are performed, the magnitude of the resulting calibration bias -if left uncorrected -represents a significant fraction of, or can even exceed, the σ vis systematic-uncertainty budget.
The methodology first developed in Refs. [4,5] to quantify the impact of beam-beam effects on the absolute luminosity scale has been adopted since 2013 by all LHC Collaborations. With as input the beam separation dialed-in at each scan step, plus the measured bunch currents and convolved beam sizes, the orbit shift induced by beam-beam deflections was calculated analytically [6], and the optical distortion associated with the dynamic-β effect was evaluated in the linear approximation [4] using the MAD-X package [7]. The two effects impacted the luminosity scale in opposite ways, resulting at the time in a net upwards correction to σ vis of 1-1.5% for vdM calibrations at √ s = 13 TeV.
A recent reevaluation of this methodology, using a new beam-beam simulation package specifically developed for vdM-scan studies, revealed that the linear approximation used in the MAD-X simulation results in a significant underestimate of the optical-distortion correction [8]. These findings motivated the studies reported in the present paper, which is organized as follows.
After a brief overview of the vdM-calibration methodology at the LHC and of the impact thereon of beam-beam effects (Sec. 2), the relevant simulation tools are presented in Sec. 3: the MAD-X package used in the original implementation [4,5]; the B*B package presented in Ref. [8], that models beam-beam dynamics in the transverse plane under the weak-strong approximation; and the long-established COMBI [9] multiparticle code, that is more CPU-intensive but can simulate beam-beam effects in the strong-strong regime with the optional inclusion of longitudinal dynamics. Following a systematic crossvalidation of the latter two packages, B*B is used in Sec. 4 to develop an easy-to-use parameterization of beam-beam corrections to vdM scans in the limit of round, initially Gaussian bunches of equal brightness that collide at a single interaction point (IP) with zero crossing angle. Deviations of the colliding-bunch configuration from this idealized limit: non-Gaussian tails, elliptical transverse bunch profiles, non-zero crossing angle, collisions at multiple IPs, or unequal-brightness beams, are then either accounted for using simulation-based adjustments to the idealized parameterization, or found to be small enough to be treated as a systematic uncertainty. These and other sources of systematic uncertainty, such as tune or β * settings, that may affect the beam-beam correction in the vdM regime are consolidated in Sec. 5. An overall summary and a terse outlook are offered in Sec. 6.

Luminosity-calibration methodology at the LHC
The vdM-scan formalism [2,3] that underpins the determination of the absolute luminosity scale at the CERN ISR, RHIC and the LHC is summarized, for the simplest case, in Sec. 2.1 below; generalizations of this formalism can be found in Refs. [8,10,11]. At the LHC, because of both instrumental and accelerator-physics reasons [1], vdM scans are not performed during normal physics operation, but rather under dedicated beam conditions (Sec. 2.2). Their two fundamental ingredients, the transverse beam separation and the measured collision rate, are both affected by the beam-beam interaction (Sec. 2.3), at a level that is significant on the scale of the precision goals outlined in Sec. 1.

Absolute luminosity scale from measured beam parameters
In terms of colliding-beam parameters, the bunch luminosity L b is given by where the opposing bunches are assumed to collide head-on and with zero crossing angle, f r is the LHC revolution frequency, n 1 n 2 is the bunch-population product andρ 1(2) (x, y) is the normalized particle density in the transverse (x-y) plane of beam 1 (2) at the IP. 1 With the standard assumption that the particle densities can be factorized into independent horizontal and vertical component distributions,ρ(x, y) = ρ x (x) ρ y (y), Eq. (1) can be rewritten as where Ω x (ρ x1 , ρ x2 ) = ρ x1 (x) ρ x2 (x) dx is the beam-overlap integral in the x direction (with an analogous definition in the y direction). In the method proposed by van der Meer [2] at the ISR (Fig. 1, top), the overlap integral (for example in the y direction) can be calculated as where R y (δ y ) is the collision rate, or equivalently the luminosity in arbitrary units, measured during a vertical scan at the time the two beams are separated vertically by the distance δ y . Because the collision rate R y (δ y ) is normalized to that at zero separation R y (0), any quantity proportional to the luminosity can be substituted in Eq. (3) in place of R.
Defining the vertical convolved bunch size Σ y [1] as and similarly for Σ x , the bunch luminosity in Eq. (2) can be rewritten as which allows the absolute bunch luminosity at zero separation to be determined from the revolution frequency f r , the bunch-population product n 1 n 2 , and the product Σ x Σ y which is measured directly during a pair of orthogonal scans. If the transverse density profile of each beam B (B = 1, 2) can be described by a single Gaussian of width σ iB (i = x, y), the convolved widths are given by where ϵ iB is the geometrical emittance of beam B in plane i, and β * iB is the corresponding value of the β function at the IP. 2 In such a case, the beam-separation dependence of the collision rate is given by The luminosity curve R i (δ i ) is also Gaussian, and Σ i coincides with the standard deviation of that distribution. It is important to note, however, that the vdM method does not rely on any particular functional form of R i (δ i ): the quantities Σ x and Σ y can be determined for any observed luminosity curve from Eq. (4) and used with Eq. (5) to determine the absolute luminosity at δ x = δ y = 0.
In the more general case where the factorization assumption breaks down, i.e. when the particle densities cannot be factorized into a product of uncorrelated x and y components, Eq. (2) no longer holds, and a single pair of horizontal and vertical scans is no longer sufficient to measure the overlap integral in Eq. (1). One must then generalize the formalism to the two-dimensional case [3], and scan over a grid in the (δ x , δ y ) beamseparation space to measure the product of the convolved bunch widths [1,8]: [Σ x Σ y ] = 1 2π R x,y (δ x , δ y ) dδ x dδ y R x,y (0, 0) .
Here the square brackets highlight the fact that in the presence of non-factorization, the quantity [Σ x Σ y ] can no longer be broken down into a product of two independent quantities. Equation (5), however, remains formally unaffected, as do Eqs. (9)-(10) below.
In terms of luminometer observables, the bunch luminosity can be written as where µ vis is the average number of inelastic collisions per bunch crossing detected by the luminometer considered, and σ vis is the associated visible cross-section. Since µ vis is a directly measurable quantity, the calibration of the absolute luminosity scale amounts to determining the visible cross-section σ vis . Equating the absolute luminosity computed from beam parameters using Eq. (5) to that measured according to Eq. (9), yields: where µ vis,pk is the visible interaction rate per bunch crossing reported at the peak of the scan curve ( Fig. 1, bottom). Equation (10) provides a direct calibration of the visible cross-section σ vis in terms of the peak visible interaction rate µ vis,pk , the product of the convolved bunch widths Σ x Σ y , and the bunch-population product n 1 n 2 .
(11) Here θ c is the full crossing angle, σ zB (B = 1, 2) are the RMS bunch lengths of beams 1 and 2, and σ cB the transverse single-beam sizes in the crossing plane. The peak luminosity is reduced by the same factor. The corresponding increase in the measured value of Σ x or Σ y is exactly compensated by the decrease in µ vis,pk , so that Eqs. (4)-(10) remain valid, and no correction for the crossing angle is needed in the determination of σ vis .

Beam conditions during van der Meer scans
The strength of the beam-beam interaction is traditionally quantified by the linear beam-beam parameter, defined as [15,16]: 2π A ion,2 γ 2 σ x1 (σ x1 + σ y1 ) (12) Here ξ x2 is the horizontal beam-beam parameter experienced by beam 2 (B2), the "witness beam"; n 1 is the bunch population of beam 1 (B1), the "source beam"; r 0 = e 2 /4πϵ 0 m p c 2 is the classical radius of the proton; A ion,B and Z B (B = 1, 2) are the atomic mass number and charge number of the beam-B particle type (proton or fully stripped ion); β * x2 is the value of the B2 horizontal β function at the IP; γ 2 is the relativistic factor of the B2 particles, and σ x1 (σ y1 ) is the horizontal (vertical) transverse RMS size of B1. Formulas for the other beam and the other plane are obtained by interchanging B1 and B2, and/or x and y.
For the most frequent case of pp collisions, Equation (12) takes the more familiar form: In the case of equally populated, equally sized round beams, this expression becomes much simpler: where n is the bunch population, and σ 0 = √ ϵ β * is the nominal RMS beam size at the IP.
In 2018, during high-luminosity physics running in proton-proton (pp) mode, the LHC collided up to 2544 bunches with typical initial intensities of 1.1 × 10 11 p/bunch, grouped in trains of 36 to 144 bunches with a minimum interbunch spacing of 25 ns. Beams crossed with a half-angle θ c /2 of ±130 µrad in order to mitigate the impact of the long-range beam-beam interaction at parasitic crossings. At the start of stable beams, the emittance was typically 2 µm·rad [17], the singlebunch luminosity around 7.8 × 10 30 cm −2 s −1 at IP1 and IP5, and the total luminosity close to 1.9×10 34 cm −2 s −1 at each of these two IPs. These values correspond to a pile-up parameter µ of around 55 inelastic pp collisions per bunch crossing, and to a bunch-averaged, head-on beam-beam parameter ⟨ξ⟩ of approximately 0.005. The brightness, however, varied significantly along the bunch string, occasionally resulting in ξ values as high as 0.007 for some of the bunches.
In contrast, during pp vdM scans (Table 1), the injected emittance is deliberately blown up and the bunch population significantly lowered, in order to reduce the impact of beam-beam effects as well as minimize the unbunched-beam fraction and the intensity of satellite bunches. The bunches are isolated rather than in trains, and their number is limited to 152 at most, in order to eliminate parasitic crossings and collide with zero crossing angle in the interaction regions where the beam-line layout so permits. The β function at the IP is increased such as to bring the pile-up parameter µ down to around 0.5, i.e. in a regime where luminometers are free of instrumental nonlinearities; this carries the additional advantage that it significantly increases the transverse luminous size, allowing a more precise measurement of its beam-separation dependence [1,13].

Beam-beam-induced biases and their correction
The mutual electromagnetic interaction between colliding bunches shifts their orbits, and therefore modifies their transverse separation at the IP (Sec. 2.3.1); it also distorts their transverse density distributions (Sec. 2.3.2). These two effects depend on the nominal separation ∆ dialed-in at each scan step. Their combination impacts both the normalized integrals (Σ x , Σ y ) and the peak (µ vis,pk ) of the luminosity-scan curves used in determining the absolute luminosity scale. The strategy for correcting the resulting biases is outlined in Sec. 2.3.3; its detailed implementation is developed in later chapters, in particular in Secs. 4.2.3 and 4.6.5.

Orbit shift
When two positively charged bunches collide with a non-zero impact parameter, they experience a mutually repulsive angular kick equivalent to that of a dipole located at the collision point, the strength of which depends on the beam separation. In the round-beam limit and for pp collisions, the angular kick experienced by a B2 bunch during a horizontal beam-separation scan is given by [18]: and similarly for a vertical scan. Here δ x (resp. δ) is the horizontal (resp. total) beam separation, and Σ R = Σ x = Σ y is the transverse convolved beam size. This formalism has been extended by Bassetti and Erskine [6] and by Ziemann [19] to the case of elliptical beams. This angular deflection, first observed at the Stanford Linear Collider in e + e − collisions [18], can be measured using beam-position monitors (BPMs) installed both upstream and downstream of the IP, as illustrated in Fig. 2 for the LHC [20].
In a circular collider, the beam-beam angular kick experienced by each beam B (B = 1, 2) in the i-plane (i = x, y) results in a shift of its position at the IP, given by [21] Fig. 2: Total beam-beam deflection angle (B1-B2) as a function of the nominal separation ∆ (denoted here by "knob setting"), during a vertical vdM scan at the ATLAS IP in pp collisions at √ s = 8 TeV (LHC fill 3316). For each beam separately, the deflection angle is obtained from the difference between the outgoing-and incomingbeam angles measured by BPMs located in the LHC arcs outside the closed-orbit bump used for the scans. The zero of the horizontal axis is arbitrary, since only relative beam displacements matter. The vertical convolved beam size, denoted here by Σ, is extracted from a fit to Eq. (14), shown by the red curve, with as input the bunchaveraged population N (Figure reproduced from Ref. [20]  where β * is the value of the β function at the IP and Q i is the betatron tune. The actual beam separation at the IP δ i therefore differs slightly from the nominal separation ∆ i : where the beam-beam-induced change in beam separation, hereafter denoted by "orbit shift", is given by δ bb i = δ bb i1 − δ bb i2 . Since δ bb iB is, to first order, proportional to the corresponding beambeam parameter ξ iB (Eqs. (12) and (14)), the orbit shift varies from one bunch to te next. The ∆-dependence of δ bb i mirrors that displayed in Fig. 2 for the deflection angle, with a peak-to-peak swing of ±(1-2) µm under typical vdM-scan conditions, to be compared to typical Σ i values in the 110-160 µm range. The value of δ bb i changes from scan step to scan step, thereby expanding in a non-linear fashion the beam-separation scale, i.e. the horizontal axis of scan curves such as that illustrated in Fig. 1b. As a result, the overlap integral of Eq. (4) increases by typically 0.7-1.4% per plane, corresponding to a positive correction to σ vis in the 1.4-2.8% range.
In practice, the correction for the beam-beam orbit shift is implemented as follows. At each scan step and for each colliding-bunch pair separately, δ bb i is calculated using the Bassetti-Erskine formula [6], with as input the measured bunch populations (n 1 , n 2 ) and uncorrected convolved bunch widths (Σ x , Σ y ), as well as the beam energy, the β * setting and the tunes. The separationdependent collision rate R i (δ i ) is then integrated according to Eq. (4) to obtain the beam-beam corrected values Σ c,Orb x and Σ c,Orb y , using the beambeam corrected separation δ (Eq. (16)) instead of the nominal separation ∆. The agreement of this simple analytical procedure with the predictions of self-consistent multi-particle simulations will be addressed in Sec. 3.5.

Optical distortions
Not only does the electromagnetic field of the B1 bunch deflect the B2 bunch as a whole: it also acts as a non-linear lens that perturbs the trajectory of the individual particles in that bunch, thereby modifying the transverse density distribution of both bunches in a separation-dependent manner.
For small-amplitude particles and beams in head-on collision (Fig. 3, short blue arrows), the force is rather linear and resembles that of a quadrupole (green line), resulting in a tune shift proportional to the beam-beam parameter ξ and in the subsequent dynamic-β effect [4]. This "quadrupole strength" is proportional to the derivative of the beam-beam force; it is largest, and repulsive, for small-amplitude particles, changes sign around 1.6σ, and becomes weakly attractive at larger amplitude (dotted brown line).
If for simplicity one assumes that for a given beam separation, all particles are subject to the same quadrupolar-like force (the strength and sign of which depend on the beam separation), then the value of β * at the scanning IP is modulated by the linear component of the beam-beam force. This results in a modulation of the transverse beam size, and therefore in a beam-separationdependent modulation of the actual luminosity; however the actual shapes of the transverse density distributions projected on the x and y axes Fig. 3: Beam-beam force exerted by the B1 (B2) source bunch as a whole (red curve), as a function of the betatron amplitude of a B2 (B1) test particle, for equally sized round beams. The amplitude is in units of the RMS beam size. The green and brown lines correspond to the linear component of the force at amplitudes of zero and 4σ respectively, and are akin to the effect of a quadrupole. The short (long) double arrows illustrate the range of transverse kicks experienced by a small (large) amplitude particle for beams either in head-on collision (solid blue lines), or transversely separated by 4σ (dashed magenta lines). remain unaffected by the quadrupolar-like force. This is the approximation that was adopted in the first implementation of the optical-distortion correction [4,5], and that will be further discussed in Sec. 3.1.
While for small amplitudes (short arrows) the force remains approximately linear, at amplitudes larger than 1σ (long arrows) it includes significant non-linear contributions. Large-amplitude particles, therefore, experience a tune shift and a β-beating that depend both on the particle amplitude (short vs. long arrows) [22], and on the beam separation (blue vs. magenta arrows). The resulting optical distortions include not only a change in optical magnification as in Ref. [4], but also distortions of the shape of the transverse density distributions. Describing their beam-separation dependence requires numerical simulations, that are detailed in Sec. 3.
In practice, the correction for optical distortions is implemented as follows. Separately for each scan step in a horizontal and vertical vdMscan pair, the luminosity-bias factor [L/L 0 ] Opt associated with beam-beam-induced optical distortions is extracted, as a function of the nominal separation ∆ i , from one of the multiparticle simulations described in Sec. 3. Here L refers to the luminosity that would be measured in the presence of beam-beam optical-distortion effects, and L 0 is the corresponding luminosity if beam-beam effects were turned off altogether, all other conditions remaining unchanged. The quantity L 0 is dubbed the nominal luminosity; it is akin to "Monte Carlo truth", and is accessible in the simulation only. The beam-beam corrected collision rate R c i (δ i ) is then computed by dividing the measured collision rate by this simulation-based luminosity-bias factor: and used instead of R i (δ i ) in computing the beambeam corrected convolved bunch sizes Σ c,Opt x and Σ c,Opt y (Eq. (4)), the peak rate µ c = µ vis,pk R c i (0)/R i (0), and from these quantities the visible cross-section σ c,Opt vis (Eq. (10)).

Beam-beam correction strategy
Conceptually, the principle of the beam-beam correction to vdM calibrations is to determine the visible cross-section σ c vis from the convolved bunch sizes and peak collision rates corrected both for the orbit shift (Sec. 2.3.1) and for optical distortions (Sec. 2.3.2), i.e. corrected to the values (Σ c x , Σ c y , µ c ) that these observables would take if the beam-beam interaction could be turned off during the scan. The quantity σ c vis is then the proportionality constant that translates a measured visible interaction rate µ vis into the corresponding bunch luminosity L b . It is important to note that even though the actual luminosity is always modified by the beam-beam interaction, including during head-on collisions typical of routine physics running, beam-beam corrections are only needed during scans, basically because the beamseparation dependence of the beam-beam effects distorts the vdM-scan curves. Once σ c vis has been determined as specified above, it can always be used to translate the measured collision rate into luminosity units, irrespective of the extent to which this collision rate has been enhanced by the beam-beam interaction.
In practice, orbit-shift and optical-distortion corrections must be applied on a bunch-by-bunch basis, with as inputs the measured bunch populations (n 1 , n 2 ) and uncorrected convolved bunch widths (Σ x , Σ y ), as well as the beam energy, the β * setting and the tunes. In many cases, these corrections can be extracted from a one-time parameterization of the simulation results in terms of the bunch-specific beam-beam parameter value and of the nominal tunes Q x , Q y . This procedure avoids the repeated use of CPU-intensive, timeconsuming multiparticle simulations; it is detailed in Sec. 4.

Linear approximation with MAD-X
Since the linear part of the beam-beam force is similar to a quadrupolar field (Fig. 3), one expects the beam-beam interaction to contribute additional focusing or defocusing, thereby shifting the tunes and affecting the optical functions all around the ring, including at the IP itself: this is known as the "dynamic-β" effect. In the specialized case of equally sized round beams colliding head-on at a single location, the resulting change in the β function at the IP is given by [4]: where β * 0 is the value of the unperturbed β function at the IP, β * its value in the presence of the beam-beam interaction and Q the tune. The beam-beam parameter ξ is proportional to the derivative of the beam-beam force. 3 Equation (18) implies that the beam-beam induced change in the IP β function, and therefore in the IP beamsize squared and in the luminosity, depends only on ξ and on Q. In addition, this dynamicβ effect is, to first order, proportional to ξ. This is the physical motivation underlying the parameterized-correction approach developed in Sec. 4.
In Ref. [4], the general-purpose optics code MAD-X [7] was used to model the dynamic-β effect during simulated vdM scans. In this software package, beam-beam elements can be inserted at one or several IPs, and their impact on the tunes and on the single-particle optical functions computed as a function of the beam separation ∆. The procedure effectively assumes that for a given beam separation, all particles in the bunch experience the same beam-beam kick, equal to that applied to a zero-amplitude particle for that particular value of ∆. Unperturbed bunches are implicitly supposed to be strictly Gaussian, and to remain so in the presence of the beam-beam interaction: only the change in optical magnification between the LHC arcs and the IP is accounted for in this method.
This study was carried out for the reference parameter set listed in Table 1. The corrections were then adapted to the beam conditions of different vdM sessions, using the assumption that ∆β/β * 0 = β * /β * 0 − 1 scales linearly with the value of ξ inferred from the beam parameters measured during each scan. The simulated beam-separation dependence of the transverse beam size squared, i.e. the value of (L/L 0 ) −1 , is illustrated in Fig. 4. Under head-on conditions (∆ = 0), the intrinsically defocusing beam-beam force results in a ∼ 0.7% reduction of the beam size squared at the IP, i.e. to an increase of the luminosity L. This apparent contradiction results from the numerical value of the fractional tunes the LHC optics was designed for (Eq. (18)). As ∆ increases, the tune shift and the dynamic-β effect weaken, change sign (in the scanning plane only) around ∆/σ 0 ∼ 1.6, then peak and finally vanish asymptotically at very large separation. The change in optical magnification at ∆ = 0 is slightly different in the x and y planes, because the corresponding fractional tunes q x and q y differ by design; the difference in beam-separation dependence also reflects the fact that the vertical kick never changes sign during a horizontal scan, while the horizontal kick does.
Historically, the optical-distortion correction to σ vis predicted by MAD-X under typical Run-1 and Run-2 vdM conditions lay in the 0.2-0.4% range, much smaller than that associated with the orbit effect. For a long time, therefore, it was considered small enough that even if imperfect, it remained sufficiently accurate in view of the systematic uncertainty assigned at the time to the overall beam-beam correction, as documented e.g. in Refs. [5,13].
In hindsight, however, the limitations of applying the MAD-X approach to beam-beam corrections may not have been fully appreciated. Since at zero beam separation, the slope of the beam-beam force is steepest for zero-amplitude particles (Fig. 3, green line), and since in MAD-X all particles in the witness bunch are assumed to experience the same linearized force, the dynamicβ effect predicted by MAD-X at ∆ = 0 is likely to be an overestimate; this is confirmed analytically below. The situation is reversed at large beam separation: here the derivative of the force is  Fig. 4: Beam-separation dependence of the transverse RMS beam sizes during a simulated horizontal vdM scan. The unperturbed beams are assumed to be round and perfectly Gaussian; their parameters are listed in the right column of Table 1. The vertical axis is the ratio squared of the actual beam size σ to its unperturbed value σ 0 ; the horizontal axis is the nominal beam separation in units of σ 0 . The MAD-X calculation of the single-particle dynamic-β effect (blue curves) is compared to the results of the B*B (square markers) and COMBI (red curves) multiparticle codes. The dark-and light-colored curves and markers illustrate the evolution, during the horizontal scan, of the horizontal and the vertical beam size respectively.
larger at large amplitude (green line) than at small amplitude (brown line), suggesting that MAD-X underestimates the optical distortions at large beam separation. This conjecture, however, can only be confirmed using multiparticle simulations.

Analytical estimate of optical distortions at zero beam separation
The particle phase space at any location s along a storage ring is described by an ellipse, the shape of which depends on s in a manner described by the well known Courant-Snyder parameters α, β and γ [16]: Here u = x, y is the deviation of the singleparticle orbit from the reference trajectory, p u = du/ds, and J u is the action variable that represents the invariant of the motion for each single particle when the reference energy is not changing. The particle position in transverse phase space is fully described by the action J u and by the corresponding phase variable ϕ u defined as: In addition to distorting the closed orbit, the beam-beam interaction acts as a non-linear electromagnetic lens that includes a quadrupolar term; the latter affects the optical functions β and the betatron tunes of the individual particles [22,23]. During a vdM scan and for each colliding bunch, the quadrupolar component of the electromagnetic field of the opposing bunch changes as a function of the relative transverse separations δ i (i = x, y) between the beams centroids, and of the single-particle transverse actions J x and J y [23,24]. This results in a shift in the betatron tune, the so-called detuning with amplitude, given by [25]: Here ξ x is the beam-beam parameter, I 0 and I 1 are modified Bessel functions of the first kind, ϵ is the emittance (assumed equal in x and y), and t is a bound variable. This integral can only be solved when one action is zero, obtaining (for J y = 0): For a zero-amplitude particle (J x = J y = 0) and beams in head-on collision (δ x = δ y = 0), the change in tune is equal to the beam-beam parameter ξ x . Head-on beam-beam β-beating can be derived from the tune shift 4 in Eq. (20) as: 4 The overall minus sign in Eqs. (20) and (21) expresses the fact that when the charges of the colliding bunches have the where is the linear β-beating or, equivalently, the βbeating at J x = J y = 0.
In the zero-separation case, the β-beating averaged over the particle distribution (assumed Gaussian) can be calculated as: For bunches colliding with zero transverse separation, therefore, the particle-action distribution is modified such that the average beam-beam beating is reduced to about 63% of the single-particle estimate computed in Ref. [4]. The corresponding impact on the head-on luminosity can be derived analytically, as follows.
The RMS single-beam size σ x = √ β x ϵ, including the action-dependent, beam-beam-induced βbeating given by Eq. (21), is computed using the following relation: Assuming Gaussian particle-density distributions, this yields where, for simplicity, we used an emittance value of ϵ = 1. The triple integral has an exact solution: same sign, the tune shift is negative, even though the beambeam parameter remains positive by definition (Eq. (12)).
This calculation implies that in head-on collisions, beam-beam-induced linear β-beating of magnitude ∆β 0 /β in both the horizontal and the vertical plane, causes a relative luminosity change of For collisions with zero transverse separation, in other words, the beam-beam interaction changes the luminosity by only half of what is expected from linear β-beating: this is consistent with the predictions of respectively COMBI and MAD-X displayed in Fig. 4.
Using the reference parameter set in Table 1, for instance, the beam-beam-induced linear βbeating amounts to about -0.6%, leading to a head-on luminosity enhancement of about 0.3%. During high-luminosity physics running typical of LHC Run 2, the effect is computed to be two to three times larger; it is both ξ-and tune-dependent (see Eq. (18)), and therefore its magnitude changes as beam conditions evolve.
When beams are transversely separated, analytical computations become impossible, forcing one to resort to numerical simulations such as those described below.

Weak-strong limit: B*B
The B*B package [8] is a multiparticle-simulation program developed specifically for assessing beambeam biases in vdM scans, that was optimized for speed by adopting several simplifying assumptions. It aims at predicting the corresponding beam-beam corrections with better than 0.1% accuracy (for a given set of input parameters), so as not to contribute significantly to the overall vdM-calibration uncertainty. The code is written in C++; it can be used as a standalone application, or as a library available in C, C++, Python and R.
The initial transverse particle density distributionsρ B (x, y) (B = s, w) are modeled by a (linear combination of) two-dimensional elliptical Gaussian(s). Since during typical vdM scans at LHC, beam-beam-induced bunch-shape deformations remain small, they are taken as negligible when computing the electromagnetic field of the source bunch s, which therefore remains unperturbed as a function of the transverse beam separation. This makes it possible to precompute this field over a two-dimensional grid in the transverse plane at initialization time, and to use only fast interpolations on all subsequent machine turns. The opposing, "witness" bunch w is represented by a set of O(1000) macroparticles, that are transported around the ring using linear maps, with as input the nominal (x, y) phase advance between consecutive collision points. The B*B simulation, therefore, falls in the "weak-strong" category in that it models the transverse deformation of the density distribution of the witness bunch (ρ w → ρ w + δρ w ) caused by the electromagnetic field of an unperturbed source bunch, the density distribution of which remains unaffected. Effects such as coherent bunch oscillations, or the distortion of the shape (and therefore of the field) of the source bunch induced by the deformation of the witness bunch, are therefore implicitly neglected.
The macroparticles are selected from a twodimensional (x, y) grid of betatron amplitudes, and assigned weights that are precalculated from the initial Gaussian density distributionρ w . Their initial betatron phase is chosen randomly; the phase then samples the full [0, 2π] interval over the next O(100 − 1000) turns or so.
The overlap integral of the perturbed and unperturbed bunches (Eq. (1)) is computed as the sum over macroparticles, of the source-bunch densityρ s (x w , y w ) evaluated at the current location (x w , y w ) of the macroparticle considered, and multiplied by the weight of that same macroparticle. If the populations and initial transverse-density distributions of the two colliding bunches are identical (n 1 = n 2 andρ 1 (x, y) =ρ 2 (x, y)), the simulation needs to be run once only; in the presence of any initial B1-B2 asymmetry, it needs to be run twice, with the roles of the source and the witness bunch swapped between the two beams. The overall beam-beam bias affecting the overlap integral is approximated by (ρ 1 + δρ 1 ) (ρ 2 + δρ 2 ) dx dy − ρ 1ρ2 dx dy ≈ (δρ 1 ·ρ 2 +ρ 1 · δρ 2 ) dx dy, up to some physical constants and where the second-order term δρ 1 δρ 2 dx dy is neglected.
The uncertainties arising from the manner in which the beam-beam force is switched on in the calculation (gradually or instantaneously), from the granularity of the simulation (finite number of macroparticles, largest sampled transverse amplitude, random choice of the initial phases), and from other simulation-control parameters such as the number of accelerator turns, are discussed in detail in Ref. [8] and found to lie well below 0.1%.
Finally, even though the density distributionŝ ρ s,w are only two-dimensional, and therefore represent the projection of the full six-dimensional distribution onto a plane perpendicular to the beam axis, the B*B package is capable of simulating the geometrical effects that arise in the presence of a non-zero crossing angle in a plane of arbitrary orientation [8]. Simulating longitudinal dynamics, however, and in particular the potential impact of a finite crossing angle on beam-beam corrections to the vdM calibration, requires a fully six-dimensional treatment such as that outlined below.

Strong-strong model: COMBI
The COMBI code [9] has been developed over the years for simulating, in a self-consistent manner, the coherent beam-beam interaction between multiple bunches coupled by head-on and/or longrange beam-beam encounters [26,27,28]. It includes a first level of parallelization based on the Message Passing Interface (MPI) [29], and a second one sharing several CPUs per node using OpenMPI [30,31].
The code has been optimized to handle simultaneously multiple bunches at several interaction points, thereby allowing flexible collision patterns. The circumference of each accelerator ring is modeled by a number of equally spaced slots that define the possible bunch positions. At each location one can assign an action (e.g. head-on or long-range beam-beam interaction, linear or nonlinear magnetic element, Landau octupole, linear or non-linear map,...) that will be executed when a bunch is present. The actions corresponding to head-on or long-range beam-beam interactions require one bunch from each beam in order to be performed. Each macroparticle is tracked individually under the effect of the preassigned actions, in either four or six dimensions. In the LHC arcs, the macroparticles are transported by applying a linear transfer map to their coordinates, using phase advances precomputed by MADX with, as input, the nominal optical configuration used during vdM sessions.
In B*B, the transverse density distribution of the source bunch, and therefore its electromagnetic field, remain unchanged from turn to turn; only the witness bunch contains macroparticles, that are tracked over thousands of turns in the transverse plane. Coherent effects, therefore, cannot be modeled, and neither can longitudinal dynamics. COMBI, in contrast, describes the two partners in a colliding-bunch pair as independent sets of macroparticles. In this paper, the initial, unperturbed density distributions are uncoupled single Gaussians by default; however, arbitrary distributions, such as a linear combination of Gaussians, can be used instead.
The beam-beam interaction can be described by different models: a four-dimensional Gaussian lens [26,27], a six-dimensional Gaussian lens [32,33], or a field computed from the actual charge distributions by the HFMM method [34]. Multiple beam-beam encounters, and therefore the evolution of bunch parameters such as emittance or transverse barycenter position, as well as coherent beam-beam effects, are treated in a self-consistent manner: the particles trajectories affected by the beam-beam interaction modify the density distribution of the corresponding bunch, and the fields produced by both partners are updated turn by turn from these modified distributions. Whether these fields are estimated in the Gaussian approximation, or by the HFMM method, yields effectively identical results for the full range of beam-beam parameter values considered in this paper (see Appendix A).
At a given IP, the luminosity per collidingbunch pair can be computed either analytically, using the Gaussian formalism of Ref. [14], or by evaluating numerically the actual overlap integral of the two colliding distributions. The first method evaluates the luminosity from Eqs. (5)-(7) by substituting the single-beam sizes σ iB with the RMS transverse widths of the macroparticle distributions, and the separation δ with the distance between the barycenters of these distributions. In the second method, the overlap integral of the macroparticle distributions is computed using functionalities developed specifically for this purpose. This is the approach adopted throughout this paper, except where explicitly specified otherwise; the numerical-integration procedure and the associated convergence studies are detailed in Appendix B.

Cross-validation of simulation codes
The full impact of the beam-beam interaction on the beam-separation dependence of the luminosity-bias factor [L/L 0 ] Full BB (∆) can be broken down as follows.
The impact of the orbit shift can be expressed either in terms of the beam-beam induced distortion of the actual beam separation δ i , as discussed in Sec. 3.5.1 below, or in terms of an orbit-related luminosity-bias factor, denoted by [L/L 0 ] Orb (∆) for a given nominal separation ∆. In physically intuitive terms, and in the context of a weakstrong model such as B*B, the orbit shift results from applying to all particles in the witness bunch the same electromagnetic kick, computed from the field produced by the source bunch as a whole and averaged over all particles in the witness bunch. Since this is tantamount to a dipole kick, only the orbit of the witness bunch is affected; the size and shape of the macroparticle distribution remain invariant as |∆| increases.
The impact of the optical distortions is quantified in terms of the luminosity-bias factor [L/L 0 ] Opt (∆) first introduced in Sec. 2.3.2. It results from the combination of beamseparation and amplitude-dependent β-beating effects that modulate the RMS transverse beam size (Sec. 3.5.2), and of ∆-dependent bunch-shape distortions that affect the beam-beam overlap integral (Sec. 3.5.3).
The B*B and COMBI packages have been mutually benchmarked, and their results compared to those obtained either analytically (where possible) or using MADX. The cross-package comparisons of the beam-beam induced orbit shift and of the predicted optical distortions are based on the reference parameter set of Table 1, and assume that in the absence of the beam-beam interaction the transverse density distributions are strictly Gaussian. The consistency of B*B and COMBI results at higher beam-beam parameters is quantified in Sec. 3.5.4. All the results presented in this Section assume in addition that the beams collide only at the IP where the vdM scan is taking place; multiple-IP effects will be discussed in Sec. 4.6. Figure 5 displays the beam-beam-induced, singlebeam orbit shift (δ bb iB ) at the IP, as simulated by MADX, B*B and COMBI, and compares it to the analytical prediction. The three models agree extremely well among themselves, and their results are indistinguishable from those of the analytical prediction. During a luminosity-calibration session, the bunches typically collide at more than one IP (Sec. 4.6); however performing simultaneous beam-separation scans at different IPs is carefully avoided. Under these conditions and to an excellent approximation, the orbit shifts associated with slightly misaligned collisions at a non-scanning IP remain static during a beamseparation scan at another IP. Their impact on the beam separation at the scanning IP is therefore expected to remain constant during the scan, and has been neglected in the simulations reported in this paper.

Impact of amplitude-dependent
β-beating effects in the Gaussian-bunch approximation The beam-separation dependence of the transverse RMS beam-size ratio (squared for easier interpretation in terms of luminosity) is presented in Fig. 4. The results of B*B and COMBI agree to better than 2 × 10 −4 . At ∆ = 0, both predict a decrease in beam-size squared (or equivalently an increase in head-on luminosity) about half of that obtained using MAD-X. This is remarkably consistent with Eq. (22): in MAD-X, all particles in the witness bunch are subject to the same quadrupole-like force as the zero-amplitude particle, while in B*B and COMBI, the slope of the beam-beam force decreases and even changes sign as the betatron amplitude increases (Fig. 3), resulting in a smaller overall change in optical demagnification at the IP.
Amplitude detuning also explains the different ∆-dependence of the beam sizes, with that in MAD-X being more pronounced than in B*B and COMBI. At very large beam separation, where the electromagnetic force becomes similar to that of a distant, point-like charge, all three curves tend towards a common asymptote. As for the different evolution of the horizontal and vertical beam sizes during a scan, it simply reflects the tune-dependence apparent in, for instance, Eq. (18).
If one assumes that the optical distortions discussed above modify the transverse RMS bunch sizes at the IP without significantly affecting their initial, purely Gaussian shape, Eqs. (5)-(7) apply. The combination of the beam-beam induced orbit shift and of the optical distortions then yields the following expression for the luminosity-bias factor during a beam-separation scan in plane i (i = x, y): where σ x,y , Σ x,y , L (resp. σ 0 , Σ 0 x,y , L 0 ) are the RMS single-beam sizes, the inverse of the overlap integrals and the luminosity in the presence (resp. absence) of the beam-beam interaction.
Combining Eq. (24) with the beam-size ratios shown in Fig. 4 yields the beam-separation dependence of the luminosity-bias factor [L/L 0 ](∆) in the Gaussian-bunch approximation, as illustrated in Fig. 6 for MAD-X (dashed blue) and COMBI (dotted grey). These curves are not simply the inverse of the beam-size ratios displayed in Fig. 4, because they include in addition the impact, at a given nominal separation, of the beam-beam induced orbit shift. At zero and moderate separation (∆/σ 0 < 1). the luminosity bias [L/L 0 ] − 1 is positive and dominated by the dynamic-β effect; as the separation increases, the orbit shift, which is represented by the exponential term in Eq. (23), progressively takes over.

Impact of optical distortions on the beam-beam overlap integral
If the transverse-density distributions are sufficiently distorted by the non-linearity of the beam-beam force, the Gaussian-bunch approximation encapsulated in Eqs. (23)- (24) is no longer valid, and the luminosity bias must be calculated numerically from the beam-separationdependent overlap integrals of the macroparticle distributions. To this effect, an integrator module, detailed in Appendix B, has been developed for COMBI; B*B provides an equivalent functionality [8]. The resulting beam-separation dependence of the luminosity-bias factor is shown by the red curve and the black markers in Fig. 6. At each simulated scan step (indicated by the markers), B*B and COMBI agree to better than 2 × 10 −4 ; the difference is not systematic, but fluctuates around zero. The difference between the Gaussian-bunch approximation and the numerically calculated overlap integral demonstrates that the non-Gaussian tails induced by the beam-beam interaction have a significant impact on the luminosity as soon as ∆/σ 0 ≥ 0.2. The excellent agreement between the two multiparticle simulations also demonstrates that strong-strong beam-beam effects, that are modeled by COMBI but not by B*B, remain negligible in the low beam-beam parameter regime considered here.
Since the predicted orbit effect is identical in all simulations (Fig. 5), the differences between MAD-X on the one hand ( Fig. 6, dashed blue curve), and B*B/COMBI on the other (black markers/red curve), is entirely associated with optical distortions. At zero separation, the difference is entirely explained by amplitude detuning (Eq. (22)); once the separation increases, beambeam-induced non-Gaussian tails play a growing role, as illustrated by the difference between the red and grey curves. 5 In the context of the beam-beam correction strategy outlined in Sec. 2.3, it is convenient to separate the full luminosity bias presented in Fig. 6 into its optical-distortion and orbit-shift 5 The very small difference, at zero separation, between the grey and red curves in Fig. 6 may grow when the bunches collide not only at the scanning IP, but also at additional IPs (Sec. 4.6). In this case, and depending on the phase advance between IPs, beam-beam-induced non-Gaussian tails at these additional IPs may contribute noticeably to the overlap integral at the scanning IP, even for zero beam separation.
components. The luminosity-bias factor associated with optical distortions and denoted by [L/L 0 ] Opt is defined as the ratio, scan step by scan step, of the full luminosity-bias factor [L/L 0 ] Full BB and of the bias factor [L/L 0 ] Orb associated with the orbit shift: Its beam-separation dependence is presented in Fig. 7. While MAD-X overestimates the dynamicβ effect at zero separation, it strongly underestimates the optical distortions at medium and large beam separations. Since all models predict the same orbit effect, and since the optical distortions and the orbit effect impact the overlap integrals in opposite ways, their mutual cancellation is stronger in B*B and COMBI than in MAD-X, resulting in a net overall beam-beam correction of significantly smaller magnitude. This will be discussed quantitatively in Sec. 4.

Beam-parameter dependence of beam-beam corrections
The results in Figs. 6 and 7 were obtained using the reference parameter set in the right column of Table 1. This corresponds to a beam-beam parameter ξ = 2.59 × 10 −3 , at the low end of the ξ range explored during vdM scans at √ s = 13 TeV. The cross-validation of B*B and COMBI was therefore repeated with a larger bunch current and a smaller emittance, both typical of routine physics running and corresponding to a beambeam parameter value well beyond the vdM range. The beam-separation dependence of the opticaldistortion luminosity-bias factor at these two ξ values is presented in Fig. 8. While in perfect agreement at low ξ, the two codes exhibit hints of a small but systematic difference in the high-ξ regime when the beam separation becomes large enough. This discrepancy is attributed to the fact that B*B is intrinsically a weak-strong model, and therefore cannot account for coherent bunch oscillations, contrarily to COMBI. This interpretation is confirmed by the comparison of the tune spectra predicted by the two packages: even in the vdM regime, a π-mode peak is apparent in the COMBI spectrum, but is missing (as it should be) from the B*B spectrum.
Quantitatively, it turns out that throughout the vdM regime and even in the high-ξ case considered here, the discrepancy is of no practical significance. Since it manifests itself only in the tails of the beams (∆ > 2σ 0 ), where the particle density is low, the difference, between B*B and COMBI, in beam-beam-induced distortions has only a very small impact on the integral of the scan curves, i.e. on the perturbed transverse convolved beam sizes (Eq. (4)). This can be quantified using the methodology that will be presented in Sec. 4.2.4. For the high-ξ setting illustrated in Fig. 8, the difference between B*B and COMBI translates into a systematic uncertainty of 0.04% on the absolute luminosity scale. In typical vdM scans, where ξ is lower by a factor of about 1.5, the inconsistency becomes negligible; there the two simulation codes remain in more than adequate agreement, validating the use of the less   Table 1. The horizontal axis is the nominal beam separation in units of the unperturbed transverse beam size σ 0 . The MAD-X (dashed blue line with markers), weak-strong B*B (black triangles) and strong-strong COMBI (red curve) optical-distortion curves are obtained by dividing the corresponding full beam-beam bias factor shown in Fig. 6 by that predicted analytically for the orbit effect alone (green curve).  • the interaction-region (IR) configuration (beam energy, nominal or measured β * and crossingangle values); • the unperturbed tunes Q x and Q y (or equivalently the unperturbed fractional tunes q x and q y ), i.e. the values of the horizontal and vertical tunes with the beam-beam interaction switched off at the IP where the simulated scans are taking place. Physically, these correspond to the tune values that would be measured, for the bunch pair under study, before the beams are brought into collision at the scanning IP considered; • the measured parameters of the bunch pair under study: bunch population, horizontal and vertical beam sizes or emittances, as well as the bunch length in case of a non-zero nominal crossing angle.
Since vdM calibrations must be performed on a bunch-by-bunch basis, with possibly over 100 colliding-bunch pairs and typically five to ten xy scan pairs per scan session and per IP, the approach sketched above can become unwieldy from the computing viewpoint. This motivated the development of a much lighter technique, based on the fact that at least in the vdM regime and under some simplifying assumptions, beam-beam biases, to a very good approximation, scale almost linearly with the beam-beam parameter ξ. This makes it possible to construct a simple polynomial parameterization of the luminosity-bias curves, that is extracted from B*B or COMBI simulations over a grid in beam-parameter space and is applicable to most cases of practical interest for pp vdM calibrations at the LHC. The impact of violating the underlying assumptions is accounted for either by a contribution to the systematic uncertainty that is associated with the beam-beam correction procedure, or, in the case of multi-IP effects, by a simulation-guided adjustment to the parameterized correction.
Some cases do not lend themselves to the parameterization technique, such as off-axis vdM scans, diagonal scans, or vdM calibrations performed with a large crossing angle (which is unavoidable at the LHCb IP).
An off-axis vdM scan is a horizontal (or vertical) beam-separation scan where the beams are partially separated in the non-scanning plane, i.e. in the vertical (or horizontal) direction. A diagonal scan is one in which the beams are scanned transversely along an inclined straight line in the x-y plane, rather than only along either the x or the y axis. Such "generalized", one-dimensional vdM scans are sometimes used to measure and correct for non-factorization effects [13,35,36]. While the orbit-shift correction can still be calculated analytically using a more general form of the Bassetti-Erskine formula, the increased dimensionality of the parameter space makes it impractical to invest in general enough and precise enough a parameterization of optical-distortion effects. In such a case, deriving the corrections from a large set of fully simulated scans, labor-and CPU-intensive as it may be, appears more tractable by comparison.
Allowing non-zero crossing angles also increases the dimensionality of the parameter space, with similar implications. Therefore, except where explicitly stated otherwise, the results offered in the remainder of this report are restricted to the zero (or moderate) crossing-angle case, which covers the needs of the ATLAS and CMS experiments (as well as those of ALICE, at the cost of a slight increase in systematic uncertainty).
This Section is organized as follows. The parameterization approach mentioned above is described in Sec. 4.2. The associated "fully symmetric Gaussian-beam configuration" assumes that in the absence of any beam-beam interaction, the colliding bunches in beam 1 and beam 2: • are beam-beam symmetric, i.e. equally populated (n 1 = n 2 ), of the same transverse size (σ i,1 = σ i,2 , i = x, y), and therefore subject to the same beam-beam parameter (ξ i,1 = ξ i,2 , i = x, y). The impact of violating these assumptions is evaluated in Sec. 4.7.

Beam-beam correction procedure in the fully symmetric Gaussian-beam configuration
The parametrization strategy relies on the fact that at fixed tunes, beam-beam induced biases to the visible cross-section depend only on the beambeam parameter ξ (Sec. 4

Validation of the scaling hypothesis
To characterize the scaling properties (or lack thereof) of beam-beam biases during vdM scans, pairs of horizontal and vertical beam-separation scans were generated under the assumptions above using B*B. 6 The beam-beam parameter spanned the range ξ ∼ 0.002 − 0.008, thereby covering from vdM scans with low-brightness proton beams at the low end, to conditions slightly above routine physics running at the high end. The input beam energies, bunch populations, emittance and β * values were chosen to be representative of vdM scans during LHC Runs 1 and 2 (Sec. 2.2); also included was a high-ξ setting representative of Run-2 high-luminosity physics running. The unperturbed fractional tunes were constrained to satisfy q y = q x +0.010 so as to reflect routine LHC collision settings. The primary output of the simulation is the beam-separation dependence of the luminositybias factor [L/L 0 ] Full BB . An example is presented  in Fig. 9, with a value of ξ chosen to lie roughly in the middle of the vdM range ( Table 1). The orbit-shift bias factor [L/L 0 ] Orb (green squares) is computed analytically: where ∆ is the nominal separation in the scanning plane, δ i is computed using Eqs. (14) to (16), The optical-distortion bias factor [L/L 0 ] Opt (red triangles) can then be calculated using Eq. (25).
For given values of q x and q y , beam-beaminduced distortions scale with the beam-beam parameter, in the sense that they depend only on ξ. This can be proven mathematically, both for the orbit shift (Eqs. (14)- (15)) and for the dynamicβ effect at zero separation (Eq. (18)). That this remains true at any beam separation, within the constraints of the fully symmetric beam configuration and over the ξ range specified above, can only be demonstrated by simulation. To this effect, luminosity-bias curves such as those presented in Fig. 9 were compared for different combinations of bunch populations, beam energies, β * and emittance values that all correspond to a given value of ξ; they were found to be identical within the statistics of the simulation. An example is presented in Fig. 10, which displays the difference, at each scan step, between the luminosity-bias factor computed for three distinct sets of beam parameters, and that associated with the parameters used in Fig. 9; all four parameter sets correspond to ξ = 0.004. The error bars represent the statistical uncertainty associated with the randomization of initial conditions in the B*B code [8]; they are computed as the error on the mean over eight different runs, with different random seeds, for each beam separation and each parameter set. The differences in luminosity-bias values between the four beam-parameter sets and the reference set are statistically consistent with zero, and never exceed 10 −4 . The exercise was repeated for a range of ξ values, leading to the conclusion that at fixed Q x and Q y , the luminosity-bias curves indeed depend only on ξ, to better than 10 −4 on the absolute luminosity scale.

Beam-beam parameter and tune dependence of the optical-distortion correction
The combined beam-separation and ξ dependence of the luminosity-bias factor [L/L 0 ] Opt is illustrated in Fig. 11  Since both the orbit shift (Eq. (15)) and the dynamic-β effect (Eq. (18)) depend on the unperturbed tunes, and as part of investigating the systematic uncertainties associated with beambeam corrections, their sensitivity to the input tune values was characterized by extending the above-described simulations to an area in the (q x , q y ) plane that encompasses the full range of operational tune settings and of beam-beam tune shifts expected during vdM sessions. Horizontal and vertical beam-separation scans were simulated, for several values of ξ, over a grid 7 bounded by 0.2975 < q x < 0.3100, with the vertical unperturbed fractional tune set to q y = q x + δ q where δ q = 0.0075, 0.0100, 0.0125.
The combined beam-separation and tune dependence of the luminosity-bias factor [L/L 0 ] Opt is illustrated in Fig. 12 for ξ = 7.83 × 10 −3 . For a given nominal separation, the lower the tunes, the more they approach the quarter integer, and therefore the smaller (i.e. the closer to unity) the optical-distortion bias factor becomes [4]. Similarly, the lower the tunes, the smaller (i.e. the further away from unity) the orbit-shift bias factor [L/L 0 ] Orb (not shown), albeit with a weaker tune dependence. As a result, the tune dependence of the full beam-beam-bias factor [L/L 0 ] Full BB is slightly steeper than that associated with optical distortions only: the lower the tunes, the faster [L/L 0 ] Full BB drops in magnitude (Fig. 13), and therefore the larger the overall beam-beam correction to the luminosity scale. In addition, the lower the tunes, the less the optical distortions contribute to the overall beam-beam bias to σ vis .

Practical implementation of parameterized beam-beam corrections
Concretely, given a measured vdM-scan curve such as that displayed in Fig. 1 (bottom), the beambeam correction procedure can be implemented as follows.
1. At each scan step in the i scanning plane (i = x, y), correct the nominal separation ∆ i for the deflection-induced separation shift δ bb i , as detailed in Sec. 2.3.1, using Eqs. (14)-(16). 2. At each scan step in the i scanning plane, correct the measured rate R i (∆ i ) for the opticaldistortion bias as detailed in Sec. 2.3.2, using Eq. (17). The beam-separation dependence of the luminosity-bias factor [L/L 0 ] Opt is either extracted directly from a set of B*B or COMBI simulations, or, if the assumptions listed in Sec. 4.1 can be considered valid, simply computed using the polynomial parameterization described below. The latter is valid only for on-axis, one-dimensional vdM scans with zero crossing angle; other configurations, and in particular off-axis one-dimensional scans as well as two-dimensional grid scans in the (δ x , δ y ) beam-separation plane, require dedicated simulations.
The simulations described in Sec. 4.2.2 were used to parameterize the optical-distortion luminosity-bias factor [L/L 0 ] Opt in bins of normalized nominal separation over the range 0 ≤ ∆ i /σ 0 ≤ 6, separately for x and y scans, by polynomials of the form Here the index k (k = 1, 25) labels the nominalseparation bin in the scanning plane considered. The round-beam equivalent beam-beam parameter ξ R is defined as where the single-beam transverse RMS beam size σ 0 in Eq. (13) is approximated by its measured 8 8 It will be shown in Sec. 4.2.4 that in the vdM regime, beam-beam effects, if left uncorrected, lead to a percent-level underestimate of Σx and Σy, and therefore of σ R . The use of the measured transverse beam size -as opposed to that of the experimentally inaccessible unperturbed beam size -in round-beam equivalent σ R = Σ x Σ y /2, thereby partially accounting for a potential ellipticity of the beams at the IP. The input fractional-tune values (q x , q y ) are those of the unperturbed tunes defined in Sec. 4.1. In the case where bunches collide only at the IP where the scan is taking place, i.e. remain fully separated, transversely and/or longitudinally, at the other three IPs, these unperturbed tunes are identical to the nominal LHC tunes, or equivalently to the values one would measure before the beams are put in collision at the scanning IP. If, however, some bunches also collide at one or more other IPs, the unperturbed tunes input to the parameterization for those bunches must take into account the beam-beam tune shift arising from collisions at non-scanning IPs; a prescription to this effect is offered in Sec. 4.6.
The parameterization above, that amounts to a second-order Taylor expansion in ξ R , q x and q y , was found sufficient to approximate the exact simulation results to better than 10 −3 on [L/L 0 ] over most of the grid of simulated points in (ξ, q x , q y , ∆) space.. The function defined by Eq. (27) is linear in all 10 parameters p lk . In each nominalseparation bin k, the set of parameters p 0k , . . . , p 9k is determined by a weighted linear least-square fit. A smooth dependence of [L/L 0 ] on the nominal separation, for given values of ξ R , q x and q y , can be achieved by, for instance, cubic spline interpolation between adjacent separation bins. Parameterizations based on the same functional form were also constructed for the orbit-shift and full beambeam bias factors [L/L 0 ] Orb and [L/L 0 ] Full BB , and achieved similar numerical accuracy.
Tabulated values of the parameters p lk in separation steps of 0.25 in ∆/σ R , for both vertical and horizontal scans and over the full range of ξ R and tune values defined earlier, have been made available to all LHC experimental Collaborations. The optical-distortion parameters, which are the most estimating ξ R therefore leads to a slight overestimate of the beam-beam parameter input to the simulation. However, since the overall beam-beam bias on the effective cross-section is typically less than 1% (Fig. 16, black curve), a sub-percent overestimate of ξ R biases the overall beam-beam correction to σvis by 1-1.5% of itself. This can usually be neglected in view of the much larger fractional systematic uncertainties detailed in Sec. 5. Alternatively, the correction can be iterated upon, by using as input to the second iteration a ξ R value based on the corrected effective single-beam size obtained from the first iteration.
24 useful, are documented in Appendix C, and are publicly accessible in computer-readable form [37].

Separation-integrated estimates of beam-beam corrections to vdM calibrations
The beam-beam correction procedure detailed in the preceding Section is designed to be applied to measured vdM-scan curves, one scan step at a time. One then extracts the beam-beam corrected peak rate and convolved beam sizes needed to calculate the visible cross-section (Eq. (10)) using a carefully chosen fit function. More often than not, the latter must diverge from a perfect Gaussian in order to faithfully model the data over the full beam-separation range.
In order to provide consistency checks on this intricate analysis chain, as well as a physically intuitive, albeit approximative, breakdown of the impact of individual beam-beam effects on the luminosity calibration under study, the same procedure can be applied to an x-y pair of hypothetical beam-separation scan curves, that suffer from no statistical fluctuations and that, in the ξ R → 0 limit would be perfectly Gaussian: Here and in the remainder of this paper, the "0" index indicates a "nominal" quantity, i.e. one for which beam-beam effects have been fully turned off in the simulation. At a given nominal separation ∆ i in scanning plane i (i = x, y), and with the beams transversely centered on each other in the non-scanning plane, the luminosity in the presence of, say, the full beam-beam effect is then calculated as One then defines "figures of merit" (FoMs) that characterize the impact of beam-beam effects on vdM observables.
1. The peak-rate bias factor characterizes the dynamic-β effect at zero beam separation: and is sensitive to optical distortions only. 2. The beam-size bias factors characterize, in each plane, the beam-beam impact on the convolved transverse beam size: where the overlap integrals are calculated numerically, with their integrands smoothed by cubic-spline interpolation. 3. The σ vis -bias factor characterizes the beambeam impact on the absolute luminosity scale, and is given by the product of the previous three FoMs: As an illustration, Fig. 14 displays the beambeam parameter dependence of the peak-rate bias factor for the full beam-beam effect (black circles), the orbit shift only (green squares), and the optical distortion only (red triangles). Since at zero separation there is no orbit shift, only optical distortions matter in this case.
The evolution of the horizontal beam-size bias factor is presented in Fig. 15. While orbit shifts, if left uncorrected, result in an underestimate of Σ x and Σ y , optical distortions have the opposite effect, resulting in a partial cancellation of the bias. This is even more apparent in the visiblecross-section bias curves (Fig. 16). While the full beam-beam effect on the visible cross section remains below 1% over the whole ξ range (black where j labels the type of correction (orbit shift only, optical distortion only, or full beam-beam). Even though the linear term is clearly dominant, a quadratic term is needed to describe the evolution of all three FoMs to satisfactory precision over the full range of beam-beam parameter values. The impact of beam-beam effects on the absolute-luminosity scale can be expressed equivalently by either the σ vis bias factor σ vis / σ 0 vis , the visible cross-section bias σ vis / σ 0 vis − 1 (typically a fraction of a percent), or a multiplicative correction factor ( σ vis / σ 0 vis ) −1 to the raw visible cross-section. This FoM approach, that lends itself to a simple polynomial parameterization in terms of ξ R , q x and q y , offers the advantage that its results are easy to interpret, as well as insensitive both to the quality of the fit to the measured scan curves and to rate fluctuations, from one scan step to the next, due for instance to counting statistics or beam-position jitter. This technique also simplifies the evaluation of some of the systematic uncertainties discussed in the remainder of this chapter. It is, however, not recommended for determining the central value of the beam-beam correction on real data, since it ignores the deviations of the actual scan curves from a perfect Gaussian.

Impact of non-Gaussian unperturbed transverse beam profiles
All the simulation results presented in this report so far make the explicit assumption that in the absence of the beam-beam interaction at the scanning IP, the unperturbed transverse beam profiles, i.e. the particle density functionsρ B (x, y) (B = 1, 2) in Eq. (1), can be perfectly modeled by the uncorrelated x-y product of two single, one-dimensional Gaussians. Should this not be the case, the beam-separation dependence of the luminosity-bias curves will deviate from that presented in Sec. 4.2. This is due both to the modified spatial dependence of the field generated by the source bunch, and to the fact that the fraction of particles in the witness bunch that experience a given electromagnetic kick is different from that in the pure-Gaussian case. Both B*B and COMBI accept as input unperturbed transverse-density distributions, functional forms that are more general than a single Gaussian. In order to assess the potential impact, on beam-beam corrections, of non-Gaussian beams, a realistic model ofρ B (x, y) is needed, or at least one that represents the "worst-case" deviation from the ideal Gaussian shape while remaining representative of actual beam conditions during vdM-calibration sessions at the LHC (Sec. 4.3.1). Taking into account potentially non-factorizable unperturbed density distributions requires a minor generalization of the beam-beam correction formalism (Sec. 4.3.2). While the influence of non-Gaussian tails on the beam-separation dependence of the luminosity-bias factors appears sizeable, the resulting bias on the beam-beam corrections to the visible cross-section remains moderate enough to be treated as a systematic uncertainty (Sec. 4.3.3).

Single-bunch models
In a hadron collider, transverse-density distributions cannot be calculated from first principles, and existing beam-profile monitors, such as wire scanners or synchrotron-light telescopes, have too limited a dynamic range to provide sufficiently precise beam-tail measurements. The only remaining experimental handle, therefore, is that provided by non-factorization analyses such as those described in, for instance, Refs. [13,38,39,40,41]. In this approach, the bunch-density distributions are modeled by the sum of two or three threedimensional Gaussians. A simultaneous fit to the beam-separation dependence, during a vdM-scan pair, not only of the luminosity but also of the position, size, shape and orientation of the luminous region, makes it possible to estimate the single-bunch parameters of each colliding-bunch pair.
In order to quantify the largest plausible impact of non-Gaussian beam profiles on the beam-beam corrections calculated in Sec. 4.2, single-bunch parameters were extracted from the results of non-factorization analyses of a 2012, a 2017 and a 2018 vdM session. The former is representative of Run-1 vdM conditions at the ATLAS IP, up to and including July 2012 [13], the latter two of Run-2 vdM sessions at the CMS IP, for which "bunch tailoring" in the injector chain [13,42] significantly reduced non-Gaussian tails. In both cases, a "worst-case" parameter set (in terms of deviations from a perfectly Gaussian shape) was selected from the fitted parameters of the analyzed colliding-bunch pairs. In order to separate the impact of non-Gaussian profiles discussed in this Section, from that of elliptical beams (Sec. 4.4) and of beam-beam imbalance (Sec. 4.7), as well as to maximise the sensitivity of the study, in each parameter set the most non-Gaussian transverse profile is assigned to both planes and both beams.
In a first step, the particle density distribution input to B*B is chosen to be factorizable by construction, and modeled by the product of two uncorrelated double Gaussians: where B = 1 or 2, the labels "n" and "w" refer to the narrow and wide components of the distribution, and their relative population is constrained by w w = 1 − w n . The functional form reflects the assumptions of unperturbed round beams and of equal particle-density distributions (ρ 1 (x, y) =ρ 2 (x, y)). The three parameter sets input to the B*B simulation are listed in Table 2.
The 2012 (2017) parameter set clearly results in the most (the least) non-Gaussian shape, as evidenced qualitatively by the combination of the largest (smallest) beam-size ratio σ w /σ n and the largest (smallest) weight of the wide component w w , and as quantified by the single-beam kurtosis computed from the parameters listed in the Table. 9 The functional form chosen forρ B (x, y) makes it possible to compute analytically the unperturbed convolved transverse beam size Σ 0 from the values of σ n , σ w , w n and w w , using Eq. (4): these are listed in the top half of Table 3 (rows 3 to 6). The Gaussian-equivalent single-beam size σ 0 SG , i.e. the transverse R.M.S. width of perfectly Gaussian bunches that would yield the same unperturbed convolved beam size as the non-Gaussian bunches studied here, is given by σ 0 SG = Σ 0 / √ 2 (Eq. (6)). These definitions naturally lead to using the "Gaussian-equivalent" beam-beam parameter ξ, inferred from σ 0 SG using Eq. (13), as the common metric in which to quantify the impact of non-Gaussian tails.

Beam-beam-correction formalism in the presence of non-factorization
In a second step, and in order to explore the potential impact of non-factorization on beambeam corrections, the simplest two-dimensional, non-factorizable double Gaussian: is chosen as an alternative density-distribution model to be input to B*B. Since by construction such a distribution cannot be expressed as 9 Deviations of a distribution from the strictly Gaussian shape can be characterized by its kurtosis. This statistic is defined as Γ2 = m4/m 2 2 −3, where m4 and m2 are, respectively, the fourth and the second moment of the distribution [43]. The kurtosis is zero for a Gaussian, positive for a leptokurtic distribution with longer tails, and negative for a platykurtic distribution with tails that fall off more quickly than those of a Gaussian. the product of uncorrelated horizontal and vertical components, the concept of one-dimensional convolved beam size defined by Eq. (4) ceases to be meaningful, and only the convolved-beam-size product [Σ x Σ y ] 2D introduced in Eq. (8) remains amenable to physical interpretation. Accordingly, the two beam-size bias factors Σ i /Σ 0 i (i = x, y) introduced in Eqs. (31)-(33) must be replaced by the single quantity where the scan and the integration are carried out over a grid in the two-dimensional nominal-separation space (∆ x , ∆ y ) that extends over ±5 σ 0 SG in both directions. The σ vis -bias factor defined by Eq. (34) then takes the more general form where the subscript "2D" indicates that the twodimensional formalism introduced in Eq. (8) is used throughout the analysis. Using ξ as a common metric for both the factorizable and the non-factorizable case also requires a generalization of the formula that expresses, in terms of single-beam parameters, an observable that is akin to the transverse convolved bunch size and that can be measured in a standard, one-dimensional vdM scan. For instance, for a vertical scan with the beams centered on each other in the horizontal plane, the measured apparent convolved width is given by: which reduces to Eq. (4) whenρ(x, y) = ρ x (x) ρ y (y). The subscript "F" emphasizes that while measurable, this observable cannot be used in Eqs. (5) and (10) unless the density distributions are factorizable. The importance of this distinction is illustrated by the significant differences between the top half (rows 3-6) and the bottom half (rows 7-10) of   Table 3: Scan parameters representative of strongly non-Gaussian beam profiles during ATLAS and CMS pp vdM scans in 2012, 2018 and 2017, calculated analytically from the single-bunch parameters listed in Table 2.  (1)). These differences, the size of which appears correlated with the single-beam kurtosis, are indicative of the magnitude of the non-factorization biases that can be introduced, even in the absence of the beam-beam interaction, by applying the one-dimensional, factorizable FoM formalism of Sec. 4.2.4 to non-factorizable beams. In contrast, the analysis strategy encapsulated in Eq. (39) guarantees that non-factorizable correlations embedded in the unperturbed density distributions cancel in the σ vis2D ratio; the latter, therefore, reflects the impact of beam-beam biases only. The study of non-factorization biases proper lies beyond the scope of the present report.

Impact of non-Gaussian tails on beam-beam-induced biases
In the case of factorizable density distributions, the impact of non-Gaussian tails on beam-beam corrections can be visualized by comparing the beam-separation dependence of the luminositybias factors [L/L 0 ] for a given set of non-Gaussian beam parameters, to that computed for perfectly Gaussian beams that in the absence of beambeam effects, would yield the same value of the unperturbed convolved beam size Σ 0 , and therefore the same head-on luminosity. This is illustrated in Fig. 17 for the orbit-shift effect. In the single-Gaussian case (red curve), the steeper drop of [L/L 0 ] Orb with increasing separation, and its asymptotic behavior for ∆ y ≥ 300 µm, indicates that at such large separations (∆ y > 4 σ 0 SG ), the number of particles contributing to the source field becomes negligible, so that the latter increasingly appears to the witness bunch as a point source. In the non-Gaussian case (black curve), the charge distribution is more spread out, and the interplay between the source-field distribution and the density distribution in the witness bunch produces a less intuitive separation dependence.
The situation is reversed for optical distortions (Fig. 18): the Gaussian beams are subject to larger optical distortions at large beam separation (red curve), presumably because in such configurations the tail particles in the witness bunch experience a more intense source-field gradient when close to the center of the source bunch. The tails of the non-Gaussian beams, being subject to a more diffuse source-field distribution, experience weaker gradients and correspondingly smaller optical distortions (black curve).
In the case of non-factorizable distributions, a physically intuitive graphical representation of beam-beam biases becomes impractical, since [L/L 0 ] depends on both ∆ x and ∆ y , and the curves in Figs. 17 and 18 must be replaced by surfaces.
The impact of non-Gaussian tails on the corresponding vdM calibrations is best summarized by the σ vis -bias factors presented in Table 4. These are all computed using Eqs. (38)- (39), in order to ensure that factorizable and non-factorizable models are treated identically; it has been verified that for the factorizable models, the one-dimensional FoM formalism detailed in Sec. 4.2.4 yields the same results within 0.01% of σ vis , well within the expected accuracy of the overlap-integral calculation (Sec. 3.3). Listed in this table is the full beam-beam bias for either factorizable or non-factorizable double-Gaussian models, for the corresponding single-Gaussian model (that is factorizable by construction), as well as for their difference DG−SG. This difference is also listed after scaling it linearly to a common, arbitrary reference value of ξ = 5.6·10 −3 , that lies at the upper end of the beam-beam parameter range covered by Run-2 vdM scans. This simple-minded procedure, inspired by the roughly linear dependence of the σ vis bias factor on ξ apparent in Fig. 16 Globally, non-Gaussian tails modify the σ vis bias by an amount that -for the models chosen here -varies from less than 10% to about half of the bias itself. Predictably enough, the larger the single-beam kurtosis, the larger the relative impact of the tails. The non-factorizable models seemingly predict a smaller influence of the tails than their factorizable counterpart; this may however depend on the particular form of x-y correlation embedded into the specific functional form chosen forρ B (x, y). A practical proposal to translate the results listed in Table 4 into a realistic systematic uncertainty is offered in Sec. 5.2.1. Repeating the above studies using more refined single-beam models extracted from future vdM sessions may help refine the understanding of such biases.  Table 4: Visible cross-section bias associated with the full beam-beam effect, calculated using B*B, for double-(DG) and single-(SG) Gaussian beams yielding the same Gaussian-equivalent beam-beam parameter ξ, and difference between the two (DG−SG). The unperturbed DG transverse-density distributions are described either by a factorizable or by a non-factorizable function; their parameters are listed in Table 2. The bunches are assumed to collide at the scanning IP only.

Impact of beam ellipticity at the interaction point
In the fully symmetric beam configuration (Sec. 4.1), the beams are round at the IP, and the beam-beam interaction can be described by a single parameter ξ (Eq. 13). In the more general case where the transverse beam profile is an upright 10 ellipse, for instance because of unequal horizontal and vertical emittances or β * values, the horizontal and vertical beam-beam parameters ξ x and ξ y (Eq. 12) no longer coincide, and the parametric approach described in Sec. 4.2.3 is, strictly speaking, no longer applicable.
During vdM scans at the LHC, however, the round-beam hypothesis is only moderately violated: the measured horizontal and vertical IP-β functions typically differ by 10% or less, the corresponding emittances by 30% at most, and the increase in projected transverse beam size associated with a non-zero crossing angle at the ALICE and LHCb IPs does not exceed 15%. To assess the impact of beam ellipticity on beam-beam corrections, B*B is used to simulate vdM scans over a square grid in (Σ 0 x , Σ 0 y ) space (Fig. 19), with as input unperturbed beams that satisfy the assumptions listed in Sec. 4.1 except for the round-beam hypothesis. The chosen Σ 0 x range is slightly wider than that observed during Run-2 vdM scans at √ s = 13 TeV; the vertical/horizontal aspect ratio Σ 0 y /Σ 0 x is varied from 0.7 to 1.4, corresponding to sampling the corresponding emittance ratio from 0.5 to 2.0, well beyond the range observed during vdM-calibration sessions. The horizontal and vertical beam-beam parameters cover the range 2.6 × 10 −3 < ξ x , ξ y < 5.1 × 10 −3 (Fig. 20). The grid is chosen such that for a constant value of the round-beam equivalent beam-beam parameter ξ R , the horizontal and vertical beam-beam parameters sample the full span allowed by the range of chosen Σ 0 x and aspect-ratio values. For example, the group of five points at ξ R ≈ 3.6 × 10 −3 in Fig. 20 corresponds to the five diagonal points in Fig. 19; the residual curvature is due to the fact that the grid is built from equally spaced (Σ 0 x , Σ 0 y ), rather than (ξ x , ξ y ), values. The resulting peak-rate, beam-size and visible cross-section bias factors (Eqs. (30)-(34)) for elliptical configurations are then compared to the same quantities calculated using the round-beam parameterization (Eq. 35) at the corresponding value of ξ R . The results are illustrated in Fig. 21 for the visible cross-section. On the scale shown here, the bias ( σ vis / σ 0 vis − 1), and therefore the corresponding beam-beam correction, seem to depend on ξ R only; stated differently, the beam-beam-induced biases appear rather insensitive to the aspect ratio. Upon closer inspection however, all four FoMs do reveal an actual dependence on the Σ 0 y /Σ 0 x aspect ratio, albeit a weak one. Since the orbit-shift effect can be accurately calculated, even for elliptical beams, using the Bassetti-Erskine formalism [6], only opticaldistortion biases matter. The difference between the bias calculated by B*B for a given set of elliptical-beam parameters, and that extracted from the round-beam parameterization using the corresponding value of ξ R , is displayed in Fig. 22. Its dependence on the vertical/horizontal aspect ratio is close to linear, and remains below 0.03% over the full range of beam parameters considered in this study. Its dependence on ξ R is illustrated   by the spread over the three parameter sets clustered at Σ 0 x /Σ 0 y ≈ 0.8 and 1.2. The beam-beam parameter values considered in this study, and that are associated with elliptical-beam configurations, span the range 3.0×10 −3 < ξ R < 4.2×10 −3 : this is representative of, but does not quite cover, the full ξ range listed in Table 1. Even if the mean value of ξ R was doubled, however, the impact of potentially elliptical beam profiles on the beambeam correction to σ vis would still remain well below 1 permil, and can therefore be treated as a systematic uncertainty.

Impact of non-zero crossing angle
In the presence of a significant crossing angle, the analysis becomes considerably more involved, and a full exposition thereof would exceed the scope of the present report. The discussion below is therefore focussed on assessing the beam-beam correction uncertainty associated with vdM scans at zero nominal crossing angle, during which the actual, residual crossing angle does not exceed few tens of microradians: this covers most luminositycalibration sessions at IP1 (ATLAS) and IP5 (CMS). The same treatment, albeit at the cost of a small additional uncertainty, remains applicable to scans at IP2 (ALICE) (and to some non-standard sessions at IPs 1 and 5) that involve nominal crossing angles θ c in the 140-400 µrad range. At IP8 however, the LHCb requirement for an unambiguous separation, in the crossing plane, of beam-gas vertices from beam 1 and beam 2, leads to θ c settings of up to 1100 µrad. Such large crossing angles would require simulation studies more extensive than can be reported here.
In the absence of any beam-beam interaction, the geometric factor F defined in Eq. (11), and generalizations thereof [8,10], provide an accurate description of the combined impact of a non-zero crossing angle and of a finite bunch length on the unperturbed apparent transverse beam size in the crossing plane, as well as of the beam-separation dependence of the unperturbed luminosity.
In the presence of the beam-beam interaction, the implications of a non-zero crossing-angle for the calculation of beam-beam biases fall in two conceptually distinct categories.
• Geometrical effects in the transverse plane that depend on the crossing angle and on the bunch length, such as modifications of the transverse beam-beam overlap and the calculation of the effective strength of the source field experienced by the particles in the witness beam, are adequately modeled [8] by the B*B package and well reproduced by COMBI simulations. B*B simulations, however, are intrinsically limited to four dimensions (horizontal and vertical positions and angles), since the longitudinal variables are "projected out" and the associated dynamics implicitly neglected. • Longitudinal dynamics, such as synchrotron motion or beam-beam-induced transverseshape variations along the bunch, can be modeled by COMBI using the 6D beam-beam strong-strong model [33]. This model divides particle distributions into slices along the length of each bunch. The kick experienced by each macroparticle in one beam is computed based on the statistical moments of the charge distributions in each slice of the opposing bunch; the computation is repeated for the macroparticles in the other beam.
In order to properly account for longitudinal effects and assess their importance, COMBI is chosen here to characterize the interplay of crossingangle and beam-beam effects. This first study is exploratory in nature and restricted to the vdM regime.
As predicted analytically [8,44], the magnitude of the orbit shift extracted from COMBI simulations for a purely horizontal crossing angle of several hundred microradians, is found consistent with that computed using the Bassetti-Erskine formula with, as input: • in the crossing (horizontal) plane, the effective convolved beam size given by Σ x,ef f = F × Σ x , where F is the geometric factor defined in Eq. (11) and Σ x is the horizontal convolved beam size at zero crossing angle, given by Eq. (6); • in the non-crossing (vertical) plane, the vertical convolved beam size Σ y .
The dependence of the beam-beam-induced visible cross-section bias on the full crossing angle is illustrated in Fig. 23 for a beam-beam parameter setting typical of Run-2 vdM scans. For the configurations shown in this figure, the geometric factor spans the range 1.0 < F < 1.09. Since this factor impacts both σ 0 vis and σ vis in an almost identical manner 11 , it cancels in the crossing-angle dependence of the ratio σ vis / σ 0 vis . The absolute value of the σ vis bias increases with the angle, even though the "effective" beam-beam parameter in the crossing plane (i.e. that calculated from the effective transverse beam size) decreases. The difference between the horizontal-and verticalcrossing configurations is due to the difference between the x and y tunes. Also notable is the sensitivity to β * : the crossing-angle dependence of the beam-beam bias is weaker for β * = 24 m than for β * = 19 m. The results are insensitive to the assumed value of the synchrotron tune, as well as to the chromaticity of the ring lattice (reducing the chromaticity to zero from its nominal setting of 15 units impacts the results by less than 0.01 %). This suggests that the effect is mostly geometrical in nature, as confirmed by the qualitative agreement of the COMBI results with those of a preliminary B*B study that covered a similar parameter space.
In practical terms, the impact of the crossing angle on the estimated beam-beam correction can be taken as the difference in visible cross-section bias between θ c = 0 and the actual crossing-angle value. For a nominally zero crossing angle, the actual angle, as measured by the BPM system or by K-modulation in the final-triplet quadrupoles, is of the order of 10 µrad: the associated shift in the beam-beam bias, as read off Fig. 23   angle would underestimate the beam-beam correction by less than 0.02%, a value small enough to be neglected. Crossing angles in excess of 500 µrad, however, induce still-to-be understood biases of 0.1% or more (Fig. 23), and may require further study.

Motivation and methodology
All results presented in this paper so far assume that the beams collide at the scanning IP only, i.e. that they are fully separated, either transversely or longitudinally, at the other three IPs of the LHC. In practice however, this is often not the case: to make optimal use of beam time during luminosity-calibration scans at, say, the CMS IP, the experimental detectors at some or all of the other IPs typically collect collision data at fixed (most often zero) beam separation, thereby significantly affecting the unperturbed-tune spectra and therefore the magnitude of beam-beam induced biases. For a given colliding-bunch pair at a given scanning IP, the number of non-scanning interaction points N NSIP is defined as the average number of head-on collisions experienced, at interaction points other than the scanning IP, by the two members of the pair. The parameter N NSIP ranges from zero when beams are fully separated at all but the scanning IP, to N NSIP = 3 when both bunches in the pair collide at three other IPs; it can assume half-integer values when the two members of the pair experience different numbers of additional collisions.
Because IP1 (ATLAS) and IP5 (CMS) are separated by half a ring circumference, all bunches that collide in IP1 also collide in IP5, unless fully separated in the transverse plane at the non-scanning IP. Each member of such a collidingbunch pair can optionally also collide in IP2 (ALICE) and/or IP8 (LHCb), both of which are located 1/8th of the ring away 12 from IP1; however its collision "partner" at these IPs cannot be the same as that at IPs 1 and 5, because of the ring layout (Fig. 24). Injected-bunch patterns are optimized, on a fill-by-fill basis, to achieve the desired sharing of collisions (or absence thereof) among the four experimental detectors.
The parameter N NSIP , therefore, depends on the scanning IP considered, and, for a given scanning IP, can vary along the bunch string in a manner that depends on the injected-bunch pattern. From 2010 to 2017, the most frequent configuration during pp vdM scans at the ATLAS and CMS IPs corresponded to N NSIP = 1 at IPs 1 and 5, i.e. to bunches that collide in IP1 and IP5 only; from 2018 onwards, all four scanning IPs had to deal with multiple non-scanning IPs, for reasons of operational efficiency.
In order to characterize the impact of multiple collision points on beam-beam corrections to luminosity calibrations, COMBI was used to simulate vdM scans at one IP (e.g. IP1), with additional head-on collision points optionally inserted at one or more other IPs (in this example IP5, IP2 and/or IP8). The bunch parameters are listed in Table 5. The beam transport between consecutive IPs is modeled by linear maps that reflect the nominal phase advance in the LHC rings during the Run-2 vdM sessions at √ s = 13 TeV (  Table 5: Beam conditions assumed for the multi-IP simulations described in Sec. 4.6. The crossing angle θ c is set to 0 at all IPs so as to fully decouple the multi-IP study from potential longitudinal-dynamics effects. This Section is organized as follows. Multi-IP collisions increase the beam-beam tune spread, thereby shifting the means of the unperturbedtune spectra to lower values (Sec. 4.6.2). They also affect the peak-rate and beam-size bias factors in a manner that depends on the phase advance between consecutive collision points around the rings, but that leaves the σ vis bias almost totally insensitive to these phase-advance values. The analysis can be greatly simplified by an astute choice of the reference "no beam-beam" configuration (Sec. 4.6.3), up to a small residual ambiguity that depends on the phase advance between the non-scanning IP(s) and that where the scans are taking place (Sec. 4.6.4). This approach naturally leads to an effective parameterization of the impact of multi-IP collisions on beam-beam corrections to vdM calibrations (Sec. 4.6.5). Candidate beam-beam correction strategies in the presence of multi-IP collisions are summarized in Sec. 4.6.6.

Impact of multiple interaction points on the tune spectra
A comparison of simulated vertical-tune spectra for N NSIP = 0 and N NSIP = 1 is presented in Fig. 25. The single-IP beam-beam parameter ξ is chosen to be typical of vdM scans in pp collisions. For beams colliding at the scanning IP only (N NSIP = 0), two narrow peaks appear: one near the nominal tune (q x ≈ 0.320), and the other near q x − Y * ξ, where Y ≈ 1.1 is the Yokoya factor in the soft-Gaussian approximation [45]. These peaks correspond to, respectively, the so-called σ and π modes of coherent bunch oscillation. When beams collide in addition at one non-scanning IP (N NSIP = 1), the peak of the σ mode is unaffected, but the beam-beam tune spread roughly doubles and the π mode is shifted further down, near q x − 2Y * ξ. In both cases, the mean tunes lie roughly half-way between the σ and π peaks, and therefore differ by roughly −ξ/2 between the two configurations. Adding a second non-scanning IP results in a further widening of the tune spread, to roughly three times that for N NSIP = 0, and in a correspondingly larger downward shift of the mean tune.

Comparison of reference "no beam-beam" configurations in the presence of one non-scanning IP
In Secs. 2, 3 and 4.1 to 4.5, it is implicitly assumed that the beams collide at the scanning IP only. The beam-beam-induced peak-rate, beam-size and luminosity-bias factors defined in Sec. 4.2.4 are all expressed with respect to an unambiguous reference configuration that corresponds to "no beam-beam anywhere", in which the "nominal " reference variables (σ 0 , µ 0 pk , Σ 0 x , Σ 0 y , L 0 , σ 0 vis ,...) are identified by a "0" index. In the presence of collisions at non-scanning IPs, however, an ambiguity arises: should beambeam-induced biases be referenced to a configuration in which the beam-beam interaction is switched off at every single IP, or to one in which the beam-beam interaction is switched off at the scanning IP only, but remains active at the non-scanning IPs (if any)?
Confronting this question requires a more precise nomenclature than that used so far: • when the reference configuration is one in which the beam-beam interaction is switched off at all IPs, i.e. if vdM observables are corrected to "no beam-beam anywhere", the corresponding reference variables carry a "0" index as before. In particular, the reference luminosity L no−bb at the scanning IP is equal to the nominal luminosity defined in Sec. 2.3.2 (L no−bb = L 0 ). This is referred to below as the "L 0 normalization"; • when the reference configuration is one in which vdM observables are corrected to "no beambeam at the scanning IP only" even when the bunches considered also collide at one or more other IPs, the corresponding reference variables carry a "u" index instead. In particular, the reference luminosity L no−bb at the scanning IP is given by the unperturbed luminosity L u (L no−bb = L u ), i.e. that which would be measured at the scanning IP if the beam-beam interaction could be turned off at that IP only. This quantity is accessible only in simulations; its beam-separation dependence is in general notably different from that of L 0 . The same notation, and physical meaning, are applied to the other reference variables. This is referred to below as the "L u normalization"; • in the limiting case of collisions at the scanning IP only (N NSIP = 0), L u and L 0 are one and the same, as are all the other reference variables: there is no difference between "nominal" and "unperturbed" quantities.
The beam-separation dependence of the L u /L 0 ratio, during a vertical scan at IP1 and in the presence of head-on collisions at IP5, is presented in Fig. 26 (blue lozenges). Even though the unperturbed luminosity L u is by definition unaffected by beam-beam effects at the scanning IP, it exhibits a monotonic increase with beam separation when compared to the reference luminosity L 0 . This seemingly counter-intuitive observation is explained by the fact that the horizontal axis in Fig. 26 is the nominal transverse beam size σ 0 , not the actual beam size during the scan; equivalently, the actual relative separation, i.e. that normalized to the actual unperturbed beam size, differs from the nominal relative separation, which is normalized to the nominal beam size. The monotonic increase of the L u /L 0 ratio therefore suggests that in this particular configuration, the collisions at IP5 result in an enlarged vertical convolved beam size at IP1 (Σ u y > Σ 0 y ). This interpretation is confirmed by the evolution, during the scan, of the individual unperturbed single-beam sizes σ u yB (B = 1, 2) that is illustrated in Fig. 27 (open lozenges): both are larger than σ 0 , they remain constant as a function of the beam separation at the scanning IP, and their difference reflects the fact that the vertical phase-advance values from IP1 to IP5 differ significantly between the two beams ( Table 6, top two rows).
Also shown in Fig. 26 is the luminosity-bias factor [L/L 0 ] Full BB in the L 0 normalization (red squares), that displays an s-like shape superposed on the monotonic beam-separation dependence of the L u /L 0 curve. The ratio of these two curves: represents the luminosity-bias factor in the L u normalization in the presence of one non-scanning IP ; it is displayed as filled green circles, and is highly reminiscent of the black curve in Fig. 9.
This same luminosity-bias factor that is presented in Fig. 26 for N NSIP = 1, is compared in Fig. 28 to the equivalent quantity for N NSIP = 0. Up to a small systematic difference that represents the impact of the additional head-on collisions that occur at the non-scanning IP, the two curves are essentially identical. This observation suggests that the beam-separation dependence of the beam-beam-induced orbit-shift and optical-distortion effects is encapsulated in that of [L/L u ] Full BB , even in the presence of collisions at the non-scanning IP(s). This interpretation is confirmed by comparing the separation dependence of the actual single-beam transverse RMS sizes σ iB (i = x, y; B = 1, 2) with, and without, additional collisions at the non-scanning IP. When N NSIP = 1 (Fig. 27), σ y1 > σ y2 ; these two quantities exhibit a very similar beam-separation dependence with respect to the corresponding unperturbed beam size σ u yB . A direct comparison of the σ yB /σ u yB ratio with and without a non-scanning IP (Fig. 29) shows that the two configurations result in an almost identical beam-separation dependence, the shape of which closely parallels that computed [4] in the linear approximation using MADX, once the qualitative impact of amplitude detuning is taken into account.
The contrast in beam-separation dependence between the nominal (L 0 ) and the unperturbed (L u ) luminosity-bias curves apparent in Fig. 26, and the "universality" of the [L/L u ] luminositybias curves (Fig. 28), can be further characterized by comparing the separation dependence of the unperturbed luminosity L u , computed numerically as the overlap integral over the unperturbed macroparticle distributions, with that of its Gaussian-equivalent counterpart L u Geq . The latter is defined as 13 : 13 Some of the constants appearing in the complete formulas (Eqs. (5)-(10)) are irrelevant here, and omitted for the sake of clarity. and similarly for L u Geq (δ y ). It is calculated under the assumption that the unperturbed macroparticle distributions are Gaussian, and adequately described (for this particular purpose) by the transverse RMS single-beam sizes σ u iB (i = x, y; B = 1, 2): The comparison between L u and L u Geq is presented in Fig. 30. During the vertical scan (filled symbols), the unperturbed luminosity (already shown in Fig. 26) and its Gaussian-equivalent counterpart are almost exactly equal, except at large separation (δ y /σ 0 > 4); both increase with separation relative to L 0 . The same occurs during the horizontal scan (open symbols), except that here both curves drop with increasing separation, because the horizontal phase advances are such that the dynamic-β effect at IP5 translates into smaller horizontal beam sizes at IP1. As already documented [4] in the linear approximation, the dynamic-β effect at the non-scanning IP(s) propagates to the scanning IP, resulting in either larger or smaller unperturbed transverse beam-sizes depending on the phase advance.
The fundamental observable of the van der Meer method is the beam-separation dependence of the measured interaction rate. Intuitively therefore, one can argue that a separation-independent change in unperturbed transverse beam size at the scanning IP (Fig. 27, open symbols) should not bias the σ vis measurement, provided the shift in unperturbed tunes is properly accounted for. The same reasoning is offered in Ref. [4], where the argument is more straightforward: in the linear approximation (Sec. 3.1), additional collisions are tantamount to inserting a quadrupole-like perturbation at the non-scanning IP(s), that in turn results in a small tune shift accompanied by a β-beating wave around each ring. In the more general case treated here, the excellent agreement between L u and L u Geq over most of the scanning range (Fig. 30) indicates that the dominant effect is that of a beam-and plane-dependent shift of the effective β * value, and of the corresponding unperturbed RMS single-beam sizes. Non-Gaussian tails associated with the non-linearity of the beambeam force at the non-scanning IP(s) have but a very minor impact on the unperturbed transversedensity distributions at the scanning IP.
As a first illustration of the magnitude of multi-IP effects and of the quantitative impact, on beam-beam corrections, of the difference between the green circles and the red squares in Fig. 26, Table 7 compares the separation-integrated FoMs (Sec. 4.2.4) extracted from a simulated vdM-scan pair at IP1, with and without additional head-on collisions at IP5.
• In the L u normalization, and as expected from Fig. 28, the hierarchy of the FoMs is very similar with (N NSIP = 1) and without (N NSIP = 0) collisions at the non-scanning IP. The beambeam bias is systematically more negative for N NSIP = 1: this is qualitatively consistent with the combination of a decrease of the mean unperturbed tunes (Sec. 4.6.2), and of the tune-dependence of the luminosity-bias curves (Fig. 13). As a result, the magnitude of the overall beam-beam bias on σ vis (row 7, columns 3 and 4) increases substantially in the presence of a non-scanning IP.
• For N NSIP = 1, the hierarchy of the FoMs is very different in the L 0 and L u normalizations, because the dynamic-β effect at the non-scanning IP has a significant impact on the transverse unperturbed beam sizes at the scanning IP (Fig. 27). The fact that the Σ x (Σ y ) bias is more negative (positive) in the L 0 normalization reflects the fact that the unperturbed horizontal (vertical) beam size at the scanning IP is smaller (larger) than its nominal value (Fig. 30). • In contrast, the σ vis -bias values for N NSIP = 1 agree to better than 0.005% between the L 0 and L u normalizations, indicating that in practice, the two normalizations yield equivalent beam-beam corrections. This can be understood by examining the bottom half of Table 7, that summarizes the impact of collisions at IP5 on the unperturbed scan variables at IP1: Σ u x decreases, Σ u y increases by a little more, resulting in a small decrease of the unperturbed peak interaction rate µ u pk ; the net change in σ vis is at the 0.004% level. Physically, this expresses the fact that the changes in the unperturbed transverse convolved beam sizes and the unperturbed head-on luminosity are correlated such that the visible-cross-section remains invariant, at least in the configuration simulated here.

Sensitivity of beam-beam corrections to the phase advance between IPs
The results presented in the preceding paragraphs indicate that even though the separation dependence of the luminosity-bias factors are very different in the L 0 and L u normalizations (Fig. 26), these two reference configurations appear to yield almost identical beam-beam corrections to σ vis (Table 7). What remains to be quantified is how robust this conclusion is against large variations in the assumed phase advance between consecutive IPs, as well as in the presence of more than one non-scanning IP. The study reported in Sec. 4.6.3 assumes the phase-advance values between IP1 and IP5 that are listed in the top two rows of Table 6; these correspond to one point in four-dimensional phaseadvance space (µ B1 x , µ B1 y , µ B2 x , µ B2 y ). To address the first question above, the analysis detailed in Sec. 4.6.3 was repeated while varying the horizontal and vertical phase advance between IP1 and IP5 over most of the [0, 2π] range, while keeping the nominal tunes constant; in order to keep the parameter space manageable, the B1 and B2 phase-advance values were in addition constrained to be the same, separately in the horizontal and the vertical plane (µ B1 . vdM scans were simulated at IP1 over a grid in (µ x , µ y ) space, with the bunches colliding headon at IP5, and using the bunch parameters listed in Table 5.
Inspection of the simulated tune spectra reveals that their detailed shape depends on the phase advance between IPs; the same applies to the transverse luminosity distributions. The interpretation is that particles at different amplitudes are subject to different combinations of beambeam kicks at the scanning and the non-scanning IP, that depend on the phase advance between the two. Beam-beam-induced bunch-shape distortions at the non-scanning IP(s) propagate around the ring in a manner that depends on that phase advance, thereby making the detailed shape of the L u distribution at the scanning IP phase-advancedependent.
At each point in the (µ x , µ y ) space, the luminosity-bias curves and the corresponding separation-integrated FoMs were computed in both the L 0 and the L u normalization. The beambeam-induced σ vis bias, or equivalently the magnitude of the overall beam-beam correction to the vdM-based absolute luminosity scale, was found to depend slightly on the choice of the reference "no beam-beam" configuration, i.e. on whether the beam-beam correction is calculated in the L 0 or in the L u normalization. The difference is phaseadvance dependent, and exhibits a clear periodic structure (Fig. 31); a periodic structure is also observed in such two-dimensional maps of the luminosity-bias factor (not shown), with different patterns for [L/L 0 ] and [L/L u ]. These periodicities have been shown not to be of statistical or numerical origin. The interplay between beambeam effects at different IPs, and its sensitivity to the phase advance between IPs, was already observed at LEP-II, where it could strongly affect the instantaneous-luminosity balance between the four experimental detectors [46]. One tentative interpretation here is cross-talk between collision points: beam-beam kicks at one IP propagate to  Table 7: Fractional peak-rate, beam-size and σ vis beam-beam-induced bias on measured (rows 3-7) and unperturbed (rows 8-12) scan variables, for simulated vdM scans at IP1, with (columns 2-3) and without (column 4) head-on collisions at IP5, using the nominal bunch parameters shown in Table 5 and the phase-advance settings listed in Table 6.
the other IP(s) in a phase-dependent manner, thereby modifying very slightly the corresponding density distributions as well as the associated beam-beam distortions. These perturbations then propagate back to the original IP, resulting in a small modulation of the σ vis bias around the central value controlled by the unperturbed tunes. This and other potential explanations remain under study. Fortunately, these subtle effects have only a minor impact on beam-beam corrections to vdM calibrations. The largest amplitude of the L 0 -L u discrepancy over the entire plane amounts to 0.06% of σ vis , which is smaller than, or at most comparable to, the absolute accuracy of the B*B or COMBI simulations. The reason why in Table 7 the L 0 and L u normalizations yield essentially identical results, is that the Run-2 phase-advance between IP1 and IP5 happens to lie in a region of the µ x − µ y plane where the inconsistency is the smallest. In this case therefore, and for bunches that collide only in these two IPs, the difference in σ vis bias between the L 0 and L u normalizations is negligible.
To assess the impact of more than one nonscanning IP, the above phase scans were repeated for collision patterns with up to three nonscanning IPs. To deal with the increased dimensionality of the problem, the fractional phase advance from IP1 to IP5 was fixed to ∆µ x = ∆µ y = 0.405, a choice that in in Fig. 31 corresponds to the largest absolute difference in σ vis bias between the L u and L 0 normalizations. With that constraint, and with the overall tunes and the IP5-IP8 phase advances fixed at their nominal values, the horizontal and vertical fractional phase-advance values from IP1 to IP2 were scanned over the range [0, 2π]. The study was then repeated with the roles of IP2 and IP8 interchanged. In all these configurations, the difference in σ vis bias between the L u and L 0 normalizations displays the same periodic pattern as that in Fig. 31, and the largest absolute difference in σ vis bias between the two reference configurations does not exceed 0.15% (for ξ = 0.006). While strictly speaking, the ensemble of these scans does not constitute an exhaustive sampling of the 12dimensional phase-advance space (2 beams × 2 planes × 3 IP combinations), it suggests that the largest possible disagreement between the L 0 and L u normalizations amounts to 0.1-0.2% of σ vis when considering the full parameter space. Since the phase advance between LHC IPs can be measured to a few degrees (∼ 0.015 × 2π) per beam and per plane, the "operating point" in phaseadvance space is known to very good precision. In some cases therefore, the ambiguity on the σ vis  Fig. 31: Phase-advance dependence of the difference, between the L 0 and the L u normalization, of the beam-beam induced σ vis bias. The horizontal and vertical axes represent the betatron phase advance between IP1 and IP5, here assumed to be the same for B1 and B2. The color scale quantifies the fractional inconsistency in σ vis bias between the two reference configurations; it ranges from −0.06 to +0.05% of σ vis .
bias associated with the choice of reference configuration can be significantly reduced, as already demonstrated for bunches that collide in IP1 and IP5 only.

Parameterization of multi-IP effects
To gain further insight into multi-IP effects and account, in a pragmatic fashion, for their impact on beam-beam corrections to luminosity calibrations, COMBI was used to simulate scans over the full range of beam-beam parameter values typical of vdM sessions at √ s = 13 TeV, using the Run-2 phase-advance settings (Table 6), both without (N NSIP = 0) and with (1 ≤ N NSIP ≤ 3) head-on collisions at one or more non-scanning IP(s).
The results are summarized in Fig. 32 for three representative values of the single-IP beam-beam parameter ξ, using the L u normalization. Since both the orbit-shift and the optical-distortion bias are sensitive to the unperturbed-tune values (Sec. 4.2.2), the overall beam-beam bias on the visible cross-section associated with a given value of ξ is significantly affected by collisions at nonscanning IPs. Taking as an example the case of ξ = 6 × 10 −3 , the σ vis bias when colliding only at the scanning IP amounts to −0.55% (orange circle). It grows to −1.0% for one non-scanning IP (pale-brown lozenge), to −1.3% for two nonscanning IPs (medium-brown lozenge), and to a little over −1.5% for the worst case of N NSIP = 3 (darkest-brown lozenge). 14 As N NSIP increases by one unit, both the horizontal and the vertical mean unperturbed tune decrease by approximately ξ/2; as a consequence, the luminosity-bias curves shift downwards (Figs. 12, 13) and the σ vis bias becomes accordingly more negative. The impact of this "multi-IP tune shift" is comparable in magnitude to the σ vis bias for N NSIP = 0; it varies from one bunch pair to another in a manner that depends on the single-bunch parameters and on the value of N NSIP associated with eachcolliding bunch pair at the IP where the scans are taking place. Fortunately however, this effect lends itself to a simple parameterization.
The diagonal curves in Fig. 32 display the visible cross-section bias for scans simulated at IP1 in the absence of collisions at any non-scanning IP, but with both the horizontal and the vertical unperturbed tunes shifted by a quantity ∆Q mIP that scales with ξ. For a given value of ξ, the intersections of the corresponding diagonal curve (orange in the example above) with the horizontal lines of the same hue corresponding to N NSIP = 1, 2, ..., determine the magnitude of the multi-IP equivalent tune shift ∆Q mIP that needs to be applied simultaneously to the unperturbed horizontal and vertical tunes input to the COMBI simulation to mimic the effect of collisions at non-scanning IPs. The corresponding values of ∆Q mIP /ξ are indicated by the vertical dashed lines in the figure.
A more telling representation of the same results is offered in Fig. 33, which displays the multi-IP equivalent tune shift as a function of the number of non-scanning IPs. The results are consistent with a linear dependence that differs only slightly from the naive scaling law inspired by the qualitative evolution of the tune spectra. The purely empirical, two-parameter linear fit represented by the purple solid line: provides a slightly better parameterization than Eq. (41), and appears accurate at the level of ∆Q mIP ∼ ±0.05 ξ. Figure 33 demonstrates that for a given scanning IP and for all values of ξ, the normalized multi-IP equivalent tune shift scales linearly with the number of non-scanning IPs. To what extent this scaling is "universal", i.e. how strongly it depends on the choice of the IP where the scans are taking place, is addressed by Fig. 34. Irrespective of the scanning IP considered, the overall dependence of the σ vis bias on the number of scanning IPs is reasonably well represented by the naive scaling law of Eq. (41), demonstrating that the dominant effect is the downward shift of the mean unperturbed tunes. There are, however, significant differences between marker families, most notably between the red crosses (IP5) and the cyan squares (IP8): these reflect the variety of phase-advance patterns between the scanning IP  Tables 5 and 6 respectively. and the other collision points. Similarly, within a given family, there exist several bunch pairings between B1 and B2 that yield the same number of additional collisions, but that correspond to different phase-advance configurations. When IP5 is the scanning IP, for instance, the configuration N NSIP = 2 can be achieved with additional collisions at IP1 and at either IP2 or IP8; these two cases can therefore yield slightly different σ vis biases. To minimize modeling errors, therefore, it is advisable that scaling parameterizations such as that illustrated by Fig. 33 be developed for each scanning IP separately, preferably taking into account the actual collision pattern experienced by each bunch that collides at the scanning IP under consideration.
Pursuing the same strategy in the L 0 normalization meets with limited success (Fig. 35): for a given value of N NSIP , the spread across different collision patterns is noticeably larger than in the L u normalization, even when restricted to one given scanning IP. This observation reflects the same underlying physics as the periodic structure in Fig. 31; its implications are addressed below.

Beam-beam correction strategies in the presence of multi-IP collisions
In summary, collisions at non-scanning IPs manifest themselves in two main ways. The widening of the tune spectra and the resulting downward shift of the mean unperturbed tunes (Sec. 4.6.2) significantly increase the beam-beam-induced bias on the absolute luminosity scale (Fig. 34), mainly because of the reduced amplitude of the opticaldistortion bias (Sec. 4.2.2). This tune shift is accompanied by a beam-and plane-dependent βbeating wave around each ring [4], that translates into a sub-percent rescaling of the unperturbed transverse single-beam sizes at the scanning IP (Fig. 27). The latter causes significant distortions of the luminosity-bias curves when using the L 0 normalization (Fig. 26), that depend strongly on the phase advance between consecutive collision points. The shift in unperturbed tunes (Sec. 4.6.2) has, by far, the dominant impact on beam-beam corrections to vdM calibrations. It can be accounted for by simulating collisions simultaneously at all relevant IPs, with the input tunes set to their nominal value, and by correcting the luminosityscan curves to "no beam-beam anywhere". This L 0 normalization is the approach adopted in Ref. [8]. The alternative is to use the L u normalization. It is inspired by the fact that the fundamental observable of the vdM method is the beam-separation dependence of the measured interaction rate. What should matter, therefore, is the variation of the beam-beam force during the scan: underlying, separation-independent modifications of the unperturbed transverse-density distributions are expected not to affect the visible cross-section. In the L u approach, the luminosityscan curves extracted from the same multi-IP simulation as above are corrected to "no beambeam at the scanning IP only" (Sec. 4.6.3); the resulting luminosity-bias functions become much less sensitive to the presence of non-scanning IPs (Fig. 28).
The L 0 and L u normalizations yield beambeam corrections that, depending on the phase advance between consecutive IPs, either are indistinguishable (Table 7), or in the worst case differ by at most 0.1-0.2% of σ vis (Sec. 4.6.4). These observations demonstrate that no phase-advance configuration exists in which beam-beam induced β-beating would significantly reduce the beambeam bias on σ vis .
The results of the full multi-IP simulations can be very closely approximated by those of the much faster single-IP simulation, with as input the same value of the beam-beam parameter ξ, but with the unperturbed tunes shifted downwards by an amount ∆Q mIP that in the L u normalization and for a given scanning IP, scales with ξ and depends linearly on N NSIP (Fig. 33). The scaling parameters are not fully universal, in the sense that they depend slightly on the choice of scanning IP (Fig. 34), or more fundamentally on the phase advance between the scanning and non-scanning IPs.
Clarifying the origin of the periodic structures in Fig. 31, and explaining why beam-beam corrections appear less sensitive to phase-advance effects in the L u (Fig. 34) than in the L 0 (Fig. 35) normalization, are topics of active study that lie beyond the scope of the present report. An ambiguity remains, therefore, as to which of the two reference "no beam-beam" configurations (i.e. which of the L 0 or the L u normalization) provides the most accurate beam-beam correction in absolute terms. Until this can be resolved, it seems prudent to keep the question open and to cover the difference by assigning a phase-dependence-related systematic uncertainty to the absolute magnitude of the beam-beam corrections.
On the practical side, the L u normalization offers several advantages. First, by reducing the multi-IP problem to that of a single collision point, it greatly simplifies the interpretation of simulation results and the estimate of systematic uncertainties. The L u approach also automatically accounts for the static (in the sense of separationindependent) change in unperturbed beam size associated with the dynamic-β effect at the nonscanning IPs, preventing it from blurring the distinction between orbit-shift and optical-distortion corrections. Finally, combining the parameterization of single-IP beam-beam biases offered in Sec. 4.2.3 with the scaling prescriptions of Sec. 4.6.5 provides a beam-beam-correction procedure that is straightforward to implement, requires no or little investment by individual users in complex, CPU-intensive simulations, and is in many cases sufficient to calculate absolute bunchby-bunch corrections to vdM calibrations, or to estimate systematic uncertainties. This is particularly valuable when dealing with five to ten vdM-scan pairs that involve up to 150 colliding bunches with different populations and emittances, arranged in multiple collision patterns that each require a different multi-IP correction [40]: the computing overhead that would be associated with simulating each of those configurations in detail is considerable.
The scaling prescriptions offered in Sec. 4.6.5 in the L u normalization account for collisions at non-scanning IPs to better than 0.1-0.2% on σ vis . It should be stressed, however, that these simple-minded models reproduce only the σ vis bias extracted from the multi-IP simulation (to an accuracy limited by the phase-advance-dependent inconsistency between the L 0 and L u normalizations), but that they are not meant to reproduce the corresponding tune spectra. An additional limitation is that they assume that collisions at non-scanning IPs occur with zero transverse separation, and with zero crossing angle. These ad hoc recipes are therefore no substitute for full-fledged 46 simulations that account in detail for every distinct multi-IP collision pattern, separately for each scanning IP. For many purposes, however, these simple scaling laws appear sufficiently accurate to cover most beam conditions encountered during Run-2 vdM sessions at the LHC. The potential interplay between multi-IP effects and vdM scans in the presence of a large crossing angle, however, may require a dedicated study.

Impact of beam-beam asymmetry
The results presented so far in this chapter, with the exception of those in Sec. 4.4, assume beams that are round in the transverse plane (β * xB = β * yB = β * B and σ xB = σ yB = σ B , where B = 1, 2). In such a case, the distinction between horizontal and vertical beam-beam parameters disappears, and the beam-beam parameter of beam 2, defined in Eq. (12), simplifies to where σ 1 is the transverse RMS size of beam 1. The same formula, mutatis mutandis, is used to define ξ 1 . The parameterization described in Sec. 4.2, and the results derived therefrom, further assume that the beam-beam forces are exactly balanced between the two beams, i.e. that ξ 1 = ξ 2 = ξ R (Eq. (28)). Beam-beam asymmetry occurs when this assumption is violated, i.e. when ξ 1 ̸ = ξ 2 , which can arise under three different scenarios (Sec. 4.7.1). Even though a parametric approach can provide physically intuitive insight (Sec. 4.7.2), dedicated simulations are required to quantify the impact of such asymmetries on the accuracy of the beam-beam correction procedure (Sec. 4.7.3).

Beam-beam asymmetry scenarios
Consider first the case of a bunch-population asymmetry (n 1 ̸ = n 2 , with β * 1 = β * 2 and σ 1 = σ 2 ), which can be simulated in B*B by assigning different populations to the source beam and to the witness beam. Since the source beam in B*B remains unaffected by the beam-beam interaction, and since the two beams are subject to different beam-beam parameters, the simulation must be run twice, first with n 1 and n 2 assigned to, respectively, the source and the witness beam, and then with the inverse assignment. The predicted luminosity-bias curves, and the corresponding corrections to the FoMs, are then given by adding the two sets of results 15 , implicitly assuming that strong-strong effects do not play a significant role even in the presence of a large bunch-population asymmetry (at least in the vdM regime). This procedure was validated by confronting the B*B results with those from COMBI for several asymmetric configurations: even in the worst case considered (n 2 /n 1 = 2, ξ 1 = 5.30 × 10 −3 , ξ 2 = 2.65 × 10 −3 ), the beam-beam corrections calculated by the two packages agreed to better than one part in 1600.
The B*B and COMBI results above are accurately reproduced by the parameterization detailed in Sec. 4.2.3, with as input the beamaveraged beam-beam parameter: This is a consequence of the fact that the ξdependence of the FoMs is close enough to linear (Figs. 14-16) for the averaging of the FoM values at two different ξ settings to be tantamount to computing the FoM once, with ⟨ξ⟩ as input.
The second scenario is that of a β-function asymmetry (β * 1 ̸ = β * 2 , with n 1 = n 2 and σ 1 = σ 2 ). This case is equivalent to that of the bunchpopulation asymmetry above, since only the product n k β * l (k, l = 1, 2 with l ̸ = k) plays a role in Eq. (42). The bunch population n k of beam k controls the scale of the beam-beam force beam k exerts on beam l, while β * l determines how sensitive beam l is to this force.
In contrast to the above two scenarios, which involve only scale effects, a beam-size or emittance asymmetry (σ 1 ̸ = σ 2 , with n 1 = n 2 = n and β * 1 = β * 2 = β * ) modifies in different ways the spatial distributions of the source field and of the particles in the witness bunch, the trajectories of which are affected by this field.

Parametric characterization of beam-beam asymmetries
A full characterization of beam-beam asymmetries requires the two parameters ξ 1 and ξ 2 . A physically intuitive parameterization can however be provided by a single, "effective" beam-beam parameterξ, defined as on the basis of the following argument.
Under the naive assumption that the optical distortion of B1 induced by the field of B2 is proportional to ξ 1 , and that it can be adequately represented by a small fractional change λ 1 in the value of σ 1 : where λ 1 ∝ ξ 1 , one can show, using Eqs. (5) and (6) and after a bit of algebra, that to first order in ξ 1 , the luminosity-bias factor [L/L 0 ] Opt,1 associated with the optical distortion of B1 alone obeys Applying the same argument to the optical distortion of B2 induced by the field of B1, and combining the impact of the optical distortions of the two beams, it follows that In the limit of equal beam sizes (σ 1 = σ 2 ),ξ reduces to ⟨ξ⟩, and as such properly characterizes bunch-population or β * asymmetries. If instead the bunch populations are balanced (n 1 = n 2 = n) but the transverse beam sizes differ (σ 2 ̸ = σ 1 ), then for a given value of Σ R = σ 2 1 + σ 2 2 , the effective beam-beam parameter can be written as thereby incorporating the effect of beam-size asymmetries in a single-parameter, but approximate, characterization of the strength of the beam-beam interaction. In the limit of σ 2 /σ 1 = 1,ξ reduces to the symmetric, round-beam equivalent beam-beam parameter ξ R (Eq. (28)).

Predicted impact of transverse beam-size asymmetries
Since the beam-beam induced orbit shift depends only on the convolved beam size Σ R (Sec. 2.3.1), but not on the σ 2 /σ 1 ratio, only the beam-beaminduced optical distortion may be affected by a beam-size asymmetry. It is natural, therefore, to characterize this effect by varying, in the simulation, the common bunch intensity and the beam-size ratio, while keeping Σ R constant. When using COMBI, this can be done in a single step for each parameter set, since the spatial density distributions of both bunches evolve simultaneously. When using B*B, however, the source bunch remains unaffected. Therefore the simulation must be run twice, first with σ 1 and σ 2 assigned to, respectively, the source and the witness beam, and then with the inverse assignment, after which the contributions of the two passes are combined. Such a study is presented in Fig. 36, for a highly asymmetric example in which ξ 1 /ξ 2 = 2 andξ = 4.42 × 10 −3 . When the wider of the two bunches (ξ 1 = 5.30 × 10 −3 ) plays the role of the witness beam (triangles and squares), the orbitshift and optical-distortion contributions are of opposite sign and of comparable magnitude; when the narrower bunch (ξ 2 = 2.65 × 10 −3 ) becomes the witness beam (dashed curves), the beam-beam parameter drops by a factor of two, the orbit-shift contribution remains unchanged as it should, but the optical-distortion effect becomes almost negligible. The combination of the two contributions, obtained by multiplying the luminosity-distortion factors associated with the two witness beams, is presented in Fig. 37. The two beams contribute equally to the orbit effect, but the opticaldistortion is dominated by that experienced by the wider (i.e. the weaker) of the two beams.
To quantify the impact of a transverse beamsize asymmetry on beam-beam corrections to vdM calibrations, the simulation protocol detailed above was repeated for several configurations in which the bunch populations were symmetric (n 1 = n 2 = n), the beam-size ratio σ 2 /σ 1 was varied from 1.0 down to 0.65, and the convolved beam size was kept constant at Σ R = 122 µm. Two values of the common bunch population were considered (n = 0.78 × 10 11 and n = 1.21 × 10 11 p/bunch), near the lower end, and beyond the upper end, of the bunch-intensity range used during vdM scans; these translate into roundbeam equivalent beam-beam parameter values of, respectively ξ R = 3.53×10 −3 and 5.50×10 −3 . The visible cross-section bias σ vis / σ 0 vis − 1 corresponding to each of these configurations, and computed from luminosity-bias curves such as those displayed in Fig. 37, is presented in Figs. 38 and 39. As the beam-size ratio deviates more and more from unity,ξ grows (see Eq. (45)), and so does the σ vis bias. Its dependence onξ is almost linear (Fig. 38), with only a small quadratic component. The scaling law, however, is approximate only, as evidenced by the slight relative misalignment of the B*B curves at the interface between the two groups (ξ R ∼ 5.2). The dependence on the beam-size ratio (Fig. 39) is well modeled by a third-order polynomial of σ 2 /σ 1 , the coefficients of which depend on the n β * product. A systematic difference between COMBI ands B*B becomes apparent at large asymmetry, suggesting that coherent oscillations may start to play a role in that regime. Quantitatively however, this discrepancy is small enough in the asymmetry range of interest at LHC (at most 0.01% on σ vis for σ 2 /σ 1 > 0.9, i.e. for an emittance imbalance of 20%) that it can be absorbed in a systematic uncertainty.
The impact of the beam-size asymmetry on the absolute luminosity scale, i.e. the difference in the σ vis bias between a given asymmetric configuration (in this example σ 2 /σ 1 = 0.8) and the fully symmetric case (σ 1 = σ 2 ), is indicated by the arrows. The same results are presented in Fig. 40 (top three curves), on an expanded scale and after subtracting the σ vis bias corresponding to the symmetric configuration.
The results displayed in Figs. 36-39 assume that the beams collide at the scanning IP only. Additional collisions at other IPs lower the unperturbed tunes (Sec. 4.6.2): as a result, the impact of the beam-beam-induced optical distortion is significantly reduced (Fig. 12), and so is the impact of the transverse beam-size asymmetry. This is illustrated by the two bottom curves in Fig. 40, for the case of bunches colliding at two non-scanning IPs (N NSIP = 2), that was modeled using B*B to simulate single-IP collisions but with the unperturbed tunes shifted as prescribed by Eq. (41). At the largest asymmetry considered (σ 2 /σ 1 = 0.65), the sensitivity of the σ vis bias to the beam-size asymmetry drops by more than a factor of two. In addition, the bunch-current dependence of this bias, manifested by the difference between the filled and the open circles, is several times smaller for the N NSIP = 2 case. This occurs because the shift ∆Q mIP in unperturbed tunes is proportional to both ξ and N NSIP : when both are large enough, Dependence on the transverse beamsize ratio σ 2 /σ 1 of the visible cross-section bias σ vis / σ 0 vis − 1 associated with optical distortions, as predicted by B*B (diamonds) or COMBI (triangles) simulations, for low-intensity (red and blue bottom curves) and high-intensity (green top curve) colliding-bunch configurations that cover a wide range of beam-size asymmetries at constant Σ R (see text). The beams collide at the scanning IP only. The σ vis bias in the fully symmetric configuration (σ 1 = σ 2 ) is indicated by a horizontal black line in each group. The curves are third-order polynomial fits to the points, and are intended to guide the eye. The arrows indicate the impact on the absolute luminosity scale of a 20% beam-size imbalance (σ 2 /σ 1 = 0.8).
this tune shift more than compensates the increase in optical distortion at the scanning IP associated with the increased beam-beam parameter experienced by the weaker of the two beams.
During Run-2 pp vdM scans at √ s = 13 TeV, the bunch population n p rarely exceeded 10 11 (Table 1), and the typical beam-size asymmetry σ 2 /σ 1 hovered around 0.94, as extracted from one of the non-factorization analyses outlined in Sec. 4.3.1. In practice therefore, for these vdM sessions the impact of beam-beam asymmetries on beam-beam corrections to vdM scans does not exceed 0.02% (Fig. 40); such a small number can be safely accounted for as a systematic uncertainty in the absolute luminosity scale.  scans at the LHC. When an uncertainty is assumptiondependent, the value flagged by an asterisk is that used in computing the total uncertainty; the latter is compared to the overall beam-beam correction itself in the bottom two rows of the Table. The rightmost column indicates the chapter(s) where the corresponding issues are discussed in detail.  Fig. 40: Impact of the beam-size asymmetry on the visible cross-section bias σ vis / σ 0 vis − 1 associated with optical distortions, as simulated using either B*B (diamonds and circles) or COMBI (triangles). Shown is the change is σ vis bias between an asymmetric beam-beam configuration characterized by a given setting of σ 2 /σ 1 < 1, and the corresponding symmetric configuration (σ 2 = σ 1 ). The top three curves (green diamonds, red triangles and blue diamonds) correspond to the data displayed in Fig. 39, when the beams collide at the scanning IP only. The bottom two curves (circles) show the same results when the beams also collide head-on at two non-scanning IPs. The open (filled) symbols correspond to the low-(high-) intensity bunches (see text). The curves are third-order polynomial fits to the points, and are intended to guide the eye.

Systematic uncertainties
The sources of systematic uncertainty that may affect the absolute magnitude of beam-beam corrections during vdM scans, and that are detailed in Sec. 4, can be regrouped in three categories: optical configuration of the LHC rings, deviations from the fully symmetric Gaussian-beam configuration, and collisions at multiple IPs. The associated uncertainties are briefly elaborated upon in Secs. 5.1, 5.2 and 5.3 below; additional sources of potential bias are addressed in Sec. 5.4. All the uncertainties are then added in quadrature, and their combined magnitude compared in Sec. 5.5 to that of the overall beam-beam correction itself.
The quantitative impact of each uncertainty on the σ vis bias, or equivalently on the magnitude of the beam-beam correction, is summarized in Table 8 for scans simulated using three different collision patterns (N NSIP = 0, 1, 2). For each uncertainty source, the beam conditions input to the simulation are detailed in the chapters listed in the last column of the table. Not all simulations, however, use exactly the same beam conditions. In order to provide a self-consistent picture of the relative magnitude of the various uncertainties, and since the magnitude of the beam-beam bias, as well as most of the uncertainties, scale with ξ, the uncertainties listed in in Table 8 are all expressed at a common value ξ sim of the beam-beam parameter; the details of this procedure can be found either in the table itself, or in Secs. 5.1 to 5.3. Because the chosen value of ξ sim lies at the upper end of the range covered during Run-2 vdM scans, Table 8 illustrates something of a "worst-case" scenario with respect to the magnitude of the beam-beam correction and of the associated systematic uncertainty. The nmerical results listed in that table should therefore not be applied blindly to estimating actual vdM-calibration uncertainties, but rather used as guidance for a case-by-case and bunch-by-bunch error analysis of actual vdM-calibration data.

β * uncertainty at the scanning IP
Since the σ vis bias scales with the beam-beam parameter (Sec. 4.2), any measurement error on ξ directly translates into an error of similar relative magnitude on beam-beam corrections to the luminosity scale. Existing accelerator instrumentation unfortunately does not allow accurate enough a determination of the transverse single-beam sizes that enter the definition of ξ (Eq. (12)). The horizontal and vertical beam-beam parameters must therefore be redefined in terms of two-beam observables that can be reliably measured: Here, β * xAv = (β * x1 + β * x2 )/2 is the beam-averaged horizontal β-function at the scanning IP, and similarly for β * yAv . The quantity n = (n 1 +n 2 )/2 refers to the beam-averaged bunch population associated with the colliding-bunch pair under study, γ = γ 1 = γ 2 is the common relativistic factor (assumed here to be the same for B1 and B2), and Σ x and Σ y are the convolved transverse bunch sizes. During vdM scans, the observables n, Σ x and Σ y are measured on a bunch-by-bunch basis to sub-percent accuracy [38,39,40,41], and the uncertainty on the absolute beam energy is totally negligible for the purposes of this paper. The main uncertainty affecting the absolute scale of the beam-beam parameter, therefore, is that associated with potential deviations of β * from its nominal or assumed value. The facts that the ξ dependence of the σ vis bias is roughly linear (Fig. 16), and that the measured value of ξ xj (Eq. (12)) depends linearly on n k and β * xj (j, k = 1, 2; j ̸ = k), justify the use of the beam-averaged variables defined in Eqs. (46)-(47) above. Furthermore, the very weak dependence of the σ vis bias on the Σ 0 y /Σ 0 x aspect ratio (Fig. 22) indicates that "averaging" over x and y, in the sense of substituting the round-beam equivalent parameter ξ R (Eq. (28)) for the pair of horizontal and vertical beam-beam parameters (ξ xAv , ξ yAv ), results in an acceptably small bias. What finally matters, therefore, is the uncertainty on the beam-and plane-averaged value of β * at the scanning IP.
The β * values at the LHC IPs are extracted from optical-function measurements that are carried out for every colliding-beam configuration as part of routine accelerator commissioning. During the first few years of LHC operation, these values were typically obtained using the phase-advance method; k-modulation is now considered a more promising approach [47]. It targets an absolute accuracy of a few percent on β * iB (i = x, y; B = 1, 2), provided sufficient beam time can be dedicated to the measurement.
The absolute β * -accuracy requirements associated with beam-beam-correction uncertainties are rather modest. When beams collide at the scanning IP only, and for ξ = 5.6 · 10 −3 , a ±10% uncertainty in β * at the scanning IP, uncorrelated between beams (because the optics of the two rings are almost completely independent) but correlated between planes (because a focusing error presumably affects both the horizontal and the vertical β-function) results in ±0.06% uncertainty in σ vis (Table 8). Treating this uncertainty as fully uncorrelated (correlated) between beams and between planes would decrease (increase) the corresponding σ vis uncertainty by a factor of √ 2. The absolute uncertainty increases to 0.10% (0.13%) for N NSIP = 1(2), but its magnitude relative to the overall beam-beam correction remains about the same.
The ±10% β * -uncertainty assumed in the example above is very conservative, and is representative of the typical difference between the nominal, i.e. the dialed-in, β * setting and the measured values. Smaller uncertainties can be achieved by calculating ξ based on the measured β * values themselves, while accounting rigorously for the systematic β * -measurement uncertainties and for their correlation between beams and between planes. Figures 12 and 13 illustrate the sensitivity of the optical distortion and of the full beam-beam bias to the unperturbed fractional tunes. Any deviation of (q x , q y ) from their nominal setting results in an inaccurate prediction of the σ vis bias.

Nominal collision tunes
For the purposes of the present discussion, a conservative uncertainty of ±0.002 is assigned to each of q x and q y , and is assumed to be correlated between beams (as could be caused by, for instance, instrumental effects) and between planes (because a focusing error presumably affects both the horizontal and the vertical betatron frequency). The magnitude of this uncertainty is based on LHC operational experience 16 , and is meant to cover setting errors, as well as the stability of the unperturbed tunes, over the course of the vdM-calibration fill(s). The corresponding uncertainty in σ vis (Table 8) dominates the overall beam-beam uncertainty. It decreases slightly when the number of non-scanning IPs increases, because the additional collisions lower the unperturbed tunes, thereby reducing the sensitivity of the optical distortions to the exact values of the fractional tunes. The same effect is responsible for the N NSIP -dependent sensitivity of σ vis to beam-beam asymmetries (Fig. 40).

Deviations from the fully symmetric Gaussian-beam configuration
Except where specifically stated otherwise, most of the results presented in this paper, and in particular the parameterizations detailed in Secs. 4.2.3 and 4.6.5, are based on the assumptions listed in Sec. 4.1. Each of these will be lifted, one at a time, in the following Sections, and their individual impact on beam-beam-correction uncertainties estimated using the studies detailed in Secs. 4.3 to 4.7.

Non-Gaussian unperturbed transverse-density profiles
Of the three single-bunch models described in Table 2, the 2012 parameter set represents an extreme case, both because this is the vdM session that revealed the largest non-factorization biases at LHC ever [13,40], and because the parameters in that Since the factorizable double-Gaussian configuration systematically yields a larger deviation from the corresponding single-Gaussian model than its non-factorizable counterpart (Table 4), it is favored for the purpose of estimating systematic errors. The uncertainty listed in Table 8 for N NSIP = 0 is therefore taken from the top half of Table 4. It is extrapolated to the cases of one and two non-scanning IPs by assuming that the fractional impact of additional collisions on the σ vis bias is dominated by the multi-IP tune shift (Sec. 4.6.2), and thus insensitive to the shape of the unperturbed transverse-density distributions.
The choice of the 2018 factorizable model to estimate the impact of non-Gaussian tails is not devoid of some arbitrariness. The uncertainties listed in Table 8, therefore, should be considered as a first indication, albeit a prudent one. A refined evaluation would require a more extensive characterization of the transverse-density distributions extracted from non-factorization analyses, that could then be used to produce a wider palette of more realistic single-bunch models for input to B*B or COMBI. In order to explore potential correlations between systematic biases, one could in addition extend the existing studies to simulate bunches that are not purely Gaussian and that in addition collide at multiple interaction points.

Beam ellipticity at the scanning IP
The dependence of the optical-distortion bias on the Σ y /Σ x aspect ratio is summarized by Fig. 22, and does not exceed ±0.025% for ξ R ≤ 4.2 · 10 −3 . Scaling this linearly to the value of ξ sim shown in Table 8 yields a beam-beam-correction uncertainty of 0.03%, under the condition that the orbit-shift correction take the actual ellipticity into account, i.e. that it be computed using the Bassetti-Erskine formalism [6] with, as input, the measured values of Σ x and Σ y .

Non-zero crossing angle
The crossing-angle dependence of the beam-beam bias in the vdM regime is summarized by Fig. 23. The associated systematic uncertainty is taken as the difference in σ vis bias between θ c = 0 and the actual crossing-angle value; it is listed in Table 8 for two values of the full crossing angle. It should be emphasized that this small uncertainty is applicable only at moderate crossing angles (θ c ≤ 150 µrad), and only in the vdM regime (ξ ≤ 0.006, β * ≥ 10 m). For larger crossing angles in the vdM regime, additional effects may come into play. This may also be the case for scans performed in the physics regime, even at moderate crossing angle, when the beam-beam parameter becomes large enough, and/or when β * becomes small enough ( β * ≤ 30 cm) for the hourglass effect to no longer be negligible. Of particular concern are the so-called "emittance scans", that are vdM-like beam-separation scans that are used to monitor luminometer performance and transverse emittances during routine LHC operation [48]: the study of their sensitivity to beam-beam-induced biases remains mostly unexplored territory.

Beam-beam imbalance
Of the three beam-beam-asymmetry scenarios considered in Sec. 4.7, only the transverseemittance imbalance (Sec. 4.7.3) requires special attention, since bunch-current and β *asymmetries are accounted for by the use of the beam-averaged beam-beam parameters defined in Eqs. (46)- (47). The emittance-asymmetry dependence of the σ vis bias is controlled by the combination of ξ R and of the transverse beam-size ratio σ 2 /σ 1 (Eq. (45)); it also is sensitive to the number of non-scanning IPs (Fig. 40).
The beam-beam-imbalance uncertainty is taken as the difference in σ vis bias between an asymmetric beam-beam configuration characterized by a given value of σ 2 /σ 1 < 1, and the corresponding symmetric configuration (σ 2 = σ 1 ); it is listed in Table 8 for three values of the transverse beam-size ratio. The ratio closest to unity (σ 2 /σ 1 > 0.95), that emerged as the most representative of typical beam conditions during the 2015-2018 pp vdM sessions [40], is the one used when computing the total beam-beam uncertainty (penultimate row of Table 8). The N NSIP dependence is discussed in Sec. 4.7.3, and is closely related to the tune-dependence of the optical-distortion bias (Sec. 5.1.2).

Multiple interaction points
Two sources of uncertainty affect the beam-beam correction strategies in the presence of multi-IP collisions, that are detailed in Sec. 4.6.6: the phase advance between consecutive IPs (Sec. 4.6.4), and the parameterization of the multi-IP equivalent tune shift ∆Q mIP (Sec. 4.6.5). These uncertainties affect only those colliding-bunch pairs in which at least one of the two opposing bunches experiences collisions at one or more IPs other than the scanning IP (N NSIP > 0).
The phase-advance uncertainty arises from the currently unresolved ambiguity in the choice of the reference "no beam-beam" configuration, or, in more technical terms, from the ambiguity between the L 0 and the L u normalizations (Sec. 4.6.3). In the absence of any phase-advance information, e.g. for a bunch pattern that mixes different collision configurations that are all treated on the same footing, a blanket systematic uncertainty of at most 0.2% appears adequate (Sec. 4.6.4). This uncertainty can, in some cases, be significantly lowered by simulating the actual collision pattern associated with each colliding-bunch pair, with as input the actual IP-to-IP phase advance per beam and per plane, and by then comparing the results between the two reference no-beam configurations. An illustration of the potential gain provided by such a refined analysis is offered by the results presented in Table 7 for the N NSIP = 1 case, a configuration that applies to the 13 TeV pp vdM scans at IP1 and IP5 in 2015, 2016 and 2017 [40].
The uncertainty associated with the equivalent multi-IP tune shift ∆Q mIP needs to be considered only when invoking the parameterization of single-IP simulations proposed in Sec. 4.6.5 and encapsulated in Eq. (41), to compute beam-beam corrections in the presence of additional collisions. This parametrization, that is illustrated by the magenta curve in Fig. 34, is in good agreement with the simulation results when IP5 is the scanning IP. The discrepancies that are apparent for scans at some of the other IPs, reflect the variety of phase-advance patterns between the scanning IP and the other collision points. Varying the p 1 parameter in Eq. (41) by ±15% of its central value covers all the collision patterns shown in Fig. 34, thereby providing a measure of the uncertainty associated with this simplified procedure (Table 8). More precise results, and a smaller uncertainty, could again be obtained by simulating separately each of the collision patterns associated with the B1 and B2 bunch strings circulating in the LHC ring during the vdM-calibration session under analysis.

Long-range encounters
During pp vdM scans at the LHC, the bunches are longitudinally isolated, in the sense that consecutive bunches are separated by more than 500 ns, rather than grouped in trains with only 25 ns spacing. In addition, the fill pattern is systematically optimized to avoid parasitic crossings, at least at the scanning IP. This strategy eliminates separation-dependent, long-range beambeam kicks that might distort the luminosity-scan curves. 17 In contrast, during high-luminosity pp physics running, consecutive bunches are typically separated longitudinally by only 25 ns, the emittance is smaller than during vdM scans, and the bunch current significantly larger. As a result, and in spite of the non-zero crossing angle, long-range beambeam effects are known to significantly impact emittance scans in several ways, most of which remain to be understood.

Lattice non-linearities
The strongest non-linearity is introduced by the beam-beam interaction itself. However, non-linear magnetic elements in the LHC rings also affect beam dynamics, potentially leading to additional distortions of the transverse particle-density distributions. Because the LHC relies on the combination of high chromaticity (Q ′ ≈ 10 − 15 units) and of high Landau-octupole currents to mitigate some strong impedance-driven instabilities, the next largest optical non-linearities at beam energies above 1 TeV are associated with the corresponding sextupoles and octupoles [49]. As demonstrated experimentally [50,51,52], machine imperfections contribute even less by comparison, and are therefore neglected in the present study.
In order to quantify the impact of these magnets on the σ vis bias, both the chromaticity and the linear detuning with amplitude introduced by the Landau octupoles were added to the COMBI model. The chromaticity dependence of the beam-beam bias is presented in Fig. 41 for two different beam energies. As expected, the impact of the octupoles is more important in the 2011 case because of the lower beam energy: at zero chromaticity, powering the octupoles changes the beam-beam bias from −0.155% to −0.130%, a reduction in relative magnitude of about one sixth. In contrast, in the 2018 configuration the effect of the octupoles is three times smaller, and almost negligible. As for the chromaticity, its impact is negligible in both configurations.
The very low sensitivity of beam-beam corrections to these non-linear elements can be understood as follows. Octupoles act more strongly mode with a non-zero nominal crossing angle. However, the bunch charge is considerably smaller than during pp vdM scans, and the beam-beam parameter about an order of magnitude smaller. Parasitic crossings, therefore, are no issue from the viewpoint of beam-beam corrections to HI vdM calibrations. on larger-amplitude particles. These barely contribute to the luminosity when the beams collide head-on; the fractional contribution of the tails to the overlap integral is highest at large beam separation, but there the luminosity is the lowest. This is why modifying the trajectory of the tail particles barely influences the separation-integrated FoMs, and therefore the σ vis bias. Similarly, high chromaticity potentially enhances synchro-betatron coupling; however since vdM scans are carried out with relaxed β * settings and at zero crossing angle, this is of no concern.
Even though lattice non-linearities are almost negligible in vdM scans (Table 8), it remains important to check systematically the potential influence of the octupoles under operational scenarios where they are pushed to maximum power to ensure beam stability: their impact on beambeam uncertainties may turn out to be no longer negligible, for instance during emittance scans.

Numerical accuracy of polynomial parameterization
Propagating to the σ vis bias the parameterization uncertainties that affect the separation-dependent luminosity-bias factors (Eq. (27)) and that are discussed in Sec. 4.2.3, results in an upper limit of 0.1% on the associated beam-beam uncertainty. This contribution should be ignored when resorting to a full B*B or COMBI simulation at each scan point, rather than to the parameterization.

Summary of systematic uncertainties
Adding in quadrature all the above uncertainties yields, for the example examined here (ξ sim = 5.6 · 10 −3 ), a total beam-beam uncertainty of 0.32 to 0.46%, depending on the number of nonscanning IPs ( Table 8). The corresponding value of the beam-beam correction varies from 0.52 to 1.17%. The absolute magnitude of the total uncertainty grows when N NSIP increases from 0 to 2, but its magnitude relative to the correction itself decreases from about 60% to roughly 40% through the mechanism dicussed in Secs. 5.1.2 and 5.2.4. Under operating conditions typical of actual vdM-calibration sessions, the number of nonscanning IPs is at least one. For N NSIP = 1, the four dominant uncertainty contributions, in order of decreasing importance, are associated with the nominal tunes, non-Gaussian tails, multi-IP effects (in the absence of collision-patternspecific studies), and the β * uncertainty. Estimating the uncertainty associated with non-Gaussian tails, and verifying whether beam-size asymmetries might contribute at a significant level, imperatively requires input from the non-factorization analysis of luminous-region data.
It is worth noting that the larger N NSIP , the larger the contribution of the orbit-shift correction relative to that of the optical distortion. Qualitatively therefore, estimating precisely the effects that contribute only to the opticaldistortion uncertainty, viz. the beam ellipticity and the beam-size asymmetry, becomes less critical as the number of non-scanning IPs increases. In contrast, the β * and tune uncertainties affect the magnitude of both the orbit shift and the optical distortion, in a fully correlated manner.
The uncertainties detailed in Table 8 have been deliberately estimated in a conservative fashion, and for the largest value of the beam-beam parameter encountered during Run-2 vdM scans in pp collisions at √ s = 13 TeV. Most uncertainties scale with ξ, either linearly or with a slightly positive quadratic component. For the purposes of a rough estimate, therefore, they can be safely extrapolated linearly to smaller ξ values, that are more representative of typical vdM conditions. A more precise determination of these uncertainties requires a case-by-case analysis, possibly augmented by some of the additional simulation studies suggested in some of the preceding sections in this chapter.

Summary and outlook
Under conditions typical of vdM luminositycalibration scans in pp collisions at the LHC, beam-beam-induced orbit shifts and optical distortions, if left uncorrected, bias the absolute luminosity scale by an amount comparable to, or larger than, the total luminosity-uncertainty budget. The magnitude of the correction, which ranges from approximately 0.2% to 1.2%, mainly depends on the beam-beam parameter ξ and on the number of non-scanning IPs. The contributions to the associated systematic uncertainty are listed in Table 8 and elaborated upon in Sec. 5.5. The total uncertainty amounts to roughly half of the correction itself, and is dominated by tune-related effects, either directly (accuracy and reproducibility of nominal-tune settings) or indirectly (beam-beam tune shift at non-scanning IPs, phase advance between consecutive IPs). The next most important source of uncertainty is the potential deviation of the unperturbed transverse bunch-density distributions from a perfectly Gaussian shape. Simulating the impact of non-Gaussian tails cannot be based on luminosity data alone: it requires, as input, single-bunch parameters that can only be extracted from the measured beam-separation dependence, during the scan, not only of the luminosity but also of the position, shape and orientation of the luminous region. These singlebunch parameters are also crucial to quantify transverse beam-size differences between the two beams, which, if large enough, may become a significant source of uncertainty.
The scope of the present paper is limited to beam-separation scans in the vdM regime (ξ ≤ 6 · 10 −3 ), under beam conditions that deviate only moderately from round, initially Gaussian bunches of equal brightness that collide with a zero or small nominal crossing angle. In the course of investigating the impact of actual departures from that idealized limit, several areas of further study have been identified, most notably: • improved single-bunch models of non-Gaussian tails (Sec. 4.3), both factorizable and nonfactorizable, that are more sophisticated than the simple double-Gaussian functions considered in Eqs. (36) and (37); • the interplay between non-factorization, non-Gaussian tails, and beam-beam corrections, and in particular the application of beam-beam corrections to two-dimensional grid scans; the physics regime (ξ ≤ 10 −2 ), with emittance scans as the primary use case; • the combined impact, during emittance scans at low β * in the physics regime, of higherbeam-beam parameter values and of the hourglass effect on the crossing-angle dependence of beam-beam biases; • the impact of parasitic crossings on the absolute accuracy of emittance scans during routine physics running; • the impact of lattice non-linearities on the absolute accuracy of emittance scans during routine physics running.

Acknowledgments
We thank our colleagues in the LHC Luminosity Monitoring and Calibration Working Group (LLCMWG), and in particular V. Balagura, X. A Impact of non-Gaussian distortions of the source field Because B*B is fundamentally a weak-strong model, beam-beam-induced shape-distortions of the witness bunch, i.e. of the multiparticle distribution, influence neither the shape of the source bunch nor the spatial dependence of the electromagnetic field it generates. The latter is computed once and for all at initialization time, and remains static both turn-by-turn and as a function of beam separation. COMBI, in contrast, is based on a strongstrong model: both bunches are represented by macroparticles; both play, with respect to their partner, the dual roles of source and witness bunch; and the electromagnetic fields they generate are updated on a turn-by-turn basis, and therefore also at each scan step. In order to estimate the impact of separation-dependent, non-Gaussian field distortions caused by the mutual interaction of the opposing bunches, the σ vis bias (Eq. (34)) has been computed using COMBI in two different modes: • for each witness bunch in turn, the field that its macropraticles are subjected to, is computed under the assumption that the transverse macroparticle distribution of the source bunch is a two-dimensional Gaussian, with the position of its centroid and its RMS widths computed on every turn from the actual spatial distribution of the "source" macropraticles. This procedure tracks the evolution of the relative position and of the width of the opposing multiparticle distributions, but the field calculation ignores potential deviations of the bunches from their initial, perfectly Gaussian shape. This is the mode in which all the COMBI results in this paper have been obtained. • for each witness bunch, the field it is subjected to is computed by the HFMM method [34] from the actual charge distribution in the source bunch. In this mode, each macroparticle in the witness bunch is subjected to an electromagnetic field computed as the vector sum of the contributions of all the macroparticles in the source bunch. Since this computation must be repeated on every turn and for the opposing bunch, it is much more CPU-intensive than the Gaussian-field approximation above. Figure 42 displays the beam-beam-parameter dependence of the σ vis bias over a ξ range that extends all the way up to 25 × 10 −3 , more than four times the highest ξ value encountered in vdM scans, and significantly higher than the largest per-IP beam-beam parameter that is sustained operationally in the actual LHC rings. The two calculations agree within 1% or better; for the ξ range and the nominal tunes considered in Table 1, this translates into a σ vis difference of less than 0.005%. Computing the electromagnetic field as that of an equivalent Gaussian, the parameters of which evolve turn by turn, therefore constitutes a fully justified approximation whose impact on beam-beam corrections can be safely neglected.

B Luminosity Integration in COMBI
In COMBI, the calculation of the luminosity, or more precisely of the overlap integral (Eq. (1)), is based on a complete description of the macroparticle distribution in the two colliding bunches: no assumption is made about the shape of the density distributions when evaluating this integral. The latter is calculated on a turn-by-turn basis, both at the scanning IP and at the non-scanning IP(s), if any. The positions of all the macroparticles are updated at each IP on every consecutive turn, and each step in a beam-separation scan is initialized and simulated separately. The transverse distribution of the macroparticles, called H B x,y , is discretized, separately for the two beams B (B = 1, 2), on a two-dimensional grid in the x-y plane. The number of cells in each grid is given by N cells = n × m, where n and m represent the number of bins in the x and y directions respectively. The grid boundaries are located at a distance k i × σ 0 i (i = x, y) from the center of the grid, where σ 0 i is the initial, unperturbed nominal transverse beam size inferred from the input emittance and β * values in the i plane; the scale factor k i is typically set to 12 in both the positive and the negative direction along the x and y axes. The cell area is thus given by ∆S = ∆x × ∆y = 2k x σ 0 x /n × 2k y σ 0 y /m .
The separation between the two beams is taken into account when filling two-dimensional histograms of the macroparticle distributions. With these definitions, the discretized macroparticle density distribution is given by where N part is the total number of tracked macroparticles. The overlap density of the two bunches, i.e. the productρ 1 (x, y)ρ 2 (x, y) in Eq. (1), is therefore represented by: The bunch luminosity can be calculated from the overlap integral I ovlp as: where n 1 and n 2 are the bunch populations and f r the revolution frequency. The overlap integral is estimated by the two-dimensional trapezoidal method: The instantaneous luminosity in Eq. (48) is calculated at each turn, and the result is averaged over the total number of selected turns. Typically a few hundred turns are necessary for the result to stabilize. The reliability of this method was confirmed by benchmarking it against analytical calculations, with the beam-beam effects turned off. Figure 43 shows the evolution of the σ vis bias as a function of the number of macroparticles used in the simulation. The result converges to within 0.01% of its asymptotic value for 5×10 6 macroparticles. Most of the simulation results presented in this paper are based on 10×10 6 macroparticles per bunch, implying that the results are numerically stable at the 0.001% level.
The uncertainty associated with statistical fluctuations in the discretization of the transversedensity distributions was evaluated separately, and typically cancels out when computing luminosity ratios such as [L/L 0 ]. It becomes significant only at large beam separation, when the overlap integral is computed from a small number of macroparticles in the tails. At these scan points, however, the luminosity values are very small, and therefore have a negligible impact on the estimation of the beam-beam bias factors.

C Parameterization of optical-distortion corrections
The polynomial parameters p 0k , ..., p 9k from which to calculate the optical-distortion luminosity-bias factor [L/L 0 ] Opt (ξ R , q x , q y |∆ k ) = p 0k + p 1k ξ R + p 2k q x + p 3k q y + p 4k ξ 2 R + p 5k q 2 x + p 6k q 2 y + p 7k ξ R q x + p 8k ξ R q y + p 9k q x q y (49) defined in Eq. (27) of Sec. 4.2.3, have been tabulated, for normalized nominal-separation values covering the range 0 < |∆ k /σ 0 | < 6 in steps of 0.25. Here the index k (k = 1, 25) labels the nominal-separation bin ∆ k in the scanning plane considered. This parameterization is valid in the fully symmetric Gaussian-beam configuration with zero non-scanning IPs, over the threedimensional space of fractional tunes (q x , q y ) and round-beam-equivalent beam-beam parameter (ξ R ) values defined in Table 9. The values of β * , γ, σ 0 and n in Eq. (13) were chosen such that they correspond to typical settings for pp vdM-calibration sessions during Run 2 of the LHC. The coefficients p 0k , ..., p 9k are listed in Tables 10 and 11 for horizontal and vertical vdM scans, respectively; they are also publicly accessible, in computer-readable form, in Ref. [37]. In the unabridged tables, i.e. for the studies   Opt luminosity-bias factor as a function of the normalized nominal separation, for vertical vdM scans. The functional form is defined in Eq. (49).