Extensions of Perturbation Theory for Spin Dressed Neutrons and 3He Comagnetometers in Superfluid 4He

1. Spin Dressing Theory

1.1. Spin Dressing Introduction

Understanding the dynamics of spin dressing in regards to the nEDM experiment is crucial for reliable measurements at never before reached sensitives. From the Golub, Lamoreaux Physics Reports paper [1], the Hamiltonian in Equation 1 contains everything I need for a basic spin dressing framework. Given a spin 1/2 particle in a static magnetic field $B_0$ along +z and an oscillating magnetic field along x with frequency $\omega$, the quantized Hamiltonian is

$$\begin{array}{c}H=\omega a^{*}a+\gamma \sqrt{2\pi \hslash \omega /L^3}{s}_x(a+a^{*})+{\omega }_0{s}_z,\end{array}$$

(1)

where $\gamma$ is the gyromagnetic ratio of the particle, ${\omega }_0=-\gamma B_0$ and $L^3$ is the volume of the trap. I will treat the last term in Equation 1, the static field term, as a perturbation. The Hamiltonian can be transformed [2] such that the unperturbed eigenstates can be defined as

where $|m_x>$ is an eigenstate of $s_x$ with eigenvalue $m_x$ and for a spin 1/2 particle, $m_x=\pm 1/2$. Going forward, I will express $|m_x>$ states at $|+{>}_x$ or $|-{>}_x$. |n> is an eigenstate of $a^{*}a$ with eigenvalue n. $\lambda =\gamma \sqrt{2\pi \hslash \omega /L^3}$ and its definition starts from the energy density of magnetic field $B_1$ given as $B_1^2/2=4\pi n_{av}\hslash \omega /L^3=2n_{av}{\lambda }^2/{\gamma }^2$ such that ${\omega }_1^2=4n_{av}{\lambda }^2$ where ${\omega }_1=\gamma B_1$.

After the transformation, the energy eigenvalues, which are the degenerate in $m_x$, present themselves as

$$\begin{array}{c}{E}_{n,m_x}=n\omega -{\displaystyle \frac{m_x^2{\lambda }^2}{\omega }}.\end{array}$$

(3)

With the Hamiltonian now transformed, expectation values will require calculation of the inner product of two bar states. The calculation and proof originally done in [2] is shown in Appendix A. The Physics Reports [1] has the $m_z$ eigenstates as

where $m=2m_z$ and is an equivalent form of the eigenvectors shown in Appendix A for spin along +z and -z ($m=\pm 1$). Going forward, I will write $m_z$ as m and any spin states with a + or - is a state in the x basis.

1.2. Expectation Values

Given Equation 4 are eigenstates that lift the degeneracy in the Hamiltonian, I can now derive expressions for spin expectation values as done in [1]. For each component that follows, I will first show the matrix elements using Equation 4 and then incorporate the elements into a time dependent wave function.

The matrix elements for ${\sigma }_y$ are

so that the expectation value is

$$\begin{array}{c}<{\sigma }_y\left(t\right)>=\sum _{m^{\prime },q}\sum _{n,n^{\prime }}{\displaystyle \frac{1}{2}}{a}_{n^{\prime }}^{*}{a}_n{e}^{-i{n}^{\prime }\omega t}{e}^{-i{m}^{\prime }{\omega }_{d}t/2}i\\ \left(m{J}_q({\omega }_1/\omega ){(-1)}^q-m^{\prime }{J}_q({\omega }_1/\omega )\right)e^{im{\omega }_{d}t/2}{e}^{in\omega t}.\end{array}$$

(6)

$a_n$ are Glauber coefficients which describe the probability amplitudes of photons of a harmonic field. They can be expressed as

$$\begin{array}{c}{a}_n=e^{-\lambda /2}{\displaystyle \frac{{\lambda }^{n/2}}{{(n!)}^{1/2}}},\end{array}$$

(7)

where $\lambda =<n>$, the average photon number. Then for large n, ${\sum }_n{\left|a_n\right|}^2=1$. Given selection rules on m and $m^{\prime }$ and summing over these values, the expectation value for ${\sigma }_y$ becomes

$$\begin{array}{c}<{\sigma }_y\left(t\right)>=\sum _{m^{\prime },q}{\displaystyle \frac{1}{2}}{e}^{iq\omega t}{e}^{-i/2(m^{\prime }-m){\omega }_{d}t}i\\ \left(m{J}_q({\omega }_1/\omega ){(-1)}^q-m^{\prime }{J}_q({\omega }_1/\omega )\right)\\ =-J_0({\omega }_1/\omega )sin{\omega }_{d}t\\ -\sum _{q>0,even}{J}_q({\omega }_1/\omega )\left[sin({\omega }_d+q\omega )t+sin({\omega }_d-q\omega )t\right].\end{array}$$

(8)

For ${\sigma }_z$, the matrix elements are

So that the expectation value is

$$\begin{array}{c}<{\sigma }_z\left(t\right)>=\sum _{m^{\prime },q}\sum _{n,n^{\prime }}{\displaystyle \frac{1}{2}}{a}_{n^{\prime }}^{*}{a}_n{e}^{-i(n^{\prime }-n)\omega t}{e}^{-i/2(m^{\prime }-m){\omega }_{d}t}\\ \left[m{J}_q({\omega }_1/\omega ){(-1)}^q+m^{\prime }{J}_q({\omega }_1/\omega )\right]\\ =\sum _{q>0,odd}{J}_q({\omega }_1/\omega )\left[cos({\omega }_d-q\omega )t-cos({\omega }_d+q\omega )t\right].\end{array}$$

(10)

And finally, for ${\sigma }_x$

So that the expectation value is

$$\begin{array}{c}<{\sigma }_x\left(t\right)>=\sum _{m,q}\sum _{n,n^{\prime }}{\displaystyle \frac{1}{2}}{a}_{n^{\prime }}^{*}{a}_n{e}^{-i{n}^{\prime }\omega t}{e}^{-i{m}^{\prime }{\omega }_{d}t/2}\\ {\delta }_{n,n^{\prime }}(1-m^{\prime }m)e^{in\omega t}{e}^{i^{\prime }{\omega }_{d}t/2}\\ =\sum _{n,m,m^{\prime }}{\displaystyle \frac{1}{2}}{a}_n^{*}{a}_n{e}^{-i/2(m^{\prime }-m){\omega }_{d}t}(1-m^{\prime }m)=cos{\omega }_{d}t.\end{array}$$

(12)

Thus, I have expressions for how spin expectations values evolve through time.

When it comes to my simulation work, the comparison of simulation data versus these expressions is important to benchmark the simulations before extending them to higher complexities. After making a fourth order Runge-Kutta (RK4) solver in Matlab to solve the spin Bloch equations, I saw there were frequencies matching with theory as well as additional frequencies which have not been accounted for in the Physics Reports. Furthermore, I saw confluence of the RK4 results with the near-exact Hamiltonian method shown in a later section. Given how the RK4 results vary with the parameter $y={\omega }_0/\omega$, perturbation theory became my focus to try to account for the additional frequencies. Also, the experiment uses $y=0.10$ and so examination of whether this parameter is sufficiently small to only use the first order corrections to the eigenvalues is quite important.

