Introduction

The Electron-Ion Collider in China (EicC) [1] is a proposed highly polarized electron-ion collider based on the high-intensity heavy-ion accelerator facility (HIAF) [2] upgrade. It provides a platform for frontier research in nuclear physics with a center of mass energy of 16.7 GeV and luminosity of 2×10³³ cm^-2s^-1. The primary parameters are presented in Table 1. The proton ring (pRing), one of the major accelerators of the EicC project, is both the accelerator ring and the collider ring of the proton beam. In the pRing, an operation mode with multi-bunches and high single-bunch intensity is essential to achieve the design goal of high luminosity. However, this leads to enormous beam dynamics challenges in the pRing. Intensity-dependent effects, such as the beam coupling impedance [3, 4], space charge effect [5], electron cloud instability [6, 7], beam–beam interaction [8-10], and their interplay [11-14], are the most severe limitations of EicC performance.

To thoroughly investigate the machine performance, beam instabilities, and associated mitigation methods subjected to complex high-intensity effects, macroparticle tracking is the only way to model and optimize the beam dynamics in pRing. Various beam dynamics simulation codes have been developed by CERN according to specific requirements. For example, PyHEADTAIL [15] is developed for electron cloud instability and impedance induced single-bunch collective effects, PyECLOUD [16] is exploited for electron cloud establishment simulation, and the longitudinal collective effect and beam manipulation are implemented in BLonD [17]. Beam–beam interactions in colliders can be studied using BeamBeam3D [18] and Athena [19] developed by LBNL and IMP, respectively. These codes have been well benchmarked with different existing codes or beam-based measurements and have become powerful tools for various beam dynamics studies. However, as the machine performance is continuously pushed to a higher level, different effects can no longer be treated independently and their complex interplay should be considered in any realistic attempt to study the high-intensity beam dynamics [12, 20]. The combination of different codes for cross-talk studies is impracticable because of differences in coding languages and common interfaces. Almost no code can simulate multiple high-intensity beam dynamics in hadron accelerators in a self-consistent manner. Therefore, as part of the EicC project, a new beam dynamics simulation code, GOAT, has been developed in the past few years and is capable of studying different single-bunch intensity-dependent beam dynamics, their complex interplay, and possible mitigation techniques simultaneously in the pRing.

GOAT is a multiparticle tracking code developed based on the Python [21] language using object-oriented programming (OOP) technology. Python significantly improves the efficiency of code development. Some of the core computing modules of the code are written using the Cython [22] package because of the slow execution speed of the interpreted language. Taking advantage of OOP technology, diversified elements and functions can be easily inserted into the code without affecting the existing ones. This also facilitates further parallelization to achieve better performance. Thus far, based on the most basic beam transportation and manipulation methods in the transverse and longitudinal planes, implementations of the impedance induced single-bunch collective instability, space charge effect, electron cloud effect, and beam-beam interaction are included in the code.

The remainder of this paper is organized as follows. The code architecture and numerical model are given in sec2. In sec3 to sec6, the detailed implementations of the four different high-intensity effects are explained separately, and the benchmark results against other codes or theories are presented and discussed. An application of the code for the cross-talk between two different effects is presented in sec7. sec8 presents the summary and outlook. As an important aspect of simulation software, the code performance is preliminarily tested and the results are summarized in appdix.

Code architecture and numerical model

GOAT consists of a beam module, data management module, infrastructure element module, and physical element module. The code architecture is illustrated in Fig. 1. All accelerator beam particles are defined in the beam module. In this module, different particle distributions, such as KV and Gaussian distributions, can be generated according to the parameters given by the user [23]. Various methods for transforming the coordinates of the particles and other parameters are implemented in the beam module, which can provide suitable parameters for specific studies. In addition, the beam module includes functions for the initialization and dynamic management of the electrons distributed in the vacuum chamber used for electron cloud simulation.

Fig. 1Full-screen View

Code architecture of GOAT.

The data management module can read external input files, such as beam distribution and optics files according to user commands, and convert them into data structures compatible with GOAT. The data, such as turn-by-turn particle statistical information and instantaneous particle distribution at specified time steps during the simulation process, can be stored by this module and written to the disk with time stamps. This module also includes data post-processing functions, such as reading the output files from the simulation, picking up the required data, and automatically generating charts.

Infrastructure element module and physical element module constitute the physical kernel of the GOAT code. The former typically contains elements that provide auxiliary functions for simulations, such as beam slicing elements and PIC method-based Poisson solvers. Intensity-independent beam transportation and manipulation can also be achieved by combining infrastructure elements. The latter specifically refers to the four types of high-intensity beam dynamics elements implemented in GOAT. Each of the four elements is used to describe a particular physical process. Based on the OOP technology, it is convenient to build new elements in both the infrastructure and physical element modules. Instead of using a command file to control the simulation procedure [16], the user can interactively customize the beamline in a specific order by using different combinations of infrastructure and physical elements to complete the desired simulation [15, 17].

In GOAT, a macroparticle beam of intensity N, energy E₀, and particle charge q is described by the 6D set of coordinates (x, px, y, py, z, pz), where x denotes the horizontal offset from the reference orbit, y the vertical offset, px and py the corresponding transverse normalized momenta, z the longitudinal offset from the synchronous particle, and pz the relative momentum deviation [24]. As for an electron macroparticle used in the electron cloud study, a 7D vector (x, vx, y, vy, s, vs, ne) is constructed [6, 16], where x and y are defined in the same way as a beam particle, s denotes the distance to the transverse reference plane, vx, vy, and vs the absolute velocity component in the horizontal, vertical, and longitudinal direction, respectively. Further, ne denotes the amount of charge carried by the macroparticle. Hereafter, for convenience, “beam” and “electron” are used to refer to an accelerator beam and electron in the electron cloud study, respectively.

Because of the distinctive motion characteristics of the beam particles and distributed electrons (the former always fulfills ultra-relativistic conditions, whereas the latter is usually non-relativistic), different methods are required to integrate the corresponding equation of motion. GOAT models the transverse beam dynamics by tracking macroparticles linearly through a transfer matrix between interaction points around the ring [24]. Longitudinal beam motion can be modeled by either linear tracking or nonlinear drift-kick integration through RF elements with a given voltage, phase, and frequency [25]. However, for an electron, the corresponding coordinates are advanced by integrating the Lorentz equation that it follows directly. A second-order symplectic Boris algorithm is implemented in the code [25], and in certain specific cases, the integration algorithm can be simplified for better performance [7].

