Traveling wave method for simulating geometric beam coupling impedance of a beamscreen with pumping holes

In particle accelerators, pumping holes in a vacuum chamber can be a source of unwanted broadband coupling impedance, leading to beam instabilities. Analytical methods have been previously developed to estimate the impedance of holes in circular-like chambers e.g. the beamscreen of the Large Hadron Collider (LHC). More sophisticated chamber designs like that of the High Energy LHC (HE-LHC) and the Future Circular Collider (FCC-hh) call for a different way to calculate the impedance. We propose using decomposition of the wakefield into synchronous traveling waves and employing a numerical solver to find the impedance of each wave. This method is compared to the direct time domain wakefield calculation method and its greater sensitivity to small impedances is shown.


I. INTRODUCTION
In hadron colliders such as the Large Hadron Collider (LHC), the proposed High Energy LHC (HE-LHC) and the Future Circular Collider (the hadron-hadron option FCC-hh), the beamscreen separates the particle beam from the magnet cold bore to avoid having the synchrotron radiation heat load on the cold mass (Fig. 1). The pumping holes connect the space inside the beamscreen to the outer region from where the air is pumped out. In this paper, we propose a new method for calculation of the geometrical beam coupling impedance of the pumping holes. In particular, we focus on the imaginary part of the longitudinal and the transverse impedances Im(Z || (f )) and Im(Z ⊥ (f )). In the mentioned colliders, estimating these quantities at the frequencies of the bunch spectrum (up to ∼ 1 GHz) is critical to ensure single-bunch beam stability. Power dissipation through the holes, related to Re(Z || (f )), is another important aspect of the impedance of the holes, but is beyond the scope of this paper.
For a beamscreen with a circular cross-section, analytical models based on Bethe's theory can be applied [1][2][3][4]. However, purely analytical models are not applicable if the cross-section of the beamscreen is far from a circle, which is the case for the FCC-hh and the HE-LHC beamscreens due to the additional shielding (Fig. 1). A semianalytical approach suggested in [5] employs Bethe's formalism together with 2D simulations to estimate the fields at the holes. However, even this method is not applicable to the FCC-hh and the HE-LHC beamscreens because the holes cannot be considered small and the fields inside a hole vary by two orders of magnitude. More than that, even for a circular beamscreen cross-section, theoretical models can be wrong by a factor of a few as will be shown in section III B. This error comes from neglecting the interference between consecutive holes and from not perfectly accounting for the influence of the outer * sergey.arsenyev@cern.ch chamber.
FIG. 1. FCC-hh beamscreen (green) inside a magnet cold bore (grey). Only the features important to the study are shown, which excludes the cooling channel and the surface coatings.
Therefore, computational tools are needed.
Often, time domain wakefield simulations with codes like CST [6], GdfidL [7], ACE3P [8] are used to calculate impedance of complex 3D structures. However, in the case of a beamscreen with pumping holes, the impedance per unit length can be so low that time domain simulations become not suitable due to the limited numerical accuracy.
The traveling wave method solves this problem by relating the impedance of the holes to a sum over waves traveling in the periodic structure. This technique by itself is not new and was applied before to estimate the impedance of a periodically corrugated beam pipe [9], [10]. However, to the best knowledge of the authors, impedance calculation through decomposition into traveling waves was never done for pumping holes (the method was briefly mentioned without a proof in an earlier study by the authors [11]). The method can only be applied to beamscreens where the positions of the holes follow a regular pattern. This is the case for the FCC-hh and the HE-LHC beamscreens, where the positions of the holes are not randomized (like in the LHC) in order to simplify the manufacturing and reduce the cost. The pattern repeats uninterrupted for the length of a cryo-dipole and part of the interconnect (15.3 m in total) with about 900 periods, making the beamscreen well suited for the approximation of an infinite periodic structure.
An important difference between the implementation of the traveling wave method in [9], [10] (for the case of a corrugated pipe) and the implementation in this paper is that the search for the synchronous waves cannot be purely analytical, but requires eigenmode simulations. Nevertheless, as a simulation tool, such method offers a valuable advantage over time domain simulations, as it allows to calculate much smaller impedances. This advantage is related to the fact that the method only requires solving for eigenmodes in a single period of the structure, making the computation much lighter.

