Dynamic aperture optimization with diffusion map analysis at CEPC using differential evolution algorithm

Nuclear Inst. and Methods in Physics Research, A 959 (2020) 163517

Contents lists available at ScienceDirect

Nuclear Inst. and Methods in Physics Research, A journal homepage: www.elsevier.com/locate/nima

Dynamic aperture optimization with diffusion map analysis at CEPC using differential evolution algorithm Jin Wu a,b ,∗, Yuan Zhang a,b ,∗, Qing Qin a,b , Yiwei Wang a , Chenghui Yu a,b , Demin Zhou c a

Key Laboratory of Particle Acceleration Physics and Technology, Institute of High Energy Physics, Chinese Academy of Sciences, Beijing 100049, China University of Chinese Academy of Sciences, Beijing 100049, China c KEK, Tsukuba, Ibaraki 305-0801, Japan b

ARTICLE

INFO

Keywords: Evolution algorithm Multi-objective optimization CEPC Dynamic aperture Beam lifetime Diffusion map analysis

ABSTRACT In storage rings, the dynamic aperture is a key concept of nonlinear beam dynamics to evaluate the overall machine performance. Therefore, its optimization becomes very essential, especially during the design stage of a machine. In recent years, with the progress in the field of parallel computing, various multi-objective optimization algorithms have been widely applied to the dynamic aperture optimization of storage ring light sources. In this paper, an alternative algorithm, so-called differential evolution algorithm is introduced to perform dynamic aperture optimization of a very high energy 𝑒+ 𝑒− storage ring collider, e.g. Circular Electron Positron Collider (CEPC). A multi-objective optimization code based on the differential evolution algorithm has been developed for this purpose and proved to be effective in increasing the dynamic aperture. In our code, a method called diffusion map analysis, in which the diffusion may come from quantum fluctuation of synchrotron radiation, beamstrahlung effect and nonlinearity in the lattice, is used to set the constraints in the dynamic aperture optimization, which can also help us find solutions with good beam lifetime.

1. Introduction The dynamic aperture of a storage ring is usually defined as a region of a given n-dimensional phase space, such that particles with initial conditions inside the region survive in tracking a certain number of turns. One can see that this is a computational definition of dynamic aperture since there is so far no consolidated theory or formulation about it. With given lattice configuration of the storage ring, the dynamic aperture can be computed using various accelerator simulation codes. Usually the dynamic aperture optimization has to be routinely done when the lattice configuration is changed. In the past years, following the evolution of parallel computing, various multi-objective optimization methods have been successfully applied to dynamic aperture optimization of storage ring based light sources. At APS [1], the genetic optimization algorithm was used to maximize the dynamic aperture and Touschek lifetime. Though the method was only based on particle tracking, the experimental tests validated the method, which brought significant improvements to APS operations. In NSLS-II [2], the multi-objective optimization of dynamic aperture was done with objective functions including both tracking data and analytical estimates of selected nonlinear driving terms. A strong correlation between the area of dynamic aperture and normalized driving terms was presented and it is found to be effective in

fastening the search of optimum solutions. As an alternative of multiobjective genetic optimization algorithm (MOGA), the multi-objective particle swarm optimization (MOPSO) was also proved to be very effective in nonlinear dynamics optimization. For the low emittance upgrade lattice of SPEAR3 [3], both methods were applied to direct optimization of machine performance using results from beam based measurements. The performance of these two algorithms were compared, and the results gave the advantages of MOPSO due to its independency on the initial distribution of solutions and faster convergence than MOGA. Multi-objective beam dynamics optimization is also used to minimize the final transverse emittances and maximize the final peak current of the beam in photoinjector design [4]. In contrast to genetic algorithms on which most previous studies were based, J. Qiang proposed a new parallel optimizer based on the differential evolution algorithm that performs much better in photoinjector beam dynamic optimization. In this paper, we show that the differential evolution algorithm is also effective in dynamic aperture optimization of a storage ring. In CEPC, the dynamic aperture is limited by a combination of various complicated effects. To maximize the dynamic aperture, a huge number of sextupole families (nearly 250) is always necessary. In FCCee, which has similar design goals as CEPC, nearly three hundreds of sextupole families per half ring are routinely used for dynamic aperture

