Geometric phases in electric dipole searches with trapped spin-1/2 particles in general fields and measurement cells of arbitrary shape with smooth or rough walls

The important role of geometric phases in searches for a permanent electric dipole moment of the neutron, using Ramsey separated oscillatory field nuclear magnetic resonance, was first noted by Commins and investigated in detail by Pendlebury et al. Their analysis was based on the Bloch equations. In subsequent work using the spin density matrix Lamoreaux and Golub showed the relation between the frequency shifts and the correlation functions of the fields seen by trapped particles in general fields (Redfield theory). More recently we presented a solution of the Schr\"odinger equation for spin-$1/2$ particles in circular cylindrical traps with smooth walls and exposed to arbitrary fields [Steyerl et al.] Here we extend this work to show how the Redfield theory follows directly from the Schr\"odinger equation solution. This serves to highlight the conditions of validity of the Redfield theory, a subject of considerable discussion in the literature [e.g., Nicholas et al.] Our results can be applied where the Redfield result no longer holds, such as observation times on the order of or shorter than the correlation time and non-stochastic systems and thus we can illustrate the transient spin dynamics, i.e. the gradual development of the shift with increasing time subsequent to the start of the free precession. We consider systems with rough, diffuse reflecting walls, cylindrical trap geometry with arbitrary cross section, and field perturbations that do not, in the frame of the moving particles, average to zero in time. We show by direct, detailed, calculation the agreement of the results from the Schr\"odinger equation with the Redfield theory for the cases of a rectangular cell with specular walls and of a circular cell with diffuse reflecting walls.


I. INTRODUCTION
The experimental search for a permanent electric dipole moment (EDM) of the neutron or other particles is motivated by its potential to search for deviations from the Standard Model. The most sensitive experimental approach so far, for neutrons, uses ultracold neutrons (UCNs) trapped in a measurement cell together with comagnetometer atoms like 199 Hg. They are exposed to a weak static magnetic field B 0 and a strong electric field E applied parallel or anti-parallel to B 0 . Employing the Ramsey method of separated oscillatory field magnetic resonance the experiment is designed to search for small changes of the Larmor frequency due to the presence of the electric field.
At the high level of sensitivity already achieved in previous EDM experiments [1][2][3] and expected to be improved in current EDM projects (e. g., [5][6][7][8][9]) a significant source of potential systematic error appears to be geometric phases. Their role in EDM experiments was first noted by Commins [4] and analyzed by Pendlebury et al. [3] on the basis of the Bloch equation for spin evo- * Electronic address: asteyerl@mail.uri.edu lution in a time-dependent magnetic field. Subsequent work by Lamoreaux and Golub [10] and Barabanov et al. [11] used the Redfield approach [12][13][14] based on the spindensity matrix and extended the scope of such studies to general NMR physics. The original magnetic field model with a constant vertical gradient [3] of the Larmor field was extended in Refs. [15][16][17][18] to include more general macroscopic fields and the microscopic field of a magnetic dipole.
In Ref. [19] we showed that the frequency shifts due to geometric phases can also be determined by directly solving the Schrödinger equation to second order in the perturbation. Solutions were presented for circular cylindrical traps with smooth walls and arbitrary macroscopic magnetic fields as well as for the microscopic field of a magnetic dipole.
The direct solution of the Schrödinger equation includes vertical spin oscillations and the non-stationary, transient spin dynamics, i.e., the gradual development of the shifts for an arbitrary number of wall reflections subsequent to the start of free spin precession in the Ramsey scheme.
In the present article we extend this approach to include diffuse reflectivity from rough walls and to vertical cylindrical measurement cells with arbitrary cross section. Their top and bottom are assumed to be perfectly flat. In Sec. II we develop an analytic solution of the Schrödinger equation for perfectly diffuse reflection, in the xy plane only, on the cylindrical walls, both for a constant gradient field in circular cylindrical geometry and for a more general field model where B 0 and E are not exactly aligned.
Our work is the first to address directly the effects of diffuse wall reflections in different geometries. In Sec. III we consider a generic cross section of the measurement cell and calculate the frequency shifts (that linear in E and those quadratic in electric field or magnetic field) for cells in the form of a rectangle with arbitrary aspect ratio. Rectangular cells are foreseen for EDM project [7]. Our results are valid over the full range of in-plane particle velocity v for given Larmor frequency ω 0 . We compare them with the simulations described in Sec. IV and observe consistency between the two approaches. The nonadiabatic and adiabatic limits of frequency shift agree with the universal predictions of Refs. [14,18,21]. In Sec. V we establish the equivalence of our method with other approaches [10,13,15]. The agreement is notable in view of the different sets of assumptions made by each.
The usual approach to these problems, the Redfield theory of frequency shifts and relaxation based on the von Neumann equation for the density matrix, makes use of a plethora of assumptions [12] and has resulted in a considerable literature concerning its range of validity, e. g. [22]. The density matrix represents the effects of both the quantum mechanics and statistical fluctuations of the system under study. In our approach, which involves first solving the Schrödinger equation, to second order in the perturbation, for the wave function for an arbitrary time varying field and then applying the statistical averaging to the result, there is a clear separation between the quantum and statistical effects. In fact our solution can be easily applied to cases where the Redfield result no longer holds, such as observation times shorter than the correlation time, and non-stochastic systems.

