rateeq returns pumping rate with nonzeros different from dijq matrix

[xuguangyue]

rateeq returns pumping rate different from the dijq matrix, with impossible transitions nonzero. The reason is that when diagonalize at the magnetic field, the base changes. I think it should be changed back to the original base since what we usually calculate are based on that.

[SteveEckel]

@xuguangyue: I'm not sure what you mean. If you apply a z magnetic field, the pumping rates should be identical in shape to the dijq matrices, modulo the laser polarizations in your problem. If you apply a magnetic field in any other direction, the program rotates the polarization of the lasers into a new basis where the z axis is defined by the local magnetic field direction. This is required to calculate forces in a MOT.

If you are still convinced this is a bug, perhaps you can provide an example?

[xuguangyue]

Hi SteveEckel, thanks for responding.

If you apply a z magnetic field, the pumping rates should be identical in shape to the dijq matrices, modulo the laser polarizations in your problem.

Yes, that's what I want to get. But the result is not what I expect. Here is an example.

import numpy as np
import pylcp

# Define the atomic Hamiltonian for 87Rb:
atom = pylcp.atom("87Rb")
H_g_D2, mu_q_g_D2 = pylcp.hamiltonians.hyperfine_coupled(
    atom.state[0].J, atom.I, atom.state[0].gJ, atom.gI,
    atom.state[0].Ahfs/atom.state[2].gammaHz, Bhfs=0, Chfs=0,
    muB=1)
H_e_D2, mu_q_e_D2 = pylcp.hamiltonians.hyperfine_coupled(
    atom.state[2].J, atom.I, atom.state[2].gJ, atom.gI,
    Ahfs=atom.state[2].Ahfs/atom.state[2].gammaHz,
    Bhfs=atom.state[2].Bhfs/atom.state[2].gammaHz, Chfs=0,
    muB=1)

dijq_D2 = pylcp.hamiltonians.dqij_two_hyperfine_manifolds(
    atom.state[0].J, atom.state[2].J, atom.I)

E_e_D2 = np.unique(np.diagonal(H_e_D2))
E_g_D2 = np.unique(np.diagonal(H_g_D2))

ham = pylcp.hamiltonian(H_g_D2, H_e_D2, mu_q_g_D2, mu_q_e_D2, dijq_D2)
magField = pylcp.constantMagneticField(np.array([0., 0., 1.e-4]))
laserBeams0 = pylcp.laserBeams(
    [{'kvec':np.array([1., 0, 0.]), 's':1, 'pol':0,
     'delta':0.1+E_e_D2[1]-E_g_D2[0]}]
)

rateeq = pylcp.rateeq(laserBeams0, magField, ham)
Neq = rateeq.equilibrium_populations(
                    np.array([0., 0., 0.]), np.array([0, 0., 0]), 0.)
Rijl = rateeq.Rijl['g->e']

The above code return a Rijl different to any of the dijq_D2. The (0,3) element is nonzero, and that's impossible.

I think the problem is because you calculate the eig value at every magnetic field and sort it according to energy. We know that the $gF$ of Rb $5S{1/2}, F=1$ is negative, which means after you sort, the order of your states is different from the one you define you Hilbert space. I notice that you rotate your dijq as well. But how can we understand the result after all these rotations.

[dsbarker]

Looking at this code, I see that the polarization of the laser beam is set at pol=0, which is not supported. Looking at the code, if pol is set using a float or int, then pol>=0 is interpreted as sigma+ and pol<0 is interpreted as sigma-. I suspect that this explains the strange results.

We may want to extend __parse_constant_polarization to interpret pol=0 as pi polarization, or add better error handling to catch this case (since 0 meaning pi makes intuitive sense given behavior of pol=+/-1).

[xuguangyue]

I don't think that's the problem, though. I replace the laserBeam to the following

laserBeams0 = pylcp.laserBeams(
    [{'kvec': np.array([1., 0, 0.]), 's': 1, 'pol': np.array([0., 0., 1.]), 'pol_coord': 'cartesian',
      'delta': 0.1 + E_e_D2[1] - E_g_D2[0]}]
)

The result is the same.

[xuguangyue]

Actually, we can use sigma+/- polarization too, and the Rijl always have nonzeros values of (0,3) element, which is impossible in any case.

[SteveEckel]

@xuguangyue: This is far more subtle than I initially gave it credit for. I can confirm that your code yields an Rijl[0, 0, 3] value of about 0.1. Correct me if I am wrong, but you claim that the (0,3) element should be impossible because, presumably, the second 0 index corresponds to $|F=1,m_F=-1 \rangle$ in the ground manifold and the third index of 3 should correspond to the state $|F=1,m_F=+1 \rangle$ in the excited manifold. This is the state ordering that comes out of pylcp.hamiltonians.hyperfine_coupled. You can see that by

