Simulations of Beam Envelope Dynamics in Circular Accelerators

The envelope instability of high intensity beams in circular accelerators is studied via multiparticle simulations. The space charge kicks are derived from a Gaussian potential model for an efficient tracking. The evolution of the envelope phase space coordinates are derived from the bunch distribution. We found that the envelope stop band played an important role in emittance growth. Correction schemes of the envelope stop bands are studied. Because the space charge force pushes the envelope tunes downward, harmonics less than twice the betatron tunes are also important on emittance growth. Our code is efficient and fast, it can be used to study the effect of space charge force on high power accelerators.


I. INTRODUCTION
High-intensity accelerators can tolerate little beam loss, which may be related to beam emittance growth. Emittance evolution for space charge dominated beams has been one of many important topics in high-power beam accelerators and storage rings. Possible sources of emittance growth in LINACs are synchrobetatron coupling, halo formation induced by collective envelope modes or structure resonances [1], envelope resonance, and systematic fourth order resonance [2].
Important candidates for the emittance growth in circular accelerators include the half-integer stop band that can perturb the beam envelope function, the Montague resonance, and other systematic space charge resonances [3][4][5]. Many numerical simulation codes, such as TRACE [6], PARMILA [7], ORBIT [8], WARP [9], etc., were developed to study the mechanism of emittance growth. COSY INFINITY has also been proposed to expand its code to deal with space charge force [10,11]. A few systematic emittance evolution experiments were done in LINACs [12], and fewer in circular accelerators [4,13]. Systematic measurements and analysis of the emittance evolution with different injection turns at the Fermilab Booster, see e.g., Fig. 10 in Ref. [4], show that the emittance growth arises from linear coupling resonances from random skew quadrupoles.
Cousineau and collaborators have carried out systematic emittance measurement of high-intensity beams at the Proton Storage Ring (PSR) in Los Alamos National Laboratory and made systematic benchmarking of a space charge simulation code for implying the importance of a half-integer stop band in the emittance growth at the PSR [14]. Yet, attempts to make half-integer stop band correction at the PSR and Fermilab Booster were not as successful.
Charged particles in accelerators encounter intrinsic Coulomb force with each other, called the space charge force. It is classified into intrabeam scattering and the "mean field" due to collective effect of all other particle. The typical time constant for intrabeam scattering is minutes or hours, an important slow process on beam quality for storage rings. On the other hand, the effect mean field force can cause emittance growth due to resonances. The emittance growth is fast and can be controlled. We study the mean field effect of Coulomb force on beam dynamics.
There were seminal works on space charge dominated beam by Kapchinskij, Vladimirskij, Lapostolle, and Sacherer [15][16][17]. There were several simulations of highly space charge dominated beam for periodic focusing paraxial beams in LINACs [18,19]. There is no systematic numerical simulation on the space charge envelope dynamics. This paper is intended to carry out systematic study on the envelope dynamics of space charge dominated beam in synchrotrons. We organize this paper as follows. Section II reviews space charge theories and discusses a multiparticle tracking technique. Envelope dynamics with simulation results are shown in Sec. III. The effects of envelope resonance stop band harmonics and the effect of emittance growth will be addressed in Sec. IV. Methods for possible emittance growth control are discussed in Sec. V. Conclusions of this study is given in Sec. VI.

