Van der Meer scan luminosity measurement and beam–beam correction

The main method for calibrating the luminosity at Large Hadron Collider (LHC) is van der Meer scan where the beams are swept transversely across each other. This beautiful method was invented in 1968. Despite the honourable age, it remains the preferable tool at hadron colliders. It delivers the lowest calibration systematics, which still often dominates the overall luminosity uncertainty at LHC experiments. Various details of the method are discussed in the paper. One of the main factors limiting proton–proton van der Meer scan accuracy is the beam–beam electromagnetic interaction. It modifies the shapes of the colliding bunches and biases the measured luminosity. In the first years of operation, four main LHC experiments did not attempt to correct the bias because of its complexity. In 2012 a correction method was proposed and then subsequently used by all experiments. It was based, however, on a simplified linear approximation of the beam–beam force and, therefore, had limited accuracy. In this paper, a new simulation is presented, which takes into account the exact non-linear force. Depending on the beam parameters, the results of the new and old methods differ by ∼1%\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim 1\%$$\end{document}. This needs to be propagated to all LHC cross-section measurements after 2012. The new simulation is going to be used at LHC in future luminosity calibrations.


Calibration methods
An absolute value of the luminosity or the cross section can be measured at an accelerator by separating the beams in the transverse plane and performing the so-called van der Meer scan [1]. To illustrate the idea, let us consider the collision of two bunches with N 1,2 particles moving in the opposite directions. If the first bunch is separated by − x, − y in the plane perpendicular to the beams, the average number a e-mail: balagura@cern.ch (corresponding author) of interactions μ with the cross section σ normalized by the number of particles is μ( x, y) where the subscript "2" of the coordinates x 2 , y 2 refers to the stationary second beam, L is the integrated luminosity and ρ 1,2 (x 2 , y 2 ) are the normalized transverse particle densities of the unseparated bunches when x = y = 0. For example, if ρ 1 (x, y) is a delta-function δ(x, y), the shifted density δ(x + x, y + y) peaks at − x, − y (note the minus sign). Integration over x and y drastically simplifies (1.1) since ρ 1 (x 2 + x, y 2 + y)ρ 2 (x 2 , y 2 )dx 2 dy 2 d x d y=1, (1.2) as can be easily proved by substituting x 1 = x 2 + x, y 1 = y 2 + y. In the new variables the integrals ρ 1 (x 1 , y 1 )dx 1 dy 1 and ρ 2 (x 2 , y 2 )dx 2 dy 2 decouple and reduce to unity by definition.
From (1.1) and (1.2) we obtain van der Meer formula σ = μ( x, y) The ratio is often called the "specific" number of interactions. Their total number accumulated during the scan, μ sp d x d y, gives σ . Though van der Meer method is well known, the formula (1.3) is sometimes reduced in the literature to the Gaussian bunch densities or the case when μ sp is distributed independently in x and y, discussed in the next section. Here, we present the method in its full generality. This is required, in particular, for explaining the novel two-dimensional scans [2], which will probably be in wide use at LHC in Run 3. After analyzing van der Meer scans for several years, the author tries to share the accumulated experience on various calibration details in the first half of the paper. The discussion is concentrated on the general method and its accuracy but not on the detector effects varying from one luminometer to the other. The second half of the paper is fully devoted to the beam-beam electromagnetic interaction, which is one of the main factors limiting the accuracy.
The derivation of (1.3) uses only μ sp but not μ or N 1,2 separately. Therefore, it remains valid even if μ and N 1 N 2 change arbitrarily but proportionally during the scan, e.g. due to a gradual decrease of beam currents with time. It is required, however, that ρ 1,2 remain constant.
Equation (1.3) can also be understood from other but equivalent perspective. The transverse movements of the first bunch smear and"wash out" its profile ρ 1 , so that effectively it becomes constantρ 1 . This reduces the complicated overlap integral to ρ 1 ρ 2 dx 2 dy 2 =ρ 1 , i.e. to unitȳ ρ 1 = ρ 1 d x y = 1 if ρ 1,2 are normalized. Equivalently, the scan can be viewed from the transverse position of the first beam i.e. in the coordinates x 1 , y 1 . Then the first beam is stationary while the second moves. In this case ρ 2 (x 1 − x, y 1 − y) is "washed out" by x, y movements, and the overlap integral reduces to ρ 1 dx 1 dy 1 = 1 and drops out as before.
To make a parallel translation of one beam in one transverse direction, one needs 4 magnets placed at the corners of the trapezoid-like beam trajectory. Therefore, to steer two beams in two directions one needs 4 × 2 × 2 = 16 magnets per one interaction point. It is not easy to synchronize precisely all of them and ensure a parallel translation of the beams with constant speeds. Therefore, the scan, e.g. at LHC, is performed not "dynamically" using the single continuous pass, but stepwise. I.e. the function μ sp ( x, y) is measured only in the predefined set of discrete points and interpolated between them or fitted to some analytic function to get the final integral. In moving the beams from one point to the next, one waits when the slowest magnet reaches it is desired current value, and only then the luminosity measurement starts.
At e + -e − colliders van der Meer method is biased by strong beam-beam interactions. The luminosity is usually measured "indirectly" by counting Bhabha scattering e + e − → e + e − events, which have a high cross-section precisely known from quantum electrodynamics. Contrary to that, in the collisions of non-elementary hadrons it is difficult to find a physical process with accurately predicted cross-section and convenient for the detection. In the hadron accelerators, the best accuracy is achieved by measuring the luminosity "directly" using its definition (1.1).
At LHC, which will only be discussed in the following, in addition to van der Meer scans there are two alternative direct methods of the luminosity calibration. They utilize precise vertex detectors. The first is the so-called beam-gas imaging [19]. Here, the profiles ρ 1,2 (x, y) of the bunches are "revealed" in their interactions with a tiny amount of gas in the beam pipe. One effectively records the bunch "photos" using the vertex detector as a "camera" and the gas as a "film". Unfolding the images with the vertex resolution yields ρ 1,2 densities. To improve the accuracy, one also uses the high statistics profile of the "luminous region" formed by the interactions of two bunches. This provides a powerful constraint on the product ρ 1 ρ 2 . The overlap integral is then calculated analytically from the reconstructed ρ 1,2 densities. This method is used up to now only at LHCb [8,11], which has a dedicated gas injection system [20], an excellent vertex detector and a flexible trigger suitable for recording beamgas interactions.
The second method is called the beam-beam imaging [21]. It is very similar, but the role of the gas plays another beam effectively "smeared" by the transverse movements in van der Meer scan. After sweeping e.g. the first bunch, it effectively becomes a wide and uniform "film" independently of the initial ρ 1 distribution. It allows making a "photo" of ρ 2 with high statistics. Alternatively, the same vertex distribution data can be viewed from the transverse position of the first bunch, where it is effectively stationary. The accumulated image gives a "photo" of ρ 1 "filmed" by smeared ρ 2 . In the original van der Meer approach the smearing allows integrating over x y and reducing the overlap integral ρ 1 ρ 2 dx dy to unity. In the beam-beam imaging the same integration is applied to the luminous region profile, i.e. to the product ρ 1 ρ 2 not integrated over dx dy, and allows reconstructing the individual densities ρ 1,2 .
Up to now, the beam-beam images were taken at LHCb [8] and CMS experiments [22]. Since in both imaging methods the overlap integrals are calculated from the measured vertex distributions, they have different systematic errors and are complementary to the "classical" van der Meer approach where the vertex distributions are ignored. All three methods might achieve similar levels of accuracy. The imaging methods are more complicated, however, because they require a deconvolution with the vertex resolution comparable to the transverse bunch widths. Due to the simplicity and sufficiency of the classical van der Meer technique, it remains the main tool of the luminosity calibrations at LHC.

X-Y factorization
It is impossible to guide particles exactly parallel to a beam axis. Therefore, in any accelerator the optic elements are designed such that the particles going away from the axis are sent back and in the end "oscillate" in the transverse plane. This creates the transverse bunch widths and determines the density profiles ρ 1,2 . In more detail this will be discussed in Sect. 3.1. To ensure stable operation, the accelerator is designed such that the oscillatory motions are separately stable in the transverse coordinates x and y and are almost independent of each other. Any "coupling" between the coordinates could create extra resonances in the oscillatory motions and, therefore, should be avoided. The bunch densities can often be factorized into x-and y-dependent parts: (1.5) From (1.1) it follows that the specific number of interactions then also factorizes, μ sp ( x, y) = μ x sp ( x) · μ y sp ( y). This is sufficient to simplify the two-dimensional integral μ sp ( x, y)d x d y and to reduce it to a product of one-dimensional integrals along the lines x = x 0 and y = y 0 . Indeed, . (1.6) The integrals in the enumerator can be measured in two onedimensional scans over x at fixed y 0 and vice versa. Note that the formula is valid for any point ( x 0 , y 0 ). This is rarely stressed in the literature. It might be advantageous to choose ( x 0 , y 0 ) not far from the point of maximal luminosity to collect sufficient statistics of interactions. There might be another advantage if the beam coordinates are not accurately measured. The potential slow drifts of the beam orbits from their nominal positions might affect both scanned and not scanned coordinates and bias the luminosity measurement. The bias from not scanned coordinate is minimized at the maximum of μ where the derivative of e.g. μ sp ( x 0 , y) on x 0 is zero. Performing a pair of one-dimensional scans instead of an expensive two-dimensional scan allows saving the beam time. Reducing the time also helps to minimize the influence of the slow drifts of the beam orbits if they are not accurately measured. Therefore, at LHC the cross-sections are usually calibrated using (1.6) instead of (1.3). This approach, however, relies on the x-y factorizability of μ, which is good at LHC but not perfect. It can be violated by many factors, essentially, by any imperfection in the accelerator leading to an x-y coupling.
The remaining non-factorizability is usually studied using the distributions of the interaction vertices. After unfolding with the vertex resolution they give the products of two bunch densities. Ideally, the shape of their projections to one coordinate should remain invariant when scanning another coordinate. The deviations are interpreted as the non-factorizability and are propagated to the cross-section corrections. This procedure is complicated because it requires the characterization of two unknown bunch shapes using only one luminous region profile. One can use the imaging methods to measure directly ρ 1,2 densities [23]. However, the beam-gas interactions have limited statistics while the beam-beam imaging suffers from the uncertainties in the beam positions and the beam-beam systematics discussed later. Because of the complexity of the x-y non-factorization studies, usually the bunch shapes are fitted assuming some smooth bunch shape model. This makes them model-dependent. The cross-section corrections due to x-y non-factorizability and the associated systematic errors are typically at the level 1%.
The accuracy can be improved further by performing twodimensional scans over the central region giving the dominant contribution to the integral in (1.3). This approach was pioneered at LHCb in 2017 [2]. Scanning only the central region was relatively fast but allowed evaluating ∼ 90% of the integral. The method is model-independent and allows reducing the non-factorization systematics by an order of magnitude.

