Modeling of nonlinear effects due to head-on beam-beam interactions

The beam-beam interaction is one of the most severe limitations on the performance of circular colliders, as it is an unavoidable strong nonlinear effect. As one aspires for greater luminosity in future colliders, one will simultaneously achieve stronger beam-beam interactions. We study the limitations caused by strong incoherent head-on beam-beam interactions, using a new code (CABIN) that calculates on a graphics processing unit (GPU), allowing for a detailed description of the long-term particle trajectories in 6D phase space. The evolution of the beam emittance and beam intensity has been monitored to study the impact quantitatively, while frequency map analysis has been performed to understand the impact qualitatively. Results from CABIN have shown good quantitative agreement with dedicated experiments in the Large Hadron Collider (LHC). For large beam-beam tune shifts, alternatives to the LHC tunes have been found to improve the beam quality. Schemes devised to cancel beam-beam driven resonances, by use of specific intermediate phase advances between the interaction points, work very well with zero crossing angle, and the accuracy required is achievable. Due to lack of symmetry, these schemes have an almost negligible impact with a significant crossing angle. The hourglass effect has been found to reduce the detrimental effects caused by the chromaticity and vice versa. The optimal level of the hourglass effect has been achieved when β 1⁄4 1.5σs. The ultimate parameters of the Future Circular Hadron Collider (FCC-hh) seem within reach, in absence of residual odd resonances.


I. INTRODUCTION
Concurrently with the hard collisions at an interaction point in a circular collider, multiple small angle deflections occur.This phenomenon, known as the beam-beam interaction, is caused by the electromagnetic fields from the opposing beam [1].The interaction force is strongly nonlinear, resulting in the possibility of betatron resonances [2].These resonances can cause a strong diffusion of the particles, or even make the beam unstable and lost within a short amount of time.The beam-beam interactions cause one of the most important limits on the performance of circular colliders, and they will possibly become even more important with the envisaged upgrades in the future, including the high-luminosity upgrade of the Large Hadron Collider (HL-LHC) and the larger Future Circular Hadron Collider (FCC-hh) [3][4][5].
In the FCC-hh the emittance will decrease during the fill, causing gradually stronger head-on beam-beam interactions [6].The usage of crab cavities is envisaged in order to mitigate the luminosity reduction due to the crossing angle, while keeping the impact of long-range beam-beam interactions under control [7].This paper therefore studies the impact of strong head-on beam-beam interactions with a crossing angle.
Although some effects of the beam-beam interaction can be derived analytically, numerical tools are necessary to assess the complete impact of the beam-beam interaction on the motion of incoherent particles, in particular their long-term behavior.The use of simulations has been an important tool in the understanding of the effects that cause instabilities [8][9][10].The design of new machines and optimization of existing ones requires a detailed knowledge of all mechanisms involved in the deterioration of the beam quality in the presence of strong beam-beam interactions.The weak-strong approach assumes that the bunch generating the electromagnetic fields stays unaltered [11].
In previous studies, it has been applied to study singleparticle stability, and the impact of the strength of the beambeam interaction [12,13].While already advanced, several mechanisms and mitigation schemes were not detailed, mainly due to the heavy computing load compared to the available resources.Here we propose an efficient numerical scheme to address these effects in a reliable way using modern computational resources.
In this paper, the interaction is studied by the weakstrong approach using a newly developed code (CABIN).A full parametric study of relevant parameters has been conducted.The discovered trends will be explained by the underlying mechanisms through frequency analysis [14], and compared to relevant previous studies [10,13,15].A realistic limit on the beam-beam interaction will be presented, relevant for the FCC-hh study.
Theoretical aspects of the beam-beam interaction are detailed in Sec.II.The numerical implementation of the physics model as well as the analysis tools are explained in Sec.III.Numerical results are presented in Sec.IV, along with a discussion of said results.Results from the LHC will be compared to results acquired with CABIN in Sec.V.The paper is concluded in Sec.VI.

II. THEORY
The beams in colliders as the LHC and the FCC-hh are crossing at multiple interaction points (IPs) with a nonzero full crossing angle θ xing , as illustrated in Fig. 1.Head-on (HO) electromagnetic interactions occur between the two beams at the IPs.In this paper, HO denotes interactions including a small crossing angle and a transverse separation.In an interaction region (IR) around each IP, each bunch experiences multiple additional parasitic long-range interactions with bunches from the opposing beam.The crossing angle is nonzero to avoid production of luminosity outside of the detector center and to reduce the effect of these long-range interactions on the beam quality.
The suggested use of crab cavities may significantly alter the interactions in the IR [7].Crab cavities work by tilting each bunch so that the two HO bunches in Fig. 1 overlap better, without reducing the crossing angle.It is proposed as a method to increase the event frequency without increasing the strongly distorting effects of the long-range interactions.In principle, the long-range interactions can be made arbitrarily small by increasing the crossing angle, limited eventually by the physical aperture of the final focusing magnets.In a high energy collider profiting from a significant radiation damping such as the FCC-hh, the effect of long-range interactions tends to decrease with the shrinking emittance, as opposed to the HO interactions which become stronger.This is assuming a fixed crossing angle.In such conditions it is clear that the HO interactions will become the main limitation rather than the long-range interactions.
The goal of a collider is to produce as many events as the physicists can analyze.The luminosity, L, is defined as the ratio of number of events detected per time per cross section.For two beams of equal size, σ i;1 ¼ σ i;2 , the luminosity is [16] where the subscript q denotes either transverse coordinate, x or y.H is a correction factor for the hourglass effect [17], being 1 for small σ s =β Ã q , and decreasing towards zero as this ratio increases.β Ã q is the transverse β function at the IP.S is a correction factor for the crossing angle, being 1 for zero crossing angle, and decreasing towards zero for increasing angles.N i is the number of particles per bunch in beam i. f rev is the revolution frequency.n b is the number of bunches in each beam.σ q is the rms transverse beam size, and σ s is the rms bunch length.Both an increase of beam emittance and a loss of particles will reduce the event frequency, assuming that the crossing angle remains fixed.
The dependence of luminosity on the crossing angle is given by Eq. (2).From this equation, a more relevant value is the Piwinski angle, ϕ PIW;q , which is defined as which is equal to the second term in the denominator in Eq. ( 2) for small angles.The crossing angle is typically small, θ xing ∼ 300 μrad in the LHC.

A. Beam-beam interaction in 4D
The HO beam-beam interaction will be studied by considering only one bunch per beam in the weak-strong regime.Beam 1 is taken to be weak and beam 2 is strong.The strong bunch will, in good agreement with real bunches, be taken to be a 3D Gaussian charge distribution.The distribution of the weak beam is not considered at the moment, because the interaction modeled in the weakstrong model is independent of it.
For nonzero crossing angle, θ xing , the bunches travel on different closed orbits around the IP.Furthermore, the transverse beam sizes, σ q , are dependent on the distance to the IP due to the hourglass effect.These effects change the beam-beam interaction as the bunches pass each other.In such a scenario, a full 6D treatment of the interaction is necessary.For θ xing ¼ 0 and σ s =β Ã q ≪ 1, the relation between the opposing bunches changes negligibly close to the IP, and the beam-beam interaction can be derived only dependent on the transverse dimensions.
The kick can be found by calculating the fields from the charge distribution, performing a Lorentz transformation of the fields, and integrating over the longitudinal direction [18].Other approaches also exist [1].Assuming round beams colliding in the 4D regime, one reaches the incoherent beam-beam kick where N is the number of particles in the strong beam, r 0 is the classical proton radius and γ is the Lorentz factor.At small radii, the kick is approximately linear.At higher radii, the force is strongly nonlinear.
If the bunches are flat, σ x ≠ σ y , the exponential function is exchanged with the Faddeeva function, sometimes referred to as the complex error function [1].

B. Beam-beam tune shift
The beam-beam parameter ξ q is a measure of the strength of the beam-beam interaction [16].It is equal to the beambeam tune shift for particles of zero transverse amplitude, jrj ¼ 0. For flat beams, the beam-beam parameter is different for the two transverse planes: The beam-beam parameter and beam-beam tune shift are negative for beams of the same charge.This paper refers to the negative of both values after this section.The combined tune shift from the machine tunes ðQ x0 ; Q y0 Þ for particles oscillating with different amplitudes, called the tune footprint, is displayed in Fig. 2. The maximum tune shift in either plane is equivalent to ξ q;Tot , neglecting the effect of any resonances, and occurs for particles of zero transverse amplitude.The tune shift then decreases for particles oscillating at larger amplitudes [19].