2. The 4th Order Runge-Kutta Solver for the Spin Bloch Equations

The spin Bloch equations can be solved in Matlab by constructing a fourth-order Runge-Kutta solver. Given the Bloch equations

$$\begin{array}{c}{\displaystyle \frac{d\overrightarrow{S}}{dt}}=\gamma \left(\overrightarrow{S}\overrightarrow{B}\right),\end{array}$$

(13)

I can specify field parameters, time steps, simulation time and initial conditions as inputs into my RK4. For spin dressing, I simply have to use the appropriate fields and solve for the spin components. The fields will be

$$\begin{array}{c}\overrightarrow{B_x}\left(t\right)=B_{1}cos\omega t\widehat{i}\\ B_y=0\\ \overrightarrow{B_z}=B_0\widehat{k}\\ y\equiv {\omega }_0/\omega ,x\equiv {\gamma }_n{B}_1/\omega ={\omega }_1/\omega .\end{array}$$

(14)

An example of a three dimensional spin dressing simulation between UCN and ³He is shown in Figure 1.

Using initial conditions on the time dependent spin components $S_i\left(t\right)$ of $S_x\left(0\right)=1$, $S_y\left(0\right)=0$ and $S_z\left(0\right)=0$, time steps of $\mathrm{tstep}=0.0005$ seconds for the spin components and $\mathrm{tstep}/2$ for magnetic field time steps, the RK4 will solve for each spin time series.

For the rest of this chapter, the simulation time is $T=10.0s$ (the actual computation can take anywhere from minutes to hours). For the number of particles, I will use $N=1$. Finally, the first results I will show have field parameters of

$$\begin{array}{c}{B}_0=1\mathrm{mG}\\ y=0.10,x=1.189.\end{array}$$

(15)

The results are spins as a function of time: ${\sigma }_x\left(t\right)$, ${\sigma }_y\left(t\right)$ and ${\sigma }_z\left(t\right)$. I then take the fast Fourier transform (FFT) of these time series and to see relevant frequencies. From the theory in the previous section, I expect certain frequencies for each spin component. The plots in Figure 2 to 4 are the FFTs of the spins for the parameters in Equation 15. In light green is simulation data and in black is theory using the expectation values from last section. The time steps are very small and therefore the frequency range is large so the plots are zoomed in for better inspection of the relevant frequencies.

The theory in Figure 2-4 are the predicted results from Equations 8, 10 and 12, respectively. As seen, there are frequencies and harmonics consistent with the expectation values in that section. Furthermore, the ratio of the heights of the expected peaks are in-line with ratios of amplitudes of Bessel functions of the first kind of corresponding order. But what is quite interesting are the peaks that the theory doesn’t predict. In Figure 4 for ${\sigma }_x$, I see these ’unexpected’ frequencies of $\omega$ and ${\omega }_d\pm 2\omega$ are almost negligible. But in Figure 2 showing ${\sigma }_y$, the unexpected peaks at $\omega$ and $3\omega$ are not negligible. And for ${\sigma }_z$ in Figure 3, the unexpected frequencies of $0\omega$ and $2\omega$ have meaningful intensities as well. Figure 5 and 6 show what Figure 3 and 4 showed but for $y=0.01$ meaning $\omega$ is ten times larger.

What’s interesting is that the unexpected frequencies decreased in height substantially between $y=0.10$ and $y=0.01$. Given y is a measure of the perturbing static field relative to the rate of oscillation of the dressing field, intensity dependence on y points to missing higher order perturbation theory terms.

The question then is: are these unexpected frequencies due to the RK4 numerical solver issues in Matlab or is it from higher order perturbation theory terms? The reduction in intensity when y is smaller would suggest the latter. Some theoretical models use y anywhere between 0.10 and 0.01 while the experiment will use $y=0.10$. Figure 7 to 9 show how these unexpected frequency peaks change as the perturbation strength changes. I chose two frequencies from ${\sigma }_x$ ($\omega$ and ${\omega }_d-2\omega$) and one each from ${\sigma }_y$ ($\omega$) and ${\sigma }_z$ ($0\omega$). Then stepped y by 0.0001 between 0.01 and 0.10.

The oscillatory behavior in Figure 7 to 9 most likely comes from numerical issues within the solver. From Figure 9, the unexpected frequencies I chose from ${\sigma }_y$ and ${\sigma }_z$ show linear in y dependence. And from Figure 7 and 8, unexpected frequencies for ${\sigma }_x$ show $y^2$ and y dependence respectively. Going forward, I will call these unexpected frequencies $\mathit{anomalies}$.

Given the dependence on y, the anomalies most likely come from perturbation theory. The next section shows a different method for solving for expectation values and frequencies of spin components to provide further evidence that the anomalies are not artifacts. If it can give confluence to the RK4 method, I will use higher order perturbation theory to more accurately represent this system.

3. The Hamiltonian Method

3.1. Theory

Diagonalizing the Hamiltonian in Equation 1 is a method that for high photon number $N_p$, converges towards an exact solution. However, computation time is a limiting factor. For $N_p$ very large, $n_{av}$ refers to an addition or subtraction away from $N_p$. This also means the $\sqrt{n_{av}}$ in $\lambda$ can be ignored. I will write $n_{av}$ as n going forward. Given this Hamiltonian and how the creation and annihilation operators effect adjacent states, the Hamiltonian matrix in the z basis up to $n=\pm 2$ can be written as

where $a=y/2$ and $b=x/4$. Additionally, I can write the Hamiltonian matrix in the x or nbar basis but for my simulations, I used the z basis.

3.2. Solving for Eigenvalues, Eigenvectors and Time Evolution Operators

Solving for the eigenvalues and eigenvectors will be the large majority of the computation. For calculating the expectation values of time dependent spin components, initialization of the state and spin matrices are needed. Specifically, in the z basis and an initial state ${\psi }_0$ along +x, ${\psi }_0$ will take the form of an Nx1 column vector of all ones and a normalization factor. The spin component operators will be NxN matrices with 2x2 Pauli matrices along their respective diagonals.

With the eigenvalues known, I move to calculating time evolution operators

$$\begin{array}{c}{U}_i=e^{-iw{t}_i\lambda },\end{array}$$

(17)

where $t_i$ is just an element in the array t representing time and since I will use $T=10.0s$ and tstep$=0.0005$s, t is a 1x20001 array. $\lambda$ is the diagonal eigenvalues matrix so that U is an NxN diagonal matrix at a given time.

$\psi \left(t\right)$ is constructed by transforming U into the eigenvector basis, letting it evolve, then transforming back and multiplying it onto ${\psi }_0$

$$\begin{array}{c}{\psi }_i\left(t\right)=K^{+}{U}_i\left(t\right)K·{\psi }_0,\end{array}$$

(18)

where K is the NxN eigenvector matrix. Finally, the expectation values over time take the usual form of

$$\begin{array}{c}{S}_i\left(t\right)=\psi {\left(t\right)}^{+}{S}_i\psi \left(t\right).\end{array}$$

(19)

