Spin dynamics investigations for the electric dipole moment experiment

Precision experiments, such as the search for a deuteron electric dipole moment using storage rings like COSY, demand for an understanding of the spin dynamics with unprecedented accuracy. In such an enterprise, numerical predictions play a crucial role for the development and later application of spin-tracking algorithms. Various measurement concepts involving polarization effects induced by an rf Wien filter and static solenoids in COSYare discussed. The matrix formalism, applied here, deals solely with spin rotations on the closed orbit of the machine, and is intended to provide numerical guidance for the development of beam and spin-tracking codes for rings that employ realistic descriptions of the electric and magnetic bending and focusing elements, solenoids, etc., and a realistically modeled rf Wien filter.


I. INTRODUCTION
The Standard Model (SM) of particle physics is not capable to account for the apparent matter-antimatter asymmetry of the Universe. Physics beyond the SM is required and it is either probed by employing high energies (e.g., at LHC), or by striving for ultimate precision and sensitivity (e.g., in the search for electric dipole moments). Permanent electric dipole moments (EDMs) of particles violate both time-reversal ðT Þ and parity ðPÞ invariance, and are via the CPT theorem also CP violating. Finding an EDM would be a strong indication for physics beyond the SM, and pushing upper limits further provides crucial tests for any corresponding theoretical model, e.g., SUSY.
Up to now, EDM searches mostly focused on neutral systems (neutrons, atoms, and molecules). Storage rings, however, offer the possibility to measure EDMs of charged particles by observing the influence of the EDM on the spin motion in the ring. These direct searches of e.g., proton and deuteron EDMs bear the potential to reach sensitivities beyond 10 −29 e cm. Since the cooler synchrotron COSY 1 at the Forschungszentrum Jülich provides polarized protons and deuterons up to momenta of 3.7 GeV=c, it constitutes an ideal testing ground and a starting point for such an experimental program.
The investigations presented here, carried out in the framework of the JEDI collaboration, 2 are relevant for the preparation of the deuteron EDM measurement [3]. A radio-frequency (rf) Wien filter (WF) [4][5][6] makes it possible to carry out EDM measurements in a conventional magnetic machine like COSY. The idea is to look for an EDM-driven resonant rotation of the deuteron spins from the horizontal to the vertical direction and vice versa, generated by the rf Wien filter at the spin precession frequency [7,8].
The search for EDMs of protons, deuterons, and heavier nuclei using storage rings [2,9,10] is part of an extensive world-wide effort to push further the frontiers of precision spin dynamics of polarized particles in storage rings. In this context, the JEDI results prompted the formation of the new Charged Particle Electric Dipole Moment (CPEDM) collaboration, 3 which aims at the development of a purely electric prototype storage ring, with drastically enhanced sensitivities to the EDM of protons and deuterons, compared to what is presently feasible at COSY [3,11].
Precision experiments, such as the EDM searches, demand for an understanding of the spin dynamics with unprecedented accuracy, keeping in mind that the ultimate aim is to measure EDMs with a sensitivity up to 15 orders in magnitude better than the magnetic dipole moment (MDM) of the stored particles.
The description of the physics of the applied approach, called rf Wien filter mapping, is discussed further in a forthcoming separate publication. The theoretical understanding of the method and its experimental exploitation are prerequisites for the planned EDM experiments at COSY [2], and will also have an impact on the design of future dedicated EDM storage rings [11].
This paper discusses various polarization effects that are induced by the rf Wien filter and static solenoids in the ring. The approach taken here strongly simplifies the machine lattice, and deals solely with spin rotations on the closed orbit [12,13], described by the SOð3Þ formalism. One aim of the work is to obtain a basic understanding about the interplay of spin rotations in a magnetic ring equipped with an rf Wien filter and solenoid magnets, under the simplifying assumption mentioned above. In an ideal machine with perfect alignment of the magnetic elements, the spin rotations on the closed orbit are generated primarily by the dipole magnets, therefore, for the time being, spin rotations in the quadrupole magnets are not considered.
As we shall demonstrate below, even with an idealized ring, the parametric rf resonance-driven spin rotations reveal quite a reach pattern of spin dynamics. Our results set the background for more realistic spin tracking calculations, based on recent geodetic surveys of COSY that make available position offsets, roll, and inclination parameters for the quadrupole and dipole magnets. The treatment of the spin transport through these individually misaligned magnetic elements, can, however, be readily incorporated in the applied matrix formalism. Besides that, the spin dynamics simulations carried out in the framework of the present paper will serve as a valuable cross-check of the analytic approximate treatment of the parametric spin resonance, based on the Bogolyubov-Krylov-Mitropolsky averaging technique [14].
The JEDI collaboration is presently implementing a beam-based alignment scheme at COSY, which aims at providing optimized beam-transfer properties of the quadrupole and dipole magnets in the ring, intending to make the beam orbit as planar as possible [15]. Once this is accomplished, the spin dynamics in the ring will be largely governed by the misaligned dipoles alone. Thus effectively, the approach described here will appropriately describe an EDM experiment using an rf Wien filter in a beam-based aligned ring.
The paper is organized as follows. In Sec. II, the effect of an EDM on the spin evolution in a ring is discussed in terms of the Thomas-BMT equation [12]. The inclusion of an rf Wien filter in an otherwise ideal ring is treated in Sec. III, while the polarization evolution with an rf Wien filter and additional solenoids is discussed in Sec. IV. The main findings are summarized in the conclusions in Sec. V. A brief outlook into additional aspects planned to be investigated using the simulation approach taken here in the near future is also given.