II. SPACE CHARGE FORCE AND TRACKING METHODS
The electric scalar potential, a solution of the Poisson equation, of the Coulomb force for a tri-Gaussian beam moving at the velocity v along the s-direction in Frenet-Serret ðx; s; zÞ coordinates is [20] ϕðx; z; sÞ ¼ where N B e is the bunch charge; γ is the relativistic Lorentz factor; ϵ 0 is the permittivity of free space; σ x , σ s , and σ z are the rms bunch lengths; and λ is a dummy variable. In addition to the scalar potential, the vector potential component A s responsible for the magnetic field of the moving bunch is given by A s ¼ γvϕ=c 2 , which is proportional to the scalar potential. The synchrotron motion is usually relatively slow compared to the transverse motion. The longitudinal space charge can induce bunch lengthening, and result in electron cloud effects [21,22]. Passive or active longitudinal focusing has been implemented to alleviate this effect [23,24].
For fast ramping synchrotrons, the longitudinal space charge effect is normally small due to substantial rf voltage in the acceleration process. Thus the bunch length change is small [4]. The effect of transverse space charge effect can be approximated by a two-dimensional potential, characterized by the beam sizes σ x and σ z which also depend on s. Combining both the scalar and vector space charge potentials, we find the space charge potential in the Frenet-Serret coordinate system as [20] The space charge potential is proportional to a dimensionless quantity called space charge perveance K sc ≡ 2Nr 0 ffiffiffiffiffi ffi 2π p σ s is the charge per unit length, β and γ are the relativistic Lorentz factors, and r 0 ≡ e 2 4πϵ 0 mc 2 is the particle classical radius. The classical radii are 2.82 × 10 −15 m for electrons and 1.535 × 10 −18 m for protons. The space charge effect is especially severe in low energy proton or ion beams because their β and γ are small.

A. Particle tracking technique
The Vlasov-Maxwell equations can self-consistently describe the evolution of particle distributions under the effect of the space charge force. Unfortunately there are no analytic solutions. A numerical simulation with Monte-Carlo analysis is necessary. A self-consistent calculation requires evaluation of mean field of beam distributions at each stage of numerical simulation [25]. This is very time consuming, and impractical in the study of essential emittance growth mechanisms. There are a few techniques and hypotheses used in beam space charge simulations. For example, TRACE [6] uses linear space charge force, and PARMILA [7] uses the particle-in-cell (PIC) method to track representatives called macroparticles in a mesh. Usually, thousands of macroparticles are used to extract the statistical properties to study the space charge beam dynamics. We will use the Gaussian distribution to evaluate space charge kick. From previous studies (see e.g., [16,17,26]), our study of the envelope dynamics based on multiparticle simulation should provide essential physics on the emittance growth associated with envelope dynamics.
The space charge force is given by the partial derivative of the 2D potential in Eq. (1). The space charge kicks become [27] The complex error function erfcð·Þ has been tabulated and fast to compute. It is more computational efficient than the calculation of the improper integral of the space charge potential in Eq. (1). For example, tracking of 5000 particles for 1000 turns and 100 space charge kickers per turn can be done in 40 minutes in an Intel Xeon 3.1 GHz CPU. There is no singularity at σ x ¼ σ z in Eq. (2), thus we can use Taylor series expansion and keep the terms up to the required accuracy for σ x ≈ σ z [27].
We integrate the space charge tracking code into a symplectic tracking code SIMTRACK [28]. The standard format of the SIMTRACK can easily be implemented into any lattice design files. Although our space charge potential of Eq. (2) is based on an initial Gaussian distribution, envelope dynamics of other initial distribution function may also be studied in this integrated tracking code. The beam distribution may change during the evolution in the accelerator, the Gaussian space charge force approximation should provide essential physics of resonance effects during the accumulation and acceleration. A more general approach proposed by Lapostolle [29] is to approximate any distribution in Hermite series expansion. Since we are studying the envelope dynamics in this paper, our calculation is an efficient approximation for extracting essential physics of space charge in beams.

B. Incoherent space charge tune shift
The space charge force is defocusing and the space charge tune shift is always negative in both transverse directions. The space charge self-force depends on the particle's transverse position in the beam, each particle will have its own tune shift. For a Gaussian beam, the incoherent tune shift is [30] where Z 1 ðxÞ ¼ e −x ½I 0 ðxÞ − I 1 ðxÞ, Z 0 ðxÞ ¼ e −x I 0 ðxÞ, R is the average radius of the synchrotron, I 0 and I 1 are the modified Bessel functions, r ¼ σ z =σ x is the beam-size ratio, and u is a dimensionless dummy variable. The tune shift dependent on particle's action is plotted in Fig. 1(a) for various beam emittance ratio with the bunch line density N ¼ 1 × 10 8 m −1 , ϵ x ¼ 1 μm, and β is the average beta function.
The maximum space charge tune shift, called Laslett or the linear tune shift parameter, corresponds to the particle at the center of the bunch with actions J x ¼ 0 and J z ¼ 0. The maximum tune shifts are The linear tune shift parameter is normally denoted by ξ sc . It is a measure of how strong the space charge force could be. It is proportional to the space charge perveance K sc . It is also a figure of merit related to the beam envelope oscillation.
Assuming a round beam with equal bare tunes in both directions, the incoherent tune spread distribution of a bi-Gaussian beam, scaled with the Laslett tune shift parameter, is shown in Fig. 1(b) [31]. The average is at 0.633Δν sc with an rms spread 0.168Δν sc .