C. Resonances and resonance canceling
Nonlinearities in the machine, due to field imperfections of the magnets, misalignments between elements, vibrations and the beam-beam interaction, cause resonances.A general expression for transverse betatron resonances can be written as [2] where Q q are the transverse tunes.In this paper, Q q refers to the fractional tune.The tunes that fulfill this relation make lines in tune space of resonance order The motion can also be coupled between the longitudinal and transverse planes, named synchrobetatron motion.The coupling enables additional resonances, where Q s is the synchrotron tune.The additional resonances appear as sidebands to the betatron resonances.This makes it increasingly difficult to avoid all resonances in tune space.As it is impossible to avoid all of them, it is of interest to understand which of the resonances deteriorate the beam the most.
There are multiple approaches to study the nonlinearity of the beam-beam interaction quantitatively.The results presented here are based on Lie transfer maps of a 2D transverse phase space [1,10,20].We first consider a single IP in the ring.The potential of the beam-beam force can be expanded in a Fourier series.Because the beam-beam force is an odd function of the transverse position, only the even coefficients of its potential are nonzero.This is not true if there is a transverse separation or a nonzero crossing angle between the beams, breaking the symmetry.Furthermore, the coefficient strength increases with the action J, and it decreases with increasing order n.Put differently, beambeam resonances of higher order n have smaller resonance coefficients, and all nonzero resonance coefficients are larger at larger transverse amplitudes.
Consider next a similar setting with two IPs instead of one, splitting the lattice in two separate parts.In the LHC as in the FCC-hh, the main experiments are diametrically FIG. 2. Tune footprint caused by the beam-beam tune shift in both transverse planes simultaneously.The lines correspond to x=σ x ∈ f0; 1; 2; 3; 4; 5; 6g and y=σ y ∈ ½0; 6, or vice versa.The markers on the corners of the footprint correspond to the normalized amplitude in horizontal and vertical phase space.The red cross marks the working point.
positioned.The beam performs betatron motion, with phase advances μ 1 and μ 2 in the two parts of the lattice.In total μ 1 þ μ 2 ¼ μ in the entire lattice combined.By carefully adjusting the intermediate phase advance between the two IPs, μ 1 , it is possible to cancel resonances of order n.This can be achieved when [10] The resonance coefficients are still zero for odd n.Another suggestion for a phase advance that improves the beam quality is

Resonances of order
splitting the phase advance equally between the two separate parts of the lattice.This phase advance was found from a symmetry perspective [10], and confirmed numerically through trial to give the best performance [13].The suggested intermediate phase advances in Eqs. ( 9) and ( 10) will be tested in this paper.

III. NUMERICAL MODEL
CABIN (Cuda-accelerated beam-beam interaction) is a code developed to track particles of a single bunch through two interaction points, separated by two individual stretches of the magnetic lattice [21].The code is implemented using PYCUDA, exploiting the massively parallel architecture of a graphics processing unit (GPU).A comparison of the code on a Central Processing Unit (CPU) vs a GPU, timing of the different beam-beam implementations, and other requirements are detailed in Appendices B and C of [18].The GPU implementation has sped up the code by 3 orders of magnitude, dependent on the chosen hardware.As all particles are considered incoherently, the GPU is exploited maximally.Despite the high performance of the GPU, the numerical modeling of all effects impacting the beam quality, such as field errors of all magnets around the lattice, external noise sources, Coulomb scattering within the beam and with the opposing beam and quantum excitations, remains out of reach.The approximated physics models are described in the following.
The relevant machine and beam parameters used in the model are given in Table I, alongside the equivalent parameters in the LHC (design) [3], LHC (2017) [22], HL-LHC [23] and FCC-hh for bunches at 25 ns spacing [4,5], which will be used in the following unless stated otherwise.When performing beam-beam parameter scans, the intensity of the strong bunch is adjusted to obtain the desired beam-beam parameter using Eq.(5).

A. Lattice
A set of normalized coordinates with reduced units has been applied in CABIN.A hat refers to normalized values, as x in comparison to x.The phase space coordinates are given by The maximum beam-beam tune shift for the FCC-hh is not achieved with the initial parameters but rather once the emittance has decreased due to synchrotron radiation [5,6].
where σ is the initial rms beam size and momentum variation of the weak beam.The transverse actions are where ϵ q;0 is the transverse beam emittance at t ¼ 0. The beam emittance is calculated as ϵ q ¼ 2 • hJ q i.In the normalized coordinates, ðx i ; pi Þ, the transfer due to the lattice from one IP to the next is modeled as a rotation, xi pi where μ 1i is the intermediate phase advance in the three phase space planes through that part of the lattice.Going from the second IP to the first is done equivalently using μ 2i to complete the turn.The zero chromaticity phase advance is divided in any way desirable with the constraint μ 1q þ μ 2q ¼ μ q ¼ 2πQ q , where Q q is the transverse tune in that plane.
The tune of the particles depends on the momentum deviation.The linear chromaticity, Q 0 , changes the tune of the individual particles to Q q þ Q 0 • δ.The chromaticity can be treated in a symplectic manner, where also the longitudinal coordinates depend on the transverse oscillation [24].The symplectic approach is necessary for simulations of lepton machines.This study is focused solely on hadron machines that have a significantly smaller synchrotron tune, and the use of the simpler implementation is considered sufficient.In CABIN, the linear chromaticity is implemented and kept equal in the two transverse planes, Q 0 x ¼ Q 0 y , and it is divided evenly over the two sections as Since the synchrotron tune is small in hadron colliders, the smooth approximation is used in the longitudinal plane and therefore Q s is also divided evenly over the two sections.

B. Noise
Intrabeam scattering (IBS) is the small angle Coulomb scattering between particles in the same bunch.This effect is therefore stronger for denser beams [25].In a circular accelerator, the emission of synchrotron radiation in bending magnets continuously leads to an exponential damping of the emittances.The emission is however not smooth, on the quantum level, this radiation is emitted as discrete photons [26].The number of photons emitted per turn is high.Due to the large number of Coulomb scatterings as well as the large number of photons emitted per turn, the effect of both the IBS and quantum emission of the synchrotron radiation can be approximated as a Gaussian random noise on the momentum of the individual particles.
External sources of noise such as power converter ripple or ground vibration contribute to the random motion of the particles.Here we do not differentiate the internal and external sources of noise, we rather model the overall effect with a random variable.The damping due to synchrotron radiation occurs on larger time scales and is therefore not taken into account in the numerical model.
Consider one normally distributed kick per turn, κ q , on the normalized momentum pq , with rms amplitude Δ ¼ 1 × 10 −4 as discussed in the Appendix.Applying Eq. (A4), the diffusivity per turn is D ¼ 5 × 10 −9 .In result, the normalized beam emittance, εq , increases with 1 × 10 −8 =turn.This noise amplitude is compatible with the observations in the LHC and is dominated by external sources of noise [27].The amplitude is also comparable to the expected excitation amplitude due to internal sources of noise in the FCC-hh configuration when the large beam-beam parameters are reached [6].
Since for the FCC-hh the internal sources of noise are unavoidable, the noise will be enabled with the amplitude mentioned above in most of the following simulations.The noise amplitude is not adjusted to the machine and beam parameters, and in particular to the phase space density, such as to improve the comparability of the different effects of the beam-beam interaction in different configurations.

C. Beam-beam interaction
The 4D beam-beam interaction was justified in the limit of zero crossing angle, θ xing ¼ 0, and negligible hourglass effect, σ s =β Ã q ≪ 1.If either of these effects are present, the longitudinal dependence of the beam-beam kick is no longer negligible and a 6D kick is needed.
The SBM is conducted by first splitting up the beam in N S slices, whereupon the kicks from the individual slices in the boosted frame are calculated as in the 4D approach, with some additional modifications.
Because the beam-beam interaction is the most timeconsuming effect implemented, there is effectively a constraint on how many slices that can be used in the 6D simulations.It was found that the number of slices required to achieve a certain tolerance on the kicks depends on the configuration that is being simulated [18].For a small crossing angle, θ xing , and a negligible hourglass effect, σ s =β Ã q ≪ 1, the convergence is fast and a single slice may be sufficient.The required number of slices grows without bound as θ xing and σ s =β Ã q increase.For a certain tolerance, σ s =β Ã q ¼ 0.2 and ϕ PIW ¼ 0 requires N S ¼ 3.With a strong hourglass effect σ s =β Ã q ¼ 3.2 and ϕ PIW ¼ 0, N S ¼ 37 is required for the same tolerance.With σ s =β Ã q ¼ 0.2 and a large crossing angle ϕ PIW ¼ 2, N S ¼ 48 is required [18].The code finds for each simulation the required number of slices to achieve a certain tolerance.No simulations for round beams have been run with less than N S ¼ 15 slices.
Based on the mapping L in [28], one can already predict a few effects of the crossing angle.The effective transverse beam sizes in the boosted frame for a horizontal crossing are given by where the effective boosted horizontal beam size is divided by the factor S in Eq. ( 2), which explains the luminosity reduction for a nonzero crossing angle.
The maximum beam-beam tune shift for head-on collisions was found to be equal to the beam-beam parameter ξ q in Eq. (5).By insertion of the effective head-on beam sizes in Eq. ( 17) into that expression, one gets for a horizontal crossing where the subscript h signifies the horizontal crossing angle.Here η q;h=v is introduced as the effective beam-beam parameter, compared to the nominal jξ q j.A common crossing scheme in an accelerator is to have two IPs of equal strength jξ q j, with one horizontal and one vertical crossing.For round beams, this leads to a total tune shift where ξ Tot ¼ 2ξ is the sum of the nominal beam-beam parameters in the two independent IPs.A previous study has followed a different approach and found the impact of both the crossing angle and the hourglass effect on the tune shift [29].