H_g_D2, mu_q_g_D2, g_basis = pylcp.hamiltonians.hyperfine_coupled(
    atom.state[0].J, atom.I, atom.state[0].gJ, atom.gI,
    atom.state[0].Ahfs/atom.state[2].gammaHz, Bhfs=0, Chfs=0,
    muB=1, return_basis=True)
H_e_D2, mu_q_e_D2, e_basis = pylcp.hamiltonians.hyperfine_coupled(
    atom.state[2].J, atom.I, atom.state[2].gJ, atom.gI,
    Ahfs=atom.state[2].Ahfs/atom.state[2].gammaHz,
    Bhfs=atom.state[2].Bhfs/atom.state[2].gammaHz, Chfs=0,
    muB=1, return_basis=True)

print(g_basis.T)
print(e_basis.T)

which produces

[[ 1. -1.]
 [ 1.  0.]
 [ 1.  1.]
 [ 2. -2.]
 [ 2. -1.]
 [ 2.  0.]
 [ 2.  1.]
 [ 2.  2.]]
[[ 0. -0.]
 [ 1. -1.]
 [ 1.  0.]
 [ 1.  1.]
 [ 2. -2.]
 [ 2. -1.]
 [ 2.  0.]
 [ 2.  1.]
 [ 2.  2.]
 [ 3. -3.]
 [ 3. -2.]
 [ 3. -1.]
 [ 3.  0.]
 [ 3.  1.]
 [ 3.  2.]
 [ 3.  3.]]

where the first column is $F$ and the second is $m_F$.

BUT, when you solve the equilibrium populations of the rate equations, pylcp rediagonalizes your hamiltonian at the mangetic field of interest. When that happens, the Hamiltonian states are automatically sorted according to energy. In many cases, the diagnoalization does not reorder the states, but it does for this Hamiltonian because $|F=1\rangle$ state in the ground manifold has a negative $g_F$ factor. Thus index 0 in the second dimension of the Rijl matrix actually corresponds to $|F=1,m_F=+1\rangle$ state, as that is the state with lowest energy at even small, non-zero magnetic fields.

You can see the transformation in rateeq.hamilotonian.U, which contains the transformation matrices from the original basis to the new, diagonalized basis. rateeq.hamilotonian.U[0] is the transformation matrix for ground manifold, and rateeq.hamilotonian.U[1] is the transformation matrix of the excited manifold. To be specific,

np.set_printoptions(precision=1)
print(np.real(rateeq.hamiltonian.U[0]))

produces

[[-0.0e+00  2.9e-16  1.0e+00  0.0e+00  7.7e-08 -1.6e-23  0.0e+00  0.0e+00]
 [-0.0e+00  1.0e+00 -0.0e+00  0.0e+00 -0.0e+00  8.9e-08  0.0e+00  0.0e+00]
 [ 1.0e+00  0.0e+00 -0.0e+00  0.0e+00 -0.0e+00  0.0e+00  7.7e-08  0.0e+00]
 [-0.0e+00  0.0e+00 -0.0e+00  1.0e+00 -0.0e+00  0.0e+00  0.0e+00  0.0e+00]
 [-0.0e+00 -1.2e-16 -7.7e-08  0.0e+00  1.0e+00 -2.0e-16  0.0e+00  0.0e+00]
 [-0.0e+00 -8.9e-08 -0.0e+00  0.0e+00 -0.0e+00  1.0e+00  0.0e+00  0.0e+00]
 [-7.7e-08  0.0e+00 -0.0e+00  0.0e+00 -0.0e+00  0.0e+00  1.0e+00  0.0e+00]
 [-0.0e+00  0.0e+00 -0.0e+00  0.0e+00 -0.0e+00  0.0e+00  0.0e+00  1.0e+00]]

which you can see flips the position of the $F=0,m_F=\pm1$ states. The small, non-zero components come from quadratic and higher-order Zeeman effects.

As for your choice of polarizations, with linear polarization along $\hat{z}$ (along the z axis), I find that Rijl[0,0,3] produces a value of about 0.2. For circular polarizations, it is half that. Circular polarization is defined with respect to the k-vector of the light. Because $$\hat{k} = (1, 0, 0)$$ in your example, the circular polarizations have Cartesian representation $$\epsilon^\pm = \frac{1}{\sqrt{2}}(0, 1, \pm i)$$ so the projection along $\hat{z}$ is $1/\sqrt{2}$, and that projection is squared to determine the pumping rate for a total factor of 1/2 relative to that of linear polarization along $\hat{z}$.

To summarize: this is expected behavior. The reason that Rijl does not match the dijq matrix is because of a Hamiltonian rotation necessary for solving the rate equations.

[xuguangyue]

Yes, I totally agree with you about your analysis. However, when the magnetic field is stronger, there will be not only reordering, but also state mixing between different $m_F$ levels, even between different $F$ levels. How can we know the contribution from each $$\lvert F, m_F \rangle$$ state? I think it's better to rotate back to the original base at some point.