A. Fundamentals
In the EDM experiments the cell walls are never perfect mirror surfaces, thus over the course of some 10 1 to 10 2 consecutive reflections randomization of trajectories is established for any initial distribution. The total number of wall reflections occurring during the measurement period in EDM experiments is much larger.
In the present article we neglect scattering on gas molecules and consider perfectly diffuse surface scattering in the xy plane perpendicular to the vertical symmetry axis of the cell (Fig. 1), thus it suffices to analyze the particle motion as projected onto the xy plane [3]. The in-plane velocity v is constant and the scattering probability distribution is ∝ cos ρ (Lambert's law with 000 000 000 000 000 000 000 000 000 000 000 000 000 in-plane scattering angle ρ), independent of incident angle ρ inc and, therefore, independent of the path history. Even a single UCN reflection will establish isotropy of flight direction. This is also a reasonable model for 199 Hg co-magnetometer atoms which stick to the wall briefly before re-emission in a random direction. Summarizing previous results [19] we start from the Hamiltonian for a classical field, where σ x , σ y and σ z are the Pauli matrices, ω 0 = −2µB 0 / is the Larmor frequency in the vertical static magnetic field B 0 and µ is the magnetic moment of the particle. The measurement cell is also exposed to a uniform vertical electric field E which gives rise to the motional magnetic field E × v/c 2 . In Eq. (1), Σ(t) = ω 0 (B x (t) + iB y (t)) /B 0 is the timedependent perturbation due to the small horizontal magnetic field (B x , B y ) seen by the particle moving from point a to point b (Fig. 1) at constant v. In the case with small uniform vertical gradient ∂B z /∂z, first analyzed in [3], the perturbing field is given as B x = −B 0 (ηΩ + ζx/R), B y = −ζB 0 y/R, where we use the dimensionless parameters Ω = v/(Rω 0 ), ζ = (R/(2B 0 )) ∂B z /∂z and η = Rω 0 E/(B 0 c 2 ).
For convenience we define a dimensionless time scale τ = ω 0 t, reset the clock at the center of each chord and solve the Schrödinger equation for Hamiltonian (1) in a coordinate system rotating about the vertical z axis at the Larmor frequency ω 0 . For B 0 pointing up, ω 0 > 0 (< 0) for µ < 0 as for neutrons (for µ > 0 as for 199 Hg). Using a second-order perturbation approach we obtain the spin-up and spin-down components, α r (τ b , τ a ) and β r (τ b , τ a ), of the wave function for a particle starting at point a at time τ a = −δ i with spin up and measured at b at τ b = +δ i : where δ i = (sin α g,i )/Ω and (see Eq. (9) of [19]) For circular cylindrical symmetry, vertical electric and magnetic fields E and B 0 and uniform vertical gradient ∂B z /∂z, the functions ν and β r in (2) are given by as shown in [19], Eq. (42). In (4), u 1 = u + 1, u = (η/ζ) + (cos α g,i )/Ω. We write the wave function at the end of n consecutive chords in the form where the initial spinor ψ (0) points in an arbitrary direction given by angles θ 0 (polar) and Φ 0 (azimuthal). In the experiments we have θ 0 π/2. The recursion relation linking chains with n−1 and n segments is (Eq. (33) of [19]) where R n = e i(αg,n−1+αg,n)/2 0 0 e −i(αg,n−1+αg,n)/2 is the rotation matrix, around the z axis, for angular change α g,n−1 + α g,n between the coordinate systems of segments n−1 and n. For each chord the y-axis is chosen to lie in the direction of motion. In (5) and (6) we have left the starting position on chord 1, given by τ a,1 , arbitrary. The asymptotic frequency shift, for n 1, does not depend on this position.
Choosing point a (Fig. 1) as the starting point of chord 1, we set τ a,1 = −δ 1 and obtain for the initial transfer matrix M (1) The unitary evolution matrix along any chord i, in the rotating system, reads and the transformation matrix from the laboratory to the rotating system at point a, is the same as that for the reverse transition at point b.
When we perform the matrix multiplications in (6) repeatedly for consecutive chords we obtain a sequence of matrices M (n) r (δ n , −δ 1 ) whose general form can be inferred by induction. Just as for a smooth wall the result is a unitary transfer matrix of the form for any n ≥ 2. As shown in [19], Eq. (34), the cumulative phase advance at the end of chord n, relative to the rotating frame and averaged over initial azimuthal spinor angle Φ 0 relative to initial flight direction, is given solely by the argument of g r : where the last form in (13) holds since g r differs from unity only slightly. For any n ≥ 2 we obtain from recursion relation (6) where In general, reversing the order of two successive chords leads to a different phase advance. This property is reflected in the form of g * r . Expression (14) is valid for any n 0 , n 1 in the range 1 ≤ n 0 < n 1 ≤ n.
Using (14) the overall phase shift (13) reads with × e i αg,n 1 −δn 1 +αg,n 0 −δn 0 +2 Dividing δϕ 0→n by the elapsed time T 0→n we obtain the average frequency shift δω 0→n . The averaging is subtle as described below in Sec. II B after a closer look at the definition of angular advance 2α g around the cell per chord.
In Eq. (7) we have assumed that, going from chord n−1 to n, the net rotation angle to the new direction of travel can be written α g,n−1 + α g,n . This is justified if, as in [3,19], we define α g as positive for passage on the right of the center O ( Fig. 1) (corresponding to counterclockwise rotation around the cell), as negative for passage on the left (for clockwise rotation), and if both chords pass O on the same side (which is always true for specular reflection and circular geometry).
However, this definition for α g (which we will later call α g,old ) fails if in case of diffuse scattering the consecutive chords pass on opposite sides of O. If so the final direction of travel is opposite that given by α g,n−1 +α g,n . This is evident for the special case α g,n = −α g,n−1 where the particle reverses its direction at the wall reflection point and thus has experienced a net rotation by angle π rather than by α g,n,old + α g,n−1,old = 0.
We avoid this inconsistency as follows: For passage on the right we define α g as before: for circular geometry, α g = ρ + π/2, for ρ < 0.
(18) can also be written α g = α g,old + π. With this choice of a continuous function for α g , instead of the old function which changes abruptly from +π/2 to −π/2 at ρ = 0, the net rotation angle can be expressed as α g,n−1 + α g,n also for passage on opposite sides. With this choice there is an offset by 2π if the next chord n passes on the left, irrespective of the geometry of chord n−1. As a result of this offset, over the course of multiple consecutive chords a cumulative offset angle 2mπ may appear, with integer m, but this offset is irrelevant both for the analysis and the simulations presented below since these depend only on the sine and cosine of angles. The choice of a continuous function for α g throughout the entire range from 0 to π, as in (17), (18), is also appropriate for the general cell geometries discussed below in section III.

B. Two averaging procedures
As shown in [3] for circular cylindrical geometry and specular reflection, the frequency shift X(α g , Ω) is averaged over all possible chord angles α g using a weight factor sin 2 α g , thus X(α g , Ω) 1 = 4 π π/2 0 dα g (sin 2 α g )X(α g , Ω). (19) We use subscript 1 to distinguish this type of averaging from a second type introduced below in (22). First we show that the averaging scheme (19) can also be derived differently and adapted to arbitrary cell geometry and rough surfaces satisfying Lambert's law.
For circular geometry the take-off angle ρ shown in Fig. 1 is complementary to α g : −ρ = π/2 − α g , implying cos ρ = sin α g . The quantity to be averaged is a frequency, i.e. the phase δϕ accumulated along the chord divided by the traversal time L/v, where L is the chord length. Since v = const., the frequency is determined by the phase per unit length, δϕ/L. A longer chord represents a larger sample for averaging, thus the chord's contribution to the mean frequency is proportional to its length L. Including the Lambert factor cos ρ we obtain the weight factor p 1 (L, ρ) ∝ L cos ρ. For circular geometry L = 2R cos ρ, thus p 1 (ρ) ∝ cos 2 ρ = sin 2 α g in agreement with the weight factor in (19).
Changing from integration over α g in (19) to integration over ρ we have where the integral covers the full scattering range −π/2 < ρ < +π/2. The normalization constant is N 1 = π/2. The measured phases are averages over both senses of particle circulation around the cell, forward (fw) and backward (bw). Eq. (20) represents the part for forward motion (in Fig. 1, from a to b). The corresponding contribution for backward motion from b to a is where we replace ρ by ρ − π and Ω by −Ω to transform to the opposite direction of motion. The fwbw even and odd contributions to the mean values are ... even = (1/2) ... fw + ... bw and ... odd = (1/2) ... fw − ... bw , respectively.
Certain averages in (15) do not involve relative quantities like the frequency (given by the phase per unit length) but properties of the chord as a whole, such as its length L or angles α g and ρ or the traversal time L/v = 2δ/ω 0 . The contribution of each chord to these averages is independent of its length L, thus the weight factor changes to p 2 (ρ) ∝ cos ρ = sin α g . For these quantities the equivalence of averaging procedure (20), (21) is, for circular geometry, Y (ρ, Ω) 2,fw = Y (ρ, Ω) 2,even + Y (ρ, Ω) 2,odd with N 2 = 2. For instance, δ 2 = δ 2,even = π/(4Ω).
In summary, we use averaging 2, given by (22) for circular geometry, to average properties of the chord as a whole, such as δ or β r (δ, −δ) or e 2i(αg−δ) (as in Eq. (27)). We average relative quantities such as frequency or (sin δ)/δ using method 1 which is given by (20), (21) for circular geometry.
In Sec. III we generalize these averaging prescriptions to a generic cylindrical cell geometry with arbitrary cross section.
C. Determining the mean frequency shift
In (14) the functions β r appear in the form of products β r,n0 β * r,n1 where n 0 and n 1 > n 0 refer to different chords. The following analysis is significantly simplified by noting that for completely diffuse scattering the chords are statistically independent of one another. There is no memory of the previous path. This is an essential point and argued as follows. For any trajectory emerging from a given surface element, such as at point a in Fig. 1, the 2D Lambert law ensures isotropy of take-off direction ρ, independent of the incident angle. While this is true for any cell geometry, the circular geometry is unique in ensuring a second type of uniformity which is also required to render any two generations of chords (labeled by indices i and j) independent of one another: a circular surface is illuminated uniformly by a diffuse point source positioned anywhere on the surface. Rays leaving point a at angle ρ within dρ intercept a surface strip of width ds = L dρ/ cos ρ where L = 2R cos ρ is the chord length. Thus ds = 2R dρ is independent of impact position. In consequence, at any step in any sequence of chords the wall is a uniform source of next generation chords. The orientation of the previous chord does not matter. We will see in Sec. III that this is not true for geometries different from circular and in those cases we will have to resort to approximations. We obtain an exact solution only for circular geometry.
In the latter case, the average of a product of chord characteristics such as β r , α g , δ for different chords i and j, or of any functions F i , G j thereof, can be factorized: In the case at hand, where n 1 = n 0 , this applies to the product of functions F n0 = β * r,n0 e i(αg,n 0 −δn 0 ) and G n1 = β r,n1 e i(αg,n 1 −δn 1 ) in (14). Since all chords are equivalent we obtain the same averages for any n 0 and n 1 : where we use averaging scheme 2 since these functions refer to a chord as a whole, as defined before (22).
The same product rule applies to the terms in the sum over k in (14) where n 0 < k < n 1 .
. Expressed in these terms the second contribution (16) to the overall phase (15) is given by the imaginary part, Q, of Each average in the last sum of terms is equal to D. Therefore this sum is geometric: where we have used Eq. (26). The second term in the square brackets is negligible for n 1 since the magnitude of D is always significantly less than 1. A similar negligible correction appears if the particle starts on chord 1 at a point different from a ( Fig. 1). In the EDM experiments using UCNs [1][2][3] the starting point is given by the beginning of the period of free spin precession in the Ramsey scheme and the particle can be at any point of the initial chord at this time.
In (31) we use the arithmetic mean value of the two contributions in the numerator of Q 1 . To show that this is justified we employ the following relation between the two averaging procedures 1 and 2: For any function F of chord variables we have For circular cell geometry this follows directly from the averaging definitions (20) - (22). Since δ = (cos ρ)/Ω, we have with normalization constants N 1 = π/2, N 2 = 2, and both for fw and bw motion. This proves (32). Thus, the two contributions in the numerator of Q 1 in (31) are equal and it was justified to use their arithmetic mean value. Relation (32) also holds for the generalization of averaging procedures 1 and 2 to general cell geometry introduced below in (44) -(47) of Sec. III A 2.
Using (32) the expression for the fw-bw averaged frequency shift (31) becomes Eq. (35) can also be understood as: mean frequency shift = mean net phase shift divided by mean net time.

Comparison between specular and diffuse reflection
Equation (35) is valid for perfectly diffuse 2D wall reflectivity. Comparing it with expression (24) for specular reflectivity and noting identity (32) we see that the first term on the right hand side (RHS) of (35), which represents the single chord contribution, is the same as in (24), thus independent of surface roughness.
The second term in (35) ensures that the spinor remains unchanged at each wall reflection. In the limit of specular reflection it reduces to the second term in (24) as follows: In this case α g is constant along an orbit, thus D = e 2i(αg−δ) . The product of two averages in (35) is replaced by the average of the single function β * r β r e 2i(αg−δ) of α g over a random distribution of α g -values. Writing , applying (32) and averaging over the fw and bw directions we recover (24) from (35).
Prior to comparing numerical results for the frequency shifts we extend the constant gradient field, considered so far, to include the possibility of a slight misalignment between the vertical E field and the Larmor field B 0 . This type of imperfection of order 10 −3 rad or more is hardly avoidable in the experimental situation.

D. Extension to general gradient field
In Refs. [15,18,19] the uniform gradient field has been extended to more general macroscopic fields B = ∇χ derived from a power series expansion of a magnetic potential χ(x, y, z), including terms of higher than second order, e.g. xy, x 4 , x 3 y, etc. Since the functions ν(δ, −δ) and β r (δ, −δ) can be determined analytically as in [19], we can use Eq. (35) to find analytic solutions also for 2D diffuse reflection.
We discuss briefly the effect of tilting the Larmor field slightly away from vertical by adding to B 0 a small uniform horizontal component B ρ0 = (B x0 , B y0 ). Slight misalignment between the volume averaged Larmor field and the symmetry axis z is inevitable in the experiments and has been taken into account in [10,23] with the result that the precession about the tilted z axis at frequency ω 0 = ω 0 1 + B 2 ρ0 /B 2 0 gives rise to the same frequency shift relative to ω 0 as for a vertical Larmor field with mag-nitude given by the field magnitude B L = B 2 0 + B 2 ρ0 . We can derive this result by exactly solving the Schrödinger equation for Hamiltonian (1) in the xyz laboratory frame. The two spinor components α(t), β(t) are determined from [19] The solutions , for initial spin up, describe uniform precession about the z axis in the xyz frame. In the limit of long time of free precession, ω 0 t/(2π) 1, which is always realized in the experiments, the precession angles about the z and about the z axis become equal and are given by the field magnitude As a result, all results for frequency shifts relative to ω 0 , presented below for purely vertical volume average B 0 of the field, also hold for a field with slightly tilted volume average, except that in this case the shift is measured relative to ω 0 . The redefinition of Ω in terms of ω 0 rather than ω 0 , Ω → v/(Rω 0 ), implies a negligible scale change in all plots of frequency shift vs. Ω.
In Figs. 2 to 4 the specular data are the same as those of [3,19] except that in Fig. 2 of [19] the dependence on finite n had been shown explicitly and the B 2 shift had not been averaged over the fw-bw senses of circulation. The sharp resonances for specular reflection appear blurred since we used a slight imaginary offset of α g , away from the poles, in the numerical integration over α g .
The limiting behavior for Ω 1 is shown more clearly in the right panels of Figs shifts ∝ EB and ∝ E 2 have been divided by Ω 2 . The lower panels of Fig. 2 show the agreement, within numerical uncertainty, of E-odd shift values calculated from Eq. (35) with the simulations described below in Sec. IV. We observed full agreement also for the second-order shifts ∝ E 2 and ∝ B 2 . In Fig. 4 the asymptotic form at large Ω is shown more clearly in the right panels, where the functions have been multiplied by Ω 2 .
Comparing the data for smooth walls and rough walls, the sharp resonances for specular reflection are largely washed out by roughness, as expected and previously shown by computer simulations in [3].

III. EXTENSION TO GENERAL CELL SHAPE
A. Generic cell cross section Figure 5 shows a cell of generic cross section. As in Fig. 1 we show a trajectory of length L extending from point a to point b where a is at position s measured along the perimeter from an arbitrary starting point s = 0. The chord has the perpendicular distance σ from the cross sectional center-of-mass point CM, thus σ is the variable replacing the closest distance x = R cos α g for circular geometry. σ is positive (negative) for trajectories passing the CM on the right (left).
As in Fig. 1, 2α g is the angle subtended by the chord L as seen from the CM. It equals the angular advance around the cell when the particle travels from a to b.
For circular geometry we had arbitrarily chosen the radius R as a convenient reference length for the definition of Ω, ζ and η: Ω = v/(Rω 0 ), ζ = R(∂B z /∂z)/(2B 0 ) and η = Rω 0 E/(B 0 c 2 ). For the rectangular cell shape ana- lyzed below (and also for the equilateral triangle analyzed in [20]) it is convenient to choose one half of the average side length as the unit of length. Thus, σ and L will denote normalized lengths which, for circular geometry, are given as σ = cos α g = − sin ρ and L = 2 sin α g = 2 cos ρ. As before we measure the wall intersection angle ρ from the surface normal and define it as positive for the counterclockwise direction (which corresponds to clockwise circulation of the particle around the cell). The trajectory at angle ρ CM , which passes through the center point, divides the region on the right of CM from that on the left. We transform from forward to backward direction of motion by replacing +Ω by −Ω and ρ by ρ − π. This operation reverses the direction of travel from a to b (in Figs. 1 and 5) to that from b to a, retaining the sign of σ and δ, but reversing the sign of α g (and adding an offset angle π, see (18)). The elapsed time is positive for either direction of motion. These symmetries are the same as for circular geometry. In general, the times at the chord endpoints, τ b = δ b (> 0 for the chord shown) and τ a = −δ a (< 0), will not add to zero. As for the circular case we define the average as δ = (δ b + δ a )/2. The imbalance is δ ba = (δ b − δ a )/2. We also distinguish between two angles subtended from σ s s + ds FIG. 5: (Color online) Measurement cell of generic shape subject to the restriction that any chord intersects the perimeter at two points only. The initial point a of chord ab is at location s which is measured along the perimeter. As for circular geometry, 2αg designates the angle subtended by the chord as seen from the areal center of mass CM but the triangle a b CM is, in general, asymmetric and the wall intersection angle ρ is no longer directly determined by αg. L and σ are the length and "lever arm" of the chord. L is divided into two sections of, generally, different length: from a to the epicenter c and from c to b with associated angles αg,a and α g,b .
the center: α g,a for the chord segment from a to c, and α g,b = 2α g − α g,a for that from c to b.
The corresponding change of the transfer matrix element g * r is taken into account in Eqs. (16), (29), (30), (31) and (35) by the replacements In the limit n 1 Eq. (30) becomes and the fw-bw averaged frequency shift takes the form (43) As in Eq. (35) for circular geometry, the first term in (43) is the single chord contribution. The second term ensures that the spinor remains unchanged at wall reflections. The two terms are averaged using the adaptation of methods 1 and 2 to general cell geometry described in the next section. We will also describe method 3 for averaging the product F G in (43), with F = β * r e iY0 and G = β r e iY1 . For circular geometry this product had been averaged as in (35). Its extension to general cell geometry is denoted (F )(G) 3 .

Averaging procedures for general 2D cell geometry
For circular cell geometry the surface coordinate s is irrelevant and does not appear in the averaging procedures, (20), (21) for method 1 and (22) for method 2. In an extension to general geometry we make use of the fact that, at equilibrium, all wall reflection points s and chord directions are equally probable, as for a 2D gas. Therefore, at any time the probability, per unit length and unit time, for a trajectory to emanate from surface position s is the same and the take-off angle ρ obeys Lambert's distribution law (cos ρ)dρ. where A is the cell area and P is its perimeter. Expression N 1 = 2πA holds since, for reversed order of integration in (46), the integral ds (cos ρ) L(s, ρ) represents the total area A for all parallel strips of width dσ = ds cos ρ, extending from the boundary on the left to that on the right in Fig. 5. The factor 2π in (46) is the result of integration over a uniform distribution of strip directions.
To average the product β * r e iY0 β r e iY1 in Eq. (43) we have to take account of the fact that consecutive chords i and i + 1 are not independent of one another. The endpoint of chord i (which is at wall position s) coincides with the starting position of chord i + 1.
We can quantify this correlation as follows: Consider trajectories starting from a given point s in random direction. Except for the singular case of circular geometry discussed in Sec. II C 2, the next wall impact point will be unevenly distributed over the wall surface and this nonuniformity extends over a number of consecutive reflections. This number increases with higher cell asymmetry, as for a rectangle with larger aspect ratio.
As a result, the averaging method used in Eq. (35) will be exact only for circular geometry. This is borne out by comparison of results for non-circular cells, based on (35), with the simulations described below in Sec. IV. It turned out that a closer approximation was obtained as follows. We average the product of functions F (s, ρ, Ω) = β * r e iY0 and G(s, ρ, Ω) = β r e iY1 adopting a third method where we treat all pairs of chords as if they were consecutive with a common wall collision point s. This model leads to with N 3 = 4P . In Eq. (48) we average over the fw-bw directions using the fw→bw transformations Ω → −Ω, ρ → ρ − π, σ → σ, (α g , α g,a , α g,b ) → (−α g + π, −α g,b + π, −α g,a +π), (δ, δ a , δ b ) → (δ, δ b , δ a ) as for a time-reversed sequence of chords. As always, we obtain the frequency shift by dividing the net phase by the elapsed time and obtain Eq. (43) via identity Eq. (32).
Rectangular measurement cell with side lengths 2(1 + w) and 2(1 − w), thus aspect ratio (1 + w)/(1 − w). Among four groups of paths originating from side IV we show a path of group 1 with endpoint on side I (see the text). Figure 6 shows a rectangular cell with side lengths 2(1 + w) and 2(1 − w), where w is in the range 0 ≤ w < 1 and determines the aspect ratio (1 + w)/(1 − w). We use one eighth of the circumference as the unit of length. Referring to the analysis of this system in Appendix A we use Eqs. (A1) -(A9) in approximation (43) for the frequency shifts and in the simulations of Sec. IV and plot the results in Figs. 7 and 8 for w = 0 (square cell) and w = 0.5 (rectangle with aspect ratio 3). In Fig. 7 the values obtained from Eq. (43) for w = 0 and diffuse wall scattering are seen to agree with the simulation data within a few percent.

B. Rectangular measurement cell
In Fig. 8, for w = 0.5, we also compare simulations and analysis for specular walls with those for rough walls. While the adiabatic and non-adiabatic limits of the simulations are independent of roughness [3,10,11,18,21], the resonances seen in the smooth-wall data are reduced by roughness. For the E-odd shift, the diffuse wall approximation (43) deviates by 6% (33%) for w = 0 (w = 0.5) from the exact result (38) of [18] which takes the form for arbitrary aspect ratio (1 + w)/(1 − w).
For specular reflection we can compare the E-odd frequency shift data with a calculation based on an extension of the analysis of relaxation times by Swank et al. [15] to frequency shifts, which is described in Sec. V C. In

C. Further common results for any cell geometry
In the adiabatic limit Ω 1 we can compare the second-order frequency shifts ∝ E 2 and ∝ B 2 , for any cell geometry, with the universal expression (18) of [21]. As seen in Figs. 3, 4, 7 and 8 we obtain, within statistics, for all cell geometries and surface reflection laws investigated. For circular resp. rectangular resp. triangular geometry: ς = 1/4 resp. (1 + w 2 )/3 resp. 1/6, in agreement with Eq. (18) of [21]. We also recover the universal adiabatic limit for the E-odd frequency shift. This limit is independent of the particle magnetic moment and can be described [3] as a purely geometric Berry phase. Further, we have verified the direct relationship between the shift ∝ EB and that ∝ E 2 , which is based on the general identity (81) derived in Sec. V A 3. This relationship translates into the identity which is clearly seen in Figs. 2, 3, 7 and 8.

IV. SIMULATIONS
We use the Bloch equation [25] in the form of the two coupled ODE's (44) and (45) of Ref. [3], where is the vertical (horizontal) spin component and Φ is the expectation value of its azimuthal angle measured in the system rotating about the z axis at the Larmor frequency ω 0 . As in [3] we can assume J z = 0 at the start t = 0 of the Ramsey measurement period and, since the perturbation is small, that Φ remains small while J xy remains almost exactly constant. J z oscillates with small amplitude about the equatorial plane. As in Ref. [3] we adopt the field model with constant vertical gradient. Including the motional field, the net magnetic field "seen" along any chord is, for any cell geometry: where time τ = ω 0 t is measured from the epicenter c (Fig. 5) and the reference length (given by the radius R in case of circular geometry) is set equal to 1. Setting S z (τ ) = J z (τ )/J xy and denoting the cumulative time elapsed since the start of the Ramsey period by τ cum we obtain from (53) and (54) the coupled ODE's In (55) we take account of the change −∆Φ of spin angle, as seen in the new reference system after each reflection. ∆Φ equals the sum of all phase jumps ∆Φ i−1,i = α g,b,i−1 + α g,a,i between consecutive chords.
For each chord, we integrate the coupled first-order ODE's (55) numerically, with S z and Φ remaining unchanged at the impact points. To select the E-odd phase shift Φ ζη this process is carried out separately for positive η (corresponding to E pointing up) and for η replaced by −η (for E down). The difference equals 2Φ ζη and the Eodd frequency shift is obtained by division by τ cum . We verified that the asymptotic frequency shifts are independent of initial spin angle Φ 0 and of starting position at any point of the cell area, as expected.
In Figs. 7 and 8 we show also the results for the second order shifts ∝ E 2 and ∝ B 2 .

V. EQUIVALENCE BETWEEN THE SCHRÖDINGER EQUATION AND THE CORRELATION FUNCTION APPROACHES
Note that in this section we mainly use dimensional quantities instead of dimensionless parameters like ζ, η and Ω in order to facilitate reference to some previous work. The symbol τ will represent the independent variable of the correlation functions and the sign convention for the Larmor frequency ω 0 is different.

A. General proof of equivalence
As the novel method of calculation used above differs from the Redfield theory usually applied in this field we now turn to the relation between the two methods.
In the following we use the direct solutions of the Schrödinger equation for confined spin-1/2 particles in arbitrary fields in [19] and [27] to investigate the conditions necessary for the Redfield theory to be valid.
We consider the solution of the Schrödinger equation for arbitrary time varying fields correct to second order, Eqs. (34) and (35)   wave functions g r (t, t 0 ) and h r (t, t 0 ) in the rotating system: As above, Σ(t) = ω x (t) + iω y (t) is the arbitrary time dependent transverse magnetic field.
Eqs. (56) and (57) are reformulations of Eqs. (17) - (19) of [19] with times t, t , t no longer restricted to a single chord but considered a common variable. We have adapted the terminology of [27] to that used in [19] and in the present article. In [27], symbols α r and β r were used for the wave functions at any time, not only for a single chord. We replace them by g r and h r as in (12).

Frequency shifts
The wave function of the spin system, ψ (t) is given by the matrix (12) acting on the initial state. We take that initial state to be the state corresponding to the spin pointing in the x direction and evaluate σ + r = σ x + iσ y r in the state ψ (t) . σ + r represents the rotating component of the spin as seen in the rotating coordinate system. The time dependence of its phase gives the frequency shift and the decay of its magnitude gives T 2. Thus neglecting rapidly varying h 2 r (2ω 0 ) term. We then evaluate the ensemble-averaged frequency shift, using where we have set t 0 = 0. T is the total elapsed time which, in the experiments, is given by the long period of free spin precession in the Ramsey (free induction decay, fid) scheme; thus we are interested in the limit T → ∞.
We rewrite (62) as where the second step is a well known result [28,29].
We have made the usual assumption that the field fluctuations have zero time average and are stationary, i.e. Σ * (t ) Σ(t ) = Σ * (0) Σ(τ = t − t ) = R(τ ) = R(−τ ). This also applies to the case of specular reflections when we take the ensemble average over different trajectories (different α g 's) and different starting positions as has been shown in [11].
In the usual physical situation we have for the observation time T : T ∼ T relax >> τ corr where T relax is the spin relaxation time and R (τ ) → 0 for τ τ corr so that in this limit which is the Redfield result [13].

Relaxation rate 1/T2
We now calculate the relaxation rate 1/T 2 . Since We now specialize to the case of a stationary stochastic system where Σ * (t ) Σ (t ) is a function of (t − t ) only. Consider a square region of the t , t plane between (t 0 , t 0 ), (t 0 , t) , (t, t 0 ) and (t, t). (See ( [28]) for a discussion of this argument). Then the double integral over the top half (t > t ) is seen to be the complex conjugate of the integral over the bottom half (t < t ), so the last term is given by the integral over the entire square.
As a result of this we have (again putting t − t = τ ) where the result is valid for times t > τ corr . Substituting we obtain the usual Redfield result: The second term is absent in the usual treatments as it is normally assumed that the cross correlation between the components of the fluctuating field vanishes but this does not hold in the presence of the E × v field. Thus we have shown that the Redfield expressions for the frequency shift and T 2 relaxation rate, see also [14], follow directly from the second order solution to the Schrödinger equation for arbitrary time dependent field with zero time average.
To calculate the longitudinal relaxation time, T 1 , we calculate ψ (t) starting in the spin up state and and then evaluate σ z r = |g r | 2 − |h r | 2 with the well known result 1/T 1 = 2/T 2 where 1/T 2 is the transverse relaxation rate in the absence of fluctuations in B z , [27].

Expressions for the frequency shifts
Using (70) we find where, as above, we are considering a magnetic field B in the presence of a motional magnetic field (E × v)/c 2 produced by an electric field in the z direction parallel to B (we use b = γE/c 2 and γ = ω 0 /B 0 is the gyromagnetic ratio of the particle) where we have assumed that the motions in the x and y directions are uncorrelated. Compare [18], equations (8) - (10). Specializing to the case of a linear magnetic gradient, ∂B x /∂x = ∂B y /∂y = G, we have Thus according to equations (65) and (75), see also [10], the frequency shift ∼ E 2 is given by while the linear in E shift is given by and similarly,S xx Substituting (78) in (77) we obtain and comparing (80) with (76) we find i. e. the frequency dependence of the two effects differs only by the factor ω 0 . This was noted by Pendlebury et al. [3] for the case of a circular cylinder with specular walls but our result is universal. It holds for all geometries, field configurations and collisional regimes and is clearly shown in Figs. 2, 3, 7 and 8. Using (78) we can also show that the result given in Eq. (64) of [32] does indeed agree with previously published results contrary to the doubt expressed in [32].

B. Correlation functions for a circular cylinder with diffuse wall collisions
We calculate the velocity correlation function for a circular vessel with diffuse reflecting walls as follows.
We consider a trajectory that starts on a chord with α g,0 (Fig. 1) at a distance x from the end of the chord. The particle reaches the end of the chord at time t n=0 . It then travels on a chord α g,1 until a time t 1 . Between t n and t n+1 the chord angle is α g,n+1. At each wall collision at t n the velocity vector rotates through an angle δθ = α g,n+1 + α g,n . For negative ρ n+1 (motion in the positive CCW direction), δθ = π +ρ n +ρ n+1 for both positive and negative ρ n while for positive ρ n+1 (motion in the CW direction) we have δθ = −π +ρ n +ρ n+1 . Since π = −π so long as angles are concerned all that matters is whether we have an odd or even number of π s. This will result in a varying sign when we take the sine or cosine of the angle. We have that between t n and t n+1 the velocity makes an angle θ n with its original direction with The wall collisions occur at times with t 0 = x/v. So the correlation function, after τ = t 0 , consists of a series of pulses between times t n and t n+1 with heights (normalized to v 2 ) cos θ n .
The Fourier transformS vv (ω) (76) will be a sum of terms each corresponding to one of the rectangular pulses. Each term is given by the Fourier transform of a rectangular pulse, S n (ω) = 2 cos θ n sin ω(δt)n where δ n = ωR v cos ρ n and We have used, from (83), We have to average the result over all possible starting positions and directions of motion. That means we average over all possible initial chord angles, α g,0 , and all positions along the chord given by the distance x from the end. x varies between x = 0 and x = 2R cos ρ. Aver-aging (84) over x gives S n (ω) = 2 sin δ 0 δ 0 e −iδ0 cos θ n sin δ n+1 ω e −i(δn+1+2 n k=1 δ k ) .
We write cos θ n = e iθn + e −iθn /2 and consider only the first term. We will then add in the second term at the end.
SummingS n (ω) from (87) over all chords except the first (labeled 00) we have, using (82), For j = 1 ... n the angles ρ j are independent of ρ 0 , ρ n+1 and of each other, so the average of (88) can be written To this we add the contribution of the second term, e −iθn , in the expansion of cos θ n and the contribution of the first chord which adds a term equal to v 2 for a time from t = 0 to t = t 0 = x/v and will have a spectrum Its average over x is and the complete expression for the spectrum is × sin δ n+1 e i(ρn+1−δn+1) 1 1 + e 2i(ρj −δj ) The average has to be taken over all possible starting positions and directions of velocity (x, ρ 0 ) . Since the number of particles located on a given chord is proportional to the length of the chord (l j = 2R cos ρ j ) we use averaging method 1 for the first chord. However because following each wall collision the angles ρ j are distributed by a cosine law (Lambert's law in 2 dimensions) we use method 2 for all other chords. Then using Eq. (76) we can write in agreement with the Eq. (36), taking into account a different sign convention for ω 0 . We have included the normalization factor v 2 in (92) and, on the last term, used (32) and again (76).

C. Correlation functions and frequency shifts for a rectangular cell with specular walls
Substituting expression (75) for Im Σ * (0) Σ(τ ) in equation (65), noting that and integrating by parts we find for the linear in E frequency shift [30]: a result first obtained in [18]. Swank et al. [15] (Eq. (17)), using the form of transport propagator derived by Masoliver et al. [31], have given an analytic result for the Fourier transform, S xx (ω) , of the correlation function x (0) x (τ ) in the case of a rectangular cell with specularly reflecting walls and gas collisions, which together with (94) yields: where l x,y = vτ c /L x,y = λ c /L x,y , L x,y are the lengths of the sides, v is the velocity, τ c is the average time between gas collisions and λ c is the mean free path. This is valid for all values of the mean free path, so taking the ballistic limit l x,y 1 we obtain the result plotted in Fig. 8 with reference to Swank et al. [15]. More details can be found in [26].

VI. SUMMARY AND CONCLUSIONS
Solving the Schrödinger equation directly we have analyzed the frequency shifts in searches for an electric dipole moment of the neutron using ultracold neutrons and comagnetometer atoms confined in measurement cells. The cells are assumed to have a perfectly flat top and bottom and can have rough cylinder walls (Sec. II) and a cross section of generic shape (Sec. III A). For perfectly diffuse 2D roughness scattering on the walls and circular cell cross section we derived expression (35) which is valid for any in-plane particle velocity v. Its generalization to arbitrary cell geometry in the form of Eq. (43) is an approximation.
We have specifically investigated the geometry of a general rectangle in Sec. III B) and compared all results with the simulations described in Sec. IV. So far, general results valid in a broad range of particle velocity had been reported only for circular and rectangular geometry and specular reflection [3,11,15,19].
In Sec. II D we showed that a Larmor field slightly tilted away from the vertical symmetry axis of the cell, due to misalignment inevitable in the experiments, can be taken into account as in Ref. [10]: The shifts remain unchanged but refer to the Larmor frequency ω 0 given by the field magnitude, rather than to ω 0 .
The adiabatic limit (Ω 1) of the E-odd shift agrees with the purely geometric Berry phase while those for the second-order frequency shifts, ∝ E 2 and ∝ B 2 , agree with the general results of [21]. This consistency and, generally, the equivalence between previous approaches and our analysis are elaborated in Sec. V. We show by direct calculation the agreement of the results from the Schrödinger equation with the Redfield theory for the cases of a circular cell with diffuse reflecting walls and of a rectangular cell with specular walls.
These comparisons serve to highlight the conditions of validity of the Redfield theory which makes use of a number of assumptions [12] and has resulted in a considerable literature concerning its range of validity, e. g. [22]. Our results can be applied to cases where the Redfield result no longer holds, such as non-stochastic or partially stochastic systems. An example of the latter is the inevitable situation where the static fields E and B 0 are not exactly aligned or where, for non-spherical cells, B 0 is not exactly aligned with the symmetry axis.
Our method also applies to cases where the observation time is on the order of or shorter than the correlation time, and thus we can describe the transient spin dynamics, i.e. the gradual development of the shift with increasing time subsequent to the start of the free precession in the Ramsey scheme. This information is more detailed than the transverse relaxation time T 2 derived from the Redfield theory or directly from the Schrödinger equation [27] or the Bloch equation [32]. Similarly, studying the vertical spin oscillations as in [19] can provide information not contained in T 1 . Our analysis in the present article does not include gravitational depolarization [33,34] due to the effect of gravity on the UCN spectra in a specific apparatus.
To cover all possible trajectories for the rectangular geometry shown in Fig. 6 it suffices to analyze the four groups of chords listed in Table I. These start from side IV (or equivalently, from II) in any direction (−π/2 < ρ < π/2). The paths of group 1 (4) end on side I (III). We have divided those ending on side II into groups 2 and 3. They are separated by the angle ρ CM (x, w) = arctan [x/(1 − w)] for the path passing through CM. Angle ρ 1 (x, w) = − arctan [(1 + w − x)/(2(1 − w)] is directed towards the upper right corner. For trajectories originating from sides I or III we replace w by −w.
The signs appropriate to the different ranges of ρ and to backward vs. forward motion are shown explicitly. For the arccosine function we use the range from 0 to π.