The consideration of transient scenarios requires describing EM fields under the quasi-static (QS) assumption [21]. In such a case, E _z is expressed as E z r , t = Re E z ̃ r , t e j ω t , (1) with j = − 1 representing the imaginary unit, ω representing the RF frequency, and E z ̃ representing the so-called quasi-static phasor. This extension of the standard phasorial formalism admits that the amplitude of the harmonic EM fields depends on time, and its variation is much slower than the RF oscillation itself.

Let E ₁ = E ₁(z, t) and E ₂ = E ₂(z, t) be two electric fields whose time-dependent quasi-static phasors are E ̃ 1 and E ̃ 2 , respectively. Furthermore, let f ̄ be the time average of a quantity f over an RF period T = 2 π ω , defined as f ̄ ( t ) = 1 T ∫ t t + T f ( τ ) d τ . Then, the following properties hold: E 1 z , t E 2 z , t = 1 2 Re E ̃ 1 z , t E ̃ 2 z , t * + Re E ̃ 1 z , t E ̃ 2 z , t e 2 j ω t , (2) E 1 z , t E 2 z , t ̄ = 1 2 T ∫ t t + T Re E ̃ 1 z , τ E ̃ 2 z , τ * d τ ≃ 1 2 Re E ̃ 1 z , t E ̃ 2 z , t * , (3) ∂ E 1 z , t E 2 z , t ∂ t ̄ = ∂ ∂ t 1 2 Re E ̃ 1 z , t E ̃ 2 z , t * ̄ ≃ ∂ ∂ t 1 2 Re E ̃ 1 z , t E ̃ 2 z , t * = ∂ E 1 z , t E 2 z , t ̄ ∂ t . (4)

Eq. 2 is obtained using the standard phasorial formalism, Eq. 3 is a direct consequence of the QS approximation [22], and Eq. 4 is obtained using Eqs 2, 3.

A first implication of the QS assumption is the temporal dependency of the following quantities defined in terms of the electric ( E ) and magnetic ( H ) fields: Volumetric density of stored EM energy  u r , t ≡ 1 2 ε 0 ‖ E r , t ‖ R 3 2 + 1 2 μ 0 ‖ H r , t ‖ R 3 2 , (5) Time - averaged stored EM energy density  u ̄ r , t ≡ 1 2 ε 0 ‖ E ̃ r , t ‖ C 3 2 + 1 2 μ 0 ‖ H ̃ r , t ‖ C 3 2 , (6) Stored energy per unit length in a cell w z k , t ≡ 1 L ∫ V k u ̄ r , t d V , (7) Time - averaged dissipated power in a cell P ̄ diss z k , t ≡ 1 2 ∫ V k R s ‖ H ̃ tan r , t ‖ C 3 2 d V , (8) Time - averaged dissipated power per unit length p diss z k , t ≡ P ̄ diss z , t L , (9) Time - averaged power flow through  S k P ̄ flow z k , t ≡ v g z k w z k , t . (10) where ɛ ₀ and μ ₀ represent the electric permittivity and magnetic permeability of vacuum, respectively; ‖ x ‖ R 3 ( C 3 ) 2 = ∑ i = 1 3 x i 2 represents the squared norm of an arbitrary real (complex) vector x = ( x 1 , x 2 , x 3 ) ∈ R 3 ( C 3 ) ; v _g represents the group velocity in the cell [23]; R _s represents the surface resistance; and H ̃ tan represents the magnetic field tangential to the surface. z _k denotes the beginning of the k–th cell, L denotes its length, and V _k denotes its volume, as shown in Figures 1A and B.

From Eqs 4–6, the following identity emerges as a consequence of the QS approximation: ∂ u r , t ∂ t ̄ = ∂ u ̄ r , t ∂ t . (11)

The definition of the effective accelerating gradient in a cell (G _eff) is also affected by the QS assumption: G eff z k , t ≡ 1 L ∫ z k z k + L Re E z ̃ z , t e j ω t q z d z , (12) where t _q (z) is the time of flight of a particle with charge q that enters the cell at time t ₀. In a single cell within the RF period, particles are assumed to travel with a constant velocity v, expressed in terms of the speed of light c, through the factor β as v = βc. Under this consideration, the expression of the time of flight is [24] t q z ≡ t 0 + ∫ z k z d z ′ v = t 0 + z − z k β c . (13)