III. ENVELOPE DYNAMICS
The envelope dynamics are governed the equation first derived by Kapchinskij and Vladimirskij (KV) in 1959 [15]. This is self-consistent theory if the beam is in KV distribution. A generalization was made by Lapostolle [16] and Sacherer [17] by extending the KV beam to beam with elliptical symmetry. In their generalization, the corresponding KV beam envelopes are replaced by the beam's second momentsx ≡ ffiffiffiffiffiffiffiffi hx 2 i p andz ≡ ffiffiffiffiffiffiffiffi hz 2 i p , and the emittance is replaced by the rms emittance. The mathematical form of the generalized rms envelope equation is the same as that of the KV equation.
One can study the envelope dynamics by solving the envelope equations. Another treatment is by constructing the envelope Hamiltonian [3,32]: Both of the above methods are perturbative because the particle motion is also perturbed by the envelope of the beam, and the beam second moments depend on the motion of all particles.

A. Envelope oscillation
Instead of solving the KV envelope equation numerically or using the Hamiltonian approach, we use multiparticle tracking to study the envelope dynamics. Through multiparticle tracking, the evolution of beam envelope radii can be dynamically evaluated. The beam envelope phase space coordinates are defined statistically by , ε x = 1 μm, β = 13 m where y denotes both transverse directions x and z. The envelope phase space coordinates can be obtained from macroparticle simulation. Consider a high power accelerator made of 18 FODO cells with circumference C ¼ 324 m, i.e., it has 18 superperiods. Systematic space charge resonances are located at betatron tunes 4.5 and 3 for the fourth and the sixth order resonances respectively. The bare betatron tunes are chosen to be 2.60 horizontally and 2.95 vertically. The typical quadrupole strength is To provide realistic simulation, the physical aperture is set to AE2.5 cm in both transverse directions. Particles falling out of the aperture are considered lost and will not be included in the statistics of the rms envelope calculation. The turn-by-turn space charge perveance is adjusted by the particle loss during the tracking and the injection filling rate to simulate multiturn accumulation. After 30-turn beam accumulation to the desired density, it is then tracked for another 1024 turns for FFT analysis.
The envelope oscillation tune can be extracted by the FFT analysis. The FFT uses the periodic data from the middle of every focusing quadrupole. There are totally 18 equally spaced focusing quadrupoles. The sampling rate is 18f 0 , therefore the Nyquist frequency is 9f 0 , where f 0 is the revolution frequency. The FFT data can extract envelope tune up to 9. An example of the FFT spectrum of the horizontal envelope is shown in Fig. 2(a) for a beam with line density N ¼ 1 × 10 8 particles=m. The spectrum peaks at an envelope tune of around 5, which is the horizontal envelope tune at 2ν x − ξ sc;x . A smaller peak at 5.7 comes from the weak coupling from the vertical envelope oscillation. At the same time, a low frequency peak shows the difference of horizontal and vertical envelope tunes.
The envelope phase space ellipses with different line densities are shown in Fig. 2(b), where the top left plot shows the phase space ellipse with no space charge force. Without space charge force, the envelope oscillation tune is exactly twice the betatron tune. Because the initial bunch distribution does not match the beam emittance, the envelope phase space is an ellipse around the matched value of β x ϵ x .
As the space charge force increases, the equilibrium envelope radius increases because the matched betatron amplitude function increases. The envelope phase space ellipse is smeared and the envelope tune spread increases with space charge force. The smear effect arises from envelope space charge coupling, which is small but visible resulting an envelope tune spread. The α function and the center of the matched σ 0 are not zero, as shown in the bottom plots of Fig. 2(b), because the space charge kick is before the calculation of these covariance quantities.
Collecting top ten peaks in the envelope spectrum with different line densities gives Fig. 3(a), where the envelope tune shift and broadening of peaks are shown, and the predictions are given by 2ν x − Δν sc;x and 2ν z − Δν sc;z of Eq. (4). The envelope tune can cross the integer stop band ν env ¼ 5 without particle loss and no emittance growth because there is no driving source for the envelope resonance. The peaks of the envelope tune spread depends on the beam distribution and the envelope tune spread, independent of the distribution, increases linearly with the space charge parameter.