D. Macroparticle distributions
The number of particles per bunch, N, in circular colliders like the LHC, is in the order of 10 11 .To reduce computation time and storage requirements, CABIN models the weak beam by a smaller number of macroparticles, N mp , of which each is tracked in six dimensions, at double precision.We require that the N mp macroparticles manage to represent well a Gaussian distribution, at low transverse amplitudes in the core to measure emittance growth, and at high transverse amplitudes in the tail to measure losses.It has been required that the bunch can represent the distribution up to 6σ, because the LHC collimation, which is usually set at 6σ, is set as a basis [3].The σ's in this section refer to the weak beam.
Multiple initial conditions (ICs) for distributing the macroparticles have been attempted.The distributions have been compared on how well they fill 6D phase space evenly through histograms of each coordinate separately, and the convergence of measurements of emittance growth rate and beam loss rate with N mp .The chi-squared distribution is defined as where N j are independently normally distributed variables, as the components of ðx; px ; ŷ; py ; ŝ; δÞ.The normalized 6D radius of a macroparticle will be referred to as How well the sums of the square of the coordinates of the tested ICs fit to the χ 2 6 distribution is visualized in Fig. 3.The Gaussian distribution is most equal to how the particles are distributed in a true bunch.Each macroparticle is weighted equally to represent N=N mp real particles.The core is modeled well, but the tail is not.
In the hollow Gaussian distribution, one half of the macroparticles have been distributed normally below 4σ in the 4D transverse phase space, and the second half have been distributed normally, requiring that they are above 4σ in the 4D transverse phase space [30].This has been done to achieve a better representation of the transverse tail.The longitudinal coordinates have been distributed as in the regular Gaussian distribution.The N reg macroparticles in each region separately are weighted equally, but the weighting is different between the different regions.The weight of all macroparticles in a given region depends on the cumulative density function (cdf) for the χ 2 4 distribution as where the limits on the 4D normalized radius are as mentioned above, f0; 4; ∞g.The core is modeled well.The tail is still less accurately represented than required for correct modeling of the beam losses.
In the 6D grid, the macroparticles are distributed linearly in each coordinate, whereupon these are meshed together in a grid.If the sixth root of N mp is not an integer, this value is rounded up to the closest integer, and the N mp macroparticles of lowest normalized 6D amplitude are kept.Macroparticle j is weighted individually based on its coordinates as whereupon the weights are normalized to give a total weight of N.This distribution does manage to represent particles in the tail, but it is far from smooth, causing multiple problems [18].The grid concept works well in lower dimensions, and a 2D grid in ðx; yÞ space will be used to perform frequency analysis.The uniform distribution is made by uniformly distributing each coordinate separately to a maximum amplitude of 6σ, and demanding that χ 6 ≤ 6.The individual macroparticles are then weighted by Eq. ( 21), making each macroparticle in the core represent more real particles than the macroparticles in the tail.The tail is much better modeled than with the previous distributions.The core is however not modeled accurately by this distribution, emphasized by the linear-scaled plot, which results in inaccurate modeling of the emittance growth.
The regionally uniform distribution is a modification of the uniform distribution, attempting to improve the modeling of the core of the bunch.The number of macroparticles is divided into three approximately equally large groups.The macroparticles in each group are distributed uniformly, requiring that they have a 6D radius within a given interval.The intervals are limited by the radii χ 6 ∈ f0; 2; 4; 6g.The number and size of the intervals were obtained by measuring the convergence of the observables of the code [18].Each macroparticle within a group is weighted as in Eq. ( 21).In addition, the macroparticles in each group are weighted proportional to the volume of the region, V reg , and to the inverse of the number of macroparticles in the region, N reg , making the weight of macroparticle j proportional to whereupon the weights are normalized to give a total weight of N. The good representation of the tail at high amplitudes is preserved from the uniform distribution.The core is much better modeled.The regionally uniform IC with N mp ¼ 1 × 10 5 macroparticles will be applied to study the evolution of the beam distribution, unless stated otherwise.

E. Beam quality simulations
This section provides information on how the beam quality is monitored, based on the modeled physics and initial conditions presented earlier.All values necessary to describe the individual particles are stored in snapshots every T mid ¼ 2 × 10 4 turns.This interval can be changed.The stored values are the six dynamical degrees of freedom (d.o.f.), whether a particle is lost or not, and the horizontal, vertical and longitudinal actions.The action of each particle is averaged over the next T ϵ ¼ 2048 turns, to avoid the angle dependent fluctuations.
During a normal physics fill in a high-energy circular hadron collider, the beam can circulate for hours.In the LHC, a fill time of 20 h equals T Tot ¼ 8 × 10 8 turns.In the FCC-hh, an estimated 3 h of run-time would equal T Tot ¼ 3 × 10 7 turns.Multiple reasons make tracking for that long using the weak-strong model inappropriate.The results presented later have been acquired from tracking particles in the weak beam through usually T Tot ¼ 2 × 106 turns, equivalent to 3 min in the LHC.

Emittance growth
We distinguish between two types of emittance growth mechanisms.Close to resonances, the effective Hamiltonian diverges, and subsequently the actions are strongly perturbed.This leads to an almost immediate and possibly large emittance growth.A different important measure is the longterm emittance growth rate, _ ϵ, after the initial readjustment.To monitor the long-term growth rate, a linear regression curve of the beam emittance is computed, starting at turn T sim ¼ 5 × 10 5 .The evolution of the emittance for the example configuration is displayed in Fig. 4. The initial emittance growth of a few percent is visible, and the regression curve for each transverse emittance growth rate is plotted on top of the evolution.The emittance growth rates in both planes are larger than _ ϵ 0 ¼ 1 × 10 −8 =turn, expected from the noise alone.The average transverse emittance growth rate will be referred to as The convergence of the emittance growth rate with N mp has been calculated for the Gaussian IC and the regionally uniform IC with noise only.The average and the standard deviation from _ ϵ 0 has been calculated based on 40 independent growth rates for each N mp .A 1σ trend line for the relative errors has been calculated.For N mp ¼ 1 × 10 5 particles, the Gaussian distribution has a 1σ relative error of 5%, while the regionally uniform distribution has a 1σ relative error of 15% [18].
The convergence study was repeated for a configuration at the LHC working point ðQ x ; Q y Þ ¼ ð0.31; 0.32Þ, also including a large beam-beam parameter ξ Tot ¼ 0.03 divided over two IPs, and a large chromaticity, Q 0 ¼ 15, in addition to the noise.The emittance growth rate is larger in this configuration.The relative errors of the emittance growth rates are for this configuration calculated relative to the average emittance growth rate calculated with N mp ¼ 1 × 10 6 particles.For N mp ¼ 1 × 10 5 particles, the Gaussian distribution has a 1σ relative error of approximately 16%, while the regionally uniform distribution has a 1σ relative error of 5%.Both ICs are reasonably accurate in calculating the emittance growth rate for both considered configurations.