II. DESCRIPTION OF THE METHOD
To employ the method, we first find the traveling waves synchronous with the beam. For this, one period of a structure is simulated using an Eigenmode solver, with the phase advance over the period φ given as a parameter and the eigenmode frequencies ω read as the output. The phase advance over one period is scanned over the Brillouin zone to obtain dispersion curves ω n (φ) (blue lines in Fig. 2). We then find intersections (φ syn n , ω syn n ) of these curves with the synchronous line of the particle beam. In this paper the beam is assumed to be ultra-relativistic (β rel ≈ 1) as is the case in FCC-hh and HE-LHC, although this is not a necessary condition to apply the method. For an ultra-relativistic beam the synchronous line on the dispersion diagram is ω = φc/L, appearing as multiple branches in the Brillouin zone (red lines in Fig. 2). Here c is the speed of light and L is the length of one period of the structure.
Then, for all the synchronous waves (φ syn n , ω syn n ), the electromagnetic fields are integrated along the beam path to obtain the longitudinal and the transverse voltage kicks V where s is the longitudinal coordinate, t is the time, and F ||,n and F ⊥,n are the complex longitudinal and trans-−π 0 π Phase advance per period, ϕ verse Lorentz forces experienced by a test particle of charge e. These voltages are then used to compute the shunt impedances (R/Q) of the waves. In this paper we use the "circuit" definition for shunt impedances: where U syn n is energy stored in one period of the n-th synchronous traveling wave.
Finally, the computed (R/Q) are used to find the longitudinal and the transverse beam coupling impedances per period of the structure Z || (ω)/N and Z ⊥ (ω)/N : where N is the number of periods, factors α are the corrections due to non-zero group velocities v g of the waves = 0 appears when the structure allows propagation of TEM-like waves with zero cut-off frequency (n = 1 band in Fig. 2). It is crucial for a calculation of the longitudinal impedance of the pumping holes, as will be demonstrated in section III B.
The relations (3) resemble those that can be obtained for ordinary standing wave cavities [12], with the exception of the additional factors α n . The factors α n are sometimes ignored for similar problems (see, for example, [9]). In many cases, this does not result in significant errors because of a small group velocity v g c. However, in the case of the pumping holes, a consistent disagreement in simulations prompted a strict derivation of the factors (see Appendix). For the case of non-zero synchronous crossing (ω syn n = 0) the factors were also experimentally observed for CLIC PETS structure [13] and later derived in [14] and [10]. While taking a different path, the derivation in Appendix leads to the same result for ω syn n = 0 and a new result for the previously not considered case ω syn n = 0.

III. APPLICATION AND BENCHMARKING OF THE METHOD
Below we consider three different types of periodic structures: • simple bellows • pumping holes in a circular pipe for which analytical models can be applied

• FCC-hh beamscreen
For each structure, the traveling wave analysis is done as described in section II. One period of a structure is independently simulated using the Eigenmode solvers of Ansys Electronics Desktop [15] and the CST Microwave Studio [16], making sure that the results agree. In both solvers, the tetrahedral meshing method is used to better approximate the curved geometries. Since we focus only on the geometrical impedance, all metal walls are considered to be made of a perfect conductor.
A MATLAB [17] script is used to externally launch the eigenmode solver and to read out the results. The script scans points along the synchronous line (red line in Fig. 2) and finds all the intersections with the mode dispersion curves (φ syn n , f syn n ) in the desired frequency range. Finding the intersection points consumes most of the computational time. Once the points are found, the longitudinal and the transverse shunt impedances are calculated by integrating the electromagnetic force along the beam path. Since all the considered structures possess mirror symmetries in X and Y , symmetry planes are applied: the magnetic symmetries in both XZ and Y Z for the longitudinal impedance, or the magnetic symmetry in XZ and the electric symmetry in Y Z for the transverse impedance in the X-direction. In all considered structures the transverse impedance in the Y -direction is lower than (or equal to) the transverse impedance in the X-direction and is not discussed.
The traveling wave method is compared to the Wakefield solver of CST [6] -a well-established time domain solver for wake and impedance calculations. The impedance computed in the CST Wakefield solver is a result of a direct Fourier transformation of the wake function, with both the longitudinal and the transverse impedances measured in Ω. We normalize the transverse impedance by the transverse offset to convert the result to Ω/m which allows for the comparison. For the wakefield computation, a 10-period long structure is terminated with open boundary conditions at the ends. Since the open boundary does not represent an infinitely repeating structure, reflected waves from the ends produce an unwanted contribution to the impedance. To cancel out this contribution, a 20-period long structure is simulated and the difference between 10 and 20 periods was considered for impedance calculation. Care is taken to make sure that this difference is reproduced when 30 periods are compared to 20 periods.