B. Envelope resonance
When quadrupole errors were added as the driving source to excite the envelope resonance, the envelope resonances appear. We consider the error ONE single quadrupole with integrated strength ΔK 1 L ¼ 0.01 m −1 , located at about 1.5 m from a defocusing quadrupole, where β x ¼ 17.5 m and β z ¼ 21.5 m. This error causes beta-beat, e.g., the perturbed beta functions at the error quadrupole become β x ¼ 10.0 m and β z ¼ 17.4 m. More importantly, the FFT spectrum of σ x (similarly for σ z ) exhibits integer harmonics. Figure 2(c) shows the FFT spectrum of σ x for the line density of N ¼ 10 8 m −1 .
The envelope phase space for beams with different line densities are shown in Fig. 2(d). The peaks of the FFT spectrum of rms envelopes with different line densities are shown in Fig. 3(b). The trend agrees with the study in PSR [33]. The horizontal envelope tune approaches an integer 5 when the line density reaches 7.5 × 10 7 m −1 . When the envelope tune falls onto an envelope resonance stop band, the resonance causes particle loss and emittance growth. Particle loss takes place during the injection-accumulation process.
Particle loss and emittance growth caused by the envelope stop band provide a limit of the Laslett tune shift. No matter how high the injected beam intensity is, the resulting beam intensity is limited by the space charge. Particle loss will reduce the beam intensity so that the envelope tune is bounded by the stop band.

A. Stop band harmonics
Since the envelope coupling is small as shown in Fig. 2(a), we can examine the envelope dynamics for a paraxial beam. Applying the Floquet transformation to the paraxial envelope equation, we find where R ≡ R b ffiffiffiffiffiffiffiffi βðϕÞϵ p is the normalized envelope radius; R b is the rms beam envelope; ϵ is the rms emittance, and ϕ ≡ 1 ν R s 0 ds βðsÞ . The new time coordinate ϕ ranges from 0 to 2π and the overdot notation stands for the derivative with respect to ϕ. For a beam with small space charge force, R ¼ 1 is the equilibrium radius. When space charge force becomes important, we will find R > 1.
The last term in Eq. (7) has the exact superperiodicity of the accelerator lattice. Expanding the space charge term in Fourier series, we find where ξ sc is the linear space charge parameter: It is proportional to the space charge perveance and the circumference C, and inverse proportional to the emittance. The ξ sc q n and χ n are the amplitude and the phase of the nth harmonic with

SIMULATIONS OF BEAM ENVELOPE DYNAMICS IN …
Phys. Rev. ST Accel. Beams 18, 024202 (2015) The ϕ-independent term increases the matched envelope radius by ξ sc 2ν , i.e., the betatron amplitude function is increased by a factor of 1 þ ξ sc =ν. Substituting R ¼ 1 þ ξ sc =ð2νÞ þ r into (7), we havë The envelope tune is ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4ν 2 − 4νξ sc p and the harmonics in the right-hand side are the driving terms. When ξ sc ≪ ν, the envelope tune is 2ν − ξ sc . For an accelerator without quadrupole field error, the envelope harmonics n is an integer multiple of the accelerator superperiod. These envelope harmonics are called systematic resonances. Accelerators are normally designed to avoid these systematic stop bands.
In the presence of quadrupole error ΔkðsÞ, the beta-beat and betatron phase perturbation will appear. The additional stop bands are called random stop bands. The beta-beat due to the quadrupole errors is [3] Δβ The corresponding phase advance perturbation is Thus the perturbation of the envelope stop band harmonics of Eq. (10) is ðe jðnϕþχ n Þ − e −jðnϕþχ n Þ Þðe jpϕ − 1Þds : Equation (10) can be used to calculate stop band harmonics arising from nonuniform spatial distribution of betatron phase advance. The effect of beta-beating and phase advance shift are linear with quadrupole errors. When the quadrupole field perturbation is small, the strength of each random envelope stop band harmonics is linearly proportional to the quadrupole perturbation [34]. Figure 4 shows the linearity of the betatron tune shifts and the dominating harmonics with the quadrupole errors. The quadrupole errors are defined by scaling a particular set of relative quadrupole errors, where one unit of scaling factor corresponds to 0.1% rms relative quadrupole errors with cutoff at three rms value. The red bar in Fig. 5 shows a set of the quadrupole errors of all 36 quadrupoles along the ring.

