Enhanced beam-beam modeling to include longitudinal variation during weak-strong simulation

Beam-beam interactions pose substantial challenges in the design and operation of circular colliders, significantly affecting their performance. In particular, the weak-strong simulation approach is pivotal for investigating single-particle dynamics during the collider design phase. This paper evaluates the limitations of existing models in weak-strong simulations, noting that while they accurately account for energy changes due to slingshot effects, they fail to incorporate longitudinal coordinate changes ($z$-variation). To address this gap, we introduce two novel transformations that enhance Hirata's original framework by including both $z$-variation and slingshot effect-induced energy changes. Through rigorous mathematical analysis and extensive weak-strong simulation studies, we validate the efficacy of these enhancements in achieving a more precise simulation of beam-beam interactions. Our results reveal that although $z$-variation constitutes a higher-order effect and does not substantially affect the emittance growth rate within the specific design parameters of the Electron-Ion Collider (EIC), the refined model offers improved accuracy, particularly in scenarios involving the interaction between beam-beam effects and other random diffusion processes, as well as in simulations incorporating realistic lattice models.


I. INTRODUCTION
Luminosity, the primary metric in accelerator physics for a collider, measures the rate of physics events per unit cross section per second during beam collisions.Achieving higher peak and integrated luminosity is pursued in colliders across energy [1], precision [2], and quantum chromodynamics froniters [3].
In circular colliders, beam-beam interactions-the electromagnetic forces between opposing beams-significantly limit performance.For lepton colliders, this interaction introduces a "beam-beam limit," where luminosity increases below the expected relationship with beam intensity beyond a certain threshold [4].In hadron colliders, beam emittance growth, primarily driven by these interactions, curtails luminosity lifetime, as demonstrated in empirical studies [5,6].The "beam-beam parameter," quantifying the maximum tune shift during collisions, serves as a key metric for assessing interaction strength.The Electron-Ion Collider (EIC), to be constructed at Brookhaven National Laboratory, aims to achieve an unprecedented luminosity of 10 34 cm −2 s −1 by colliding 10 GeV electrons and 275 GeV protons, requiring large electron and proton beam-beam parameters [7].Throughout the design phase of the EIC, a critical objective is to minimize proton emittance growth amidst beam-beam interaction.
Hadron-Elektron Ring Anlage (HERA) pioneered the collision of leptons and protons, serving as the first operational collider of its kind.However, in comparison to the main parameters achieved during routine operation at HERA [8], the EIC seeks to escalate the peak luminosity objective by two orders of magnitude, accompanied by a fourfold increase in both proton and electron beam-beam parameters.Given the unprecedented nature of the EIC's design objective, beam-beam simulations emerge as the sole method to validate the efficacy of the parameter combination, with accurate modeling of beam-beam interaction being crucial for assessing proton emittance growth in the simulation.
In beam-beam simulations, two principal approaches are distinguished: strong-strong and weak-strong.The strongstrong concept was introduced by Peggs in the 1980s [9].In this paper, the strong beam distribution is fixed during tracking, and the role of the strong beam were switched between the interacting beams over successive tracking turns.Modern strong-strong simulations have advanced to a point where the beam distributions are dynamically updated based on the electromagnetic fields calculated from a two-dimensional Poisson equation.Subsequently, these distributions evolve according to the Vlasov equation, ensuring a self-consistent simulation process.
In contrast to the strong-strong approach, weak-strong simulations simplify the modeling of beam-beam interactions by assuming the strong beam maintains a rigid Gaussian distribution.This simplification facilitates the calculation of electromagnetic forces exerted on test particles within the weak beam using the Bassetti-Erskine formula [10], a method that, despite lacking self-consistency and the ability to address coherent effects, remains crucial for collider design due to its lower computational demands.Particularly in hadron storage ring design, where effective cooling methods are absent and particles must be tracked over millions of turns to ascertain stability, the weak-strong model proves indispensable.Moreover, it introduces less numerical noise compared to the strong-strong method, which employs the particle-in-cell technique to expedite Poisson equation solutions [11][12][13].This technique, while efficient, often leads to significant numerical noise by projecting particle distributions onto a finite twodimensional grid, as observed during the EIC design phase [14][15][16][17][18].Such noise can obscure the particle diffusion effects attributed to beam-beam interactions, making them discernible only through weak-strong simulations.
For the EIC, the adoption of a flat hadron beam configuration-characterized by a vertical emittance that is an order of magnitude smaller than horizontal emittance-maximizes luminosity while introducing sensitivity to various real-world fluctuations in the vertical plane.The weak-strong simulation has emerged as a cornerstone in addressing EIC design challenges, including the dynamic aperture reduction from interaction region (IR) magnetic field errors [19], emittance growth due to electron orbit perturbations from dipole magnet power supply variability [20], and crab cavity phase noise [21,22].The necessity for a precise weak-strong model transcends the EIC, underscoring its critical importance across the collider physics community for ensuring the integrity of collider design.
This paper is organized as follows.Section II revisits the concept of synchro-beam mapping, highlighting the z variation effect, and the energy changes resulting from slingshot effects.Section III introduces two symplectic approaches designed to encompass both the z variation effect and the slingshot effect-induced energy changes.Section IV provides a comparative analysis of the simulation outcomes using these three distinct methodologies.A conclusive summary is presented in Sec.V.