A. Bellows
The bellows structure shown in Fig. 3 is one of the simplest geometries to simulate, hence a good agreement between the different methods is expected. Dispersion curves for the first five longitudinal and transverse modes are shown in Fig. 4. Results of the traveling wave method are shown in Tables I and II for the longitudinal and the transverse cases, respectively. The corresponding total impedances given by Eq. 3 agree within 5% with the results of the CST wakefield solver (also shown in Tables I and II). For this simple structure, the impedances are dominated by the lowest frequency traveling wave in both the longitudinal and the transverse cases. Notice also that if the correction factors α were not taken into account, the results would have been underestimated by 40%.

B. Pumping holes in a circular pipe
We now consider a simple structure with pumping holes, shown in Fig. 5. It consists of a circular pipe with two circular holes per cross section, opening up to a cylindrical outer region. This simple geometry allows comparing calculated impedances to an analytical model by S. Kurennoy [1]. According to this model (and confirmed by the later study accounting for the outer coaxial region [4]), for circular holes the imaginary part of where N is the number of periods, Z 0 = 377 Ω is the impedance of vacuum, R h is the hole radius, b is the inner radius of the pipe, and M is the number of holes per cross section (M = 2 in the considered structure). For the considered geometry, the effect of the wall thickness on the polarizabilities can be neglected (see Fig. 5 of [18]). An important complication relative to the case of bellows arises from the fact that the outer region allows propagation of a TEM-like mode, similar to a coaxial line. This mode is represented by the lowest blue line in the corresponding dispersion diagram (Fig. 6), almost matching the synchronous line. The exact intersection with the synchronous line appears at the point φ = 0, f = 0 and is treated as a zero synchronous crossing as described in section II.
Results of the traveling wave method are shown in Ta cases, respectively. In case of the longitudinal impedance, the contribution of the first synchronous wave is much greater than that of all the other modes combined. This is due to the very high factor α caused by the dispersion curve of the first band almost matching the synchronous line. This shows the importance of treating the case ω syn = 0 separately. Indeed, if we had mistakenly used the non-zero crossing expression for α (Eq. (4), top), the final answer would have been a factor of 2 off. In case of the transverse impedance, the contributions are spread over a large number of traveling waves. Modes up to the frequency of 20 GHz were taken into account, with 143 modes in total. To check that the number of modes is sufficient, the evolution of the impedance is plotted as a function of the number of modes in Fig. 7. After summing more than 50 traveling waves (f > 15 GHz), the impedance converges to a value which is taken as the total impedance.
In Tables III and IV, the results of the traveling wave calculation are also compared to those obtained by the wakefield solver and by the analytical model of Eq. 6. All three methods agree within 15%. However, while the agreement between the traveling wave method and the wakefield solver is always to be expected, the agreement with the analytical model is only achieved thanks to the careful choice of the geometry parameters. In particular, a large distance between the holes L and a large radius of the outer wall R out are necessary for Eq. 6 to correctly  To prove this statement, we consider an example of a geometry that better reflects the features of a real beamscreen, and call it "Geometry B" (as opposed to the already analyzed "Geometry A"). In particular, we reduce the period length L to 20 mm to make the holes more tightly packed, and reduce the outer chamber radius R out to 22 mm to reflect the fact that the beamscreen occupies most of the cold bore space. As far as the analytical approach is concerned, the changes in L and R out have no impact on the impedance. This, however, is not the case, as is shown in Table V, where the results of the two computational methods are compared to the analytical model for both geometries. One can see that the analytical approach works well only for Geometry A, but significantly overestimates the impedance for Geometry B. Thus, even if the FCC-hh beamscreen had a circularlike cross-section, the theory could have been only used to get a rough (within a factor of a few) estimate. We consider the 2018 version of the FCC-hh beamscreen design (Fig. 1), which is different from the previous designs considered, for example, in [11]. An important difference from the case of the circular pipe (section III B) comes from the addition of the hole shielding. Because of the shielding, the holes are not directly visible to the beam and instead are coupled to the beam region through a slit. The narrower the slit, the stronger the shielding of the holes, and, consecutively, the lower the hole impedance. The low impedance per period of the structure, however, still needs an exact estimate, as with 4.8 million periods in the FCC-hh it can reach a critical level. When estimating low impedances, the traveling wave method is expected to have a significant advantage over the wakefield solver. This is due to the fact that only one period of the structure needs to be simulated, and due to the more accurate conformal tetrahedral mesh available in eigenmode solvers.
In order to calculate very low impedances with the traveling wave method, special care is taken to reduce the numerical noise floor of the eigenmode simulations. For that, the mesh nodes are forced to lie on the voltage integration line, and the mesh steps along the line are made much smaller than in the rest of the volume (Figure 8). The synchronicity condition is strictly enforced, as the phase advance is adjusted until the frequencies of the traveling waves to lie within 10 kHz from the synchronous line. To measure the noise level, we simulate a modified geometry for which the impedance is known to be strictly zero. For that, the consecutive holes are connected together creating an infinitely long continous hole. The cross-section of such structure is a constant of the longitudinal coordinate, leading to the zero geometrical impedance. Therefore, for each traveling wave, a small numerically obtained impedance gives the noise floor. The mesh and eigenmode solver parameters were chosen to minimize this noise floor. For the actual FCC-hh beamscreen design, the wakefield solver fails to give an impedance estimate as the impedance is far below the computational noise. Nevertheless, the wakefield solver can be compared to the traveling wave method for geometries with wider slits, where the impedance is not low. Such comparison is shown in Fig. 9 for the slit width w ranging from 7.5 mm (the actual design) to 24.44 mm (the screen is completely removed). For the case of w = 24.44 mm the two methods agree within 20%, but the corresponding lines start diverging for narrower slits due to the wakefield results affected by the computational noise. In the range w < 15 mm the wakefield solver becomes unusable, while the traveling wave solver gives converging results even for w = 7.5 mm. Thus, the sensitivity of the traveling wave method to low impedances exceeds the one of the wakefield solver by three orders of magnitude.
For narrow slits w < 15 mm, the curves for Im(Z || )(w) and Im(Z x )(w) approach straight lines on the log-log scale (Fig. 9). These lines correspond to a very steep dependence on the slit width that can be approximated by Im(Z) ∝ w 9.5 , confirming that the shielding is very effective at reducing the impedance of the holes. For the actual design (w = 7.5 mm) the traveling wave method gives Im(Z || ) per period /f = 4.49 × 10 −7 Ω/GHz Im(Z x ) per period = 4.48 × 10 −4 Ω/m.
The corresponding dispersion diagram and the convergence check are plotted in Figures 10 and 11, and the eigenmode data is presented in Tables VI and VII for the first 40 modes. The total impedance of 4.8 million periods of the holes is Im(Z || ) total /n = 6.7 × 10 −6 Ω Im(Z x ) total = 2.2 × 10 3 Ω/m, where, as often used for stability studies, the longitudinal impedance is normalized over n = f /f rev , and the revolution frequency f rev = 3.067 kHz. In the FCChh, the instability threshold for Im(Z || )/n is in the mΩ range [19] and the instability threshold for Im(Z x ) is in the MΩ/m range [20]. Therefore, as far as the beam stability is concerned, the impedance of the holes in the FCC-hh beamscreen is completely negligible.

