The spin Bloch equations can be solved in Matlab by constructing a fourth-order Runge-Kutta solver. Given the Bloch equations 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 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 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.

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 For the primed states

where I define $A^{\prime }$ as 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 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 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 The first few terms are 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 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 For ${\sigma }_z\left(t\right)$, I saw the anomalies are for even q. Then I have with the first 4 terms being 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 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) 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 And the fourth case Adding together the results from the third and fourth cases gives the full linear in y dependence of ${\sigma }_x$ as 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.

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: 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}$ 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.

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 so that the result is 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 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 The new critical dressing condition is 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.