In addition, GOAT is equipped with two beam slicing methods: one is the uniform-length slicing method, in which all slices have the same length, and the other is the uniform-charge slicing method, in which each slice contains the same number of macroparticles [15]. In each method, slicing is performed for the beam rather than the ring, and the upper and lower limits of the longitudinal extent are given by the user. Both methods have their own advantages and can be used for different simulations. Additionally, to meet the requirements for solving the Poisson equation, two types of solvers based on the particle-in-cell (PIC) method are implemented. One uses the finite-difference (FD) time-domain method with perfect electric conducting boundary conditions [6]. The vacuum pipe is set as the boundary, which could be either elliptical or rectangular. The other method uses the integrated Green function (IGF) with open space boundary conditions [9]. The boundary is assumed to be at infinity, which is valid when the beam size is much smaller than the vacuum chamber. This can save computational resources and ensure computational speed while maintaining accuracy if the ratio of the pipe size to the beam size is large. Furthermore, the linear chromaticity, the thin nonlinear elements [24], and the bunch-by-bunch feedback system [26] are available in GOAT to explore possible mitigation measures for high-intensity effects.

Relying on the abundant infrastructure element module and flexible numerical model, the impedance induced collective instability, space charge effect, electron cloud effect, beam-beam interaction, and their interplay can be simulated in the GOAT code using the kick approximation. Together with the benchmark results against other frequently-used codes, their physical models and numerical approaches are introduced in the following sections.

Impedance

In GOAT, the simulation of impedance induced collective instability is simple. Only beam transportation and impedance elements are required to form the beamline. The transverse dipolar impedance, transverse quadrupolar impedance, and longitudinal impedance are available in the code. Instead of the impedance, the wake function, which is the equivalent expression of the impedance in the time domain, is used to calculate the wake kick experienced by the beam particles [4, 24]: $\begin{array}{c}\text{Δ}{p}_{u\mathrm{,}i}=-\frac{q^2}{{\beta }^2{E}_0}{\displaystyle \sum _{j=i+1}^{N_{\text{Slice}}}{N}_j[\langle u_j\rangle W_{\text{D}}(z_i-z_j)}\\ +u_i{W}_{\text{Q}}\left(z_i-z_j\right)]\mathrm{,}\\ \text{Δ}{p}_{z\mathrm{,}i}=-\frac{q^2}{{\beta }^2{E}_0}{\displaystyle \sum _{j=i+1}^{N_{\text{Slice}}}{N}_j{W}_0}\left(z_i-z_j\right),\end{array}$ (1) where u denotes x or y of the ith particle, $\langle u_j\rangle$ the centroid of all Nj particles in the jth slice, pz the momentum deviation, β the Lorentz velocity, and W_D, W_Q, and W₀ the transverse dipolar wakefield, transverse quadrupolar wakefield, and longitudinal wakefield, respectively.

For the broadband resonator (BBR) impedance model, a built-in method is used to calculate the wake function in the space domain. For other impedance models, coarse wake functions are obtained by reading external input files with two columns of data: one for the position and the other for the corresponding wake function. The linear interpolation method is adopted for finer wake function calculations. The slicing technique can be employed to achieve a better computational performance than calculating the wake force between two macroparticles [27]. The wake force between slices is calculated using the transverse and longitudinal centroids. All the macroparticles contained in each slice experienced the same wake force. If only the transverse impedance is included, the beam is transported linearly in the ring using a transfer matrix for both the transverse and longitudinal planes. Only an RF element is required for the longitudinal impedance study. It is noteworthy that the RF element models particle motion in the (Δ E, θ) phase space, whereas the (z, pz) phase space is used in the impedance element. The particle coordinates are transformed using the beam module. In both cases, the wake kick for the beam is integrated at a single interaction point [24, 27].

3.1

Transverse mode coupling instability

In transverse mode coupling instability (TMCI), the frequency shift of each azimuthal satellite increases with impedance or beam intensity. When two adjacent modes collide and merge, the beam becomes unstable, and the oscillation of the bunch center starts to grow exponentially [4]. This type of instability usually requires a detailed study because it is destructive and can cause beam loss. Based on the beam parameters in Table 1, the impedance threshold for pRing is studied. BBR is used in the simulation to approximate the transverse impedance model in pRing, which is a reasonable estimation when the pRing’s impedance model is not fully built. The following values are considered for other parameters: quality factor Q = 1 and resonant frequency fr = 1 GHz. The parameters are chosen mainly because the cutoff frequnecy of the vacuum chamber with an average radius is approximately 1 GHz.

The proton beam is initialized with transverse and longitudinal Gaussian distributions. Taking the vertical plane as an example, the impedance is scanned from 0 to 8 MΩ/m at intervals of 0.04 MΩ/m. Each set of simulations is performed for 2¹⁵ turns. Fig. 2a shows the spectrum obtained via fast Fouirer transformation (FFT) using the turn-by-turn vertical bunch centroid recorded in the simulation. As shown in the figure, with the increase of impedance, the frequency of 0 mode and -1 mode moves down and up, respectively. These two modes are coupled at an impedance of 7 MΩ/m, after which beam loss occurs rapidly and the spectrum is no longer accurate. As a benchmark, the same simulation is performed using PyHEADTAIL, and the spectrum is shown in Fig. 2b. As observed, the behaviors of the 0 and -1 modes are exactly the same as those predicted by GOAT, but the instability threshold of 7 MΩ/m is slightly higher. The prediction of the instability threshold differs by only 0.6% between the two codes. Comparing the two spectra, -2 mode and +1 mode are also clearly visible apart from the two strongest modes. Meanwhile, several lines appeared for each azimuthal satellite, which can be interpreted as different radial modes of the same azimuthal mode. From the perspective of the spectra, the results of the two codes are consistent.

Fig. 2Full-screen View

(Color online) Real part of the normalized mode-frequency shifts in the transverse plane given by (a) GOAT and (b) PyHEADTAIL via the macroparticle tracking with the BBR impedance model.

In addition, in Fig. 3, the instability growth rates obtained from the two codes are compared by fitting the envelope of the bunch centroid oscillation. The fitting window is in the interval of 20–80% of the total data selected before beam loss. Again, the TMCI growth rates predicted by the two codes matches perfectly. All benchmark results against PyHEADTAIL show that the beam transportation element and impedance element implemented in GOAT function very well.

Fig. 3Full-screen View