II. SYNCHRO-BEAM MAPPING
A single-particle dynamics is described by the set of canonical coordinates: x ¼ ðx; p x ; y; p y ; z; p z Þ; ð1Þ where x and y represent the transverse positions in the horizontal and vertical planes, respectively, and p x;y are their associated momenta, normalized to the design momentum, P 0 , of the reference particle.The longitudinal position is given by z ¼ s − l, where s is the designed path length for the reference particle, and l is the actual path length.The term p z ¼ ðP − P 0 Þ=P 0 denotes the deviation in momentum relative to the reference particle.For the purposes of clarity and focus, this paper restricts its scope to the case of high-relativistic dynamics.
The choice of canonical coordinates and the Hamiltonian is discussed in Appendix A.

A. Hirata's original approach
The synchro-beam mapping, formulated by Hirata, Moshammer, and Ruggiero, is extensively utilized for simulating beam-beam interactions in the presence of synchrotron motion, providing a symplectic mapping within a six-dimensional phase space [23].This method involves longitudinally splitting the strong bunch into multiple slices and employing a drift-kick-drift model to calculate the particle-slice interaction.In this approach, the particle engages in a virtual drift from the interaction point (IP) to the collision point (CP), receives a beam-beam kick from the opposing slice at the CP, and then returns back to the IP.
Denoting z Ã the longitudinal coordinate of the strong slice in its own axis, the separation between the IP and CP is represented as Throughout this paper, we will use S to represent Sðz; z Ã Þ without causing any ambiguity.The virtual drift from the IP to the CP is characterized by an exponential Lie operator D 0 , defined as This operator's action on the canonical coordinates results in the transformation: The beam-beam interaction at the CP is represented by another Lie operator B, expressed as where U denotes the beam-beam potential induced by the strong slice.The application of the operator B to the coordinates results exclusively in modifications to the momenta, described by The spatial coordinates remain unaffected by the beambeam interaction, given the assumption that the slice is sufficiently "thin."Transverse momentum changes, Δp x;y , are derived from the well-established Bassetti-Erskine formula, and the formulation for the longitudinal momentum change, Δp z , follows Hirata's proposal in [23].
Appendix B provides detailed formulas for momentum changes by a strong slice characterized by a bi-Gaussian distribution.
Upon experiencing the beam-beam kick at the CP, the test particle is subjected to another virtual drift, leading it back from the CP to the IP.This sequential process is encapsulated by the overall mapping equation: which represents the combined effect of the initial drift, the beam-beam interaction, and the subsequent return drift.The resultant transformation of the particle's coordinates and momenta is given by where the new positions and momenta are adjusted according to the shifts induced by the beam-beam interaction.