From a physical point of view, G _eff is the mean electric field that affects a particle during its flight. The design of TW structures often requires on-crest acceleration, achieved when t _q is such that the field affecting the particle is | E ̃ z ( z , t ) | . In this case, G _eff matches the accelerating gradient (G), defined as G z k , t ≡ 1 L ∫ z k z k + L | E ̃ z z , t | d z . (14)

The magnitude relating G and G _eff is the time-transit factor ( T ) : T z k , t ≡ G eff z k , t G z k , t = ∫ z k z k + L Re E z ̃ z , t e j ω t q z d z ∫ z k z k + L | E z ̃ z , t | d z . (15)

To relate G with the stored EM energy and the dissipated power, the quality factor (Q) and the normalized shunt impedance per unit length ( r Q ) should be examined. These magnitudes, which depend only on the dimensions of the cell and its material [25], lead to the following relations [23, 24]: p diss z k , t Q z k = ω w z k , t , (16) G z k , t 2 = ω r Q z k w z k , t . (17)

The starting point for the derivation of a PDE describing gradient reduction in an accelerating structure is the Poynting theorem. Applied to the kth cell (see Figure 1B), it is expressed as [2] − ∂ u r , t ∂ t = ∇ ⋅ S r , t + E r , t ⋅ J r , t , r ∈ V k , t ∈ 0 , ∞ , (18) where S is the Poynting vector and J is the current density. Under the paraxial approximation (v _z ≫ v _x , v _y ), the current density can be expressed as J r , t ≡ ρ r , t v r , t ≃ ρ r , t c β z z k , t z ̂ , (19) where ρ is the volumetric charge density and β _z is assumed to be constant within a cell and period.

Eq. 18 is time-averaged over a period. Using Eqs 1, 11, and 19, considering only the on-axis longitudinal electric field, we obtain − ∂ u ̄ r , t ∂ t = ∇ ⋅ S ̄ r , t + Re E ̃ z z , t β z z k , t 1 T ∫ t t + T c ρ r , τ e j ω τ d τ . (20)

Then, Eq. 20 is integrated over the cell volume (see Figure 1B) and divided by L. Using the definition from Eq. 7 and denoting the longitudinal charge density as λ _q , this results in − ∂ w z k , t ∂ t = 1 L ∮ S ̄ r , t ⋅ d S + 1 L ∫ z k z k + L Re E ̃ z z , t β z z k , t 1 T ∫ t t + T c λ q r , τ e j ω τ d τ d z . (21)

For the second term, it should be recalled that the Poynting vector relates to the power balance in the cell as [2] ∮ S S r , t ⋅ d S = P diss z k , t + P flow z k + L , t − P flow z k , t . (22)

Regarding the time integral of the last term of Eq. 21, if a Gaussian bunch of charge q _bunch and bunch-size σ _z are assumed, then λ q ( z , t ) = q bunch 2 π σ z exp ( − ( t − t q ( z ) ) 2 2 σ z 2 / c 2 ) . Considering t _q (z) ∈ [t, t + T] and σ _z /c ≪ T, it can be simplified as β z z k , t 1 T ∫ t t + T c λ q r , τ e j ω τ d τ = e j ω t q z k β z z k , t q bunch T e − ω 2 σ z 2 / 2 c 2 = e j ω t q z k I ̃ z k , t . (23)

Substituting Eqs 22, 23 into Eq. 21 and using the definitions in Eqs 7–14, we obtain − ∂ w z k , t ∂ t = p diss z k , t + P flow z k + L , t − P flow z k , t L + G eff z k , t I ̃ z k , t . (24)

For structures exhibiting high-resonant frequencies, their cells are short enough to assume the limit L → 0 Eq. 24. This leads to the following expression: − ∂ w z , t ∂ t = p diss z , t + ∂ P ̄ flow z , t ∂ z + G eff z , t I ̃ z , t . (25)

Finally, substituting Eqs 10, 15–17 into Eq. 25, after manipulation, one obtains − ∂ G z , t ∂ t = v g z ∂ G z , t ∂ z + ∂ v g z ∂ z − v g z r Q z ∂ r Q z ∂ z + ω Q z G z , t 2 + ω r Q z T z , t I ̃ z , t 2 . (26)