The results are spin components as time arrays just like in the RK4 and I proceed to take fast Fourier transforms of them to extract frequencies. What I saw was agreement between this method and the RK4. At first, for low N, the Hamiltonian method was all over the place. But as N got larger, the results converged towards what the RK4 showed.

3.3. Comparing the Hamiltonian Method to the RK4

The Hamiltonian method showed reasonable results for $n=\pm 150$ which is $N=602$. However, to be more accurate, I will use $n=\pm 200$ in the following plots. Figure 10 to 13 show ${\sigma }_y$ and ${\sigma }_z$ for $y=0.10$ and $y=0.01$ and as with the RK4 plots, I have zoomed in to show relevant frequencies.

Figure 10. ${\sigma }_y$ frequencies for y=0.10 using the Hamiltonian method, N=802.

Figure 10. ${\sigma }_y$ frequencies for y=0.10 using the Hamiltonian method, N=802.

Figure 11. ${\sigma }_y$ frequencies for y=0.01 using the Hamiltonian method, N=802.

Figure 11. ${\sigma }_y$ frequencies for y=0.01 using the Hamiltonian method, N=802.

Figure 12. ${\sigma }_z$ frequencies for y=0.10 using the Hamiltonian method, N=802.

Figure 12. ${\sigma }_z$ frequencies for y=0.10 using the Hamiltonian method, N=802.

Figure 13. ${\sigma }_z$ frequencies for y=0.01 using the Hamiltonian method, N=802.

Figure 13. ${\sigma }_z$ frequencies for y=0.01 using the Hamiltonian method, N=802.

The immediate conclusion is that the Hamiltonian method is showing very similar results to the RK4 method. The difference in heights between the two methods mostly comes from the difference in time stepping and therefore total intensity. The ratio of heights between the two methods as well as the ratio of heights within a method and compared to the other would converge for n much larger. And again, I am seeing the anomalous frequency heights decreasing as y gets smaller. For comparison to Figure 9, I will loop through four values of y: 0.10, 0.05, 0.02 and 0.01 and show the results for $n=\pm 200$ in Figure 14. The reason I’m only using four values whereas the RK4 plots used 901 y values is because each Hamiltonian method simulation takes a long time. The plot shows a near linear relationship for these anomalies’ heights and y, as seen before.

4. First Order Perturbation Theory, Corrections to the Eigenstates

In the first section, the degeneracy was lifted by finding first order energy shifts. The $m_z$ eigenvectors were then used for calculating expectation values. Given the evidence laid out in the sections showing results from the RK4 and Hamiltonian methods, it’s clear extensions of perturbation theory are needed to account for the anomalous frequencies. More specifically, corrections to the eigenstates are required rather than corrections to the eigenvalues. This is because the only way for the peak heights to have a dependence on $y^n$, where n is the order of perturbation theory, is from corrections to the eigenstates. Whereas corrections to the eigenvalues introduce changes to energy levels and therefore frequencies.

An expansion of an eigenstate using higher order corrections is written as

The first term on the right hand side are the $m_z$ eigenstates whereas the second term is the 1st order correction.

The first order correction is generally written as and so first order corrections to unprimed and primed states are

where I define A as

$$A(m,m^{\prime },q^{\prime })={\displaystyle \frac{\left[m{J}_{q^{\prime }}\left(x\right){(-1)}^{q^{\prime }}+m^{\prime \prime }{J}_{q^{\prime }}\left(x\right)\right]}{4q^{\prime }}}.$$

(24)

For the primed states

where I define $A^{\prime }$ as

$$A^{\prime }(m,m^{\prime },q,q^{\prime })={\displaystyle \frac{\left[m^{\prime }{J}_{q^{\prime }-q}\left(x\right){(-1)}^{q^{\prime }-q}+m^{\prime \prime }{J}_{q^{\prime }-q}\left(x\right)\right]}{4(q^{\prime }-q)}}.$$

(26)

From before, $q=n-n^{\prime }$ but now $q^{\prime }=n-n^{\prime \prime }$ such that $n^{\prime }-n^{\prime \prime }=q^{\prime }-q$. States cannot be equal in the formulation of the first order corrections which then requires $q^{\prime }\ne 0$ and $q^{\prime }-q\ne 0$ for certain values of $m,m^{\prime },m^{\prime \prime }$ but also to avoid singularities.

Now a state $\overline{|n,m>}$ up to first order is

And similarly, a state $\overline{|n^{\prime },m^{\prime }>}$ up to first order is

The matrix elements for ${\sigma }_y$ are

where the first two terms in the brackets are the usual zeroth order, unperturbed terms and $O\left(y^2\right)$ refers to terms proportional to $y^2$ that I will ignore. Then, the expectation value is

$$\begin{array}{c}<{\sigma }_y\left(t\right)>\\ =\sum _{m,m^{\prime },m^{\prime \prime }}\sum _{q,q^{\prime }}{\displaystyle \frac{1}{2}}{a}_{n^{\prime }}^{*}{a}_n{e}^{iq\omega t}{e}^{i(m-m^{\prime }){\omega }_{d}t/2}[im{J}_q\left(x\right){(-1)}^q-i{m}^{\prime }{J}_q\left(x\right)\\ +iy{\displaystyle \frac{\left[m^{\prime }{J}_{q^{\prime }-q}\left(x\right){(-1)}^{q^{\prime }-q}+m^{\prime \prime }{J}_{q^{\prime }-q}\left(x\right)\right]}{4(q^{\prime }-q)}}\left[m{J}_{q^{\prime }}\left(x\right){(-1)}^{q^{\prime }}-m^{\prime \prime }{J}_{q^{\prime }}\left(x\right)\right]\\ +iy{\displaystyle \frac{\left[m{J}_{q^{\prime }}\left(x\right){(-1)}^{q^{\prime }}+m^{\prime \prime }{J}_{q^{\prime }}\left(x\right)\right]}{4q^{\prime }}}\left[-m^{\prime }{J}_{q^{\prime }-q}\left(x\right)+m^{\prime \prime }{J}_{q^{\prime }-q}\left(x\right){(-1)}^{q^{\prime }-q}\right]].\end{array}$$

(30)

Next I will choose $m=m^{\prime }$ so that ${\omega }_d$ no longer contributes to any frequency. Then, I examine the four cases when $m=m^{\prime }=\pm 1$ and $m^{\prime \prime }=\pm 1$. These calculations are located in the appendix in B. I know from the RK4 and diagonalizing the Hamiltonian that for ${\sigma }_y$, the anomalous peaks are for odd q. Choosing odd q, then varying $q^{\prime }$ (which dictates whether $q^{\prime }-q$ is even or odd), I have