Crossing angle between the beams
In van der Meer scans at LHC the beams are not always opposite. They may collide at a small angle of the order O(100 µrad). This separates the beams outside the interaction region and suppresses possible parasitic collisions between nominally not colliding bunches. It is not immediately obvious how van der Meer formalism should be extended to the case of not parallel beams. In addition, the particles in the two beams can be different. Up to now, LHC has performed van der Meer scans with the proton and lead ion beams. It might be not clear whether in the asymmetric proton -ion collisions the cross-section can be calibrated in the laboratory frame, or it is necessary to make a transformation to the center-of-mass system.
These questions were answered in [21] using two alternative derivations. In the first, the simple two-dimensional integral in (1.1) was extended to four dimensions and taken over x, y in the same way as (1.3) was obtained from (1.1). In the second, equivalent derivation the direct mathematical calculation was substituted by simple physical arguments. They will be elaborated in more detail below.
Let's denote the velocities of the beam particles by v 1,2 as shown in Fig. 1. They can be decomposed into two parallel (v 1,2 ) and one common perpendicular There exist infinitely many relativistic frames where the beams are parallel. In any of them, (1.3) is valid, for example, in the center-of-mass or the rest frame of one of the particles. It is even valid in the frame where the particles move in the same direction one running after the other. All such frames can be obtained by boosting the laboratory frame first with the velocity v ⊥ and then with an arbitrary velocity parallel to v. Indeed, in the laboratory frame boosted by v ⊥ , or, equivalently, after the "active" boost of the beam particles by −v ⊥ , their perpendicular momentum transforms to p ⊥ = γ ⊥ p ⊥ −γ ⊥ β ⊥ E/c, which is zero since by construction β ⊥ = p ⊥ c/E. Here, c is the speed of light, E is the energy in the laboratory system, β ⊥ = v ⊥ /c and γ ⊥ = (1 − β 2 ⊥ ) −1/2 are the beta-and gamma-factors of the boost, respectively. Therefore, after the active boost by −v ⊥ the beams become parallel to v and any further boost along v preserves this parallelism.
Note that although the perpendicular velocity after the boost by −v ⊥ becomes a simple difference v ⊥ − v ⊥ = 0, the remaining velocity v is, of course, not equal to the difference v − v ⊥ = v . This would be the case for the Galilean but not for the Lorentz transformation. The correct formula follows most easily from the opposite transformation from the primed to the laboratory frame: Let's find out how van der Meer formula (1.3), proved in the primed coordinates with the parallel beams, modifies in the laboratory frame. The quantities μ, N 1,2 and σ are Lorentz-invariant, only the transverse area d x d y is not. Following notations of this subsection, the latter will be denoted in the primed coordinates as d x d y . Its transformation to the laboratory system requires some explanations given below. This material complements the discussion in [21].
The space-time coordinates x = (t, x) of any particle moving with the four-momentum p = (E, p) satisfy the equation where λ is a free parameter and x 0 = (t 0 , x 0 ) is an arbitrary point on the particle trajectory corresponding to λ = 0. Four equations in (1.7) with one free parameter define a line in the four-dimensional space. The values of λ uniquely label the line points and can be expressed, for example, via the time coordinate: λ = (t − t 0 )/E. The beam displacements = (0, ) during van der Meer scan change the x 0 parameter: Note that the definition of the scan implies that the beams are moved only spatially, so there is no time component in . In the following, it will be convenient to decompose spatial vectors into three components: parallel to v ⊥ , to the vectorproduct [ v × v ⊥ ] and to v. They will be denoted by the subscripts ⊥, × ⊥ and , respectively, e.g. (1.9) The × ⊥ component is perpendicular to the plane of Fig. 1.
Ideally, the beam displacements should be orthogonal to v but we consider below the general case = 0. After the active boost by − v, the line defined by (1.8) transforms to Note that the solutions x, x of (1.8) and (1.10) for the same λ correspond to the same four-dimensional point in the laboratory and primed frames, respectively. After the boost, the beam displacement acquires the time component t = −β ⊥ γ ⊥ ⊥ . In other words, it is impossible to transform a spatial scan in the laboratory to merely a spatial scan in the primed frame. To resolve this complication one can use the following argument. Any system with the parallel beams is special in the sense that the number of interactions created by two particles is determined only by the transverse distance between their lines. It does not depend on the initial positions of the particles x 0 + along the lines or on the initial time t 0 + t since the particles are anyway assumed to travel from t = −∞ to t = +∞. Contrary to the transverse initial coordinates, the time and longitudinal shifts do not change the particle line. Therefore, the scan with from (1.11) is equivalent to the one with where the time and longitudinal coordinates are simply set to zero. Comparing this equation with (1.9) one can see that the beam displacements ×⊥ transform from the laboratory to the primed system without any changes while the displacements ⊥ are "extended" by the γ ⊥ -factor. Therefore, the area element d x d y in the primed coordinates is larger than in the laboratory system by γ ⊥ : Note that one can invert the arguments and make the opposite boost from the primed scan with the displacements (0, ⊥ , ×⊥ , 0) to the laboratory system: Here, the γ ⊥ -factor again appears in the transverse components but now in the laboratory system. Concluding from this formula that γ ⊥ d x d y = d x d y with γ ⊥ on the opposite side would be a mistake, however, since the time and longitudinal components can be freely changed and zeroed only in the primed frame with the parallel beams. Setting t = β ⊥ γ ⊥ ⊥ to zero in the formula above is not allowed and spoils the equivalence of the scans in the laboratory and the primed frames.
Using (1.13), the cross-section formula in the laboratory system can finally be written as (1.14) Note that the relativistic correction γ ⊥ depends on the velocities but not on the momenta or masses of the particles. For example, it coincides for proton and lead ion beams if their velocities are the same. The formula is relativistically invariant and is valid in any frame. The area element d x d y by definition lies in the plane perpendicular to v. Let's denote this plane by P. If the scan planeP is inclined with respect to P at an angle α, an area element d˜ x d˜ y onP should be projected to P, i.e.
The longitudinal translations along v do not matter. The beam crossing angles at LHC are small, so β ⊥ < 10 −3 and γ ⊥ − 1 < 10 −6 . Therefore, the relativistic correction at LHC can be safely neglected.
The luminosity, however, is modified as it can be seen from the following. The typical longitudinal sizes of the LHC bunches σ i L along the beam i = 1, 2 in the laboratory frame are 5-10 cm. They are much larger than the transverse sizes σ i T ∼ 100 µm. Since γ ⊥ ≈ 1, the Lorentz and Galilean transformations from the laboratory to primed coordinates are almost equivalent. The latter preserves the bunch shapes. Therefore, in the primed system the longitudinal spread σ i L gets projected to v ⊥ at the angles α 1,2 shown in Fig. 1. For the Gaussian bunches this increases the primed transverse width σ i⊥ in the beam crossing plane to (1. 16) This formula is valid up to the second-order α i -terms enhanced by σ i L /σ i T 1 factor. Note that often in the literature σ i⊥ is expressed in the form of a rotation as (σ i T cos α i ) 2 + (σ i L sin α i ) 2 . Writing cos α i = 1 − α 2 i /2 + · · · instead of unity as in (1.16) implies its validity at least up to the terms ∝ α 2 i not enhanced by σ i L /σ i T . At this level of accuracy one can not neglect the difference between Lorentz and Galilean transformations, however, and should take into account the relativistic corrections.
The exact formula can be obtained using the same formalism as for van der Meer scan above. Let's interpret the parameter = (0, ) in (1.8) not as the beam displacement but as a stochastic variable describing the spatial spread of the particles in the bunch in the laboratory frame. For example, can be a zero-mean random variable normally and independently distributed along the transverse x, y and longitudinal axes with the standard deviations σ x,y,L , respectively. Let's consider the general case when neither x nor y lies in the crossing plane of Fig. 1. Then, the transverse bunch widths in the laboratory frame along ⊥ and × ⊥ directions are Here, the notation like cos(x i , ⊥) denote the cosine of the angle between the x direction of the i-th bunch and v ⊥ . According to (1.12), the bunch spread σ ⊥ is multiplied by the γ ⊥ -factor in the primed frame, while σ ×⊥ remains unchanged: If one of the transverse axes, e.g. x, lies in the crossing plane, this simplifies to For ultrarelativistic beams with v 1,2 ≈ c, like at LHC, and arbitrary α 1,2 , which should be approximately equal in this case, α 1 ≈ α 2 ≈ α, the beta-and gamma-factors can be expressed via the angle α: (1.20) This gives The smallness of α i at LHC is partially compensated by the large σ i L /σ i T ratio, so due to the crossing angle the effective transverse widths σ i⊥ increase by 5-20%. The luminosity reduces by the same amount. Once again, this does not modify van der Meer formula (1.14), since the normalization ρ i (x, y) dx dy = 1 remains invariant. For broader bunches, one just needs to enlarge proportionally the region of integration.
Note that if the transverse bunch densities ρ 1,2 factorize in the primed directions x and y but none of them lies in the crossing plane of Fig. 1, σ i L has non-zero projections on both x and y . This makes x and y distributions slightly correlated and to some extent breaks the x -y factorizability.