B. z variation and slingshot effects
Hirata's approach provides a symplectic mapping for the beam-beam interaction, which includes the energy change and the bunch-length effect.This approach employs exponential Lie operators to represent both the virtual drift and the beam-beam kick, ensuring that each component of the synchro-beam mapping maintains symplecticity.Consequently, the cumulative mapping that results from sequentially applying these steps preserves the symplecticity.
However, it is imperative to acknowledge a critical distinction: the transformation characterizing the transition between the IP and the CP does not retain symplecticity when considered in the context of the standard accelerator coordinates ðx; p x ; y; p y ; z; p z Þ.This discrepancy arises from the dynamic nature of the CP location, which varies among particles based on their initial coordinates at the IP.Specifically, some particles engage in collisions at CP1, others at CP2, and so forth, necessitating the computation of particle coordinates at their respective CPs.Once all particles have transitioned to their respective CPs, the overarching transformation from the IP to these dynamically defined CPs does not conform to symplectic principles.
To illustrate this concept, we can incorporate Eq. ( 2) into an actual drift map.Consider the simple drift Hamiltonian given by H 0 ¼ ðp 2 x þ p 2 y Þ=2.The transformation corresponding to a real drift over a distance S is described by The corresponding Jacobian matrix is given by M ¼ ∂ðx; p x ; y; p y ; z; p z Þ ∂ðx 0 ; p x;0 ; y 0 ; p y;0 ; z 0 ; p z;0 Þ ¼

:
It is straightforward to verify that where J is the 6-by-6 symplectic form matrix.This deviation from symplecticity can be intuitively understood through an analogy to photography in phase space.According to Liouville's theorem, in a closed and isolated system, where dynamics are governed by Hamiltonian principles without any dissipation or diffusion, the volume occupied by the system in phase space remains invariant over time.When all particles are captured in a single "photograph," the phase volume should remain constant.However, when we capture only a subset of particles in sequential "photographs," the volume of phase space represented in each image may vary depending on the distribution evolution.Consequently, when summing the volumes from these sequential photographs taken at different times, the total phase volume may not remain constant.
Easing off the symplectic constraint, Fig. 1 presents the actual drift of a particle between the IP and the CP.The respective path lengths are quantified as follows: where x 0 ¼ dx=ds and y 0 ¼ dy=ds denote the derivatives of x and y with respect to the path length s, while Δx 0 and Δy 0 represent the deflections due to the beam-beam interaction at the CP.This difference in path length accounts for the alteration in the z coordinate as the test particle returns to the IP: In this paper, the longitudinal coordinate change is referred as z variation effect.When every particle undergoes a beam-beam kick at the CP and returns to the IP, the entire process as a whole must retain symplecticity.However, given the actual drift between the IP and the CP deviates from symplectic behavior, the beam-beam interaction at the CP must also be nonsymplectic in nature to ensure the overall process conserves symplecticity.
In addition to the momentum change or energy change induced by distribution variation, as shown in Eq. ( 6), there is an additional mechanism of energy alteration reminiscent of the gravitational slingshot effect encountered in orbital mechanics, as explained in [24].The comprehensive derivation of this effect is provided in Appendix C. In terms of the accelerator coordinates, the momentum change due to the slingshot effect is By comparing Eqs. ( 11) and ( 12) with Eq. ( 8), we observe that while Hirata's approach encapsulates the energy change attributable to the slingshot effect, it does not account for z variation.In the original formulation by Hirata, the z coordinate is presumed to remain unchanged throughout the interaction.This assumption, however, does not fully capture the dynamics of the system.The effect of z variation may have significant implications in the design, especially given that the hadron particles are often tracked over millions of turns.Recognizing and addressing this z variation is crucial for accurately modeling the long-term behavior of particles within accelerators, pointing to areas where existing theoretical frameworks might be enhanced to better reflect physical realities.