II. SPIN ROTATIONS IN THE RING A. Thomas-BMT equation
Below, the basic formalism to describe the spin evolution in electric and magnetic fields is briefly reiterated. The generalized form of the Thomas-BMT equation describes the spin motion of a particle with spin ⃗ S in an arbitrary electric ( ⃗ E) and magnetic field ( ⃗ B). Including EDMs (in SI units), it reads [16] where Here m, γ, and ⃗ β are the mass, Lorentz factor, and the velocity of a particle in units of the speed of light c in vacuum. ⃗ S (in units of ℏ) is given in the particle rest frame, and the fields ⃗ E and ⃗ B are in the laboratory system. The magnetic dipole moment ⃗ μ (MDM) and the electric dipole moment ⃗ d (EDM) are defined via the dimensionless Landéfactor g and η EDM The magnetic anomaly is given by B. EDM tilt angle ξ from the Thomas-BMT equation In an ideal machine without unwanted magnetic fields, the axis about which the particle spins precess is given by the purely vertical magnetic field ⃗ B¼ ⃗ B ⊥ ¼B ⊥ · ⃗e y . Equating the COSY orbital angular velocity Ω rev ¼ 2πf rev and the relativistic cyclotron angular velocity yields, for ⃗ E ¼ 0 with the parameters given in Table I, a vertical magnetic field of F. RATHMANN, N. N. NIKOLAEV, and J. SLIM PHYS. REV. ACCEL. BEAMS 23, 024601 (2020) This value can be considered as the magnetic field that corresponds to an equivalent COSY ring where the magnetic fields are evenly distributed. Inserting ⃗ B from Eq. (6) and ⃗ E ¼ 0 into Eq. (2), yields for the angular velocity in the particle rest system In the laboratory system, however, we observe with the parameters of Table I the angular velocity with respect to the cyclotron motion of the momentum, where ⃗ Ω rev denotes the COSY angular velocity along ⃗ e y . The spin-precession frequency yields the familiar value of which is also listed in Table I. The angle by which the stable spin axis is tilted, i.e., the angle between ⃗ Ω Lab and ⃗ e y is obtained by evaluating Inspecting Eq. (8), the effect of an EDM in a magnetic machine can be expressed by the tilt of the stable spin axis away from the vertical orientation in the ring, given by 4 For an assumed EDM of d ¼ 1 × 10 −20 e cm, and for deuterons at a momentum of 970 MeV=c, Eqs. (3) and (11) yield η EDM and ξ EDM , as listed in Table I.   TABLE I. Parameters of the deuteron kinematics, the COSY ring, the deuteron elementary quantities, the electric dipole moment (EDM) assumed, and the field integrals of the idealized rf Wien filter (to eight decimal places). The deuteron momentum P is used to specify the deuteron kinetic energy T, and the Lorentz factors β and γ. The COSY circumference l COSY is used to specify the COSY revolution frequency f rev and the spin-precession frequency f s . The deuteron mass m and the deuteron g factor, taken from the NIST database [17] (not from the most recent one), are used to specify G. The deuteron EDM d is used to quantify η EDM and ξ EDM .

C. Rotation matrices
Our description of the spin dynamics is based on the SOð3Þ formalism. A rotation by an angle θ around an arbitrary axis given by the unit vector ⃗ n ¼ ðn 1 ; n 2 ; n 3 Þ is described by the matrix [19] Rð⃗ n; θÞ ¼

D. One turn spin rotation matrix with EDM
With a nonvanishing EDM, in the rotation matrix of Eq. (12), the spins do not precess anymore around the vertical axis ⃗ e y , but rather around the direction given by Therefore, the ring rotation matrix can be obtained by inserting into Eq. (12) the coefficients c 1 , c 2 , c 3 from Eq. (14), and by setting Here, the time t is defined by the number of momentum revolutions n in the ring, The spin-precession frequency f s , related to ⃗ Ω Lab introduced in Eq. (8), can be expressed also via where f rev denotes the revolution frequency, and ν s the spin tune, i.e., the number of spin precessions per turn in the ring. A negative G factor indicates that the precession proceeds opposite to the orbit revolution. Thus, a one-turn matrix for the ring including the EDM effect is obtained by inserting θðtÞ from Eq. (15) into Eq. (12) at t ¼ T rev ¼ 1=f rev . For comparison with numerical simulations, the ring matrix is explicitly given below (to four decimal places) using the parameters listed in Table I

E. Polarization evolution in the ring
The evolution of the polarization vector ⃗ S 1 as a function of time in the ideal bare ring is then described by where ⃗ S 0 denotes the initial polarization vector. Throughout the present paper, the single-particle spin evolution is described by unitary transformations which preserve the magnitude of the polarization. Figure 1 shows the situation when the spin rotation axis ⃗ c, defined by Eq. (14), is tilted with respect to the normal to the ring plane ⃗ n (y axis in the figure). 5 In Fig. 2 for ten turns. It is visible that the polarization evolution occurs counterclockwise with respect to the clockwise rotation of the particles in the ring, since the deuteron G factor is negative.

A. Electric and magnetic fields of the rf Wien filter
The rf Wien filter, described in [4], has been designed in order to be able to manipulate the spins of the stored particles, avoiding as much as possible effects on the beam orbit. To this end, great care was taken to minimize the unwanted field components of the Wien filter and to characterize them via the polynomial chaos expansion [5]. In the EDM mode, the main component of the magnetic induction ⃗ B WF is oriented along the y axis, and the main component of the electric field ⃗ E WF along the x axis. In order to avoid betatron oscillations in the beam, the magnetic and electric fields must be matched to each other to provide a vanishing Lorentz force ⃗ F L [see Eq.
(3) of [4] ], According to a full-wave simulation (FWS), 6 including the ferrite cage (see label 6 in Fig. 1 of [4]), for an input power of 1 kW, a field integral of ⃗ B WF along the beam axis of is obtained. Here, the active length of the rf Wien filter [4], denoted by is defined as the region, where the fields are nonzero. Under these conditions, the corresponding integrated electric field components with ferrites are The design and construction of the rf Wien filter includes a ferrite cage surrounding the electrodes, which improves the field homogeneity and increases the magnitude of the fields [4]. However, in order to simplify the installation, the rf Wien filter was installed at COSY without ferrites. In addition, it was decided to proceed without ferrites until a first direct deuteron EDM measurement is available.
For this situation without ferrites, and for an input power of 1 kW [ignoring the unwanted components of the field integrals (B WF x , B WF z , and E WF y , E WF z )], one obtains from the full-wave simulation (FWS) The ratio of electric and magnetic field integrals from the FWS yields FIG. 2. Polarization evolution during idle precession for ten turns in an ideal ring using Eq. (19) and the parameters listed in Table I. Panel (a) shows p x ðtÞ, p z ðtÞ and p y ðtÞ for an initially longitudinal polarization, and panel (b) the same for initial sideways polarization. The bunch revolution is indicated as well. The magnitude of the p y oscillation amplitude is equal to the tilt angle ξ EDM [see also Eq. (14) and Fig. 1].
and should ideally be equal to unity. The subsequent calculations use the field integrals of an idealized Wien filter with vanishing Lorentz force ⃗ F L , given in the last column of Table II. A field amplification factor is applied in the simulations to increase the field integrals of the ideal rf Wien filter (last column Table II) in the simulations, so that The field amplification allows one to speed up the simulation calculations accordingly, without affecting other aspects of the spin dynamics of the polarization evolution in the ring. The reason is that the involved spin rotations scale with the field integrals applied by the rf Wien filter. It should be noted that if instead the magnitude of the EDM would be scaled up, the spin dynamics would be affected because ξ EDM becomes larger (see Fig. 1). In the description of the spin evolution via spin rotations on the closed orbit, momentum and position kicks are not considered.