Beam loss
One unavoidable source of particle losses is the luminosity burn-off, which CABIN does not take into account.The other main source of beam losses is particles that drift to large transverse amplitudes, whereupon they are extracted by the collimation system of the accelerator [31].Some particles are transported beyond the limit within a few couple of turns, following the same process as the initial emittance growth.After the initial losses, the longterm diffusion can also make particles drift out.These losses are partly unavoidable without any damping mechanisms, and keeping them low is a main objective in the design and operation of a collider.The long-term loss rate of particles is calculated based on the conserved intensity for each snapshot.The loss rates are calculated as the relative decrease from the initial intensity, by a linear regression starting at turn T ¼ 5 × 10 5 .
A particle in a real machine is lost at large amplitudes because it hits the collimation system.In the LHC, this is typically at a transverse radius of R M ¼ 6σ ϵ , where σ ϵ is the transverse sigma corresponding to the operational emittance of the beams.In CABIN, the transverse radius is checked once per turn, and if it is ever beyond a given limit R, it is considered lost forever.The particle is still tracked to study the emittance growth to larger amplitudes.The intensity below the machine limit, R M , and a beam limit, R B ¼ 5σ Σ , are displayed in Fig. 5.Here σ Σ is the beam size of the strong beam in the simulation.Many particles are lost beyond R B that appear to stop between the two limits.This is due to the fact that the strength of the beambeam interaction decreases fast for larger amplitudes.It has been assumed that particles being transported to R B are likely to be lost within a small number of turns due to other effects such as the lattice nonlinearities which, as opposed to head-on beam-beam interactions, become stronger at higher amplitudes.Hence, the results presented later refer to the loss rate from the beam (LR Beam ) beyond R B .The values are scaled to represent loss rates per hour in a ring with revolution frequency f rev ¼ 11.245 kHz, such as the LHC.If the beams are flat, the limits are taken proportional to the relevant σ, giving an elliptical collimation.In order to fully describe the diffusion of the particles from R B to the actual collimator position, a complete nonlinear model of the magnetic lattice should be implemented, but that is beyond the scope of this study.
The beam loss rate is mostly dependent on the particles in the tail of the bunch.The convergence of the beam loss rate with N mp has been calculated for the Gaussian IC and the regionally uniform IC, for a configuration only including a noise of rms amplitude Δ ¼ 1 × 10 −4 .The average loss rates have been calculated based on 20 independent values for each N mp .The standard deviations have been calculated relative to the average loss rate calculated for N mp ¼ 1 × 10 6 particles.For N mp ¼ 1 × 10 5 particles, the Gaussian distribution has a 1σ relative error of approximately 86%, while the regionally uniform distribution has a 1σ relative error of 8%.
The convergence of the beam loss rate was repeated for a configuration at the LHC working point ðQ x ; Q y Þ ¼ ð0.31; 0.32Þ, also including a large beam-beam parameter, ξ Tot ¼ 0.03, divided over two IPs, and a large chromaticity, Q 0 ¼ 15, in addition to the noise.The error was also for this configuration calculated relative to the average loss rate for N mp ¼ 1 × 10 6 particles.For N mp ¼ 1 × 10 5 particles, the Gaussian distribution has a 1σ relative error of 8%, while the regionally uniform distribution has a 1σ relative error of 5%.The Gaussian IC is strongly inaccurate for configurations that would give small loss rates.That is the main reason why the regionally uniform distribution is used for the beam quality simulations in this paper.

F. Frequency map analysis
The beam quality simulations introduced in Sec.III E enable CABIN to quantitatively study the evolution of measurable quantities.Frequency map analysis (FMA) is implemented in CABIN to study the underlying mechanisms driving the beam quality degradation.FMA is a method that can qualitatively visualize the detrimental effects of beambeam interactions, based on the evolution of the tunes of individual particles [14].
For a configuration only including the linear lattice for zero chromaticity, all particles have tunes ðQ x ; Q y Þ. Due to the nonlinear beam-beam interaction, the tunes of each particle j are shifted to ðQ xj ; Q yj Þ separately, as was detailed in Fig. 2. The action of a particle may increase close to a resonance, leading to a change of the tune.Hence, it is of interest to monitor the tunes of individual particles closely, both to know which resonances may be excited, and because the change in tune can indicate the severity of the diffusion close to a given resonance.One can calculate a tune diffusivity, D j , for each particle, based on how much the tunes change per turn where T is the turn number.The change per turn can be calculated as a linear regression based on the different tune measurements.The tune diffusivity is then used as a color scale for each particle, presented as a footprint in tune space, ðQ x ; Q y Þ, or in initial amplitude space, ðA x ; A y Þ.In the results that will be presented later, the diffusivity constant, D j in Eq. ( 24), is calculated for N mp ¼ 4 × 10 4 particles.These are initially distributed in a quadratic, linear IC in the first quadrant of the ðx; yÞ plane, with sðT ¼ 0Þ ¼ 1, unless stated otherwise.The FMA simulations will be run with zero noise.The tunes may be drifting, or they may oscillate around a mean, as would be expected e.g. for a nonzero chromaticity.The tune oscillation with chromaticity has a period of 1=Q s turns.By measuring each tune based on data from more than 1=Q s turns, the periodic tune variations due to the chromaticity will be reduced.By calculating the diffusivity based on multiple tunes, the calculation is less sensitive to periodic variations caused by either chromaticity or other sources.On the other hand, the tune calculations are slow, favoring less points.As a compromise, the tune diffusivity will be calculated as the slope of a linear regression of n Q ¼ 10 tunes spaced by T mid ¼ 1 × 10 4 turns.This combination has been found through a convergence study [18].
A FFT can only calculate tunes of resolution 1=T Q where T Q is the number of turns used for the calculation.For the FMAs presented later, the tunes are calculated based on T Q ¼ 4096 turns.The individual tunes are calculated using an interpolated FFT method implemented in the code HARPY [32].

IV. NUMERICAL RESULTS
The beam-beam interaction is calculated in the weakstrong approach, assuming that the beam distributions change negligibly during the simulation.Self-consistent models overcome this limitation, usually at the cost of a large increase of the computing load [8,9].Nevertheless, the simulation results showing only a weak distortion of the beam, corresponding to potential operational configurations, remain well approximated by the weak-strong model [33].

A. Search for working point
The total beam-beam tune shift is relatively small in the LHC today, and a working point close to the nominal ðQ x ; Q y Þ ¼ ð0.31; 0.32Þ preserves the beam quality sufficiently well.Future upgrades will have a larger beambeam tune shift, as this is necessary to achieve a higher luminosity.
To look for other working points that are possibly more suitable for larger beam-beam tune shifts, a working point scan is performed in ðQ x ; Q y Þ space.The scans are performed with a tune resolution of 0.0025.Because there is complete symmetry between the tracking of horizontal and vertical phase space in CABIN, only the upper half of the diagonal has been calculated.All simulations are run for a total beam-beam tune shift of ΔQ Tot ðϕ PIW Þ ¼ 0.03, combined with a large chromaticity, Q 0 ¼ 15, and a considerable Piwinski angle, ϕ PIW ¼ 0.9, corresponding to a full crossing angle of θ xing ¼ 300 μrad.This is considered similar to the worst-case scenario for the FCC-hh, based on the values presented in Table I.The crossing angle is chosen nonzero to also include odd resonances in order to remain in a pessimistic scenario.
The tune scan on the interval Q q ∈ ½0.255; 0.345, in both transverse tunes, is presented in Fig. 6.There is generally better preservation of the beam quality with a working point on or close to the coupling resonance, Q x ¼ Q y .While it is not excluded that such working points can be used, a proper demonstration is required, since for example coherent instabilities [34] and skew magnetic errors are not included in the numerical model.To be conservative, it is assumed in the following that a finite tune separation of jQ x − Q y j ≥ 0.01 is needed, as for the LHC.The working point with the lowest beam loss rate and emittance growth rate in this tune area, at least this far from the diagonal in tune space, is ðQ x ; Q y Þ ¼ ð0.315; 0.325Þ.
The tune scan on the interval Q q ∈ ½0.45; 0.505, in both transverse tunes, is presented in Fig. 7.The working point with best preservation of beam quality in this tune area, requiring again that jQ x − Q y j ≥ 0.01, is ðQ x ; Q y Þ ¼ ð0.475; 0.485Þ.Close to the half-integer resonance, it can be difficult to correct for dipole errors and nonideal optics.However, modern correction techniques suggest that such working points may be attainable [35].
The footprint of the bunch in ðQ x ; Q y Þ space is wider, in addition to longer, with a larger beam-beam parameter.It could therefore be argued that the difference between the tunes in the alternative working points should be larger FIG. 6. Scan in ðQ x ; Q y Þ space, for ΔQ Tot ðϕ PIW Þ ¼ 0.03, β Ã q ¼ 40 cm, ϕ PIW ¼ 0.9 and Q 0 ¼ 15.The black cross (×) marks the LHC working point, ðQ x ; Q y Þ ¼ ð0.31; 0.32Þ.The green diamond ( ) marks the suggested working point for this configuration, ðQ x ; Q y Þ ¼ ð0.315; 0.325Þ.

FIG. 7. Scan in ðQ
q ¼ 40 cm, ϕ PIW ¼ 0.9 and Q 0 ¼ 15.At the working points marked by blue squares with white crosses, the entire beam was lost before the long-term loss rate could be calculated.The green diamond ( ) marks a suggested working point in this tune area, ðQ x ; Q y Þ ¼ ð0.475; 0.485Þ.than jQ x − Q y j ¼ 0.01.To compare to the LHC working point, we keep the tune separation equal.