(Color online) Growth rates of the transverse mode coupling instability as a function of the impedance given by PyHEADTAIL (black dots) and GOAT (red dots).

3.2

Longitudinal microwave instability

In the longitudinal plane, the induced voltage generated by the longitudinal impedance is superimposed on the RF voltage, causing potential well distortion. As the impedance increases further, the microwave instability region is reached [28]. When the instability threshold is exceeded, the energy spread and bunch length start to grow drastically. Although the microwave instability can self-stabilize and will not lead to beam loss, the beam-beam performance in colliding beams is likely to be affected. Similar to the transverse plane, the BBR is used to estimate the longitudinal impedance model in pRing. The quality factor is Q = 1. For the resonant frequency, f_r = 2 GHz is chosen such that the rms bunch length is longer than the oscillation period of the BBR model, and makes the bunch particles see mostly the inductive part of the impedance. Thus, the extent of bunch lengthening can be investigated, as it is an important factor for peak luminosity.

The simulations are conducted using BLonD and GOAT. The proton beam parameters used in the simulations are listed in Table 1. The beam size, momentum spread, and emittance are simulated by increasing the longitudinal impedance from 0 to 80 kΩ with an interval of 5 kΩ. Each set of simulations is performed for 50000 turns. In order to avoid the filamentation caused by the induced voltage at the initial stage, the impedance is elevated adiabatically to maximum within the first 8000 turns, which is more than 100 times the synchrotron period. As shown in fig: Fig. 4, the instability threshold predicted by both codes is 35 kΩ. All the results plotted in the figures are obtained by averaging the bunch statistics after reaching a dynamic balance. In the potential well distortion region, the predictions of the longitudinal beam parameters from the two codes are identical. Note that the bunch length increased and the momentum spread shrank with an increase in impedance. This is because of the absence of synchrotron radiation in proton beams. Below the instability threshold, the evolution of the bunch length and momentum spread must ensure that the longitudinal emittance is conserved [29], as shown in subfig: Fig. 4c.

Fig. 4Full-screen View

(Color online) (a) Bunch length, (b) momentum spread, and (c) emittance as a function of the impedance obtained by GOAT (black dots) and BLonD (red dots).

However, as the strength of the microwave instability is enhanced with higher impedance, divergences gradually appeared in the simulation results. Taking the momentum spread as an example, the maximum relative error between different codes is approximately 9%. Similar to results reported in Ref. [30], bifurcations arise between the different simulation codes when the microwave instability become sufficiently strong. Nevertheless, both BLonD and GOAT provides the same instability threshold, and the tendencies of the beam parameters obtained through the different codes are not significantly different. The correctness of the GOAT code in modeling the longitudinal collective effect is verified. In addition, other information, such as the bunch distribution and detailed evolution of the beam parameters, can also be extracted from the simulations. Fig. 5a shows the longitudinal bunch shape in the potential well distortion region with an impedance value of 25 kΩ, and Fig. 5b shows the evolution of the bunch length and momentum spread, where the microwave instability has already occurred with an impedance value of 45 kΩ.

Fig. 5Full-screen View

(Color online) (a) Longitudinal bunch shape normalized to unit with the longitudinal impedance value of 25 kΩ in the PWD region. (b) Evolution of the bunch length and momentum spread with the impedance value of 45 kΩ in the microwave instability region.

Space charge effect

The simulation of the space charge effect is straightforward. However, the particle motion is continuous under the action of the space charge self-field. Applying an integrated space charge kick to the beam once per turn leads to incorrect results [27]. Therefore, the ring is divided into several segments using an external optical input file. These segments can be uniform or nonuniform. A space charge element is placed at each node. Within each space charge element, the particle coordinates are first Lorentz boosted from the co-moving laboratory frame to the beam frame, and the FD Poisson solver is used to compute the electrostatic field slice-by-slice in the beam frame. After considering the magnetic field when converted to the laboratory frame, the space charge kick is applied to each particle. The transverse tracking between two adjacent space charge elements is linear. When synchrotron motion is considered, longitudinal dynamics can be considered as either linear or performed by the RF element.

Considering the space charge self-field, the envelope equations can be written as [31] $\begin{array}{l}{r^{″}}_x+K_x\left(s\right)r_x-\frac{{ϵ}_{x\mathrm{,}\text{geo}}^2}{r_x^3}-\frac{K^{\text{SC}}}{2\left(r_x+r_y\right)}=\mathrm{0,}\\ {r^{″}}_y+K_y\left(s\right)r_y-\frac{{ϵ}_{y\mathrm{,}\text{geo}}^2}{r_y^3}-\frac{K^{\text{SC}}}{2\left(r_x+r_y\right)}=\mathrm{0,}\\ K_{\text{SC}}=\frac{q\lambda }{2\pi {ϵ}_0{\beta }^2{\gamma }^2{E}_0}\mathrm{,}\end{array}$ (2) where r denotes the rms beam widths, $K\left(s\right)$ the focusing function of lattice, ϵ_geo the geometry emittance, ϵ₀ the permittivity of vacuum, γ the Lorentz factor, and λ the beam line density. The coherent quadrupolar tune that describes the envelope oscillation can be derived analytically using the coupled envelope equations in a smoothly approximated ring [31]: $\begin{array}{c}{Q}_{\pm }=2\left(Q_x^2+Q_y^2\right)-\frac{K^{\text{SC}}}{\left(r_x+r_y\right)}[1+\frac{3}{4}\left(\frac{r_y}{r_x}+\frac{r_x}{r_y}\right)\\ \mp \frac{\sqrt{1+D^2}}{2}]\mathrm{,}\\ D=4\frac{K_y-K_x}{K^{\text{SC}}}\left(r_x+r_y\right)+\frac{3}{2}\left(\frac{r_y}{r_x}-\frac{r_x}{r_y}\right),\end{array}$ (3) where D denotes the coupling parameter, and Qx,y are the bare tunes in the horizontal and vertical, respectively.