B. Rotations induced by the rf Wien filter
The effect of the rf Wien filter on the polarization evolution in the ring is implemented by an additional rotation matrix. The spin rotation in the Wien filter depends on the applied field integrals (right column of Table II), multiplied by the factor f ampl .

Spin rotation angle in the Wien filter
In the following, the spin rotation angle ψ WF in the rf Wien filter is calculated numerically using the Thomas-BMT equation of Eqs. (1) and (2) with ⃗ Ω EDM ¼ 0. We start with an initial spin vector and we compute the final polarization vector ⃗ S fin via Electric and magnetic field vectors for ⃗ Ω MDM in Eq.
(2) are obtained by computing the average fields from the idealized field integrals of the rf Wien filter (last column of Table II), given by where the effective length of the Wien filter is taken from Eq. (22). These conditions provide for a vanishing Lorentz force ⃗ F L [see also Eq. (20)]. After passing the rf Wien filter once, the final polarization vector is given by and, after normalizing ⃗ S fin to unity, the angle between S in and ⃗ S fin is determined from the four-quadrant inverse tangent using the parameters listed in Table I.   TABLE II. Values for the main electric and magnetic field integrals from the full wave simulation with and without ferrites for an input power of 1 kW where ⃗ B WF k⃗ e y . The last column lists the electric and magnetic field integrals of an idealized Wien filter used in the simulations. In this case, the unwanted field components vanish, i.e., The spin-rotation angle in the rf Wien filter, divided by the idealized transverse magnetic field integral from Table II, yields Validating the numerical result for the spin-rotation angle ψ WF in the rf Wien filter obtained in Eq. (32) against the analytic expression, given in Eq. (13) of [18], yields where the time interval Δt in the Wien filter has been expressed through the length l WF . The spin rotation angle in the rf Wien filter, given in Eq. (34), constitutes an upper limit, which corresponds to a situation when a sharp δ-function-like bunch passes through the device. Realistically, the bunch distribution has to be folded in, and the average spin-rotation angle will be reduced correspondingly.

rf Wien filter rotation matrix
The spin-rotation angle of the rf Wien filter changes as a function of time according to where The parametric spin resonance condition for the Wien filter frequency is given by the sum of a harmonic multiple (K) of the revolution frequency f rev and the spin-precession frequency f s [Eq. (17)], The rf Wien filter was designed to allow for the operation at different frequencies f WF [4]. Figure 3 shows the available harmonics (closed circles) for protons (at T ¼ 135 MeV) and for deuterons (at P ¼ 970 MeV=c).
The rf Wien filter rotation matrix is given by where in the generic case, ⃗ n WF is a unit vector along the magnetic field of the Wien filter. The case for instance, denotes the Wien filter EDM mode. The rf Wien filter matrix U WF ðtÞ is only evaluated once per turn when the condition is met stroboscopically, otherwise, the implemented function returns the I 3 unit matrix. When the Wien filter is rotated around the beam axis (z) by some angle ϕ WF rot , and the oscillations also receive a contribution from the rotation of the MDM in the horizontal magnetic field.

C. Polarization evolution in the ring with rf Wien filter
The evolution of the polarization vector ⃗ S as a function of time t in the ring with rf Wien filter can be numerically evaluated via The corresponding situation of the ring is illustrated in Fig. 4. The spin rotations in the ring can be described by U ring . A turn begins with the revolution in the ring at t ¼ 0, T rev ; …, n · T rev , and it ends with one pass through the rf Wien filter. Between two successive points in time at which a particle encounters the rf Wien filter, its spin is just idly precessing in the machine. Sequential unitary rotations in Eq. (42) manifestly preserve the magnitude of the polarization. According to Eq. (42), the spin motion is stroboscopic in the sense that the spin rotation follows the angle ψðtÞ of the rf Wien filter [Eq. (35)] turn by turn. The rf Wien filter therefore induces a stroboscopic turn-by-turn conversion of the transverse in-plane polarization into a vertical one (or vice versa). In the spirit of Eq. (30), Eq. (42) can be cast in the form of a system of finite difference equations. In the continuous limit, one obtains a system of homogeneous first order differential equations with time-dependent coefficients which are periodic functions of ω WF · t. The latter property is imposed by the time-dependent spin rotation in the Wien Filter, see Eq. (35).
The transition to the continuous limit is furnished by the Bogolyubov-Krylov-Mitropolsky (BKM) averaging [14], which approximates the turn-by-turn evolution by a continuous dependence on the revolution number, given by n ¼ f rev · t [Eq. (16)]. Then, Eq. (37) provides the condition for exact parametric resonance [20][21][22] that generates the up-down oscillation of the polarization. For the generic orientation of the rf Wien filter, the BKM averaged buildup of the vertical polarization proceeds with the resonance tune (or strength) [18] The above formula is universal for rf spin rotators of all kinds, rf Wien filters, electric and magnetic rf dipoles, and rf solenoids.
In the EDM mode [see Eq. (39)], the above equation yields using the parameters of Table I. Therefore, in the absence of other perturbing spin rotations in the ring, a measurement of the resonance tune constitutes a direct measurement of the EDM. 7 The generic case when the stable spin axis ⃗ c acquires an additional tilt from in-plane magnetic fields will be discussed in Sec. IV.
The direct simulations using Eq. (42), discussed below, will furnish important cross-checks with respect to the accuracy of the analytic approximations based on the BKM averaging.
D. Radial magnetic rf field in the Wien filter 1. Driven oscillations and resonance strength ε MDM Driven oscillations of the vertical polarization p y ðtÞ can also be induced by the horizontal magnetic field of the rf Wien filter that couples to the deuteron MDM. As an illustration of the principal features of the polarization evolution under exact parametric resonance [see Eq. (37)], we discuss below the situation where the rf Wien filter is operated in the so-called MDM mode with magnetic field along −⃗ e x , i.e., for ϕ WF rot ¼ 90°, and where the initial polarization ⃗ In this case ⃗ n WF ¼ −⃗ e x , and Eq. (43) predicts the resonance strength where ψ WF from Eq. (32) was used, and ξ EDM from Table I. Using the function for ⃗ S 2 ðtÞ, given in Eq. (42), for the conditions of Table I, driven oscillations for the rf Wien filter with magnetic field aligned along −⃗ e x [see Eq. (41)] were simulated. One example for K ¼ −1 is shown in Fig. 5. Subsequently, the simulated oscillations were fitted using the function The quality of the fit to the numerical data is evaluated in terms of squared deviations via where the weight factors are w i ¼ 1, and p y ðtÞ ¼ ⃗ e y · ⃗ S 2 ðtÞ. In the last row of Table III, the reduced χ 2 ¼ SSE=ndf is given, where n points ¼ 101, and ndf ¼ n points − 4 ¼ 97, since the fitted function in Eq. (46) has four parameters.
The fit results are summarized in Table III. The angular velocity Ω driven ¼ b was obtained using the field integrals, listed in the right column of Furthermore, the angular velocity normalized to the real magnetic field integral yields The induced driven oscillations, shown in Fig. 5, correspond to a resonance strength of where the factor f ampl in the denominator corrects for the field enhancement used in the simulation. The resulting value agrees nicely within errors with jε EDM j of Eq. (45).