IV. CONCLUSIONS
We have shown that the traveling wave method is well suited for calculation of the longitudinal and transverse impedances of pumping holes. The method was benchmarked against the CST time domain wakefield solver and the analytical formulas for three different types of geometries. The method agrees well with both approaches in the cases when they are applicable. More than that, the method gives converging results even in the cases when one or both of the compared appoaches cannot be applied. In particular, the analytical approach only works for simple geometries, and only for some sets of parameters (section III B). The wakefield solver, on the other hand, fails to calculate very low impedances due to the numerical noise (section III C). As a consequence, the traveling wave method is the only option for the FCChh beamscreen, where the geometry is complex and the impedance per unit length is very low. In this case, the method provides a three orders of magnitude better sensitivity than the wakefield solver.
The transverse impedance of the pumping holes in the FCC-hh beamscreen is found to be negligible for single bunch instabilities. This is a consequence of the hole shielding that reduces the impedance by more than a factor of 30000.

V. ACKNOWLEDGEMENTS
This work was supported by the European Union's Horizon 2020 research and innovation programme under grant No 654305. We also thank Walter Wuensch for the useful discussions.

APPENDIX: CORRECTION FACTORS DUE TO THE GROUP VELOCITY
Below we will show the validity of Eq. (3), and, specifically, derive the correction factors α(v g ) listed in Eq. (4). To do that, we represent an infinite structure as a limiting case of an N -periods long structure, with N → ∞. In a finite structure, instead of a continous dispersion curve only discrete points φ p = pπ/N exist (the circles in Fig. 2). Each point corresponds to a standing wave, and for each standing wave, the resonator impedance model [12] can be applied: where ω res is the resonant frequency, Q is the quality factor, and R || and R ⊥ are defined according to Eq. (2). In order to make the transverse resonator formula consistent with the adopted definition for R ⊥ (in Ω), the R ⊥ (measured in Ω/m in [12]) was multiplied by the additional factor (ω res /c). By decomposing Eq. (9) into the powers of ω and summing over all standing waves (without any correction factors α n ), we get Here the symbol Σ means that the shunt impedances are not per-period, but the total impedances of the structure. We will show that the sum over the indexes n × p of the standing waves can be reduced to the sum over indexes n of the synchronous traveling waves with the correction factors α n (v g ), as in Eq. 3.
Let us consider only one mode n and omit the index n below. To calculate the shunt impedances (R/Q) Σ p , we first write down the total voltage kick V Σ p (either longitudinal or transverse). The electromagnetic force F of a standing wave can be expressed as a sum of the forward and the backward waves F 0 (s)e iωps/c e −iφps/L + e iφps/L ds (11) where e is the elementary charge and F 0 (s) = F 0 (s+L) is the periodic envelope function for both the forward and the backward waves, and the index p can be any integer in the range 0 ≤ p ≤ N − 1. Using the periodicity, we write where the first summand corresponds to the +p branch and the second summand corresponds to the −p branch in Fig. 2. Futhermore, we can assume that for each p, only one of the branches matters. Indeed, if the considered wave (φ p , ω p ) lies close to the synchronous line ω = φc/L on the dispersion diagram, its counterpart (−φ p , ω p ) will be far from synchronous and will not amount to a significant impedance. Therefore, we can simplify Eq. (12) by leaving only the first summand, and exdending the range of p to −(N − 1) ≤ p ≤ N − 1.
The integral in the first square brackets of Eq. (12) is equal to V p -the voltage of the p-th traveling wave in a structure consisting of only one period. Therefore where θ p = Lω p /c − φ p is the phase slippage per period. A geometrical representation of Eq. 13 is depicted in Fig. 12. The farther away the point p is from the synchronous point (a cross in Fig. 2), the lower is its total voltage. Mathematically it can be expressed as |V Σ p | = N |V p |g N (θ p ), where the function g N (θ) describes the attenuation of the voltage sum due to asynchronicity. From Eq. (13) it can be shown that This compact result for |V Σ p | is only possible thanks to the assumption that only one of the forward and the backward waves matters. If both waves were considered, |V Σ p | would have to depend on the two voltages as |V + + V − |, where V + and V − correspond to the terms in Eq. (12).
In one period of the structure, the energy stored in a standing wave is twice that stored in the corresponding traveling wave. Over N periods, the energy of the standing wave becomes U Σ p = 2N U p . Given |V Σ p | and U Σ p , the standing wave shunt impedance can be related to the per-period impedance of p-th traveling wave: We then plug this expression in Eq. 10 and take the limit N → ∞ to obtain the contribution of the n-th resonance to the total impedances As N increases, the function g N (θ) becomes a sharper and sharper peak around θ = 0. Using this fact, the quantities (R/Q) The only remaining part is to prove that the coefficient α(v g ) can be written as a function of the group velocity in the form of Eq. (4). To do that, the cases ω syn = 0 and ω syn = 0 are treated separately. Let us first treat the case ω syn = 0. In the vicinity of the synchronous point, the dispersion curve can be approximated by a straight line ω p = ω syn + (φ p − φ syn )∂ω/∂φ = φ syn c/L + (φ p − φ syn )v g /L. The phase slippage per period becomes The coefficient α(v g ) can be expanded as where we have defined ξ = π 2 (1 − v g /c) and p * = φ syn N/π. As N increases, the terms in the denominator approach ξ(p * − p), and α approches We can then use a remarkable property of the sinc function ∞ k=−∞ sinc 2 (ak − b) = ∞ k=−∞ sinc 2 (ak) = π/a and rewrite α in a more compact way Now, let us consider the remaining case ω syn = 0. The sum in Eq. (17) has to be split in two parts because for the negative p the group velocity has the opposite sign v g (−φ p ) = −v g (φ p ), as shown in Figure 2. By following the same steps as for the case ω syn = 0, we arrive to which proves Eq. (4).