Luminosity calibration accuracy
Van der Meer scan is the main tool of the absolute luminosity calibration at LHC. For the given beam particles and the LHC energy, the scan is performed in every experiment at least once a year to check the stability of the luminosity detectors. One or two LHC fills with carefully optimized experimental conditions are allocated for this purpose. The accuracy is determined by various factors discussed below. The overall calibration uncertainty is typically 1-2%.
The calibration constant is then propagated to the luminosity of the whole physics sample using linear luminometers. In ATLAS and CMS operating at higher pile-up μ-values, the accurate linearity is required in larger dynamic range since van der Meer scans are typically performed with μ 1. The uncertainties caused by the deviations from the linearity due to irradiation ageing, long-term instabilities etc. are reduced by comparing several luminometers with different systematics. Therefore, the overall luminosity uncertainty is often dominated by the calibration error.
Averaging over many colliding bunch pairs and several van der Meer scans often reduces the statistical error of the calibration to a negligible level. The systematics from the bunch population N 1 N 2 measurements is typically at the level ∼ 0.2-0.3% except in the very first LHC scans in 2010 with low-intensity beams.
The cross-section has the dimensionality of a length unit square. It appears in van der Meer formula (1.14) due to the integration over x and y. Any error in the beam displacements directly affects the cross-section. An accurate measurement of x, y scale is performed in a dedicated "length scale calibration" (LSC), which always accompanies van der Meer scans. The simplest LSC is described below.
Let's assume that the true beam positions i in the laboratory frame for the beam i = 1, 2 can be written as where a i , b i are unknown vectors close but not exactly equal to the unit x, y-vectors and α i , β i are the nominally set val-ues for x, y beam movements, respectively. The latter are known exactly. Potential nonlinearities in the i dependence on α i , β i e.g. due to beam orbit drifts, are neglected here but will be briefly discussed later. The constant vectors 0 i corresponding to α i = β i = 0 are also unknown.
In the simplest LSC the beams are nominally displaced by the same amount, i.e. with α 1 = α 2 = α L SC , β 1 = β 2 = β L SC . If the bunches have equal shapes, the center of the luminous region O L SC is positioned in the middle between them, (1.23) It can be accurately measured by the vertex detectors. Most often van der Meer scans at LHC are performed such that the beams are displaced symmetrically in opposite directions corresponding to α 1 = −α 2 = 2α sym , β 1 = −β 2 = 2β sym . The advantage of the symmetric scan is that it allows to reach maximal separations in the limited allowed range of the beam movements. The distance between the beams 12 = 1 − 2 is then 12 = α sym a 1 + a 2 2 A comparison of (1.23) and (1.24) shows that the measurements of O L SC in the vertex detector are sufficient to calibrate the scales (a 1 +a 2 )/2, (b 1 +b 2 )/2 necessary for the symmetric scans. The fact that the shapes of the colliding bunches are different can normally be neglected at LHC after averaging over many colliding bunch pairs. Non-symmetric scans depend on other linear combinations than (a 1 + a 2 )/2 and (b 1 + b 2 )/2 measurable in (1.23). The simplest way to calibrate the lengths |a 1,2 | = a 1,2 and |b 1,2 | = b 1,2 individually is to use the measurements of the luminosity. Ideally it should be stationary during LSC. Any small variation indicates that the distance between the beams changes.
For example, let's consider the LSC in the x-direction. Assuming x-y factorizability and neglecting the angle between a 1,2 and the x-axis, one arrives at the scalar equations The movement of the luminous region between any two LSC points, The corresponding change of the x-distance between the beams 1 12,x − 2 12,x can be deduced from the luminosity change L 1 L SC − L 2 L SC and van der Meer scan data. For example, if one of the x-scans was performed symmetrically, the beam separation change required to modify the luminosity by a given amount can be calculated from the derivative Together with (a 1 + a 2 )/2 from (1.26) this allows to calibrate the scales a 1,2 individually.
To improve the sensitivity of the method and to increase L 1 L SC − L 2 L SC , LSC can be performed at a point close to the maximum of the derivative d L/d 12 , where the second derivative is zero. For the Gaussian bunches with the widths σ 1,2 , the luminosity dependence on the beam separation is also Gaussian with the sigma = σ 2 1 + σ 2 1 , and the optimal LSC beam separation is 12 = .
After the calibration, the length scale systematics is typically well below 1%.
Note that LSC is not needed for the "static" imaging methods, namely, for the beam-gas and also for the beam-beam imaging if in the latter the reconstructed bunch is stationary during the scan in the laboratory frame. In both cases, the reconstruction is performed in the vertex detector coordinates, which always define an accurate scale.
The remaining most important sources of systematics include x-y non-factorizability of the bunch densities, the beam orbit drifts and the so-called beam-beam effects due to the electromagnetic interaction between two colliding bunches.
As it was discussed in Sect. 1.2, the x-y non-factorizability can be circumvented by performing two-dimensional scans over the central region giving the main contribution to the integral in (1.14). Three two-dimensional scans, already performed at LHCb at the end of Run 2, allowed to reduce this uncertainty approximately by an order of magnitude [2].
The beams can drift at LHC by a few microns leading to the cross-section uncertainties at the level of 1%. Further improvements require accurate monitoring of the beam positions. The LHC Beam Position Monitors (BPMs) could not provide the necessary accuracy in Run 1 because of the temperature drifts in the readout electronics. In Run 2 they were upgraded and all interaction points were equipped with the so-called DOROS BPMs. The accuracy was significantly improved and reached a sub-micron level. This was proved by calibrating and comparing with the beam positions reconstructed with the beam-gas imaging at LHCb. The latter is relatively slow and requires one or a few minutes to reach the required accuracy even with the gas injection. This is sufficient for the calibration, however, and after that, the DOROS BPMs can accurately measure even fast beam drifts with 0.1 second time resolution. Sufficient accuracy can, possibly, be achieved also at other experiments without the beam-gas imaging using the correlations between the DOROS measurements and the positions of the luminous centers. Using well-calibrated DOROS BPM data, one can significantly reduce the scan-to-scan non-reproducibility and potentially achieve the overall calibration accuracy below 1%.
The last beam-beam systematics is caused by the electromagnetic interaction between the bunches. The electromagnetic force kicks the beam particles and modifies their accelerator trajectories and the bunch densities ρ 1,2 . If the perturbations were constant during the scan it would not produce any bias since van der Meer formula (1.3) is valid for any densities. The kick strength, however, depends on the transverse profile of the opposite bunch and the distance to it. The perturbation of the densities ρ 1,2 , therefore, depends on x, y. For example, it vanishes at large beam separations. Such ρ 1,2 dependence on x, y breaks the derivation of (1.3) from (1.1) and introduces biases in the cross-section formulas that are difficult to estimate. The beam-beam correction is the main subject of this paper and will be discussed in detail in the following sections.
The lead ion bunches in LHC van der Meer scans carry much smaller charge than the proton ones. Therefore, the beam-beam systematics is more significant for the protonproton scans.
The beam-beam interaction also affects the imaging methods. Like the classical van der Meer approach, the beambeam image is biased when ρ 1,2 densities vary during the scan. The beam-gas imaging is also biased when the beambeam interaction changes, e.g. during van der Meer scan at the same or any other LHC experiment. In the latter case, the associated beam-beam distortions propagate through the accelerator everywhere in the ring. If the beams are stationary, however, the densities ρ 1,2 are constant and the appropriate ρ 1,2 fit model can describe the bunches accurately. The required model might be complicated, however, and dependent on the constant beam-beam interactions at all LHC experiments.
In the first LHC publications, the experiments either not considered the beam-beam systematics or assigned ∼ 1% error to their luminosity measurements in the proton-proton scans without making any correction. In 2012 the method [24] was proposed for correcting the classical van der Meer technique. It was subsequently used by all LHC experi-ments, and the systematic error was reduced to 0.3-0.7% [9][10][11]13,[15][16][17][18]. However, the beam-beam force in this method was oversimplified and approximated by a linear function of the transverse coordinates. More specifically, the electromagnetic field was described by the dipole and quadrupole magnets responsible for the offset and the slope of this linear function, respectively. This was discovered by the author of this paper in January 2019 using a new, independently developed simulation. It will be presented in this paper. Instead of the linear approximation, the simulation uses the accurate formula of the beam-beam force. Unfortunately, the old and new beambeam corrections differ by ∼ 1% as shown in Fig. 11. The difference is dependent on van der Meer scan beam parameters. This requires the corresponding rescaling of all LHC cross-sections after 2012 that were based on the luminosity calibrated with the old oversimplified beam-beam model.
The new simulation is primarily oriented at the classical van der Meer method. It is optimized for calculating the luminosity but not the bunch shapes required in the imaging methods. Some limited tools for predicting the shapes are implemented in the simulation, however, and can be extended.
The beam-beam force depends on both transverse coordinates and, therefore, introduces x-y coupling and nonfactorization. The new simulation allows to correct the luminosity measurements at each point of van der Meer scan such that the beam-beam perturbation is effectively removed together with its x-y coupling. The cross-section can then be calculated from the corrected μ values using unmodified (1.3) or (1.6).
The new simulation is sufficiently general. The bunch profiles can be approximated by an arbitrary weighted sum of the Gaussians with the common center. The x-and y-widths can be different. In addition, the luminosity correction can be calculated in the presence of the beam-beam kicks at an arbitrary number of interaction points. The bunch shapes of all colliding bunches are specified individually.

Momentum kick induced by the beam-beam interaction
For LHC physics one usually considers the collision of two protons (or ions) ignoring other particles in the bunches. Contrary to that, the beam-beam electromagnetic interactions have a long-range and act simultaneously between many particles. At large distances, one may neglect quantum effects and use classical electrodynamics. Any associated electromagnetic radiation of protons or ions at LHC will be ignored. The formula of the momentum kick induced by the beambeam interaction is well known in the accelerator community. However, it might be not so easy to find in the literature its rigorous derivation with the discussion of all simplifying Fig. 2 The electrical field from the charge q 1 at rest acting on the charge q 2 from another bunch assumptions affecting the accuracy. To make the material of this paper self-contained, we present below the derivation of this formula from the first principles.