Width of the spin resonance
The detuning of the frequency f WF at which the rf Wien filter is operated away from the spin-precession frequency f s [see Eqs. (17) and (37)] can be parametrized by substituting in Eq. (36) From about 50 simulations similar to the one shown in Fig. 6, the oscillation amplitudes and the oscillation frequencies as function of Δf WF are obtained by fitting, invoking again the parametrization given in Eq. (46). In order to reduce the time required for the simulations, a field amplification factor of f ampl ¼ 10 3 was used, which leads to oscillations that are faster by the same factor. The results, shown in Fig. 7, are corrected for the field amplification factor employed in the simulations. The dependence of the oscillation amplitude a of the simulated data can be described by a Lorentz curve (Breit-Wigner function) of the form FIG. 5. Simulated driven oscillations on resonance using ⃗ S 2 ðtÞ from Eq. (42) with initial vertical polarization for the parameters given in Table I and for the harmonic K ¼ −1. The plot contains 101 points for a total of 10 000 turns. where The normalization constant and the width of the Breit-Wigner, obtained from the fit shown in the left panel of Fig. 6, amount to h ¼ 1.000 00 AE 0.000 09 and For other harmonic excitations K ¼ 0, 1, and AE2 used in the rf Wien filter, within the errors the simulations yield the same results as given above. It should be emphasized that the off-resonance parametric modulation of the spin motion exhibits a profound difference compared to the more familiar oscillations that are driven by an external harmonic force [20][21][22]. In parametric resonances, detuning causes a continuously growing phase difference between rf phase and spin-vector phase. The accumulation of the vertical polarization comes to an end as soon as Δϕ WF ∼ π=2. Evidently, the exact pattern of the buildup of the vertical polarization will depend on the relative magnitude of Δf WF and the buildup frequency f driven ¼ ϵ EDM f rev , as evidenced by comparing Figs. 5 and 6. A full-fledged discussion of the Fourier spectrum of the off-resonance polarization evolution will be reported elsewhere. Here, only the gross features will be illustrated using the previous example when the rf Wien filter is operated in MDM mode, The obtained fit results of the simulated oscillations Ω driven ðΔf WF Þ ¼ b [Eq. (46)] are summarized in the panels on the right in Fig. 7. Ω driven exhibits a linear dependence on Δf WF at large jΔf WF j (strong detuning) and a parabolic dependence for weak detuning. The linear fits to the six outermost points, indicated in the top right panel in Fig. 7, show that the transition between the two regimes occurs at about Δf WF ≈ AE0.2 Hz.
In the following, the off-resonance results of Fig. 7 are subjected to an interpretation in terms of the BKM averaging. For an exact parametric resonance, the BKM averaging precisely predicts the harmonic evolution ⃗ S y ðtÞ ¼ −⃗ e y cos ð2πϵ EDM f rev tÞ; ð56Þ with unity amplitude. In the off-resonance case, the BKM averaging dictates the substitution The above substitution, together with Eq. (50), leads to For the limit of strong detuning, Δf WF ≫ f driven , one obtains In the above expression, the higher harmonics cos ð2πNΔ WF tÞ are suppressed by powers of the small parameter ðf driven =Δf WF Þ 2N and were omitted. The gross features of S y ðtÞ can then be approximated by the harmonic expansion of Eq. (46), yielding The linear behavior of Ω driven at large Δf WF , governed by Eq. (61) and depicted in the right top panel of Fig. 7, as expected yields for both branches a slope parameter near AE2π, d 1 ¼ −6.07 AE 0.02; and The quadratic fit to Ω driven ðΔf WF Þ in the bottom right panel of Fig. 7 yields which is consistent with the result for the same quantity, given in Eq. (48).
Equating the asymptotic behavior of Eq. (60) to the one of the Breit-Wigner parametrization of Eq. (52), one obtains which agrees well with the fit result given in Eq. (54). At weak detuning, the argument of the cosine function in Eq. (58) initially rises linearly with t, then reaches a maximum and starts to decrease around At still larger t, the argument is a sine function of time with the period 4T BKM , i.e., it exhibits the beating pattern, typical of a spin echo. In our simulations, the spinevolution time is limited, t ≤ T max ¼ 0.013 s. Evidently, the finite T max introduces a new frequency scale of With decreasing Δf WF , we expect a change of the dependence of Ω driven as a function of Δf WF as soon as the argument of the cosine in Eq. (58) In this situation (EDM mode), the experimental determination of the resonance strength ε EDM amounts to the determination of the tilt angle ξ EDM and of the associated EDM, via Eqs. (11) and (3).

