On solving the Thomas Bargman-Michel-Telegdi equation using the Bogoliubov Krylov method of averages and the calculation of the Berry phases

TL;DR

This paper applies the Bogoliubov-Krylov-Mitropolski method to solve the Thomas Bargman-Michel-Telegdi equation and compute Berry phases, aiding precise spin evolution simulations for charged particle EDM measurements.

Contribution

It introduces a novel application of the Bogoliubov-Krylov-Mitropolski method to the T-BMT equation and Berry phase calculations in proton EDM storage rings.

Findings

Validated the numerical method against analytical estimates.

Demonstrated convergence of simulations with lattice imperfections.

Enhanced understanding of systematic errors in EDM experiments.

Abstract

Several proposals aimed at measuring the Electric Dipole Moment (EDM) for charged particles require very precise simulations and understanding of the systematic errors that can contribute to a spin buildup mimicking the EDM signal to be detected. In what follows, one used the Bogoliubov-Krylov-Mitropolski method of averages to solve the T-BMT equation and calculate the Berry phases arising for a proton EDM storage ring. The formalism employed proved to be particularly useful to determine the evolution of the spin at the observation point, i.e. at the location of the polarimeter. Several selected cases of lattice imperfections were simulated and benchmarked with the analytical estimates. This allowed the proof of the convergence of the numerical simulations and helped gain better understanding of the systematic errors.

Tables

Table 1. Table 1: Table of the ring parameters of the proton EDM experiment.

Total beam energy	1.171 GeV
Focusing structure	FODO
$N_{cells}$, number of cells	20
Deflector shape	cylindrical
Number of deflectors per cell	6
Bending radius $\rho$	52.3089 m
Radial E field	8.016 MV/m
Gap	3 cm
Bending voltage	$\pm$ 120 kV
Horizontal tune $Q_x$	2.42
Vertical tune $Q_y$	0.44

Table 2. Table 2: Leading terms of the geometric phases and contributing fields.

Geometric phases (leading terms)	Acting fields
	${\Omega }_l\propto (y^{\prime }{E}_r,x^{\prime }{E}_y,y^{\prime }{B}_y,x^{\prime }{B}_x,B_l)$	${\Omega }_y\propto (E_r\Delta p/p,x^{\prime }{E}_l,y^{\prime }{B}_l,B_y)$
	Leading terms $\propto y^{\prime }{E}_r^2\Delta p/p,y^{\prime }{E}_r{B}_y,B_l{B}_y$
	${\Omega }_l\propto (y^{\prime }{E}_r,x^{\prime }{E}_y,y^{\prime }{B}_y,x^{\prime }{B}_x,B_l)$	${\Omega }_r\propto (E_y,y^{\prime }{E}_l,x^{\prime }{B}_l,B_x)$
	Leading terms $\propto {\left(y^{\prime }{E}_r\right)}^2{B}_x,B_l^2{B}_x$
	${\Omega }_r\propto (E_y,y^{\prime }{E}_l,x^{\prime }{B}_l,B_x)$	${\Omega }_y\propto (E_r\Delta p/p,x^{\prime }{E}_l,y^{\prime }{B}_l,B_y)$
	Leading terms $\propto {\left(E_r\Delta p/p\right)}^2{B}_x,B_y^2{B}_x$

Equations

$$\frac{d\text{s}}{dt}=\left({\mathbf{\Omega }}_{\mathbf{BMT}}+{\mathbf{\Omega }}_{\mathbf{EDM}}\right)\times \mathbf{s}$$

$${\mathbf{\Omega }}_{\mathbf{BMT}}$$

$${\mathbf{\Omega }}_{\mathbf{EDM}}=-\frac{e}{mc}\frac{\eta }{2}\left[\mathbf{E}-\frac{\gamma }{\gamma +1}(\mathbf{\beta }.\mathbf{E})\mathbf{\beta }+c\mathbf{\beta }\times \mathbf{B}\right]$$

$$\frac{d\text{s}}{ds}=\left(\frac{dt}{ds}{\mathbf{\Omega }}_{\mathbf{BMT}}-\mathbf{\kappa }\times {\mathbf{e}}_{\mathbf{l}}\right)\times \mathbf{s}$$

$$\frac{dt}{ds}=\frac{1+x/\rho }{{\beta }_{l}c};e{E}_r^b=-\frac{{\gamma }_m{E}_0{{\beta }_m}^2}{\rho };\mathbf{\kappa }=\frac{1}{\rho }{\mathbf{u}}_{\mathbf{x}}$$

$$\left(\begin{array}{c}d{s}_r/ds\\ d{s}_y/ds\\ d{s}_l/ds\end{array}\right)$$

$$\left(\begin{array}{c}d{s}_r/dt\\ d{s}_y/dt\\ d{s}_l/dt\end{array}\right)$$

$$=\left(\begin{array}{ccc}0& -{\Omega }_l& {\Omega }_y\\ {\Omega }_l& 0& -{\Omega }_r\\ -{\Omega }_y& {\Omega }_r& 0\end{array}\right)\left(\begin{array}{c}{s}_r\\ s_y\\ s_l\end{array}\right)$$

$$\mathbf{\beta }=\frac{\mathbf{p}c}{E}$$

$$\beta =\frac{1+\frac{\Delta p}{p}}{{\left[{\left(1+\frac{\Delta p}{p}\right)}^2+{\left(\frac{E_0}{p_{ref}c}\right)}^2\right]}^{1/2}}$$

$$\gamma =\frac{p_{ref}c\left(1+\frac{\Delta p}{p}\right)}{\beta E_0}={\left[1+{\left(\frac{p_{ref}c}{E_0}\right)}^2{\left(1+\frac{\Delta p}{p}\right)}^2\right]}^{1/2}$$

$$p_{x,y}=\frac{P_{x,y}}{P_{ref}}=\frac{{\beta }_{r,y}\gamma }{{\beta }_{ref}{\gamma }_{ref}}=\frac{{\beta }_{r,y}}{\beta }\left(1+\frac{\Delta p}{p}\right)$$

$${\beta }_r$$

$${\beta }_y$$

$${\beta }_l$$

$${\Omega }_l$$

$${\Omega }_r$$

$${\Omega }_y$$

$$a+\frac{1}{\gamma +1}-\frac{1}{{\gamma }_m{\beta }_m^2}$$

$$p_m=\frac{mc}{\sqrt{a}}=\frac{938.2721\text{mbox}\text{ MeV/c}}{\sqrt{1.79285}}=700.7398\text{ MeV/c}$$

$$a+\frac{1}{\gamma +1}-\frac{1}{{\gamma }_m{\beta }_m^2}$$

$${\gamma }_m-\gamma$$

$${\Omega }_y$$