B. Crossing angle
This section studies the effect of varying the crossing angle for the main head-on interaction.Beam quality simulations for round beams have been run for each of the following studies, with Piwinski angles ϕ PIW ∈ f0; 0.1; 0.3; 0.6; 1; 1.5; 2g, and beam-beam parameters ξ Tot ∈ f0.01; 0.02; 0.03; 0.04; 0.05g, from which the effective maximum beam-beam tune shifts can be calculated by Eq. (19).
The results from a crossing angle scan for the machine tunes of ðQ x ; Q y Þ ¼ ð0.31; 0.32Þ, and zero chromaticity, are presented in Fig. 8(a).As the angle increases from ϕ PIW ¼ 0 to 0.6, the detrimental effects gradually become stronger, strongly dependent on ξ Tot .This behavior is qualitatively predictable, because the beam-beam interaction with a crossing angle drives odd resonances as well.For ξ Tot ¼ 0.03 and ϕ PIW ¼ 0, the particle trajectories are affected by the tenth order resonances and weakly by the 16th order, as seen in the FMA in Fig. 9(a).There is only weak overlap between different resonances.The 13th order resonance is active with ϕ PIW ¼ 0.1, as displayed in the FMA in Fig. 9(b), and it affects predominantly particles at large transverse amplitudes.This resonance may contribute to the large loss rates combined with the moderate emittance growth rates.This FMA was produced for a small angle, ϕ PIW ¼ 0.1, because the picture became distorted beyond clear understanding for larger angles.The additional detrimental effects caused by the crossing angle saturates around ϕ PIW ≈ 1.For larger angles, the beam quality preservation begins to improve as the crossing angle increases further.This improvement is due to the reduction of the effective beam-beam tune shift.
Chromaticity is often needed in modern colliders to mitigate collective instabilities.The crossing angle scan for the machine tunes of ðQ x ; Q y Þ ¼ ð0.31; 0.32Þ was therefore repeated with a large chromaticity, Q 0 ¼ 15, as is currently needed in the LHC [36], and the results are presented as functions of ΔQ Tot ðϕ PIW Þ in Fig. 8(b).The most evident effect is that the beam quality is worse preserved for small angles, ϕ PIW ≤ 0.1, and large beam-beam tune shifts, ΔQ Tot ðϕ PIW Þ ≥ 0.02, and that the emittance growth rate has increased substantially.The chromaticity has the effect that the tune of a particle oscillates in tune space, parallel to the diagonal, Q x ¼ Q y , as Q 0 is kept equal in both planes.Effectively, particles far from resonance lines become affected by resonances, visible in the FMA in Fig. 9(c).This oscillation is dependent on the energy variation, causing a mixing between the longitudinal and transverse planes, activating synchrobetatron resonances.The introduction of the longitudinal d.o.f. is able to increase the diffusion rate through e.g.Arnold diffusion [37].The chromaticity does not activate odd resonances, since it does not introduce a transverse asymmetry to the interaction.Therefore, there still is an increase in detrimental effects towards a saturation as the angle increases towards ϕ PIW ≈ 1.
The crossing angle scan with a nonzero chromaticity was repeated at the shifted working point, (0.315,0.325), and the results are presented as functions of ΔQ Tot ðϕ PIW Þ in Fig. 10(a).Both outputs support that this working point might be better than the LHC working point, except for the configurations with ΔQ Tot ðϕ PIW Þ ≤ 0.02 and ϕ PIW ≥ 0.3, which indeed are the most relevant to describe the configurations in the LHC today.With the shift of working point, other resonances are important to understand the dynamics.The FMA for ξ Tot ¼ 0.03, Q 0 ¼ 15 and ϕ PIW ¼ 0.3 is presented in Fig. 11(a).The odd and even resonances overlap for particles at large transverse amplitude.This explains the clear distinction between the cases of zero and nonzero crossing angle, based on the loss rate.The overlap is close to the working point, and will affect the footprint also for smaller ΔQ Tot ðϕ PIW Þ.This explains why this working point appears slightly worse than the LHC working point, for small ξ Tot and large ϕ PIW .
The results from a crossing angle scan for the machine tunes of ðQ x ; Q y Þ ¼ ð0.475; 0.485Þ, for nonzero chromaticity, Q 0 ¼ 15, are presented in Fig. 10(b) as functions of ΔQ Tot ðϕ PIW Þ.A strange effect is that, for ϕ PIW ¼ 0, the beam quality seems better preserved for a stronger beambeam interaction.This was identified as an artifact of the numerical model, due to the fact that the noise amplitude is input relatively to the unperturbed optics functions, which are actually significantly affected by the beam-beam interaction in the vicinity of the half-integer resonances [18].The extraordinary good behavior for ξ Tot ¼ 0.05 and zero crossing angle can be understood by the FMA in Fig. 11(b).There  are only a few resonance lines that affect the footprint, and there is no overlap between resonances.The maximum beam-beam tune shift is smaller than ξ Tot because of the half integer resonance.For a nonzero crossing angle, odd higher order resonances overlap with the even resonances already present with zero crossing angle.This working point gives by far the best long-term preservation of beam quality for the considered configurations.

C. Separation
By separating the beams at the IP in the plane transverse to the crossing angle, the luminosity can be reduced intentionally.The crossing scheme in these simulations is both with and without crossing angle HV (horizontal in the first IP, vertical in the other), and the separation scheme is thus VH.Because the beambeam interaction is weaker at larger amplitudes, one could naively assume that the beam quality would be preserved better with a nonzero separation.To not reduce the luminosity to zero, the relevant separations are of only a few σ Σ , which changes the nonlinear force completely.It is therefore necessary to do numerical simulations to find out if the separation has a beneficial impact on the beam quality.At larger separations, it is assumed that the impact of the beam-beam interaction will decrease.Beam quality simulations for round beams have been run for each of the following studies, with separations Δx ⊥ ∈ f0; 0.25; …; 2; 3; …; 10g • σ Σ , large chromaticity, Q 0 ¼ 15, and nominal beam-beam parameters ξ Tot ∈ f0.01; 0.02; 0.03; 0.04; 0.05g, at the LHC working point.
The results for a separation scan with zero crossing angle are presented as functions of the separation Δx ⊥ in Fig. 12(a).For a separation of 6σ Σ , the beam loss rate has converged to good values.For a separation of 4σ Σ , the emittance growth rate has converged to the value expected from noise alone.At lower amplitudes, the loss rate and emittance growth rate both oscillate with local peaks at separations of approximately 0.5σ Σ and 1.5σ Σ .This is in agreement with an earlier study [15].The worse behavior at small separations is due to the activation of odd resonances, as visualized by the FMA in Fig. 13.
The same scan was repeated with a nonzero Piwinski angle, ϕ PIW ¼ 1, and the results are presented as functions of the separation in Fig. 12(b).The loss rate and emittance growth rate converges again at 6σ Σ and 4σ Σ respectively.The worst beam quality preservation occurs on the interval Δx ⊥ ∈ ½0.25; 1.5 • σ Σ .The behavior in this interval is smoother with a nonzero crossing angle, and the increased deterioration is smaller.That is because the odd resonances are already activated by the nonzero Piwinski angle.
For a transverse separation of 5σ r , the particles of approximately zero transverse amplitude always experience the weak beam-beam force at 5σ r , independent of the betatron phase.The particles at large transverse amplitudes up to 5σ r will oscillate between 0 and 10σ r , experiencing at times the strongly nonlinear force in the core of the strong beam.This explains the convergence of emittance growth rate at a lower transverse separation than the convergence of the loss rate.