$$\begin{array}{c}\sum _{q\ge 0,\mathrm{odd}}{\displaystyle \frac{1}{2}}\left[e^{iq\omega t}+e^{-iq\omega t}\right]iy\\ \left\{\sum _{q^{\prime }\mathrm{odd},q^{\prime }-q\mathrm{even}}\right[{\displaystyle \frac{J_{q^{\prime }}\left(x\right)J_{q^{\prime }-q}\left(x\right)+J_{q^{\prime }}\left(x\right)J_{q^{\prime }-q}\left(x\right)}{q^{\prime }}}\\ +{\displaystyle \frac{-J_{q^{\prime }-q}\left(x\right)J_{q^{\prime }}\left(x\right)-J_{q^{\prime }-q}\left(x\right)J_{q^{\prime }}\left(x\right)}{q^{\prime }-q}}]\\ +\sum _{q^{\prime }\mathrm{even},q^{\prime }-q\mathrm{odd}}[{\displaystyle \frac{-J_{q^{\prime }}\left(x\right)J_{q^{\prime }-q}\left(x\right)-J_{q^{\prime }}\left(x\right)J_{q^{\prime }-q}\left(x\right)}{q^{\prime }}}\\ +{\displaystyle \frac{-J_{q^{\prime }-q}\left(x\right)J_{q^{\prime }}\left(x\right)-J_{q^{\prime }-q}\left(x\right)J_{q^{\prime }}\left(x\right)}{q^{\prime }-q}}\left]\right\}\\ \\ =iy\sum _{q\ge 0,\mathrm{odd}}cosq\omega t\\ \left\{\sum _{q^{\prime }\mathrm{odd},q^{\prime }-q\mathrm{even}}\right[{\displaystyle \frac{J_{q^{\prime }}\left(x\right)J_{q^{\prime }-q}\left(x\right)+J_{q^{\prime }}\left(x\right)J_{q^{\prime }-q}\left(x\right)}{q^{\prime }}}\\ +{\displaystyle \frac{-J_{q^{\prime }-q}\left(x\right)J_{q^{\prime }}\left(x\right)-J_{q^{\prime }-q}\left(x\right)J_{q^{\prime }}\left(x\right)}{q^{\prime }-q}}]\\ +\sum _{q^{\prime }\mathrm{even},q^{\prime }-q\mathrm{odd}}[{\displaystyle \frac{-J_{q^{\prime }}\left(x\right)J_{q^{\prime }-q}\left(x\right)-J_{q^{\prime }}\left(x\right)J_{q^{\prime }-q}\left(x\right)}{q^{\prime }}}\\ +{\displaystyle \frac{-J_{q^{\prime }-q}\left(x\right)J_{q^{\prime }}\left(x\right)-J_{q^{\prime }-q}\left(x\right)J_{q^{\prime }}\left(x\right)}{q^{\prime }-q}}\left]\right\},\end{array}$$

(31)

with $J_{q^{\prime }}\left(x\right)$ and $J_{q^{\prime }-q}\left(x\right)$ commuting. Of the four terms in Equation 31, the first and the fourth will be equal but opposite sign for certain q and $q^{\prime }$ values. Whereas the second and the third will be equal and same sign. I will show two examples of this when $q=1$ and I choose different values of $q^{\prime }$ based off the sum rules in Appendix C.

Finally, the expectation value for $m=m^{\prime }$ and odd q can be written as

$$\begin{array}{c}<{\sigma }_y{\left(t\right)>|}_{m=m^{\prime }}\\ =-4iy\sum _{q>0,\mathrm{odd}}\sum _{q^{\prime }\ne 0,\mathrm{even}}cosq\omega t{\displaystyle \frac{J_{q^{\prime }-q}\left(x\right)J_{q^{\prime }}\left(x\right)}{q^{\prime }}}.\end{array}$$

(32)

The first few terms are

$$\begin{array}{c}-4iy[cos\omega t{\displaystyle \frac{J_1\left(x\right)J_2\left(x\right)}{\left(2\right)}}+cos\omega t{\displaystyle \frac{J_{-3}\left(x\right)J_{-2}\left(x\right)}{(-2)}}\\ +cos3\omega t{\displaystyle \frac{J_{-1}\left(x\right)J_2\left(x\right)}{\left(2\right)}}+cos3\omega t{\displaystyle \frac{J_{-5}\left(x\right)J_{-2}\left(x\right)}{(-2)}}+...]\\ =-2iy{J}_2\left(x\right)[\left[J_1\left(x\right)+J_3\left(x\right)\right]cos\omega t\\ +\left[J_5\left(x\right)-J_1\left(x\right)\right]cos3\omega t+...].\end{array}$$

(33)

For all the other cases ($m\ne m^{\prime }$ with q odd or even and $m=m^{\prime }$ with even q) all terms vanish.

For ${\sigma }_z$, the matrix elements are

So that the expectation value up to y is

$$\begin{array}{c}<{\sigma }_z\left(t\right)>={\displaystyle \sum _{m,m^{\prime },m^{\prime \prime }}}{\displaystyle \sum _{q,q^{\prime }}}{\displaystyle \frac{1}{2}}{a}_{n^{\prime }}^{*}{a}_n{e}^{iq\omega t}{e}^{i(m-m^{\prime }){\omega }_{d}t/2}\\ [m{J}_q\left(x\right){(-1)}^q+m^{\prime }{J}_q\left(x\right)+my{A}^{\prime *}{J}_{q^{\prime }}\left(x\right){(-1)}^{q^{\prime }}+m^{\prime \prime }y{A}^{\prime *}{J}_{q^{\prime }}\left(x\right)\\ +m^{\prime }yA{J}_{q^{\prime }-q}\left(x\right)+m^{\prime \prime }yA{J}_{q^{\prime }-q}\left(x\right){(-1)}^{q^{\prime }-q}].\end{array}$$

(35)

This time, I’ll exclude the zeroth order terms. The calculation for the sum over spins for $m=m^{\prime }$ is similar to the ${\sigma }_y$ calculation. Then the linear in y terms are

$$\begin{array}{c}<{\sigma }_z\left(t\right)>\\ =\sum _{q,q^{\prime }}{\displaystyle \frac{1}{2}}{e}^{iq\omega t}y[{\displaystyle \frac{J_{q^{\prime }-q}\left(x\right){(-1)}^{q^{\prime }-q}{J}_{q^{\prime }}\left(x\right){(-1)}^{q^{\prime }}+J_{q^{\prime }-q}\left(x\right)J_{q^{\prime }}\left(x\right)}{q^{\prime }-q}}\\ +{\displaystyle \frac{J_{q^{\prime }}\left(x\right){(-1)}^{q^{\prime }}{J}_{q^{\prime }-q}\left(x\right)+J_{q^{\prime }}\left(x\right)J_{q^{\prime }-q}\left(x\right){(-1)}^{q^{\prime }-q}}{q^{\prime }}}].\end{array}$$

(36)

For ${\sigma }_z\left(t\right)$, I saw the anomalies are for even q. Then I have