∗ Corresponding authors. E-mail addresses: [email protected] (J. Wu), [email protected] (Y. Zhang).

https://doi.org/10.1016/j.nima.2020.163517 Received 11 October 2019; Received in revised form 31 December 2019; Accepted 22 January 2020 Available online 27 January 2020 0168-9002/© 2020 Elsevier B.V. All rights reserved.

J. Wu, Y. Zhang, Q. Qin et al.

Nuclear Inst. and Methods in Physics Research, A 959 (2020) 163517 Table 1 Main parameters of CEPC in CDR.

optimization during the optics design [5]. For comparison, there are in total 18 sextupole families in BEPCII and 54 sextupole families in SuperKEKB [6]. In a light source, the number of sextupole families is no more than 10 in nonlinear optimization, which is almost one order smaller than that in a modern collider. For example, only four families of multipoles are used for optimization in HEPS [7]. In this paper, we discuss the application of multi-objective optimization of dynamic aperture in CEPC. First, we briefly introduce the differential evolution algorithm, and show its potential to optimize the performance of a machine like CEPC, where hundreds of parameters are to be varied. To manage the inconsistency between dynamic aperture and beam lifetime in our solutions, we developed a method called diffusion map analysis, aiming to describe the diffusion rate in transverse amplitude space by tracking only few turns. The diffusion map presents good consistency with tracking results of beam lifetime for which usually more computing time is required. By setting constraints based on this method, the optimization usually gives a good solution of both dynamic aperture and beam lifetime.

Beam energy (GeV) Circumference (km) Synchrotron radiation loss/turn (GeV) Luminosity/IP (1034 cm−2 s−1 ) Number of IPs Momentum compaction factor (10−5 ) Betatron tune 𝜈𝑥 ∕𝜈𝑦 Synchrotron tune 𝜈𝑠 Energy spread (%) Emittance x/y (nm) 𝛽 function at IP (m) 𝜎𝑥 ∕𝜎𝑦 at IP (m) Bunch length 𝜎𝑧 (mm) Damping time 𝜏𝑥 ∕𝜏𝑦 ∕𝜏𝑧 (turns) Beam–beam parameter

120 100 1.73 3 2 1.11 363.10/365.22 0.065 0.134 1.21/0.0024 0.36/0.0015 2.1 × 10−5 ∕6.0 × 10−8 4.4 140∕140∕71 0.018/0.109

2. Differential evolution algorithm Similar to genetic algorithm, differential evolution (DE) is a population based optimization algorithm proposed by Storn and Price [8] in 1995. Mutation, recombination and selection are the crucial parts of DE since they are applied to every solution in each generation to create a new solution. Compared with other algorithms, the DE has mainly three advantages: finding the true global minimum value independent on the initialization of individuals, fast convergence, and using a few control parameters [9,10]. The major steps of the DE algorithm in each generation can be stated as follows:

Fig. 1. DA in horizontal direction with radiation damping on/off and without radiation damping from quadrupoles.

(1) Randomly initialize the population of NP parameter vectors; (2) Create the offspring solutions from the parent and combine them; (3) Identify the non-dominated solution set, which are assigned the front number; (4) Sort the total solutions and select the best NP solution as the parents in the next generation; (5) Return to step 2, if stopping condition is not satisfied.