D. Intermediate phase advance
The focus in this section is the impact of the intermediate phase advance between the two IPs in the model, further explained in Sec.II C. A previous study has found that certain resonances can be canceled by setting μ 1 ¼ π=2 or μ 1 ¼ μ=2 [10].In this analytical model, a phase error of more than 0.001 could be sufficient to lose the cancellation effect.Since those constraints were not compatible with realistic optics tolerances, no conditions on the phase advance between IPs were considered for the LHC.However, the theoretical model does not allow for an estimation of the gain in terms of beam quality for a partial cancellation of the resonances, which could allow for more relaxed tolerances on the optics control.
Beam quality simulations have been run for configurations with different intermediate phase advances, Δμ 1 ∈ ½−0.03; 0.03 • 2π, close to the phase advances introduced above, nominal beam-beam parameters ξ Tot ∈ f0.03; 0.04; 0.05g, and large chromaticity, Q 0 ¼ 15.The difference, Δμ 1 , from the design phase advance, π=2 or μ=2, is kept equal in the horizontal and vertical plane.Because the goal is to study how much the beam quality preservation improves, and how accurately the intermediate phase advance must be set, only the largest beam-beam parameters have been simulated.Only the results at the LHC working point, (0.31,0.32), will be presented.Similar behaviors have been found at the alternative working points (0.315,0.325) and (0.475,0.485).
The results from a scan in intermediate phase advance close to μ=2 for zero crossing angle are presented as functions of μ 1 − μ=2 in Fig. 14(a).The best preservation of beam quality is achieved when μ 1 ¼ μ 2 .This is supported by the FMA in Fig. 15(a).The improvement can be understood from a symmetry argument.With equal phase advances, and equal modeling of the beambeam interaction in the two IPs, the model represents effectively a collider with a single IP at a working point of ðQ x ; Q y Þ ¼ ð0.155; 0.16Þ.The total beam-beam tune shift would then be half of what it is in the original configuration.The required accuracy in μ 1 seems dependent on the beam-beam parameter, larger ξ Tot requires a more accurate intermediate phase advance to achieve the full improvement.The beam quality simulations are less sensitive to the phase error than what the resonance coefficients are.
The scan in intermediate phase advance close to μ=2 has been repeated with a nonzero crossing angle, ϕ PIW ¼ 1.The results are presented as functions of μ 1 − μ=2 in Fig. 14(b).With one horizontal and one vertical crossing, the symmetry is broken and therefore the lattice is no longer equivalent to a shorter lattice with a single IP.The resonances are not canceled, and the improvement from the use of this intermediate phase advance is limited.This argumentation is supported by the FMA in Fig. 15(b).The results from a scan in intermediate phase advance close to π=2, for zero crossing angle, are presented as functions of μ 1 − π=2 in Fig. 16(a).The best preservation of beam quality is achieved when μ 1 ¼ 0.512 • π, not when μ 1 ¼ π=2.The optimum phase advance has been found to be dependent on the working point.The theory applied to achieve the resonance canceling condition, μ 1 ¼ π=2, was for a 2D transverse phase space, and assumed that the beam-beam interaction was weak.This assumption has possibly been broken.For resonances dependent on both transverse phase space planes, it seems that π=2 is not the optimal choice.This is supported by FMAs.As seen in Fig. 17(a), the cancellation of the resonances is not as good as with the symmetric phase advance.The FMA can however only give a qualitative understanding of the beam dynamics.The beam quality simulations did show an equally good improvement in preservation of beam quality, and similar required accuracy in maintaining μ 1 .
The scan in intermediate phase advance close to π=2 has been repeated with a nonzero crossing angle, ϕ PIW ¼ 1.The results are presented as functions of μ 1 − π=2 in Fig. 16(b).The improvement is modest at best.The choice μ 1 ¼ π=2 was only supposed to cancel resonances of order f2; 6; 10; 14; 18; …g, not the 7th, 13th or 16th.This is supported by the FMAs for a small crossing angle in Fig. 17(b).Even if there is an improvement compared to Fig. 9(b), there is no improvement according to the beam quality simulations.

E. Hourglass effect
The β function varies as a parabola close to the IPs, causing the transverse beam size to vary as well.This variation is weak for configurations where σ s =β Ã q ≪ 1. β Ã q is squeezed to small values in all modern colliders to increase the luminosity, making the hourglass effect important as a side effect.Beam quality simulations for round beams have been run for the following studies with β Ã q ∈ f2.5; 4; 8; 10; 12; 14; 20; 30; 40; 60; 80; 160; 400g cm, and ξ Tot ∈ f0.01; 0.02; 0.03; 0.04; 0.05g.The simulations are run with σ s ¼ 8 cm and zero crossing angle.
The results from a β Ã q scan at the LHC working point, ðQ x ; Q y Þ ¼ ð0.31; 0.32Þ, with zero chromaticity, are presented as functions of σ s =β Ã q in Fig. 18(a).A change occurs for σ s =β Ã q ≥ 1, for which the loss rates and emittance growth rates are increasing.The strong hourglass effect causes significant mixing with the longitudinal d.o.f., activating and driving the synchrobetatron resonances.This is supported by the FMA for β Ã q ¼ 2.5 cm in Fig. 19(a).The effect is similar to the one caused by chromaticity.
The β Ã q scan for the LHC working point, ðQ x ; Q y Þ ¼ ð0.31; 0.32Þ, has been redone with a large chromaticity, Q 0 ¼ 15, and the results are presented in Fig. 18(b).The  chromaticity activates synchrobetatron resonances visualized by the FMA for β Ã q ¼ 4 m and Q 0 ¼ 5 in Fig. 19(b).That is why the beam quality is worse preserved for nonzero chromaticity and σ s =β Ã q ≪ 1.This behavior was also expected for the more significant hourglass effect.Instead, there appears an optimum σ s =β Ã q of approximately 2=3, limited by the chromaticity from below and the hourglass effect from above.It seems that when the hourglass effect becomes significant, it reduces the strength of some resonances that strongly affect the dynamics in the presence of chromaticity, visualized by the FMA for β Ã q ¼ 12 cm and Q 0 ¼ 5 in Fig. 19(c).The hypothesis is supported by both the beam quality simulations and FMAs.
At the first alternative working point, ðQ x ; Q y Þ ¼ ð0.315; 0.325Þ, we observe the same behavior, with a wider optimal interval of β Ã q ∈ ½10; 30 cm.The beam quality is preserved better for all tested configurations at this working point.At the second alternative working point, ðQ x ; Q y Þ ¼ ð0.475; 0.485Þ, there is no noticeable impact of the chromaticity for zero crossing angle, as was found in Sec.IV B. The loss rate increases fast for σ s =β Ã q ≥ 0.4 [18].
The results from a chromaticity scan with negligible hourglass effect, σ s =β Ã q ¼ 0.02, and zero crossing angle, ϕ PIW ¼ 0, are presented as functions of Q 0 in Fig. 20(a).The scan was repeated with a significant hourglass effect, σ s =β Ã q ¼ 2=3, and the results are presented in Fig. 20(b).The most important difference is that the onset of detrimental effects begins at a larger chromaticity with the significant hourglass effect.For the largest beam-beam parameter, ξ Tot ¼ 0.05, there seems to be a cancellation of the detrimental effects caused by the hourglass effect at Q 0 ¼ 0, for a small nonzero chromaticity.This behavior produces an optimum chromaticity, being approximately Q 0 ¼ 3.75 for this configuration.
The chromaticity scans have been repeated with a significant crossing angle, ϕ PIW ¼ 1.The results from a scan with negligible hourglass effect, σ s =β Ã q ¼ 0.02, is presented as functions of Q 0 in Fig. 21(a).The scan with a significant hourglass effect, σ s =β Ã q ¼ 2=3, is presented in Fig. 21(b).Note that this Piwinski angle for a smaller value of β Ã q corresponds to a smaller actual crossing angle.For the largest beam-beam parameters, ξ Tot ≥ 0.04, the detrimental effects caused by the hourglass effect and crossing angle when Q 0 ¼ 0, is reduced when a small nonzero chromaticity is included.This behavior produces an optimum chromaticity, being approximately Q 0 ¼ 2.5 for this configuration.
The optimum balance between chromaticity and hourglass effect, observed both here and in Sec.IV E, is assumed to be linked to how the two effects alternate in affecting the incoherent particles.The hourglass effect is felt more strongly for a large longitudinal displacement, s, while the chromaticity affects the particles more strongly for a large energy deviation, δ.This alternation with the synchrotron motion is assumed to disturb the resonances from building up periodically over time.Exactly how this mechanism works remains to be understood.Unlike the smart choice of phase advance, the improvement from this effect has been found to remain after activation of odd resonances with a crossing angle.