$$\begin{array}{c}<{\sigma }_z\left(t\right)>=\sum _{q\ge 0,\mathrm{even}}{\displaystyle \frac{1}{2}}\left[e^{iq\omega t}+e^{-iq\omega t}\right]y\\ \left\{\sum _{q^{\prime }\mathrm{odd},q^{\prime }-q\mathrm{odd}}\right[{\displaystyle \frac{J_{q^{\prime }-q}\left(x\right)J_{q^{\prime }}\left(x\right)+J_{q^{\prime }-q}\left(x\right)J_{q^{\prime }}\left(x\right)}{q^{\prime }-q}}\\ +{\displaystyle \frac{-J_{q^{\prime }}\left(x\right)J_{q^{\prime }-q}\left(x\right)-J_{q^{\prime }}\left(x\right)J_{q^{\prime }-q}\left(x\right)}{q^{\prime }}}]\\ +\sum _{q^{\prime }\mathrm{even},q^{\prime }-q\mathrm{even}}[{\displaystyle \frac{J_{q^{\prime }-q}\left(x\right)J_{q^{\prime }}\left(x\right)+J_{q^{\prime }-q}\left(x\right)J_{q^{\prime }}\left(x\right)}{q^{\prime }-q}}\\ +{\displaystyle \frac{J_{q^{\prime }}\left(x\right)J_{q^{\prime }-q}\left(x\right)+J_{q^{\prime }}\left(x\right)J_{q^{\prime }-q}\left(x\right)}{q^{\prime }}}\left]\right\}\\ =4y\sum _{q^{\prime }\ne 0,\mathrm{odd}}{\displaystyle \frac{J_{q^{\prime }}^2\left(x\right)}{q^{\prime }}}\\ +4y\sum _{q>0,\mathrm{even}}\sum _{q^{\prime }\ne 0,\mathrm{odd}}cosq\omega t{\displaystyle \frac{J_{q^{\prime }-q}\left(x\right)J_{q^{\prime }}\left(x\right)}{q^{\prime }}},\end{array}$$

(37)

with the first 4 terms being

$$\begin{array}{c}4y[{\displaystyle \frac{J_1\left(x\right)J_1\left(x\right)}{\left(2\right)}}+{\displaystyle \frac{J_{-1}\left(x\right)J_{-1}\left(x\right)}{(-1)}}-cos2\omega t{\displaystyle \frac{J_1\left(x\right)J_1\left(x\right)}{\left(1\right)}}\\ -cos2\omega t{\displaystyle \frac{J_3\left(x\right)J_1\left(x\right)}{(-1)}}]\\ =4y\left[\right(\mathrm{oscillatory}q=0\mathrm{term})\\ -cos2\omega t{J}_1\left(x\right)\left[J_1\left(x\right)-J_3\left(x\right)+...\right]+...].\end{array}$$

(38)

Regarding ${\sigma }_x$, I saw anomalous frequencies at $q\omega$ and ${\omega }_d\pm 2\omega$. The matrix elements for ${\sigma }_x$ are

so that the expectation value is

$$\begin{array}{c}<{\sigma }_x\left(t\right)>=\sum _{m,m^{\prime },m^{\prime \prime }}\sum _{n,n^{\prime },n^{\prime \prime }}{\displaystyle \frac{1}{2}}{a}_{n^{\prime }}^{*}{a}_n{e}^{iq\omega t}{e}^{i(m-m^{\prime }){\omega }_{d}t/2}\\ [(1-m{m}^{\prime }){\delta }_{n,n^{\prime }}\\ +y{A}^{\prime }(1-m{m}^{\prime \prime }){\delta }_{n,n^{\prime \prime }}+yA(1-m^{\prime }{m}^{\prime \prime }){\delta }_{n^{\prime },n^{\prime \prime }}].\end{array}$$

(40)

For summing over all n’s and m’s in Equation 40, I have the following cases and their constraints (due to $k\ne n$ in the formulation for first order perturbation theory)

$$\begin{array}{c}n=n^{\prime },n^{\prime \prime }\ne n\to n^{\prime \prime }\ne n^{\prime }\to \mathrm{no}m\mathrm{constraints};\\ \\ n=n^{\prime },n^{\prime \prime }=n\to n^{\prime \prime }=n^{\prime }\to m^{\prime \prime }\ne m^{\prime },m^{\prime \prime }\ne m;\\ \\ n\ne n^{\prime },n^{\prime \prime }\ne n\to n^{\prime \prime }=n^{\prime }\to m^{\prime \prime }\ne m^{\prime };\\ \\ n\ne n^{\prime },n^{\prime \prime }=n\to n^{\prime \prime }\ne n^{\prime }\to m^{\prime \prime }\ne m.\end{array}$$

(41)

Of these four cases, the first leads to the zeroth order $cos{\omega }_{d}t$ term. The second case has $q^{\prime }=0$ and $q^{\prime }-q=0$ which is not allowed. Cases three and four are the ones I care about. Given I’m in a spin 1/2 system, a constraint like $m^{\prime \prime }\ne m^{\prime }$ means $m^{\prime \prime }=-m$ and so I can simply replace all $m^{\prime \prime }$’s by $-m$. Also, the third case has $q^{\prime }-q=0$ so that any $q^{\prime }$ can be replaced by q. And the fourth case has $q^{\prime }=0$ but given $n^{\prime \prime }\ne n^{\prime }$, there’s no singularity. The calculation for the third case is

$$\begin{array}{c}\sum _{m,m^{\prime },m^{\prime \prime }}\sum _{n,n^{\prime }}{\displaystyle \frac{1}{2}}{a}_{n^{\prime }}^{*}{a}_n{e}^{iq\omega t}{e}^{i(m-m^{\prime }){\omega }_{d}t/2}y\\ {\displaystyle \frac{m{J}_q\left(x\right){(-1)}^q+m^{\prime \prime }{J}_q\left(x\right)}{4\left(q\right)}}[1-m^{\prime }{m}^{\prime \prime }]\\ \\ =\sum _q\sum _{m,m^{\prime }}{\displaystyle \frac{1}{2}}{e}^{iq\omega t}{e}^{i(m-m^{\prime }){\omega }_{d}t/2}y\\ {\displaystyle \frac{m{J}_q\left(x\right){(-1)}^q-m^{\prime }{J}_q\left(x\right)}{4q}}[1+{\left(m^{\prime }\right)}^2]\\ \\ =y\sum _q{e}^{iq\omega t}{\displaystyle \frac{J_q\left(x\right)}{4q}}[\left[{(-1)}^q-1\right]+e^{-i{\omega }_{d}t}\left[-{(-1)}^q-1\right]\\ +e^{i{\omega }_{d}t}\left[{(-1)}^q+1\right]+\left[-{(-1)}^q+1\right]]\\ =y\sum _q{e}^{iq\omega t}{\displaystyle \frac{J_q\left(x\right)}{4q}}\left[{(-1)}^q+1\right]\left[e^{i{\omega }_{d}t}-e^{-i{\omega }_{d}t}\right]\\ =y\sum _{q>0,\mathrm{even}}\left[e^{iq\omega t}+e^{-iq\omega t}\right]{\displaystyle \frac{J_q\left(x\right)}{4q}}\left(2\right)(2isin{\omega }_{d}t)\\ =2iy\sum _{q>0,\mathrm{even}}{\displaystyle \frac{J_q\left(x\right)}{q}}cosq\omega tsin{\omega }_{d}t.\end{array}$$

(42)

And the fourth case