B. Emittance growth
With the random quadrupole error, we can calculate its effect on emittance growth through multiparticle tracking. The emittance growth can be categorized into two phases.  The first phase is an emittance leap in the first few turns during injection. This is caused by beam filamentation due to the mismatched betatron functions or the matched envelope radii. Due to the space charge force, the actual beta function mismatches with the original ring beta function. Higher space charge perveance causes bigger beta function mismatch therefore bigger emittance leap at the beginning. This part of emittance growth can be reduced by tweaking the beta functions of the injection beam radii to match the perturbed ones. The second phase of emittance growth is caused by the envelope dynamics. The perturbation and the driving source of envelope resonance are proportional to the linear space charge parameter and are also proportional to the quadrupole error. Therefore, the emittance growth rate, defined as the ratio of the final emittance and the initial emittance over a certain period, depends on the linear space charge parameter and the quadrupole error.
We carry out systematic simulation to study the effect of the emittance growth with different line densities and quadrupole errors. The lattice chosen is the same FODO lattice with 18 superperiods, with the bare tune designed at ðν x ; ν z Þ ¼ ð2.6; 2.9Þ. Therefore this lattice will encounter envelope resonance in the horizontal plane when the linear space charge parameter reaches ξ sc ¼ 0.2. The injected beam is matched with the unperturbed beta functions with no space charge force. The injection process takes 30 turns with increasing space charge kick. The initial beam emittance is ϵ x ¼ ϵ z ¼ 1 μm. We track 5000 macroparticles for 1000 turns.
The emittance growth vs line density and quadrupole error scaling factor in phases 1 and 2 are shown in Figs. 6(a) and 6(b). Here the quadrupole error scaling factor 1.0 corresponds to rms error of 0.1% to all FODO cell quadrupoles. By checking the color scale horizontally we find that, during the injection process (phase one), emittance growth depends on the space charge line density and is less sensitive to the quadrupole errors. The emittance growth during the storage or acceleration depends sensitively on quadrupole field error.
The number of survival particles is shown in Fig. 6(c). Few particles are lost if the line density is larger than N ¼ 1.6 × 10 9 m −1 . In phase two, the emittance growth is significant when the quadrupole error and the line density are large enough. If the quadrupole errors are controlled small, the emittance growth is less. Therefore wellcontrolled quadrupole error helps the reduction of the emittance growth.
Around N ¼ 0.1 × 10 9 m −1 and N ¼ 0.6 × 10 9 m −1 at high quadrupole errors, there is particle loss. This is due to the emittance growth caused by the envelope resonance. At these two line densities, the linear space charge tune shift parameter pushes the envelope tune to an integer or the envelope resonance.
The space charge tune shift parameters with different line densities and quadrupole errors at turn 30 and 1000 are plotted in Figs. 6(d) and 6(e), obtained from Eq. (4). Note that the linear space charge tune shift parameter is linear at low line density and will saturate when the line density is of the order ∼10 9 because the filamentation becomes very important during the injection process.
At one density there are many points plotted for different quadrupole errors. In Fig. 6(e) when the emittance grows more the tune shift is reduced. This is the reason that the curve smears. In particular the envelope resonance induced particle loss makes the spikes around N ¼ 0.1 × 10 9 m −1 and N ¼ 0.6 × 10 9 m −1 .