G. Maximum beam-beam tune shift
In lepton colliders, a limit in the beam-beam tune shift is usually observed [1].Having studied individually the different beam quality deterioration mechanisms in a hadron collider, we now ask ourselves whether a realistic beam-beam limit exists also there.Beam quality simulations have been run for different configurations at the LHC working point, (0.31, 0.32), and the alternative working points, (0.315, 0.325) and (0.475, 0.485), in search of a maximum acceptable beam-beam tune shift.The limit will be taken where there appears to be a distinct threshold on the preservation of beam quality, and not linked to a hard limit on either output.Pessimistically, the chromaticity will be set to Q 0 ¼ 15 in these simulations, assumed necessary to prevent coherent instabilities.The beta function will be set to β Ã q ¼ 12 cm, as have been found optimal to reduce the problems caused by this level of chromaticity at the LHC working point.This β Ã q is similar to the design value for the HL-LHC.
The first configuration that will be considered has round beams colliding head-on with zero crossing angle.The intermediate phase advance is set to μ 1 ¼ μ=2.This is similar to the configuration considered in [13], and the results presented as functions of ΔQ Tot ðϕ PIW Þ in Fig. 22(a) are in good agreement.Indeed, in this configuration, the maximum acceptable beam-beam tune shift from the LHC working point is approximately ΔQ Tot ¼ 0.26.The limits are slightly larger for the alternative working points.This configuration is however quite unrealistic.First of all, it is not possible to keep exactly the smart intermediate phase advance μ 1 ¼ μ=2.The acceptable phase error has been found to decrease for increasing ξ Tot , and was approaching zero for a beam-beam parameter of ξ Tot ¼ 0.05.Second of all, it was assumed that the nonlinearities of the lattice, and in particular the corresponding resonance driving terms, could be reduced arbitrarily, and were consequently neglected in the numerical model.It is obvious that this cannot be achieved for low order resonances such as the fourth order that is strongly excited for example by the sextupoles necessary for chromatic corrections.This resonance would limit the beam-beam tune shift to at least 0.06 with the LHC working point.Similarly the absence of odd resonance in this configuration seems difficult to realize in practice.
The required accuracy in the intermediate phase advance is perhaps impossible to achieve for large beam-beam parameters, ξ Tot > 0.05.The scan for zero crossing angle has therefore been repeated with zero intermediate phase advance, and the results are presented as functions of ΔQ Tot ðϕ PIW Þ in Fig. 22(b).Without the improvement from the intermediate phase advance, the maximum acceptable beam-beam tune shift from the LHC working point is reduced significantly to approximately ΔQ Tot ¼ 0.043.Above this limit the loss rate increases rapidly.The limit for the first alternative working point is approximately 0.067.The limit is similar for the second alternative working point, although it is less abrupt.
In future colliders, the strength of long-range interactions might be weaker and the Piwinski angle might be canceled by use of crab cavities.Nevertheless, it is possible that odd resonances cannot be fully suppressed due to the imperfections of the crabbing system and the remaining longrange interactions.Therefore, the beam-beam parameter scan has been repeated for a configuration with a significant Piwinski angle, ϕ PIW ¼ 1, and zero intermediate phase advance.The results of this scan are presented in Fig. 22(c).The limit at the LHC working point is approximately ΔQ Tot ¼ 0.028.For the first alternative working point, there is a slight decrease in preservation of beam quality for small beam-beam tune shifts, as discussed in Sec.IV B. The general behavior is acceptable until ΔQ Tot ¼ 0.036.
For the second alternative working point, the limit is approximately ΔQ Tot ¼ 0.06.
The last two studies have been rerun with β Ã q ¼ 30 cm, the design value for the FCC-hh.This choice reduced the limit on the tune shift further.The limits for the LHC working point, (0.31, 0.32), and the first alternative working point, (0.315, 0.325), are collected in Table II.

V. EXPERIMENTAL RESULTS
A dedicated experiment was performed in the LHC to study the limitations due to strong head-on beam-beam interactions [38,39].The experiment tested collisions at different working points between individual bunches of higher intensity and smaller normalized emittance than the ones produced for regular operation with bunch trains, but close to the bunch brightness expected for the HL-LHC project.The resulting total beam-beam tune shift was just below 0.02.The experiment found that the loss rate depended strongly on the working point, which was moved along the diagonal close to the LHC working point.The parameters of the LHC were fixed with the standard setup of the 2017 run, except for the different working points.The particles were colliding in IP1 and IP5 only, with alternating full crossing angle of θ xing ¼ 280 μrad with β Ã q ¼ 0.4 m.The chromaticity was moderate, Q 0 ¼ 7.
The configurations of the bunches that were acted upon by the strongest beam-beam interactions have been simulated by CABIN.Each working point was tested in a limited window, making the emittance measurements erroneous.Therefore, only the loss rates will be compared.This experiment was run over two fills in the LHC.The first fill began at (0.311,0.318), whereupon the tunes of both beams were decreased simultaneously parallel to the diagonal.The second fill began at (0.311,0.319), whereupon the tunes of only beam 1 where increased one at a time.To use the implementations for the round beam-beam interaction, the average of the horizontal and vertical emittances has been used.To use the 4D implementation, the beam-beam parameter has been reduced by a factor S as in Eq. ( 17).The loss rates in the LHC and calculated with CABIN are visualized in Fig. 23.
In the LHC, the loss rate is large for working points close to the tenth order resonance, Q x ≤ 0.302.The loss rate is also large in the opposite end, for the working point (0.317,0.328).This is possibly caused by the 14th order resonance, 4Q x þ 10Q y ¼ 2. Because there is a nonzero crossing angle, it could also be caused by the 7th order resonance 2Q x þ 5Q y ¼ 1.There are also quite large loss rates close to (0.309,0.316).This might be caused by a combination of the 13th and 16th order resonances.The smallest loss rate was measured at a working point of (0.315,0.324).
The configurations were simulated with the round 4D beam-beam implementation in CABIN.For the working points close to the tenth order resonance, Q x ≤ 0.305, the calculated loss rates are high, comparable to the experimental values.The code calculates high loss rates for working points further from the tenth order resonance than what was measured in the experiment.For larger horizontal tunes, Q x > 0.305, the round 4D beam-beam implementation does not calculate high loss rates for any configuration.
The configurations were also simulated with the round 6D beam-beam implementation.For the working points close to the tenth order resonance at Q x ≤ 0.305, the calculated loss rates are high, as they were for the 4D implementation.For larger horizontal tunes, Q x > 0.305, the round 6D beam-beam implementation calculates loss rates similar to the experimental values.
The configurations were also simulated with the flat 6D beam-beam implementation.With this implementation, the large loss rates close to the tenth order resonance are measured for even larger horizontal tunes, Q x ≤ 0.306.For larger horizontal tunes, Q x > 0.306, the dependence on the working point is similar to the values calculated by the FIG.23.Loss rates measured in the LHC, subtracted the estimated luminosity burn-off of 5% h −1 , and loss rates (LR Beam ) calculated by CABIN for the same configurations using the different implementations of the beam-beam interaction.The circles and stars correspond to the first and second fill respectively.Even resonance lines up to 14th order are plotted on top of the LHC measurements, with width proportional to the strength of the resonance coefficient [40].The simulations are affected by the detrimental effects from the tenth order resonances close to Q x ¼ 0.3 at larger horizontal tunes than are measured in the LHC.The code counts particles as lost at lower amplitudes than where they actually are lost.Particles at 5σ Σ have tunes further from the working point than particles at 6σ Σ , as was detailed in Fig. 2.This difference could explain the discrepancy in the loss rates for Q x ∈ ½0.303; 0.305.Collimating at the larger amplitude has been found to push the limit for large loss rates down to Q x ¼ 0.303, but this also reduces the loss rates in the configurations far from the tenth order resonance significantly.Another possible cause is that the measured emittances in the machine during the experiment were erroneous.The beam-beam tune shift was also measured, and it suggested that the emittances of the strong beam were underestimated by 20% to 30%, which leads to an overestimation of the beam-beam tune shifts of beam 1 [38].
The loss rates calculated with the 4D implementation stand out for the working points with Q x > 0.306.The difference comes from the inability of the 4D implementation to model the effect of the crossing angle and the mixing with the longitudinal d.o.f.The reduction of the beam-beam parameter for the 4D model does reduce the length of the tune footprint, but it also reduces the width of the footprint, as the crossing angle does not do to the same extent.More importantly, the seventh and ninth order resonance lines with a nonzero crossing angle are replaced by the 14th and 18th order resonance lines with the 4D model.The 4D model is not able to correctly model the behavior at these working points, for which the beams can be affected strongly by odd resonances.
There are small differences between the round and the flat 6D simulations as well.Some of these differences are within the ∼10% deviation in the simulation results.The possibly most crucial difference is in the limit on Q x for large loss rates close to the tenth order resonance.The strong beam has a smaller vertical than horizontal normalized emittance, ϵ n;x2 > ϵ n;y2 , making the beam squeezed in the vertical direction and the beam-beam tune shift larger in the vertical tune than in the horizontal tune, in agreement with Eq. ( 5).In effect, the tune footprint with the flat implementation is shifted towards the coupling resonance.The total tune shift is also slightly larger for the flat implementation than the round implementation.This is caused by taking an average of the emittances when modifying the values of the strong bunch to be able to use a round beam-beam implementation.If the measured emittances were correct, the flat beam implementation would model the physics most accurately.It is possible that the values used for the round 6D implementation are closer to the actual values.