Equation 26 shows that the gradient reduction in accelerating cavities takes place because of the power flow between cells, heat dissipation, and the beam–cavity interaction. The last term of Eq. 26 is called the beam-loading term, and it captures the influence of the beam in gradient reduction. It depends linearly on the charge and the normalized shunt-impedance per unit length, i.e., the design of the cavity.

For the case of an ultrarelativistic beam with on-crest acceleration ( T = 1 ) and point-like bunch (σ _z → 0), Eq. 26 matches the proposed expression in [17].

Standing-wave accelerating structures do not exhibit the power flow from cell to cell, meaning that the net input power entering the kth cell, P _input, has to be considered in the power balance since no flow from an external cell (coupler) feeds the structure. Therefore, the flux of the Poynting vector across the cell surface (see Figure 1B) has to be reconsidered as ∮ S S r , t ⋅ d S = P diss z k , t − P input z k , t . (27)

For the input power P _input, the transient response of a resonator to the power provided by the source has to be considered as follows: P input z k , t L = G unloaded z k 2 r z k 1 − exp ω t 2 ⟨ Q ⟩ , (28) with G _unloaded being the unloaded steady gradient that builds up in the cavity long-enough after the RF is turned on and ⟨Q⟩ being the loaded quality factor of all cells. However, r/Q is a geometrical factor that does not depend on coupling, so the following relation holds: r _loaded/Q _loaded = r _unloaded/Q _unloaded ≡ r/Q [8].

Starting from Eq. 21 and using Eqs 27, 28, and following the same steps as before, the following PDE is derived for gradient reduction in SW structures: − ∂ G z , t ∂ t = ω 2 Q z G z , t − ω 2 Q z G unloaded z 2 G z , t 1 − exp ω t 2 ⟨ Q ⟩ + ω r Q z T z , t I ̃ z , t 2 . (29)

RF-Track [19] is a tracking code developed at CERN, which allows arbitrary particle tracking under the effect of external and self-induced forces. It is written in parallel C++, and it is imported into Python and Octave as a pre-compiled library. This software program offers a flexible accelerator description, handling complex 3D fieldmaps and computing self-consistent single-particle and collective effects such as space-charge, incoherent synchrotron radiation, or short-range wakefields, among others.

A dedicated BL module has been implemented in RF-Track based on the resolution of Eqs 26, 29 using the finite-difference method, which requires the discretization of longitudinal and temporal coordinates [26]. For this purpose, let L _total be the total length of the accelerating structure where tracking is to be performed and T _end be the time at which the considered beam leaves the structure.

Then, a two-dimensional space–time mesh of N × M is created, where the longitudinal spacing is Δ z = L total N − 1 and the temporal spacing is Δ t = T end M − 1 . Each point of the mesh (n, m) refers to the position z _n = nΔz and time t _m = mΔt, with n = 0.(N − 1) and m = 0.(M − 1). As a consequence, the functions defined in Eqs 26, 29 are expressed as vectors and matrices evaluated in those mesh points.

Afterward, the derivatives are discretized by Taylor-expanding them to the first order. This leads to the following expression for the calculation of the gradient in TW structures: G n , m + 1 = Δ t B n G n , m + Δ t v g z n Δ z G n − 1 , m − C n , m , n = 1 . N − 1 m = 1 . M − 1 , (30) where B n = 1 2 d v g d z z n + d r / Q d z z n v g z n 2 r / Q z n − ω 2 Q z n − v g z n Δ z + 1 Δ t , (31) C n , m = Δ t ω r Q z n I ̃ z n , t m 2 . (32)

For a SW structure, the gradient is expressed as G n , m + 1 = 1 − ω Δ t 2 Q z n G n , m − C n , m + ω Δ t 2 Q n G n , 0 2 G n , m 1 − e ω m Δ t 2 ⟨ Q ⟩ . (33)

At each tracking step (where particles have a position r _part and time t _part), the BL effect is computed by applying an additional kick ( F _BL). For this, a TM₀₁ look-alike EM field ( E ₀₁ , B ₀₁ ) with gradient G* is calculated by solving Maxwell’s equations. G* is defined as G* = G(z _part, t _part) − G(z _part, 0) and is obtained by cubic-interpolating the result of Eqs 30, 33. Therefore, the BL kick is F BL r part , t part = q E 01 r part , t part + v r part , t part × B 01 r part , t part . (34)