III. MODIFIED VIRTUAL DRIFT
Despite the nonsymplectic nature of the actual transformation from the IP to the CP in practical scenarios, Hirata's approach retains its significance for two key reasons: first, it ensures the symplecticity of the overall synchro-beam mapping, and second, it accurately captures the slingshot effects.Building on this foundation, we aim to extend Hirata's methodology in this section by devising an approximate yet symplectic transformation for the virtual drift.This transformation will account for both z variation and slingshot effects at the lowest order, thereby refining the model to more closely mirror the complexities observed in Sec.II B.

A. Chromatic Hamiltonian
Incorporating p z into the Hamiltonian expansion allows for calculation of z variation effects.This leads to the refined Hamiltonian: Following Hirata's strategy, the associated exponential Lie operator is defined as We remind the reader that the Lie operator is a shorthand notation for the Poisson bracket.For any two functions of f and g, the Poisson bracket is defined as 1. Illustration of test particle drift between the IP and CP.The trajectory from IP to CP is depicted by the blue line, whereas the return path from CP to IP is represented by the red line.
Application of the ∶h 1 ∶ operator to longitudinal coordinates where Iterative applications of the operator reveal for n ≥ 1.Here ð•Þ n is the Pochhammer symbol of rising factorial The resulting transformations for z and p z are Utilizing the identity the transformation of the virtual drift from the IP to the CP can be expressed as The transformation from the CP to the IP reverses the effects of D 1 , achieved by inverting the sign of S and S 0 : The operator D 1 in Eq. ( 14) implies that the location of the CP is determined by the longitudinal coordinate z at the IP.During the drift from the IP to the CP, z is no longer a constant, as shown in Eq. (22).This variation means that the test particle may not engage with the counterpropagating slice at the initially calculated CP location.To effectively apply the collision map B, delineated in Eq. ( 5), we have to presume that the z variation-arising from the modified virtual drift-exerts negligible influence on the collision map B.
This assumption, despite not being entirely realistic, is not anticipated to introduce significant discrepancies.The electromagnetic field generated by a relativistic beam is characterized by a narrow angular width of 1=γ Ã , where γ Ã denotes the Lorentz factor of the strong slice.Let σ represent the transverse size encountered by a test particle interacting with this strong slice.The test particle will predominantly experience the electromagnetic field from the thin slice if: In most practical scenarios, this condition is generally met.