GOAT code is employed to study the transverse coupling phenomenon induced by the space charge effect. As indicated in eq3, the envelope tunes depend on the bare transverse tunes of the ring. Thus, seven groups of simulations with different initial bare tunes are performed by tracking the beam for 512 turns under the action of the space charge self-field. The bare vertical tune changes from Qy=21.28 to Qy=21.34, whereas the horizontal tune is fixed at Qx=21.31. The beam is initialized as a coasting beam with a transverse KV distribution represented by 1×10⁶ macroparticles. The horizontal and vertical beam sizes are 20 and 10 mm, respectively. According to numerical convergence studies, 500 space charge elements are uniformly placed along the ring. The smooth approximation is used in the optics file. In the simulations, the rms beam size is calculated and recorded at each node. Perform FFT on the turn-by-turn beam size at one node, and two envelope coherent modes are obtained. In Fig. 6, the simulated coherent tunes are compared with those predicted by eq3. The variation in the vertical beam size owing to the vertical beta function scaling with the vertical tune is considered for both the simulation and analytical calculation. As shown in the figure, GOAT matches the theory well.

Fig. 6Full-screen View

(Color online) Envelop mode tunes as a function of the bare vertical tune of the machine for the coasting pRing beam obtained by theory (lines) and macroparticle tracking (dots).

In addition, the incoherent tune spread of the beam particles with different amplitudes is a direct dynamic result of the space charge self-field. It can be used as a crosscheck for space charged elements. Under the smooth approximation, all particles in the KV distribution have exactly the same incoherent tune shift when subjected to the space charge self-field, which can be expressed by the following formula [31]: $\text{Δ}{Q}_{x\mathrm{,}y}^{\text{KV}}=-\frac{K^{\text{SC}}{R}^2}{R_{x\mathrm{,}y}\left(R_x+Ry\right)Q_{x\mathrm{,}y}}\mathrm{,}$ (4) where Rx,y are the rms beam sizes of a KV distribution in horizontal and vertical plane, respectively. Here, a KV beam with a radius of 10 mm is used, and the incoherent tune shift calculated using eq4 is -0.04. Correspondingly, the maximum incoherent tune spread of a Gaussian beam can also be obtained by analytical calculations [31]: ${\left\{\text{Δ}{Q}_{x\mathrm{,}y}^{\text{Gaussian}}\right\}}_{\text{max}}=\frac{K^{SC}}{4\pi }{\displaystyle \oint \text{d}s}\frac{{\beta }_{x\mathrm{,}y}\left(s\right)}{{\sigma }_{x\mathrm{,}y}\left(s\right)\left[{\sigma }_x\left(s\right)+{\sigma }_y\left(s\right)\right]}\mathrm{,}$ (5) where $\beta \left(s\right)$ is the beta function, the σx,y the rms beam size of a Gaussian distribution in the horizontal and vertical plane, respectively. According to the theory in accordance with the rms equivalence principle, a Gaussian beam with an rms size of 5 mm is considered, and the maximum incoherent tune spread is -0.08, which is twice the incoherent tune shift in the KV beam [5, 32]. For comparison, the numerical tracking is performed using the beam parameters employed in the analytical evaluation. In both cases, the beams are initialized with 5×10⁵ macroparticles, and 100 space charge elements are uniformly distributed along the ring. The beams are tracked linearly in the transverse plane and the longitudinal motion is frozen. As shown in Fig. 7, the black and red dots represent the incoherent tune spreads of all particles in the Gaussian and KV beams, respectively. The incoherent tune spread of the Gaussian beam particles reaches the theoretically calculated value indicated by the blue pentagram. Although the tune spread of the KV beam given by the simulation is not strictly a point marked by the green pentagram, as predicted by theory, the dispersion is very small. This is consistent with the results reported in Ref. [5]. This may be caused by fluctuations in the particle distribution or numerical errors in the field solver. Overall, the space charge element implemented in GOAT is correct.

Fig. 7Full-screen View

Incoherent tune spread obtained by the macroparticle tracking with transversely a Gaussian distribution (black dots) and a KV distribution (red dots). The incoherent tune shift of the KV distribution is denoted by the green pentagram, and the maximum incoherent tune spread of the Gaussian distribution is marked by the blue pentagram. The bare tune is represented by the purple pentagram.

Electron cloud effect

In modeling electron cloud effects, the cloud build-up process and beam-cloud interaction are treated separately, or the so-called “weak-strong” model [6]. This is because the dynamic balance of the electron cloud density in the vacuum chamber can be reached in a single or very few turns, whereas the beam-cloud instability study requires tracking the beam for many turns, at least longer than the instability growth time. In GOAT, two methods are established for beam-cloud interactions: the linearized method [31] and the self-consistent tracking method [7]. The linearized method is based on cloud dipolar and quadrupolar forces obtained from dedicated simulations. In the self-consistent tracking method, both the beam particles and electrons are characterized by macroparticles. The electromagnetic interaction between the beam and cloud is modeled by a set of thin interactions along the ring, and the forces acting on each other are solved numerically using the FD Poisson solver implemented in the code. Currently, the linearized method is preferred for studying the beam-cloud instability because the computation overhead for the tracking method is expensive.

5.1

Build-up simulation

In the build-up simulation, the very low density cold electrons are uniformly generated in the chamber as the “seed” electron at the initialization stage. When the bunch passes, primary electrons are created due to the residual-gas ionization or beam loss. These electrons, along with the surviving electrons, are accelerated to some energy under the action of the beam field. Secondary electrons are produced as energetic electrons hit the chamber wall. Only elastic scattering and real secondary emission are included in the secondary electron emission model. The energy carried by the electron determines the secondary electron yield (SEY) of an elastically scattered electron [6]: ${\delta }_{\text{elastic}}\left(E\right)=R_0{\left(\frac{\sqrt{E}-\sqrt{E+E_0^{\text{es}}}}{\sqrt{E}+\sqrt{E+E_0^{\text{es}}}}\right)}^2\mathrm{,}$ (6) where R₀ corresponds to the elastic reflection probability for an electron in the limit of zero primary electron energy, $E_0^{\text{es}}$ the fitting parameter of experimentally measured energy spectrum of elastically scattered electrons (chosen as a specific value according to the vacuum chamber’s material), and E the energy of the incident electrons (calculated via the velocity of each electron in the simulation, in unit of eV), or a true secondary electron is produced with SEY [6]: ${\delta }_{\text{true}}\left(E\right)={\delta }_{\text{max}}\frac{s\frac{E}{E_{\text{max}}}}{s-1+{\left(\frac{E}{E_{\text{max}}}\right)}^2}\mathrm{,}$ (7) where E_max is the energy of the electrons at which the SEY is maximum, δ_max the maximum value of the SEY at normal incidence, and s the fitting parameter of measured SEY curve (s=1.35 is used widely since it provides the most reasonable fit to the experimental data.). The newly generated electrons follow a logarithmic normal distribution in energy and a cosine distribution in direction. This process is repeated successively according to the user-defined beam filling pattern. It is worth noting that, due to the multipacting effect, the number of electron macroparticles and the charge carried by a single electron change constantly during the simulation [16]. For considerations of computing performance and memory, as in PyECLOUD, both the number of macroparticles and the amount of charge of single macroparticles are dynamically managed. All alive electrons are mixed and redistributed in the 6D phase space when the above two quantities exceed the pre-defined threshold. After the regeneration process, the charge of each electron is at the same level.