VI. CONCLUSION
The object of this paper has been to study the parametric dependence of detrimental mechanisms related to beambeam interactions.This has been done extensively by use of new beam quality simulations and FMAs.
In order to accurately quantify the impact on the beam quality, a high-performance tracking code named CABIN has been developed.It is designed to study the beam-beam interaction, which has been implemented by the weakstrong approach for both round and flat beams in both 4D and 6D.The modeling of the beam-beam interaction consumes in general more than 95% of the computation time of the simulations.The implementation using a GPU has reduced the tracking time by approximately 3 orders of magnitude, compared to a single CPU.A new regionally uniform initial distribution has been designed to model the evolution of both the core and the tail of the bunch accurately, while maintaining realistic computing requirements.It represents a 6D Gaussian bunch up to 6σ with the required accuracy, using only N mp ¼ 1 × 10 5 macroparticles.Especially the relative error of the loss rate measurement has been reduced from 86% to 8%, compared to a Gaussian distribution of the same number of macroparticles.Simulation results produced with this code showed good quantitative agreement with a dedicated experiment done at the LHC, and allowed to identify the role of the crossing angle in the observed loss mechanism.
Some effects have been shown to improve the beam quality for a given total beam-beam parameter.The best beam performance has been found for zero crossing angle, as the odd resonances are not activated.The crossing angle is generally nonzero in modern circular colliders, but may be forced zero for the head-on interactions by crab cavities in the HL-LHC and FCC-hh.Thus, the odd synchrobetatron resonances due to the head-on interactions could be suppressed.Increasing the crossing angle further past ϕ PIW ¼ 1 also reduces the detrimental effects.This improvement has been found to be caused by the decrease of the beam-beam tune shift and spread.As a result, the particles in the beams are affected by fewer strong resonances.Nevertheless, this regime is not preferable for operation, since the luminosity is also decreasing for larger crossing angles.
A transverse separation of more than 2σ also improves the beam quality by reducing the beam-beam tune shift.This is done to intentionally reduce the luminosity in a single experiment.However, for a nonzero separation of less than 2σ, the altered nonlinearities can deteriorate the beam quality faster than with zero separation.Using such separations to reduce the luminosity is therefore not an optimal strategy to minimize the impact of the beam-beam interaction on the beam quality.
A few specific configurations were also found to improve the beam quality without reducing the luminosity.It has been found that a significant hourglass effect reduces the detrimental effects caused by chromaticity, and vice versa.With a large chromaticity, Q 0 ¼ 15 at the LHC working point, (0.31,0.32), the optimal strength of the hourglass effect is when σ s =β Ã q ≈ 2=3, corresponding to β Ã q ¼ 12 cm in the FCC-hh, which is smaller than the current design value.This effect is only weakly dependent on the crossing angle.
Enforcing intermediate phase advances of μ 1 ¼ μ=2 or π=2 were expected from first order theory to reduce the effect of resonances.The symmetric phase advance, μ=2, improved the performance for all working points.An optimal intermediate phase advance was also found close to π=2, but it was different from the predicted value and found to depend on the working point.Both phase advances improved the beam quality best for zero crossing angle.The effect was marginal for a significant crossing angle, ϕ PIW ¼ 1.It was found that the accuracy required to keep the resonance coefficients suppressed was high.However, the effect on the beam quality remained under control with a reduced accuracy of order Δμ 1 ∼ 0.01, which is within reach using modern optics correction methods.
The LHC working point, (0.31,0.32), has been found to work well for small beam-beam parameters.With a Piwinski angle of ϕ PIW ¼ 1, a significant chromaticity, Q 0 ¼ 15, and a beam-beam tune shift equal to the maximum tune shift in the FCC-hh, ΔQ Tot ¼ 0.03, two alternative working points were found to better preserve the beam quality, (0.315,0.325) and (0.475,0.485).The beam quality was generally also found to be better closer to the coupling resonance.A realistic, pessimistic maximum beam-beam tune shift from the LHC working point, is found to be ΔQ Tot ¼ 0.043 with zero crossing angle.With a Piwinski angle of ϕ PIW ¼ 1, this limit is reduced to ΔQ Tot ¼ 0.028.The limits are larger for the alternative working point, (0.315,0.325), being ΔQ Tot ¼ 0.067 and 0.036 respectively.These values correspond to the optimum hourglass effect with β Ã q ¼ 12 cm.For the FCC-hh design parameters, β Ã q ¼ 30 cm at the LHC working point, the limits for zero and nonzero crossing angle are reduced to ΔQ Tot ¼ 0.035 and 0.018 respectively.While the ultimate beam-beam tune shift of 0.03 considered for the FCC-hh baseline with crab cavities seems within reach, the limit for nonzero crossing angle is significantly smaller indicating that the presence of odd resonances could be highly detrimental to its performance and should be taken into consideration.In particular, other sources of odd resonances, such as the imperfections of the magnetic lattice or long-range beam-beam interactions, have to be kept under control in order to avoid a further decrease of the maximum acceptable beam-beam tune shift.Further extensions of the numerical model have to be envisaged to evaluate quantitatively those effects.The assessment of the effect of the lattice nonlinearities is particularly challenging due to the extra numerical load of the lattice model.
Consistently with past studies, a working point close to the half integer was found to allow for significantly larger beam-beam tune shifts, supporting further investigations of the difficulties such as the optics correctability in this area of the tune diagram.
Consider an example relevant for simulations done in CABIN [21].Set the rms amplitude of the kick to be Δ ¼ σ p x ;0 × 10 −4 , and assume that the momentum of each particle is kicked once per turn, δt ¼ 1=f rev , where f rev is the revolution frequency.The evolution of the second moment of p x is then where T is the turn number, equal to t • f rev .For the normalized Gaussian distribution considered in this paper, σ p x ;0 ¼ 1.The growth rate of the second moment of p x is equivalent to the emittance growth rate expected from noise alone, _ ϵ 0 ¼ 1 × 10 8 =turn.ϵ 0 is the emittance based on normalized coordinates, equal to 2 for a perfect Gaussian distribution.

FIG. 3 .
FIG. 3. Comparison of initial distributions of N mp ¼ 1 × 10 5 macroparticles, with the actual χ 2 6 distribution.The sum of the weight of all macroparticles is N ¼ 1 × 10 11 .The bottom right plot displays the same distribution as the plot above, but with a linear scale.

FIG. 4 .
FIG. 4. Example of evolution of beam emittance in both transverse planes, ϵ x and ϵ y .The beam emittances are given based on normalized coordinates, being ϵ ¼ 2 for a perfect initial Gaussian distribution.The straight lines represent the linear regressions giving the emittance growth rates.

FIG. 11 .
FIG.11.FMAs illustrating activation of synchrobetatron resonances at the alternative working points, related to the beam quality study in Fig.10.Possible resonance lines up to the 16th order have been plotted on top of the tune footprint.A few labels were added to guide the reader.

FIG. 15 .
FIG.15.FMA illustrating canceling of resonances using an intermediate phase advance μ 1 ¼ μ=2, for a working point of ðQ x ; Q y Þ ¼ ð0.31; 0.32Þ, Q 0 ¼ 0 and ξ Tot ¼ 0.03.Possible resonance lines up to the 16th order have been added, in addition to a few labels meant to guide the reader.

FIG. 17 .
FIG. 17. FMAs illustrating canceling of resonances using an intermediate phase advance μ1 ¼ 0.512 • π, for a working point of ðQ x ; Q y Þ ¼ ð0.31; 0.32Þ, Q 0 ¼ 0 and ξ Tot ¼ 0.03.Possible resonance lines up to the 16th order have been added, in addition to a few labels meant to guide the reader.

FIG. 22 .
FIG.22.Beam quality reduction for increasing beam-beam parameter, ξ Tot , until a threshold is found.Simulations are run for a large chromaticity, Q 0 ¼ 15, and a significant hourglass effect, β Ã q ¼ 12 cm, changing the crossing angle, intermediate phase advance and working point.The outputs are presented as functions of ΔQ Tot ðϕ PIW Þ calculated by Eq. (17).
1Þ, where m is an integer, may be canceled if μ 1 ¼ fπ=2; 3π=2; …g.That is, resonances of order n ¼ f2; 6; 10; …g can be canceled simultaneously.There is no initial phase advance that can cancel resonances of order n ¼ 4 • m.

TABLE I .
Values for parameters in present and future circular colliders, in addition to the default settings of b May be countered by use of crab cavities.cAssumed collisions in two IPs, zero Piwinski angle in HL-LHC and FCC-hh.dSet for each simulation separately.These are the default values.e

TABLE II .
Maximum acceptable beam-beam tune shift, for pessimistic, realistic configurations.
The loss rates for large vertical tunes, Q y ≥ 0.327, are larger for the flat implementation than for the round 6D implementation and the experimental values.