$$\begin{cases}\Omega_{y}=K-\dfrac{e}{mc}\left[a+\dfrac{1}{\gamma+1}\right]\dfrac{\beta_{l}}{1+x/\rho}x^{\prime}E_{l}\\ \\ \Omega_{l}=\dfrac{e}{mc}\left[a+\dfrac{1}{\gamma+1}\right]\dfrac{\beta_{l}}{1+x/\rho}(x^{\prime}E_{y}-y^{\prime}E_{r})\\ \\ \Omega_{r}=\dfrac{e}{mc}\left[a+\dfrac{1}{\gamma+1}\right]\beta_{l}\left(\dfrac{y^{\prime}}{1+x/\rho}E_{l}-E_{y}\right)\\ \\ K\approx\left\{\begin{array}[]{@{}ll@{}}\dfrac{-e/mc}{a\left(\gamma_{m}+1\right)^{2}\gamma_{m}}\beta_{l}E_{r}^{b}\dfrac{\Delta p}{p}+\dfrac{e}{mc}\left[a+\dfrac{1}{\gamma+1}\right]\beta_{l}(E_{r}-E_{r}^{b})-\beta_{l}c\dfrac{x}{\rho^{2}}&\text{if bend}\\ \\ \dfrac{e}{mc}\left[a+\dfrac{1}{\gamma+1}\right]\beta_{l}E_{r}&\text{otherwise}\end{array}\right.\end{cases}$$

$$\frac{d\mathbf{S}}{dt}=\mathbf{M}(t)\mathbf{S}(t)=ϵ{\mathbf{M}}_1(t)\mathbf{S}(t)=ϵ\mathbf{Y}(t,\mathbf{S})$$

$$⟨\mathbf{Y}(t,\mathbf{S})⟩=\underset{T\to \infty }{\lim }\frac{1}{T}{\int }_0^T\mathbf{Y}(t,\mathbf{S})dt$$

$$\mathbf{S}(t)=\mathbf{\xi }(t)+ϵ\tilde{\mathbf{Y}}(t,\mathbf{\xi })$$

$$\frac{d\mathbf{\xi }}{dt}=ϵ⟨\mathbf{Y}(t,\mathbf{\xi })⟩=⟨\mathbf{M}(t)⟩\mathbf{\xi }$$

$$\tilde{\mathbf{Y}}(t,\mathbf{\xi })=\int \left[\mathbf{Y}(t,\mathbf{\xi })-⟨\mathbf{Y}(t,\mathbf{\xi })⟩\right]dt$$

$$\mathbf{S}(t)=\left[1+\tilde{\mathbf{M}}(t)\right]\mathbf{\xi }(t)$$

$$\frac{d\mathbf{\xi }}{dt}=⟨\mathbf{M}(t)⟩\mathbf{\xi }$$

$${\lambda }_1$$

Full text

††preprint: APS/123-QED

I Introduction

A storage ring represents a new attractive method to search and measure the electric dipole moment of the proton (or other charged particles) by using polarized particles at their magic momentum farley . However, to reach the sensitivity level required for such a measurement, it is crucial to understand and mitigate the systematic errors that can yield a fake signal mimicking the EDM one. The objective of this paper is to contribute some better understanding regarding that matter: starting from the spin precession equation, we will establish the formalism and all necessary quantities to compute the spin evolution in a storage ring. Then, using a perturbation method, an approximate solution to this equation will be derived and benchmarked with numerical simulations in MATHEMATICA. Then, using the tracking code BMAD, we will investigate the impact of several imperfections on the spin buildup for a fully electrostatic ring lattice.

II Frozen Spin

The Thomas Bargman-Michel-Telegdi (T-BMT) equation gives the precession rate of the angle between the spin and momentum vectors of a relativistic particle in the presence of electromagnetic fields thomas bargman :

[TABLE]

where

[TABLE]

is the precession vector due to the particle’s magnetic moment and

[TABLE]

is the precession vector due the particle’s finite electric dipole moment edm . The charged particle considered throughout this paper is the proton.

The above equation describes the spin changes with the time $t$ in the laboratory frame. s is the classical spin vector expressed in the rest frame of the particle while $\bm{B}$, $\bm{E}$ are the magnetic and electric field vectors expressed in the laboratory frame. $\bm{B}_{\perp}$ and $\bm{B}_{\|}$ are the components perpendicular and parallel to the particle’s momentum. $a$ is the particle’s anomalous gyro-magnetic factor and $\gamma$ is the relativistic gamma factor.

In accelerator coordinates, i.e. in a frame co-moving with the reference particle (curvilinear coordinates system), the above equation can be re-written in the following way hoff :

[TABLE]

where $s$ is the arc length of the reference particle, $\bm{e_{l}}$ or $\bm{e_{z}}$ is the unit vector in the longitudinal direction and $\bm{\kappa}$ is the curvature vector as illustrated in fig 1.

For simplicity, we will first assume that there are no magnetic fields and no EDM effects, i.e. $\eta=0$ and $\bm{B}=\bm{0}$.

Defining $\rho$ as the radius of curvature of the reference orbit, one can write:

[TABLE]

where $\beta_{m}$ and $\gamma_{m}$ are the Lorentz factors at the magic energy, $E_{0}$ the energy of mass of the proton and $E_{r}^{b}$ is the radial electric field of the electrostatic bends for the ideal particle. Thus, in the following formalism, the spin buildup is measured with respect to the planar reference orbit silenko . This will be the reference frame that we will use both for the analytical estimates as well as the numerical simulations.

In matrix notation, Eq. (4) writes:

[TABLE]

where the subscripts $r$ or $x$ denote equivalently the horizontal transverse coordinate pointing towards outer radii, $y$ denotes the vertical transverse coordinate orthogonal to the median plane of the unperturbed particle motion and $l$ or $z$ the longitudinal coordinate in the Frenet Serret frame. Differentiating with respect to time, one obtains:

[TABLE]

which writes in the general form $d\bm{S}/dt=\bm{M}\bm{S}$. Now, we compute the spin precession components ($\Omega_{r}$, $\Omega_{y}$, $\Omega_{l}$) in order to express them as a function of the accelerator coordinates $(x,x^{\prime},y,y^{\prime},z,\Delta p/p)$ in the Frenet Serret system (where ${}^{\prime}=d/ds$). For that, one has to express $(\beta_{r},\beta_{y},\beta_{l})$ as a function of $(x,x^{\prime},y,y^{\prime},z,\Delta p/p)$.

II.1 Calculation of the Lorentz factors

Given that $\bm{p}=\gamma mc\bm{\beta}$ and $E=\gamma mc^{2}=\gamma E_{0}$ where $E_{0}$ is the mass energy, then the normalized particle velocity is obtained:

[TABLE]

Recalling that $E^{2}=p^{2}c^{2}+E_{0}^{2}$, it results that for a particle with a momentum offset $\Delta p$ and slope $\Delta p/p=(p-p_{ref})/p_{ref}$:

[TABLE]

In electrostatic elements, the change in energy implies a change in the particle momentum, i.e. a variation of $\Delta p/p$. It results from Eq. (9) that:

[TABLE]

In accordance with BMAD bmad , which has been used for the simulation presented in this paper, the phase space momenta are defined as:

[TABLE]

and

[TABLE]

II.2 Calculation of the spin precession vectors

Injecting Eqs. (13) and (14) into the expression of $\Omega_{l}$ and $\Omega_{r}$, one obtains:

[TABLE]

and

[TABLE]

By means of Eq. (5), the horizontal spin precession vector can be calculated for an electrostatic bend:

[TABLE]

where $E_{r}^{b}$ is the design radial electric field of the electrostatic deflector.

If $\gamma=\gamma_{m}$, i.e. for a particle at the reference energy and with no transverse displacements ($x=x^{\prime}=y=y^{\prime}=0$),

[TABLE]

The above quantity vanishes at the momentum $p=mc/\sqrt{a}$ which is generally referred to as the magic momentum since then, this guarantees that $\Omega_{r}=\Omega_{y}=\Omega_{l}=0$, i.e. the spin precession vanishes. This is referred to as the frozen spin condition and is considered as a robust way to improve the sensitivity of the EDM measurement for charged particles. The main idea is to keep the spin aligned along the momentum vector of the particle and to provide a strong enough radial electric field acting on the EDM for the duration of the beam storage time farley .

For protons, this is achieved for,

[TABLE]

corresponding to a proton kinetic energy of $E_{m}\sim 232.79$ MeV.

In reality though, the particle momentum deviates from the magic value. Thus, it is important to evaluate the impact of any momentum offset on the spin precession rate.

Back to Eq. (18). Let’s assume that $\gamma\neq\gamma_{m}$. Then, one can write:

[TABLE]

Substituting Eq. (11) into the first order term of the above expression yields:

[TABLE]

Injecting the latter into the expression of $\Omega_{y}$, one finally obtains:

[TABLE]

In the above derivation, it is important to note that the spin precession vector is calculated with respect to the momentum of the reference particle orbit. In summary, the expressions of the spin precession components are as follows:

[TABLE]

III On solving the Thomas BMT equation using perturbation method

Eq. (8) can also be written in the matrix form:

[TABLE]

where $\epsilon>0$ is a small parameter. Thus, the Bogoliubov-Krylov-Mitropolski (BKM) method of averages can be invoked BKM . The main idea of the averaging method is based on the assumption that the derivative $d\bm{S}/dt$ is small, which is typically the case for the frozen spin method considered here (see (Eq. 20)). Thus, the evolution of $\bm{S}$ contains two contributions (obeying two timescales): a slowly varying term $\bm{\xi}$ due to the presence of $\epsilon$ in front of $\bm{M}_{1}$, and small rapidly oscillating terms, due to the presence of $t$ in $\bm{M}_{1}$, i.e. describing the spin precession changes within the elements.

The basic idea of this approach was first developed by Krylov and Bogoliubov (1934). Later on, in 1958, Bogoliubov and Mitropolski established the general scheme and a more rigorous treatment for this method. Finally, in 1969, Perko almost completed the theory with error estimates for the periodic and quasi-periodic cases perko .

It is assumed that the vector function $\bm{Y}(t,\bm{S})$ possesses an average value:

[TABLE]

This condition is satisfied for all cases considered here given that $\bm{Y}$ is almost periodic in time. The asymptotic solution of Eq. (25) can be obtained in the following way: in a first approximation in terms of $\epsilon$, one can write:

[TABLE]

where $\bm{\xi}(t)$ satisfies the equation:

[TABLE]

The average is to be taken only with respect to the explicit variable $t$ and the tilde $\bm{\sim}$ is the integrating operator defined by:

[TABLE]

From Eq. (25), (27) and (29), it results that a first order approximate solution to the T-BMT equation where the frozen spin condition is quasi-satisfied ($\epsilon$ small) is given by:

[TABLE]

To simplify the analysis, let’s solve for the slowly varying term in $\bm{S}$. From what preceded (Eq. 28), the zero order approximation of $\bm{S}$ is a solution of:

[TABLE]

In order to apply the above formalism, one has to compute the average of the matrix $\bm{M}(t)$ by tracking the particle over a sufficient number of turns, or in a simpler way by searching for the closed orbit which satisfies the periodicity of $\bm{M}(t)$ and thus the averages are defined over one full turn.

NB: In the above definitions, $\left<\right>$, $\bm{\sim}$ represent operators acting on all the elements of the matrix.

$\left<M\right>$ is a skew symmetric matrix. So all eigenvalues are imaginary. Solving for the eigenvalues $\lambda_{i}$, one obtains:

[TABLE]

and the corresponding eigenvectors are:

[TABLE]

A general solution to this equation writes in the form: $\bm{\xi}(t)=c_{1}\bm{v_{\lambda_{1}}}e^{\lambda_{1}t}+c_{2}\bm{v_{\lambda_{2}}}e^{\lambda_{2}t}+c_{3}\bm{v_{\lambda_{3}}}e^{\lambda_{3}t}$. Now, solving for $\xi$ such that $\xi$(0)=(0,0,1), one obtains the zero order solution:

[TABLE]

Assuming that the argument of the $\cos$ and $\sin$ functions is small (EDM signal level), then Eq. (35) may be rewritten:

[TABLE]

To verify the validity of the above solutions, one can compute the norm of $\bm{\xi}(t)$ and it can be shown that for any given $t$, $\parallel\bm{\xi}(t)\parallel=({\xi_{r}}^{2}+{\xi_{y}}^{2}+{\xi_{l}}^{2})^{1/2}=\parallel\bm{\xi}(0)\parallel=1$. The conservation of the Hermiticity of the zero order solutions is thus verified. From Eqs. (35) and (30) it results that, to the first order in the method of averages, the solution of the T-BMT equation is given by:

[TABLE]

where

[TABLE]

In the above derivation, it is assumed that $\bm{S}$ changes slowly. This is formulated in such a way that the right hand side term of Eq. (25) is proportional to $\epsilon$, therefore $\bm{S}\sim\epsilon$. Thus, $\bm{S}$ changes appreciably only in a large time interval $t\sim 1/\epsilon$. Of course, the above argument is not rigorous. However, the foundation of the method of averaging was established by N. N. Bogoliubov who demonstrated that one can construct a system of first approximation and higher orders, the solutions of which approximate the solutions of the initial system within an arbitrary degree of accuracy. It is important to note that the tilde parameters in Eq. (39) give rise only to small vibrations around $\xi$ and do not influence the systematic change of the spin.

Now, let’s characterize the evolution of the tilde parameters on the closed orbit:

Property

*“On the closed orbit, the tilde parameters vanish after the completion of each turn.”

*To establish this, let’s consider a quantity $x$. $x$ refers to the field or the coordinates of the particle. Such a quantity defined on the closed orbit is periodic in $T$. Thus, without any loss of generality, $x$ can be written in the following way:

[TABLE]

where $g_{x}(t)$ is a periodic function in $T$ such that $\left<g_{x}\right>=0$. Now,

[TABLE]

Thus, given that $\left<g_{x}\right>=0$, one obtains

[TABLE]

III.1 Spin rotation matrix

We start from

[TABLE]

Then, assuming that $\xi(t)=R(t)\xi(t_{0})$ where $R(t)$ is the spin rotation matrix satisfying $R(t_{0})=I$, one obtains:

[TABLE]

The solution is obtained using the Euler-Rodrigues formula:

[TABLE]

Now, assuming that $\xi(0)=(0,0,1)$, one obtains:

[TABLE]

which is the same result as obtained in Eq. (35) above.

III.2 2nd order approximation

In reality, when solving for the vertical spin component, the buildup is determined by the product of the rapidly oscillating terms of the particle motion (velocities and experienced fields) with the rapidly oscillating spin components. In certain configurations where the average value of the rapidly oscillating terms of the spin precession vanishes i.e. $\left<\Omega\right>=0$, a spin buildup still exists. This is due to the non-commutativity of spin rotations around different axes and is generally referred to as geometric phase. Such a phenomena, also known as Pancharatnam-Berry phase geom1 or often Berry phase geom2 is defined as the phase difference acquired over the course of a cycle, when a system is subjected to cyclic adiabatic processes. Therefore, to account for such a contribution to the vertical spin buildup, a second order approximation is employed in which one can integrate the Thomas BMT equation by re-injecting the first order solution based on the method of averages. First, one re-writes the radial spin buildup from the exact T-BMT equation:

[TABLE]

where $s_{l}$ and $s_{y}$ are solved using the method of averages (Eq. 39):

[TABLE]

and $\phi_{l}$, $\phi_{y}$ are the rapidly oscillating terms that vanish after each turn completion (on the closed orbit) as established in the property above and defined by:

[TABLE]

Next, one injects this improved first approximation into the exact equation (48). This yields:

[TABLE]

Assuming that $\xi_{y}(t),\xi_{r}(t)\ll\xi_{l}(t)$, and keeping the terms up to the second order, one can integrate the above equation (assuming $\xi_{l}\approx 1$ and $\xi_{y}\approx-\left<\Omega_{r}\right>t$) to compute the radial spin component:

[TABLE]

Now, one expresses $\Omega_{l}$ in the following way:

[TABLE]

Thus, by means of an integration per parts, the expression of the integral can be simplified:

[TABLE]

so that Eq. (54) re-writes in the following way:

[TABLE]

where the subscript 2 denotes the second order approximation for which the frozen solution is given by:

[TABLE]

and the rapidly oscillating terms that vanish after each turn completion are:

[TABLE]

Thus, the second order approximation yielded additional terms for the frozen solution in comparison with the $1^{st}$ order approximation:

• •

The first additional term $\left<\Omega_{l}\tilde{\Omega}_{r}\right>t$ accounts for the geometric phases induced by the rapid oscillations of the vertical spin component (although its average is neglected) and longitudinal spin precession. As will be established later on, such a contribution becomes dominant at the proximity to the magic energy ($\left<\Omega_{y}\right>\approx 0$).

• •

The second additional term $-\left<\Omega_{r}\right>\left<\tilde{\Omega}_{l}\right>t$ is due to the slowly linear varying term of the vertical spin component coupled with the rapidly oscillating terms of the longitudinal spin precession.

Next, one repeats the same steps to compute the vertical spin buildup from the exact T-BMT equation:

[TABLE]

where $s_{r}$ and $s_{l}$ are solved using the method of averages (Eq. 39):

[TABLE]

and $\phi_{r}$, $\phi_{l}$ are the rapidly oscillating terms that vanish after each turn completion (on the closed orbit) as established in the property above and defined by:

[TABLE]

Next, one injects this improved first approximation into the exact equation (60). This yields:

[TABLE]

where the third (Right Hand Side) term can be neglected. Now, injecting the approximate expressions of $\xi$ ($\xi_{l}\approx 1$, $\xi_{r}\approx\left<\Omega_{y}\right>t$), keeping the terms up to the second order and integrating over time, one obtains:

[TABLE]

and injecting Eq. (56) into the above expression, one can write a second order approximation to the vertical spin component:

[TABLE]

where the subscript 2 denotes the second order approximation for which the frozen solution is given by:

[TABLE]

and the rapidly oscillating terms that vanish after each turn completion are:

[TABLE]

In summary, the second order frozen solution yields 4 terms contributing to the vertical spin buildup:

The first contribution is due to the longitudinal spin component coupled with the average contribution of the radial spin precession. This accounts for any net average radial magnetic field for instance. 2. 2.

The second contribution is due to the rapidly oscillating terms of the radial spin components coupled with the longitudinal spin precession. This accounts for the geometric phases as well for any net average longitudinal spin precession. 3. 3.

The third contribution is due to the slowly linearly varying term of the radial spin component coupled with the rapidly oscillating terms of the longitudinal spin precession. 4. 4.

The fourth and last contribution is due to the slowly linearly varying term of the radial spin component coupled with the average contribution of the longitudinal spin precession (quadratic increase).

Finally, one repeats the same steps to compute the second order approximation of the longitudinal spin buildup:

[TABLE]

Assuming $\xi_{l}\approx 1$, $\xi_{r}\approx\left<\Omega_{y}\right>t$, $\xi_{y}\approx-\left<\Omega_{r}\right>t$, one can integrate the above equation. First, one can establish the following identity:

[TABLE]

so that one can solve for the longitudinal spin component which writes in the following way:

[TABLE]

Remarks:

• •

The second order approximation above is modified from the original BKM method of averages: if one follows the original method, the second order approximation is obtained from the exact T-BMT equation by re-injecting the first order approximation and taking the averages with respect to the time as an explicit variable while considering $\bm{\xi}$ constant. Such an assumption is not rigorous since it neglects the impact of the product of the linearly varying terms of $\bm{\xi}$ with the rapid oscillations of the spin precession components $\bm{\Omega}$.

• •

If the average value of $\Omega_{l}$ or $\Omega_{y}$ is zero, then the quadratic evolution with time of the vertical spin vanishes. Therefore, clearing the frozen radial spin component by means of a feedback system can eliminate the quadratic buildup. However, this is only valid up to the second order. Higher order terms yield additional quadratic time dependence.

• •

The frozen solution $\xi$ (see Eq. (68)) is only valid at time $t$ which is an integer number of the revolution period. In other terms if the spin components are recorded after each turn completion. Therefore it is a natural quantity to describe the signal measured at the location of the polarimeter.

• •

If $\left<\Omega_{r}\right>=\left<\Omega_{y}\right>=\left<\Omega_{l}\right>=0$ then $\xi_{y,2}(t)\approx\left<\Omega_{l}\tilde{\Omega}_{y}\right>t$ which accounts for the geometric phases to the second order.

• •

The above calculation can be easily extended to the next order. Although there seems to be no need in general to pursue this, one might be interested though to determine the high order terms of the Berry phases. The final result of the third order frozen spin solution is given by:

[TABLE]

It results that the only contribution of the geometric phases up to the third order is given by:

[TABLE]

since $\left<\Omega_{r}(\tilde{\Omega}_{r})^{2}\right>=0$ if $\left<\Omega_{r}\right>=0$.

This will be discussed later on in this paper.

III.3 Error analysis

The above second (and third) order approximation to the T-BMT equation is based on the assumption that the average spin precession component is negligible on the timescales of the EDM experiment. If the spin coherence time is $T_{coh}=1000\text{ s}$ as is generally assumed to reach the sensitivity required to measure the EDM, then the previous assumption can be formulated as follows:

[TABLE]

From the above scheme one can infer that the general frozen solution to the T-BMT equation can be classified into three main regimes depending on the value of the average spin precession:

If $0<\left<\Omega\right>\ll 1/T_{coh}$ then any non-linear increase with time of the spin can be neglected on the timescales of the EDM experiment. And using the $N^{th}$ order approximation based on the Bogoliubov method of averages, it can be seen that:

[TABLE]

For instance, the second order approximation satisfies the T-BMT equation to an accuracy of the order of $\left<\Omega\right>^{3}t$.

In general though, the non-linearities arising from the solution can be neglected for short timescales in comparison with $T_{coh}$ ($\sim$ milliseconds). 2. 2.

If $\left<\Omega\right>\gtrsim 1/T_{coh}$ then the spin evolution is mainly governed by the averages of the spin precession components if the same property holds for all components i.e. $\left<\Omega_{y}\right>,\left<\Omega_{r}\right>,\left<\Omega_{l}\right>\gtrsim 1/T_{coh}$. Therefore, the first order solution (Eq. (35)) is the most accurate one. 3. 3.

In the limit where $\left<\Omega\right>=0$, i.e. $\left<\Omega_{r}\right>=\left<\Omega_{l}\right>=\left<\Omega_{y}\right>=0$, the geometric phases are the only contribution to the spin buildup. The latter is governed by the non-commutativity of the rotation around different axes. There appears no obvious way to estimate the error of the solution in this particular case, since the spin buildup is mainly determined by how the errors are distributed in the ring with respect to each other. However the third order approximation (see Eq. (73)) contains the dominant contributions. Note that the geometric phases which arise initially as a linear function of time can be transferred from one plane to the other therefore yielding additional terms.

It is important to note that the Hermiticity of the approximate frozen solution is only conserved for the $1^{st}$ order terms of the approximation. To show this, let’s consider the second order approximation in the particular case where $\left<\Omega\right>=0$ and calculate the norm of the frozen solution only at times $t=kT$ i.e. after each turn completion:

[TABLE]

Such an effect is negligible for the timescales of the EDM experiment. However, the Hermiticity can be improved by keeping the higher order terms in the expansion of the sinusoidal functions of the $1^{st}$ order solution. This will not be pursued here since the EDM experiment requires slowly varying parameters of the average spin precession.

III.4 Benchmarking of the analytical solution with numerical simulations

To establish the validity of the analytical solutions, one simulated several cases by solving the T-BMT differential equations using explicit Runge Kutta tracker in MATHEMATICA. The spin precession components are defined as a linear combination of sinusoidal functions. For all cases considered, one assumed $\left<\Omega_{r}\right>=0$ to simplify the analysis. For more details regarding the numerical solution, see Appendix A. Note that, for all cases considered here, unless otherwise specified, one relies on the second order approximation since the first order approximation is particularly inaccurate when the averages of the spin precession components vanish. This is generally the case at the proximity to the magic energy.

III.4.1 Geometric phases

In the first case considered, one assumes $\left<\Omega\right>=0$. The result is that the vertical spin buildup is determined by the geometric phases such as $\xi_{r,2}(t)\approx\left<\Omega_{l}\tilde{\Omega}_{r}\right>\approx 0$ and $\xi_{y,2}(t)\approx\left<\Omega_{l}\tilde{\Omega}_{y}\right>t=-1.1*10^{-6}t$ rad. Such an effect is always present and is strongly related to how the errors are distributed in the ring. A comparison of the numerical simulation with the analytical result is shown in fig 2.

We define $\epsilon_{abs}$ (resp. $\epsilon_{rel}$) as the absolute error (resp. relative error) between the convergent numerical solution and the analytical one. Fig 3 shows the absolute error which grows linearly with time (the relative error is therefore constant).

III.4.2 Geometric phases and growing radial spin component

In this case, $\left<\Omega\right>=\left<\Omega_{y}\right>=10$ rad/s which typically occurs when the particle is not at the magic energy. This yields a non-vanishing frozen radial spin component $\xi_{r,2}(t)\approx\left<\Omega_{y}\right>t=-10*t$ rad as well as an additional term contributing to the vertical spin buildup such as $\xi_{y,2}(t)\approx\left<\Omega_{l}\tilde{\Omega}_{y}\right>t-\left<\Omega_{y}\right>\left<\tilde{\Omega}_{l}\right>t=-1.1*10^{-6}t-2.2*10^{-5}t$ rad.

A comparison of the analytical results with the numerical simulations is shown in fig 4 where one can observe a good agreement. However, as illustrated in fig 5, the error grows cubically with time. This is explained by the fact that in the above derivations, one only considered the first order linear term in the expansion of $\xi_{r}$ (Eq. 35). This is valid for very small timescales. However, as the time grows, the approximation may no longer be valid, and this is particularly enhanced here by the fact that $\left<\Omega\right>=10$ rad/s is non negligible. Thus, one needs to include the next higher order term in the expansion, which yields the cubical time dependence since $\sin(t)\approx t-t^{3}/6$.

Note as well that the vertical spin oscillates with an amplitude growing over time. As predicted by Eq. (69), this is due to $\left<\Omega_{y}\right>\neq 0$. Thus, it is advantageous to use more than one polarimeter located at specific azimuthal locations around the ring to probe such imperfections.

III.4.3 Geometric phases, growing radial spin component and quadratic increase

In this case, one added a non-vanishing longitudinal spin precession such as $\left<\Omega_{l}\right>=10^{-2}$ rad/s in comparison with the previous case. This yields an additional quadratic spin buildup so that:

[TABLE]

Figure 6 displays the tracking results as well as a comparison with the analytical formula for the first 20 revolution periods where one can observe the impact of the quadratic spin buildup in comparison with the linear one. The error (see fig 7) is similar to the previous case although one can observe the existence of additional higher order terms.

III.4.4 No vertical spin precession and zero averages

A question that arises frequently is whether a beam injected at the magic energy so that the vertical spin precession vanishes everywhere can still yield a vertical spin buildup? To answer this, one assumes no vertical spin precession as well as zero averages of the spin precession components in the ring, i.e. $\Omega_{y}(t)=0$ and $\left<\Omega_{l}\right>=\left<\Omega_{r}\right>=0$. It results from the second order approximation that:

[TABLE]

so that the radial spin increase results from the net effect of the spin rotations around the longitudinal and radial axes. However, given that $\Omega_{y}(t)=0$, the second order approximation yielded no vertical spin buildup which is not in agreement with the numerical simulations. This is a sign that higher order terms are missing. To derive the 3rd order approximation of the vertical spin component in this particular case, one re-injects Eq. (79) into (60). This yields:

[TABLE]

or, after integration:

[TABLE]

so that the frozen spin solution is approximated by:

[TABLE]

Thus, in the above system, the vertical spin buildup arises first by means of radial spin precession. This leads to a linear radial spin buildup via rotations around the longitudinal axis that do not commute with the radial one. Such a radial spin buildup acts back on the vertical one by means of rotations around the longitudinal axis. These rotations transfer the geometric phases from the horizontal plane to the vertical one. The following diagram summarizes this scheme:

[TABLE]

and comparison of the tracking simulations with the analytical results are shown in fig 8.

Therefore, in a fully electrostatic ring, and for a beam injected at the magic energy, vertical spin buildup can occur by means of consecutive rotations around the longitudinal and radial axes induced by the electric fields and the non-vanishing slope of the particle trajectories.

Furthermore, assuming $\Omega_{l}(t)=0$ and $\left<\Omega_{r}\right>=\left<\Omega_{y}\right>=0$, it can be verified (although it will not be discussed here) that the vertical spin buildup can occur due to the rapidly oscillating terms of the longitudinal spin component coupled with the radial spin precession such as $s_{y}(t)\approx 1/2\left<\Omega_{r}(\tilde{\Omega}_{y})^{2}\right>t$.

III.5 System oscillating with several frequencies

In reality, the spin precession components can be oscillating with different frequencies, not necessarily equal to the particle revolution frequency. In this case, the question arises of whether there exists an efficient way to compute the averages of the previous equations without having to track a particle for a consequent number of turns? The question can also be formulated differently: can we define a period $T$ associated to the system under which the spin oscillations are represented as a sum of slowly varying terms and rapidly rotating terms that vanish at times $t=kT,k$ integer? To simplify the analysis, let’s assume that the system is oscillating with two frequencies: the revolution frequency of the particle since the stable motion is always near the closed orbit, and a second frequency perturbing the spin motion. The latter can be due to any external or internal force of the beam, mismatching issue, etc. For instance, the vertical spin precession component can be written in the following way:

[TABLE]

where $f_{1}$ and $f_{2}$ are two periodic functions in $t$ with periods $T_{rev}$ and $T_{per}$ respectively. Now, assuming there exists a rational number $p/q$ ($p,q$ mutually prime numbers) such as:

[TABLE]

this implies that $T$ is the frequency of the system as defined above. In reality, $T_{rev}/T_{per}$ is a real number. Therefore, the idea is to obtain its closest rational approximation. This can be done with the help of continued fractions.

III.6 Accuracy of the numerical simulations

In general, the higher the order of the integration method, the more accurate is the numerical scheme. Accuracy here refers to the smallness of the local error, i.e. the error introduced after each integration step of the differential equation. Therefore, to improve the accuracy of the solution, it is generally assumed that smaller step size is a sufficient condition. However, when stiff systems are considered, some of the solutions obtained can be highly unstable. Besides, it appears that the accuracy of the obtained solutions is highly sensitive to the step size chosen. Stiffness is mainly encountered when hard edge models are assumed, i.e. due to an abrupt change in the particle field or its coordinates. In addition, the total energy conservation is a fundamental criterion to assess the accuracy of the numerical simulations as is discussed in selcuk . Using the explicit Runge Kutta tracker available in MATHEMATICA, one investigated the impact of the stiffness of the T-BMT equaton on the error of the calculation (see appendix B).

IV Spin tracking simulations in an all-electric ring lattice

As established earlier, when the proton is at the magic energy (Eq. (20)), the spin is frozen in the horizontal plane along the momentum direction. Such a condition is met at one specific energy for a purely electric ring.

In what follows one simulated the all-electric Lebedev ring lattice in BMAD based on the storage ring EDM collaboration proposal lebedev sredm . The ring structure consists of 4 superperiods, each including 5 FODO cells with 3 cylindrical deflectors per half cell. An overview of the simulated ring is shown in fig 9 in the ideal case where the particle momentum vector is perfectly aligned with the spin vector. The main ring and beam parameters are summarized in table 1.

Fig 10 shows the Twiss parameters in the ring. In particular, the weak vertical focusing, $\beta_{y}^{max}=216$ m, was chosen in order to achieve enough sensitivity on the radial magnetic field. The aim of the remaining sections of this paper is to benchmark the previously established analytical formula against tracking simulations. The approach relies on computing the averages of the spin precession components on the closed orbit (no betatron motion considered) and comparing the results with the tracking simulations.

V Case with no longitudinal and no magnetic field components

In order to undertake the benchmarking analysis, one needs to reduce the complexity of the problem. For that reason, one first considers a lattice which is free of magnetic fields as well as any longitudinal electric field imperfections $E_{l}$ = 0. Thus, the spin precession components may be written:

[TABLE]

From the first order approximation (Eq. (39)), one can write for very short timescales:

[TABLE]

so that one can compute the rapidly oscillating terms for the vertical spin component:

[TABLE]

Now, for a fully eletrostatic ring, the particle equation of motion in the vertical plane may be written:

[TABLE]

which yields after integration

[TABLE]

Thus, on the closed orbit, $\left<E_{y}\right>=0$ and consequently $\left<\Omega_{r}\right>=0$. Injecting Eq. (91) into (89) yields:

[TABLE]

To the first order, the rapidly oscillating terms of the vertical spin components are mainly determined by the vertical particle slope.

V.1 Benchmarking of the analytical solution with numerical simulations for very small timescales

Now, let’s compare the analytical solution with the results of tracking simulations in a fully electrostatic ring in the case where $t\to 0$, i.e. for a few turns from injection: a particle is launched with an initial momentum offset $\Delta p/p=10^{-5}$ and tracked over several 1000 of turns to compute its trajectory. One defocusing quadrupole is vertically misaligned in the lattice by \Delta y=10\text{ \mum} thereby yielding a vertical orbit motion. The spin is computed by solving the Thomas BMT equation in the Frenet Serret coordinate system using a Runge Kutta integrator available in BMAD. From the orbit tracking results, one can compute the averages as well as the tilde of all quantities needed in Eq. (86).

Comparison of the tracking results with Eqs. (87) and (92) is shown in figures 11 and 12 respectively. The agreement is quite good for both quantities which demonstrates that the perturbation approach based on the method of averages is an adequate way to obtain an approximate solution of the T-BMT equation since all quantities evolve adiabatically in time. In particular, it is shown that the rapid oscillating terms for the vertical spin component are determined by the vertical angles of the particle. From Eq. (92), it follows that a first order approximation for very short timescales is given by:

[TABLE]

Note that the above analysis still holds if the particle is off the closed orbit. However, in that case, the averages should be calculated over several turns in order to remove the fluctuations of the spin precession components.

V.2 Study of misalignment errors

Next, one compares the previously established analytical formula with the tracking results without restriction to very short timescales. The Lebedev lattice lebedev is modified by introducing a vertical misalignment of one or several quadrupoles. The spin is only recorded after each turn completion. Thus, the benchmarking analysis with the analytical formula is based on the frozen spin solution. We discuss two cases, depending on the value of the average spin precession. Some of the following simulation results were discussed in previous publications selcuk_quad .

V.2.1 Case of a beam injected off the magic momentum

If the particle is injected with an initial momentum offset, then a vertical spin buildup will occur in presence of vertical motion. Assuming a net vertical misalignment of one quadrupole, one can search for the closed orbit ensuring the periodicity of the ring as illustrated in fig 13.

For instance, let’s consider the case of a particle injected with a momentum offset $\Delta p/p=10^{-5}$. Given that the particle is traveling through electrostatic elements, the total energy conservation is a crucial aspect of the simulation. The effect is most dominant inside the electrostatic bends mane leading to a non negligible change of the momentum offset of the particle as shown in fig 14.

The latter can be calculated in the following way:

[TABLE]

where the subscript $out$ denotes the value outside the electrostatic elements. From what preceded, one can calculate the average spin precession components and one obtains:

[TABLE]

so that the spin precession is dominated by the average of its vertical component $\left<\Omega_{y}\right>$. It results that the non-linear terms cannot be neglected, even for timescales of the order of few milliseconds. Thus, relying on the first order approximation (35), one can compute an approximate solution to the vertical spin buildup:

[TABLE]

Then, several particles are tracked with different initial momenta. Figs 15 shows a comparison of the longitudinal and radial spin results (obtained from BMAD tracking) with the analytical formula (96) and (97). As can be seen, the agreement is good.

Next, one varies the initial momentum offset as well as the quadrupole misalignment to determine its impact on the vertical spin buildup. From Eq. (98), a Taylor expansion yields:

[TABLE]

so that the vertical spin buildup is proportional to the momentum offset from the magic energy as well as to the vertical slope in the bends and therefore to the net vertical misalignment of the quadrupole. This is confirmed by tracking simulations as shown in fig 16. Nevertheless, the first order approximation fails at the magic energy, i.e. when $\Delta p/p\rightarrow 0$ since then, $\left<\Omega\right>\rightarrow 0$ and one needs to exploit the higher order terms in order to explain any spin buildup. This is discussed in the next section.

NB: In order to evaluate the above average quantities properly, one has to calculate the averages of all quantities by integrating over the closed orbit. For instance, a particle off-momentum yields a closed orbit that is different from the ideal one. Thus, a search for the closed orbit is necessary in order to avoid any fluctuations of the calculated quantities, since the above formalism is based on the assumption that all quantities are periodic so that the average value is unchanged from one turn to the next.

V.2.2 Case of a beam injected at the magic momentum

If a beam is injected at the magic momentum, then, in presence of misalignment errors, spin precession will occur. To show this, let’s consider a lattice that has two misaligned quadrupoles as follows: in the first quarter of the ring, a defocusing quadrupole is misaligned vertically and horizontally by $(+\Delta x,+\Delta y)$. In the third quarter, i.e. 180 degrees out of phase, a second defocusing quadrupole is misaligned by $(-\Delta x,-\Delta y)$. Thus, the average misalignment vanishes in this configuration. For a beam injected at the magic energy, an energy spread will appear at the entrance of the misaligned quadrupole. This results in a change of the closed orbit and consequently in a change of the kinetic energy within the electrostatic bends:

[TABLE]

This consequently yields a radial spin buildup within the electrostatic bends where the effect is most important, which generates a vertical spin as well by means of a longitudinal spin precession. Assuming $\Delta x=\Delta y=10\mu m$, one obtains:

[TABLE]

and

[TABLE]

Thus, the vertical spin buildup is mainly due to the geometric phases given by:

[TABLE]

which is proportional to the product of the displacements of both quadrupoles: the horizontal displacement of the quadrupoles yields larger radial spin oscillations due to the variation of the kinetic energy in the electrostatic bends while the vertical displacement of the quadrupoles yields a vertical slope inside the electrostatic bends, therefore a longitudinal spin precession which rotates the radial spin into the vertical plane. Such an effect yields a non-vanishing average value, therefore the frozen spin is proportional to both displacements as verified by tracking simulations in fig 17.

VI Magnetic field imperfections and EDM effect

In the case where the magnetic fields as well as the EDM signal are taken into account, the previous formalism holds, the only difference being the coefficients of the transfer matrix. The spin vector is still propagated with respect to the local reference frame of the accelerator. In other words, the presence of magnetic fields is considered as a perturbation that does not modify the reference frame: the latter remains attached to the ideal orbit in the median plane of the accelerator and although the horizontal restoring force contains a perturbation induced by the transverse magnetic fields, the accelerator reference frame does not rotate with such a perturbation. Such definition shall be taken with care though: in some derivations, one subtracts the angular velocity induced by the vertical magnetic fields $eB_{y}/(\gamma m)$. Given that the bending is only achieved by the electrostatic deflectors, we do not make such a correction in our calculations. Besides, successive quadrupoles in the lattice provide alternate bending in opposite directions to keep the beam stable so that, on average, such a correction may be neglected. An excellent note discussing this matter can be found in luccio .

Keeping terms up to the second order only, Eq. (24) re-writes:

[TABLE]

In what follows, one focuses on investigating the impact of the systematic errors originating from magnetic field imperfections. To simplify the analysis, one assumes a vanishing EDM signal, i.e. $\eta=0$. Unless otherwise specified, one relies on the second order approximation for the benchmarking of the analytical results with the tracking simulations.

VI.1 Geometric phases contributing to the vertical spin buildup

If the beam is injected at the magic energy, the geometric phase effects may become important. It is extremely difficult though to reach such configuration where $\left<\Omega_{r}\right>=\left<\Omega_{y}\right>=\left<\Omega_{l}\right>=0$. In general, one may reach a configuration where only the averages of two components of the spin precession vector vanish. The remaining component can however play an important role to transfer the spin buildup to another plane as is discussed in section VI.4. Table 2 summarizes the leading terms of the geometric phases as established in Eq. (73). Some of these contributions, particularly the third order terms cannot be eliminated by using counter-rotating beams. However, if the ring is properly shielded such as the integrated magnetic fields imperfections are below the nT.m level, then one can reasonably neglect such contributions to the vertical spin buildup in comparison with the aimed sensitivity of the experiment ($\sim nrad/s$ EDM signal). Precise orbit control is a condition sine qua non to eliminate such contributions mimicking the EDM signal.

VI.2 Horizontal magnetic field imperfections

The particle equation of motion in the vertical plane may be written:

[TABLE]

which yields after integration

[TABLE]

which vanishes on the closed orbit. This implies:

[TABLE]

Now, integrating the expression of the radial spin precession, and assuming no longitudinal electric fields one obtains:

[TABLE]

given that, at the magic energy $\left(a+\dfrac{1}{\gamma+1}\right)\beta_{l}^{2}\gamma=1$. Now, injecting Eq. (107) into (108) to compute the average on the closed orbit yields:

[TABLE]

where the magnetic fields are expressed in units of Tesla. For an aimed EDM sensitivity of $d_{p}=10^{-29}\text{e.cm}$ corresponding to $\eta=1.9*10^{-15}$, the resulting vertical spin buildup is

[TABLE]

so that an average radial magnetic field of $\approx 10\text{ aT}$ will generate the same signal as the smallest spin precession to be identified sredm .

In addition, the presence of any longitudinal magnetic fields can cause a vertical spin buildup and needs to be taken into account in the analysis. Therefore clearing the average radial magnetic fields is not a sufficient condition to eliminate the first order linear terms of the systematic errors arising from the radial spin precession.

VI.3 Longitudinal and vertical magnetic field imperfections

VI.3.1 Case with $\dfrac{\Delta p}{p}=0$

In the configuration where the particle is at the magic energy, the quadratic evolution with time is negligible since the average quantities of the spin precession vectors are negligible. Assuming that the imperfection is due to the presence of both vertical and longitudinal magnetic fields which are $90$ degrees out of phase as illustrated in fig 18 and for which the integrated fields are $\pm 1\text{ nT.m}$, one obtains:

[TABLE]

The third term is non vanishing due to an average radial spin precession such as $s_{r}(t)\approx\left<\Omega_{y}\right>t\approx 7.36*10^{-7}t$ which is mainly due to the variation of the energy within the electrostatic deflectors and is proportional to the horizontal orbit offset, therefore to the amplitude of the vertical magnetic fields. However, the second term is the dominant contribution to the vertical spin buildup. The order of the spin rotations is particularly important to explain the latter as shown below: $\tilde{\Omega}_{y}$ represents the integral of $\Omega_{y}-\left<\Omega_{y}\right>$ therefore accounts for the presence of any rapidly oscillating radial spin component. Such alternating magnetic field imperfections in the ring yield alternating spin precession vectors as illustrated in fig 19.

Next, one computes the contribution of the Berry phases to the vertical spin buildup:

[TABLE]

where $\tilde{\Omega}_{y}$ is defined by;

[TABLE]

so that

[TABLE]

The spin precession is a result of several localized imperfections in the ring. Using the hard edge model, one can express the perturbations as a sum of box distributions in the following way:

[TABLE]

where $L_{i}$ is the length of the perturbation located at $s_{i}$ and $\Pi(x)$ is the box distribution defined by:

[TABLE]

The integral of the box distribution is the stepwise function defined by:

[TABLE]

so that $\tilde{\Omega}$ can be expressed in the following way:

[TABLE]

Injecting this into Eq. (114) yields:

[TABLE]

where one assumed that both longitudinal and vertical perturbations do not overlap. $H$ is the Heaviside function which accounts for the fact that a vertical spin buildup arises after the radial one. In the case considered above this yields:

[TABLE]

which is proportional to the amplitude of the magnetic field perturbations. The tracking simulation results are finally summarized in fig 18 and 20 where one obtained a good agreement.

VI.3.2 Case with $\dfrac{\Delta p}{p}\neq 0$

In the case where $\Delta p/p=10^{-9}$, i.e. the particle is slightly above the magic energy, one obtains:

[TABLE]

A comparison of the tracking simulations with the analytical estimate (Eq. 119) is shown in fig 21. The dominant contribution is due to the third term which accounts for the linear slowly varying term of the radial spin component coupled with the rapidly oscillating term of the longitudinal spin precession.

In particular, such an effect sets a limit on the difference of energy of the two counter-rotating beams that can be employed to extract the EDM signal: for the case considered here with alternating magnetic field imperfections of 1nT.m, one obtains:

[TABLE]

In reality though, the amplitude of the vertical spin buildup is dependent on the location of the perturbation in the lattice, in particular on the shape of the integral of the longitudinal spin precession function, i.e. $\tilde{\Omega}_{l}$. To show this, let’s calculate the latter based on the Heaviside function model:

[TABLE]

Now, assuming that $\Omega_{l}$ consists of two alternating perturbations at locations $s_{1}$ and $s_{2}$ with the same length $L$, one can compute the integral:

[TABLE]

which shows that the amplitude of the vertical spin buildup due to the $\left<\Omega_{y}\right>\left<\tilde{\Omega}_{l}\right>$ is proportional to the distance between the two alternating perturbations of the longitudinal magnetic fields. To verify this result, one constructs several cases where the location of the perturbation is varied as shown in fig 22. Tracking simulations are performed for each case and one can see in fig 23 good agreement between the numerical simulations and the analytical estimates.

VI.4 Longitudinal and radial magnetic field imperfections

Assuming that the only perturbation for a particle at the magic energy is due to the radial and longitudinal magnetic fields such that $\left<\Omega_{y}\right>,\left<\Omega_{r}\right>\approx 0$ and $\left<\Omega_{l}\right>=4.3$ rad/s, it results from the second order approximation that the radial and vertical spin buildup are mainly due to the geometric phases such as:

[TABLE]

The simulated imperfections are shown in fig 24 and the comparison with the analytical results are summarized in fig 25. As can be seen the second order approximation is sufficient to explain the radial spin buildup. However, it fails to explain the quadratic spin buildup in the vertical plane. Thus, the third order approximation is invoked which yields:

[TABLE]

In summary, the geometric phases explain in part the radial spin buildup. However, the latter are transferred to the vertical plane by means of an average longitudinal spin precession.

VI.5 Longitudinal magnetic field imperfections

Assuming that the only perturbation is due to longitudinal magnetic field imperfections such that $\left<B_{l}\right>=0$, a vertical spin buildup can arise only if the particle is not at the magic energy which yields two contributions: a non-vanishing radial spin precession due to the horizontal orbit displacements coupled with the longitudinal magnetic fields, therefore $ds_{y}/dt=\Omega_{r}s_{l}\propto\left<x^{\prime}B_{l}\right>s_{l}$, and a non vanishing longitudinal spin precession coupled with the radial spin component. For instance, assuming an initial momentum offset $\Delta p/p=10^{-8}$ and $B_{l}L=\pm 10^{-7}\text{ T.m}$, it follows from Eq. (68) that:

[TABLE]

where $\tilde{\Omega}_{l}$ is given by:

[TABLE]

Such an effect is illustrated in fig 26 where one can see that, starting from point A, a linear radial spin component grows linearly with time. The latter yields a vertical spin buildup which depends on the location of the perturbations in the ring. This is verified by simulating two different lattices with different alternating field imperfections as illustrated in fig 27.

If the average contribution of the longitudinal magnetic fields is non-vanishing, then a quadratic increase of the vertical spin buildup with time is expected. Assuming $\Delta p/p=10^{-12}$ and an average longitudinal magnetic field of $1$ nT.m, one evaluates the vertical spin buildup by computing the averages:

[TABLE]

Comparison with the tracking simulations is shown in fig 28 where one can observe a good agreement between the two calculations.

VII Conclusion

Several formula were established and benchmarked against numerical simulations for the purpose of computing the systematic errors in an all-electric proton EDM storage ring. In particular, it appears that the established formalism is very useful to calculate the spin precession rates at the location of the polarimeter, i.e. at specific locations of the storage ring. To reach the sensitivity level required for such an experiment, the second order approximation based on the BKM method of averages is generally sufficient. In addition, it is very practical to benchmark different simulation codes for precision tracking by simply comparing the average quantities of the spin precession components and their periodic tilde functions.

Acknowledgements

We acknowledge useful discussions with members of the CPEDM and JEDI collaboration, specifically Mike Lamont, Sig Martin, Selcuk Hacıömeroğlu, Andreas Lehrach, Yannis Semertzidis, Ed Stephenson, Hans Stroeher and Richard Talman. Special thanks to David Sagan and Yann Dutheil for helping with BMAD and useful discussions as well.

Appendix A Spin Tracking example using Mathematica

\boldsymbol{\text{ClearAll}[\text{{``}Global\grave{}*{''}}]}

\boldsymbol{\text{SetOptions}[\text{$FrontEndSession},\text{PrintPrecision}\to\text{Ceiling}[\text{$MachinePrecision}]]}

$\boldsymbol{T=2.7872*10^{-6};}$

$\boldsymbol{\text{Norder}=7;}$

\boldsymbol{\Omega_{l}[\text{t_}]=-20*(\text{Sin}[8*\text{Pi}/T*t])-10^{-2};}

$\boldsymbol{\Omega_{\text{la}}=\text{Integrate}\left[\left.\Omega_{l}[t]\right/T,\{t,0,T\}\right]}$

\boldsymbol{\Omega_{\text{lt}}[\text{x_}]=\text{Integrate}\left[\Omega_{l}[t]-\Omega_{\text{la}},\{t,0,x\},\text{Assumptions}\to x\text{>=}0\&\&x<0.01\right];}

$\boldsymbol{\Omega_{\text{lta}}=\text{Integrate}\left[\left.\Omega_{\text{lt}}[t]\right/T,\{t,0,T\}\right]}$

\boldsymbol{\Omega_{\text{ltt}}[\text{x_}]=\text{Integrate}\left[\Omega_{\text{lt}}[t]-\Omega_{\text{lta}},\{t,0,x\},\text{Assumptions}\to x\text{>=}0\&\&x<0.01\right];}

\boldsymbol{\Omega_{r}[\text{t_}]=-2*10^{-10}*\text{Sin}[6*\text{Pi}/T*t];}

$\boldsymbol{\Omega_{\text{ra}}=\text{Integrate}\left[\left.\Omega_{r}[t]\right/T,\{t,0,T\}\right]}$

\boldsymbol{\Omega_{\text{rt}}[\text{x_}]=\text{Integrate}\left[\Omega_{r}[t]-\Omega_{\text{ra}},\{t,0,x\},\text{Assumptions}\to x\text{>=}0\&\&x<0.01\right];}

\boldsymbol{\Omega_{y}[\text{t_}]=20*\text{Sin}[10*\text{Pi}/T*t]+\text{Cos}[8*\text{Pi}/T*t]-10;}

$\boldsymbol{\Omega_{\text{ya}}=\text{Integrate}\left[\left.\Omega_{y}[t]\right/T,\{t,0,T\}\right]}$

\boldsymbol{\Omega_{\text{yt}}[\text{x_}]=\text{Integrate}\left[\Omega_{y}[t]-\Omega_{\text{ya}},\{t,0,x\},\text{Assumptions}\to x\text{>=}0\&\&x<0.01\right];}

$\boldsymbol{\Omega_{\text{lyta}}=\text{Integrate}\left[\Omega_{l}[t]*\left.\Omega_{\text{yt}}[t]\right/T,\{t,0,T\}\right]}$

\boldsymbol{\Omega_{\text{lytt}}[\text{x_}]=\text{Integrate}\left[\Omega_{l}[t]\Omega_{\text{yt}}[t]-\Omega_{\text{lyta}},\{t,0,x\},\text{Assumptions}\to x\text{>=}0\&\&x<0.01\right];}

************* second order frozen solution for the radial spin **************

*************        $\xi_{r,2}(t)=(r1+r2+r3)*t+r4*t^{2}$

$\boldsymbol{\text{r1}=\Omega_{\text{ya}}}$

$\boldsymbol{\text{r2}=\text{Integrate}\left[\Omega_{l}[t]*\left.\Omega_{\text{rt}}[t]\right/T,\{t,0,T\}\right]}$

$\boldsymbol{\text{r3}=-\Omega_{\text{ra}}*\Omega_{\text{lta}}}$

$\boldsymbol{\text{r4}=\Omega_{\text{ra}}*\left.\Omega_{\text{la}}\right/2}$

************* second order frozen solution for the vertical spin **************

*************        $\xi_{y,2}(t)=(t1+t2+t3)*t+t4*t^{2}$

$\boldsymbol{\text{t1}=-\Omega_{\text{ra}}}$

$\boldsymbol{\text{t2}=\Omega_{\text{lyta}}}$

$\boldsymbol{\text{t3}=-\Omega_{\text{ya}}*\text{Integrate}\left[\left.\Omega_{\text{lt}}[t]\right/T,\{t,0,T\}\right]}$

$\boldsymbol{\text{t4}=\Omega_{\text{ya}}*\left.\Omega_{\text{la}}\right/2}$

$\boldsymbol{\Omega_{\text{tot}}=\left(\Omega_{\text{ra}}{}^{2}+\Omega_{\text{ya}}{}^{2}+\Omega_{\text{la}}{}^{2}\right){}^{0.5}}$

$\boldsymbol{\text{sol}=\text{NDSolve}[\{s_{r}^{\prime}[t]==\Omega_{y}[t]*s_{l}[t]-\Omega_{l}[t]*s_{y}[t],s_{y}^{\prime}[t]==\Omega_{l}[t]*s_{r}[t]-\Omega_{r}[t]*s_{l}[t],}$

$\boldsymbol{s_{l}^{\prime}[t]==\Omega_{r}[t]*s_{y}[t]-\Omega_{y}[t]*s_{r}[t],s_{y}[0]==s_{r}[0]==0,s_{l}[0]==1\}},$

$\boldsymbol{\{s_{r}[t],s_{y}[t],s_{l}[t]\},\{t,0,400.1*T\},\text{MaxStepSize}\to 10^{-9},\text{MaxSteps}\to 8000000,}$

$\boldsymbol{\text{Method}\to\{\text{{``}ExplicitRungeKutta{''}},\text{{``}DifferenceOrder{''}}\to\text{Norder}\}];}$

************* solution for the vertical spin using explicit Runge Kutta tracker **************

\boldsymbol{\text{sytrack}[\text{t_}]=s_{y}[t]\text{/.}\text{sol};}

\boldsymbol{\text{frozsytrack}[\text{t_}]=\text{sytrack}[\text{IntegerPart}[t/T]*T];}

************* second order analytical solution for the vertical spin **************

*************        $\xi_{y,2}(t)+\phi_{y,2}(t)$

\boldsymbol{\text{syanal}[\text{t_}]=(\text{t1}+\text{t2}+\text{t3})*t+\text{t4}*t^{2}-\Omega_{\text{rt}}[t]+\Omega_{\text{ya}}*\left(\Omega_{\text{lt}}[t]*t-\Omega_{\text{ltt}}[t]\right)+\Omega_{\text{lytt}}[t];}

\boldsymbol{\text{frozsyanal}[\text{t_}]=(\text{t1}+\text{t2}+\text{t3})*t+\text{t4}*t^{2};}

$\boldsymbol{\text{plotsy}=\text{Plot}[\{\text{sytrack}[t],\text{syanal}[t],\text{frozsytrack}[t],\text{frozsyanal}[t]\},\{t,0,1.1*T\}]}$

************* compute and plot the error of the analytical approximation **************

\boldsymbol{\text{error}[\text{t_}]=\text{sytrack}[t]-\text{syanal}[t];}

\boldsymbol{\text{frozerror}[\text{t_}]=\text{frozsytrack}[t]-\text{frozsyanal}[\text{IntegerPart}[t/T]*T];}

$\boldsymbol{\text{errorplot}=\text{LogLogPlot}[\{\text{Abs}[\text{error}[t]],\text{Abs}[\text{frozerror}[t]]\},\{t,1.1*T,400.1*T\}]}$

Appendix B Stiffness studies using explicit Runge Kutta tracker

When hard edge models are encountered, stiffness becomes the bottleneck of precision tracking. When solving the equation of motion and/or the T-BMT equation, some singularities may arise at the entrance/exit of the hard edge elements due to discontinuities in the fields calculated or their derivatives. In order to assess its impact on the convergence of the Runge Kutta tracker, one simulated several cases in Mathematica where one started from a smooth model with continuous functions and approximate it by stepwise functions for which the stiffness can be controlled. For instance, increasing the stiffness can be achieved by reducing the width of the stepwise function, therefore increasing the slope at the entrance/exit of the element (the hard edge model corresponds to an infinite slope). This is better illustrated in fig 29. The idea is to investigate the convergence of the Runge Kutta tracker by slowly increasing the stiffness of the system. One first starts from smooth functions describing the spin precession components as sinusoidal oscillations. Then, using the above model, one can vary the stiffness of the system: from one case to the next, one reduces by half the rising time of the spin precession components. The T-BMT equation is then solved for each case using a specific order and a fixed step size of the explicit Runge Kutta tracker in Mathematica. The relative error after 10 revolution periods is computed for each case corresponding to a fixed integration step size, a specific order of the Runge Kutta method and a given spin precession function. In one case the assumed reference solution corresponds to the $9th$ order Runge Kutta method with the smallest step size. In the second case, the assumed reference solution corresponds to the second order analytical solution. The results are displayed in fig 30 where one can observe some oscillatory behavior of the relative error due to instability of the Runge Kutta method. However, in general, the higher the stiffness, the larger the error becomes for a given step size. Nevertheless both the analytical and the Runge Kutta solutions converge towards each other. Using implicit Runge Kutta tracker does reduce the error. This will be discussed in later publication.