e1e-mail<EMAIL_ADDRESS>11institutetext:  Laboratoire Leprince-Ringuet,
CNRS/IN2P3, École polytechnique, Institut Polytechnique de Paris

# Van der Meer Scan Luminosity Measurement and Beam–Beam Correction

Vladislav Balagurae1,addr1

(Received: date / Accepted: date)

###### Abstract

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 $\sim 1\%$. 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.

††journal: Eur. Phys. J. C

##  1 Van der Meer scan

###  1.1 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 vanderMeer:296752 . 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 $-\Delta x$, $-\Delta
y$ in the plane perpendicular to the beams, the average number of interactions
$\mu$ with the cross section $\sigma$ normalized by the number of particles is

$\displaystyle\frac{\mu(\Delta x,\,\Delta
y)}{N_{1}N_{2}}=\sigma\frac{L}{N_{1}N_{2}}$
$\displaystyle=\sigma\iint\rho_{1}(x_{2}+\Delta x,\,y_{2}+\Delta
y)\rho_{2}(x_{2},y_{2})dx_{2}\,dy_{2},$ (1.1)

where the subscript “2” of the coordinates $x_{2}$, $y_{2}$ refers to the
stationary second beam, $L$ is the integrated luminosity and
$\rho_{1,2}(x_{2},y_{2})$ are the normalized transverse particle densities of
the unseparated bunches when $\Delta x=\Delta y=0$. For example, if
$\rho_{1}(x,y)$ is a delta-function $\delta(x,y)$, the shifted density
$\delta(x+\Delta x,y+\Delta y)$ peaks at $-\Delta x,-\Delta y$ (note the minus
sign). Integration over $\Delta x$ and $\Delta y$ drastically simplifies
section 1.1 since

$\int\rho_{1}(x_{2}+\Delta x,\,y_{2}+\Delta
y)\rho_{2}(x_{2},y_{2})dx_{2}\,dy_{2}\,d\Delta x\,d\Delta y=1,$ (1.2)

as can be easily proved by substituting $x_{1}=x_{2}+\Delta x,\
y_{1}=y_{2}+\Delta y$. In the new variables the integrals
$\iint\rho_{1}(x_{1},\,y_{1})dx_{1}\,dy_{1}$ and
$\iint\rho_{2}(x_{2},\,y_{2})dx_{2}\,dy_{2}$ decouple and reduce to unity by
definition.

From section 1.1 and 1.2 we obtain van der Meer formula

$\sigma=\iint\frac{\mu(\Delta x,\,\Delta y)}{N_{1}N_{2}}d\Delta x\,d\Delta y.$
(1.3)

The ratio

$\mu_{sp}=\mu/N_{1}N_{2}.$ (1.4)

