questions about units

[Yuna Park]

Hello!

Thank you for distributing the laser cooling scripts !

I was trying the real atoms example with the magnetic forces turned on when computing the rate equation. I see that the output force profile is discontinuous when alpha = 1. The real Hamiltonian is more discontinuous and oscillating, while the 2->3 Hamiltonian is also discontinuous at the center (attached a screen shot).  If this is an expected behavior, could you explain why this is the case? 

I also wanted to try using real units to understand the script better (positions in cm, time in seconds, field in gauss). I tried running the following commands but encountered an error when plotting the force that says "unable to find solutions." 

atom_name = '7Li';

BGcm=10; #G/cm

atom0 = atom(atom_name)

klab = 2*np.pi*atom0.transition[1].k # Lab wavevector (without 2pi) in cm^{-1}

tau = atom0.state[2].tau #seconds

gamma = 1/tau;

k=klab; t0=1; x0=1; #x0 in cm, t0 in s

mass = atom0.mass #kg

alpha = BGcm #G/cm

muB = cts.value('Bohr magneton')*1e4*1e-4 # J/T to cm^2/s^2/G

s=1.3; det=-1e3

H_g_23, mu_q_g_23 = pylcp.hamiltonians.singleF(F=2, gF=1/2, muB=muB)

H_e_23, mu_q_e_23 = pylcp.hamiltonians.singleF(F=3, gF=2/3, muB=muB)

dijq_23 = pylcp.hamiltonians.dqij_two_bare_hyperfine(2, 3)

hamiltonian_23 = pylcp.hamiltonian(H_g_23, H_e_23, mu_q_g_23, mu_q_e_23, dijq_23,mass=mass,k=k,gamma=gamma,muB=muB)

# Define the laser beams for 2->3

laserBeams_23 = pylcp.conventional3DMOTBeams(s=s, delta=det,k=k)

magField = pylcp.quadrupoleMagneticField(alpha)

# Make the trap for 2->3

eqn = pylcp.rateeq(

    laserBeams_23, magField, hamiltonian_23, include_mag_forces=True)

Thank you!

Yuna

[Barker, Daniel S. (Fed)]

Hello Yuna,

 

Thanks for using the package!

 

I think both issues that you’re seeing is caused by mismatched units. The https://python-laser-cooling-physics.readthedocs.io/en/latest/examples/MOTs/06_real_atoms_3D_MOT.html example uses non-standard (for pylcp) length units: setting alpha = 1 and muB = 1 puts lengths in units of 1/(muB*B’/hbar*gamma). Those length units work fine for force profiles using the rate equations, but many pylcp functions expect distances to be specified in units of 1/k. I believe that the magnetic forces are calculated assuming distances are in 1/k, so setting alpha = 1, muB = 1 makes a gradient that shifts the transition by 1 linewidth per wavelength of distance (absurdly large!).

 

In scripts where I want magnetic forces (or to allow atoms to move), I always use 1/k distance units to avoid this kind of issue. Looking at your physical units version of the example, I think you’ll get reasonable results if you set muB = 1 in the Hamiltonian, use the cts value of muB in Hz/T (*1e-4 to convert to G), and then set alpha = BGcm*100*muB/(gamma/(2*np.pi)*k).

 

Let us know if you have any more questions (or if you’re still getting strange results).

 

Best,

 

Daniel

--
You received this message because you are subscribed to the Google Groups "pylcp" group.
To unsubscribe from this group and stop receiving emails from it, send an email to pylcp+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/pylcp/a021fdb5-ca30-475a-a0ef-fa6813477790n%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Eckel, Stephen P. (Fed)]

Thanks Shang-wey,

I concur with your changes and I am adding them to the forthcoming 1.0.3 bug fix release.

-Steve

--

Dr. Stephen Eckel
Physicist
Sensor Sciences Division
Fundamental Thermodynamics Group
100 Bureau Drive, Stop 8364
Gaithersburg, MD 20899-8364
(301) 975-8571

On Apr 8, 2024, at 4:33 PM, Sheng-wey Chiow <swc...@gmail.com> wrote:

Hi Daniel and all,

To follow up on the units issue: I tried to run and modify real_atoms_3D_MOT using different units. Initially it didn't work as expected. Here are some modifications I made to have consistent results irrespective of units. I only worked with rateeq.

1. Modify line 378 of rateeq.py from Neq = np.compress(S <= self.svd_eps, VH, axis=0).T to  Neq = np.compress(S/np.max(S) <= self.svd_eps, VH, axis=0).T, so that the script is robust against large values of gamma.

2. Modify line 516 of hamiltonian.py from block.matrix to  block.matrix,k=block.parameters['k'],gamma=block.parameters['gamma'], so that the rotated_hamiltonian inherits the correct gamma and k.

With these, the script can run with numerical values in SI units, or other units. Any feedback is welcome.

Cheers,

Sheng-wey

To view this discussion on the web visit https://groups.google.com/d/msgid/pylcp/f61ef05f-a2f7-405b-8aba-fdad9cdf10c0n%40googlegroups.com.