3. Optimization of dynamic aperture in CEPC using differential evolution algorithm based code 3.1. Dynamic aperture of CEPC At the collision energy of 240 GeV, the CEPC will generate millions of Higgs particles in 7 years as a Higgs factory. The Conceptual Design Report (CDR) was published in 2018, which summarized the work of years [13]. The main parameters for Higgs of CEPC in CDR are listed in Table 1 [1]. At this energy, the dynamic aperture is significantly influenced by the strong synchrotron radiation. In addition to the radiation at dipoles, the energy loss from quadrupoles for particles at large amplitude is also important and consequently sets strong limit to dynamic aperture. With the radiation damping, the off-momentum dynamic aperture becomes larger, but the peak at 𝛿𝑝 = 0 disappears when considering the radiation from quadrupoles, compared with radiation only from the dipoles (see Fig. 1). The particles with large amplitude could not stay on-momentum due to the synchrotron motion caused by radiation from quadrupoles, such that the –I transportation between the non-interleaved sextupoles breaks down and nonlinear driving terms arise. Similar results could be seen in FCC-ee at 175 GeV [5], they gave a more detailed explanation about this figure. From tracking results, A. Bogomyagkov introduced a new mechanism called ‘‘self-inducing parametric resonance’’ by solving the equations of motion with radiation from quadrupoles and gave an analytical estimation of DA limit in FCC-ee at 45.6 GeV [14]. Consideration of quantum fluctuations shrinks dynamic aperture further, particularly in the vertical direction (Fig. 2), where the synchrotron radiation is mainly from the final focus quadrupoles.