Polarization evolution with development of p y ðtÞ
In the following, the polarization buildup in the machine is addressed. The interplay of the different frequencies involved is illustrated in Fig. 8.
The same situation as in Fig. 8 is depicted in Fig. 9, the only difference is the larger turn number. The graph illustrates the experimental evidence for an EDM, namely a nonvanishing slope of the vertical polarization component p y ðtÞ. This slope describes the steady out-of-plane rotation of the polarization vector on the background of oscillations shown in the bottom panels of Fig. 2.
The slope can be determined by fitting using where f s is not a fit parameter, but taken from Eq. (17). The oscillation amplitude A in Fig. 9 perfectly matches the angle ξ EDM , used in the simulation (see Table I). Using the above parametrization, the initial slope is given by 2. p y (t) dependence on the phases ϕ rf and ϕ S x 0 The rf phase ϕ rf is introduced in Eq. (35). During a real experiment, this phase needs to be maintained by a phase-locking system (for details see [23]). Another way to parametrize the same effect is via the angle ϕ S x 0 ¼∠ð ⃗ S 0 ; ⃗e x Þ, as illustrated in Fig. 10(a), keeping ϕ rf constant.
Within the formalism described in [18], it is the interplay between the stable spin axis ⃗ c at the rf Wien filter and its magnetic axis ⃗ n WF (k ⃗ B WF Þ that controls via ½⃗ c × ⃗ n WF the dependence on the orientation of ⃗ S 0 . On the other hand, one could start by fixing the orientation of ⃗ S 0 by picking some angle ϕ S x 0 , which amounts to shifting the spin phase while keeping the rf phase fixed. The resulting evolution of p y ðtÞ, however, must be the same, except for a possible constant shift between the two phases ϕ rf and ϕ S x 0 . The buildup of a vertical polarization component, which is equivalent to a rotation of the polarization vector out of the ring plane due to the EDM for a set of random azimuthal angles ϕ S x 0 and ϕ rf has been computed. The results are shown in Fig. 11. The fit results are listed in Table IV. Within the given uncertainties, the two simulated data sets for ϕ S x 0 and ϕ rf , as expected, yield the same results. The only difference is a phase shift of π=2 between fðϕ S x 0 Þ and gðϕ rf Þ, as evidenced by the difference of the parameters b 1 − b 2 from Table IV. The weights that are used to find the optimum parameters are all equal for each point of the two data sets. Correcting the initial slope parameter a in Table IV for the employed field amplification factor used in the simulation yields a prediction for the initial slope that one would expect in an ideal ring in the presence of an EDM of   Table I. The red line is a fit to the data using Eq. (72) that yields an initial slope of dp y ðtÞ=dtj t¼0 ¼ B ¼ ð4305.059 AE 5.268Þ × 10 −6 s −1 (for f ampl ¼ 10 3 ). the parameters for the idealized rf Wien filter, given in the last column of Table II, one obtains The comparison of _ p y ðtÞj t¼0 with experiment requires knowledge about the magnitude of ⃗ SðtÞ. The approach taken in [24] appears convenient, because the out-of-plane rotation angle α is independent of the magnitude of the beam polarization. The quantity of interest, indicated in Fig. 10(b), in that case is _ αðtÞj t¼0 . The polarimeter measures p y ðtÞ, irrespective of the in-plane polarization p xz ðtÞ, given by From this it follows that 3. Initial slope versus slow oscillation Figure 12(a) shows the initial slopes for four different assumed EDMs, for an ideal ring and an idealized Wien filter, based on the conditions listed in Table I. The EDMs manifest themselves twofold, namely in different slopes and in larger amplitudes of the fast oscillation. The linear slopes in Fig. 12(a) of course reflect just the very beginning of a sinusoidal oscillation that becomes visible only when the EDM is large. Such a situation is depicted in Fig. 12(b), where has been used in the simulation. The initial slope of the vertical polarization component is related to the strength of the EDM spin resonance. Another way to obtain this information is to vary the rf phase ϕ rf , as indicated in Fig. 11. The initial slope can of course also be obtained from the slow oscillation. The slope can be described by   Table IV. Each data point is obtained from a graph like the one shown in Fig. 9, but for 10 000 turns and 501 points. TABLE IV. Summary of parameters obtained (for K ¼ −1) via fitting the oscillatory patterns of the initial slopes shown in Fig. 11 as function of ϕ S x 0 and ϕ rf , still including the factor f ampl ¼ 10 3 . For the other harmonics (K ¼ 0, 1, and AE2), within the given uncertainties, the same values are obtained.
which respects the property that for any ϕ rf , p y ðtÞj t¼0 ¼ 0.
The derivative of p y ðtÞ with respect to time is where the value given corresponds to the situation shown in Fig. 12(b). Numerically, the red curve in Fig. 12(b) has been parametrized by the function fðtÞ ¼ p y ðtÞ ¼ a sinðω · t þ ϕÞ: The amplitude of the averaged oscillation [red curve in Fig. 12 The envelope b osc ðtÞ of the fast oscillations is consistent with the law b osc ðtÞ ¼ sin½ξ EDM ðdÞ · cosðωtÞ: ð83Þ According to [18], the EDM induced angular velocity ω in Eq. (78) can be expressed through the EDM resonance strength ε EDM and the orbital angular velocity frequency ω rev , via Now we can apply Eq. (84) to interpret the result for the initial slope for d ¼ 10 −20 e cm, given in Eq. (74). Equations (80) and (81) entail _ p y ðtÞj t¼0 ¼ cos ξ EDM · jε EDM j · ω rev ¼ 4.305 033 × 10 −6 ; ð85Þ where ε EDM from Eq. (44) was used with ξ EDM from Table I. Nice agreement is obtained with the value of the initial slope from the simulations [Eq. (74)].
In terms of the initial slope, the resonance strength is given by While the slopes can be easily determined as a function of ϕ rf , the latter method using Eq. (86) clearly also requires knowledge about the oscillation amplitude a. Knowing the initial slopes alone does not allow one to determine the resonance strength ε EDM .
Using the technique of variation of ϕ rf , as shown in Fig. 11, Fig. 13 yields an initial slope of

Determination of the running spin tune, based on the polarization evolution ⃗ S 2 ðtÞ
The standard definition of the spin tune as a rotation around the local stable spin axis ⃗ n s at every point in the machine holds for a static machine. It does not involve a time dependence of the polarization evolution, like the one generated by the rf Wien filter. If the polarization is time dependent, the term running or instantaneous spin tune will be used in the following. In case there is a timedependent or instantaneous spin tune, the direction of ⃗ n s also changes as a function of time, i.e., ⃗ n s ≡ ⃗ n s ðtÞ (see further Sec. III E 5).
Using the numerical simulations for ⃗ S 2 ðtÞ, or any other spin-evolution function, one can numerically determine the running spin tune in the following way. For this one needs three spin vectors from the spin-evolution function, say Using these three vectors, two more vectors are constructed, The in-plane angle between ⃗ dðtÞ and ⃗ eðtÞ can be used to determine the running, time-dependent spin tune ν s ðtÞ. To this end, we define the normal vector ⃗ N of the plane that contains ⃗ d and ⃗ e, which corresponds to the instantaneous spin axis. Using ⃗ N, we find the in-plane components of ⃗ b and ⃗ c, via The normalized versions of these vectors are called and the running spin tune is determined from The factors in front of arctangent take care that ν s ðtÞ generates the correct sign based on the G-factor and the number of spin precessions per turn. As a cross-check of the algorithm, with the rf Wien filter switched off, for the beam conditions given in Table I which is very close to the achievable machine precision (see footnote 8). During a revolution in the machine, as prescribed by ⃗ S 2 ðtÞ using Eq. (42), the spin tune remains constant during each turn [see Fig. 14(a)]. When the rf Wien filter is switched on, due to the additional spin rotation in the timevarying rf field, the instantaneous spin tune jumps from turn to turn. As depicted in Fig. 14(b), the oscillation amplitude of the spin tune variation due to the rf Wien filter using a power of 1 kW (see Table I) is well consistent with the expectation from the spin rotation formalism: . Each data point is obtained from a graph like the one shown in Fig. 9, but for ten turns and 1001 points.
The average spin tune, however, remains constant.

Instantaneous spin orbit determination based on ⃗ S 2 ðtÞ
The running spin orbit vector ⃗ n s can be easily determined from the procedure described in the previous section, using the normal vector ⃗ N, defined in Eq. (90), Similarly to the running (instantaneous) spin tune, the instantaneous spin orbit exhibits oscillating in-plane polarization components.

A. Evolution equation with additional static solenoids
In the course of this paper, with the rf Wien filter in EDM mode ( ⃗ B WF k⃗ e y ), the EDM interaction with the motional electric field in the ring was the only source of up-down spin oscillations.
In the following, two static solenoids in the straight sections will be added to the ring. Besides that, we shall make an allowance for rotations of the rf Wien filter around the longitudinal ⃗ e z (momentum) direction. Such rotations induce a radial magnetic rf field, and, in conjunction with the solenoidal magnetic fields, we start mixing the EDM and MDM induced rotations. The idea, common to all EDM experiments, is to disentangle the EDM signal using an extrapolation to a vanishing MDM contribution [25,26].
With two static solenoids added to the ring, the resulting sequence of elements is depicted in Fig. 15. The one-turn ring matrix can be split into two arcs, one arc made of the dipole magnets D 1 to D 12 , and the second arc made of dipoles D 13 to D 24 . Since U ring ð⃗ c; T rev Þ ¼ U arc 2 ring ð⃗ c; T rev =2Þ × U arc 1 ring ð⃗ c; T rev =2Þ; ð98Þ the two additional solenoids can be inserted before and behind arc 2, leading to invoking again the generic rotation matrix Rð⃗ e z ; χ rot Þ from Eq. (12).
In a similar fashion as in Eq. (42), one can then write for the polarization evolution,  Table I s for each turn. It should be noted that the initial spin vector ⃗ S 0 is not in the ring (xz) plane (see Fig. 1). ⃗ S 3 ðtÞ ¼ U ring ð⃗ c; t − n · T rev Þ |fflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl} rest of last turn

B. Spin-rotation angle in a static solenoid
In a solenoidal magnet with a field integral BDL ¼ R B k dl, the spins are rotated around the longitudinal direction ⃗ e z , and the rotation angle is given by The spin-rotation angle in the solenoid for deuterons at a momentum of P ¼ 970 MeV=c, normalized to the magnetic field integral, amounts to C. Spin tune and spin closed orbit with solenoids using ⃗ S 3 ðtÞ In the following, the abbreviation, e.g., χ Sol 1 rot ¼ χ 1 is used. For an ideal ring, free of magnetic imperfections, the spin tune change Δν s ðχ 1 ; χ 2 Þ, due to solenoids S 1 and S 2 in the ring (see Fig. 15), the left side of Eq. (30) of Ref. [18] can be approximated by πΔν s ðχ 1 ; χ 2 Þ · sinðπν 0 s Þ, where ν 0 s denotes the unperturbed spin tune in the machine. For small spin-rotation angles in the solenoids, Eq. (30) of [18] can thus be approximated by In order to validate the spin-evolution equation for ⃗ S 3 ðtÞ, given in Eq. (100), in Fig. 16 the spin-tune changes Δν s are compared to the approximation of Eq. (103) for four different cases.

D. Spin-closed orbit in a nonideal lattice
The static solenoids or magnetic imperfections in the ring affect the spin-closed orbit vector ⃗ n s ¼ ⃗ c in the machine. The situation is similar to the one depicted in Fig. 1, but there, only the tilt due to the EDM was taken into account. The presence of static solenoids in the ring can be numerically evaluated using Eq. (97) with ⃗ S 3 ðtÞ from Eq. (100).
Since the time t begins to count right behind the rf Wien filter (see Fig. 15), evaluation of Eq. (97) at t ¼ T rev {or integer multiples of T rev [see Eq. (40)]}, yields the orientation of the spin-closed orbit vector ⃗ c at the rf Wien filter Figure 17 shows how the stable spin axis ⃗ c ¼ ðc x ; c y ; c z Þ at the rf Wien filter is affected by the two solenoids S 1 and S 2 , and the presence of an EDM d. For comparisons, a number of special cases was numerically evaluated and is listed in Table V.

E. Strength of the EDM resonance
As depicted in Fig. 14, and already discussed in Sec. III E 4, the operation of the rf Wien filter modulates the spin tune. While the average spin tune is equal to the one obtained when the rf Wien filter is switched off, solenoids and magnet misalignments in the ring, however, affect the spin tune. Therefore, the spin-precession frequency and thus the frequency at which the rf Wien filter should be operated, differs from the unperturbed spin tune. The spin tune ν s must be determined anew for every  Fig. 15) under the conditions of Table I  solenoid setting to ensure that the resonance frequency for the rf Wien filter is given by and this frequency needs to be used in ψðtÞÞ [Eq. (35)], as it controls the rf Wien filter spin-rotation matrix R½⃗ n WF ; ψðtÞ [Eq. (38)]. The EDM resonance strength ε EDM , actually a resonance tune, is defined as the ratio of the angular velocity of the vertical polarization oscillation Ω p y induced by the EDM relative to the orbital angular velocity Ω rev , Since Ω p y corresponds to ω [first line in Eq. (79)], the resonance strength can in principle be determined from the sole observation of Ω p y . Alternatively, the resonance strength can be determined from the last line in Eq. (79) via   Fig. 15). but this requires that the initial slopes need to be determined as a function of, e.g., ϕ ¼ ϕ rf . The statistical aspects of this will be further elucidated in Sec. IV E 2.
1. Evolution of p y ðtÞ as a function of ϕ WF rot and χ Sol 1 rot The EDM resonance strength ε EDM [Eq. (106)] manifests itself in the oscillation frequency, as illustrated in Fig. 18 for two pairs of Wien filter rotation angle and spin-rotation angle in the solenoid S 1 , ðϕ WF rot ; χ Sol 1 rot Þ, while χ Sol 2 rot ¼ 0. The resulting oscillation pattern of p y is fitted using amplitude a and frequency ω are given in each panel, together with various other parameters. The calculation for the ideal ring situation in panel (b) uses a 1000 times larger assumed EDM value of d ¼ 10 −17 e cm and a larger number of turns n turns ¼ 100 000, in order to make the oscillations of p y ðtÞ visible when the rf Wien filter does not contribute a sideways magnetic field.
2. Comparison of ε EDM from Ω p y and _ p y ðtÞj t = 0 by variation of ϕ rf One would expect that the variation of the rf phase ϕ rf will affect the resulting oscillation amplitudes a and offsets b of Fig. 18, while the oscillation frequencies ω, and thus the resonance strengths ε EDM , remain unchanged.
In the panels of Fig. 19, for the same combinations of ðϕ WF rot ; χ Sol 1 rot Þ, shown in Fig. 18, _ p y ðtÞj t¼0 and the oscillation frequency ω are computed for 36 randomly picked values of ϕ rf . The graph illustrates that in the presence of solenoid fields and rf Wien filter rotations, the determination of _ p y ðtÞj t¼0 by variation of ϕ rf , making use of Eq. (107) yields results comparable to the direct determination of the resonance strength from the oscillation frequency Ω p y via Eq. (106). The oscillation amplitudes a and _ p y j t¼0 exhibit an identical dependence on ϕ rf , while the obtained resonance tune ε EDM remains constant over the whole range of ϕ rf .
The resonance strengths extracted from _ p y ðtÞj t¼0 and Ω p y make use of the very same simulated data. The results are summarized in Table VI, where for the numbers that should match, the markers (A), (B), and (C) are used. Although the different extraction methods show good overall agreement, the uncertainties of ε EDM ðΩ p y Þ, however, are substantially smaller than those from ε EDM ð _ p y j t¼0 Þ by a factor of at least 20. The reason for this is that in general frequencies can be measured more accurately than other quantities, and the determination of ε EDM ðΩ p y Þ involves fewer uncertainties in the error propagation. The most accurate determinations are obtained from Ω p y when χ Sol 1 rot ¼ 0. In the following, we briefly comment on some features of the results obtained so far (Fig. 18, Table VI). We observe that numerically 2 sin πν s ¼ 1.0041 ≃ 1. Then, according to the Appendix, we expect in good agreement with the results shown in Fig. 19. The resonance tunes determined from _ p y j t¼0 and from Ω p y are identical. For the above reason of 2 sin πν s ≃ 1 and small EDM contributions, the following equalities hold: F. Resonance strength ε EDM for random points ðϕ WF rot ;χ Sol 1 rot Þ The resonance strengths shown in Fig. 20 are obtained using the fit function of Eq. (108) (ω ¼ Ω p y ) and then Eq. (106) for a set of randomly chosen pairs of ðϕ WF rot ; χ Sol 1 rot Þ and χ Sol 2 rot ¼ 0. For all points, ϕ rf ¼ 0 and ⃗ S 0 ¼ ð0; 0; 1Þ,  Table I apply. The rf Wien filter is operated at harmonic K ¼ −1. The extracted resonance strengths are summarized in Table VI.  Fig. 20(b)]. The shift amounts to about 0.18°¼ sinðξ EDM ¼ 0.3 mradÞ. The relative uncertainties of the points shown in Fig. 20, obtained from the fits, range from Δε EDM =ε EDM ¼ 2.0 × 10 −5 to 4.1 × 10 −2 .
For the set of points ðϕ WF rot ; χ Sol 1 rot Þ shown in Fig. 20, the initial spin tunes ν s , i.e., before the rf Wien filter is turned on, are shown in Fig. 21. The result indicates the familiar quadratic dependence Δν s ðχ 1 ; χ 2 ¼ 0Þ ∝ χ 2 1 , as described by Eq. (103).

