Discrepancy in clock-qubit fields

[RHanley1]

I have noticed that the magic-fields for the clock-qubits don't quite match the numbers in TPH thesis.

I have calculated the fields using the ion_phys package example "clock_qubits.py", as well explicitly using the Breit-Rabi formula. Both the ion_phys package and the BR formula agree when using the same optimisation procedure. However, if one plots the transition frequency and manually looks at the minimum field, it does agree with TPH.

We currently calculate the gradient by adding a small positive change to B (see below), but this skews the result as it is not symmetric about the original B.

def df_dB(ion, lower, upper, eps=1e-6):
    f = ion.delta(lower, upper)
    ion = deepcopy(ion)
    ion.setB(ion.B+eps)
    fpr = ion.delta(lower, upper)
    return (fpr - f)/eps

We should calculate the gradient symmetrically about a point. I have checked both methods for calculating the magic fields using the BR formula and ion_phys, and if you do it symmetrically you get agreement with TPH thesis. Not sure what the best way to implement this change, but something like

def df_dB(ion, lower, upper, eps=1e-6):
    ion = deepcopy(ion)
    ion.setB(ion.B-eps)
    fl = ion.delta(lower, upper)
    ion.setB(ion.B+2*eps)
    fu = ion.delta(lower, upper)
    return (fu - fl)/(2*eps)

maybe?

[pathfinder49]

I was suprised, when I saw the ionphys was implementing numerical differentiation. I'm always inclined to use a package. Numpy has somehting like this.

[hartytp]

The numbers in my thesis are based on the analytic derivative of the BR formula. If we want something accurate then that is what we should do in this package as well. I think there is a TODO comment in the code to that effect.

[hartytp]

IMO this code should be explicitly special cased (and documented) to only handle J=1/2 systems and use the BR formula rather than diagonalising H. As an added bonus, that will make it much faster.

Might be nice to add a test that compares the BR to the diagonalisation for the energy levels.

[RHanley1]

Yes you could do that, but if one wanted to find out df/dB or any other derivative for any other transitions, you'd get the wrong answer. You could add in the special case of analytic BR, but the function we currently have shouldn't numerically give the wrong answer.

[hartytp]

Yes you could do that, but if one wanted to find out df/dB or any other derivative for any other transitions, you'd get the wrong answer

Is that something that comes up? I've never wanted to know that for a transition outside the ground level. (If I only want a rough idea I'd tend to reach for pen and paper rather than bothering to open python and write code).

But, yes, if we use the BR then we either need to (a) raise an exception if J isn't 1/2 or (b) revert to a numerical method. Whichever we do obviously needs docs and tests.

The concern I have with relying on numerics is accuracy. e.g. for g factor measurements you really want to be confident that there are no rounding errors anywhere in the calculation. Going through the whole diagonalisation process as well as numerical minimization makes that somewhat hard to guarantee. In practice, you'd always end up double checking against the BR formula (as you've done) to make sure the more complex numerics are doing what you expect, which kind of undermines the point of having them.

[pathfinder49]

Is that something that comes up?

Not for our usecase. However, quadrupole clock transitions are sometimes used.

[hartytp]

Yes, but only really at 0 field. I can't think of an experiment that has used an intermediate field quadrupole clock transition, so I'm not convinced that this is something we need to cater for.

As a piece of philosophy/design goals, I think it's important to remember that we're building a framework that people can use to solve all kinds of problems, but it's okay if they have to do some of the legwork/write glue code themselves. It's fine to have some utility code to reduce boilerplate in really common cases, but we have to make sure that it's sufficiently useful to justify its place in the library as opposed to user code (or a user library).

With that in mind, it's not even clear to me that the current code should be in the utils module rather than an example. I'd argue that having some BR routines (maybe in ion_phys.breit_rabi) add extra value and justify a place in the library, but it's less clear to me if we're just sticking the main solver into a python minimizer...

[hartytp]

Not for our usecase. However, quadrupole clock transitions are sometimes used.

And, again, if you're building a clock then you really want to use the analytic routines to make sure you understand the numerics.