B. Exact Hamiltonian
By applying an exponential Lie operator to obtain the virtual drift transformation, we have extended Hirata's methodology for the chromatic drift Hamiltonian.
However, it proves to be impractical for the exact drift Hamiltonian due to the resulting expression's complexity and unwieldiness.Therefore, it is imperative to explore an alternative strategy to effectively tackle this challenge.
In Sec.III A, it is demonstrated that the CP location is not solely determined by the initial longitudinal coordinates.Here, we continue to use S to represent the drift length from the IP to the CP.To account for geometric effects, S must be contingent on the transverse momenta p x;0 and p y;0 , and to consider the chromatic effect, S is also influenced by p z;0 .The superscript "0" means these quantities are evaluated at the IP.As a result, the relationship governing S can be expressed as follows: The variable z Ã remains constant when the thin slice traverses between the IP and the CP.At the CP, the test particle has a longitudinal coordinate z.Although the expression of z is still unknown yet, S can be expressed by z and z Ã .The elapsed time for the reference particle from the IP to the CP is t r ¼ S=v r .The superscript "r" indicates that v r is the velocity of the reference particle.Using the definition of z in Eq. (A7), where t is when the test particle arrives at CP.For the opposite slice, the arrival time is obtained by substituting S with −S, where the negative sign comes from the opposite s axis: The test particle will collide with the slice when t ¼ t Ã : Taking the relativistic limit v ¼ v Ã ¼ c, Eq. ( 28) will turn into Eq.( 2).It is crucial to note, however, that the interpretation of the variable z varies between the two equations.In Eq. ( 28), z represents the longitudinal coordinate at the CP, whereas in Eq. ( 2), the same symbol z actually denotes the quantity z 0 determined at the IP.As a result, S has a form of S ¼ Sðz 0 ; z Ã ; p x;0 ; p y;0 ; p z;0 Þ ¼ Sðzðz 0 ; p x;0 ; p y;0 ; p z;0 Þ; The spatial coordinate transformation governed by the exact drift Hamiltonian is detailed as where The variable z on the right-hand side in Eq. ( 30) is evaluated at the CP, and z can be resolved by combining Eqs. ( 28) and (30).
We introduce the generating function to find out the transformation of p x;y;z : The spatial coordinates at the IP are given by This relationship confirms the alignment of the spatial coordinate transformation equations with those previously established in Eq. (30).The transformation of momentum coordinates at the CP is similarly derived from the partial derivatives of G 3 with respect to the spatial coordinates, yielding: Similar to the Hirata's approach, the transverse momenta remain unchanged, and the correction is applied to p z to preserve the symplectic structure of the system.The reverse transformation, from the CP back to the IP, is achieved by representing the initial coordinate and momentum set ðx 0 ; p x;0 ; y 0 ; p y;0 ; z 0 ; p z;0 Þ in terms of the final set ðx; p x ; y; p y ; z; p z Þ, utilizing a combination of Eqs.(30) and (34).
For relativistic case, the transformation from the IP to the CP is summarized as and the reverse transformation from the CP to the IP: p s;0 S; p y;0 ¼ p y ; y 0 ¼ y − p y p s;0 S; Section III details the derivation of the enhanced virtual drift, incorporating z variation and energy adjustments due to slingshot effects.Weak-strong simulations are conducted to validate the enhanced map's accuracy and assess its implications, ensuring the theoretical modifications are both substantiated and practically relevant.

A. Model difference demonstration
To demonstrate the distinctions among the three models, simulations are performed using artificially constructed parameters.These "man-made" parameters, which do not correspond to real-world machines, are specifically designed to highlight the differences between the models.The beam-beam parameters are set at 0.01, with the associated beam characteristics presented in Table I.
Figure 2 illustrates the variation in the longitudinal coordinate z after the test particle undergoes a virtual drift from the IP to the CP, receives a beam-beam kick at CP, and drifts back to IP.The initial coordinate of the test particle is set at ðσ; 0; σ; 0; 0; 0Þ T , where σ ¼ 70 μm.As the location of the thin slice changes, so does Δz.For Hirata's model, Δz ¼ 0. In both the chromatic and exact models, Δz is aligned with each other but remains small with a magnitude of 10 −11 m.
The z variation is too small to impact beam dynamics observably when the beam-beam kick is not sufficiently large.Figure 3 depicts the frequency diffusion over the course of tracking 1000 turns.The thin slice is positioned at z Ã ¼ −β=2 ¼ 30 cm.The initial particle distribution spans AE5σ across the x-y plane, with the longitudinal coordinates and p x;y set to zero.The diffusion index D is determined using the formula: Here, ν 1 and ν 2 are the frequencies computed from the first and second 500 turns, respectively.The results across all three models show similar frequency maps.A seventh order resonance 4ν y ¼ 7 can be observed.
To enhance the z variation, one effective approach is to amplify the beam intensity.The beam-beam parameters can be preserved by reducing the β functions and maintaining the slice position at z Ã ¼ −β=2.Suppose the beam intensity is increased by a factor of A, and concurrently, the β functions are reduced by the same factor.This adjustment results in an A-fold increase in both beam divergence and the beam-beam kick.According to Eq. ( 11): where x 0 , y 0 , Δx 0 , and Δy 0 also increase by a factor of A, while the distance S between the CP and IP is reduced by a factor of A. Consequently, it leads to an A-fold increase in Δz.
In our case, the initial longitudinal coordinates are set to 0. Therefore, the longitudinal action is derived from the energy kick due to beam-beam interactions.When the beam intensity is scaled up, the associated longitudinal energy kick also increases proportionally.Consequently, the change in the longitudinal coordinate, Δz, remains negligible compared to the longitudinal action.As the β functions decrease during this scaling process, the hourglass effect becomes relevant when Δz approaches the magnitude of the β function.This variance in Δz across different models will influence long-term particle diffusion.
Figure 4 displays the frequency maps obtained from different models when the beam intensity is increased by a factor of 100.The reduction in β functions leads to fewer particles, indicated in red color, being affected by the seventh order resonance.In the chromatic and exact virtual models, particles close to the center exhibit a higher diffusion index.Therefore, in long-term tracking, the z variance effect should be included.