Two different BL algorithms have been implemented in RF-Track depending on the velocity of the bunch. For ultrarelativistic bunches (v ≃ c) in TW structures, Eq. 32 can be precomputed prior to tracking, speeding up the overall process. If not, Eq. 32 has to be re-examined at each step of tracking. The flowcharts and differences of both approaches are shown in Figure 2.

The performance of the new RF-Track BL modules has been assessed by simulating two well-studied scenarios of the compact linear collider (CLIC), a TeV-scale high-luminosity linear electron–positron collider under development by international collaborations hosted by CERN [9].

The two situations where beam-loading effects have been reproduced are CLIC accelerating structures (CLIC ASs), which consist of an X-band TW structure where the main beam is accelerated, and CLIC power extraction and transfer structures (CLIC PETS), which are passive cavities where CLIC’s drive beam is decelerated due to the BL effect and power is extracted to feed CLIC AS.

The structural specifications required for the simulations are shown in Table 1.

The CERN Linear Accelerator for Research at CERN is a versatile user facility providing a 200-MeV electron beam for accelerator R&D, irradiation studies for space, and medical applications [29]. It consists of a 25-m-long electron beam injector, called CALIFES [30], and a 16-m-long user beamline. The BL measurements have been performed in CALIFES, whose schematic representation is shown in Figure 5.

The electron beam is generated on a Cs₂Te photo-cathode upon laser incidence and is captured by an RF gun [31]. Thereafter, three 4.5-m-long TW structures are placed to bunch and accelerate the beam. Both the RF-gun and the TW structures are fed by two independent 3-GHz modulators/klystrons, and a rich set of beam diagnostics is set along the beamline to monitor the beam. The main parameters of the accelerating structures are displayed in Table 1.

A beam of 150 bunches with f _b = 1.5 GHz was injected into the CLEAR RF gun. For different q _bunch, the beam induced energy losses have been measured as exposed in Eq. 37. The results are shown in Figure 7A. The same scenario was simulated with RF-Track, and the beam-induced energy loss (E _{loss, BL}) dependency with q _bunch is illustrated in Figure 7B.

A linear dependency of the beam-induced energy loss with q _bunch is found, as reflected by the high value of the r ² coefficient from the linear fit of both experimental and simulation data. This is consistent with the beam-loading term in Eq. 29, which depends linearly on I ̃ and therefore on charge.

The positive sign of the slope (m) confirms that the higher the charge of the injected bunch, the greater the gradient reduction due to BL. A good agreement is found between the slopes obtained from the experimental data and RF-Track’s simulation, as shown in Table 2, where a slope deviation (δ _m ) of 3.4% is obtained for this case.

The validity of the BL model in non-ultrarelativistic SW structures relies on experimental results as to the knowledge of the authors, no tracking code incorporates this feature yet for β < 1. For that reason, additional measurements were taken in the CLEAR gun. In this case, a beam with a fixed charge per bunch of q _bunch = 294 ± 24 pC was injected into the structure.

Figure 8 shows the dependency of E _{loss, BL} on the number of bunches. A similar conclusion can be drawn from this comparison: larger charges lead to greater energy losses as a consequence of beam-loading-induced gradient reduction, and accordance between simulations and measurement is obtained, as reflected by the large r ² value shown in Table 2.

Finally, Figure 9 shows the evolution of the bunch-averaged time-transit factor from bunch to bunch. A slight variation in the time-transit factor from bunch to bunch is observed due to the BL-induced gradient reduction. Lower energy gains lead to lower velocities along the structure, which result in particles arriving at different times to the starting position of the cell, thus seeing a different phase of the accelerating electric field.

A beam of 50 bunches is injected into CLEAR, and its energy is measured in VESPER, as shown in Section 3.3.2. Several measurements have been taken for different values of the charge per bunch, ranging from 100 to 400 pC. Figures 10A, C show the dependency of ⟨E⟩ and E _range with the charge per bunch. Figures 10B, D show the analogous simulations in RF-Track.

For both ⟨E⟩ and E _range, an agreement is found between measurements and simulations since results deviate only up to 5.9% and high r ² values are reported (see Table 2). In addition, the results indicate the same dependency of energy loss with beam charge, which is consistent with Eqs 26, 29.