Electromagnetic interaction of two particles
As it will be shown in a moment, at LHC the beam-beam force changes the transverse particle momentum in the laboratory frame by a few MeV, i.e. negligibly compared to the total momentum. Therefore, one can assume that all particles creating the electromagnetic field move without deflections with constant velocity. Since the equality and constancy of velocities are preserved by boosts, this approximation can be used in any other frame.
The electromagnetic field of a particle with the charge q 1 is simplest in its rest frame where it reduces to the electrical Coulomb component E 1 shown in Fig. 2. The momentum kick exerted on a particle q 2 from another bunch can be calculated in this frame as an integral of the infinitesimal momentum changes along the trajectory. According to our assumption, the velocity of q 1 is constant, so its position is fixed. Since q 2 in the rest frame of q 1 has even larger momentum than in the laboratory, while the transverse kick is the same, one can safely assume that its speed dz/dt = β 0 c is also constant and q 2 moves along the straight line denoted as the z-axis in Fig. 2.
If the velocities of the colliding particles in the laboratory system are ±βc, β 0 can be expressed as which is the double-angle (or double-rapidity) formula for the hyperbolic tangent, the analog of tan(2φ) = 2 tan φ/(1− tan 2 φ). Of course, at LHC one can safely assume β ≈ β 0 ≈ 1. The momentum kick received by q 2 is given by the line integral Because of the reflection invariance z → −z, the zprojection of the kick should vanish, so only the perpendicular component of the electric field is written in (2.2).
Its integral E ⊥ 1 (x, y, z) dz depends only on the transverse coordinates. It will be denoted in the following simply by E 1 (x, y). According to Gauss's flux theorem applied to the cylinder with the radius R shown in Fig. 2, it is equal to where 0 is the electric constant. Here, we are using the system of units where the first Maxwell's equation is written as where α = e 2 /(4π c 0 ) is the fine-structure constant, e is the proton charge, Z 1,2 = q 1,2 /e are either 1 or 82 for proton or lead ion LHC beams, respectively, and is the reduced Plank constant.
Since the momentum kick p 2 and R are perpendicular to the z-axis, they are conserved by the boosts along this line. Therefore, the corresponding terms in (2.4) remain invariant and equal in all frames where q 1,2 velocities are parallel. The perpendicular electric field E 1 (x, y), however, depends on the boosts and acquires a magnetic counterpart. In (2.4) it should be used only in the q 1 rest frame.
If q 1 and q 2 collide in the laboratory system with the small crossing angles α 1,2 as in Fig. 1, (2.4) receives corrections of the order α 2 1,2 . For example, they can be calculated by boosting the kick (2.4) from the primed to the laboratory system with the velocity v ⊥ and the subsequent projection to the planes transverse to the beams. At LHC these corrections are negligible.
Equation (2.4) is sufficient for a rough estimation of the momentum kick induced by the whole bunch. If its charge ∼ 10 11 e is condensed into q 1 , the kick at, for example, one bunch sigma ∼ 100 µm is equal to p = 3 MeV/c. This is, indeed, negligible compared to LHC energies. The exact formulas of the Gaussian bunch fields are presented later in (2.10) and (2.12) and lead to the same conclusion. Equation (2.4) was derived for the stationary charge q 1 . Due to the momentum conservation, however, it receives the opposite kick p 1 = − p 2 and starts moving. Our formulas do not depend on the mass of q 1 and, therefore, should be valid even when the mass is much less than p 1 ∼ 3 MeV/c. But then after the kick q 1 becomes ultrarelativistic, and its field can not be described by the simple electrostatics.
To solve this seeming contradiction one should recall that the field from the ultrarelativistic charge q 2 in Fig. 2 is concentrated only in a thin "pancake" perpendicular to z and travelling together with the particle. So, q 1 receives the kick p 1 when q 2 passes z = 0. Only after that and almost instantaneously q 1 becomes ultrarelativistic. Let's denote this moment by t (z = 0). It then takes some time to propagate this information back to q 2 , which is escaping almost at the speed of light. Namely, q 2 "sees" the initial stationary field from q 1 until it passes the forward light cone emitted from q 1 at t (z = 0). This moment t (z) can be found from From this consideration we see again that the assumption of constant q 1,2 velocities is valid only in the ultrarelativistic limit. Therefore, the dependence 1/β 0 in (2.4) is not justified and should be dropped. The formula should be written as under the explicit condition β 0 ≈ 1. The corresponding angular kick for the particles with the momentum p is Note that in the present literature the angular kick and the related parameters traditionally and most often are expressed via the classical particle radius r c = q 2 /(4π 0 mc 2 ) determined by the particle mass m. For example, one can refer to the so-called beam-beam parameter. In the recent Particle Data Group review [25] it is defined in Eq. (31.13) as for the y-direction, where β y2 is the beta-function discussed later, m 2 , γ 2 are the mass and γ -factor of q 2 and m e , r e are the electron mass and classical radius. Writing mass in the formulas is misleading, since, as it was explained above, the kick does not depend on m. In the ultrarelativistic case β 0 ≈ 1 the momentum change p depends only on the electric charges q 1,2 and R, since it is determined by Gauss's or Coulomb's law. If the mass m is introduced in the beam-beam equations, it should necessarily cancel as in the combination m e r e above. For example, the kick is the same for protons and electrons if their momenta are the same, despite the difference in their classical radii. To stress this invariance, the fine-structure constant α should be used instead of the classical radius because only the electromagnetic interaction is relevant here. It is better to drop completely the mass m from the formulas.

Simplifying assumptions in the particle interaction with the opposite bunch
Let's demonstrate that for calculating the kick from the whole bunch one can assume that all bunch particles move in the same direction with the same speed. As it will be discussed in Sect. 3.1, the angular spread in the laboratory frame is of the order δα = σ T /β O(10 −5 ) where σ T is the transverse bunch size (40-100 µm) and β is the beta-function in the range 1-20 m during van der Meer scans in the interaction points at LHC. Note also that the kick is perpendicular to the velocity difference v = v 1 − v 2 from Fig. 1, so the kick angular variation is of the same negligible order δα.
There is one effect where this spread is enhanced. An angular deviation of one particle changes its crossing angle with respect to the opposite bunch. This affects the transverse bunch width σ T visible from the particle according to (1.16). Therefore, the angular variation δα leads to the effective smearing of σ T : In spite of the large enhancement factor σ L /σ T 1000, the values of σ L /β O(10 −2 ) and α O(10 −4 ) are so small that in van der Meer scans the variations of the transverse width δσ T /σ T 0.01 can be neglected in the beam-beam kick calculations.
The longitudinal momentum spread δp/ p of the beam is completely negligible for our purposes, since the beam-beam kick is determined by the velocities that are close to the speed of light at LHC and almost insensitive to the momentum change. Namely, if v is the velocity corresponding to the rapidity φ, v = tanh φ, its change is δv = δφ/ cosh 2 φ = δ(sinh φ)/ cosh 3 φ = (δp/ p)·β/γ 2 ∝ 1/γ 2 . The associated angular variation of v = v 1 − v 2 due to the crossing angle is additionally suppressed by the smallness of α < 10 −3 .
Finally, the angular variation due to the beam-beam kick itself is also small, p/ p ∼ 10 −6 . Since the typical longitudinal bunch length σ L is 5-10 cm, the kick has no time to develop to a sizable displacement during the interaction. The particles should travel freely much longer distances of the order σ T · p/ p ∼ 100 m before their transverse displacements reach σ T . However, the accelerator elements controlling the transverse movements correct the trajectories and bring the particles back. In Sect. 3.8 Eq. (3.37) it will be shown that in the end the beam orbit is shifted by less than 1% of the bunch width. The angular distribution shifts by p/2 p 10 −6 . The beam-beam luminosity bias typically does not exceed 1%. Therefore, to achieve the required overall calibration accuracy of 0.1% and to estimate the bias with the relative uncertainty < 0.1%/1% = 10%, it is sufficient to calculate the momentum kicks using the electromagnetic fields of the unperturbed densities ρ 1,2 .
If one can assume that the particles in the bunches move with constant and opposite velocities, this greatly simplifies our four-dimensional electromagnetic problem and reduces it to the two-dimensional electrostatics. Indeed, in (2.3) one can easily recognize the Coulomb's law in two dimensions. The circle circumference 2π R in the denominator substitutes the sphere area 4π R 2 in the three-dimensional Coulomb's law in accordance with the Gauss's electric flux theorem. Therefore, from (2.5) is just the Coulomb's two-dimensional force between q 1 and q 2 . The calculation of p 1,2 kick, i.e. the problem of the electromagnetic interaction of the ultrarelativistic laboratory bunches, reduces to the calculation of the two-dimensional electrostatic forces between the transversely projected static charges in the frame with the parallel beams.
As it was already discussed, the longitudinal bunch distributions do not matter in this frame. Indeed, the accumulated kick remains invariant if particles in the opposite bunch are arbitrarily displaced longitudinally as long as they follow the same lines and traverse the whole interaction region.

Electrostatic field from two-dimensional Gaussian distribution
In this subsection we present the formulas of the electrostatic field from the two-dimensional Gaussian density For the round bunch with σ x = σ y = σ , the azimuthally symmetric field can be determined from the charge Q(1 − e −R 2 /2σ 2 ) inside the disk x 2 + y 2 < R 2 and the Gauss's flux theorem: Therefore, the beam-beam angular kick of the particle with the charge Z 1 e induced by the round Gaussian bunch with N 2 particles with the charges Z 2 e is The field from an elliptical bunch with σ x = σ y is more complicated. It was derived by Bassetti and Erskine in [26]: where the path-independent integral is taken in the complex plane between the points It can be expressed via the complex error function erf(z) = 2 z 0 e −ζ 2 dζ / √ π or its scaled version named Faddeeva function as . (2.15) Note that in [26] the sign in front of y 2 /2σ 2 y was misprinted as plus. A simplified proof of Bassetti-Erskine formula found by the author will be published in a separate paper.
The Faddeeva function w(z) grows exponentially when the imaginary part I m(z) of its argument tends to −∞. In this case, calculating the difference between two large numbers in the enumerator of (2.15) becomes numerically unstable. In practice, to ensure the positiveness of I m(z 1,2 ), the calculation can be performed in the following way. In the case σ x < σ y the formulas might be applied with the swapped x-and y-directions. The obtained components E x and E y should be swapped back. This ensures that the square roots in (2.13) are always taken with σ x > σ y and, therefore, are real. Then I m(z 1,2 ) becomes negative only if y < 0. Since the field is centrally symmetric, E(x, y) = −E(−x, −y), this case can be circumvented by calculating the field at the opposite point (−x, −y) and by inverting the signs of the obtained components E x , E y .

Average kick of bunch particles
Up to now, we have discussed the kicks of individual particles. In this subsection we present a simple formula for the kicks averaged over the bunches. For the Gaussian shapes, it was derived in Appendix A in [27]. Here we give an alternative proof based on simple arguments and extend the formulas to arbitrary ρ 1,2 .
Let's denote the momentum kick of the i-th particle in the first bunch exerted by the j-th particle in the second by p i j . It can be calculated as p i j = F i j /c where is the two-dimensional Coulomb's force between the charges q i , q j separated by r i j = r i − r j in the transverse plane. The full force on the first bunch is the sum i, j F i j that can be approximated by the integral 17) where N 1,2 are the number of particles in the bunches. Since F depends only on the difference r = r 1 −r 2 , it is convenient to use r as the integration variable where ρ also depends only on r: This is the cross-correlation ρ 2 ρ 1 or, equivalently, the convolution ρ 1 * ρ 2 whereρ 2 (r) = ρ 2 (−r) i.e. ρ 2 (r) with an opposite argument. Equations (2.17) and (2.18) show that r spread in the integral can be equivalently represented either by the two bunch densities ρ 1,2 or by only one ρ. In the latter case the second bunch effectively collapses to the point-like charge at the origin. Indeed, (2.18) follows from (2.17) if ρ 1 and ρ 2 are substituted by the artificial bunch density ρ = ρ 1 * ρ 2 and by the delta-function at zero, respectively. This is illustrated schematically in Fig. 3. The left picture shows the overall electrostatic force j F i j exerted by the second bunch on the i-th particle. The origins of F i j vectors are varied according to the density ρ 2 (x, y). To obtain the full force, one needs to sum over i, i.e. to vary the ends of F i j vectors. Figure 3b shows such variations for the j-th particle of the second bunch. Since parallel translations do not change vectors F i j , the variations of their ends can be obtained by the opposite variations of the origins, as shown in the picture. The full sum i, j F i j in Fig. 3c, therefore, can be calculated by smearing the origins with both probability densities ρ 1 (−r) and ρ 2 (r). Equivalently, it can be calculated by varying the ends of r by ρ = ρ 1 (r) * ρ 2 (−r) in agreement with (2.19). For the Gaussian bunches the convolution ρ 1 * ρ 2 is again Gaussian with the sigmas x,y = σ 2 1x,y + σ 2 2x,y . Since the momentum is conserved, the full momentum kicks of two bunches are opposite, F 1 = −F 2 . The average kicks of the particles are equal to F 1 /N 1 , F 2 /N 2 and are different if N 1 = N 2 .
(a) (b) (c) Fig. 3 The force between the bunches does not change when one bunch density ρ i (r) (i = 1, 2) is collapsed to the point charge at the origin while the other is convolved with ρ i (−r). For the Gaussian ρ 1,2 densities the convolution has sigma = σ 2 1 + σ 2 2 effect. Below we describe a new numerical simulation developed for this purpose named "B*B" or "BxB" (pronounced "B-star-B") [28,29]. Before going into details, let's briefly remind the transverse dynamics of the particles in an idealized accelerator, the so-called "betatron motion".