B. Realistic case: Simulation for EIC
The parameters used in simulations are shown in Table II.The crossing angle is as large as 25 mrad, and the Lorentz transformation is applied to accurately address the crossing angle within the simulation [26].In all subsequent simulations, the weak proton beam is represented by one million macroparticles.A second order harmonic crab cavity is used to flatten the proton bunch [27].The one-turn map at IP is represented by the linear betatron map.There is no momentum dispersion and crab dispersion at IP [28].The electron beam, functioning as the strong beam, exhibits a rigid bi-Gaussian distribution.Parameters for the electron beam, as listed in Table II, FIG. 3. Frequency map analysis in amplitude space (top) and tune space (bottom).FIG. 4. Frequency map analysis in amplitude space when the beam intensity is scaled by 100 times.are anticipated to be achieved upon reaching equilibrium, accounting for the effects of beam-beam interactions.
Figure 5 illustrates the evolution of proton emittance through three distinct mapping approaches.The blue curve represents the original Hirata's mapping, detailed in Eq. ( 8).Meanwhile, the orange curve displays the results obtained from a modified virtual drift model incorporating the chromatic Hamiltonian, as derived from Eqs. ( 22) and (23).Finally, the green curve showcases the outcomes derived from employing a modified virtual drift model with the exact Hamiltonian, as outlined in Eqs. ( 35) and (36).
The emittance growth rate is determined through a linear fit of the final 50% of the tracking data.To more clearly illustrate the trend of emittance evolution, the data are averaged over every 1000 turns.In the context of the EIC design, vertical emittance growth is of particular concern.For the baseline design, it is imperative that the vertical emittance growth does not surpass 20%/h, a limit set by the Strong Hadron Cooling requirements.As depicted in Fig. 5, the vertical emittance growth observed across the three examined approaches shows negligible differences.This outcome is expected, given that the z variation effect constitutes a higher-order effect, and the parameters listed in Table II have been fine-tuned to minimize vertical emittance growth.
Figure 6 presents the normalized relative error derived from tracking data, offering valuable insights into our enhanced model.The linear progression of relative error over time suggests that the proton beam is free from significant resonance phenomena.Furthermore, the error's magnitude is notably small, on the order of 10 −5 , affirming the effectiveness of Hirata's original mapping across a TABLE II.Simulation parameters from EIC-CDR [7]."H" stands for horizontal and "V" denotes vertical below.The electron beam, acting as the strong beam, is characterized by design parameters anticipated to be realized upon achieving equilibrium.