$$\begin{array}{c}\sum _{m,m^{\prime },m^{\prime \prime }}\sum _{n,n^{\prime }}{\displaystyle \frac{1}{2}}{a}_{n^{\prime }}^{*}{a}_n{e}^{iq\omega t}{e}^{i(m-m^{\prime }){\omega }_{d}t/2}y\\ {\displaystyle \frac{m^{\prime }{J}_{q^{\prime }-q}\left(x\right){(-1)}^{q^{\prime }-q}+m^{\prime \prime }{J}_{q^{\prime }-q}\left(x\right)}{4(q^{\prime }-q)}}[1-m{m}^{\prime \prime }]\\ =\sum _q\sum _{m,m^{\prime }}{\displaystyle \frac{1}{2}}{e}^{iq\omega t}{e}^{i(m-m^{\prime }){\omega }_{d}t/2}y\\ {\displaystyle \frac{m^{\prime }{J}_{-q}\left(x\right){(-1)}^{-q}-m{J}_{-q}\left(x\right)}{4(-q)}}[1+m^2]\\ =-y\sum _q\sum _{m,m^{\prime }}{e}^{iq\omega t}{e}^{i(m-m^{\prime }){\omega }_{d}t/2}{\displaystyle \frac{m^{\prime }{J}_q\left(x\right)-m{J}_q\left(x\right){(-1)}^q}{2q}}\\ =-{\displaystyle \frac{1}{2}}y\sum _q{e}^{iq\omega t}{\displaystyle \frac{J_q\left(x\right)}{2q}}[\left[1-{(-1)}^q\right]+e^{-i{\omega }_{d}t}\left[1+{(-1)}^q\right]\\ +e^{i{\omega }_{d}t}\left[-1-{(-1)}^q\right]+\left[-1+{(-1)}^q\right]]\\ =-y\sum _q{e}^{iq\omega t}{\displaystyle \frac{J_q\left(x\right)}{4q}}\left[1+{(-1)}^q\right]\left[e^{-i{\omega }_{d}t}-e^{i{\omega }_{d}t}\right]\\ =-y\sum _{q>0,\mathrm{even}}\left[e^{iq\omega t}+e^{-iq\omega t}\right]{\displaystyle \frac{J_q\left(x\right)}{4q}}\left(2\right)(-2isin{\omega }_{d}t)\\ =2iy\sum _{q>0,\mathrm{even}}{\displaystyle \frac{J_q\left(x\right)}{q}}cosq\omega tsin{\omega }_{d}t.\end{array}$$

(43)

Adding together the results from the third and fourth cases gives the full linear in y dependence of ${\sigma }_x$ as

$$4iy\sum _{q>0,\mathrm{even}}{\displaystyle \frac{J_q\left(x\right)}{q}}cosq\omega tsin{\omega }_{d}t.$$

(44)

From the results of Equation 44, it can be seen that the ${\sigma }_x$${\omega }_d\pm q\omega$ for even q, anomaly has be resolved. Also just as important is that the $q\omega$ anomaly did not show itself since data in the earlier section suggested it would depend on $y^2$. The $y^2$ terms from the first order corrections that I ignored plus the $y^2$ terms that would come about in second order corrections should show the frequencies $q\omega$ for odd q. The plots with the first order corrections added to the theory, analogous to Figure 2 to 4 are shown below in Figure 15 to 17.

In Figure 15, the ${\omega }_d\pm 2\omega$ theoretical heights in black from the first order corrections are too large. In Figure 16, the $\omega$ theory peak is too small while the $3\omega$ is too large. Figure 17 has discrepancies between RK4 heights and the new terms added to the theory as well. I also simply used 0.10 as a placeholder for the $q=0$ oscillatory term. The differences are most likely due to either factors of 2 missing for certain frequency amplitudes, missing Bessel function terms since I stopped the sums after two terms or issues with normalization for Fourier transforms in Matlab.

What has been demonstrated is that the $J_0$ approximation is not valid for $y=0.10$ considering the non-negligible intensities in the ${\sigma }_y$ and ${\sigma }_z$ FFTs. The next section will show that again, when using $y=0.10$, perturbation theory needs to be extended for the eigenvalues as well.

5. The Critical Dressing Condition

Critical dressing is defined by when the precession rate about the z-axis is matched between the two particle species. In regards of UCN and ³He, this means matching their effective Larmor frequencies. From earlier, the normal Larmor frequencies ${\omega }_0=-\gamma B_0$ has been shifted by a factor of $J_0\left(x\right)$. Then the critical dressing condition is:

$$\begin{array}{c}-{\gamma }_n{B}_0[J_0\left(x_c\right)-\alpha J_0\left(\alpha x_c\right)]=0,\end{array}$$

(45)

where $x_c={\omega }_1/\omega$ is the critical dressing parameter and has been chosen to satisfy Equation 45 by adjusting the magnitude of the oscillating field $B_1$. Thus, the condition is defined using the first order corrections to the eigenvalues. When I mentioned the $J_0$ approximation last section, this is precisely what the approximation says: the first order corrections to the eigenvalues are sufficient. However, I found in my simulations that the critical dressing parameter was varying drastically from the solution to Equation 45. Not only that but $x_c$ also had $y^2$ dependence. Furthermore, the reason I investigated further into the critical dressing condition is because I wanted to perform simulations where a ³He pseudomagnetic field would be modulated via the dressing field and its affect nullified. If the critical dressing parameter changes, the point where the modulation centers on changes. A plot of $J_0\left(x\right)$ is shown in Figure 18 and the absolute value of the critical dressing condition is shown in Figure 19.

The first zero is at the critical dressing point the experiment will use because the nEDM term gets diluted by the Bessel function correction so the largest value of $J_0\left(x\right)$ should be used. This solution in the $J_0$ approximation is $x_c=1.189018$. When I performed simulations where I solved for $x_c$ using the same field parameters as the experiment ($B_0=30$ mG, $y=0.10$), critical dressing was at $x_c=1.18246...$. This 0.55% difference between my simulation value and the analytical $J_0$ approximation value was the first indication that something was missing. The second indication was when I tried modulating the dressing field. When the dressing field is set for critical dressing, the neutrons and Helium precess about the z-axis at the same rate. Since dressing is controlling the relative rate, a change in x has the effect, when averaged over, of producing a $\delta \omega$ about the z-axis between the UCN and ³He. The Bessel function is asymmetric about critical dressing so I will label the $\delta \omega$’s as ${\omega }_{z1}$ and ${\omega }_{z2}$

$$\begin{array}{c}{B}_{rf1}=(x_c+x_m)\omega /{\gamma }_n\\ B_{rf2}=(x_c-x_m)\omega /{\gamma }_n\\ {\omega }_{z1}=-\gamma B_0[J_0(x_c+x_m)-\alpha J_0\left(\alpha (x_c+x_m)\right)]\\ {\omega }_{z2}=-\gamma B_0[J_0(x_c-x_m)-\alpha J_0\left(\alpha (x_c-x_m)\right)]\end{array}$$

(46)