Recurrence relation
Every particle in the beam oscillates around a stable orbit with a constant amplitude. Ideally, the oscillations in x and y are independent. They are described by the Hill's equation where u is the transverse coordinate (x or y), s is the circular coordinate along the ring, u = ∂ 2 u/∂s 2 and K u (s) is a function defined by the quadrupole accelerator elements whose field is proportional to u. The solutions of (3.1) are where u is a constant defining the oscillation amplitude and called "emittance" in the accelerator language, β u (s) is the so-called "beta-function" determined by the equation .
is the "phase advance" whose value at s = 0 is denoted by φ 0,u . From (3.2) and (3.4) one can calculate u = ∂u/∂s, i.e. the tangent of the angle between the particle and the orbit: At the interaction point the beams are maximally focused to reach the maximal luminosity. As it follows from (3.2), β u is then minimal and β u = 0, so that (3.5) simplifies to After every turn in the accelerator the particle phase advance increases by the constant called the "tune". Here, the integral is taken over the whole accelerator length L. As it follows from (3.2) and (3.6), at the interaction point it is convenient to merge the phase-space coordinates (u, u ) to one complex variable Its evolution is described by the simple rotation in the complex plane where z n+1,u and z n,u = √ u β u e i(φ n,u −φ 0,u ) are the complex coordinates at the turns n+1 and n, respectively. Note that the minus sign in (3.8) is chosen according to the minus in (3.6), so that the rotation is counter-clockwise by definition. Since there are two transverse axes x and y, there are independent rotations in the two complex planes z x = x − i x β x , z y = y − iy β y and the full phase-space is four-dimensional.
The bunch transverse shapes are typically approximated by Gaussians. The distributions in every complex plane then have a form of the two-dimensional Gaussians with the same standard deviations along u and −u β u . They are invariant under rotations around the origin and transform to themselves after every accelerator turn.
As it was discussed, the beam-beam kick changes the angle u while the instantaneous change of u is negligible. Equation (3.9) then modifies to the recurrence relation z n+1,u = (z n,u − iβ u u )e 2πi Q u . (3.10) According to (2.8), the angular kicks in the first bunch, for example, are determined by ( x , y ) = q 1 E 2 (x, y)/ pc. The electrostatic field is given by (2.10) or (2.12) for round and elliptical bunches, respectively. The beam-beam deformations of the bunch creating the field are neglected, as it was explained in Sect. 2.2. The strategy of the B*B simulation is, therefore, the following. In the beginning, the particles are distributed in the phase-space according to the given initial density. Then, they are propagated through the accelerator turn-by-turn using (3.10) and the change of the luminosity integral (ρ 1 + δρ 1 )ρ 2 dx dy (3.11) is calculated. To take into account the beam-beam perturbation of the second bunch shape, the simulation is repeated with the swapped bunches yielding ρ 1 (ρ 2 + δρ 2 ) dy. The full luminosity change with respect to the unperturbed value is approximated as The second-order term δρ 1 δρ 2 dx dy is neglected. If the bunches are identical, two terms in (3.12) coincide and it is sufficient to perform one simulation and to double the correction.
The main challenge of the numerical calculation of the beam-beam modified luminosity (3.11) is the required accuracy. It should be negligible compared to other systematic uncertainties, i.e. less than 0.1% at LHC. Reaching this level with the Monte Carlo integration in the four-dimensional phase space requires simulating many particles and too much CPU time. Therefore, several optimizations are implemented in B*B that are described below.

Particle weights
In (3.12) only the integrals of perturbed and unperturbed densities are required. Contrary to the unperturbed profile defined in B*B by a continuous analytic formula, the other density, e.g. ρ 1 + δρ 1 in (3.11), should be represented by the point-like particles. This can be achieved by splitting the full phase-space into volumes V i , i = 1, 2, . . . . Each of them can be assigned to one "macro-particle" with the weight w i equal to the phase-space density integrated over In this way any continuous density can be approximated by a weighted sum of delta-functions placed at the macro-particle coordinates. The integral from (3.11) can then be expressed as the discrete sum (3.13) Here, the density ρ 2 is considered constant in x-y projection of each V i and substituted by its value ρ 2 (x i , y i ) at the particle position (x i , y i ).
Let's assume that N 1 real particles in the first bunch are approximated by N B * B 1 N 1 macro-particles. Let's define that the association of the real and macro-particles does not change, so that each V i always contains the same N 1 V i (ρ 1 + δρ 1 )d 4 V = N 1 w i particles. With this definition, the volumes V i deform due to the beam-beam force but the weights w i are conserved as it follows from the conservation of particles. To simplify notations, the simulated macro-particles in the following discussion will also be called "particles" as the meaning will be clear from the context.
Since w i are conserved, in B*B they are defined using the simplest unperturbed density ρ 1 as explained below. To ensure a good sampling of the four-dimensional space, B*B adopts a two-step approach. Initially, the macro-particles are distributed at the circles whose radii r x,y form an equidistant grid where n i x,y = 1, 2, . . . , n max . The index i runs across all 1, 2, . . . , n 2 max radii pairs (r x , r y ), where is the configurable parameter of the simulation. Each pair receives one macro-particle placed randomly at the circles in the z x , z y planes. In this way the sampling of the absolute values |z x |, |z y | is realized. The sampling of z x,y phases is performed by the accelerator simulation itself. After every turn the particle gets rotated by 2π Q x,y angles in the corresponding planes around the origin and slightly shifted vertically by the beam-beam kick according to (3.10). Since the beam-beam interaction is small at LHC, the particle trajectories remain approximately circular. After N B B turn ∼ O(100−1000) accelerator turns the particle well samples its circular trajectory and z x , z y phases. It was proved that the initial choice of random phases has negligible impact on the final integral. Equation In principle, it is possible to assign initially not one but several particles to the i-th pair of z x,y -circles with the radii (r i x , r i y ). This would reduce the phase sampling dependency on the tunes and the associated evolution in the accelerator. The choice of one particle per circle in B*B was made to simulate more accelerator turns using the same number of calculations. This allows checking that no new effects appear after a very large number of turns.
For the normally distributed density in the complex plane where u = x, y, and the equidistant radii r u from (3.14), the macro-particle weights w i are given by the integrals over the rings n i u u < r u < (n i u + 1) u : The final weight of the particle placed at the radii (r i x , r i y ) is the product In van der Meer scans at LHC, the Gaussian bunch approximation is not always sufficient. The x, y-projections sometimes can be better described by the sum of two Gaussians. For flexibility, B*B simulation allows defining ρ 1,2 as the sum of an arbitrary number of Gaussians with configurable weights and widths and independently in x and y. The field E(x, y) is calculated from each Gaussian individually using (2.10) or (2.12). After weighing, all contributions are summed. To speed up the simulation, the field map is precalculated in the beginning, and then the interpolations are used. This is especially important in the case σ x = σ y when the field should be computed using the complicated Bassetti-Erskine formula (2.12). Similarly, both the initial weights w i and the densities ρ 2 (x i , y i ) at every turn appearing in (3.13) are calculated as weighted sums of the contributions from all Gaussians.
The radii limits n max x , n max y in (3.14) in the multi-Gaussian case are chosen in the B*B simulation such that they have the same efficiency as r x < N σ σ x and r y < N σ σ y cuts for the simple single Gaussian shape. The parameter N σ is configurable. Its default value N σ = 5 is usually sufficient to reach the required accuracy. For the single Gaussian bunch the excluded weight, i.e. the density integral of the not simulated region (r x > 5σ x or r y > 5σ y ), is 1 − (1 − e −5 2 /2 ) 2 = 7.5 × 10 −6 . To additionally reduce the CPU time, the largest (r x , r y ) pairs are not simulated, namely, it is required that r x n max x 2 + r y n max y 2 < 1. (3.19) This increases the lost weight to (5 2 /2 + 1)e −5 2 /2 = 5.0 × 10 −5 . The remaining weights are normalized. The default value of the configurable N part parameter from (3.15) is 5 000. The real number of generated particles is 23% less due to (3.19). This default value is used in the simulation examples discussed in the following. The other beam parameters are listed in Table 1. They are taken from [24] to compare its old, biased method used at LHC in 2012-2019 and the results of the B*B simulation. Increasing N part to 100 000 with the default B*B settings changes relatively the luminosity integral by ≤ 4 × 10 −5 in the full simulated range of bunch separations from 0 to 200 µm.