Electron cloud build-up simulations are performed in the drift, dipole, and quadrupole regions using PyECLOUD and GOAT software. The numerical parameters in Table 2 are used. Fig. 8 shows the evolution of the electron line density within the vacuum chamber simulated by the two codes. The line density in the quadrupole region shown in the figure is halved for comparison. Both codes yield the same results. The time scale shown here is only 1/8 of the revolution period. Nevertheless, owing to the extremely short bunch spacing, the electron cloud reaches equilibrium after dozens of bunches pass through. As the cloud becomes saturated, the electron densities in the field-free and quadrupole regions are relatively close, at approximately 4×10¹⁰/m. Although the horizontal motion is frozen in the presence of a dipole magnetic field, the electron density also reaches 2.5×10¹⁰/m and tends to increase slowly. In addition, the oscillation amplitude of the electron density is smaller in the presence of an external magnetic field than in the field-free region because most electrons are trapped by the magnetic field lines.

Table 2Full-screen View

Numerical parameters used for electron cloud simulations.

Bunch intensity, Nb (1011)	1.25
Beam size, σx, y (mm)	3.0/2.3
Chamber shape	Elliptical
Horizontal chamber aperture, a (mm)	45
Vertical chamber aperture, b (mm)	35
Dipole field (T)	1.0
Quadrupole gradient (T/m)	20.0
Electron reflectivity R0	0.7
$E_0^{\text{es}}$ for elastic scatters (eV)	150
s for true secondaries	1.35
Energy of δmax, Emax (eV)	332
Maximum SEY, δmax	2.0

Show more

Fig. 8Full-screen View

(Color online) Comparisons of the evolution of the electron line density in the drift, dipole, and quadrupole regions. The results are obtained by PyECLOUD (solid lines) and GOAT (dashed lines). The electron line density in the quadrupole region is halved.

Figure 9 shows the electron spatial distribution exported from GOAT in the dipole region just before the bunch arrival. Similar to other studies, vertical stripes are formed in the horizontal direction in the presence of a dipolar magnetic field [6, 7]. In fig: Fig. 10, comparisons of the transverse phase space and corresponding histograms of the two codes are shown. In the horizontal direction (subfig: Fig. 10a), the electron distribution is figure 8 shaped. The central density is relatively low. Almost all the electrons gather in the range of two stripes. The velocity distribution of the cloud is Gaussian-like, and the electrons are nearly static. In the vertical direction (subfig: Fig. 10b), however, the electrons are distributed over the chamber since the motion is unconstrained. At the same time, it can be noted that the electrons are highly concentrated near the chamber wall, which is a significant feature of a saturated cloud. The electrons in the vicinity of the chamber wall form a potential barrier to prevent the energetic electrons from hitting the pipe and producing secondary electrons. The vertical velocity component is sharper than the horizontal one, and the value is approximately an order of magnitude higher. The comparison shown in the histograms further confirms the consistency between PyECLOUD and GOAT.

Fig. 9Full-screen View

(Color online) Snapshot of the electron distribution in the dipole region before the bunch arrival at the cloud sturation stage.

Fig. 10Full-screen View

(Color online) The horizontal (a) and vertical (b) phase space distributions and statistical histograms of the electron cloud in the presence of a dipole magnetic field. The histograms are normalized to unit.

5.2

Beam-cloud interaction

Two quantities are required due to the electron pinch as the bunch passes through in modeling the beam-cloud interaction using the linearized method: the generalized 2D dipolar wakefield [33-35] and the longitudinally varied betatron tune shift [32, 36]. The dipolar term is first dicussed. When an on-axis bunch passes through a cloud, the electrons are attracted toward the axis and the centroid of the cloud does not move. However, the situation changes if a small transverse offset is introduced to the bunch. When the displaced bunch passes through a cloud, the cloud is redistributed and begins to oscillate. Then, the subsequent portions of the bunches are deflected. Therefore, to obtain the dipolar wakefield, the bunch is sliced longitudinally, and a transverse offset is added to the driving slice. The dipolar wakefield generated by the driving slice can be computed for the subsequent testing slices via [35] $W_{j\mathrm{,}x\mathrm{,}y}\left(z_i\mathrm{,}{z}_j\right)=\frac{E_{j\mathrm{,}x\mathrm{,}y}}{e{\text{Δ}}_{i\mathrm{,}x\mathrm{,}y}{N}_i}L\mathrm{,}$ (8) where L denotes the length of the cloud, Ex, y the electric field averaged over the particle distribution, Δ the transverse offset attached to the driving slice containing N particles, and the subscript i and j are respectively the driving and the testing slices. Instead of depending on the relative distance between the ith driving slice and the jth testing slice only, the dipolar wakefield produced by the electron cloud depends on two longitudinal positions simultaneously as a result of the electron pinch. For example, in the field-free region, by changing the position of the driving slice along the slice sequence, the generalized dipolar wakefield can be obtained, as shown in Fig. 11. z > 0 corresponds to the bunch head. The electron distribution used in the simulation is extracted from the dedicated buildup simulation described in the previous section. The transverse offset is 10% of the vertical beam size. Longitudinally, a Gaussian distribution with nominal bunch length is used. The bunch is cut into 50 slices within the range of ($-4{\sigma }_z$, $4{\sigma }_z$), where σz is the rms bunch length. It is worth noting that for the purpose of reducing the numerical noise in the simulation, the electron distribution is created by merging ten realistic distributions at the same time step just before the bunch arrival. Moreover, a set of reference “driving slice” groups without attached initial transverse offset is established in each simulation, which has the same distribution as the beam and is used as the base of the numerical noise.

Fig. 11Full-screen View

(Color online) The generalized 2D dipolar wakefield generated by the electron cloud in the field-free region. z > 0 means the head of the bunch. The realistic distribution of electrons from dedicated build-up simulation are used.