is often called the “specific” number of interactions. Their total number
accumulated during the scan, $\iint\mu_{sp}d\Delta x\,d\Delta y$, gives
$\sigma$.

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 $\mu_{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 twoD , 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 $\mu_{sp}$ but not $\mu$ or $N_{1,2}$
separately. Therefore, it remains valid even if $\mu$ and $N_{1}N_{2}$ change
arbitrarily but proportionally during the scan, eg. due to a gradual decrease
of beam currents with time. It is required, however, that $\rho_{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
$\rho_{1}$, so that effectively it becomes constant $\bar{\rho_{1}}$. This
reduces the complicated overlap integral to
$\iint\bar{\rho_{1}}\rho_{2}dx_{2}\,dy_{2}=\bar{\rho_{1}}$, ie. to unity
$\bar{\rho_{1}}=\int\rho_{1}d\Delta x\,\Delta y=1$ if $\rho_{1,2}$ are
normalized. Equivalently, the scan can be viewed from the transverse position
of the first beam ie. in the coordinates $x_{1},\ y_{1}$. Then the first beam
is stationary while the second moves. In this case $\rho_{2}(x_{1}-\Delta
x,y_{1}-\Delta y)$ is “washed out” by $\Delta x,\,\Delta y$ movements, and the
overlap integral reduces to $\int\rho_{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\times 2\times
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, eg. at LHC, is performed not
“dynamically” using the single continuous pass, but stepwise. Ie. the function
$\mu_{sp}(\Delta x,\,\Delta 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.

This beautiful method was invented by van der Meer more than 50 years ago for
the ISR accelerator Carboni:156499 . It was proposed for Sp$\bar{\mathrm{p}}$S
SppS , successfully applied with various modifications at RHIC Drees:1307872 ;
Drees2013AnalysisOV  and LHC White:2010zzc ; lhcb12 ; CMS-PAS-LUM-13-001 ;
Aad:2013ucp ; lhcb14 ; Abelev:2014epa ; Aaboud:2016hhf ; alice:2160174 ; cms15
; cms16 ; cms17 ; atlas19 . The method delivered a record accuracy between
0.7% and a few percent.

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^{-}\to 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 section 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 Ferro-Luzzi:844569 . Here, the profiles $\rho_{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 $\rho_{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 $\rho_{1}\rho_{2}$. The overlap integral is then calculated
analytically from the reconstructed $\rho_{1,2}$ densities. This method is
used up to now only at LHCb lhcb12 ; lhcb14 , which has a dedicated gas
injection system SMOG , an excellent vertex detector and a flexible trigger
suitable for recording beam–gas interactions.

The second method is called the beam–beam imaging Balagura:2011yw . 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 eg. the first
bunch, it effectively becomes a wide and uniform “film” independently of the
initial $\rho_{1}$ distribution. It allows making a “photo” of $\rho_{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 $\rho_{1}$
“filmed” by smeared $\rho_{2}$. In the original van der Meer approach the
smearing allows integrating over $\Delta x\,\Delta y$ and reducing the overlap
integral $\int\rho_{1}\rho_{2}\,dx\,dy$ to unity. In the beam–beam imaging the
same integration is applied to the luminous region profile, ie. to the product
$\rho_{1}\rho_{2}$ not integrated over $dx\,dy$, and allows reconstructing the
individual densities $\rho_{1,2}$.

Up to now, the beam–beam images were taken at LHCb lhcb12  and CMS experiments
Zanetti:1357856 . 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.

###  1.2 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 $\rho_{1,2}$. In more detail this will be discussed in
Sec. 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:

$\rho_{1,2}(x,y)=\rho_{1,2}^{x}(x)\cdot\rho_{1,2}^{y}(y).$ (1.5)

From section 1.1 it follows that the specific number of interactions then also
factorizes, $\mu_{sp}(\Delta x,\Delta y)=\mu^{x}_{sp}(\Delta
x)\cdot\mu^{y}_{sp}(\Delta y)$. This is sufficient to simplify the two-
dimensional integral $\iint\mu_{sp}(\Delta x,\Delta y)d\Delta x\,d\Delta y$
and to reduce it to a product of one-dimensional integrals along the lines
$\Delta x=\Delta x_{0}$ and $\Delta y=\Delta y_{0}$. Indeed,

$\displaystyle\sigma$ $\displaystyle=\iint\mu_{sp}(\Delta x,\Delta y)d\Delta
x\,d\Delta y$ $\displaystyle=\int\mu^{x}_{sp}(\Delta x)d\Delta
x\int\mu^{y}_{sp}(\Delta y)d\Delta y\,\frac{\mu^{y}_{sp}(\Delta
y_{0})\mu^{x}_{sp}(\Delta x_{0})}{\mu^{y}_{sp}(\Delta
y_{0})\mu^{x}_{sp}(\Delta x_{0})}$ $\displaystyle=\frac{\int\mu_{sp}(\Delta
x,\Delta y_{0})\,d\Delta x\times\int\mu_{sp}(\Delta x_{0},\Delta y)\,d\Delta
y}{\mu_{sp}(\Delta x_{0},\Delta y_{0})}.$ (1.6)

The integrals in the enumerator can be measured in two one-dimensional scans
over $\Delta x$ at fixed $\Delta y_{0}$ and vice versa. Note that the formula
is valid for any point $(\Delta x_{0},\Delta y_{0})$. This is rarely stressed
in the literature. It might be advantageous to choose $(\Delta x_{0},\Delta
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 $\mu$ where the derivative of eg.
$\mu_{sp}(\Delta x_{0},\Delta y)$ on $\Delta 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 $\mu$, 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 $\rho_{1,2}$ densities bbicms . 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 $\lesssim 1\%$.

The accuracy can be improved further by performing two-dimensional scans over
the central region giving the dominant contribution to the integral in 1.3.
This approach was pioneered at LHCb in 2017 twoD . Scanning only the central
region was relatively fast but allowed evaluating $\sim 90\%$ of the integral.
The method is model-independent and allows reducing the non-factorization
systematics by an order of magnitude.

###  1.3 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~{}\mu 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 Balagura:2011yw  using two alternative
derivations. In the first, the simple two-dimensional integral in section 1.1
was extended to four dimensions and taken over $\Delta x$, $\Delta y$ in the
same way as 1.3 was obtained from section 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 $\mathbf{v}_{1,2}$ as
shown in Fig. 1.1. They can be decomposed into two parallel
($\mathbf{v}_{1,2}$) and one common perpendicular component
$\mathbf{v}_{\perp}$ with respect to their difference
$\mathbf{v}_{1}-\mathbf{v}_{2}=\Delta\mathbf{v}$.

Figure 1.1: Decomposition of the beam particle velocities $\mathbf{v}_{1,2}$
to the perpendicular $\mathbf{v}_{\perp}$ and the parallel
$\mathbf{v}_{1\parallel}$, $\mathbf{v}_{2\parallel}$ components with respect
to their difference $\mathbf{v}_{1}-\mathbf{v}_{2}=\Delta\mathbf{v}$.

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 $\mathbf{v}_{\perp}$ and then with an arbitrary velocity parallel to
$\Delta\mathbf{v}$. Indeed, in the laboratory frame boosted by
$\mathbf{v}_{\perp}$, or, equivalently, after the “active” boost of the beam
particles by $-\mathbf{v}_{\perp}$, their perpendicular momentum transforms to
$\mathbf{p}_{\perp}^{\prime}=\gamma_{\perp}\mathbf{p}_{\perp}-\gamma_{\perp}\mathbf{\beta}_{\perp}E/c$,
which is zero since by construction
$\mathbf{\beta}_{\perp}=\mathbf{p}_{\perp}c/E$. Here, $c$ is the speed of
light, $E$ is the energy in the laboratory system,
$\mathbf{\beta}_{\perp}=\mathbf{v}_{\perp}/c$ and
$\gamma_{\perp}=(1-\beta_{\perp}^{2})^{-1/2}$ are the beta- and gamma-factors
of the boost, respectively. Therefore, after the active boost by
$-\mathbf{v}_{\perp}$ the beams become parallel to $\Delta\mathbf{v}$ and any
further boost along $\Delta\mathbf{v}$ preserves this parallelism.

Note that although the perpendicular velocity after the boost by
$-\mathbf{v}_{\perp}$ becomes a simple difference
$\mathbf{v}_{\perp}-\mathbf{v}_{\perp}=\mathbf{0}$, the remaining velocity
$\mathbf{v}_{\parallel}^{\prime}$ is, of course, not equal to the difference
$\mathbf{v}-\mathbf{v}_{\perp}=\mathbf{v}_{\parallel}$. 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: $\mathbf{p}_{\parallel}=\mathbf{p}_{\parallel}^{\prime}$,
$E=\gamma_{\perp}E^{\prime}$, so
$\mathbf{v}_{\parallel}^{\prime}=\gamma_{\perp}\mathbf{v}_{\parallel}$.

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
$\mu$, $N_{1,2}$ and $\sigma$ are Lorentz-invariant, only the transverse area
$d\Delta x\,d\Delta y$ is not. Following notations of this subsection, the
latter will be denoted in the primed coordinates as $d\Delta
x^{\prime}\,d\Delta y^{\prime}$. Its transformation to the laboratory system
requires some explanations given below. This material complements the
discussion in Balagura:2011yw .

The space-time coordinates $x=(t,\,\mathbf{x})$ of any particle moving with
the four-momentum $p=(E,\,\mathbf{p})$ satisfy the equation

$x-x^{0}=\lambda p,$ (1.7)

where $\lambda$ is a free parameter and $x^{0}=(t^{0},\,\mathbf{x}^{0})$ is an
arbitrary point on the particle trajectory corresponding to $\lambda=0$. Four
equations in 1.7 with one free parameter define a line in the four-dimensional
space. The values of $\lambda$ uniquely label the line points and can be
expressed, for example, via the time coordinate: $\lambda=(t-t^{0})/E$.

The beam displacements $\Delta=(0,\mathbf{\Delta})$ during van der Meer scan
change the $x^{0}$ parameter:

$x-x^{0}-\Delta=\lambda p.$ (1.8)

Note that the definition of the scan implies that the beams are moved only
spatially, so there is no time component in $\Delta$.

In the following, it will be convenient to decompose spatial vectors into
three components: parallel to $\mathbf{v}_{\perp}$, to the vector-product
$[\Delta\mathbf{v}\times\mathbf{v}_{\perp}]$ and to $\Delta\mathbf{v}$. They
will be denoted by the subscripts $\perp$, $\parallel\times\perp$ and
$\parallel$, respectively, eg.

$\Delta=(0,\,\Delta_{\perp},\Delta_{\parallel\times\perp},\Delta_{\parallel}).$
(1.9)

The $\parallel\times\perp$ component is perpendicular to the plane of Fig.
1.1. Ideally, the beam displacements should be orthogonal to
$\Delta\mathbf{v}$ but we consider below the general case
$\Delta_{\parallel}\neq 0$.

After the active boost by $-\Delta\mathbf{v}$, the line defined by 1.8
transforms to

$x^{\prime}-x^{0\prime}-\Delta^{\prime}=\lambda p^{\prime}.$ (1.10)

Note that the solutions $x$, $x^{\prime}$ of 1.8 and 1.10 for the same
$\lambda$ correspond to the same four-dimensional point in the laboratory and
primed frames, respectively.

After the boost, the beam displacement

$\Delta^{\prime}=(-\beta_{\perp}\gamma_{\perp}\Delta_{\perp},\,\,\gamma_{\perp}\Delta_{\perp},\,\Delta_{\parallel\times\perp},\,\Delta_{\parallel})$
(1.11)

acquires the time component
$\Delta_{t}^{\prime}=-\beta_{\perp}\gamma_{\perp}\Delta_{\perp}$. 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\prime}_{\parallel}+\Delta_{\parallel}^{\prime}$ along the lines or on
the initial time $t^{0\prime}+\Delta_{t}^{\prime}$ since the particles are
anyway assumed to travel from $t^{\prime}=-\infty$ to $t^{\prime}=+\infty$.
Contrary to the transverse initial coordinates, the time and longitudinal
shifts do not change the particle line. Therefore, the scan with
$\Delta^{\prime}$ from 1.11 is equivalent to the one with

$\tilde{\Delta}^{\prime}=(0,\,\gamma_{\perp}\Delta_{\perp},\,\Delta_{\parallel\times\perp},\,0),$
(1.12)

where the time and longitudinal coordinates are simply set to zero. Comparing
this equation with 1.9 one can see that the beam displacements
$\Delta_{\parallel\times\perp}$ transform from the laboratory to the primed
system without any changes while the displacements $\Delta_{\perp}$ are
“extended” by the $\gamma_{\perp}$-factor. Therefore, the area element
$d\Delta x^{\prime}\,d\Delta y^{\prime}$ in the primed coordinates is larger
than in the laboratory system by $\gamma_{\perp}$:

$d\Delta x^{\prime}\,d\Delta y^{\prime}=\gamma_{\perp}d\Delta x\,d\Delta y.$
(1.13)

Note that one can invert the arguments and make the opposite boost from the
primed scan with the displacements
$(0,\Delta^{\prime}_{\perp},\,\Delta^{\prime}_{\parallel\times\perp},\,0)$ to
the laboratory system:

$\Delta=(\beta_{\perp}\gamma_{\perp}\Delta^{\prime}_{\perp},\,\,\gamma_{\perp}\Delta^{\prime}_{\perp},\,\Delta^{\prime}_{\parallel\times\perp},\,0).$

Here, the $\gamma_{\perp}$-factor again appears in the transverse components
but now in the laboratory system. Concluding from this formula that
$\gamma_{\perp}d\Delta x^{\prime}\,d\Delta y^{\prime}=d\Delta x\,d\Delta y$
with $\gamma_{\perp}$ 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
$\Delta_{t}=\beta_{\perp}\gamma_{\perp}\Delta^{\prime}_{\perp}$ 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

$\sigma=\gamma_{\perp}\int\frac{\mu(\Delta x,\Delta y)}{N_{1}N_{2}}d\Delta
x\,d\Delta y.$ (1.14)

Note that the relativistic correction $\gamma_{\perp}$ 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\Delta x\,d\Delta y$ by definition lies in the plane perpendicular
to $\Delta\mathbf{v}$. Let’s denote this plane by $P$. If the scan plane
$\tilde{P}$ is inclined with respect to $P$ at an angle $\alpha$, an area
element $d\tilde{\Delta}x\,d\tilde{\Delta}y$ on $\tilde{P}$ should be
projected to $P$, ie.

$\sigma=\gamma_{\perp}\int\frac{\mu(\tilde{\Delta}x,\tilde{\Delta}y)}{N_{1}N_{2}}\cos\alpha\,d\tilde{\Delta}x\,d\tilde{\Delta}y.$
(1.15)

The longitudinal translations along $\Delta\mathbf{v}$ do not matter.

The beam crossing angles at LHC are small, so $\beta_{\perp}<10^{-3}$ and
$\gamma_{\perp}-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 $\sigma_{iL}$ along the beam
$i=1,2$ in the laboratory frame are 5–10 cm. They are much larger than the
transverse sizes $\sigma_{iT}\sim 100\,\mu m$. Since $\gamma_{\perp}\approx
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 $\sigma_{iL}$ gets
projected to $\mathbf{v}_{\perp}$ at the angles $\alpha_{1,2}$ shown in Fig.
1.1. For the Gaussian bunches this increases the primed transverse width
$\sigma_{i\perp}^{\prime}$ in the beam crossing plane to

$\sigma_{i\perp}^{\prime}\approx\sigma_{iT}\sqrt{1+(\alpha_{i}\,\sigma_{iL}/\sigma_{iT})^{2}}.$
(1.16)

This formula is valid up to the second-order $\alpha_{i}$-terms enhanced by
$\sigma_{iL}/\sigma_{iT}\gg 1$ factor.

Note that often in the literature $\sigma_{i\perp}^{\prime}$ is expressed in
the form of a rotation as
$\sqrt{(\sigma_{iT}\cos\alpha_{i})^{2}+(\sigma_{iL}\sin\alpha_{i})^{2}}$.
Writing $\cos\alpha_{i}=1-\alpha_{i}^{2}/2+\ldots$ instead of unity as in 1.16
implies its validity at least up to the terms $\propto\alpha_{i}^{2}$ not
enhanced by $\sigma_{iL}/\sigma_{iT}$. 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 $\Delta=(0,\mathbf{\Delta})$ 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, $\mathbf{\Delta}$ can be a zero-mean random variable normally and
independently distributed along the transverse $x$, $y$ and longitudinal axes
with the standard deviations $\sigma_{x,y,L}$, respectively. Let’s consider
the general case when neither $x$ nor $y$ lies in the crossing plane of Figure
1.1. Then, the transverse bunch widths in the laboratory frame along $\perp$
and $\parallel\times\perp$ directions are

$\displaystyle\sigma_{i\perp}$
$\displaystyle=\sqrt{\left[\sigma_{ix}\cos(x_{i},\perp)\right]^{2}+\left[\sigma_{iy}\cos(y_{i},\perp)\right]^{2}+\left[\sigma_{iL}\sin\alpha_{i}\right]^{2}}$
$\displaystyle\sigma_{i\parallel\times\perp}$
$\displaystyle=\sqrt{\left[\sigma_{ix}\cos(x_{i},\parallel\\!\times\\!\perp)\right]^{2}+\left[\sigma_{iy}\cos(y_{i},\parallel\\!\times\\!\perp)\right]^{2}}.$
(1.17)

Here, the notation like $\cos(x_{i},\perp)$ denote the cosine of the angle
between the $x$ direction of the $i$-th bunch and $\mathbf{v}_{\perp}$.
According to 1.12, the bunch spread $\sigma_{\perp}$ is multiplied by the
$\gamma_{\perp}$-factor in the primed frame, while
$\sigma_{\parallel\times\perp}$ remains unchanged:

$\sigma_{i\perp}^{\prime}=\gamma_{\perp}\sigma_{i\perp},\quad\sigma_{i\parallel\times\perp}^{\prime}=\sigma_{i\parallel\times\perp}.$
(1.18)

If one of the transverse axes, eg. $x$, lies in the crossing plane, this
simplifies to

$\sigma_{i\perp}^{\prime}=\gamma_{\perp}\sqrt{\left[\sigma_{ix}\cos\alpha_{i}\right]^{2}+\left[\sigma_{iL}\sin\alpha_{i}\right]^{2}},\
\ \ \sigma_{i\parallel\times\perp}^{\prime}=\sigma_{y}.$ (1.19)

For ultrarelativistic beams with $v_{1,2}\approx c$, like at LHC, and
arbitrary $\alpha_{1,2}$, which should be approximately equal in this case,
$\alpha_{1}\approx\alpha_{2}\approx\alpha$, the beta- and gamma-factors can be
expressed via the angle $\alpha$:

$\beta_{\perp}\approx\sin\alpha,\quad\gamma_{\perp}\approx
1/\sqrt{1-\sin^{2}\alpha}=1/\cos\alpha.$ (1.20)

This gives

$\sigma_{i\perp}^{\prime}\approx\sqrt{\sigma_{ix}^{2}+(\sigma_{iL}\tan\alpha)^{2}}.$
(1.21)

One can see that $\sigma_{ix}$ is not multiplied by $\cos\alpha$. Instead of
the rotation formula, for $\sigma_{i\perp}^{\prime}$ one should use either
exact section 1.3, 1.18, 1.19 or one of the approximate equations 1.16, 1.21.

The smallness of $\alpha_{i}$ at LHC is partially compensated by the large
$\sigma_{iL}/\sigma_{iT}$ ratio, so due to the crossing angle the effective
transverse widths $\sigma_{i\perp}^{\prime}$ 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 $\int\rho_{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 $\rho_{1,2}$ factorize in the
primed directions $x^{\prime}$ and $y^{\prime}$ but none of them lies in the
crossing plane of Figure 1.1, $\sigma_{iL}$ has non-zero projections on both
$x^{\prime}$ and $y^{\prime}$. This makes $x^{\prime}$ and $y^{\prime}$
distributions slightly correlated and to some extent breaks the
$x^{\prime}$-$y^{\prime}$ factorizability.

###  1.4 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 $\mu$-values, the accurate linearity is required in larger dynamic
range since van der Meer scans are typically performed with $\mu\lesssim 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 $\sim 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 $\Delta x$ and
$\Delta y$. Any error in the beam displacements directly affects the cross-
section. An accurate measurement of $\Delta x$, $\Delta 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 $\mathbf{\Delta}_{i}$ in the
laboratory frame for the beam $i=1,2$ can be written as

$\mathbf{\Delta}_{i}=\alpha_{i}\mathbf{a}_{i}+\beta_{i}\mathbf{b}_{i}+\mathbf{\Delta}_{i}^{0},$
(1.22)

where $\mathbf{a}_{i}$, $\mathbf{b}_{i}$ are unknown vectors close but not
exactly equal to the unit $x$, $y$-vectors and $\alpha_{i}$, $\beta_{i}$ are
the nominally set values for $x$, $y$ beam movements, respectively. The latter
are known exactly. Potential nonlinearities in the $\mathbf{\Delta}_{i}$
dependence on $\alpha_{i},\beta_{i}$ eg. due to beam orbit drifts, are
neglected here but will be briefly discussed later. The constant vectors
$\mathbf{\Delta}_{i}^{0}$ corresponding to $\alpha_{i}=\beta_{i}=0$ are also
unknown.

In the simplest LSC the beams are nominally displaced by the same amount, ie.
with $\alpha_{1}=\alpha_{2}=\alpha_{LSC}$, $\beta_{1}=\beta_{2}=\beta_{LSC}$.
If the bunches have equal shapes, the center of the luminous region
$\mathbf{O}_{LSC}$ is positioned in the middle between them,

$\mathbf{O}_{LSC}=\alpha_{LSC}\,\frac{\mathbf{a}_{1}+\mathbf{a}_{2}}{2}+\beta_{LSC}\,\frac{\mathbf{b}_{1}+\mathbf{b}_{2}}{2}+\frac{\mathbf{\Delta}_{1}^{0}+\mathbf{\Delta}_{2}^{0}}{2}.$
(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
$\alpha_{1}=-\alpha_{2}=2\alpha^{sym}$, $\beta_{1}=-\beta_{2}=2\beta^{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
$\mathbf{\Delta}_{12}=\mathbf{\Delta}_{1}-\mathbf{\Delta}_{2}$ is then

$\mathbf{\Delta}_{12}=\alpha^{sym}\,\frac{\mathbf{a}_{1}+\mathbf{a}_{2}}{2}+\beta^{sym}\,\frac{\mathbf{b}_{1}+\mathbf{b}_{2}}{2}+\frac{\mathbf{\Delta}_{1}^{0}-\mathbf{\Delta}_{2}^{0}}{2}.$
(1.24)

A comparison of 1.23 and 1.24 shows that the measurements of
$\mathbf{O}_{LSC}$ in the vertex detector are sufficient to calibrate the
scales $(\mathbf{a}_{1}+\mathbf{a}_{2})/2$,
$(\mathbf{b}_{1}+\mathbf{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
$(\mathbf{a}_{1}+\mathbf{a}_{2})/2$ and $(\mathbf{b}_{1}+\mathbf{b}_{2})/2$
measurable in 1.23. The simplest way to calibrate the lengths
$|\mathbf{a}_{1,2}|=a_{1,2}$ and $|\mathbf{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 $\mathbf{a}_{1,2}$ and the
$x$-axis, one arrives at the scalar equations

$\displaystyle\Delta_{i,x}$
$\displaystyle=\alpha_{LSC}\,a_{i}+\Delta_{i,x}^{0},$
$\displaystyle\Delta_{12,x}$
$\displaystyle=\alpha_{LSC}(a_{1}-a_{2})+(\Delta_{1,x}^{0}-\Delta_{2,x}^{0}),$
$\displaystyle O_{LSC,x}$
$\displaystyle=\alpha_{LSC}\,\frac{a_{1}+a_{2}}{2}+\frac{\Delta_{1,x}^{0}+\Delta_{2,x}^{0}}{2}.$
(1.25)

The movement of the luminous region between any two LSC points,
$O^{1}_{LSC,x}-O^{2}_{LSC,x}$, allows to measure the average $x$-scale

$\frac{a_{1}+a_{2}}{2}=\frac{O^{1}_{LSC,x}-O^{2}_{LSC,x}}{\alpha_{LSC}^{1}-\alpha_{LSC}^{2}}.$
(1.26)

The corresponding change of the $x$-distance between the beams
$\Delta_{12,x}^{1}-\Delta_{12,x}^{2}$ can be deduced from the luminosity
change $L_{LSC}^{1}-L_{LSC}^{2}$ 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 $d\Delta_{12,x}/dL=d\alpha^{sym}/dL^{sym}\cdot(a_{1}+a_{2})/2$.
This leads to the equation

$\displaystyle\Delta_{12,x}^{1}-\Delta_{12,x}^{2}$
$\displaystyle=(\alpha_{LSC}^{1}-\alpha_{LSC}^{2})(a_{1}-a_{2})$
$\displaystyle=\frac{d\alpha^{sym}}{dL^{sym}}\cdot\frac{a_{1}+a_{2}}{2}(L_{LSC}^{1}-L_{LSC}^{2}),$
(1.27)

which allows to obtain $(a_{1}-a_{2})/(a_{1}+a_{2})$ from the measurable
values $d\alpha^{sym}/dL^{sym}$, $L_{LSC}^{1}-L_{LSC}^{2}$ and the set
difference $\alpha_{LSC}^{1}-\alpha_{LSC}^{2}$. 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_{LSC}^{1}-L_{LSC}^{2}$, LSC can be performed at a point close to the
maximum of the derivative $dL/d\Delta_{12}$, where the second derivative is
zero. For the Gaussian bunches with the widths $\sigma_{1,2}$, the luminosity
dependence on the beam separation is also Gaussian with the sigma
$\Sigma=\sqrt{\sigma_{1}^{2}+\sigma_{1}^{2}}$, and the optimal LSC beam
separation is $\Delta_{12}=\Sigma$.

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 Sec. 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 twoD .

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 $\rho_{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 $\rho_{1,2}$, therefore,
depends on $\Delta x$, $\Delta y$. For example, it vanishes at large beam
separations. Such $\rho_{1,2}$ dependence on $\Delta x$, $\Delta y$ breaks the
derivation of 1.3 from section 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 proton–proton scans.

The beam–beam interaction also affects the imaging methods. Like the classical
van der Meer approach, the beam–beam image is biased when $\rho_{1,2}$
densities vary during the scan. The beam–gas imaging is also biased when the
beam–beam interaction changes, eg. 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 $\rho_{1,2}$ are constant and the
appropriate $\rho_{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 $\sim 1\%$ error to their luminosity
measurements in the proton–proton scans without making any correction. In 2012
the method madx12  was proposed for correcting the classical van der Meer
technique. It was subsequently used by all LHC experiments, and the systematic
error was reduced to $0.3-0.7\%$ CMS-PAS-LUM-13-001 ; Aad:2013ucp ; lhcb14 ;
Aaboud:2016hhf ; cms15 ; cms16 ; cms17 ; atlas19 . 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 beam–beam corrections
differ by $\sim 1\%$ as shown in Fig. 3.8. 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 non-factorization. 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
$\mu$ 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.

##  2 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 beam–beam 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 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.

###  2.1 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, ie.
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.

Figure 2.1: The electrical field from the charge $q_{1}$ at rest acting on the
charge $q_{2}$ from another bunch.

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
$\mathbf{E}_{1}$ shown in Fig 2.1. 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=\beta_{0}c$ is also constant and $q_{2}$ moves along the
straight line denoted as the $z$-axis in Fig. 2.1.

If the velocities of the colliding particles in the laboratory system are
$\pm\mathbf{\beta}c$, $\beta_{0}$ can be expressed as

$\beta_{0}=2\beta/(1+\beta^{2}),$ (2.1)

which is the double-angle (or double-rapidity) formula for the hyperbolic
tangent, the analog of $\tan(2\phi)=2\tan\phi/(1-\tan^{2}\phi)$. Of course, at
LHC one can safely assume $\beta\approx\beta_{0}\approx 1$.

The momentum kick received by $q_{2}$ is given by the line integral

$\Delta\mathbf{p}_{2}=\int q_{2}\mathbf{E}_{1}^{\perp}(x,y,z)\,dt=\int
q_{2}\mathbf{E}_{1}^{\perp}(x,y,z)\,\frac{dz}{\beta_{0}c}.$ (2.2)

Because of the reflection invariance $z\to-z$, the $z$-projection of the kick
should vanish, so only the perpendicular component of the electric field is
written in 2.2. Its integral $\int\mathbf{E}_{1}^{\perp}(x,y,z)\,dz$ depends
only on the transverse coordinates. It will be denoted in the following simply
by $\mathbf{E}_{1}(x,\,y)$. According to Gauss’s flux theorem applied to the
cylinder with the radius $R$ shown in Fig. 2.1, it is equal to

$E_{1}(x,\,y)=\int E_{1}^{\perp}(x,y,z)\,dz=\frac{q_{1}}{2\pi R\epsilon_{0}},$
(2.3)

where $\epsilon_{0}$ is the electric constant. Here, we are using the system
of units where the first Maxwell’s equation is written as
$\mathbf{\mathbf{\nabla}}\cdot\mathbf{E}=\rho/\epsilon_{0}$, like the
International System of Units (SI).

Substituting 2.3 to 2.2 gives

$\Delta p_{2}=\frac{q_{2}E_{1}(x,\,y)}{\beta_{0}c}=\frac{q_{1}q_{2}}{2\pi
R\epsilon_{0}\beta_{0}c}=\frac{2\alpha\hslash Z_{1}Z_{2}}{R\beta_{0}},$ (2.4)

where $\alpha=e^{2}/(4\pi c\hslash\epsilon_{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 $\hslash$ is the reduced Plank
constant.

Since the momentum kick $\Delta 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 $\alpha_{1,2}$ as in Fig. 1.1, 2.4 receives corrections of the
order $\alpha_{1,2}^{2}$. For example, they can be calculated by boosting the
kick 2.4 from the primed to the laboratory system with the velocity
$\mathbf{v_{\perp}}$ 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 $\sim\\!10^{11}e$ is condensed into $q_{1}$,
the kick at, for example, one bunch sigma $\sim 100\,\mu m$ is equal to
$\Delta 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
$\Delta\mathbf{p}_{1}=-\Delta\mathbf{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 $\Delta p_{1}\sim 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.1 is concentrated only in a thin
“pancake” perpendicular to $z$ and travelling together with the particle. So,
$q_{1}$ receives the kick $\Delta 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

$(t(z)-t(z=0))^{2}=(z/c\beta_{0})^{2}=(R^{2}+z^{2})/c^{2},$

so $q_{2}$ travels the distance $z=R\beta_{0}\gamma_{0}\gg 1$ where
$\gamma_{0}=(1-\beta_{0})^{-1/2}$. Therefore, most of $\Delta p_{2}$ kick is
created by the stationary $q_{1}$.

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/\beta_{0}$ in 2.4 is not justified and should be dropped. The
formula should be written as

$\Delta p=\frac{q_{2}E_{1}(x,\,y)}{c}=\frac{2\alpha\hslash Z_{1}Z_{2}}{R}$
(2.5)

under the explicit condition $\beta_{0}\approx 1$. The corresponding angular
kick for the particles with the momentum $p$ is

$\Delta\phi=\frac{q_{2}E_{1}(x,\,y)}{pc}=\frac{2\alpha\hslash
Z_{1}Z_{2}}{Rp}.$ (2.6)

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\pi\epsilon_{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 10.1093/ptep/ptaa104  it
is defined in Eq. (31.13) as

$\xi_{y2}=\frac{m_{e}r_{e}q_{1}q_{2}N_{1}\beta_{y2}}{2\pi
m_{2}\gamma_{2}\sigma_{y1}(\sigma_{x1}+\sigma_{y1})}$

for the $y$-direction, where $\beta_{y2}$ is the beta-function discussed
later, $m_{2}$, $\gamma_{2}$ are the mass and $\gamma$-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 $\beta_{0}\approx 1$ the
momentum change $\Delta 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
$\alpha$ 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.

###  2.2 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 Sec. 3.1, the angular spread in the
laboratory frame is of the order $\delta\alpha=\sigma_{T}/\beta\lesssim
O(10^{-5})$ where $\sigma_{T}$ is the transverse bunch size ($40-100\,\mu m$)
and $\beta$ 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
$\Delta\mathbf{v}=\mathbf{v}_{1}-\mathbf{v}_{2}$ from Fig. 1.1, so the kick
angular variation is of the same negligible order $\delta\alpha$.

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 $\sigma_{T}^{\prime}$ visible from the
particle according to 1.16. Therefore, the angular variation $\delta\alpha$
leads to the effective smearing of $\sigma_{T}^{\prime}$:

$\frac{\delta\sigma_{T}^{\prime}}{\sigma_{T}^{\prime}}\approx\left(\frac{\sigma_{L}}{\sigma_{T}}\right)^{2}\alpha\
\delta\alpha=\left(\frac{\sigma_{L}}{\sigma_{T}}\right)\left(\frac{\sigma_{L}}{\beta}\right)\alpha.$
(2.7)

In spite of the large enhancement factor $\sigma_{L}/\sigma_{T}\lesssim 1000$,
the values of $\sigma_{L}/\beta\lesssim O(10^{-2})$ and $\alpha\lesssim
O(10^{-4})$ are so small that in van der Meer scans the variations of the
transverse width $\delta\sigma_{T}^{\prime}/\sigma_{T}^{\prime}\lesssim 0.01$
can be neglected in the beam–beam kick calculations.

The longitudinal momentum spread $\delta 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 $\phi$, $v=\tanh\phi$, its change is $\delta
v=\delta\phi/\cosh^{2}\phi=\delta(\sinh\phi)/\cosh^{3}\phi=(\delta
p/p)\cdot\beta/\gamma^{2}\propto 1/\gamma^{2}$. The associated angular
variation of $\Delta\mathbf{v}=\mathbf{v}_{1}-\mathbf{v}_{2}$ due to the
crossing angle is additionally suppressed by the smallness of
$\alpha<10^{-3}$.

Finally, the angular variation due to the beam–beam kick itself is also small,
$\Delta p/p\sim 10^{-6}$. Since the typical longitudinal bunch length
$\sigma_{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 $\sigma_{T}\cdot p/\Delta p\sim 100$ m before
their transverse displacements reach $\sigma_{T}$. However, the accelerator
elements controlling the transverse movements correct the trajectories and
bring the particles back. In Sec. 3.8 equation 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 $\Delta p/2p\lesssim 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
$\rho_{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\pi R$ in the denominator substitutes the sphere
area $4\pi R^{2}$ in the three-dimensional Coulomb’s law in accordance with
the Gauss’s electric flux theorem. Therefore,

$\mathbf{F}_{12}=\Delta\mathbf{p}_{1}c=q_{1}\mathbf{E}_{2}=-q_{2}\mathbf{E}_{1}=-\Delta\mathbf{p}_{2}c$
(2.8)

from 2.5 is just the Coulomb’s two-dimensional force between $q_{1}$ and
$q_{2}$. The calculation of $\mathbf{\Delta}p_{1,2}$ kick, ie. 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.

###  2.3 Electrostatic field from two-dimensional Gaussian distribution

In this subsection we present the formulas of the electrostatic field from the
two-dimensional Gaussian density

$\rho=\frac{Q}{2\pi\sigma_{x}\sigma_{y}}\exp\left(-\frac{x^{2}}{2\sigma_{x}^{2}}-\frac{y^{2}}{2\sigma_{y}^{2}}\right).$
(2.9)

For the round bunch with $\sigma_{x}=\sigma_{y}=\sigma$, the azimuthally
symmetric field can be determined from the charge
$Q\left(1-e^{-R^{2}/2\sigma^{2}}\right)$ inside the disk $x^{2}+y^{2}<R^{2}$
and the Gauss’s flux theorem:

$E=\frac{Q}{2\pi\epsilon_{0}R}\left(1-e^{-R^{2}/2\sigma^{2}}\right).$ (2.10)

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

$\Delta\phi=\frac{Z_{1}eE}{pc}=\frac{2\alpha\hslash
Z_{1}Z_{2}N_{2}}{pR}\left(1-e^{-R^{2}/2\sigma^{2}}\right).$ (2.11)

The field from an elliptical bunch with $\sigma_{x}\neq\sigma_{y}$ is more
complicated. It was derived by Bassetti and Erskine in Bassetti:122227 :

$E_{x}-iE_{y}=\frac{Q\cdot
e^{-z_{2}^{2}}}{\pi\epsilon_{0}\sqrt{2\left(\sigma_{x}^{2}-\sigma_{y}^{2}\right)}}\int_{z_{1}}^{z_{2}}e^{\zeta^{2}}\,d\zeta,$
(2.12)

where the path-independent integral is taken in the complex plane between the
points

$z_{1}=\frac{x\frac{\sigma_{y}}{\sigma_{x}}+iy\frac{\sigma_{x}}{\sigma_{y}}}{\sqrt{2\left(\sigma_{x}^{2}-\sigma_{y}^{2}\right)}},\quad
z_{2}=\frac{x+iy}{\sqrt{2\left(\sigma_{x}^{2}-\sigma_{y}^{2}\right)}}.$ (2.13)

It can be expressed via the complex error function
$\mathrm{erf}(z)=2\int_{0}^{z}e^{-\zeta^{2}}\,d\zeta/\sqrt{\pi}$ or its scaled
version named Faddeeva function

$w(z)=e^{-z^{2}}\left(1+\frac{2i}{\sqrt{\pi}}\int_{0}^{z}e^{\zeta^{2}}d\zeta\right)$
(2.14)

as

$E_{x}-iE_{y}=-iQ\frac{w(z_{2})-w(z_{1})\exp\left(-\frac{x^{2}}{2\sigma_{x}^{2}}-\frac{y^{2}}{2\sigma_{y}^{2}}\right)}{2\epsilon_{0}\sqrt{2\pi(\sigma_{x}^{2}-\sigma_{y}^{2})}}.$
(2.15)

Note that in Bassetti:122227  the sign in front of $y^{2}/2\sigma_{y}^{2}$ 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
$Im(z)$ of its argument tends to $-\infty$. 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
$Im(z_{1,2})$, the calculation can be performed in the following way. In the
case $\sigma_{x}<\sigma_{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 $\sigma_{x}>\sigma_{y}$ and, therefore, are real. Then $Im(z_{1,2})$
becomes negative only if $y<0$. Since the field is centrally symmetric,
$\mathbf{E}(x,y)=-\mathbf{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}$.

###  2.4 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
Hirata:1987cn . Here we give an alternative proof based on simple arguments
and extend the formulas to arbitrary $\rho_{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 $\Delta\mathbf{p}_{ij}$. It
can be calculated as $\Delta\mathbf{p}_{ij}=\mathbf{F}_{ij}/c$ where

$\mathbf{F}_{ij}=\frac{q_{i}q_{j}}{2\pi\epsilon_{0}}\frac{\Delta\mathbf{r}_{ij}}{|\Delta\mathbf{r}_{ij}|^{2}}$
(2.16)

is the two-dimensional Coulomb’s force between the charges $q_{i}$, $q_{j}$
separated by $\Delta\mathbf{r}_{ij}=\mathbf{r}_{i}-\mathbf{r}_{j}$ in the
transverse plane. The full force on the first bunch is the sum
$\sum_{i,j}\mathbf{F}_{ij}$ that can be approximated by the integral

$\mathbf{F}_{1}=\sum_{i,j}\mathbf{F}_{ij}=N_{1}N_{2}\int\mathbf{F}(\mathbf{r}_{1}-\mathbf{r}_{2})\rho_{1}(\mathbf{r}_{1})\rho_{2}(\mathbf{r}_{2})\,d\mathbf{r}_{1}\,d\mathbf{r}_{2},$
(2.17)

where $N_{1,2}$ are the number of particles in the bunches. Since $\mathbf{F}$
depends only on the difference
$\Delta\mathbf{r}=\mathbf{r}_{1}-\mathbf{r}_{2}$, it is convenient to use
$\Delta\mathbf{r}$ as the integration variable

$\mathbf{F}_{1}=N_{1}N_{2}\int\mathbf{F}(\Delta\mathbf{r})\rho(\Delta\mathbf{r})\,d\Delta\mathbf{r},$
(2.18)

where $\rho$ also depends only on $\Delta\mathbf{r}$:

$\rho(\Delta\mathbf{r})=\int\rho_{1}(\Delta\mathbf{r}+\mathbf{r}_{2})\rho_{2}(\mathbf{r}_{2})\,d\mathbf{r}_{2}.$
(2.19)

This is the cross-correlation $\rho_{2}\star\rho_{1}$ or, equivalently, the
convolution $\rho_{1}*\tilde{\rho}_{2}$ where
$\tilde{\rho}_{2}(\mathbf{r})=\rho_{2}(-\mathbf{r})$ ie.
$\rho_{2}(\mathbf{r})$ with an opposite argument.

Equations 2.17 and 2.18 show that $\Delta\mathbf{r}$ spread in the integral
can be equivalently represented either by the two bunch densities $\rho_{1,2}$
or by only one $\rho$. 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 $\rho_{1}$ and $\rho_{2}$ are substituted by the artificial bunch
density $\rho=\rho_{1}*\tilde{\rho}_{2}$ and by the delta-function at zero,
respectively.

This is illustrated schematically in Fig. 2.2. The left picture shows the
overall electrostatic force $\sum_{j}\mathbf{F}_{ij}$ exerted by the second
bunch on the $i$-th particle. The origins of $\mathbf{F}_{ij}$ vectors are
varied according to the density $\rho_{2}(x,y)$. To obtain the full force, one
needs to sum over $i$, ie. to vary the ends of $\mathbf{F}_{ij}$ vectors. Fig.
2.2b shows such variations for the $j$-th particle of the second bunch. Since
parallel translations do not change vectors $\mathbf{F}_{ij}$, the variations
of their ends can be obtained by the opposite variations of the origins, as
shown in the picture. The full sum $\sum_{i,j}\mathbf{F}_{ij}$ in Fig. 2.2c,
therefore, can be calculated by smearing the origins with both probability
densities $\rho_{1}(-\mathbf{r})$ and $\rho_{2}(\mathbf{r})$. Equivalently, it
can be calculated by varying the ends of $\Delta\mathbf{r}$ by
$\rho=\rho_{1}(\mathbf{r})*\rho_{2}(-\mathbf{r})$ in agreement with 2.19. For
the Gaussian bunches the convolution $\rho_{1}*\tilde{\rho}_{2}$ is again
Gaussian with the sigmas
$\Sigma_{x,y}=\sqrt{\sigma_{1x,y}^{2}+\sigma_{2x,y}^{2}}$.

Since the momentum is conserved, the full momentum kicks of two bunches are
opposite, $\mathbf{F}_{1}=-\mathbf{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}\neq
N_{2}$.

Figure 2.2:  The force between the bunches does not change when one bunch
density $\rho_{i}(\mathbf{r})$ ($i=1,2$) is collapsed to the point charge at
the origin while the other is convolved with $\rho_{i}(-\mathbf{r})$. For the
Gaussian $\rho_{1,2}$ densities the convolution has sigma
$\Sigma=\sqrt{\sigma_{1}^{2}+\sigma_{2}^{2}}$.

##  3 Beam–beam numeric simulation B*B

Unfortunately, there is no known analytic method that can predict the
luminosity change caused by the beam–beam effect. Below we describe a new
numerical simulation developed for this purpose named “B*B” or “BxB”
(pronounced “B-star-B”) BxB1 , BxB2 . Before going into details, let’s briefly
remind the transverse dynamics of the particles in an idealized accelerator,
the so-called “betatron motion”.

###  3.1 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

$u^{\prime\prime}+K_{u}(s)u(s)=0,$ (3.1)

where $u$ is the transverse coordinate ($x$ or $y$), $s$ is the circular
coordinate along the ring, $u^{\prime\prime}=\partial^{2}u/\partial 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

$u=\sqrt{\epsilon_{u}\beta_{u}(s)}\cos(\phi_{u}(s)-\phi_{0,u}),$ (3.2)

where $\epsilon_{u}$ is a constant defining the oscillation amplitude and
called “emittance” in the accelerator language, $\beta_{u}(s)$ is the so-
called “beta-function” determined by the equation

$\frac{1}{2}\beta_{u}\beta_{u}^{\prime\prime}-\frac{1}{4}\beta_{u}^{\prime
2}+\beta_{u}^{2}K_{u}=1,$ (3.3)

while

$\phi_{u}(s)=\int_{0}^{s}\frac{d\zeta}{\beta_{u}(\zeta)}.$ (3.4)

is the “phase advance” whose value at $s=0$ is denoted by $\phi_{0,u}$. From
3.2 and 3.4 one can calculate $u^{\prime}=\partial u/\partial s$, ie. the
tangent of the angle between the particle and the orbit:

$u^{\prime}=\sqrt{\frac{\epsilon_{u}}{\beta_{u}}}\left(\frac{\beta_{u}^{\prime}}{2}\cos(\phi_{u}-\phi_{0,u})-\sin(\phi_{u}-\phi_{0,u})\right).$
(3.5)

At the interaction point the beams are maximally focused to reach the maximal
luminosity. As it follows from 3.2, $\beta_{u}$ is then minimal and
$\beta_{u}^{\prime}=0$, so that 3.5 simplifies to

$u^{\prime}=-\sqrt{\frac{\epsilon_{u}}{\beta_{u}}}\sin(\phi_{u}-\phi_{0,u}).$
(3.6)

After every turn in the accelerator the particle phase advance increases by
the constant

$Q_{u}=\frac{1}{2\pi}\int_{0}^{L}\frac{d\zeta}{\beta_{u}(\zeta)}$ (3.7)

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^{\prime})$ to one
complex variable

$z_{u}=u-iu^{\prime}\beta_{u}.$ (3.8)

Its evolution is described by the simple rotation in the complex plane

$z_{n+1,u}=\sqrt{\epsilon_{u}\beta_{u}}e^{i(\phi_{n+1,u}-\phi_{0,u})}=z_{n,u}e^{2\pi
iQ_{u}},$ (3.9)

where $z_{n+1,u}$ and
$z_{n,u}=\sqrt{\epsilon_{u}\beta_{u}}e^{i(\phi_{n,u}-\phi_{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-ix^{\prime}\beta_{x}$, $z_{y}=y-iy^{\prime}\beta_{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^{\prime}\beta_{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^{\prime}$ 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\beta_{u}\Delta u^{\prime})e^{2\pi iQ_{u}}.$ (3.10)

According to 2.8, the angular kicks in the first bunch, for example, are
determined by $(\Delta x^{\prime},\,\Delta
y^{\prime})=q_{1}\mathbf{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

$\int(\rho_{1}+\delta\rho_{1})\rho_{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
$\int\rho_{1}(\rho_{2}+\delta\rho_{2})\,dy$. The full luminosity change with
respect to the unperturbed value is approximated as

$\displaystyle\int\left(\rho_{1}+\delta\rho_{1}\right)\left(\rho_{2}+\delta\rho_{2}\right)dx\,dy-\int\rho_{1}\rho_{2}dx\,dy$
$\displaystyle\approx\int\left(\delta\rho_{1}\cdot\rho_{2}+\rho_{1}\cdot\delta\rho_{2}\right)dx\,dy.$
(3.12)

The second-order term $\int\delta\rho_{1}\delta\rho_{2}dx\,dy$ is neglected.
If the bunches are identical, two terms in section 3.1 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, ie. 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.

###  3.2 Particle weights

In section 3.1 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, eg. $\rho_{1}+\delta\rho_{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,\ldots$ . Each of
them can be assigned to one “macro-particle” with the weight $w_{i}$ equal to
the phase-space density integrated over $V_{i}$ ,
$w_{i}=\int_{V_{i}}(\rho_{1}+\delta\rho_{1})d^{4}V$, where
$d^{4}V=dx\,dx^{\prime}\,dy\,dy^{\prime}$. 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

$\displaystyle\int(\rho_{1}+\delta\rho_{1})\rho_{2}dx\,dy\approx$
$\displaystyle\sum_{i=1}^{N}\rho_{2}(x_{i},y_{i})\int_{V_{i}}(\rho_{1}+\delta\rho_{1})d^{4}V$
$\displaystyle=$ $\displaystyle\sum_{i=1}^{N}w_{i}\rho_{2}(x_{i},y_{i}).$
(3.13)

Here, the density $\rho_{2}$ is considered constant in $x$-$y$ projection of
each $V_{i}$ and substituted by its value $\rho_{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_{1}^{B*B}\ll 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}\int_{V_{i}}(\rho_{1}+\delta\rho_{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 $\rho_{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

$r_{x}^{i}=(n_{x}^{i}-0.5)\Delta_{x},\quad
r_{y}^{i}=(n_{y}^{i}-0.5)\Delta_{y},$ (3.14)

where $n_{x,y}^{i}=1,2,\ldots,n_{max}$. The index $i$ runs across all
$1,2,\ldots,n_{max}^{2}$ radii pairs $(r_{x},r_{y})$, where

$n_{max}^{2}=N_{part}\sim O(1000)$ (3.15)

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\pi 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_{turn}^{BB}\sim 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 3.10 is applied approximately $N_{part}\times
N_{turn}^{BB}=O(10^{6}-10^{7})$ times and all calculated coordinates
contribute to the sampling and the final integral. The luminosity sum over all
particles section 3.2 is calculated in B*B after every turn. The average gives
the final result.

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_{x}^{i},r_{y}^{i})$. 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
$z_{u}=u-i\beta_{u}u^{\prime}$,

$\rho_{1}^{u}(z_{u})=\frac{1}{2\pi\sigma_{u}^{2}}\exp\left(-\frac{|z_{u}|^{2}}{2\sigma_{u}^{2}}\right),$
(3.16)

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_{u}^{i}\Delta_{u}<r_{u}<(n_{u}^{i}+1)\Delta_{u}$:

$w_{i}^{u}=\iint_{i}\rho_{1}^{u}d\phi_{u}\,dr_{u}\approx\frac{1}{\sigma_{u}^{2}}\exp\left(-\frac{(r_{u}^{i})^{2}}{2\sigma_{u}^{2}}\right)r_{u}^{i}\Delta_{u}.$
(3.17)

The final weight of the particle placed at the radii $(r_{x}^{i},r_{y}^{i})$
is the product

$w_{i}=w_{i}^{x}w_{i}^{y}.$ (3.18)

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 $\rho_{1,2}$
as the sum of an arbitrary number of Gaussians with configurable weights and
widths and independently in $x$ and $y$. The field $\mathbf{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 $\sigma_{x}\neq\sigma_{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
$\rho_{2}(x_{i},y_{i})$ at every turn appearing in section 3.2 are calculated
as weighted sums of the contributions from all Gaussians.

The radii limits $n_{max}\Delta_{x}$, $n_{max}\Delta_{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_{\sigma}\sigma_{x}$ and
$r_{y}<N_{\sigma}\sigma_{y}$ cuts for the simple single Gaussian shape. The
parameter $N_{\sigma}$ is configurable. Its default value $N_{\sigma}=5$ is
usually sufficient to reach the required accuracy. For the single Gaussian
bunch the excluded weight, ie. the density integral of the not simulated
region $(r_{x}>5\sigma_{x}\ \text{or}\ r_{y}>5\sigma_{y})$, is
$1-(1-e^{-5^{2}/2})^{2}=7.5\cdot 10^{-6}$. To additionally reduce the CPU
time, the largest $(r_{x},r_{y})$ pairs are not simulated, namely, it is
required that

$\left[\frac{r_{x}}{n_{max}\Delta_{x}}\right]^{2}+\left[\frac{r_{y}}{n_{max}\Delta_{y}}\right]^{2}<1.$
(3.19)

This increases the lost weight to $(5^{2}/2+1)e^{-5^{2}/2}=5.0\cdot 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 madx12
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 $\leq 4\cdot 10^{-5}$
in the full simulated range of bunch separations from $0$ to $200\ \mu$m.

Table 1: Bunch parameters from madx12  describing one of van der Meer calibration scans in ATLAS. p | $Z_{1,2}$ | $\beta_{x,y}$ | tune $Q_{x},Q_{y}$ | bunch $\sigma$ | $N_{1,2}$
---|---|---|---|---|---
3500, GeV | 1 | 1.5 m | 64.31, 59.32 | 40 $\mu$m | $8.5\cdot 10^{10}$

###  3.3 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 section 3.2 partially cancels
in the ratio

$R=\frac{\int(\rho_{1}+\delta\rho_{1})\rho_{2}dx\,dy}{\int\rho_{1}\rho_{2}dx\,dy}$
(3.20)

if both integrals are taken numerically in the same way. Therefore, B*B starts
from simulating $N_{turn}^{no\ BB}$ accelerator turns without the beam–beam
interaction ie. with the beam–beam kick $\Delta u^{\prime}=0$ in 3.10. The
unperturbed luminosity $\int\rho_{1}\rho_{2}dx\,dy$ is estimated turn-by-turn
using section 3.2. The default value of $N_{turn}^{no\ BB}$ is conservatively
chosen to be 1000. In Sec. 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_{turn}^{no\ BB}=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
$\int\rho_{1}\rho_{2}dx\,dy$ integrals less than $2\cdot 10^{-4}$ in the
practically important region of the beam separations contributing $\sim
99.9\%$ to the cross-section integral in 1.3. The resulting bias of the ratio
$R$ in 3.20 should, therefore, be smaller.

After the first $N_{turn}^{no\ BB}$ turns the beam–beam kick is switched on.
The user has two options: either instantaneously apply the nominal kick
$\Delta u^{\prime}$ in 3.10 or increase it slowly or “adiabatically”, namely,
linearly from zero to the nominal value during $N_{turn}^{adiab}$ 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. 3.1. 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
1000-th 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^{\prime}$. This will be discussed in more detail in section 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.

Figure 3.1: The trajectories of the particles with $r_{x}=20$, $r_{y}=30\
\mu$m in $z_{x}$ projection for the bunches separated in $x$ by 40 $\mu$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\ \mu$m is shown by
the green dashed line.

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. 3.1. 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 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 $\mu$m in Fig. 3.1. Therefore, the contribution to the
overlap integral in section 3.2 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. 3.2 by the blue
and red points for various beam separations.

Figure 3.2: The deviation of the ratio $R$ defined in 3.20 from unity in per
mille in the $x$-scan versus the separation of the bunches expressed in the
bunch widths. The beam parameters are taken from Table 1. The error bars show
the standard deviation of the results when the simulation was repeated with
different random seeds 100 and 25 times for instantaneous (blue) and adiabatic
(red points) switch of the beam–beam force, respectively. The number of
simulated turns are $[N_{turn}^{no\
BB},N_{turn}^{adiab},N_{turn}^{BB}]=[100,0,500]$ (blue), $[100,100,500]$
(red), $[10\,000,10\,000,1\,000\,000]$ (black). The red and blue points are
displaced horizontally to the left and right, respectively, to reduce
overlapping.

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. 3.2 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. 3.2,
simulating only 500 turns already gives sufficient accuracy. The default B*B
values, $N_{turn}^{no\ BB}=1000$, $N_{turn}^{adiab}=1000$ and
$N_{turn}^{BB}=5000$, are, therefore, quite conservative and may be reduced by
a factor of 10 in practice.

After the beam–beam force is fully switched on and before calculating the
overlap integral over $N_{turn}^{BB}$ turns in the final stage, the B*B
simulation has an option to run $N_{turn}^{stab}$ 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_{turn}^{stab}$ parameter exists only for flexibility and for experimenting
with the B*B simulation.

###  3.4 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 $\phi^{0}_{x,y}$ is described by 3.9:

$x_{n}=r_{x}\cos(2\pi Q_{x}n+\phi^{0}_{x}),\ \ y_{n}=r_{y}\cos(2\pi
Q_{y}n+\phi^{0}_{y}),$ (3.21)

ie. appears to be a Lissajous curve. The beam–beam force leads to the
dispersion of the trajectory as shown in Fig. 3.3 left by the black points for
the artificial values $Q_{x}=0.3$, $Q_{y}=1/3$ and 5000 accelerator turns.

Figure 3.3: Left: The trajectory of one particle with $r_{x}=20$, $r_{y}=30$
$\mu$m followed during 5000 turns in the presence of the beam–beam force. The
bunch parameters are from Table 1 except the tunes. The black points are
obtained with $Q_{x}=0.3$, $Q_{y}=1/3$. The red points, scattered in the
background, are simulated with $Q_{x}$, $Q_{y}$ shifted by
$+e^{-8.5}/4,-e^{-8.5}/4$, respectively. Right: the probability density of
$x$-$y$ localization of the particle followed during $10^{6}$ turns.

In general, for the rational $Q_{x,y}$ and a small beam–beam force the
trajectory “cycles” after $LCD(Q_{x},Q_{y})$ turns where $LCD$ denotes the
lowest common denominator. In the above example $LCD(3/10,1/3)=30$, so after
one “cycle” with 30 turns the phase increments $30\cdot 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 beam–beam 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 $\delta Q_{x,y}$ randomly distributed in the
interval $[-\epsilon,\epsilon]$ are added to $Q_{x}$ and $Q_{y}$ by default,
where $\epsilon$ is chosen to be $e^{-8.5}/2=1.01734\ldots\times 10^{-4}$. As
an example, the points obtained with the shifts $\delta Q_{x}=\epsilon/2$,
$\delta Q_{y}=-\epsilon/2$ in 5000 turns are shown in Fig. 3.3 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. 3.3 right for
$N_{turn}^{BB}=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$,
$\rho_{r_{x},r_{y}}(x,y)=\rho_{r_{x}}(x)\rho_{r_{y}}(y)$. The terms
$\rho_{r_{u}}(u)$, where $u=x,y$, can be derived by projecting the uniform
density $d\phi/2\pi$ of the $z_{u}$-circle with the radius $r_{u}$ to the
$u$-axis:

$\rho_{r_{u}}(u)=\frac{1}{2\pi}\left|\frac{d\phi}{du}\right|=\frac{1}{2\pi
r_{u}\left|\frac{d\cos\phi}{d\phi}\right|}=\frac{1}{2\pi\sqrt{r_{u}^{2}-u^{2}}}.$
(3.22)

Therefore,

$\rho_{r_{x},r_{y}}(x,y)=\frac{1}{4\pi^{2}\sqrt{(r_{x}^{2}-x^{2})(r_{y}^{2}-y^{2})}}.$
(3.23)

It is interesting that the smooth two-dimensional $x$-$y$ Gaussian shape of
any bunch is always intrinsically composed from such peaking
$\rho_{r_{x},r_{y}}$ densities. As follows from 3.17, they should be taken
with the weights
$\propto\exp(-r_{x}^{2}/2\sigma_{x}^{2}-r_{y}^{2}/2\sigma_{y}^{2})dr_{x}^{2}\,dr_{y}^{2}$.
Conversely, only the distributions decomposable to $\rho_{r_{x},r_{y}}$ can
represent the bunch shapes $\rho_{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 $LCD(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 $\delta Q\sim 10^{-4}$ mentioned above introduces a
sizeable spread to the trajectory only after $O(10^{3})$ turns. Fig. 3.2
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 $\delta
Q_{x,y}\sim 10^{-4}$ does not matter. The luminosity bias due to modifications
of $Q_{x,y}$ is also negligible, it is less than $5\cdot 10^{-5}$ for all
bunch separations. Therefore, $\delta Q_{x,y}$ plays no role here but might be
useful in the cases when $LCD(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 $\sim O(100)$, it is important to choose them as
multiples of $LCD(Q_{x},Q_{y})$. For the LHC tunes, they should be multiples
of 100, like $N_{turn}^{no\ BB}=100$, $N_{turn}^{BB}=500$ in Fig. 3.2. 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 3.4:  $(R-1)\times 10^{3}$ in the $x$-scan versus the “cycle” number,
where $R$ is defined in 3.20 and averaged per “cycle”, ie. within non-
intersecting windows of 100 accelerator turns. Green and yellow bands
delimited by the two vertical dashed lines correspond to $N_{turn}^{no\
BB}=1000$ first turns without the beam–beam force and next
$N_{turn}^{adiab}=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.

Fig. 3.4 shows the luminosity $R$-correction from 3.20 in the $x$-scan
calculated “cycle-by-cycle”, ie. 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. 3.2 at higher bins. They show the level of statistical
fluctuations of $R$-averages over $N_{turn}^{BB}=500$ turns, ie. over the
first 5 beam–beam “cycles” to the right from the second vertical dashed line
in Figure 3.4.

###  3.5 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 $\beta^{i}_{x,y}$ and the constant phase advances normalized by
$2\pi$

$Q^{i}_{x,y}=\frac{1}{2\pi}\int_{s_{1}}^{s_{i}}\frac{d\zeta}{\beta_{x,y}(\zeta)}$
(3.24)

between the first and $i$-th points. By definition, the first point has
$Q^{1}_{x,y}=0$. The emittance $\epsilon_{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_{u}^{i}=\sqrt{\epsilon_{u}\beta_{u}^{i}}\propto\sqrt{\beta_{u}^{i}}$. The
recurrence relations then take the form

$\displaystyle z_{n+1,u}^{2}$
$\displaystyle=\left(z_{n,u}^{1}-i\beta_{u}^{1}\Delta u^{\prime
1}\right)e^{2\pi
i\left(Q_{u}^{2}-Q_{u}^{1}\right)}\sqrt{\beta_{u}^{2}/\beta_{u}^{1}},$ (3.25)
$\displaystyle z_{n+1,u}^{3}$
$\displaystyle=\left(z_{n,u}^{2}-i\beta_{u}^{2}\Delta u^{\prime
2}\right)e^{2\pi
i\left(Q_{u}^{3}-Q_{u}^{2}\right)}\sqrt{\beta_{u}^{3}/\beta_{u}^{2}},$
$\displaystyle\ldots$ $\displaystyle z_{n+1,u}^{1}$
$\displaystyle=\left(z_{n,u}^{N}-i\beta_{u}^{N}\Delta u^{\prime
N}\right)e^{2\pi
i\left(Q_{u}-Q_{u}^{N}\right)}\sqrt{\beta_{u}^{1}/\beta_{u}^{N}}.$

Let’s consider, for example, the last equation. The kick $\Delta u^{\prime N}$
is determined by the electrostatic field at the last interaction point $N$ at
the position $(x^{N},y^{N})=(Re(z_{x}^{N}),Re(z_{y}^{N}))$. The factor
$e^{2\pi i\left(Q_{u}-Q_{u}^{N}\right)}$, where $Q_{u}$ is the full tune,
rotates the phase, while the term $\sqrt{\beta_{u}^{1}/\beta_{u}^{N}}$ 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\to(N-1)\to\ldots\to 1\to 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 $\sqrt{s}=13$ TeV are listed in Table 2.

Table 2: The phase advances $Q_{x,y}^{i}$ 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 $\sqrt{s}=13$ TeV. Last column lists the beta-functions $\beta_{x}=\beta_{y}$ from madx12  used in Figure 3.5 simulation. | $Q_{x}^{\text{beam}\ 1}$ | $Q_{x}^{\text{beam}\ 2}$ | $Q_{y}^{\text{beam}\ 1}$ | $Q_{y}^{\text{beam}\ 2}$ |  $\beta_{x,y}$, m
---|---|---|---|---|---
ATLAS | 0 | 0 | 0 | 0 | 1.5
ALICE | 8.2960 | 8.2728 | 7.6692 | 7.9577 | 10
CMS | 31.9757 | 31.9844 | 29.6486 | 29.7613 | 1.5
LHCb | 56.0648 | 55.7990 | 51.0171 | 51.7158 | 3

Figure 3.5 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.

Figure 3.5: $(R-1)\times 10^{3}$ versus the beam separation expressed in bunch
widths. The facets correspond to the $x$-scans at four LHC interaction points
$i=1,2,3,4$. The solid black (colored dotted) lines with the smaller open
circles correspond to the correction $R_{i,i}$ ($R_{i,j}$, $j\neq i$) when the
beam–beam force is switched on only at the scanned LHC interaction point $i$
(one of three other points $j\neq i$). The case when the beam–beam force is
switched on at all four points, $R_{i,1234}$, is shown by the blue solid
curves with the larger solid circles. The purple dashed lines with the larger
open circles are their approximations by the four-products $R_{i,1}\cdot
R_{i,2}\cdot R_{i,3}\cdot R_{i,4}$. The beams at not scanned interaction
points collide head-on. The simulation parameters are listed in Tables 1 and
2.

The $\beta$-function values for Figure 3.5 simulation are listed in Table 2.
Together with the ATLAS beam parameters from Table 1 they are taken from
madx12 . The number of particles and emittances of all bunches are the same.
The emittances are determined from the bunch width and $\beta$-value at ATLAS.
Four plots in Fig. 3.5 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 Figure 3.5, 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 exceed $4\cdot
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 $\beta_{x,y}$, $p$, $\epsilon$ and $N$ scalings

As one can see from Table 2, the beta-functions $\beta^{i}_{x,y}$ and,
therefore, the bunch widths $\sigma_{i}\propto\sqrt{\beta^{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. 3.5 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. 3.2. Let’s understand why this happens despite different bunch widths and
kicks.

Scaling of both $\beta_{x}$ and $\beta_{y}$ by a constant factor
$\beta_{x,y}\to\alpha\beta_{x,y}$ changes the linear scale in the $z_{x,y}$
planes by $\sqrt{\alpha}$. The two-dimensional electrostatic field from a
point charge $q$ drops with the distance as $\propto 1/R$, so the angular kick
changes as $\Delta u^{\prime}\to\Delta u^{\prime}/\sqrt{\alpha}$. According to
the definition 3.8, in the $z_{u}$ complex planes the angular kick is
additionally multiplied by $\beta$, so $\Delta z_{u}=-i\beta_{u}\Delta
u^{\prime}$ scales as $\Delta z_{u}\to\sqrt{\alpha}\Delta z_{u}$. Therefore,
the distributions of two bunches and the kicks $\Delta z_{x,y}$ scale
proportionally to $\sqrt{\alpha}$. New particle trajectories can be obtained
by the simple $\sqrt{\alpha}$-scaling of the whole $z_{x,y}$ planes. This
modifies the luminosities by $\alpha$ but keeps constant their ratios, eg. $R$
in 3.20.

To stress this invariance, it is better to rewrite 3.25 via the ratios
$z_{u}/\sqrt{\beta_{u}}$:

$\frac{z_{n+1,u}^{"i+1"}}{\sqrt{\beta_{u}^{"i+1"}}}=\left(\frac{z_{n,u}^{i}-i\beta_{u}^{i}\Delta
u^{\prime i}}{\sqrt{\beta_{u}^{i}}}\right)e^{2\pi
i\left(Q_{u}^{"i+1"}-Q_{u}^{i}\right)},$ (3.26)

where $"i+1"$ index denotes the next interaction point, eg. the first after
the last, while $Q_{u}^{"N+1"}$ is the full tune $Q_{u}$.

Variables in the complex planes $z_{u}/\sqrt{\beta_{u}}$ have another
advantage: they do not change across the interaction points in the ideal
linear accelerator. Indeed, both the phase $Arg(z_{u})=\phi_{u}$ and the
absolute value $\left|z_{u}/\sqrt{\beta_{u}}\right|=\sqrt{\epsilon_{u}}$ are
invariant, since the emittance $\epsilon_{u}$ is conserved. So, the complex
planes $z_{x}/\sqrt{\beta_{x}}$ and $z_{y}/\sqrt{\beta_{y}}$ are common to all
interaction points. As we have seen, the beam–beam kick $-i\beta_{u}^{i}\Delta
u^{\prime i}/\sqrt{\beta_{u}^{i}}$ in 3.26 happen to be invariant to the
scaling $\beta_{x,y}^{i}\to\alpha^{i}\beta_{x,y}^{i}$ with arbitrary
$\alpha_{i}$ because of $\propto 1/R$ Coulomb’s law in two dimensions.
Therefore, the trajectories in the $z_{u}/\sqrt{\beta_{u}}$ planes drawn from
the fixed initial points $z_{u}^{1}/\sqrt{\beta_{u}}$ using 3.26 do not depend
at all on $\beta_{u}^{i}$ values. This is true under the condition that
$\beta_{x}^{i}$ and $\beta_{y}^{i}$ scale proportionally but independently in
different experiments. If $\beta_{x}^{i}/\beta_{y}^{i}$ ratio changes at any
point $i$, this also changes the relative $x^{i}$ and $y^{i}$ scales in
$R^{i}=\sqrt{(x^{i})^{2}+(y^{i})^{2}}$ and in the Coulomb’s law, and breaks
the kick invariance. At LHC, however, typically at all interaction points the
two beta-functions are equal, $\beta_{x}^{i}=\beta_{y}^{i}$.

The horizontal axes in Figures 3.2 and 3.5 are also chosen in the form of the
scale-invariant ratio. Therefore, all $R_{i,i}$-correction curves shown in
these figures do not depend on the specific $\beta_{x}^{i}=\beta_{y}^{i}$
values from Tables 1 and 2. They are determined only by the distribution of
the emittances $\epsilon_{u}$ or the areas of the $z_{u}$-circles $\pi
r_{u}^{2}=\pi\epsilon_{u}$. For example, $\beta^{i}=1$ m at all four
experiments leads to the same figures.

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\neq j$ curves in
Figure 3.5 are different.

There is another interesting invariance of the beam–beam perturbations. During
acceleration of the particles, their angles $u^{\prime}$ 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

$\epsilon_{u}=\tilde{\epsilon}_{u}\cdot m/p,$ (3.27)

where $m$ is the mass of the particle and $\tilde{\epsilon}_{u}$ is the
constant normalized emittance. Therefore, when $p$ increases, the complex
plane $z_{u}/\sqrt{\beta_{u}}$ shrinks as $1/\sqrt{p}$ and the beam–beam force
increases as $\sqrt{p}$. The angular beam–beam kick $\Delta p/p$ contains $p$
in the denominator as in 2.6, and decreases as $1/\sqrt{p}$, ie. coherently
with the $z_{u}/\sqrt{\beta_{u}}$ complex plane. Therefore, in the complex
plane corresponding to the normalized emittances
$z_{u}\sqrt{p/m\beta_{u}}=\sqrt{\tilde{\epsilon}_{u}}e^{i\phi_{u}}$ the
trajectories remain invariant if the phase advances are constant and
$\tilde{\epsilon}_{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 $\alpha$ increases
the linear scale by $\sqrt{\alpha}$ and reduces the beam–beam angular kick by
$1/\sqrt{\alpha}$. In order to have the same $\sqrt{\alpha}$ change in the
linear and angular scales for the first bunch, the beam–beam force from the
second should be enhanced by $\alpha$. This can be achieved by increasing
$Z_{1}Z_{2}N_{2}$. Therefore, the simultaneous scaling of
$\tilde{\epsilon}_{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

$z_{1u}\sqrt{\frac{p_{1}}{m_{1}\beta_{u1}Z_{1}Z_{2}N_{2}}}=\sqrt{\frac{\tilde{\epsilon}_{1u}}{Z_{1}Z_{2}N_{2}}}e^{i\phi_{1u}},$
(3.28)

the results depend only on the initial distributions, the phase advances and
$\beta_{1x}/\beta_{1y}$ ratios, but do not depend on the individual values of
$\beta_{1x,1y}$, $\epsilon_{1u}$, $p_{1}$, $Z_{1,2}$ or $N_{2}$. Only their
combination

$\frac{\tilde{\epsilon}_{1u}}{Z_{1}Z_{2}N_{2}}=\frac{\epsilon_{1u}p_{1}}{Z_{1}Z_{2}N_{2}m_{1}}$
(3.29)

matters. This can also be confirmed with the B*B simulation.

###  3.7 Simulation of the beam crossing angle

If the beams collide with the crossing angle as in Figure 1.1, the betatron
oscillations occur in the planes transverse to the beam vectors
$\mathbf{v}_{1,2}$. They determine the transverse widths $\sigma_{iT}$ of the
bunches $i=1,2$. The luminosity, however, depends on the bunch widths
$\sigma_{iT}^{\prime}$ perpendicular to $\Delta\mathbf{v}$. According to 1.21,
they get additional contributions from the longitudinal widths $\sigma_{iL}$.

The beam–beam kick is also perpendicular to $\Delta\mathbf{v}$ either in the
laboratory or, for example, in the rest frame of the first bunch particle
$q_{1}$ shown in Figure 2.1. Therefore, $q_{1}$ “sees” the opposite bunch
width $\sigma_{2T}^{\prime}$ enhanced by $\sigma_{2L}$, and this value should
be used in the B*B simulation for $\rho_{2}$. The bunch creating the field is
“static” in the B*B model, its betatron transverse motion and $\sigma_{2T}$ do
not matter.

However, to simulate the betatron trajectories in the first bunch, the initial
widths $\sigma_{1x}$, $\sigma_{1y}$ should be specified in B*B without the
longitudinal component. The betatron oscillations are insensitive to
$\sigma_{1L}$ spread. Since the crossing angles are small at LHC,
$\alpha_{1,2}=\alpha<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
$\sigma_{1T},\sigma_{2T}^{\prime}$ instead of the required
$\sigma_{1T}^{\prime},\sigma_{2T}^{\prime}$. Therefore, the final density
$\rho_{1}+\delta\rho_{1}$, corresponding to $\sigma_{1T}$, should be
additionally smeared by the contribution from $\sigma_{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 $\beta_{1}$ the angle between the $x$-axis and the projection of the
crossing plane to the $x$-$y$ plane. The perturbed density
$\rho_{1}+\delta\rho_{1}$ should then be convolved with the two-dimensional
Gaussian $G(x,y)$ with the sigmas $\alpha\sigma_{1L}\cos\beta_{1}$ and
$\alpha\sigma_{1L}\sin\beta_{1}$ in $x$ and $y$, respectively. If, for
example, the $x$-axis belongs to the crossing plane and $\beta_{1}=0$, the
smearing occurs only along $x$ with the sigma $\alpha\sigma_{1L}$, while
$G(x,y)$ projection to $y$ reduces to the delta-function. Instead of smearing
$\rho_{1}+\delta\rho_{1}$ in the overlap integral, it is simpler to perform an
equivalent smearing of $\rho_{2}$:

$\displaystyle L$
$\displaystyle\propto\int((\rho_{1}+\delta\rho_{1})*G)\,\rho_{2}\,d^{2}\mathbf{r}$
$\displaystyle=\int\left(\int(\rho_{1}+\delta\rho_{1})(\mathbf{r}^{\prime})G(\mathbf{r}-\mathbf{r}^{\prime})d^{2}\mathbf{r}^{\prime}\right)\rho_{2}(\mathbf{r})d^{2}\mathbf{r}$
$\displaystyle=\int(\rho_{1}+\delta\rho_{1})(\mathbf{r}^{\prime})\left(\int
G(\mathbf{r}^{\prime}-\mathbf{r})\rho_{2}(\mathbf{r})d^{2}\mathbf{r}\right)d^{2}\mathbf{r}^{\prime}$
$\displaystyle=\int(\rho_{1}+\delta\rho_{1})\,(G*\rho_{2})\,d^{2}\mathbf{r}^{\prime},$
(3.30)

since $G(\mathbf{r}-\mathbf{r}^{\prime})=G(\mathbf{r}^{\prime}-\mathbf{r})$,
where $\mathbf{r}$ and $\mathbf{r}^{\prime}$ denote the two-dimensional
vectors in the $(x,y)$ plane.

To implement this scheme, the B*B simulation takes the values

$\Delta\sigma_{1x}^{i}=\alpha\sigma_{1L}\cos\beta_{1},\quad\Delta\sigma_{1y}^{i}=\alpha\sigma_{1L}\sin\beta_{1}$
(3.31)

as configurable parameters defining the Gaussian $G^{i}(x,y)$ for each
interaction point $i$. Its convolution with the multi-Gaussian $\rho_{2}$ is
performed analytically. The result is used in section 3.2 instead of
$\rho_{2}$ for the overlap integral calculation. Note that it receives
contributions from both $\sigma_{1L}$ and $\sigma_{2L}$. The field is
determined from $\rho_{2}$ not smeared by $G^{i}(x,y)$, where only
$\sigma_{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.

###  3.8 Cross-checks and comparison with the old model used at LHC in
2012–2019

In Section 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
$\rho_{2}$ to $\tilde{\rho}_{1}*\rho_{2}$. In case of the Gaussian bunches
with the widths $\sigma_{1,2;x,y}$ and the centers $\mathbf{r}^{0}_{1,2}$, the
cross-correlation $\tilde{\rho}_{1}*\rho_{2}$ is also Gaussian and has the
widths $\Sigma_{x,y}=\sqrt{\sigma_{1;x,y}^{2}+\sigma_{2;x,y}^{2}}$ and the
center $\Delta\mathbf{r}_{21}^{0}=\mathbf{r}^{0}_{2}-\mathbf{r}^{0}_{1}$. Its
electrostatic field $E_{2u}$, given by 2.10 or 2.12, allows to calculate the
average angular beam–beam kick $\Delta
u^{\prime}=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 one-turn phase advance by $2\pi Q=\phi$, introduce
the constant $\Delta=\beta_{1}\Delta u^{\prime}=\beta_{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

$\displaystyle z_{n+1}$ $\displaystyle=(z_{n}-i\Delta)e^{i\phi}$
$\displaystyle=z_{1}e^{in\phi}-i\Delta
e^{i\phi}\left(1+e^{i\phi}+e^{2i\phi}+\ldots+e^{i(n-1)\phi}\right)$
$\displaystyle=z_{1}e^{in\phi}-i\Delta
e^{i\phi}\frac{e^{in\phi}-1}{e^{i\phi}-1}.$ (3.32)

This can be rewritten as

$z_{n+1}-\frac{i\Delta e^{i\phi}}{e^{i\phi}-1}=\left(z_{1}-\frac{i\Delta
e^{i\phi}}{e^{i\phi}-1}\right)e^{in\phi}.$ (3.33)

A comparison with the equation $z_{n+1}=z_{1}e^{in\phi}$ without the beam–beam
interaction shows that the constant kick only shifts the center of rotation
from the origin to the point

$z_{0}=\frac{i\Delta
e^{i\phi}}{e^{i\phi}-1}=\frac{e^{i\phi/2}}{2\sin\frac{\phi}{2}}=\frac{\Delta}{2\tan\frac{\phi}{2}}+i\frac{\Delta}{2}.$
(3.34)

Its real part determines the shift of the beam orbit, while the scaled
imaginary part $-Im(z_{0})/\beta_{1}=-\Delta u^{\prime}/2$ gives the angular
shift just before the beam–beam interaction. After receiving the kick $+\Delta
u^{\prime}$, it flips from $-\Delta u^{\prime}/2$ to $+\Delta u^{\prime}/2$,
while the next accelerator turn changes it back to $-\Delta u^{\prime}/2$.

The orbit shift expressed in the bunch widths should be invariant and should
depend only on $\tilde{\epsilon}_{1}/Z_{1}Z_{2}N_{2}$. Indeed, since
$\Delta=\beta_{1}\Delta
u^{\prime}=\left(\sigma_{1}^{2}p_{1}/\tilde{\epsilon}_{1}m_{1}\right)\cdot\left(eZ_{1}E_{2u}/p_{1}c\right)$,
it can be written as

$\frac{Re(z_{0})}{\sigma_{1}}=\frac{e(E_{2u}\sigma_{1})}{Z_{2}N_{2}}\cdot\frac{1}{2m_{1}c\tan\frac{\phi}{2}}\cdot\frac{Z_{1}Z_{2}N_{2}}{\tilde{\epsilon}_{1}}.$
(3.35)

Here, $E_{2u}$ has $\propto Z_{2}N_{2}/R$ dependence, so
$e(E_{2u}\sigma_{1})/Z_{2}N_{2}$ is invariant. For example, for the round
bunches with $\sigma_{1}=\sigma_{2}=\sigma$, $\Sigma_{x,y}=\sqrt{2}\sigma$
separated in $x$ by $\Delta x$, substituting $E_{2u}$ from 2.10 or 2.11 gives

$\frac{Re(z_{0})}{\sigma}=\frac{1-e^{-\Delta x^{2}/4\sigma^{2}}}{\Delta
x/\sigma}\cdot\frac{\alpha\hslash}{m_{1}\tan\frac{\phi}{2}}\cdot\frac{Z_{1}Z_{2}N_{2}}{\tilde{\epsilon}_{1}}.$
(3.36)

The maximal value of the first term depending only on $\Delta x/\sigma$ is
$0.31908\ldots\ $. Therefore, for protons one has

$\frac{Re(z_{0})}{\sigma}\leq\frac{N_{2}}{(2.0420\ldots\times
10^{12})\cdot[\tilde{\epsilon}_{1}/\mu\mathrm{m}]\cdot\tan\frac{\phi}{2}},$
(3.37)

which is less than 1% for the typical LHC values $\tilde{\epsilon}_{1}\approx
3-4\ \mu$m, $N_{2}\approx(7-9)\cdot 10^{10}$ and $\phi_{x}/2=0.31\pi$,
$\phi_{y}/2=0.32\pi$.

In case of several interaction points with the phase advances $2\pi
Q_{i}=\phi_{i}$ with respect to the first point where by definition
$\phi_{1}=0$, and the constant kicks $\Delta_{i}$, $i=1,2,\ldots,N$, the first
equation in section 3.8 modifies to

$z_{n+1}=z_{n}e^{i\phi}-i\Delta_{1}e^{i\phi}-i\Delta_{2}e^{i(\phi-\phi_{2})}-\ldots-i\Delta_{N}e^{i(\phi-\phi_{N})},$
(3.38)

which becomes equivalent to section 3.8 after the substitution

$\sum_{i=1}^{N}\Delta_{i}e^{-i\phi_{i}}=\Delta.$ (3.39)

So, we have again a circular trajectory with the center shifted to

$\displaystyle z_{0}$
$\displaystyle=\frac{i\left(\sum_{i=1}^{N}\Delta_{i}e^{-i\phi_{i}}\right)e^{i\phi}}{e^{i\phi}-1}=\frac{\sum_{i=1}^{N}\Delta_{i}e^{i[\phi/2-\phi_{i}]}}{2\sin\frac{\phi}{2}}$
$\displaystyle=\frac{\sum_{i=1}^{N}\Delta_{i}\cos(\frac{\phi}{2}-\phi_{i})}{2\sin\frac{\phi}{2}}+i\frac{\sum_{i=1}^{N}\Delta_{i}\sin(\frac{\phi}{2}-\phi_{i})}{2\sin\frac{\phi}{2}},$
(3.40)

whose real and opposite imaginary parts define the spatial and
$\beta_{1}$-scaled angular shifts, respectively. This formula is very well
known in the accelerator physics, which uses different tools for its
derivation Chao:1490001 .

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 $\sum_{i=1}^{N}w_{i}u_{i}$ over the
simulated particle coordinates $u_{i}$, averaged over accelerator turns
similarly to section 3.2. For example, with $(N_{turn}^{no\
bb},N_{turn}^{adiab},N_{turn}^{bb})=(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]\ \mu$m has the standard deviation $\sigma\lesssim 2$
nm. The maximal predicted orbit shift is $0.28\ \mu$m. The simulated shifts
without the beam–beam force are compatible with zero with $\sigma\lesssim 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
Figure 3.1. One particle trajectory corresponds to two circles in
$z_{x,y}$-planes. The orbit shift is then approximately determined by the
beam–beam force $\mathbf{F}^{tr}$ averaged over the trajectory like shown in
Figure 3.3, with the $x$-$y$ density projection 3.23. The circle from Figure
3.1 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_{y}$-plane. So,
different $z_{x,y}$-circle pairs have different $\mathbf{F}^{tr}$ and shift
differently. However, the ratio $r_{Re/Im}$ of the real and imaginary parts in
equations 3.34 and 3.8 remains constant regardless of $\mathbf{F}^{tr}$.
Therefore, the shifted centers $z_{0}=x_{0}+iy_{0}$ lie on the line
$r_{Re/Im}=x_{0}/y_{0}$ and the spatial and angular shifts are proportional.
Averaging over the bunch should preserve this proportionality. Therefore, if
$\Delta$, $\Delta_{i}$ denote exactly known average bunch angular kicks in
3.34 and 3.8, 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 $\Delta
u^{\prime}$ and $\Delta$, the values of $z_{0}$, ie. the angular and the orbit
shifts calculated from 3.34 and 3.8, 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, Figure 3.6 shows the $x$-orbit shifts in the same simulation as
in Figure 3.5. 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 section 3.8 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 $\sigma=0.4$ nm.

Figure 3.6: Upper row: $x$-coordinate averaged over the first bunch at four
experiments, lower row: its deviation from section 3.8 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 Figure 3.5.

In the old beam–beam model used at LHC in 2012–2019, the field
$\mathbf{E}_{2}$ in the momentum kick
$\Delta\mathbf{p}_{1}=q_{1}\mathbf{E}_{2}/c$ was 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
$\partial\Delta\mathbf{p}_{1}/\partial x$,
$\partial\Delta\mathbf{p}_{1}/\partial 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
$\partial\Delta p_{1,y}/\partial x$, $\partial\Delta p_{1,x}/\partial y$
vanish, and the $x$, $y$-slopes were taken as $\partial\Delta p_{1,x}/\partial
x$, $\partial\Delta p_{1,y}/\partial 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=\Delta x$, one has

$\displaystyle\left.\frac{1}{k_{1}}\frac{\partial\Delta p_{1,x}}{\partial
x}\right|_{x=\Delta x,y=0}$ $\displaystyle=-\frac{1}{\Delta
x^{2}}+\left(\frac{1}{\sigma_{2}^{2}}+\frac{1}{\Delta x^{2}}\right)e^{-\Delta
x^{2}/2\sigma_{2}^{2}},$
$\displaystyle\left.\frac{1}{k_{1}}\frac{\partial\Delta p_{1,y}}{\partial
y}\right|_{x=\Delta x,y=0}$ $\displaystyle=\frac{1-e^{-\Delta
x^{2}/2\sigma_{2}^{2}}}{\Delta x^{2}},$ (3.41)

where $k_{1}=2\alpha\hslash Z_{1}Z_{2}N_{2}$. The full kick at $\Delta
x$-separation was approximated as the following linear function of $x$ and
$y$:

$\displaystyle\Delta p_{1,x}$ $\displaystyle\approx
x\cdot\left.\frac{\partial\Delta p_{1,x}}{\partial x}\right|_{x=\Delta
x,y=0}+k_{1}\frac{1-e^{-\Delta x^{2}/2(\sigma_{1}^{2}+\sigma_{2}^{2})}}{\Delta
x}$ $\displaystyle\Delta p_{1,y}$ $\displaystyle\approx
y\cdot\left.\frac{\partial\Delta p_{1,y}}{\partial y}\right|_{x=\Delta
x,y=0}.$ (3.42)

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^{-\Delta
x^{2}/2\sigma_{2}^{2}})/\Delta x$. It contains $\sigma_{1}^{2}+\sigma_{2}^{2}$
instead of $\sigma_{2}$ in the exponent. Therefore, section 3.8 is not a
linear expansion of $\Delta 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\cdot\left.\partial\Delta
p_{1,u}/\partial u\right|_{x=\Delta x,y=0}$ represent the quadrupole magnets
usually used to focus or defocus the beams Chao:1490001 . They modify the
Gaussian widths according to the formula

$\sigma_{u}^{\prime}=\sigma_{u}\sqrt{1+\frac{\beta_{1}}{2\tan(2\pi
Q_{u})}\left.\frac{\partial\Delta p_{1,u}}{\partial u}\right|_{x=\Delta
x,y=0}}\ \ .$ (3.43)

Another simplification of section 3.8 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 madx12 , 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
Figure 3.7 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.

Figure 3.7: $(R-1)\times 10^{3}$ correction for the bunches with the
parameters from Table 1 colliding at one interaction point versus their
separation in the $x$-scan expressed in the bunch widths. The solid points are
the results of the B*B simulation when in the recurrence relation 3.10 the
exact formula (E), dipole constant dependent on the separation (D) or their
difference (E-D) are used. The open circles show the tabulated values from the
old simulation with the quadrupole kick (Q) and their sum with the values from
the dipole analytic formula (D+Q). The curves are analytic predictions for the
dipole (D), quadrupole (Q) kicks and their sum (D+Q).

To compare these known results with B*B, its exact kick formula was replaced
by section 3.8, by only its dipole or quadrupole parts. The simulated values
$R-1$ are shown in Figure 3.7 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 Figure 3.2. 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 $\sigma$ eg. 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

$\mu_{sp}\propto\frac{e^{-\Delta
x^{2}/2(\sigma_{1}^{2}+\sigma_{2}^{2})}}{2\pi\sigma_{1}\sigma_{2}}$ (3.44)

modulated by $R^{2}$ from Figure 3.7 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\neq 0.32$. The resulting biases of van der Meer
cross-section $\sigma^{\prime}/\sigma-1$ are shown in Figure 3.8 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 emittance $\tilde{\epsilon}/N$ assuming identical bunches
with $N=N_{1}=N_{2}$ and $\beta_{1}=\beta_{2}$. Therefore, this variable is
chosen as the horizontal axis. As shown by the dashed line, the difference
between 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.

Figure 3.8: The beam–beam cross-section bias $(\sigma^{\prime}/\sigma-1)\times
10^{3}$ determined from 1.6 versus $\tilde{\epsilon}/N\cdot 10^{11}$, where
$\tilde{\epsilon}$ is the normalized emittance
$\tilde{\epsilon}=\epsilon/(p/m_{p})$, identical for both bunches, $m_{p}$ is
the proton mass and $N=N_{1}=N_{2}$ is the number of protons. The bias is
calculated for the single interaction point using B*B, analytic dipole (D),
quadrupole (Q) approximations or their sum (D+Q). The tunes are $Q_{x}=0.31$,
$Q_{y}=0.32$ and it is assumed that $\beta_{x}=\beta_{y}=\beta$. Under these
conditions, $\sigma^{\prime}/\sigma$ depends only on the combination
$\tilde{\epsilon}/N$ but not on $\epsilon$, $p$, $N$ or $\beta$ individually.
The vertical dashed line shows the value $\tilde{\epsilon}/N\cdot 10^{11}$ for
the bunch parameters from Table 1.

##  4 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, eg. 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 Section 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