Stages in the simulation
The bias associated to the phase space limit (3.19) and to the approximation of the continuous integral by the discrete sum in (3.13) partially cancels in the ratio R = (ρ 1 + δρ 1 )ρ 2 dx dy ρ 1 ρ 2 dx dy (3.20) if both integrals are taken numerically in the same way. Therefore, B*B starts from simulating N no B B turn accelerator turns without the beam-beam interaction i.e. with the beam-beam kick u = 0 in (3.10). The unperturbed luminosity ρ 1 ρ 2 dx dy is estimated turn-by-turn using (3.13). The default value of N no B B turn is conservatively chosen to be 1000. In Sect. 3.4 it will be explained, however, that for the tunes Q x,y with two digits after the comma, like at LHC, any multiple of 100 leads to identical results. Therefore, N no B B turn = 100 is sufficient in this case. With the beam parameters listed in Table 1 this value gives the relative deviation between the numerical and analytical ρ 1 ρ 2 dx dy integrals less than 2 × 10 −4 in the practically important region of the beam separations contributing ∼ 99.9% to the crosssection integral in (1.3). The resulting bias of the ratio R in (3.20) should, therefore, be smaller.
After the first N no B B turn turns the beam-beam kick is switched on. The user has two options: either instantaneously apply the nominal kick u in (3.10) or increase it slowly or "adiabatically", namely, linearly from zero to the nominal value during N adiab turn turns. In the former case, the particle trajectory instantaneously changes from the ideal circle to the one perturbed by the beam-beam force. The intersection of the two trajectories is the last point on the ideal circle. It is positioned randomly, and different points on the ideal circle lead to different perturbed trajectories. This is depicted in Fig. 4. Two initially opposite points, marked in the figure by the small open circles, create two outer blue trajectories. Their evolution is followed during 10 7 turns and every 1000th point is shown in the plot. The region in grey in the middle is filled with all other trajectories. The unperturbed ideal circle with the center at the origin is shown by the green dashed line. As one can see, the center of the blue circles is shifted due to the change of the orbit in x and x . This will be discussed in more detail in Sect. 3.8. Note that the ideal green circle is infinitely "thin", but the blue points are scattered because of the beam-beam x-y coupling and the variations in the other z y -projection. The "thickness" of the blue trajectories increases with the force strength. For much larger forces the trajectory becomes significantly non-circular.
In the adiabatic case, the trajectories change slowly and the particles have time to redistribute over them. The initial position on the circle then has little importance, and the whole initial circle transforms to approximately one final trajectory shown in red in Fig. 4. Here, the beam-beam interaction is slowly switched on during 1000 turns and then, as in the previous case, every 1000-th turn is shown out of 10 7 in total. The red trajectory has approximately the same spread as the  Fig. 4 The trajectories of the particles with r x = 20, r y = 30 µm in z x projection for the bunches separated in x by 40 µm with the parameters from Table 1 as an example. The outer blue trajectories are formed by two initially opposite particles (marked by small open circles) after the instantaneous switch of the beam-beam force. The grey band between them is composed of such trajectories from all particles. The adiabatic trajectory is in red and the initial circle with r x = 20 µm is shown by the green dashed line blue outer ones but less than the grey band composed of many trajectories initiated by the instantaneous switch. In any case, the spread is negligible compared to the width of the opposite bunch, namely, 40 µm in Fig. 4. Therefore, the contribution to the overlap integral in (3.13) is essentially determined by the infinitely "thin" average trajectory, which is the same in both instantaneous and adiabatic switch cases. The integral does not depend on the way how the beam-beam force is switched on. This is demonstrated in Fig. 5 by the blue and red points for various beam separations.
In both cases only 500 (100) accelerator turns were simulated to determine the perturbed (initial) overlap integral. In the adiabatic case the force was switched on during 100 turns. The error bars show the standard deviations of the results obtained with different random generator seeds. Clearly, the adiabatic switch is preferable. The instantaneous switch leads to the spread of the trajectories and larger statistical fluctuations of the final integral.
The central black points in Fig. 5 are the result of the simulation when the beam-beam force was switched on slowly in 10 4 turns and the luminosity integral was calculated over 10 6 turns. As one can see from Fig. 5, simulating only 500 turns already gives sufficient accuracy. The default B*B val- After the beam-beam force is fully switched on and before calculating the overlap integral over N B B turn turns in the final stage, the B*B simulation has an option to run N stab turn turns for a "stabilization". This period is not used for calculating the luminosity integral. Normally this parameter can be set to zero because, for example, in the adiabatic case the particle arrives at its final trajectory as soon as the force reaches its nominal value. The following evolution does not change the trajectory. N stab turn parameter exists only for flexibility and for experimenting with the B*B simulation.

Lissajous curves
Without the beam-beam interaction the x-y trajectory of a particle placed at the z x , z y -circles with the radii (r x , r y ) and with the initial phases φ 0 x,y is described by (3.9): x n = r x cos(2π Q x n + φ 0 x ), y n = r y cos(2π Q y n + φ 0 y ), (3.21) i.e. appears to be a Lissajous curve. The beam-beam force leads to the dispersion of the trajectory as shown in Fig. 6  Fig. 6 Left: the trajectory of one particle with r x = 20, r y = 30 µm followed during 5000 turns in the presence of the beam-beam force. The bunch parameters are from Table 1  In general, for the rational Q x,y and a small beambeam force the trajectory "cycles" after LC D(Q x , Q y ) turns where LC D denotes the lowest common denominator. In the above example LC D(3/10, 1/3) = 30, so after one "cycle" with 30 turns the phase increments 30 · Q x,y become integer and the following turns pass through approximately the same region of the phase-space. Therefore, they do not immediately improve its sampling. After many more turns the beambeam spread forces the points to fill and to sample the whole rectangle, but then the simulation takes too much CPU time.
To improve the sampling but keep the number of turns low, one can artificially shift Q x,y by a small amount and, for example, make them irrational. This "opens" the Lissajous curve, so that it fills the whole rectangle |x| < r x , |y| < r y even without the beam-beam force. The modification should be sufficiently small to ensure that the overlap integral changes negligibly. In B*B two small irrational numbers δ Q x,y randomly distributed in the interval [− , ] are added to Q x and Q y by default, where is chosen to be e −8.5 /2 = 1.01734 · · · × 10 −4 . As an example, the points obtained with the shifts δ Q x = /2, δ Q y = − /2 in 5000 turns are shown in Fig. 6 left by the red points. They much better fill the rectangle. This modification of the tunes is configurable and can be switched off if desired.
When the trajectory fills the whole rectangle, the probability density of the particle peaks at the boundaries. This is shown in Fig. 6 right for N B B turn = 10 6 . This can be understood from the simple case when the beam-beam force is absent and at least one of Q x,y is irrational, so that the Lissajous curve is open. Then, the density is factorizable in x and y, ρ r x ,r y (x, y) = ρ r x (x)ρ r y (y). The terms ρ r u (u), where u = x, y, can be derived by projecting the uniform density dφ/2π of the z u -circle with the radius r u to the u-axis: Therefore, . (3.23) It is interesting that the smooth two-dimensional xy Gaussian shape of any bunch is always intrinsically composed from such peaking ρ r x ,r y densities. As follows from (3.17), they should be taken with the weights ∝ exp(−r 2 x /2σ 2 x − r 2 y /2σ 2 y )dr 2 x dr 2 y . Conversely, only the distributions decomposable to ρ r x ,r y can represent the bunch shapes ρ 1,2 (x, y).
The functions K u (s) in (3.1) and the accelerator phase advances (3.4) can be modified by configuring the quadrupole currents. The LHC tune values in all Run 1 and 2 van der Meer scans were deliberately kept constant: Q x = 64.31, Q y = 59.32. The fractional parts of these values have two digits after the comma, and their lowest common denominator is LC D(0.31, 0.32) = 100. This is larger than 30 from the example above, so the phase-space sampling is better. The tiny tune change δ Q ∼ 10 −4 mentioned above introduces a sizeable spread to the trajectory only after O(10 3 ) turns. Figure 5 shows, however, that 500 turns are already sufficient to get the required accuracy for the bunch settings from Table 1 if the beam-beam force is switched on adiabatically. With 500 turns the small irrationality δ Q x,y ∼ 10 −4 does not matter. The luminosity bias due to modifications of Q x,y is also negligible, it is less than 5 × 10 −5 for all bunch separations. Therefore, δ Q x,y plays no role here but might be useful in the cases when LC D(Q x , Q y ) is relatively small or to improve the accuracy of the overlap integrals in experimenting with individual particles.
With a low number of turns ∼ O(100), it is important to choose them as multiples of LC D(Q x , Q y ). For the LHC tunes, they should be multiples of 100, like N no B B turn = 100, N B B turn = 500 in Fig. 5. This ensures that every particle trajectory is sampled an integer number of "cycles". Not complete "cycles" introduce a bias and a dependence on the initial phases. Figure 7 shows the luminosity R-correction from (3.20) in the x-scan calculated "cycle-by-cycle", i.e. averaged over non-intersecting 100-turn windows. Six facets correspond to different bunch separations. One can see small "oscillations" for large separations. This is the reason for larger red error bars in Fig. 5 at higher bins. They show the level of statistical fluctuations of R-averages over N B B turn = 500 turns, i.e. over = 1000 first turns without the beam-beam force and next N adiab turn = 1000 turns when it is being adiabatically switched on. The facet labels denote the bunch separation in the bunch sigmas. The simulated beam parameters are listed in Table 1 the first 5 beam-beam "cycles" to the right from the second vertical dashed line in Fig. 7.