where I will write $B_1$ as $B_{rf}$ going forward since there will be two different magnitudes of the dressing field. Intuitively, ${\omega }_{z1}$ being near the magnitude of ${\omega }_{z2}$ should give symmetry about critical dressing. This would make sense since the UCN and ³He need to alternate relative phase shifts in a symmetric manner while $J_0$ itself is asymmetric as mentioned. The data I looked at for symmetry was the scintillation signal which the experiment relies on. This signal is maximum when there’s a maximum component of anti-parallel spins between the two species (Appendix D talks about the spin dependent scattering length of of neutron-³He capture which leads to a periodic scintillation rate and a pseudomagnetic field). Therefore, observing the behavior of ${\overrightarrow{\sigma }}_3·{\overrightarrow{\sigma }}_n$ with modulation gives evidence as to whether the modulation is doing a good enough job to offset the pseudomagnetic field. But what I first saw was that ${\omega }_{z1}\approx \left|{\omega }_{z2}\right|$ was leading to large asymmetries in ${\overrightarrow{\sigma }}_3·{\overrightarrow{\sigma }}_n$ as shown in Figure 20.

Note, in most of my simulations, the initial conditions have the spin parallel to each other. However Figure 20 has the spins initially perpendicular so that the symmetry or asymmetry about 0 can be seen easier. As can be seen, the blue curve which had ${\omega }_{z1}=-{\omega }_{z2}$, was asymmetric about 0 and therefore continuously had phase shifts building. Whereas the three other curves for $x_{m1}=x_{m2}$ are mostly symmetric about 0. What this points to given the definition of ${\omega }_z$ is that the expressions for ${\omega }_{z1,2}$ are wrong. This can be conceptualized by thinking about why symmetric $x_m$ is leading to ${\overrightarrow{\sigma }}_3·{\overrightarrow{\sigma }}_n$ symmetry. Looking at Figure 19, symmetric $x_m$ ($x_{m1}\approx x_{m2}$) for $x_m$ small should lead to ${\omega }_{z1}\approx \left|{\omega }_{z2}\right|$. But as can be seen in the legend of Figure 20, ${\omega }_1$ and ${\omega }_2$ have large percentage differences. This led me to think back to whether the $J_0$ approximation is valid for $y=0.10$ and the idea of extending the eigenvalues to higher orders.

6. Corrections to the Eigenvalues

The general expression for the second order correction to the eigenvalues is

and using the bar eigenstates of $m_z$ where $m_z=\pm {\displaystyle \frac{1}{2}}=\pm {\displaystyle \frac{m}{2}}$

Next are the third order corrections to the eigenvalues which I suspected to be the source. The third order corrections to the eigenvalues have two summation terms and one is easier to calculate than the other. The general expression is

Defining $V_{nm,n^{\prime }{m}^{\prime }}=\overline{<n^{\prime },m^{\prime }|}{\displaystyle \frac{1}{2}}{\omega }_0{\sigma }_z\overline{|n,m>}$, the third order correction in this context is

The Physics Reports [1] made the approximation that the first term, which is more complicated than the second, is equal in magnitude as the second term and therefore did not calculate it. What I found is that the first term is actually larger than the second term and when included in the critical dressing condition, produces much better convergence with simulation values. The calculation for the second term in Equation 50 is

such that for different q’s, $m^{\prime }=m$ or $m^{\prime }\ne m$. The terms up to $J_3\left(x\right)$ are

$$\begin{array}{c}\sum _{m^{\prime }=\pm m}[{\displaystyle \frac{{[-m{J}_1\left(x\right)+m^{\prime }{J}_1\left(x\right)]}^2}{{\left(1\right)}^2}}+{\displaystyle \frac{{[m{J}_1\left(x\right)-m^{\prime }{J}_1\left(x\right)]}^2}{{(-1)}^2}}\\ +{\displaystyle \frac{{[m{J}_2\left(x\right)+m^{\prime }{J}_2\left(x\right)]}^2}{{\left(2\right)}^2}}+{\displaystyle \frac{{[m{J}_2\left(x\right)+m^{\prime }{J}_2\left(x\right)]}^2}{{(-2)}^2}}\\ +{\displaystyle \frac{{[-m{J}_3\left(x\right)+m^{\prime }{J}_3\left(x\right)]}^2}{{\left(3\right)}^2}}+{\displaystyle \frac{{[m{J}_3\left(x\right)-m^{\prime }{J}_3\left(x\right)]}^2}{{(-3)}^2}}+...],\end{array}$$

(52)

so that the result is

$$\begin{array}{c}=-{\displaystyle \frac{1}{4}}m{\omega }_0{J}_0\left(x\right)y^2\left[J_1^2\left(x\right)+{\displaystyle \frac{J_2^2\left(x\right)}{4}}+{\displaystyle \frac{J_3^2\left(x\right)}{9}}+...\right]\\ =-{\displaystyle \frac{1}{2}}{m}_z{\omega }_0{J}_0\left(x\right)y^2\left[J_1^2\left(x\right)+{\displaystyle \frac{J_2^2\left(x\right)}{4}}+{\displaystyle \frac{J_3^2\left(x\right)}{9}}+...\right].\end{array}$$

(53)

For the calculation of the first term in Equation 50, first I would like to consider the terms that contain $J_0\left(x\right)$ since they will be similar in magnitude to the result in Equation 53. The only way to get a $J_0\left(x\right)$ in the first term is for $m=k$ which means $n^{\prime }=n^{\prime \prime }$ (but with $m^{\prime }=\pm m^{\prime \prime }$) giving

where the sum over $m^{\prime \prime }$ made $m^{\prime \prime }=m^{\prime }$ or else the whole term vanishes. Continuing the calculation of just the sum

$$\begin{array}{c}\sum _{m^{\prime },q\ne 0}{\displaystyle \frac{m^{\prime }{[m{J}_q\left(x\right){(-1)}^q+m^{\prime }{J}_q\left(x\right)]}^2}{q^2}}\\ =\sum _{m^{\prime }}[{\displaystyle \frac{m^{\prime }{[-m{J}_1\left(x\right)+m^{\prime }{J}_1\left(x\right)]}^2}{1^2}}+{\displaystyle \frac{m^{\prime }{[m{J}_1\left(x\right)-m^{\prime }{J}_1\left(x\right)]}^2}{{(-1)}^2}}\\ +{\displaystyle \frac{m^{\prime }{[m{J}_2\left(x\right)+m^{\prime }{J}_2\left(x\right)]}^2}{2^2}}+{\displaystyle \frac{m^{\prime }{[m{J}_2\left(x\right)+m^{\prime }{J}_2\left(x\right)]}^2}{{(-2)}^2}}\\ +{\displaystyle \frac{m^{\prime }{[-m{J}_3\left(x\right)+m^{\prime }{J}_3\left(x\right)]}^2}{3^2}}+{\displaystyle \frac{m^{\prime }{[m{J}_3\left(x\right)-m^{\prime }{J}_3\left(x\right)]}^2}{{(-3)}^2}}+...]\\ =\left[{\displaystyle \frac{-16m^3{J}_1^2\left(x\right)}{1}}+{\displaystyle \frac{16m^3{J}_2^2\left(x\right)}{4}}+{\displaystyle \frac{-16m^3{J}_3^2\left(x\right)}{9}}+...\right].\end{array}$$

(55)

The summation over $m^{\prime }$ is simply a summation of $m^{\prime }=m$ and $m^{\prime }=-m$. Then substituting $m=2m_z$ as before, the final result for the first term in Equation 50 is

