Numerical Modeling of a Proof-of-Principle Experiment on Optical Stochastic Cooling at the IOTA Electron Storage Ring

,


I. INTRODUCTION
Particle-beam cooling [1,2] − the process of reducing beam emittances − has been an essential element in the success of accelerator-based science [3,4].For example, in the case of particle colliders, beam cooling is critical for increasing and preserving the collider's luminosity.A broad range of cooling techniques have been developed [2,[5][6][7][8].Of particular importance, the method of stochastic cooling (SC) [9,10] enabled the production of intense antiproton beams, which subsequently led to the discovery of the W and Z bosons in 1983 at the Super Proton Synchrotron (SPS) [11].SC uses an electromagnetic pickup to measure position information associated with short longitudinal slices of the bunch [12].This information is encoded in the temporal structure of a radio frequency (RF) signal, typically in the microwave regime, amplified, and then used to correct the slices' average positions at a downstream kicker device.The correction is applied turn-by-turn as the beam circulates in the ring.At each turn, the sampled slices contain different randomized groupings of particles.Over many turns, each particle's coherent contribution to its own cooling signal comes to dominate over the diffusive contribution of neighboring particles, and thus, the beam is "stochastically" cooled [13].
SC has been successfully implemented at many accelerator facilities for both particle accumulation and cooling [11,[13][14][15][16]. SC has also been successfully demonstrated and operationalized at the collision energy in the Relativistic Heavy Ion Collider (RHIC) for moderate-intensity bunched beams of both protons and gold ions [17,18].Unfortunately, due to its limited band- * adick1@niu.eduwidth, microwave SC becomes ineffective at the higher intensities that are typical of, e.g., proton-(anti)proton colliders.
In the case of optimal cooling, where the system gain is balanced against diffusive effects, the maximum damping rate is approximately where ∆f is the bandwidth of the integrated system, N the number of particles in the bunch, σ s the root-mean-square (RMS) bunch length, and C is the circular-accelerator circumference.For reference, with a typical proton beam in the Large Hadron Collider (σ s ≃ 10 cm, C ≃ 30 km, N ≃ 10 11 ), an SC system with a bandwidth ∆f ≃ 4 × 10 9 Hz, would have a damping time on the order of τ ≃ 2 × 10 3 hours.
Effective cooling of such beams would require an increase in damping rates by at least three orders of magnitude.A possible solution is to significantly increase the bandwidth of the integrated system.One such extension of SC to optical frequencies and bandwidths, optical stochastic cooling (OSC), was proposed in the early 1990s [19].The OSC mechanism can support an optical bandwidth of up to ∆f ∼ 100 THz, potentially increasing the achievable damping rate by four orders of magnitude.The basic principles of OSC are similar to those of SC; however, OSC replaces the conventional microwave pickups and kickers with undulator magnets that deflect the beam particles to produce and couple to optical radiation.This radiation can be amplified using an optical amplifier similar to those employed in high-power free-space lasers.The possible use of OSC to cool hadron [20,21] and lepton [22][23][24][25] beams has been proposed, and OSCbased techniques capable of manipulating the beam, e.g.halo control, have also been explored [26,27].
The OSC technique was first demonstrated in 2021 at Fermilab's IOTA storage ring using low-charge electron beams [28].The original experiment implemented a "passive" version of OSC where optical radiation emitted in the pickup interacts with the beam in the kicker without prior amplification.An amplified version of the experiment is currently under development at Fermilab.
To guide OSC R&D activities, we have developed a computational model of OSC that can simulate the full OSC mechanism on a turn-by-turn basis.The model was benchmarked against the data collected during the passive OSC experiment at the IOTA ring.This model supplements the previous analytical models of OSC and will aid in the study of various aspects of the phase-space dynamics that are inaccessible in experiments.It will also serve as a valuable tool in the development of techniques for advanced beam control.

II. TRANSIT-TIME OSC (TTOSC)
A. Theoretical Background 1.General principle OSC extends the principle of SC to optical frequencies by using undulator magnets for the pickup and kicker [19,29].For a planar undulator, the on-axis fundamental wavelength of undulator radiation (UR) is given by where λ u is the undulator period, γ is the Lorentz factor, and K ≡ eλuBu 2πmc is the undulator parameter where e and m are respectively the elementary charge and the particle's rest mass, c the speed of light, and B u the peak magnetic field in the undulator [30].
The original OSC concept included a quadupolar pickup undulator to produce emission based on the fluctuations of each sample slice, a form similar to conventional SC [19].In subsequent work, it was shown that OSC could be implemented using conventional undulators separated by a delay line that provides an energydependent path length for the particles -a configuration referred to as transit-time OSC (TTOSC) [31].Figure 1 depicts the TTOSC concept: the pickup (PU) and kicker (KU) undulators are separated by a bypass beamline which introduces a nominal pathlength s 0 for the reference particle.A particle with phase-space coordinate Z P U,i ≡ (x P U,i , x ′ P U,i , y P U,i , y ′ P U,i , s P U,i , δ P U,i ) T at the PU location (with respect to the reference particle) will experience a modified path length [31][32][33] where R 5j refers to the elements of the 6x6 transfer matrix between the PU and KU, (x P U,i , x ′ P U,i ) are the hori- (a) A TTOSC system with pickup and kicker undulators, a bypass chicane, an optical line, and (b) the position of a particle and its own undulator radiation at the exit of the pickup undulator, and at the entrance and exit of the kicker undulator.The particle slips behind the radiation one wavelength each undulator period and the bypass introduces a delay so that the particle arrives at the kicker at the front of the radiation wavepacket.The delay can be tuned such that the subsequent energy exchange in the kicker produces corrective kicks that cool the beam.
zontal phase-space coordinates and δ P U,i is the fractional momentum deviation.Equation 3 considers a bypass beamline acting on the 4D phase space (x, x ′ , s, δ) so that the motion in (y, y ′ ) is decoupled.We introduce the phase shift associated with the i-th particle as where k ≡ 2π λ = ω c and ∆t i and ∆t 0 are the time of flight for the i-th particle and the reference particle respectively (i.e.∆t i = si βc where β ≡ 1 − γ −2 ).The radiation emitted by the particle in the PU propagates to the KU through an optical transport line which includes the optical amplifier, imaging lenses, and a delay-control system [34].The optical line introduces an optical path length (OPL) ℓ.The delay-control system must maintain an OPL so that the optical delay δℓ ≡ s 0 − ℓ = λ/4.This ensures that the reference particle's periodic maxima in transverse velocity are in phase with the zero crossings of the PU-radiation field, thus resulting in zero net energy exchange.
The energy change experienced by a particle in a single pass will depend on its phase difference with respect to the reference particle ∆ϕ i as [34] where K represents the maximum possible energy change, L u is the length of the undulators, and ψ 0 is the phase offset which represents a temporal shift between the beam and radiation packet and is described further in Section III.The strength of the focused PU-radiation field is computed at the KU location.The latter equation also assumes that s i − ℓ ≤ λ.Equation 6 represents the interaction of a particle with its own PU radiation while passing through the KU; this is commonly referred to as the coherent contribution to the kick.The coherent kick is the one that ultimately contributes to phase-space cooling.The UR wave packet produced by a particle as it passes through the PU can be approximated as a truncated sinusoidal wave with a total duration N u λ/c.By the end of the PU the wave packet is located just ahead of the particle; see Fig. 1(b).This slippage length S ≡ N u λ arises because the radiation propagates faster than the particle, which "slips" by one optical wavelength (λ) per undulator period [35].In order to maximize the overlap of the interaction within the KU, the particle needs to arrive just ahead of the wave packet; see Fig. 1(b).If the unwrapped phase is larger than 2π (i.e.s i − ℓ > λ) then the kick strength K reduces as discussed in Section III B.
The cooling mechanism in TTOSC relies on the particles' longitudinal dynamics between the PU and KU.By its nature, TTOSC only cools the longitudinal phase space; however, the cooling force can be redistributed to the other phase planes by modifying the bypass and lattice optics.Generally, the cooling process leads to an emittance reduction parameterized as where m = x, y, s refers to the considered degree of freedom, ε m,0 and ε m,∞ are respectively the injection and equilibrium emittances, and τ m is the damping time.Throughout this paper, the emittance values correspond to the geometric emittances statistically computed from the macroparticle distribution [36] as Generally, τ m = J −1 m τ 0 with the partition numbers J m connected by J x + J y + J s = 4 [37] and τ 0 being the damping time associated with the considered energy-loss process (e.g.synchrotron radiation and/or OSC in the present work).Consequently, the damping decrements satisfy The ratio between damping decrements associated with the different phase-space planes can be controlled by the horizontal dispersion and the lattice coupling.Specifically, the cooling force can be shared with the horizontal plane by introducing dispersion in the undulators and modifying the PU-KU beamline settings to control the path-length dependence on horizontal coordinates (x, x ′ ) in Eq. 3. The cooling force can be shared with the vertical plane by coupling the horizontal and vertical planes using skew-quadrupole magnets or by operating the storage ring on a transverse-coupling resonance.

Multi-particle interactions
In addition to the coherent kick discussed above, a particle also experiences the field produced by other particles whose radiation wave packets overlap with the considered particle.These "incoherent" contributions of the neighboring particles can be derived similarly with an additional phase difference due to their relative separation.
The radiation wave packet produced by the j-th particle interacts with the i-th particle in the KU if |t i − t j | < N u λ/c, where t i and t j are the particles' arrival times at the PU center.Introducing the relative phase between the particles, ϕ ij = ω(t i − t j ), the effective kick produced by the j-th particle on the i-th particle can be calculated in the same way as Eq. 6 using ϕ ij as the relative phase term.The total incoherent kick experienced by the i-th particle is then obtained by summing over all particles having time coordinates within the temporal slice, i.e.
and is given by The incoherent kick for a given particle thus depends on the number of particles within a slippage length to either side and corresponds to the summation of the wave packets produced by each particle [38].

Transverse effects
So far our description of the TTOSC method has assumed the UR wavepacket to be a plane wave.We may also account for the non-uniform transverse field distribution of UR and the particles' transverse motion between the PU and KU.The UR is emitted within a cone with apex angle θ ∼ 1/γ and is then imaged in the kicker.Ignoring prefactors and assuming an undulator parameter K ≪ 1 [34,39], the far-field distribution of the electric field in the KU has the form where ξ ≡ rkθ 1+(γθ) 2 , J 0 and J 2 are Bessel functions of the first kind of order n = 0, 2 respectively, and θ is the polar angle of a given ring-shaped surface element on the focusing lens relative to the center of the PU.The horizontal component of the electric field at the focal point in the KU is then proportional to the integral over the angular acceptance θ m of the lens, where r ≡ x 2 + y 2 .In the absence of depth-of-field effects, the UR from the PU is imaged into the KU with a total transfer matrix M = −I, where I is the identity matrix [40].If the transfer matrices for the particles and light are not matched in the transverse planes, then the particles may arrive in the KU off-axis relative to their PU radiation.The relative transverse offset of a particle is determined by the difference in transverse coordinates at the pickup and kicker.The offset of a single particle is defined as ∆x i = x U,i − M x P U,i .The electric field sampled by the i-th particle in the PU is then E x (∆x i , ∆y i ).
The correction for the total energy change experienced by an off-axis particle is given by substituting K → Kϱ(x, y) in Eq. 6 where It should be noted that, for a given OSC system, the finite angular acceptance of the optical system also decreases the maximum energy change [39] where the suppression factor F T () depends on the angular acceptance of the optical system θ m as detailed in [34].Therefore the transverse effects leads to Eq. 6 being modified as The previous equation together with Eq. 10 assumes a passive-cooling configuration.When an optical amplifier is included both of these equations should include a multiplicative factor √ G, where G represents the opticalpower gain.

B. OSC proof-of-principle experiment configuration at IOTA
The OSC process was experimentally demonstrated at Fermilab's IOTA storage ring [41][42][43] and data collected during the experiment [28] are used to benchmark the computational model in the next sections.Table I summarizes the main experimental parameters and a diagram of the system appears in Fig. 2(a).The OSC bypass includes four dipole magnets arranged as a chicane, four quadrupole magnets for control of beam focusing and dispersion, a single quadrupole (QX1) to control the coupling of OSC between the longitudinal and horizontal planes, and three sextupole magnets for mitigation of nonlinear path lengthening in the bypass.
For the uncoupled lattice [with computed lattice functions displayed in Fig. 1(c)], OSC only cools the longitudinal phase space; however, the cooling force can be distributed between the longitudinal and horizontal planes by powering QX1, which is located at the symmetry point of the particle bypass with significant horizontal dispersion (see Fig. 2).Additionally, the cooling force can further be shared with the vertical plane using the skew-quadrupole magnet SQ1 that introduces coupling between the (x, x ′ ) and (y, y ′ ) transverse phase  spaces [43].This latter configuration supports threedimensional cooling.These two quadrupole magnets were the only lattice elements tuned during the bench-marking process to investigate the cooling dynamics in the three coupling configurations.
The primary metrics on which we evaluated the model were (i) the cooling rates of various lattice configurations and (ii) the equilibrium longitudinal and transverse beam distributions.A complete lattice of the IOTA ring was simulated in elegant [44] using magnet settings from the experiment lattice; see Fig. 2(c).

III. MODELING OSC IN ELEGANT
We opted to implement the OSC process in the elegant particle-tracking program owing to its large user base and open-source nature [44].Elegant has been extensively used for the design of storage rings and linear accelerators.Particle tracking in elegant can be performed using matrices of selectable order, canonical kick elements, numerically integrated elements, or any combination thereof.Canonical kick elements are available for dipole, quadrupole, and sextupole magnets, and higher-order multipoles.All of these elements support optional classical synchrotron radiation losses, as well as random variation in radiation, leading to quantum excitation.The particle-tracking implementation of most of the beamline elements is parallelized on CPU [45] and some have GPU capabilities [46].The code also supports parallel simulation of multi-bunch beams with short-and long-range wakefields [47], beam-loaded rf cavities with feedback [48], and bunch-by-bunch transverse and longitudinal feedback, all of which are of interest in modeling realistic beam behavior.The parallelization of elegant is critical to performing detailed simulations of the longterm beam dynamics and understanding the cooling process.

A. General considerations
The main components of the OSC implementation in elegant are the cpickup and the ckicker beamline elements.Both are zero-length elements.The cpickup records the internal 6-D coordinates of every macroparticle while the ckicker element applies a momentum kick.The two elements form a linked pair within elegant using an identifier string.In its simple implementation, OSC is modeled by placing a cpickup and ckicker in the center of respectively the KU and PU.Such a "single" kick approximation captures most of the features associated with the OSC mechanism.However, the kick can also be applied in a distributed fashion using multiple cpickup-ckicker pairs distributed along each of the undulators.Such a configuration allows for the simulation of angular misalignment, and for investigating the impact of optical focusing and magnification errors.Elegant uses the 6-D coordinates (x, x ′ , y, y ′ , s, δ) where s is the total, equivalent distance traveled, and δ ≡ (p − p 0 )/p 0 (where p 0 is the reference-particle momentum) is the fractional momentum deviation [49].
At present, the OSC model has three main components discussed in Section II A: (i) coherent kicks, (ii) incoherent kicks, and (iii) transverse effects.The following three sections discuss the specifics of the implementation associated with these three processes.

B. Coherent kicks
The coherent energy kick described by Eq. 6 is implemented as a relative momentum change experienced by the i-th macroparticle following where ∆ϕ i = ω(∆t i − ∆t 0 ) and the variables κ, ω, and ψ 0 , are user-defined parameters for respectively the maximum normalized kick strength, radiation frequency, and unwrapped optical-delay phase.The function Ψ() describes the temporal overlap as detailed below.The 6-D coordinates of each macroparticle are recorded in the PU using the pickup element.At the location of the kicker element, the time of flight from the location of the kicker element is computed as ∆t i = t i,kU − t i,P U ; this corresponds to the path length s i ≃ c∆t i for ultra-relativistic macroparticles.The reference particle time of flight ∆t 0 is computed either from the bunch average time of flight, ∆t 0 ≃ s 0 /c = ⟨t KU ⟩ − ⟨t P U ⟩ where ⟨...⟩ indicates the statistical averaging over all the macroparticles, or from the closed orbit trajectory found by elegant.
The delay phase, ψ 0 , is an important parameter in the model as it affects both the phasing of the OSC effect and the temporal overlap of the beam and radiation pulse.Under normal operation, ψ 0 = 0 and a particle with zero momentum deviation will arrive in phase with its radiation, as illustrated in Fig. 1(b); however, if there is an additional delay of ψ 0 = ±π in the optical line, the reference particle arrives antiphased with its radiation and will receive a kick with the opposite sign and be driven away from the design momentum.As the UR is periodic, if the delay is an entire period (ψ 0 = ±2π) then the system will return to the cooling mode but with a slightly reduced cooling force due to the imperfect temporal overlap.The energy is only transferred from the radiation to the beam for as long as the two are in the kicker together.Thus, as the delay increases, the total energy transferred decreases.We model this effect using the unwrapped delay phase as The nominal kick in Eq. 6 is multiplied by Ψ() yielding Eq. 16.This causes the kick strength to decrease linearly with the optical delay until it vanishes when the particle and its undulator-radiation wavepacket no longer overlap.

C. Incoherent kicks
As described above, the field experienced by a particle in the kicker is the superposition of the field produced by all particles in its sample slice z i ± N u λ.In order to simulate the effect of neighboring particles, we consider their individual contributions to the corrective kick instead of their contributions to the electric field in the kicker [50].Each neighboring particle will affect the i-th particle in a similar way to the coherent contributions but with an unwrapped phase term ψ ij determined by the difference in arrival time at the pickup.
The incoherent-kick contribution to the i-th particle is modeled by considering the contributions of all other particles j within the slippage interval |t i − t j | ≤ N u λ/c at the PU.The phase difference ϕ ij = ω(t i − t j ) at the PU is combined with the coherent phase for particle i at the KU to ensure the correct is kick received from each particle j, as described in Eq. 10.The kick reduction due to the total optical delay described by Ψ() also influences the incoherent kick.In practice, these contributions become more significant as the longitudinal beam density increases but are negligible for the considered low-charge beams.In our numerical implementation in elegant, this effect is parallelized and can be toggled in the lattice file using the incoherentMode parameter.

D. Transverse effects
The implemented TTOSC model accounts for the effect of the transverse motion of particles on the time of flight described in Eq. 5. Additionally, the macroparticle may enter the KU with a transverse offset thereby sampling off-axis regions of the electric field as discussed in Section II A 3. We account for this effect by calculating the relative change of the field with respect to the on-axis field, ϱ(x, y), sampled by the macroparticle with transverse offset (x, y); see Eq. 13.The user-defined maximum kick parameter κ in Eq. 16 is modified as The transverse motion of particles is calculated by subtracting the position in the kicker from the position in the pickup.The transverse displacement between UR emitted from the pickup and imaged in the kicker element and the position of the i-th macroparticle in the kicker element is given by xi where M is the optical magnification between the pickup and kicker.The magnification is a user-defined parameter (M = −1 by default, corresponding to a single-lens optical-imaging system such as employed in the IOTA experiment).The relative field is calculated for each macroparticle using Eq. 12 with the correction factor ϱ(x i , ỹi ).The model does not currently include the transverse effects for the incoherent contributions and instead uses the on-axis strength.This approximation is acceptable for modeling the IOTA experiment as the incoherent effects are negligible (owing to the low number of electrons per sampling slice).The algorithm for modeling the transverse effects is parallelized and may be toggled in the model using the transverseMode input parameter in elegant.

IV. COMPARISON OF NUMERICAL AND EXPERIMENTAL RESULTS
An important motivation for developing a precise numerical model of OSC is to gain insight into the single particle evolution usually inaccessible in experiments.Elegant incorporates several diagnostics that provide either ensemble-averaged information on the beam distribution or phase-space coordinates of each macroparticle at a given location on a turn-by-turn basis.In particular, elegant's watch element was extensively employed to export macroparticle distributions at several locations around the IOTA storage ring.This capability permitted the reconstruction of the data collection methods used at IOTA while providing insight into the dynamics of individual macroparticles in the phase space.

A. Measurements
The IOTA experiment relies on several diagnostics systems to measure the beam distribution.The transverse (x, y) distribution of the beam is captured by two complementary metal-oxide-semiconductor (CMOS) cameras that image the synchrotron radiation produced in the M1L and M2R dipole magnets; see Fig. 2(a).The required exposure time depends on the beam intensity.Typical values were between 1 and 100 ms, corresponding to many thousands of beam revolutions.The temporal distribution is measured using a streak camera installed at the M3R diagnostics station.The streak camera uses a continuous sweep voltage and is phase-locked to the 11 th harmonic (82.5 MHz) of the circulation frequency.Correspondingly, the elegant simulations are configured to record the macroparticle distributions at the center of the M3R bending magnet for direct comparisons with the measurement.The signal can be integrated over multiple turns to mimic the response time of the cameras.
Elegant allows the user to toggle most advanced effects independently, e.g.synchrotron radiation energy loss and quantum excitation effects may be applied individually.This feature is also built into the model of OSC so that incoherent contributions and transverse effects can be added to the basic transit-time model.The measurements of emittance damping in Tables III, IV, and V do not include quantum excitation or other diffusive effects.This reduces random noise in the emittance calculation and does not impact the observed cooling rate as diffusive effects only increase the equilibrium term in Eq. 7.

B. Comparison of the nominal IOTA lattice with elegant simulations
The lattice configuration and magnet strengths were taken from the IOTA experiment in the Summer of 2021.The three coupling configurations can be controlled using the quadrupole magnet QX1 located in the OSC bypass chicane and the skew quadrupole magnet SQ1 near the dipole magnet M1R.The full ring and OSC beamline is diagrammed in Fig. 2(a,b).The PU and KU are modeled as 16 periods of alternating dipole magnets with the two outermost dipoles on either side forming a standard [1/4,3/4] termination scheme [51]; see Table I.
Prior to investigating the OSC process in detail, we compared results obtained from the elegant model with available nominal values for some of the lattice parameters in the absence of OSC.The theoretical values were initially calculated for the designed lattice [43] and were recomputed, using various optics codes, for slight changes in the final lattice used in the experiment.The parameters of the final lattice are considered here and are summarized in Table II.During the experiment, the synchrotron-oscillation frequency was measured for RF-cavity voltages within [70, 140] V.It was ultimately set to 105 V corresponding to a synchrotron frequency of f s = 426.5 Hz in agreement with elegant value within 3%.Likewise, the synchrotron-radiation damping dynamics and the equilibrium emittances were simulated with elegant.This is particularly important for the OSC studies since the OSC damping rates and equilibrium emittances are directly compared to the SR damping ones.We simulate a beam initially above equilibrium and allow it to circulate for several seconds corresponding to ∼ 25 × 10 6 turns.The damping decrement and equilibria, fit from these tracking simulations, are in good agreement with theoretical estimates, indicating that the baseline lattice and performance are properly simulated in elegant; see Table III.The simulations outlined here consider the theoretical and experimental cases separately.The simulations of the theoretical design use the nominal values for angular acceptance, kick strength, and coupling terms.In the latter simulations, we attempt to recreate the conditions of the experiment by reducing the angular acceptance of the focusing lens, the total kick strength of the OSC element, and the strength of the coupling quadrupole magnet.We also include a simple scattering element to mimic residual gas scattering.

Comparison with design values
In the theoretical design, the total emittance damping rate due to OSC, in the absence of synchrotron radiation, scattering, transverse field effects, and incoherent contributions, is τ −1 tot ≃ 38 s −1 [43].These effects can each be disabled individually in the Elegant model to give a clear picture of the unaltered OSC damping.When simulating the "base" OSC effects in this way, we measure a combined transverse damping rate of τ −1 x,y = 19.1 s −1 and a longitudinal damping rate of τ −1 s = 18.0 s −1 corresponding to a total rate in all planes of τ −1 tot = 37.0 s −1 .Next we simulated OSC in each of the coupling configurations while including synchrotron radiation effects and the transverse field calculation.Beginning with the uncoupled lattice, we set the parameters of the pickup and kicker elements to their theoretical values.In this case, the total longitudinal damping rate increases to τ −1 s = 29.6 s −1 from the one initially computed with SR damping alone (τ −1 s = 2.03 s −1 ).The horizontal and vertical rates remain the same, as determined by SR damping.In the coupled cases, the nominal ratio of cooling rates between the transverse and longitudinal planes is τ −1 s : τ −1 x = 1.00 : 1.03.We then set the strength of QX1 to the theoretical value corresponding to 1:1.03 coupling.Here, the longitudinal damping rate decreases to τ −1 s = 16.35 s −1 and the horizontal damping rate becomes τ −1 x = 16.49s −1 .Finally, the SQ1 quadrupole is activated and the cooling force is seen in all three planes.The longitudinal damping rate is approximately the same, and the transverse cooling is now evenly split between the horizontal and vertical degrees of freedom with τ −1 x = 9.3 s −1 and τ −1 y = 9.2 s −1 .This demonstrates the expected redistribution of the total cooling force between the three phase-space planes as the coupling terms are introduced.Table IV summarizes the results of these simulations.The total damping rate, summed across all dimensions, is consistent between the three configurations.Subtracting the SR-damping rates in each case leaves the OSC rates, which are observed to be ∼ 25% lower than the simulated total rate discussed above (τ −1 tot = 37.0 s −1 ).This can be explained by the reduced kick strength that results from the transverse effects described in Section III D. Figure 3 presents the evolution of the projected distributions for the three coupling configurations.As expected, in the uncoupled case, only the longitudinal beam size is reduced.As the OSC force is shared with the other planes, the equilibrium size in the longitudinal plane grows larger while the equilibrium beam size in the transverse planes decreases thus confirming the redistribution of the cooling process among the coupled planes.

Comparison with IOTA experiment
The IOTA experiment serves as an important benchmark for the computational model.In the experiment, the strength of OSC was characterized by measuring the projected beam distribution along each direction (x,y, and z) and accounting for the effects of intrabeam scattering (IBS) using a simple model [28].Relative to SR damping alone, the IOTA experiment measured an increase in damping of 8.06 times in the longitudinal and of 2.94 times in the transverse.This corresponds to a total OSC emittance cooling rate of τ −1 tot = 18.4 s −1 and is approximately half of the theoretical estimate.Several sources were identified as possible explanations for this difference [28].The first was a reduced, asymmetric aperture for the pickup radiation, due to vacuum-chamber OSC is turned on instantaneously at t = 0 ms (indicated with the vertical dash lines).The "1D", "2D", and "3D" columns correspond to lattice configured for respectively longitudinalonly, longitudinal-horizontal, and three-dimensional cooling.
misalignments, with a corresponding reduction in the total cooling force as well as a modified transverse profile of the focused UR.The second was increased scattering due to residual gas in the IOTA ring.This would increase the equilibrium transverse beam sizes and single-particle amplitudes, exacerbating the transverse effects described in Section III D. Others included nonlinearities in the electron bypass and distorted trajectories in the undulators due to saturation in the steel poles.
Three-dimensional scans of the integrated experimental apparatus and correspondingly modified simulations of the UR suggest a reduction of the angular acceptance of the lens by ∼ 25% with a corresponding reduction of the maximum kick strength by ∼ 15%.The experiment also observed a weaker-than-expected coupling ratio for the nominal excitation of the QX1 quadrupole magnet, which suggests the presence of additional sources of coupling in the bypass.Reducing the strength of QX1 by half results in the observed experimental coupling ratio of 1.00:0.35.With this reduced coupling, the effects of OSC are seen most prominently in the longitudinal plane.To account for the increased residual gas scattering, we insert a scattering element in the lattice that simulates coulomb scattering events.The equilibrium transverse emittance in the experiment due to SR damping alone was measured to be approximately four times larger than the design value.The scattering strength was tuned to match this empirical value and the same simulations as before are run with these modifications.The results are outlined in Table V.
As expected, the overall rates are significantly reduced and the effect is seen mostly in the longitudinal plane.Subtracting the SR-damping rates, the total OSC damping rates in the 1D, 2D and 3D configurations are 18.7 s −1 , 17.2 s −1 and 17.7 s −1 , respectively.This is in good agreement with the experimentally measured value of τ −1 tot = 18.4 s −1 and suggests that exacerbated transverse effects were the principal source of reduced OSC force in the experiment.

Optical Delay and Equilibrium Distributions
The standard configuration of OSC is designed to reduce the deviations in the particles' positions and momenta, producing a beam with a lower 6D emittance.This relies on proper temporal alignment of the beam and UR in the kicker such that a reference particle experiences no net energy change.The optical delay system is responsible for establishing and maintaining the correct optical path length.In the IOTA experiment, the optical delay system is made up of two rotating glass plates which provide fine control over the optical delay [52].We implement this in the model using a phase term, ψ 0 , in Eq. 6.
The normal operation of OSC has a delay phase of ψ 0 = 0 while a delay phase of ψ 0 = π establishes a "heating mode" in which the small-amplitude motion is unstable and particles are driven away from the reference momentum.However, because of the periodic nature of the electromagnetic radiation, the particles are pushed toward higher-order attractors in phase space, resulting in ring-like structures in the different phase-space planes.
The IOTA experiment demonstrated the stability of the OSC bypass and the control over the optical delay by slowly sweeping the delay and observing the response of the beam.The optical delay is initially set such that δℓ ∼ +N u λ r , where there is no overlap between a particle and its corresponding radiation wave packet.The delay plates are then rotated slowly, at a rate of 0.01 deg•s −1 (∼ 0.03λ r •s −1 ) to a delay length of δℓ ∼ −N u λ r .This sweep shows the sequence of heating and cooling modes with the beam distribution reaching equilibrium at each point in the scan [28].The full envelope of the scan constitutes a measurement of the integrated system bandwidth, which was determined to be ∼ 20 THz.
In the case of passive OSC, a complete simulation of this scan was not feasible due to the slow scan rate required to ensure that the beam distribution reaches equilibrium for each delay setting [53]; therefore, we simulate such a scan for a single period δℓ ∈ [0, λ].Even so, this shorter scan still required the use of a faster scan rate than the experiment, and some resulting differences are apparent.Figure 4  tions, we used 200 macroparticles and executed a slow phase sweep through a single period of the OSC force, starting at δℓ = 0 and ending at δℓ = 1.25λ r , with a scanning speed of ≃ 1.39λ r s −1 (corresponding to a total phase shift of ∆ψ 0 ≃ 2π every 600 synchrotron periods).Due to computational restrictions, this sweep was ∼ 600 times faster than the experiment.As a result, the elegant data in Fig. 4(b) shows a lagging response to the optical-delay sweep compared to the measurement displayed in Fig. 4(a) as the beam could not fully equilibrate.However, the simulation qualitatively reproduces the experimental data with cooling occurring at ∼ 0 and δℓ ≃ λ r and maximum heating at δℓ ≃ 0.5λ r , albeit with a shift of ∼ 0.25λ r due to the rapid scan rate.
To further investigate the delay variable we simulated equilibrium beam distributions for delay settings of ψ 0 = 0 and ψ 0 = π, which correspond to the OSC cooling and heating modes respectively.The results are compared with the equilibrium distributions recorded in the experiment and are shown in Fig. 5.The transverse distributions were measured in elegant by recording particle coordinates in the center of the M2R bending magnet and binning over several thousand turns.The horizontal distributions [see Figs.The longitudinal distributions were recorded at the M3R dipole and compared with streak camera projections from the experiment.The longitudinal heating equilibrium distributions show very good agreement indicating that the computational model behaves as expected when operating in a heating mode.The simulated longitudinal distribution for the cooling case is slightly more narrow than the experiment.This is most likely due to the increased IBS in the cooling mode, as seen in the previous section.

D. Heating Mode Dynamics
In the nominal 2-D cooling configuration of OSC (shared cooling between the longitudinal and horizontal planes), the particles are cooled simultaneously in both planes towards the design orbit.However, as discussed in Section IV C 3, when the optical delay is set to ψ 0 = ±π, particles at low amplitudes are driven away from the design orbit towards stable, high-amplitude orbits in phase space.
In order to quantitatively describe this effect, it is useful to rewrite Eq. 6 to also include betatron oscillations.Following Ref. [54], Eq. 6 is parameterized in terms of normalized betatron and synchrotron amplitudes and phases, (a x , Ψ x ) and (a p , Ψ p ) respectively, as where Φ c is the phase shift between the momentum kick and betatron motion.Using this definition, integrating over the synchrotron and betatron oscillations reveals the presence of stable-orbit points (or attractors) in the (a p , a x )-parameter space [20].The synchrotron amplitude is where δ m is the maximum fractional-momentum deviation over a full synchrotron period.The normalized betatron amplitude is is the Courant-Snyder invariant, and (β x , α x ) are the usual Courant-Snyder parameters [55].In the heating mode (ψ 0 = ±π) for the 2-D configuration, the lowest-order attractors have high amplitude in one phase plane and simultaneously low amplitude in the other.This feature was observed in the elegant simulations by examining the evolution of the horizontal and longitudinal amplitudes associated with randomly sampled macroparticles within the bunch; see Fig. 6(a).The macroparticles evolve to one of the two heating attractors with either high-amplitude (blue dots) or low-amplitude (green dots) longitudinal motion.The bifurcation of the ensemble towards these attractors results in an equilibrium distribution with a "bullseye"-shaped distribution in the longitudinal and horizontal phase-space consisting of a bright core with a peripheral ring as shown respectively in Fig. 6(b) and (c).These bifurcated phase space distributions yield a multi-modal spatiotemporal (s, x) distribution as shown in Fig. 6(d).The ratio of macroparticles populating these regions of the (a p , a x ) parameter space depends on the strength of the longitudinal-tohorizontal coupling (1:0.34 in this case).For an imbalanced ratio, the particles will be preferentially driven to the stronger attractor, which results in cooling for the other phase plane.This effect was clearly observed in the OSC experiment [28].
In the experiment, the equilibrium spatiotemporal distribution was recorded at M3R for the 2-D cooling configuration and is shown in Fig. 7(a).Elegant simulations were performed for the same configuration and the macroparticle coordinates were recorded at the effective source plane of the M3R diagnostics.
The results of this simulation appear in Fig. 7(b,c).Figure 7(b) shows the bifurcation of the beam as particles move to one of the two high-amplitude attractors as already examined in Fig. 6(d).The simulated streak camera image confirms the presence of two lobes along the longitudinal axis corresponding to the highsynchrotron-amplitude, low-betatron-amplitude attractor (a x = 0, a p ≃ ν 11 ) where ν 11 ≃ 3.81 represents the first zero of the Bessel function of the first kind J 1 (x) [20].The two lobes along the ordinate of Fig. 7(b) correspond to the other attractor (a x ≃ ν 11 , a p = 0).The fine features observed in the simulated distribution shown in Fig. 7(b) are not present in the experimental data [Fig.7(a)] owing to the finite resolution and other limitations of the streak camera imaging system.To make a qualitative comparison, a Gaussian point-spread function was applied to the simulated distribution along the x dimension.This simulated "measurement" appears in Fig. 7(c) and agrees well with the experimentallymeasured distribution of Fig. 7(a) thus confirming the ability of the numerical model to accurately simulate the dynamics of heating modes as the beam evolves toward equilibrium.Where direct comparisons are possible, we find the simulated beam distributions to be in good quantitative agreement with the measured distributions.

V. ADDITIONAL FEATURES
The OSC model was developed to be a general tool, and it includes several features that could not be easily studied in the IOTA experiment.The initial experiment was designed to use beams with low bunch charge to re- duce the effects of IBS.This also has the effect of reducing the incoherent contributions to OSC as very few particles (on average N s ∈ [20,2000]) will populate each sampling slice.There are other features of the model which can be demonstrated with minimal additional setup such as amplified OSC and errors in angular alignment using multiple OSC PU-KU pairs.In this section, we present brief examples and analysis of these features.

A. Incoherent OSC
The incoherent contributions of neighboring particles have been modeled before using a statistical approach [50] but have not been experimentally studied.Our model computes this effect directly by applying a kick to the individual particle from every particle in its sample slice.Figure 8(a) shows the average kick each particle in a beam receives as a function of its average momentum deviation.The values are averaged over 50 turns, which corresponds to elapsed time of 0.67 ms (i.e.0.3% of a synchrotron period).The change in longitudinal position and momentum of a single particle is negligible over this duration, which provides a reasonable average kick for a given momentum deviation.We simulated this effect for three different beam currents.As the charge density increase, the incoherent effects become more dominant.Figure 8(b) shows the standard deviation of the incoherent kick that a particle will experience in a single turn.The variance σ 2 inc grows linearly with the longitudinal charge density.The average standard deviation scales with the number of turns N as ⟨σ inc ⟩ N = σ inc / √ N so this effect is negligible for low charge and long bunch length.

B. Amplified OSC
In future experiments, the UR produced in the pickup will be amplified before interacting with the particles in the kicker.By amplifying the PU radiation, the beam can be cooled much faster than in the passive case.In the model, amplification is included by simply increasing the kick strength parameter.The limited optical delay of the initial IOTA experiment could not accommodate an optical amplifier; however, for purposes of demonstration, Fig. 9 presents the simulated evolution of the longitudinal emittance ε t ≡ ε s /c as the beam is cooled by OSC; see Eq. 8.These calculations used the nominal uncoupled IOTA lattice and consider various opticalpower gains G ranging from 10 to 70 dB (corresponding to √ G ∈ [3,3000]).The simulations include incoherent effects but neglect IBS.When incoherent effects are negligible, the cooling rate increases linearly with the maximum energy kick (proportional to √ G).In the absence of size-dependent diffusion effects like IBS, the longitudinal emittance in equilibrium follows the relationship ε t,OSC /ε t,SR = τ OSC /τ SR [28].As the gain increases, the incoherent contributions limit the longitudinal emittance.Eventually, large gains yield an unstable particle distribution with increasing longitudinal emittance.The simulations indicate that an optimal gain for providing the lowest-emittance beam while maintaining equilibrium in the IOTA lattice is G ≃ 40 dB.

C. Multiple pairs of OSC elements
The OSC model discussed so far consists of a pair of lattice elements, the pickup and the kicker linked using an ID-string within elegant.Specifically, the results discussed in Section IV use a single OSC pair considered as thin elements.However, multiple pairs can be defined so as to introduce an OSC element for each undulator period (i.e. in a sequence of paired pickup and kicker elements).Such a configuration allows for the investigation of thick-lens effects in the OSC elements.
One such effect is angular misalignment in the PU and KU.Similar to the discussion in Section II A 3, a beam entering the kicker with at an angle will experience the off-axis electric field in the KU.This generally will reduce the cooling damping decrement.Figure 10 presents the evolution of the bunch length of a beam in the presence of OSC with increasing misalignment in the KU.The angular misalignment θ is introduced by assigning a transverse misalignment δx n = (s n − s KU ) tan θ to each kicker element (indexed by the n subscript) depending on its location within the KU.Angular misalignment decreases the cooling rate and results in a longer equilibrium bunch length.In the example shown in Fig. 10, corresponding to the IOTA passive-OSC case, a ∼ 2 fold increase in the equilibrium bunch length and a ∼ 30% reduction in cooling rate is observed for an angular misalignment θ = 500 µrad.

VI. CONCLUSION
We have developed a fast computational model of the OSC mechanism and have demonstrated its ability to accurately simulate the dynamics of the passive-OSC proof- of-principle experiment performed at the IOTA storage ring.We modeled the IOTA storage ring in elegant and observed good quantitative and qualitative agreement between the simulations and the measured performance of the OSC system for a variety of system configurations.
The current model includes the effects of the transverse profile of UR, optical delay, and incoherent contributions to OSC of neighboring particles.Other features of the model (incoherent kicks, effects of misalignment, and optical gain) were also showcased and will guide future OSC experiments while providing a microscopic understanding of OSC dynamics.Finally, although the current model implemented in elegant considers a transit-time OSC configuration, it could be easily extended to implement other stochastic cooling configurations and, more broadly, be used to investigate the control of beam distributions using self-field radiated at an earlier time.
FIG. 1. (a)A TTOSC system with pickup and kicker undulators, a bypass chicane, an optical line, and (b) the position of a particle and its own undulator radiation at the exit of the pickup undulator, and at the entrance and exit of the kicker undulator.The particle slips behind the radiation one wavelength each undulator period and the bypass introduces a delay so that the particle arrives at the kicker at the front of the radiation wavepacket.The delay can be tuned such that the subsequent energy exchange in the kicker produces corrective kicks that cool the beam.

FIG. 2 .
FIG. 2. Diagram of the IOTA storage ring at Fermilab (a) with a zoomed-in view of the OSC line (b) and horizontal and vertical betatron functions around the ring (c).In (a,b) blue, red, and greed shapes correspond respectively to dipole, quadrupole, and sextupole magnets.In (c), the shaded area corresponds to the location of the OSC section and the labels MiR (i = 1, 2, 3) indicate the location of the first three bending magnets [see (a)] where beam diagnostics are available.The beam circulates in the clockwise direction in the ring (i.e.SQ1→ M1R → M2R → M3R→ OSC → SQ1).

FIG. 3 .
FIG.3.Waterfall plot of the temporal evolution of longitudinal (a) horizontal (b) and vertical beam distribution.OSC is turned on instantaneously at t = 0 ms (indicated with the vertical dash lines).The "1D", "2D", and "3D" columns correspond to lattice configured for respectively longitudinalonly, longitudinal-horizontal, and three-dimensional cooling.

FIG. 4 .
FIG. 4. Waterfall plot showing the evolution of longitudinal projections measured (a) and simulated with elegant (b) as the optical delay varies within the range δℓ ∈ [0, 1.25λr].

FIG. 5 .
FIG. 5. Comparison of the experimental (blue) and simulated (green) equilibrium distributions in longitudinal (a,b), horizontal (c,d) and vertical (e,f) directions for the cooling (ψ0 = 0, left column) and heating (ψ0 = π, right column) modes.The shaded area in plots (c,d) represents the regions where the measurements are accurate (see text for details).The development of a tail toward negative values of x is a measurement artifact; see text for details.

FIG. 6 .
FIG.6.Evolution of macroparticles from their initial positions (black dots) in (ap-ax) amplitude space due to OSC heating in the 2D-cooling configuration (a).The blue dots labeled "1" represent macroparticles associated with the high longitudinal-amplitude attractor, while the green dots labeled "2" correspond to macroparticles in the high horizontalamplitude attractor.The equilibrium distribution is shown in longitudinal (b) and horizontal phase space (c).The spatiotemporal distribution (s, x) of particles is shown in (d), where the axes correspond to the spatial components of (b) and (c).The attractors "1" and "2" respectively correspond to (ax = 0, ap ≃ 3.81) and (ax ≃ 3.81, ap = 0) where 3.81 corresponds to the first zero of the J1 Bessel's function.

FIG. 7 .
FIG. 7. Streak camera measurement of the spatiotemporal distribution (s, x) at M3R for the 2D-cooling configuration (s − x coupling) operated in the heating mode (a) and corresponding distribution simulated with elegant (b) along with the simulated measurement obtained by smearing the simulated distribution with a Gaussian point-spread function along the x axis.The labels "1" and "2" in (b) refer to the same attractors as defined in Fig. 6.

FIG. 8 .
FIG.8.The energy kick ∆E ≡ ∆E + ∆ Ẽ, including the coherent and incoherent contributions, each particle receives relative to the maximum kick K (a) and the standard deviation of the incoherent contribution for a single turn normalized to the maximum kick strength(b).

FIG. 9 .
FIG. 9. Effect of amplification on OSC damping for gains G of 10-70 dB.The simulation used an amplitude gain, √ G.

TABLE I .
Nominal undulator and beam parameters for the passive TTOSC proof-of-principle experiment at the IOTA facility.

TABLE II .
[43]arison of lattice parameters between the nominal ("theory") values[43]and the simulated values obtained from the elegant model.These calculations do not consider the OSC process.
† Computed in elegant using lattice optics

TABLE III .
[43]arison of synchrotron-radiation (SR) damping and equilibrium-distribution parameters analytically computed ("theory")[43]with simulated values obtained from the elegant model ("elegant") .These calculations do not consider the OSC process or lattice coupling.The emittance and bunch length values are RMS values.

TABLE IV .
Total damping rates simulated with elegant for the IOTA design configuration.

TABLE V .
Damping rates simulated with elegant for the IOTA experimental ("as built") configuration.