V. SPACE CHARGE ALLEVIATION AND COMPENSATION
From the perspective of beam dynamics there are some principles to minimize the space charge effects, extend the intensity limitation and minimize emittance growth.

A. Effect of ring circumference and betatron tunes
Since the linear space charge parameter is proportional to the space charge perveance and the ring circumference, and inverse proportional to the original beam emittance, a smaller ring with larger beam emittance can accommodate higher beam intensity.
The Laslett tune shift pushes the envelope tune downward. Before it hits an integer, there are rooms for tune shift. The rule is to adjust the bare tune so that it can accommodate more tune shift, i.e., more beam intensity. An ideal tune is 2ν 0 just below an integer stop band, so that the envelope tune at 2ν 0 − ξ sc will not encounter the next integer harmonic. Figure 3(c) shows the envelope tune with different line densities when the ν x bare tune is raised to 2.85. The line density can reach as much as N ¼ 3 × 10 8 m −1 before the envelope tune hit an integer. This is much better than the case in Fig. 3(a), in which case the bare tune is ν x ¼ 2.6.
Figure 6(f) shows the same thing as Fig. 6(b), except the horizontal tune is shifted from 2.6 to 2.8. Compared to Fig. 6(b), the envelope resonance around N ¼ 0.1 × 10 9 m −1 and N ¼ 0.6 × 10 9 m −1 disappears and a weaker resonance exists around N ¼ 0.3 × 10 9 m −1 . Furthermore, there is no loss in the first 1000 turns.
Another consideration about the betatron tune is to avoid the systematic space charge resonances [35]. An ideal betatron tune is below the resonance line because the space charge shifts the tune lower. In our lattice example with 18 superperiods, the systematic fourth and sixth order resonances are ν x ¼ ν z ¼ 4.5 and ν x ¼ ν z ¼ 3 respectively. The bare tune has better to be either below 4.5 and above 3; or simply below 3 in our simulation example.

B. Effect of injection rate
It is possible to shift the envelope tune across the envelope resonance. One way is to make the injection process faster. If the beam accumulates faster, the time for staying on the resonance will be shorter. Less particles will be trapped in the resonance therefore the less loss will be expected. By inspecting the FFT spectra with different injection turns shown in Fig. 7(a) we note that the envelope tune can go through the envelope stop band. In this case the injected line density is kept N ¼ 2 × 10 8 m −1 . Figure 7(b) shows the FFT spectrum peak distributions and the particle survival rate vs the number of injection turns. When the injection turn number is less than 7, the envelope tune can go through the stop band and more particles will survive.

C. Effect of quadrupole errors
The other way to tolerate a larger space charge tune shift is to reduce the strength of the random half-integer stop bands. The stop band integral can be reduced by decreasing the quadrupole errors.
Typical sources of quadrupole error are manufacturing error, quadrupole power supply error, misalignment of sextupoles, etc. The quadrupole error can be reduced by quality magnet fabrication and the orbit control and alignment. Figure 8 shows the envelope tune peaks and the particle survival rate with different quadrupole errors. When the quadrupole error is less then ΔK 1 L ¼ 0.004 m −1 the envelope tune can across the resonance in our simulation example.