Beam-beam interactions at several interaction points
During van der Meer scan, at each LHC experiment there are colliding and not colliding bunches. Since ATLAS and CMS detectors are exactly opposite in the LHC ring, due to this symmetry they share the same colliding bunch pairs. On the other hand, these pairs never collide in ALICE and LHCb. The latter can have either "private" bunch pairs not colliding anywhere else or pairs with one or both bunches colliding in one, two or three other experiments. The beam-beam disturbance at any interaction point propagates everywhere in the LHC ring and biases van der Meer calibrations in other experiments.
B*B allows to determine the luminosity correction in the general case of one bunch colliding with an arbitrary number N of fixed bunches of specified geometries at other interaction points. One should also specify their beta-functions β i x,y and the constant phase advances normalized by 2π between the first and i-th points. By definition, the first point has Q 1 x,y = 0. The emittance x,y of the particle in the linear accelerator is conserved everywhere in the ring, therefore, the radii of the circles at different interaction points scale according to (3.2) as r i u = u β i u ∝ β i u . The recurrence relations then take the form Let's consider, for example, the last equation. The kick u N is determined by the electrostatic field at the last interaction point N at the position (x N , y N ) = (Re(z N x ), Re(z N y )). The factor e 2πi Q u −Q N u , where Q u is the full tune, rotates the phase, while the term β 1 u /β N u changes the radii scale from the last to the first interaction point. Note that to simulate the perturbation of the second bunch, the order of the interaction points should be reversed, N → (N − 1) → · · · → 1 → N , due to its opposite direction. Also note that, contrary to the overall tunes, the phase advances at LHC are different for two beams. For example, their values in Run 2 proton-proton van der Meer scans at √ s = 13 TeV are listed in Table 2. Figure 8 shows how the luminosity in a scan at one LHC interaction point is affected by the beam-beam force at other points. As an example, here the bunch of the first beam colliding at all four LHC experiments is simulated. The luminosity change due to the beam-beam perturbation of only this bunch is presented in the plot. To get the full change, a similar contribution due to the perturbation of the opposite bunch should be simulated and added. Note that each of the four opposite bunches might collide with other bunches at other interaction points and this needs to be included in its simulation.
The β-function values for Fig. 8 simulation are listed in Table 2. Together with the ATLAS beam parameters from Table 1 they are taken from [24]. The number of particles and emittances of all bunches are the same. The emittances are determined from the bunch width and β-value at ATLAS. Four plots in Fig. 8 show the R-corrections from (3.20) when the scan is performed along the x-axis at ATLAS, ALICE, CMS and LHCb, respectively. Let's denote the scanned experiment by the index i = 1, 2, 3, 4. The corrections are calculated and color-coded separately when the force is switched on only at the j-th LHC interaction point (R i, j ) or simultaneously at all (R i,1234 ). One can see that the latter, shown by the solid circles, is relatively close to the product of the four corrections R i,1 R i,2 R i,3 R i,4 each calculated in the absence of the beam-beam force at three other points.
In addition to the luminosity corrections at the scanned interaction point, B*B automatically calculates similar corrections at all other experiments where the beams remain stationary. In the simulation of Fig. 8, they collide head-on. The corresponding R-corrections are not shown because they are not needed for van der Meer calibration. They less depend on the separation and are closer to unity, the mismatches do not Table 2 The phase advances Q i x,y defined in (3.24) for two LHC beams in x and y directions with respect to ATLAS in Run 2 proton-proton van der Meer scans at √ s = 13 TeV. Last column lists the beta-functions β x = β y from [24] used in Fig. 8 4 . The beams at not scanned interaction points collide head-on. The simulation parameters are listed in Tables  1 and 2 exceed 4×10 −3 . B*B also allows to simulate the beam-beam interactions at multiple points and to calculate all luminosity corrections when the scans are performed simultaneously in several experiments.
3.6 Invariance under β x,y , p, and N scalings As one can see from Table 2, the beta-functions β i x,y and, therefore, the bunch widths σ i ∝ β i in the simulated example are different at ATLAS/CMS, ALICE and LHCb. ALICE has the largest bunches and the smallest beam-beam kick.
The solid black curves in Fig. 8 correspond to R i,i correction when the beam-beam force is switched on only at the scanned experiment. One can see that all such curves in the four facets are the same and also identical to Fig. 5. Let's understand why this happens despite different bunch widths and kicks.
Scaling of both β x and β y by a constant factor β x,y → αβ x,y changes the linear scale in the z x,y planes by √ α. The two-dimensional electrostatic field from a point charge q drops with the distance as ∝ 1/R, so the angular kick changes as u → u / √ α. According to the definition (3.8), in the z u complex planes the angular kick is additionally multiplied by β, so z u = −iβ u u scales as z u → √ α z u . Therefore, the distributions of two bunches and the kicks z x,y scale proportionally to √ α. New particle trajectories can be obtained by the simple √ α-scaling of the whole z x,y planes. This modifies the luminosities by α but keeps constant their ratios, e.g. R in (3.20).
To stress this invariance, it is better to rewrite (3.25) via the ratios z u / √ β u : where "i + 1 index denotes the next interaction point, e.g. the first after the last, while Q "N +1 u is the full tune Q u . Variables in the complex planes z u / √ β u have another advantage: they do not change across the interaction points in the ideal linear accelerator. Indeed, both the phase Arg(z u ) = φ u and the absolute value z u / √ β u = √ u are invariant, since the emittance u is conserved. So, the complex planes z x / √ β x and z y / β y are common to all interaction points. As we have seen, the beam-beam kick −iβ i u u i / β i u in (3.26) happen to be invariant to the scaling β i x,y → α i β i x,y with arbitrary α i because of ∝ 1/R Coulomb's law in two dimensions. Therefore, the trajectories in the z u / √ β u planes drawn from the fixed initial points z 1 u / √ β u using ( Contrary to the tune values, the phase advances are specific for each experiment. Therefore, with multiple interaction points this symmetry between the experiments is lost and R i,1234 and R i, j , i = j curves in Fig. 8 are different. There is another interesting invariance of the beam-beam perturbations. During acceleration of the particles, their angles u decrease due to the growth of the longitudinal momentum. If the acceleration is sufficiently slow, like at LHC, the emittance decreases according to the "adiabatic damping" formula u =˜ u · m/ p, (3.27) where m is the mass of the particle and˜ u is the constant normalized emittance. Therefore, when p increases, the complex plane z u / √ β u shrinks as 1/ √ p and the beam-beam force increases as √ p. The angular beam-beam kick p/ p contains p in the denominator as in (2.6), and decreases as 1/ √ p, i.e. coherently with the z u / √ β u complex plane. Therefore, in the complex plane corresponding to the normalized emittances z u √ p/mβ u = √˜ u e iφ u the trajectories remain invariant if the phase advances are constant and˜ u are conserved. The beam-beam perturbations of the bunches with identically distributed normalized emittances are identical.
Finally, let's formulate the scaling law for the number of particles N 2 in the opposite bunch. It will be called "second" in the following to distinguish from the "first" bunch with the studied trajectories. The bunch variables will be marked by the corresponding indices.
Multiplying the normalized emittance of the first bunch by α increases the linear scale by √ α and reduces the beambeam angular kick by 1/ √ α. In order to have the same √ α change in the linear and angular scales for the first bunch, the beam-beam force from the second should be enhanced by α. This can be achieved by increasing Z 1 Z 2 N 2 . Therefore, the simultaneous scaling of˜ 1u and Z 1 Z 2 N 2 by the same factor changes the scale but not the shapes of the first bunch trajectories. It also keeps constant the luminosity correction ratios. If the trajectories are analyzed in the complex plane the results depend only on the initial distributions, the phase advances and β 1x /β 1y ratios, but do not depend on the individual values of β 1x,1y , 1u , p 1 , Z 1,2 or N 2 . Only their combinatioñ (3.29) matters. This can also be confirmed with the B*B simulation.

Simulation of the beam crossing angle
If the beams collide with the crossing angle as in Fig. 1, the betatron oscillations occur in the planes transverse to the beam vectors v 1,2 . They determine the transverse widths σ i T of the bunches i = 1, 2. The luminosity, however, depends on the bunch widths σ i T perpendicular to v. According to (1.21), they get additional contributions from the longitudinal widths σ i L . The beam-beam kick is also perpendicular to v either in the laboratory or, for example, in the rest frame of the first bunch particle q 1 shown in Fig. 2. Therefore, q 1 "sees" the opposite bunch width σ 2T enhanced by σ 2L , and this value should be used in the B*B simulation for ρ 2 . The bunch creating the field is "static" in the B*B model, its betatron transverse motion and σ 2T do not matter.
However, to simulate the betatron trajectories in the first bunch, the initial widths σ 1x , σ 1y should be specified in B*B without the longitudinal component. The betatron oscillations are insensitive to σ 1L spread. Since the crossing angles are small at LHC, α 1,2 = α < 10 −3 , the kicks calculated in the primed frame with the parallel beams can be propagated without changes to the frame of the betatron motion. Similarly, the simulated x, y beam transverse coordinates can be propagated back for the luminosity calculation.
In this way the luminosity correction is determined for the unperturbed widths σ 1T , σ 2T instead of the required σ 1T , σ 2T . Therefore, the final density ρ 1 +δρ 1 , corresponding to σ 1T , should be additionally smeared by the contribution from σ 1L . This is performed in the B*B simulation in the following way. Let's consider the general case when neither x nor y of the first bunch lies in the crossing plane, and denote by β 1 the angle between the x-axis and the projection of the crossing plane to the x-y plane. The perturbed density ρ 1 + δρ 1 should then be convolved with the twodimensional Gaussian G(x, y) with the sigmas ασ 1L cos β 1 and ασ 1L sin β 1 in x and y, respectively. If, for example, the x-axis belongs to the crossing plane and β 1 = 0, the smearing occurs only along x with the sigma ασ 1L , while G(x, y) projection to y reduces to the delta-function. Instead of smearing ρ 1 + δρ 1 in the overlap integral, it is simpler to perform an equivalent smearing of ρ 2 : since G(r − r ) = G(r − r), where r and r denote the two-dimensional vectors in the (x, y) plane.
To implement this scheme, the B*B simulation takes the values as configurable parameters defining the Gaussian G i (x, y) for each interaction point i. Its convolution with the multi-Gaussian ρ 2 is performed analytically. The result is used in (3.13) instead of ρ 2 for the overlap integral calculation. Note that it receives contributions from both σ 1L and σ 2L . The field is determined from ρ 2 not smeared by G i (x, y), where only σ 2L contributes. The beam-beam perturbation of the longitudinal beam dynamics is neglected and the longitudinal spread is always approximated by a single Gaussian. By default the parameters (3.31) are zero.