As shown in Fig. 11, the consequences of the electron pinch are significant. The generalized dipolar wakefield does not satisfy the translational invariance. This is a basic property of the conventional electromagnetic impedance. The generalized wakefield peaks at the bunch head at a value of 8×10¹⁶ VC^-1m^-1. In contrast to the results reported in [35] with an initially static electron distribution, no significant increase in the wakefield amplitude caused by the electron pinch is observed when a realistic electron distribution is used. This illustrates that the enlargement of the wakefield amplitude is closely related to the electron velocity distribution.

To obtain the betatron tune shift variation along the bunch length, the self-consistent macroparticle tracking method is used. A symplectic expression has been derived for the beam-cloud interaction [37] $\begin{array}{l}{x}^{\text{new}}=x^{\text{old}}\mathrm{,}\\ p_x^{\text{new}}=p_x^{\text{old}}-\frac{qL}{{\beta }^2{E}_0}\frac{\partial \Phi }{\partial x}\left(x\mathrm{,}y\mathrm{,}z\right)\mathrm{,}\\ y^{\text{new}}=y^{\text{old}}\mathrm{,}\\ p_y^{\text{new}}=p_y^{\text{old}}-\frac{qL}{{\beta }^2{E}_0}\frac{\partial \Phi }{\partial y}\left(x\mathrm{,}y\mathrm{,}z\right)\mathrm{,}\\ z^{\text{new}}=z^{\text{old}}\mathrm{,}\\ p_z^{\text{new}}=p_z^{\text{old}}-\frac{qL}{{\beta }^2{E}_0}\frac{\partial \Phi }{\partial z}\left(x\mathrm{,}y\mathrm{,}z\right)\mathrm{,}\end{array}$ (9) where Φ is the electrostatic potential of the cloud. Instead of the typical approaches of including only the transverse part of this map, the complete map is implemeted in GOAT. Based on the numerical convergence studies, 20 thin electron cloud elements are uniformly placed around the ring. The beam is represented by 5×10⁵ macroparticles with 50 longitudinal bunch slices. Fresh cloud distributions with approximately 2×10⁶ macroparticles at each node are generated and saved at the initialization stage. Again, a realistic electron distribution is used. Between electron cloud elements, the beam is linearly transported in the transverse plane. And the longitudinal motion is frozen to obtain the instantaneous tune spread.

In Fig. 12, the blue dots represent the incoherent tune spread caused by the electron cloud force. There are several steps in the tune spread distribution, particularly at the bunch head with z > 0. This is mainly because the discrete equation of motion is used to integrate the electron motion, and the change in the spatial morphology of the cloud is discontinuous. The tune shift obtained by averaging over the tune spread of particles within the slices at different longitudinal positions is plotted as the red line. The electron pinch effect results in a change in the focusing forces during the bunch passage, and therefore, in the betatron tune modulation with longitudinal coordinates. One can observe that the tune shift increases by a factor of eight from the head to the center of the bunch and then slowly decays toward the tail. In addition, a line parallel to the horizontal axis at Δ Qy=0.18 is clearly visible. This is because some particles at the bunch center experiencing very strong cloud-gradient forces and crossing the half-integer resonance Qy+0.18=22.5, as indicated by the dashed green line.

Fig. 12Full-screen View

(Color online) Incoherent tune spread (blue dots) of the bunch under the action of the electron cloud force in the field-free region. The varied transversal tune shift in y direction (red line) along the longitudinal position can be obtained by averaging the incoherent tune spread of particles in the slices at each longitudinal position. The half-integer resonance line Qy+0.18=22.5 is shown in green dashed line.

Taking the above results as a wakefield file for the external input, the impedance element can be used to simulate the beam-cloud instability. The cloud-generated dipolar and quadrupolar forces can be modeled by transverse driving and detuning impedance implemented in the impedance element. The impedance element is well benchmarked in sec3. The nonlinear characteristics of the beam-cloud interaction are omitted in such a linearized method. However, this method is correct and has advantages obviously in fast and enormous parameter scans [32].

Beam-beam interaction

For colliding beams, the GOAT code describes the electromagnetic interaction of two counter-rotating beams through the 6D symplectic Synchro-Betatron Mapping method [38]: $\begin{array}{l}{x}^{\text{new}}=x^{\text{old}}-\frac{eS}{{\beta }^2{E}_0}{\tilde{E}}_x^{\text{CP}}\mathrm{,}\\ p_x^{\text{new}}=p_x^{\text{old}}+\frac{e}{{\beta }^2{E}_0}{\tilde{E}}_x^{\text{CP}}\mathrm{,}\\ y^{\text{new}}=y^{\text{old}}-\frac{eS}{{\beta }^2{E}_0}{\tilde{E}}_y^{\text{CP}}\mathrm{,}\\ p_y^{\text{new}}=p_y^{\text{old}}+\frac{e}{{\beta }^2{E}_0}{\tilde{E}}_y^{\text{CP}}\mathrm{,}\\ z^{\text{new}}=z^{\text{old}}\mathrm{,}\\ p_z^{\text{new}}=p_z^{\text{old}}+\frac{e}{2E_0}[{\tilde{E}}_x^{\text{CP}}\left(p_x^{\text{old}}+\frac{e}{2{\beta }^2{E}_0}{\tilde{E}}_x^{\text{CP}}\right)\\ +{\tilde{E}}_y^{\text{CP}}\left(p_y^{\text{old}}+\frac{e}{2{\beta }^2{E}_0}{\tilde{E}}_y^{\text{CP}}\right)+{\tilde{E}}_z^{\text{CP}}]\mathrm{,}\end{array}$ (10) where the S denotes the drift distance, ${\tilde{E}}_x^{\text{CP}}$, ${\tilde{E}}_y^{\text{CP}}$, and ${\tilde{E}}_z^{\text{CP}}$ the electric field generated by the opposite bunch slice at the collision point. The linear transfer matrix is used for the motion of beam particles in both transverse and longitudinal planes between different interaction points.

The calculation sequence of the beam-beam interaction is as follows. The two bunches are first sliced longitudinally using the uniform charge slicing method implemented in the code, and then the bunches collide slice-by-slice. At each collision step, the particles contained in the slice are drifted from the interaction point (IP) to the collision point (CP) located at $S=\left(z_1+z_2\right)/2$, where z₁ and z₂ are the longitudinal centroids of the two slices. The IGF Poisson solver is used to compute the electric fields generated by the slices according to the current distributions. Subsequently, taking into account the magnetic field and time-of-flight effects, the kick is applied to the particles in the opposite slice. Finally, the particles are transferred back to the IP via the inverse drift operator. In general, the linear interpolation method [39] is applied to calculate the kick for each particle in terms of computational convergence and speed.