D. Stop band harmonics correction by SVD
A set of trim quadrupoles can be excited to minimize the driving term in Eq. (8). For a given harmonic n, one needs to compensate the amplitude q n and phase χ n , i.e., each harmonic requires two trim quadrupoles to correct its amplitude and phase.
There were a few effective beta-beat correction methods in the literature [36]. For the space charge envelope stop band correction, these methods will not work because the beta-beat is most sensitive to the harmonics nearest to twice of the betatron tunes, while the space charge effect is sensitive to those integer envelope harmonic lower than twice the bare betatron tune. A conventional way to correct the stop band is to use the singular value decomposition (SVD) method.
The envelope tune without space charge is 2ν 0 . Therefore the dominant harmonics in our accelerator are q 5x , q 6x , q 5z , and q 6z . In the present space charge, the envelope tune is 2ν 0 − ξ sc . If the space charge perveance is large enough, the envelope tune will be shifted down below 5 or even 4. A set of ten trim quadrupoles properly distributed in the ring is used to compensate six variables q 5x , q 6x , q 5z , q 6z , ν x and ν z . They can be assigned with a x-before z-before x-after z-after (e)  x-before z-before x-after z-after  weight. First, we assume all equal weight. Figure 5 shows the quadrupole strength changes. The maximum change of quadrupole strength is less than 0.5% of the original quadrupole integrated strength, which is about Figure 9(a) shows the harmonics and the tune before and after correction. The fifth and sixth harmonics are effectively reduced. Figure 9(b) shows the horizontal and vertical beta-beat before and after the correction. They are effectively reduced too. However, the result of this SVD method is not good, as shown in Fig. 9(c), where the effect of beta-beat correction on emittance growth is compared. Note that the beta-beat correction may lead to higher emittance growth due to a higher stop band widths at lower harmonics. In our case, the linear space charge parameters are −0.48, −1.10, and −1.29 for N ¼ 2 × 10 8 m −1 , 6 × 10 8 m −1 , and 1 × 10 9 m −1 , respectively. When the linear space charge parameter is larger enough to bring the envelope tune below 5, the fourth and lower harmonics become very important. There is no more benefit from quadrupole setting for the beta-beat correction. The emittance even grows faster than before correction. To ease the emittance growth for high density the lower order harmonics correction must be considered.
The emittance leap in the first few turns is because of the filamentation from the mismatched betatron functions between the injected beam and the lattice perturbed by space charge. Higher space charge perveance causes bigger beta function mismatch therefore bigger emittance leap in the beginning.

E. Beta-beat correction
Since the beta-beat is caused by the quadrupole error along the ring, instead of stop band harmonics correction, we can also target beta-beat measurement and correction [36] so that the beta function is close to the ideal accelerator lattice. We use 18 beam position monitors along the ring to monitor beta functions and the weights of the betatron tunes are 200 times larger than those of beta functions. We also use ten trim quadrupoles for fair comparison with the previous method. The results are shown in Fig. 9(e). All stop band harmonics are effectively reduced as shown in Fig. 9(d). The emittance growth is also slightly reduced as shown in Fig. 9(f).

F. Stop band harmonics correction by other methods
Another approach is to use optimization algorithms [37] to minimize the stop band harmonics [Eq. (10)] while keeping the betatron tune the designed value. Since the lower harmonics near 2ν − ξ sc are important to space charge envelope dynamics, we target harmonics the third, fourth, and fifth harmonics in our accelerator model. In this case, we use 14 trim quadrupoles. The stop band harmonics and beta-beating before and after correction are shown in Figs. 9(g) and 9(h). Note that the emittance growth at N ¼ 2 × 10 9 m −1 has also been reduced to the value of a perfect accelerator, as shown in Fig. 9(i).

VI. CONCLUSION
This paper studies the space charge envelope dynamics by using and developing an efficient multiparticle tracking technique. Our tracking code can read accelerator design code and carry out analysis on space charge effects. With this tracking method, our numerical simulations demonstrated the importance of the envelope dynamics on beam stabilities.
The simulation code can also be used to study emittance growth due to the systematic fourth and sixth order space charge resonances. They can be avoided by carefully choosing the bare betatron tunes.
The next important space charge effect is the envelope resonances or the half-integer stop bands. We carry out a systematic study on the envelope resonance by multiparticle simulation. Our numerical simulations show good agreement with the KV-Sacherer-Lapostolle envelope equation. We carry out a systematic study on the emittance growth for different space charge perveance and quadrupole errors. The envelope tune should avoid the half-integer stop band to avoid the envelope resonance. A careful choice of the betatron tune will minimize emittance growth for space charge dominated beams.
Furthermore, we study the envelope resonance correction methods. We note that the stop band correction method at resonance near 2ν − ξ sc is the most effective method to reduce the emittance growth and beam loss. Although there are no successful experiments on stop band correction for the space charge envelope dynamics, our study can renew this effort on the space charge envelope stop band correction.