Parameter
Proton  The legend describes the virtual drift models employed in the simulations: original Hirata's method (blue), chromatic Hamiltonian (orange), and exact Hamiltonian (green).
majority of scenarios.Additionally, a comparison between the modified drift model, which integrates the chromatic Hamiltonian, and the exact Hamiltonian model reveals a remarkable degree of consistency, highlighting the robustness of our approach.
It is important to note that a relative error on the order of 10 −5 or 10 −4 per turn may not be considered negligible when compared to the magnitude of random diffusion processes, such as intrabeam scattering (IBS).In the context of the EIC, where the horizontal or longitudinal IBS growth time is approximately 2-3 h, the relative amplitude of IBS diffusion is about 10 −7 per turn.From this perspective, the inclusion of the z variation effect in beam-beam interaction becomes essential in the design of the hadron ring, when addressing the interplay between beam-beam interaction and the IBS.
In accelerator physics, it is a well-established fact that minor deviations can be exponentially magnified over time for chaotic particle motion.This scenario typically arises when there is an overlap of multiple resonance lines within the tune space.Specifically, for the EIC beam-beam simulation, high-order synchro-betatron resonances are present within the footprint, as detailed in [29].Figure 7 illustrates the tracking outcomes when the strong electron beam is artificially displaced by one σ x in the horizontal plane.This displacement excites the higher-order synchrobetatron resonances, leading to a noticeable increase in vertical emittance across three distinct simulation models.
Incorporating a realistic lattice model inevitably introduces more higher-order resonances, particularly at larger amplitudes.Figure 7 underscores the necessity of incorporating z variation in the dynamic aperture study.

V. CONCLUSION AND OUTLOOK
In weak-strong simulations for beam-beam interactions, Hirata's model is widely used, which offers an approximate, yet symplectic mapping that accounts for energy changes and bunch length effects.However, it operates under the assumption of constant longitudinal coordinates (z), thereby excluding z variation from its considerations.
In the design of hadron rings, it is essential to track particles across millions of turns, making it impractical to disregard z-variation during beam-beam interactions.Building upon Hirata's pioneering work, we propose two new models designed to refine the virtual drift process.It can be shown that these models are symplectic and accurately integrate the z-variation effect at the lowest order.
Simulations are performed to benchmark these three distinct models.Although the result indicates that z variation is a higher-order effect that does not substantially alter the emittance growth rate under the specific design considerations of the EIC, the discrepancy among the models is significant when contrasted with the amplitude of IBS diffusion.This distinction becomes even more critical when considering realistic lattice models that involve higher-order resonances, underscoring the importance of accounting for z variation.
Moreover, the approach discussed in Sec.III B can be readily adapted to include external fields, such as those from a detector solenoid, without having to assume that the z coordinate remains constant in the presence of such fields.
We have to emphasize that the enhanced models presented in Sec.III, along with Hirata's original model, should be regarded as approximations.These models presuppose that the transformations between the IP and the CP adhere to symplectic principles with respect to accelerator coordinates, a notion that deviates from physical reality as discussed in Sec.II B. Nonetheless, we also recognize that in the absence of dissipation and diffusion, real-world physical processes must inherently conform to symplectic characteristics.
An alternate approach involves standard tracking (no energy kick) between IP and CP with the distance between IP and CP determined by the longitudinal coordinates z of the test particle and the z Ã of the opposing slice along with the trajectory of the test particle.This is coupled with a beam-beam interaction at the CP that includes an energy kick due to the slingshot effect and an energy kick due to the varying shape of the slice.This strategy has been applied within the framework of Bmad [24].
Shatilov and Zobov in [30] also employed this strategy: they used the real drift map and modify the beam-beam kick accordingly.In their paper, the additional energy kick is attributed to v x B y − v y B x , where v x;y are transverse velocities and B x;y are transverse magnetic field from beambeam interaction.In Bmad and in Appendix C, this energy kick is derived by analogy to the slingshot effect.Both treatments yield consistent formulas under the highrelativistic approximation.
The advantage of this method is that it more closely mirrors what is actually happening.The disadvantage is that symplecticity of the overall map is not assured.Future work will entail a detailed comparison between the methodologies employed by Hirata and the one used by Bmad and Shatilov.