In addition, the requirements for the fast separation of the two colliding beams, the overall detector component and interaction region (IR) magnet arrangements strongly depend on a large crossing angle shown in Fig. 13. To include the crossing angle in the beam-beam simulation procedure given by eq10, which is derived without the crossing angle, the Lorentz boost [40]: $\begin{array}{l}{\tilde{p}}_x=\frac{p_x-h\tan {\theta }_{\text{c}}}{\cos {\theta }_{\text{c}}}\mathrm{,}\\ {\tilde{p}}_y=\frac{p_y}{\cos {\theta }_{\text{c}}}\mathrm{,}\\ {\tilde{p}}_z{}_z=p_z-p_x\tan {\theta }_{\text{c}}+h{{\displaystyle \tan }}^2{\theta }_{\text{c}}\mathrm{,}\\ \tilde{x}=x\left[1+\frac{{\tilde{p}}_x}{{\tilde{p}}_s}\sin {\theta }_{\text{c}}\right]+z\tan {\theta }_{\text{c}}\mathrm{,}\\ \tilde{y}=y+\frac{{\tilde{p}}_y}{{\tilde{p}}_s}x\sin {\theta }_{\text{c}}\mathrm{,}\\ \tilde{z}=\frac{z}{\cos {\theta }_{\text{c}}}-\frac{\tilde{h}}{{\tilde{p}}_s}x\sin {\theta }_{\text{c}}\mathrm{,}\end{array}$ (11) can be used to transform the beam particles from laboratory frame to head-on frame. A variable marked with and without a superscript tilde indicates the quantity in the head-on frame and laboratory frame, respectively. The variable θ_c denotes the half crossing angle, and the h hamiltonian, $\begin{array}{l}h=1+p_z-\sqrt{{\left(1+p_z\right)}^2-p_x^2-p_y^2},\\ \tilde{h}=\frac{h}{{{\displaystyle \cos }}^2{\theta }_c}\mathrm{,}\\ {\tilde{p}}_s=\sqrt{{\left(1+{\tilde{p}}_z\right)}^2-{\tilde{p}}_x^2-{\tilde{p}}_y^2}.\end{array}$ (12) After collision, the beam particles need to be transformed back to the laboratory frame via inverse Lorentz boost, $\begin{array}{l}{p}_x={\tilde{p}}_x\cos {\theta }_{\text{c}}+h\tan {\theta }_{\text{c}}\mathrm{,}\\ p_y={\tilde{p}}_y\cos {\theta }_{\text{c}}\mathrm{,}\\ p_z={\tilde{p}}_z+p_x\tan {\theta }_{\text{c}}-h{{\displaystyle \tan }}^2{\theta }_{\text{c}}\mathrm{,}\\ x=\left(\tilde{x}-\tilde{z}\sin {\theta }_{\text{c}}\right)/\left(1+\frac{{\tilde{p}}_x}{{\tilde{p}}_s}\sin {\theta }_{\text{c}}+\frac{\tilde{h}}{{\tilde{p}}_s}{{\displaystyle \sin }}^2{\theta }_{\text{c}}\right),\\ y=\tilde{y}-\frac{{\tilde{p}}_y}{{\tilde{p}}_s}x\sin {\theta }_{\text{c}}\mathrm{,}\\ z=\tilde{z}\cos {\theta }_{\text{c}}+\frac{\tilde{h}}{{\tilde{p}}_s}x\sin {\theta }_{\text{c}}\cos {\theta }_{\text{c}}\mathrm{.}\end{array}$ (13) Currently, only the ‘strong-strong’ model is available in GOAT.

Fig. 13Full-screen View

(Color online) Beams colliding at an angle in the original laboratory frame (top) and in the boosted head-on frame (bottom).

Predicting beam parameters such as luminosity, beam size, and centroid oscillation, studying the dependencies between different parameters, and investigating the role of beam-beam dynamics on the cross-talk of multiple beam dynamics are the main tasks of the beam-beam simulation. In order to check the validity of the beam-beam model implemented in GOAT, a set of beam-beam simulations is carried out with and without considering the crossing angle utilizing Athena and GOAT. The beam parameters in Table 1 are used. The luminosity and beam size are shown in Fig. 14a and fig: Fig. 15. When the two beams collide in the head-on frame, the transverse beam sizes of the proton beam are stable at the design value due to its small beam-beam parameter. For the electron beam, the beam sizes in both the horizontal and vertical directions shrink to slightly smaller values than the design value, which is also why the luminosity is higher than the calculated value. For the collision with a crossing angle of 50 mrad, although the luminosity is reduced by more than six times compared to the head-on collision, the beams are stable, and no significant beam sizes blow up are observed. The transverse beam sizes of protons and electrons are stable around the design values after reaching equilibrium. This indicates that the luminosity degradation is caused by geometric loss. Again, the simulation results given by Athena and GOAT agree very well with each other.

Fig. 14Full-screen View

(Color online) (a) The luminosity evolution in case of collision with and without considering the crossing angle obtained by Athena and GOAT. (b) The impact of the crab cavity frequency on the luminosity obtained by GOAT.

Fig. 15Full-screen View

(Color online) The corresponding horizontal (a) and vertical (b) beam sizes of colliding beams for the cases of Fig. 13(a).