Operation of rf Wien filter exactly on resonance
In this section, the contour of the surface ε EDM ðϕ WF rot ; χ Sol 1 rot Þ, shown in Fig. 20(a), is compared to the theoretical expectation, given in Eq. (A5). The functional dependence describes a quadratic surface, also know as elliptic paraboloid, and is used here in the form of where the unperturbed spin tune ν ð2Þ s for the EDM of d ¼ 10 −18 e cm, assumed in the simulation, is given by It should be emphasized that the simulations shown in Fig. 20 reflect the situation when the rf Wien filter is operated exactly on resonance. During the corresponding EDM experiments in the ring, however, a certain spin-tune feedback is imperative to maintain the resonance condition for long periods of time, i.e., the spin-precession frequency in Eq. (35), using the measured spin tune [27]. To maintain phase and frequency lock when the rf Wien filter is actively operating, turns out to be much more tricky, and more sophisticated approaches, beyond those outlined in [23], are presently being pursued by the JEDI collaboration. Only such a phase and frequency lock during a measurement cycle enables one to take full advantage of the large spin-coherence time (SCT) of τ SCT ≃ 1000 s, achieved by JEDI at COSY [28,29].
The result of a fit without weighting is shown in Fig. 22(a). The fit parameters are listed in Table VII. It should be noted that within the uncertainties obtained from the fit, A ¼ B, while C and χ 0 are compatible with zero. Here, χ 0 represents a primordial tilt of the stable spin axis at the rf Wien filter along the horizontal axis, c x . For the model ring, one would expect a property which is nicely returned by the fit shown in Fig. 22(a). In addition, the fit to the simulated data is expected to return ϕ 0 ¼ jξ EDM ðd ¼ 10 −18 e cmÞj ¼ 0.3054 mrad, given by Eq. (11), and the fitted result returns this value accurately.