Cross-checks and comparison with the old model used at LHC in 2012-2019
In Sect. 2.4 it was shown that the full beam-beam force exerted on the first bunch can be easily calculated as the force between the first bunch collapsed to the origin (0, 0) and the second one with the density "inflated" from ρ 2 toρ 1 * ρ 2 .
In case of the Gaussian bunches with the widths σ 1,2;x,y and the centers r 0 1,2 , the cross-correlationρ 1 * ρ 2 is also Gaussian and has the widths x,y = σ 2 1;x,y + σ 2 2;x,y and the center r 0 21 = r 0 2 − r 0 1 . Its electrostatic field E 2u , given by (2.10) or (2.12), allows to calculate the average angular beam-beam kick u = q 1 E 2u / p 1 c = eZ 1 E 2u / p 1 c, of the first bunch particles. This value is exact if all bunch particles move with the opposite and constant velocities close to the speed of light, so that one can calculate the kicks using the laws of electrostatics.
Let's substitute the beam-beam angular kick by its average. The x-and y-directions then decouple. To simplify notations, let's drop the coordinate subscript u, denote the oneturn phase advance by 2π Q = φ, introduce the constant = β 1 u = β 1 q 1 E 2u / p 1 c and consider for simplicity the case with only one interaction point. The recurrence relation (3.10) then defines the geometric sequence Its real part determines the shift of the beam orbit, while the scaled imaginary part −I m(z 0 )/β 1 = − u /2 gives the angular shift just before the beam-beam interaction. After receiving the kick + u , it flips from − u /2 to + u /2, while the next accelerator turn changes it back to − u /2. The orbit shift expressed in the bunch widths should be invariant and should depend only on˜ 1 Here, E 2u has ∝ Z 2 N 2 /R dependence, so e(E 2u σ 1 )/Z 2 N 2 is invariant. For example, for the round bunches with σ 1 = σ 2 = σ , x,y = √ 2σ separated in x by x, substituting E 2u from (2.10) or (2.11) gives The maximal value of the first term depending only on x/σ is 0.31908 . . . . Therefore, for protons one has (3.37) which is less than 1% for the typical LHC values˜ 1 ≈ 3 − −4 µm, N 2 ≈ (7−9) × 10 10 and φ x /2 = 0.31π , φ y /2 = 0.32π . In case of several interaction points with the phase advances 2π Q i = φ i with respect to the first point where by definition φ 1 = 0, and the constant kicks i , i = 1, 2, . . . , N , the first equation in (3.32) modifies to (3.38) which becomes equivalent to (3.32) after the substitution (3.39) So, we have again a circular trajectory with the center shifted to (3.40) whose real and opposite imaginary parts define the spatial and β 1 -scaled angular shifts, respectively. This formula is very well known in the accelerator physics, which uses different tools for its derivation [30].
For experimenting and debugging, B*B has configurable options to substitute the exact kick formula by its average and to output the x, y orbit shifts. The latter are calculated as weighted sums N i=1 w i u i over the simulated particle coordinates u i , averaged over accelerator turns similarly to (3.13). For example, with (N no bb turn , N adiab turn , N bb turn ) = (100, 100, 500), default other settings and the single interaction point with the parameters from Table 1, the mismatch between the simulated and predicted orbit shifts averaged over the beam separation range [0−200] µm has the standard deviation σ 2 nm. The maximal predicted orbit shift is 0.28 µm. The simulated shifts without the beam-beam force are compatible with zero with σ 0.3 nm.
If the beam-beam interaction is small, even with the exact, not constant kick the particle trajectory remains approximately circular as illustrated in Fig. 4. One particle trajectory corresponds to two circles in z x,y -planes. The orbit shift is then approximately determined by the beam-beam force F tr averaged over the trajectory like shown in Fig. 6, with the x-y density projection (3.23). The circle from Fig. 4 in the z x -plane shifts according to the average x-component of the force F tr x , dependent, however, on the radius in z yplane. So, different z x,y -circle pairs have different F tr and shift differently. However, the ratio r Re/I m of the real and imaginary parts in Eqs. (3.34) and (3.40) remains constant regardless of F tr . Therefore, the shifted centers z 0 = x 0 +iy 0 lie on the line r Re/I m = x 0 /y 0 and the spatial and angular shifts are proportional. Averaging over the bunch should preserve this proportionality. Therefore, if , i denote exactly known average bunch angular kicks in (3.34) and (3.40), these equations should correctly predict the average z 0 shifts. For example, the real part Re(z 0 ) predicts the spatial average or "orbit" shifts. Note that contrary to the exact average angular kicks u and , the values of z 0 , i.e. the angular and the orbit shifts calculated from (3.34) and (3.40), depend on the assumption that the trajectories are approximately circular, which is violated for strong kicks or a large number of interaction points.
As an example, Fig. 9 shows the x-orbit shifts in the same simulation as in Fig. 8. The beam-beam kick is switched on at all four interaction points. The mismatches between the shifts calculated by the B*B simulation and the analytic approximation (3.40) reach 26 nm. Although small in the absolute scale, they are significantly larger than the statistical fluctuations. The orbit shift along y-coordinate, where the beams are not separated, is compatible with zero within σ = 0.4 nm.
In the old beam-beam model used at LHC in 2012-2019, the field E 2 in the momentum kick p 1 = q 1 E 2 /c was Fig. 9 Upper row: x-coordinate averaged over the first bunch at four experiments, lower row: its deviation from (3.40) approximation, both in nm, versus the beam separation expressed in the bunch widths. The column labels denote the scanned interaction points (marked by the solid lines in the plots), at other points the beams collide head-on (the dashed lines). The simulation is the same as in Fig. 8 not taken from (2.10) or (2.12) but was approximated by a simple linear function of the x and y coordinates. The constant term was chosen to reproduce the orbit shift from (3.34). The slopes were taken from the derivatives ∂ p 1 /∂ x, ∂ p 1 /∂ y at the first bunch center. The model was limited to the case when the Gaussian bunches collided head-on in the not scanned coordinate. Under this assumption, the cross-derivatives ∂ p 1,y /∂ x, ∂ p 1,x /∂ y vanish, and the x, y-slopes were taken as ∂ p 1,x /∂ x, ∂ p 1,y /∂ y, respectively. They can be calculated from (2.10) for the round bunches. For example, in the x-scan at the center of the first bunch, where R = x, one has The full kick at x-separation was approximated as the following linear function of x and y: Note that to reproduce the orbit shift in the first formula, the constant term is not equal to the value of the kick at the bunch center k 1 (1 − e − x 2 /2σ 2 2 )/ x. It contains σ 2 1 + σ 2 2 instead of σ 2 in the exponent. Therefore, (3.42) is not a linear expansion of p 1,x .
The linear kick significantly simplifies the analysis. As known from the accelerator physics, in this case the Gaussian bunch remains Gaussian. The offset term is equivalent to the dipole magnet. It shifts the Gaussian center according to (3.34). The linear terms u · ∂ p 1,u /∂u x= x,y=0 represent the quadrupole magnets usually used to focus or defocus the beams [30]. They modify the Gaussian widths according to the formula Another simplification of (3.42) is that the x and y kicks are independent, so the beam-beam x-y coupling is neglected. Multiple interaction points also decouple and the analysis can be limited to the point where the scan is performed. Indeed, other interactions only modify bunch Gaussian widths and centers, but in any case, they are considered as free parameters in van der Meer analyses. Since the transverse positions of the beams at not scanned points are kept constant, the beam-beam modifications are also constant and do not bias van der Meer calibration. Although the beam-beam kick at multiple points have been discussed and simulated in [24], where the old model was introduced, in practice this was never used.
For the Gaussian bunches, the luminosity can be calculated from the known centers and widths. In the old model the luminosity modification due to the orbit shift, induced by the equivalent dipole magnet, was calculated analytically. The contribution from the quadrupole magnet was obtained using LHC MAD-X accelerator simulation for the beam parameters from Table 1. They were tabulated and then extrapolated analytically to other bunch settings. Two contributions were summed.
The results of the old quadrupole simulation for Table 1 settings are shown in Fig. 10 by the open circles. The curves denoted by "D", "Q" and "D+Q" correspond to the dipole, quadrupole contributions and their sum, respectively. As in other plots, the luminosity bias R − 1 from (3.20) due to the perturbation of only one bunch is shown. The analytic corrections due to the dipole orbit shifts (3.34), the quadrupole bunch width changes (3.43) or both are shown by the corresponding solid lines. One can see that the quadrupole MAD-X simulation is well described analytically and (3.43) could be used instead of the tabulated values in practice.
To compare these known results with B*B, its exact kick formula was replaced by (3.42), by only its dipole or quadrupole parts. The simulated values R − 1 are shown in Fig. 10 by the solid points. They are also in good agreement with the analytic predictions "D+Q", "D" and "Q", respectively.
The B*B results with the exact kick are shown by the solid circles connected by the dashed line "E". They were already shown in Fig. 5. These results significantly differ from the old "D+Q" simulation. Both "E" and "D+Q" contain the dipole contribution. To compare with "Q" alone, the B*B simulation was performed with the dipole constant subtracted from the exact kick in the recurrence relation. The result is shown as "E-D" curve. It differs significantly from "Q" and to a larger extent compensates the luminosity reduction by "D". As expected, the sum of "D" and "E-D" is in agreement with "E".
The knowledge of the R-correction allows determining the bias induced on the reference cross-section σ e.g. when using (1.6) for the Gaussian bunches in one-dimensional x, y van der Meer scans. The integrals in (1.6) should be taken over the unperturbed values μ sp ∝ e − x 2 /2(σ 2 1 +σ 2 2 ) 2πσ 1 σ 2 (3.44) modulated by R 2 from Fig. 10 or the corresponding figure for the y-scan. R should be squared to take into account the perturbation of the second bunch. The luminosity corrections for the x-and y-scans are slightly different because of the difference in the fractional parts of the tunes Q x , Q y , 0.31 = 0.32. The resulting biases of van der Meer crosssection σ /σ − 1 are shown in Fig. 11 separately for B*B, dipole (D), quadrupole (Q) approximations and their sum (D+Q). For the proton bunches colliding at one point, the correction depends only on the specific normalized emittancẽ /N assuming identical bunches with N = N 1 = N 2 and β 1 = β 2 . Therefore, this variable is chosen as the horizontal axis. As shown by the dashed line, the difference between  Table 1 B*B and the old model for Table 1 settings is 0.96%. The old linear kick approximation was too simple to describe accurately the beam-beam luminosity bias.

Conclusions
The main tool for the absolute luminosity calibration at LHC is van der Meer scan. Its systematics dominates the overall luminosity uncertainty, which, in turn, gives one of the main contributions to the uncertainty of the accurate cross-section measurements, for example, in the electroweak sector. Various details of van der Meer method are discussed in the paper. Currently, the main sources of systematic uncertainties are the beam orbit drifts, x-y non-factorizability and the bunch deformations induced by the beam-beam electromagnetic interaction. The first two can be significantly reduced by the accurate monitoring of the beam positions and the two-dimensional scans, respectively. The formalism of one-and two-dimensional scans is presented in the paper in detail. Other sources of systematics include the measurements of the bunch populations, the length scale calibration, luminometer detector effects, other bunch shape changes during the scan, e.g. after particle losses, and various unknown factors contributing to the scan-to-scan non-reproducibility. The alternative beam-gas and beam-beam imaging calibration methods are also briefly mentioned.
The beam-beam bias is the main subject of the paper. The derivation of the beam-beam kick is presented in Sect. 2 from the first principles together with the discussion of various approximations and the induced errors. It is shown that under the assumption of the constant and opposite velocities of the bunch particles, the calculation of the beam-beam electromagnetic force reduces to the simple electrostatics between the charges in the transverse plane. In particular, this allows deriving the average kick formulas (2.18) and (2.19) for the bunches of arbitrary shapes.
In the last section, we present the B*B simulation for calculating the beam-beam luminosity corrections. Contrary to the previous model with the linear kick used at LHC in 2012-2019, it is based on the exact nonlinear electrostatic force between the point and the Gaussian charge density, either round (2.10) or elliptical (2.12). The perturbed particle trajectories are followed in the accelerator assuming their ideal transverse betatron motion with the known phase advances. The perturbed luminosity is calculated at the interaction points with maximally focused beams, where the derivative of the β-function is zero, and the elliptical phase-space trajectories of the betatron motion become circular.
The luminosity corrections due to the perturbations of two bunches are calculated separately and then summed. It is assumed that the bunch creating the field is not disturbed by the beam-beam interaction. Therefore, the coherent oscillation modes of the two beams are neglected.
The B*B simulation allows to correct the beam-beam luminosity bias in van der Meer scan point-by-point and to remove it together with the associated x-y non-factorizability. The bunch shapes may be approximated by an arbitrary sum of Gaussians. The electrostatic field is pre-calculated and then the interpolations are used to save CPU time. An arbitrary number of interaction points is allowed. The simulation of different particles is parallelized in the processors with multiple cores. The B*B code is written in C++. It can be used as a standalone application or as a library available in four computer languages: C, C++, python and R.
In Sect. 3.6 it is shown that, for example, the first bunch phase-space trajectories drawn in the complex planes ˜ 1x,1y /Z 1 Z 2 N 2 · e iφ 1x,1y are determined only by the initial distribution of˜ 1x,1y /Z 1 Z 2 N 2 ratios, by the phase advances and β 1x /β 1y ratios, if the normalized emittances˜ 1x,1y = 1x,1y /( p 1 /m 1 ) are conserved. They do not depend on the individual values of the beta-functions β 1x,1y , the emittances 1x,1y , the momentum p 1 , the proton numbers Z 1,2 or the number of particles in the opposite bunch N 2 . This leads to the corresponding invariance of the luminosity correction ratios.
For the bunch parameters from Table 1, i.e. for the reference values of the old model with the linear kick, the B*B simulation predicts 0.96% less van der Meer cross-section correction than the old model. As one can see from Fig. 11, the latter significantly overcorrects the bias. This needs to be propagated to all LHC cross-sections after 2012 taking into account the bunch parameters in the corresponding van der Meer calibrations. For other Run 2 proton-proton van der Meer scans the discrepancies are in the range 0.8-1.4%. The B*B predictions are going to be used at LHC in the future calibration analyses, for example, they already appear in [31].