Here NP is defined as the size of population, which is usually ten or more times the number of parameters. In our strategies, a specified solution is chosen to be set in the initialization of the first population and others are initialized randomly. We can also start a new task by using the last generation population of the old task in the conditions that our objective functions are changed. The Kung’s algorithm [11] is used to find the non-dominated set of solutions. The differential evolution algorithm for multi-objective optimization in our work is similar to Qiang’s work [4]. An associated mutant vector 𝑣𝑖,𝑗 in generation j is generated using the strategy defined by ] [ ] [ (1) 𝑣𝑖,𝑗 = 𝑥𝑖,𝑗 + 𝐹 × 𝑥𝑏,𝑗 − 𝑥𝑖,𝑗 + 𝐹 × 𝑥𝑟1,𝑗 − 𝑥𝑟2,𝑗 , where 𝑥𝑖,𝑗 is the ith population of generation j and 𝑥𝑏,𝑗 is the best parameter solution. The r1 and r2 are random integers generated in the range [1, NP]. Another parameter F in the equation is a scale factor chosen to be between 0 and 1 for each generation. If the new trial solution produces a better objective function value compared with its parent, it will be put into the next generation population. Otherwise, the old one will be kept. A key issue in multi-objective optimization is the selection of best solutions from current populations. The concept of domination is always used [12]. A solution 𝑥1 is said to dominate another solution 𝑥2 , if both the following conditions are true: (1) The solution 𝑥1 is no worse than 𝑥2 in all objectives. (2) The solution 𝑥1 is strictly better than 𝑥2 in at least one objective.

3.2. Multi-objective optimization of dynamic aperture Our code named MODE [15] was developed for dynamic aperture optimization of storage rings. It is a multi-objective code based on the 2

J. Wu, Y. Zhang, Q. Qin et al.

Nuclear Inst. and Methods in Physics Research, A 959 (2020) 163517

differential evolution algorithm. SAD [16] is used for particle tracking

excitation due to synchrotron radiation, can be implemented in each element. SAD is a parallel code, but the scalability is not very good. To reduce the computing time, MODE is used, as a MPI-based parallel code, to call SAD to do the optics calculation and particle tracking, which is much more scalable and efficient. In the conceptual design report of CEPC [13], the dynamic aperture is required to be 8𝜎𝑥 × 15𝜎𝑦 , with momentum acceptance 1.35% for the purpose of injection and beam–beam effect concerning beam lifetime [17]. In our optimization, radiation fluctuation in each element is considered in the particle tracking, where beam–beam effect is also included. According to Yang’s work [2], a larger DA area may not necessarily give a better solution because of the deviation from the shape of an ellipse. In order to reduce the random noise, we do the DA tracking more than once. Since the full DA tracking is time-consuming, we only track the on-momentum DA for extra 10 times and take the minimum value. The DA result is clipped to avoid abnormal DA score at some momentum deviation, such that the DA at large momentum deviation will not be larger than that at small momentum deviation when we try to evaluate the score, as shown in Fig. 3. After clipped, the shape of DA is closer to an ellipse. All the sextupoles, in total around 250 families, could be freed in the dynamic aperture optimization. However, the improvement in dynamic aperture becomes quite small when the number of variables exceeds 50. Therefore, in our DA optimization, we used 50 variables in total, with 32 sextupole families in arc region and 10 sextupole families in interaction region. In addition, 8 families of phase advance tuning knobs in straight sections are dedicated to adjust the nonlinear coupling between arc regions. Since the lattice design work is still ongoing, we only use an example version here as a test of optimization. Dynamic aperture was tracked for 100 turns (about one and a half damping times) with two initial phase (0, 0) and (𝜋∕2, 𝜋∕2).

and DA determination, where the dynamical effect, such as random

The objective functions are listed in the following:

Fig. 2. On-momentum DA in transverse space with radiation fluctuation and with radiation damping in each element.

Fig. 3. An example of DA clipping.

Fig. 4. The x-z aperture and y-z aperture of CEPC before optimization (a) and after optimization (b). The value of gradient color represents the average survival turn of particles. Lines in each picture are the apertures in four different initial phases, 0, 𝜋∕2, 𝜋 and 3𝜋/2, respectively, with 90% survival of 100 samples. Radiation fluctuation and beam–beam interaction is included. The number of turns is 145 (about two damping times). The coupling is 0.2%.. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

3

J. Wu, Y. Zhang, Q. Qin et al.

Nuclear Inst. and Methods in Physics Research, A 959 (2020) 163517

Fig. 5. On-Momentum DA in transverse space and beam lifetime in horizontal/vertical direction of the two cases. Case (a) has a larger DA but a shorter beam lifetime, while case (b) is opposite.

4. Beam lifetime

(1) The fractional parts of tunes 𝑄𝑥 and 𝑄𝑦 in a half ring by linear optics are kept in the range of [0.52, 0.58] and [0.58, 0.64], respectively, for 𝛿 ∈ [−0.005, 0.005], with the fixed tunes 0.55/0.61 for 𝛿 = 0; (2) The objective of horizontal aperture is defined with a boundary 𝑥2 𝑧2 of 25 2 + 222 = 1, where x is the horizontal amplitude in unit of RMS size with 0.2% coupling, and z is the energy deviation in unit of RMS energy spread; (3) The objective of vertical aperture is defined with a boundary of 𝑦2 𝑧2 + 22 2 = 1, where y is the vertical amplitude in unit of RMS 302 size with 0.2% coupling.

4.1. Estimation of beam lifetime Beam lifetime is very critical for the feasibility of the collider. Here we use the method in Ohmi’s BBWS code to calculate lifetime [23, 24]. In the equilibrium, the distribution of action 𝐽 is 𝑓 (𝐽 ), with ∞ ∫0 𝑓 (𝐽 ) 𝑑𝐽 = 1. Limited by the boundary A, the estimation of lifetime is 𝜏𝑖 T𝑖 = , (2) 2𝐴𝑓 (𝐴) where 𝜏𝑖 is the synchrotron radiation damping time in the corresponding direction. For a Gaussian distribution in phase space, the distribution of action 𝐽 is

The typical convergence time in one task is about 50 h with 500 CPU cores. The optimized solution enlarges the dynamic aperture significantly, especially for the off-momentum dynamic aperture. The tracking results of dynamic aperture before and after optimization are shown in Fig. 4. We calculated the four initial betatron phases of (0, 0), (𝜋∕2, 𝜋∕2), (𝜋, 𝜋) and (3𝜋∕2, 3𝜋∕2) in the transverse space in the dynamic aperture survey.

𝑓 (𝐽 ) =

1 𝑒𝑥𝑝(−𝐽 ∕𝜖𝑖 ), 𝜖𝑖

(3)

where 𝜖𝑖 (𝑖 = 𝑥, 𝑦, 𝑧) is the emittance in each direction. The lifetime of Gaussian distribution according to Eq. (3) is estimated as T𝑖 =

𝜏𝑖 𝜖𝑖 𝐴 𝑒𝑥𝑝 . 2 𝐴 𝜖𝑖

(4)

The number of tracking turns is 50 000. We do statistics with data of particles’ coordinates in the equilibrium (after seven damping times), and get the distribution 𝑓 (𝐴). Then lifetime could be calculated according to Eq. (2). The tracking process contains the effects of synchrotron radiation, beam–beam force and lattice nonlinearity, but it is very time-consuming.

3.3. Speed up method

In practice, the direct dynamic aperture tracking is usually very time-consuming. To speed up the optimization, it is very important to simplify the objective. For example, one option to save computing time is to track the particle with less number of turns in the early of optimization or in the process of trying new objectives. Among the different objectives, some are time-consuming, and some may be much faster. We can first optimize the fast objectives as constraints, such as chromaticities in the linear optics calculation, and then do the slow tracking only when the necessary constraints are satisfied. This method is justified by Ehrlichman’s work as demonstrated in [18]. Actually, the choice of objective and consequent speeding up of computation quite depends on if the physics is well define or not. As an alternative, the square matrix method introduced by Yu and others [19–22] is another suitable choice to save computing time.

4.2. Dynamic aperture vs. beam lifetime Dynamic aperture is usually strongly correlated to beam lifetime. However, it is found that larger dynamic aperture does not ensure longer lifetime in some cases (Fig. 5). The lifetime was obtained by Eq. (2). It should be mentioned that the beam–beam interaction is on during the tracking of dynamic aperture and lifetime. Turns of tracking for dynamic aperture (about 150) is much smaller than that for lifetime (about 50 000). Dynamic aperture is a result of short-term behavior but beam lifetime is a long-term result. A shorter lifetime means that more large-amplitude particles exist in equilibrium distribution 𝑓 (𝐴) 4

J. Wu, Y. Zhang, Q. Qin et al.

Nuclear Inst. and Methods in Physics Research, A 959 (2020) 163517

beamstrahlung effect all contribute to the diffusion and may lead to the halo distribution in coordinate space that result in a short lifetime. The frequency map analysis [25] method which is usually used in light source does not work due to the nonsymplectic condition. Inspired by some previous work [26–28], we attempt to describe the diffusion caused by the complicated effects through using a method based on particle tracking. In our method, hundreds of particles with the same initial amplitude and the same phase are tracked. As the number of turns grows, motions of the particles are quite different because of randomness. The amplitude of particle is defined as √ 2𝐽𝑖 𝑎𝑖 ≡ (𝑖 = 𝑥, 𝑦, 𝑧) , 𝜀𝑖

(5)

𝜎𝑎𝑖 is represented as the RMS value of particles’ amplitude. To combine the three directions, we do the geometric sum √ def 2 + 𝜎2 + 𝜎2 . 𝜎𝑎 = 𝜎𝑎𝑥 (6) 𝑎𝑦 𝑎𝑧 In our models, 200 particles with the same initial phase space coordinates were tracked. Fig. 7 gives the evolution of 𝜎𝑎𝑥 , 𝜎𝑎𝑦 , 𝜎𝑎𝑧 and 𝜎𝑎 with turns in different models with initial amplitude (𝑎𝑥, 𝑎𝑦) = (1, 10). The amplitude was selected according to the beam tail distribution in Fig. 6. We usually initialized the particles with phases (𝜓𝑥 , 𝜓𝑦 ) = (0, 0), and the tracking results is similar with another initial phases, such as ( 𝜋2 , 𝜋2 ). In the model without radiation, we initialized the 200 particles with little variance to measure the nonlinear effect in the motion. The results indicated that the radiation fluctuation and beam–beam effect (beamstrahlung effect included) are the two dominating effects of diffusion in three directions.

Fig. 6. The beam tail distribution in x–y space of the two cases in Fig. 5.

(Fig. 6), which is a non-Gaussian distribution. None of the particles is lost because the amplitude is far from the boundary of dynamic aperture.

Fig. 8 shows the evolution of variances with two different initial amplitudes. One has a larger amplitude in horizontal direction and one has a larger vertical amplitude. It is found that 𝜎𝑎𝑦 always takes the largest proportion, even in the case with large horizontal amplitude. The result agrees well with the significant decrease of dynamic aperture in vertical direction caused by radiation excitation. However, we only care about the first few turns when the diffusion increases before damping dominates the motion.

5. Diffusion map analysis 5.1. Definition of diffusion In such a high energy collider, the nonlinearity of lattice, strong synchrotron radiation fluctuation, beam–beam interaction and strong

Fig. 7. Evolution of variance in different directions and in different models.

5

J. Wu, Y. Zhang, Q. Qin et al.

Nuclear Inst. and Methods in Physics Research, A 959 (2020) 163517

Fig. 8. Evolution of variance in the model of radiation with beam–beam interaction at two different initial amplitudes. The left one is at (𝑎𝑥, 𝑎𝑦) = (1, 10) and the right one is at (𝑎𝑥, 𝑎𝑦) = (8, 1).

Fig. 9. The diffusion maps of four cases. The green lines are corresponding to the value of the function 𝑓 = 1.52. . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

5.2. Diffusion map

which have almost the same dynamic aperture but different lifetime. (Figs. 9 and 10)

Diffusion map is defined as an image of beam diffusion rate in transverse amplitude space [29]. According to the above analysis, the variance in first few turns could be a good index to describe the diffusion rate. In a general case, particles with large amplitude could have a large diffusion. On the other hand, a large diffusion usually means the motion near the corresponding amplitude in the map is more unstable. The function in diffusion map based on tracking is defined as

One can see that a clear difference of these lines between the first two cases and the last two cases is shown in the figure, especially in the region of large vertical amplitudes. The lifetime in vertical direction of the four cases by tracking are shown in Fig. 11. The diffusion map analysis presents good consistency with beam lifetime that needs much more turns of tracking.

( 𝑓 (𝑎𝑥, 𝑎𝑦) ≡ log10

∑

In our lifetime estimation, we usually track 1000 particles with Gaussian distribution for 50 000 turns with synchrotron radiation and beam–beam interaction. Here we also give the evolution of 𝜎𝑎𝑦 of these 1000 particles in our lifetime tracking of the four cases (see Fig. 12). 𝜎𝑎𝑦 is usually smaller in the cases with good lifetime. Thus, we established the relationship between beam lifetime and diffusion map. We try to use the diffusion map analysis as a figure of merit for halo particles distribution to do the dynamic aperture optimization.

) 𝜎𝑎2

.

(7)

𝑡𝑢𝑟𝑛

Here we tracked for only 25 turns at each initial amplitude with both radiation fluctuation and beam–beam effect included. The map is produced by scanning the value of 𝑓 in transverse space, with the step of (0.2𝜎𝑥 , 0.2𝜎𝑦 ), and it is shown in graded colors. Here we give the diffusion maps of four cases with different configuration of sextupoles 6

J. Wu, Y. Zhang, Q. Qin et al.

Nuclear Inst. and Methods in Physics Research, A 959 (2020) 163517

Fig. 10. The green lines of the four cases in Fig. 9 are plotted in one graph with different colors and enlarged in a certain region. . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 12. Evolution of variance 𝜎𝑎𝑦 in beam lifetime tracking of the four cases.

Several tests have been done to optimize the dynamic aperture using this constraint. Without focusing on the slight difference of dynamic aperture, Fig. 13 gives some tracking results of beam lifetime. Almost all the results present longer lifetime than that in first two cases in Fig. 11. Our experience shows shat the beam lifetime is usually bad in most cases without constraint about diffusion.

6. Optimization with diffusion map analysis Beam lifetime is not feasible to be an additional objective function in dynamic aperture optimization, as it takes too much time by tracking to be determined. Based on the correlation between lifetime and diffusion map, the information of the latter may be a good alternative. In our optimization strategy, minimizing the chromaticity is necessary as a constraint and the dynamic aperture is our objective. Here we add a constraint of diffusion rate: { 𝑐𝑜𝑛1 = 𝜉𝑥 , 𝜉𝑦 , 𝑐𝑜𝑛2 = 𝑓 (𝑎𝑥, 𝑎𝑦) , (8) 𝑜𝑏𝑗1 = 𝐷𝐴𝑥 , 𝑜𝑏𝑗2 = 𝐷𝐴𝑦

7. Conclusion A multi-objective optimization code based on differential evolution (MODE) is developed. Its functionality has been proved by the application in CEPC. We have shown that multi-objective optimization is an effective way of optimizing dynamic aperture for both on-momentum and off-momentum particles. An ongoing task is to speed up the convergence process for saving computing time. The synchrotron radiation is an important factor that has influence on both dynamic aperture and beam lifetime in the high energy collider. Some tracking results of solutions in optimization show that a larger dynamic aperture may not have a longer lifetime. We developed a method called diffusion map analysis to describe the influence of combined effects where radiation fluctuation and beamstrahlung dominate. This method presents good consistency with beam lifetime by tracking. Constraints of diffusion rate in amplitude space are used in optimization

where 𝜉𝑥 , 𝜉𝑦 is the chromaticity, 𝑓 (𝑎𝑥, 𝑎𝑦) is the value defined by Eq. (7). According to the tail distribution in x–y space (Fig. 6), 𝑓 (1, 10) is chosen based on our experience. We try to use the constraint on 𝑓 (𝑎𝑥, 𝑎𝑦) to select good solution where exists less large amplitude particles in equilibrium distribution, in other words, a longer lifetime. We can also add another constraint at a larger horizontal amplitude, such as 𝑓 (9, 1), to avoid that the solution does not have a good lifetime in horizontal direction. Actually, more constraints on 𝑓 (𝑎𝑥, 𝑎𝑦) can be added to diminish the influence of randomness if computing time is acceptable.

Fig. 11. The vertical lifetime of 4 cases.

7

J. Wu, Y. Zhang, Q. Qin et al.

Nuclear Inst. and Methods in Physics Research, A 959 (2020) 163517

Fig. 13. Lifetime in three directions of 4 selected optimized results.

of dynamic aperture, solutions always present good results of beam lifetime.

[2] L. Yang, Y. Li, W. Guo, S. Krinsky, Multiobjective optimization of dynamic aperture, Phys. Rev. ST Accel. Beams 14 (2011) 054001. [3] X. Huang, J. Safranek, Nonlinear dynamics optimization with particle swarm and genetic algorithms for SPEAR3 emittance upgrade, Nucl. Instrum. Methods Phys. Res. A 757 (2014) 48–53. [4] J. Qiang, Y. Chao, C.E. Mitchell, S. Paret, R.D. Ryne, A parallel multi-objective differential evolution algorithm for photoinjector beam dynamics optimization, in: Proc. IPAC’13. [5] K. Oide, et al., Design of beam optics for the future circular collider 𝑒+ 𝑒− collider rings, Phys. Rev. Accel. Beams 19 (2016) 111005. [6] Y. Ohnishi, et al., Accelerator design at SuperKEKB, Prog. Theor. Exp. Phys. (2013) 03A011. [7] G. Xu, Y. Jiao, Y. Peng, Chin. Phys. C 40 (2) (2016) 027001. [8] Differential evolution, available at http://www1.icsi.berkeley.edu/~storn/code. html. [9] D. Karaboga, S. Okdem, A simple and global optimization algorithm for engineering problems: differential evolution algorithm, Turk. J. Electron. Eng. 12 (1) (2004). [10] J. Vesterstrøm, R. Thomsen, A comparative study of differential evolution, particle swarm optimization, and evolutionary algorithms on numerical benchmark problems, IEEE Congr. Evol. Comput. (2004) 980–987. [11] H.T. Kung, F. Luccio, F.P. Preparata, On finding the maxima of a set of vectors, J. ACM 22 (1975) 469–476. [12] D.W. Corne, J.D. Knowles, M.J. Oates, The Pareto envelope-based selection algorithm for multi-objective optimization, in: Parallel Problem Solving from Nature PPSN VI. 6th International Conference. Proceedings, Lecture Notes in Computer Science 1917, Paris, France, 2000, pp. 839–848. [13] The CEPC study group, CEPC conceptual design report, arXiv:1809.00285. [14] A. Bogomyagkov, S. Sinyatkin, S. Glukhov, Dynamic aperture limitation in 𝑒+ 𝑒− colliders due to synchrotron radiation in quadrupoles, Phys. Rev. Accel. Beams 22 (2019) 021001. [15] Y. Zhang, D. Zhou, Application of differential evolution algorithm in future collider optimization, in: Proc. IPAC’16, Busan, Korea, 2016. [16] SAD home page, available at http://acc-physics.kek.jp/SAD/. [17] Y. Zhang, et al., Some Issues on beam-beam interaction at CEPC, in: Proc. eeFACT2018, Hong Kong, China, 2018. [18] M.P. Ehrlichman, Genetic algorithm for chromaticity correction in diffraction limited storage rings, arXiv:1603.02459.

Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. CRediT authorship contribution statement Jin Wu: Writing - original draft, Validation, Formal analysis, Data curation. Yuan Zhang: Conceptualization, Software, Funding acquisition, Writing - review & editing. Qing Qin: Supervision. Yiwei Wang: Resources. Chenghui Yu: Project administration, Resources. Demin Zhou: Methodology, Writing - review & editing. Acknowledgments The authors would like to thank Ji Qiang (LBNL) for his introduction of the DE algorithm. We also thank K. Ohmi (KEK) and K. Oide (CERN) for their continuous help and support. This work was supported by National Key Programme for S&T Research and Development, China (Grant NO. 2016YFA0400400) and National Natural Science Foundation of China (No. 11775238). References [1] M. Borland, Parallel tracking-based optimization of dynamic aperture and lifetime with application to the APS, in: 48th ICFA Beam Dynamics Workshop on Future Light Sources, Menlo Park, California, 2010. 8

J. Wu, Y. Zhang, Q. Qin et al.

Nuclear Inst. and Methods in Physics Research, A 959 (2020) 163517 [24] K. Ohmi, F. Zimmermann, FCC-ee/CEPC beam-beam simulations with Beamstrahlung, in: Proc. IPAC’14, Dresden, Germany, 2014. [25] J. Laskar, The chaotic motion of the solar system: A numerical estimate of the size of the chaotic zones, Icarus 88 (1990) 266–291. [26] D. Shatilov, E. Levichev, E. Simonov, M. Zobov, Application of frequency map analysis to beam-beam effects study in crab waist collision scheme, Phys. Rev. ST-AB 14 (2011) 014001. [27] K. Ohmi, M. Tawada, K. Oide, S. Kamada, Study of the diffusion processed caused by the beam-beam interactions, in: Proc. APAC’04, Gyeongju, Korea, 2004. [28] G. Stancari, et al., Measurements of the effect of collisions on transverse beam halo diffusion in the Tevatron and in the LHC, FERMILAB-CONF-13-054-APC. [29] J. Wu, Y. Zhang, Q. Qin, Diffusion map analysis in high energy storage ring based 𝑒+ 𝑒− collider, in: Proc. IPAC’19, Melbourne, Australia, 2019.

[19] L.H. Yu, Analysis of nonlinear dynamics by square matrix method, Phys. Rev. Accel. Beams 20 (2017) 034001. [20] Y. Li, L.H. Yu, L. Yang, Optimize nonlinear beam dynamical system with square matrix method, arXiv:1706.02195. [21] Y. Li, L.H. Yu, Using square matrix to realize phase space manipulation and dynamic aperture optimization, in: Proc. NAPAC2016: Chicago, Illinois, USA, 2016. [22] Y. Li, L.H. Yu, Applying square matrix to optimize storage ring nonlinear lattice, in: Proc. IPAC’17. Copenhagen, Denmark, 2017. [23] K. Ohmi, D. Shatilov, D. Zhou, Y. Zhang, Beam-beam simulation study for CEPC, in: Proc. IPAC’14, Dresden, Germany, 2014.

9