Validation of the scale of ε EDM
The fit with the elliptic paraboloid, shown in Fig. 22(a), indicates that the surface is described with A ¼ B. In the following, the first fit function from Eq. (111) is slightly altered, yielding where a factor F ¼ 10 20 has been introduced to scale the resonance strength. The second fit now uses weights derived from the uncertainty of the fitted Ω p y using Eq. (106). The fit obtained is shown in Fig. 22(b), and the results are summarized in Table VIII. The agreement between the theoretical model and the simulated data is good, the χ 2 =ndf ¼ 374.4=194 ¼ 1.9. According to Eq. (A5), the factor in front of the brackets in Eq. (115) reads where the Wien filter rotation angle ψ WF from Eq. (34) is used. Inserting the numerical value of D from the fit (Table VIII), and taking into account that the results are in mrad, the ratio   Fig. 22(b). The fit parametrizes the simulated data shown in Fig. 20(a) compared to the first one, shown in Fig. 22(a), and χ 0 and C are both compatible with zero.

V. CONCLUSIONS AND OUTLOOK
The SOð3Þ matrix formalism used here to describe the spin rotations on the closed orbit, i.e., the spin dynamics of the interplay of an rf Wien filter with a machine lattice that includes solenoids, proved very valuable. The general features of the deuteron EDM experiment at COSY can be obtained rather immediately.
The polarization evolution in the ring in the presence of an rf Wien filter that is operated on a parametric resonance, in terms of the resonance tune or resonance strength ε EDM is theoretically well understood. This will allow us to investigate in the future effects of increasingly smaller magnetic imperfections, either through additional solenoidal fields in the ring, or by transverse magnetic fields via the rotation of the rf Wien filter around the beam axis.
In the near future, it is planned to incorporate into the developed matrix formalism also dipole magnet displacement and rotation parameters, available from a recent survey at COSY. This will allow us to determine the orientation of the stable spin axis of the machine at the location of the rf Wien filter, and to extract the EDM from a measurement of the resonance strengths as a function of ðϕ WF rot ; χ Sol 1 rot Þ. In addition, it shall be possible to incorporate the spin rotations from misplaced and rotated quadrupole magnets on the closed orbit into the formalism as well. Of course, the approach taken is no substitute for more advanced spin-tracking codes, but the results obtained here can be applied to benchmark those codes.
It should be noted that the JEDI collaboration is presently applying beam-based alignment techniques to improve the knowledge about the absolute beam positions in COSY [15]. Once the orbit corrections based on the results of the beam-based alignment have been implemented, the approach described here to parametrize the spin rotations solely on the basis of the closed orbit, will become even more realistic.
The collaboration devoted a considerable effort to the experimental optimization of the spin-coherence time [28,29]. For deuterons at momenta near 970 MeV=c, spin-coherence times in excess of 1000 s are routinely achieved nowadays by careful adjustment of the sextupole magnets in the ring. The collaboration is presently preparing the corresponding investigations for protons.
An obvious limitation of the analytic treatment presented here, is the implicit assumption that the beam emittance is vanishing and correspondingly, the spin-coherence time is infinitely large. In the near future, we intend to apply more sophisticated spin-tracking algorithms to better understand the relation between a finite beam emittance and the corresponding spin-coherence time, both for deuteron and proton beams. It remains to be seen, how far one can develop analytic descriptions to actually model the spin-coherence time and other spin dynamics aspects of beams with finite emittance.
A full-fledged analytic treatment of decoherence of the polarization in a beam bunch due to synchrotron oscillations was reported in [30], and these findings will be applied to the analysis of the experimental data. Yet another potential source of decoherence is intrabeam scattering, which evidently randomizes the synchrotron motion and thereby the spin phases of the stored particles. Effects from intrabeam scattering can be readily incorporated in the formalism exposed in [30].
Recently, in conjunction with the design of the rf Wien filter [5], an approach based on the polynomial chaos expansion has been successfully applied to determine a hierarchy of uncertainties. Such a methodology, in combination with the spin-tracking approach based on the matrix formalism outlined here, can be employed to efficiently generate a hierarchy of uncertainties for the EDM prototype ring [11], based on the design parameters of the machine.
The employed spin-tracking approach shall be also applied to study various aspects of the presently applied spin-tune feedback system, which is used to phase lock the spin vector to the rf of the Wien filter [23].
At the location of the polarimeter, only the vertical and radial components of the beam polarization [S y ðtÞ and S x ðtÞ] can be determined. At the rf Wien filter, the orientation of the stable spin axis is denoted by ⃗ c, and in EDM mode the direction of the magnetic field by ⃗ n WF [see Eq. (39)]. The in-plane S x ðtÞ thus obviously depends on ½⃗ n WF × ⃗ c. In an ideal all-magnetic ring under consideration, the stable spin axis is close to the vertical direction ⃗ e y , ⃗ c ¼ cos ξ EDM · ⃗ e y þ sin ξ EDM · ⃗ e x ≈ ⃗ e y þ ξ EDM · ⃗ e x : ðA1Þ In EDM mode, the magnetic axis of the rf Wien filter can be approximated by The stable spin axis ⃗ c can be manipulated by static solenoids in the ring, and the drift solenoids S 1 and S 2 of the electron coolers (or the Siberian snake instead of S 1 ) generate the spin kicks χ 1 and χ 2 . When both solenoids S 1;2 are turned on, one can write for the stable spin axis c z ¼ 1 2 sin πν s ðχ 1 þ χ 2 cos πν s Þ: ðA3Þ In case solenoid S 2 is off (χ 2 ¼ 0), one obtains Thus the resonance strength squared can be written as a sum of two independent quadratic functions, where ψ WF is defined in Eq. (34).