APPENDIX A: HAMILTONIAN AND CANONICAL VARIABLES
In accordance with Forest [31], the Hamiltonian for a relativistic particle in an external magnetic vector potential A, navigating through a curved coordinate system of radius ρ, is formulated as Here, δ ¼ ðP − P 0 Þ=P 0 represents the relative momentum deviation from the reference momentum P 0 , and p x , p y denote canonical momenta normalized by P 0 .The canonical variables in this context are represented by two equivalent sets: ðx; p x ; y; p y ; −l; δÞ or ðx; p x ; y; p y ; δ; lÞ; ðA2Þ with l denoting the particle's path length.The reference particle's design path length, s, serves as the independent variable in the Hamiltonian K expressed in Eq. (A1).
The generating function yields the transformations: and leading to a new Hamiltonian formulation: The corresponding canonical variables are x;p x ¼ P x P 0 ;y;p y ¼ P Similar to [26], a constant 1 is added to the Hamiltonian.It has no effect on beam dynamics.For a particle in a drift space, where ρ ¼ ∞; A s ¼ 0, the Hamiltonian simplifies to

APPENDIX B: BEAM-BEAM KICK FOR BI-GAUSSIAN DISTRIBUTION
For a bi-Gaussian particle distribution, the beam-beam potential U is given by Uðx; y; σ where N is the total particle number, r 0 ¼ e 2 =ð4πϵ 0 mc 2 Þ the classical radius, γ 0 the relativistic factor of the test particle, Q 1;2 the charge numbers of particles from two colliding bunches, and σ x;y are the rms beam sizes of the strong slice at the CP.The parameters σ x and σ y are not constants but vary as functions of the distance S from the IP to the CP.The first order differential was derived by Bassetti and Erskine [10] where U x ≡ ∂Uðx; yÞ ∂x ; U y ≡ ∂Uðx; yÞ ∂y x − σ 2 y Þ q : In Eq. (B2), wðzÞ is the Faddeeva function defined as Its derivative is given by The parameters σ x;y vary as functions of the longitudinal coordinate z.The derivative U z was first given by Hirata in [23]: Taking the partial derivatives on both sides of Eq. (B2), the second order derivatives can be obtained As a result, the momentum changes in Eq. ( 6) for a strong slice with bi-Gaussian distribution are fully resolved.

APPENDIX C: ENERGY CHANGE DUE TO A MOVING MAGNET
The four-momentum of a test particle in the laboratory frame is denoted as ðE=c; PÞ, where E represents the energy of the particle, c is the speed of light, and P is its momentum vector.When transitioning to the rest frame of the moving magnet, the energy of the test particle undergoes a transformation as per the principles of Lorentz transformation.This transformation is articulated as where β Ã ¼ v Ã =c is the magnet's velocity normalized by the speed of the light, and γ Ã is the corresponding Lorentz factor.θ represents the angle between the particle's momentum vector and the opposite direction of the magnet's motion.The notation Ē indicates the energy measured in the magnet's rest frame.
In this rest frame, the test particle's trajectory is altered while its energy is conserved.Therefore, For relativistic scenarios, we approximate where subscripts "1" and "2" denote the states before and after the interaction with the magnet, respectively.Substituting these approximations into Eq.(C2) yields the energy change: Utilizing the paraxial approximation where θ 1;2 ≈ 0 and considering the high relativistic limit where β Ã ≈ 1, the formula simplifies further to

FIG. 7 .
FIG. 7. Proton emittance evolution in weak-strong simulation with strong electron beam offset by one σ x in the horizontal plane: top (horizontal plane) and bottom (vertical plane).Here, σ x represents the rms horizontal size of the electron beam at the IP.The legend follows the same conventions as in previous figures.

TABLE I .
[25]ficial parameters used to illustrate z variations among different virtual drift models.The transverse working point is based on the SuperKEKB design[25].
FIG.2.Longitudinal coordinate variation after projecting back to IP.The legend describes the virtual drift models employed in the simulations: original Hirata's method (blue), chromatic Hamiltonian (orange), and exact Hamiltonian (green).The orange and green curves overlap with each other.