Issue about runtime of OBE and units in general

[Vanessa Grauer]

Hello Everyone, 

First and foremost thank you for making this nice package available! 

I am currently working with the package and ran into some issues and hoped someone here might be able to help me with that...

Firstly, my goal is to rewrite both 02_F2_to_F3_1D_molasses.ipynb and 10_recoil_limited_MOT.ipynb examples for the Rb-87 cooling transition of |F=2, mF=2> to [F'=3, mF=3>

So my first issue is a rather simple one: Suppose I want to implement a Gaussian Beam with a beam diameter of 2cm @ 1/e^2. I am using the units native to PyLCP, so x is in terms of 1/k. What would be the correct way to implement the beam radius here? I am currently doing this:

pylcp.laserBeams([
        {'kvec':np.array([0., 0., 1.]), 'pol':pol[0],
         'pol_coord':pol_coord, 'delta':delta, 's':s, 'wb':1.},
        {'kvec':np.array([0., 0., -1.]), 'pol':pol[1],
         'pol_coord':pol_coord, 'delta':delta, 's':s, 'wb':1.},
        ], beam_type=pylcp.gaussianBeam)

My 2nd issue is with the script 10_recoil_limited_MOT.ipynb:

• When looking at the example  Recoil-limited MOT — pylcp 1.0.0 documentation in cell [5], I get the first error when trying to set the starting condition of the atom: "ValueError: Npop has only 4 entries for 12 states". My fix to that was the following:

tmax = 0.05/(t0)
if isinstance(eqn1, pylcp.obe):
    eqn1.set_initial_position(np.array([0., 0., 0.]))
    eqn1.set_initial_velocity(np.array([0., 0., 0.]))
    eqn1.set_initial_rho_equally()
eqn1.evolve_motion([0, tmax], random_recoil=True, progress_bar=True, max_step=1.)

now, t0 here is defined as 1/gamma and for Rb87, this is O(10³) smaller than the example with strontium, as gamma = 2*np.pi*5.75e6 Hz (See down below, how I adjusted the units in this script for Rb-87). Consequently, if I dont adjust tmax above, I get runtimes of the single atom simulation of ~12 hours, which is of course not feasible. When I'm reducing tmax, the results are not really informative I assume... Is there anything I am doing wrong here? Help would be much appreciated!

rb87 = atom('87Rb')
taulab = rb87.state[2].tau  # Lifetime of 6P_{3/2} state (in seconds)
k = 2*np.pi/780E-7    # cm^{-1}
x0 = 1/k              # our length scale in cm
gamma = 2*np.pi*5.75e6 # 5.75 MHz linewidth
t0 = 1/gamma          # our time scale in s
gradBlab = 15 # About 15 G/cm is a typical gradient for Rb

# Magnetic field gradient parameter (the factor of 3/2 comes from the
# excited state g-factor.)
alpha = cts.value('Bohr magneton')*1e-4*gradBlab*x0*t0/cts.hbar

# The unitless mass parameter:
mass = 87*cts.value('atomic mass constant')*(x0*1e-2)**2/cts.hbar/t0

# Gravity
g = -np.array([0., 0., 9.8*t0**2/(x0*1e-2)])

• Lastly, I am thinking about how to possibly implement a Pyramid MOT with the package but have no Idea where to start.. has anyone maybe already done something like that and could give me a hint?

I hope my description of the problems is sufficient enough, if not I would be glad to provide more infos on my scripts etc.

Cheers, 

Vanessa Grauer

[daniel.barker]

Hello Vanessa,

I've attached a pruned version of one of my simulation notebooks, which should help answer the units questions. It mostly has parameter definitions and a bit of 1D modelling for F=0 to F'=1 MOTs. The main things that you'll want to look at are the definition of de_scale (in the second cell) and x_scale (in the 4th cell), which convert real atomic mass to 1/mass (in pylcp units) and convert field gradient (in G/cm) to pylcp units. The 5th cell shows how to use cell 2 and 4 to generate force profiles with distances in units of the k vector. The notebook uses a rate equation model, but the units for a OBE simulation are the same.

The Npop ValueError is due to the array that is passed to obe.set_initial_pop having only 4 elements, but you now have 12 states in the Hamiltonian. Your fix should be fine, but if you'd like to have more control, just change the argument to set_initial_pop have 12 elements (which specify the diagonal elements of the density matrix with the states ordered by energy).

As you note, runtimes for OBE simulations can be very long. There are a few approaches to keeping the runtime short. The first (as you pointed out) is just to reduce tmax. In the recoil limited MOT example, tmax is 50 ms (0.05 s divided by t0 to convert to pylcp units) and that's much longer than the equilibration timescale for an alkali system. As a first guess, I would reduce tmax by roughly the ratio of the linewidths. However, since you're just getting started, I'd suggest going much shorter (maybe set tmax = 10 to tmax = 1000) to ensure that the simulations are making sense and then go to larger tmax once you have intuition about the actual equilibration time for your system. The second approach is to switch to the DOP853 integrator, which performs a bit better than the default RK45 (this is documented in two github issues: https://github.com/JQIamo/pylcp/issues/71 and https://github.com/JQIamo/pylcp/issues/44). The third approach is to define events that terminate the simulation when, for example, an atom leaves the overlap volume of the MOT beams (see the solve_ivp docs: https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html). This prevents the simulations from wasting time integrating atoms that you've already gotten the useful information from. (Parallel integration is also important, but that's already in cell [7] of the example, so I assume you're already planning on it).

I don't know of anyone that has simulated a pyramid MOT with pylcp. However, you could start by looking at our code for masked Gaussian beams and grating MOTs to see how we've implemented non-standard beam shapes and geometries. The grating MOT code is undocumented functionality, but comes with the released version of pylcp and its functions can be called by importing the subpackage:

from pylcp import gratings

You'll need to look at the underlying code to see what's going on (either dig down into your python distribution for pip install, look in your github pylcp directory for developer install, or just pull the gratings.py file from github for reference).

I hope this is helpful and let us know if you have any more questions.

Daniel