In the EicC conceptual design, the luminosity reduction caused by the crossing angle is compensated by the crab cavity. The transfer map of the crab cavity in the x-z plane for each particle is given by [41], $\begin{array}{l}{p}_x^{\text{new}}=p_x^{\text{old}}+\frac{qV}{{\beta }^2{E}_0}\sin \left(\frac{\omega z^{\text{old}}}{c}+\varphi \right),\\ p_z^{\text{new}}=p_z^{\text{old}}+\frac{qV\omega }{{\beta }^{2}c{E}_0}{x}^{\text{old}}\cos \left(\frac{\omega z^{\text{old}}}{c}+\varphi \right),\end{array}$ (14) where V is the voltage of the crab cavity, ω the angular frequency of the cavity, and ϕ the phase of the cavity. The cavity voltage is chosen as: $V=\frac{c{E}_0\tan \left(\frac{{\theta }_c}{2}\right)}{q\omega \sqrt{{\beta }^{*}{\beta }_{\text{CC}}}\sin \text{Δ}\psi }\mathrm{,}$ (15) where β^* is the beta function at the IP, β_CC the beta function at cavity location, and Δ ψ the phase advance between the IP and the cavity. For the proton beam, the higher energy and smaller beta function at the IP lead to a much higher cavity voltage requirement than that for the electron beam. Generally, higher frequencies are preferred for reducing the cavity voltage. Fig. 14b shows the luminosity evolution obtained by scanning the crab cavity frequencies at a fixed half-crossing angle of 25 mrad. The luminosity degradation is low since the cavity frequency is below 200 MHz. As the frequency continuously increases, the luminosity degradation becomes higher. This phenomenon can be explained by the x-z distribution at the IP for different crab cavity frequencies, as shown in Fig. 16. Because of the incomplete deflection caused by the crab with a higher frequency (lower wavelength), the head and tail of the bunch are farther off from the axis at the IP, resulting in an increase in the overlap area of the two bunches and a decrease in luminosity. In addition, the simulated luminosity over 20000 turns is higher than the design value, even for a cavity frequency of 400 MHz. However, the synchrotron-betatron (SB) resonance may still be excited owing to the longitudinally position-dependent beam-beam kick, causing a reduction in the luminosity lifetime [42].

Fig. 16Full-screen View

(Color online) The x-z distributions of the bunch at IP with different crab cavity frequencies. The axis are normalized to the nominal horizontal beam size and bunch length. The crab frequency is 100 MHz (top left), 200 MHz (top right), 300 MHz (bottom left), and 400 MHz (bottom right).

Applicatiion of cross-talk

The physical elements implemented in GOAT are well benchmarked in the previous subsections. Thus, its validity is guaranteed. However, the beam in a real machine cannot be affected by a single effect, and there is cross-talk between different effects. Among the many high-intensity effects in a collider, the beam-beam interaction is one of the most important beam dynamics processes, which has a direct impact on the luminosity and integrated luminosity of the machine. Several instabilities have been observed in previous studies due to the cross-talk between the beam-beam interaction and other intensity-dependent effects [11-13, 43, 44]. It is necessary to explore the mechanisms of these instabilities and corresponding mitigation measures at the machine design stage.

In this subsection, a numerical example is presented for the cross-talk between the beam-beam interaction and the pRing’s vertical impedance in EicC. To emphasize the importance of cross-talk simulation, three sets of simulations are performed. The first purely includes the effect of the pRing’s vertical impedance, the second considers only the beam-beam interaction, and the third considers the self-consistent treatment of the two effects. The first two sets of simulations are performed in previous subsections, and the results of the evolution of beam parameters are presented here. In the self-consistent simulation, the transverse impedance element, beam-beam interaction element, and transverse and longitudinal linear transportation elements are combined to form a ring. A vertical impedance value of 6 MΩ/m is chosen, which is below the instability threshold. And the beam-beam interaction is assumed to take place in the head-on frame. Fig. 17a and Fig. 17b show the evolution of the beam vertical centroid and emittance for the three cases. It can be observed that below the instability thresholds of TMCI and beam-beam interaction considered alone, a dipole instability arises in the simulation. The vertical centroid grows exponentially, and the emittance blows up. The simulation results suggest that there is no coupling between the coherent beam-beam mode and the longitudinal sideband because the bare tunes of the two beams are different. This is very different from the mode coupling instability for symmetric collisions reported in Ref. [12]. It is also worth noting that such a dipole instability does not appear in the horizontal plane. This implies that the obseved head-tail type instability is closely related to the hourglass effect, because it is the only major difference between the two transverse planes. Further studies are still required to identify and understand the underlying mechanisms of this coherent instability.

Fig. 17Full-screen View

(Color online) Evolution of vertical bunch centroid (a) and emittance (b). Three cases are considered, first is purely beam-beam interaction, second is purely transverse impedance, and third is cross-talk of beam-beam interaction and transverse impedance.

Obviously, this numerical application indicates the necessity of simulating multiple physical processes in a self-consistent manner. This also demonstrates the flexibility, completeness, and advancement of GOAT.

In addition, the linear chromaticity and the ideal bunch-by-bunch feedback system implemented in GOAT are utilized to suppress this coherent instability. In Fig. 18, the instability growth rates for different chromaticity values are presented in the range from -10 to 10. The growth rates increase and then remain almost constant when the chromaticity values are negative. For positive chromaticity values, non-monotonic growth rate behavior is observed. When the charomaticity value is greater than eight units, the instability is fully suppressed. Then, the ideal bunch-by-bunch feedback system is employed to damp the instability. The evolution of the vertical centroid and emittance for different damping rates are presented in Fig. 19a and Fig. 19b, respectively. Compared to chromaticity, the feedback system is more effective in suppressing this instability. It can be seen that the dipole motion is suppressed, and the emittance is conserved, even though the damping rate is very small. This can be explained by Fig. 20, which illustrates the intra-bunch motion of the proton beam extracted from the simulation in successive turns. The most unstable mode is the 0-mode and the bunch-by-bunch feedback can eliminate the oscillations of the bunch as a whole.

Fig. 18Full-screen View

The impact of chromaticity on the growth rate of coherent instability.

Fig. 19Full-screen View

(Color online) The impact of ideal bunch-by-bunch feedback system on the coherent instability. The evolution of vertical bunch centroid (a) and emittance (b) with different damping rates.

Fig. 20Full-screen View

Intra-bunch motion for successive 100 turns of proton beam extracted from the cross-talk simulation.

Summary and outlook

In this paper, a simulation code, GOAT, is developed for single-bunch high-intensity beam dynamics. It can be used to simulate all intensity-dependent effects in the pRing of the EicC project. The code architecture and numerical model are introduced. Four simulation examples, including impedance induced collective instability, space charge effect, electron cloud effect, and beam-beam interaction, are conducted based on abundant elements and the flexible numerical models provided by GOAT. The results are well benchmarked with other codes and theories. In addition to separate simulations, an application is presented for cross-talk between the beam coupling impedance and the beam-beam interaction. The comprehensiveness and correctness of GOAT are verified. The Python program coded by the OPP technology ensures the scalability of GOAT as well as the independence of the modules. Different effects can easily be integrated into the code. GOAT can also be used to simulate the intensity-dependent beam dynamics in other accelerators and colliders. In the future, it is planned to improve the performance of the GOAT code, including algorithm optimization and hardware support (such as parallelization based on GPU (Graphic Processing Unit)). With the help of parallelization techniques, GOAT can be upgraded to include the multi-bunch beam dynamics.