$$\begin{array}{c}{\displaystyle \frac{1}{32}}{\omega }_0{J}_0\left(x\right)y^2\sum _{m^{\prime },m^{\prime \prime },q\ne 0}{\displaystyle \frac{{[m{J}_q\left(x\right){(-1)}^q+m^{\prime }{J}_q\left(x\right)]}^2[m^{\prime }+m^{\prime \prime }]}{q^2}}\\ =m_z{\omega }_0{J}_0\left(x\right)y^2\left[-J_1^2\left(x\right)+{\displaystyle \frac{J_2^2\left(x\right)}{4}}-{\displaystyle \frac{J_3^2\left(x\right)}{9}}+...\right].\end{array}$$

(56)

The new critical dressing condition is

$$\begin{array}{c}({\omega }_n-{\omega }_3)=0=(-\gamma B_0)[J_0\left(x_c\right)-\alpha J_0\left(\alpha x_c\right)\\ -J_0\left(x_c\right){\displaystyle \frac{y^2}{2}}\left[J_1^2\left(x_c\right)+{\displaystyle \frac{J_2^2\left(x_c\right)}{4}}+{\displaystyle \frac{J_3^2\left(x_c\right)}{9}}+...\right]\\ +{\alpha }^3{J}_0\left(\alpha x_c\right){\displaystyle \frac{y^2}{2}}\left[J_1^2\left(\alpha x_c\right)+{\displaystyle \frac{J_2^2\left(\alpha x_c\right)}{4}}+{\displaystyle \frac{J_3^2\left(\alpha x_c\right)}{9}}+...\right]\\ \\ +J_0\left(x\right)y^2\left[-J_1^2\left(x_c\right)+{\displaystyle \frac{J_2^2\left(x_c\right)}{4}}-{\displaystyle \frac{J_3^2\left(x_c\right)}{9}}...\right]\\ -{\alpha }^3{J}_0\left(x_c\right)y^2\left[-J_1^2\left(\alpha x_c\right)+{\displaystyle \frac{J_2^2\left(\alpha x_c\right)}{4}}-{\displaystyle \frac{J_3^2\left(\alpha x_c\right)}{9}}+...\right]]\\ \\ =(-\gamma B_0)[J_0\left(x_c\right)-\alpha J_0\left(\alpha x_c\right)\\ -J_0\left(x_c\right){\displaystyle \frac{y^2}{2}}\left[3J_1^2\left(x_c\right)-{\displaystyle \frac{J_2^2\left(x_c\right)}{4}}+{\displaystyle \frac{3J_3^2\left(x_c\right)}{9}}+...\right]\\ +{\alpha }^3{J}_0\left(\alpha x_c\right){\displaystyle \frac{y^2}{2}}\left[3J_1^2\left(\alpha x_c\right)-{\displaystyle \frac{J_2^2\left(\alpha x_c\right)}{4}}+{\displaystyle \frac{3J_3^2\left(\alpha x_c\right)}{9}}+...\right]].\end{array}$$

(57)

Now, going back to the other terms that do not include $J_0\left(x\right)$, they will either vanish when summed over or be quite small. An example is shown in Appendix E.

Figure 21 shows the zoomed-in view of the new critical dressing condition and how $x_c$ varies for different values of y. This demonstrates how the nEDM experiments needs to account for this shift in $x_c$ not only for the sake of critical dressing and modulating the dressing field but for an accurate expression for the scintillation rate. This third order correction also dilutes the nEDM signal.

What I now see from my simulations is much better agreement between my simulation $x_c$ and the analytical expression for $x_c$ that includes the third order corrections. Figure 22 shows the analytical and simulation critical dressing parameters for different values of y. And Figure 23 shows the percentage difference of these values at each given value of y. There is still a $y^2$ dependence in Figure 23 because the third order correction terms I’ve ignored are of course still proportional to $y^2$. Given I now have an accurate expression for the critical dressing condition, I can more accurately modulate around it.

7. Modulation and the Pseudomagnetic Field

With the third order corrections from the last section, I will briefly show what modulation looks like with a pseudomagnetic field on and off. Figure 24 shows a near symmetric scintillation signals with modulation on and an applied pseudomagnetic field off.

As seen in the legend, the corrected ${\omega }_{z1,2}$ do have magnitudes much closer to each other than before.

For the magnitude of $\overrightarrow{P}$, I wanted to use the value from the paper by Leung et al.[3] where for $x_3={10}^{-10}$ and $P_3=1$ such that $|\overrightarrow{P}|\approx 0.10\mu$G. However, this magnitude is very small compared to the fields already present and the sensitivity of the simulation so instead I looked at $|\overrightarrow{P}|\approx 10\mu G$. The relevant plots are shown in Figure 25 to 27. For a benchmark, I first show the difference in $\overrightarrow{{\sigma }_3}·\overrightarrow{{\sigma }_n}$ for the pseudomagnetic field off and on without modulation. Followed by plots where modulation will cancel out the effects from $\overrightarrow{P}$. Note that in Figure 25, I chose to use the color blue for the pseudomagnetic field on and orange for it off which differs from other plots. This was solely to give a better contrast to the naked eye.

Figure 25 showed that $\overrightarrow{P}$ was inducing a non-negligible phase shift between the ³He and UCN that grows linear in time (with no modulation). Then Figure 26 and 27 show the modulation doing a good job at canceling the phase difference that was built in the plots without modulation.

Finally, one way which Robert Golub proposed to quantitatively check how the modulations were performing was to check the relationship between the first and second harmonic signals in Fourier space. The idea is that the first harmonic at 1.25 Hz should be near zero if the modulations are symmetric about $x_c$. Whereas the second harmonic will have a much larger height and increase in height as the first harmonic signal decreases. So checking the ratio of the second to the first harmonic with modulation on pseudomagnetic field off and modulation on pseudomagnetic field on, will be a better test of whether the modulations are doing their job. I optimized symmetry by looping through $x_{m1}$ values that produced the lowest first harmonic signal. The results were $x_{m1}=0.05546$ and $x_{m2}=0.056107$, close to the values I’ve been using and I used a simulation time of 8 seconds. Figure 28 shows the FFT of $\overrightarrow{{\sigma }_3}·\overrightarrow{{\sigma }_n}$ with modulations, no pseudomagnetic field and zoomed on the modulation harmonics.

The FFTs are shows in Figure 29 where I zoomed-in on the harmonics but further zoomed-in on the first (1.25 Hz) and second (2.5 Hz) harmonics. The colors have been changed to better see where the curves differ and overlap.

The black curve for when the pseudomagnetic field is on has a larger value at the first harmonic which is expected. Due to the binning of the frequency axis which comes from the time steps and length of the simulation, there was no value exactly on 1.25 Hz or 2.5 Hz. For the 1.25 Hz, I took the average of the values at 1.1875 Hz and 1.3125 Hz. For the 2.5 Hz, I used the value at 2.4375 Hz. The results were that the ratio of the second harmonic signal to the first with the pseudomagnetic field off was 30.055. The ratio with the pseudomagnetic field on was 13.346. Thus, the pseudomagnetic field being flipped to on changed the relative ratios by about a factor of 2.25. This change is a good benchmark as to how well the modulations are doing going forward and with a field that is not 100x the experiment’s strength.