Regarding stark shifts in pylcp

188 views
Skip to first unread message

Karthick Ramanathan

unread,
Feb 5, 2022, 9:09:20 AM2/5/22
to pylcp
Hello everyone,

I am confused as to whether pylcp accounts for the light shift while dealing with multi-level atoms. 

Let's say I have a 3-level system in a lambda configuration with 2-level detuning - Delta and Raman detuning - delta. Now, for Delta>>Gamma (where Gamma is the lifetime of the excited state), we can adiabatically eliminate the excited state to give an effective 2-level system which oscillates with a generalised Rabi-frequency of :

Omega_r = sqrt(Omega_eff + (delta-delta_AC))

where Omega_eff = Omega1*Omega2/(2*Delta) and del_AC is the differential stark shift. 

Now if I construct this 3 level system with 2 lasers (similar to the construction in the three-level susceptibility example given in the documentation), I will observe Rabi flopping between the two ground states. However in practice, this rabi flopping frequency will be crucially dependent on the ratio of intensities of the 2 laser beams, as the ratio decides the differential AC stark shift. 

 My question is: if I evolve the 3-level hamiltonian with 2 laser beams using the OBEs, and plot the populations (by taking the diagonal entries of the solution rho matrix) as a function of time, will I see oscillation with a frequency of Omega_r or just Omega_eff? I am not able to see the sensitivity of the solution to the intensity ratio as is expected theoretically. Am I missing something here, or am I supposed to calculate the stark shift by hand and include it in the Hamiltonian?

Any help would be appreciated.

Thanks,
Karthick

Sheng-wey Chiow

unread,
Feb 8, 2022, 1:21:36 AM2/8/22
to Karthick Ramanathan, pylcp
Hi Karthick,

As you said, AC Stark shifts come from the adiabatic elimination of the excited state. So the familiar effect of AC Stark shifts will show up only when the detuning from the excited state is large enough that the fast oscillations at the detuning frequency can be averaged out. PyLCP doesn't handle large detunings well, since it solves the OBE by brute force. 

That said, for small detunings that PyLCP can handle, if you have each laser couples to both 1-3 and 2-3 transitions where 3 being the excited state, you should see a fuzzy Rabi oscillation at Omega_eff. The AC Stark shifts really come from the coupling to off-resonant states. If you have laser 1 only couples to 1-3 and laser 2 only couples to 2-3, then you don't have AC Stark shifts.

Hope that makes sense.

Sheng-wey

--
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/888e5ac7-d7ff-4fe2-a81a-9a859feec335n%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Eckel, Stephen P. (Fed)

unread,
Feb 8, 2022, 9:15:04 PM2/8/22
to Sheng-wey Chiow, Karthick Ramanathan, pylcp

Hi Kartick and Sheng-way,

 

Thank you for your questions and comments; this particular set got me back into running some sims.

 

Let me start off by saying the pylcp can handle large detunings fine, provided that you define the rotating frame properly.  In the three_level_susceptibility.ipynb example, the detunings are all written into the Hamiltonian (i.e., the rotating frame is set by the fields), which is the most computationally efficient.  You can set Delta (the two-level detuning) pretty large.  In what follows, I’ve set it to -100 with no serious computational effect.

 

Also, any effect that is within the context of your Hamiltonian should pop out in pylcp.  Differential Stark shifts in three-level systems are a simple consequence of the Hamiltonian in the far-detuned limit, and we can make them rather obvious pretty easily.

 

To do so, let’s modify the three_level_susceptibility.ipynb example.  First, set Delta=-100 (two-level detuning) and delta=0 (three-level detuning) in the cell that defines the Hamiltonian and execute it.

 

I then set the Ige=Ire=400, which makes Omega_re/Gamma=Omega_ge/Gamma=10, which makes the effective three level Rabi rate equal to Omega_eff/Gamma=0.5 (because Omega_eff=Omega_re*Omega_ge/2/Delta).  After executing that cell, you can print out the structure of the three-level Hamiltonian for this configuration to verify the Rabi rates:

 

Eq = {}

for key in laserBeams:

    Eq[key] = laserBeams[key].total_electric_field(np.array([0., 0., 0.]), 0.)

hamiltonian.return_full_H(Eq, np.array([0., 0., 0.]))

 

which yields

 

array([[  0.+0.j,   0.+0.j, -10.+0.j],

       [  0.+0.j,   0.+0.j, -10.+0.j],

       [-10.+0.j, -10.+0.j, 100.+0.j]])

 

If you then execute the “Damped Rabi Oscillations” cell, it shows slightly damped Rabi oscillations with a period of Gamma t = 4 pi, as expected.

 

You can modify this example to see the differential Stark shift rather easily.  Add the following code that re-runs the damped Rabi oscillation calculation but changes the relative intensity between the two lasers, keeping the product of the intensities the same.

 

fig, ax = plt.subplots(1, 1)

 

for factor, ls in zip([1., 3., 10.], ['-','--','-.']):

    laserBeams = return_three_level_lasers(factor*Ige, Ire/factor)

    obe = pylcp.obe(laserBeams, magField, hamiltonian,

                    transform_into_re_im=True)

    obe.set_initial_rho_from_populations(np.array([0., 1., 0.]))

    obe.evolve_density([0, 10*np.pi])

   

    ax.plot(obe.sol.t/2/np.pi, np.real(obe.sol.rho[0, 0]), c='C0', ls=ls ,lw=0.5, label='$\\rho_{gg}$')

    ax.plot(obe.sol.t/2/np.pi, np.real(obe.sol.rho[1, 1]), c='C1', ls=ls, lw=0.5, label='$\\rho_{rr}$')

   

ax.set_xlabel('$\Gamma t/2/\pi$')

ax.set_ylabel('$\\rho_{ii}$')

 

My plot looks like this:

 

Blue is rho_{gg}, orange is rho_{rr}; solid is equal intensities (resonant), dashed is a factor of 9 between the two, and the last dash-dot is a factor of 100 between the two.  I haven’t quantitatively compared to theory, but it certainly seems that by changing the relative intensities, I make Rabi oscillations that look more and more detuned.  That is exactly the differential Stark shift you were looking for.

 

Hope that helps!

 

-Steve

 

-- 

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

Sheng-wey Chiow

unread,
Feb 8, 2022, 11:08:35 PM2/8/22
to Eckel, Stephen P. (Fed), Karthick Ramanathan, pylcp
Hi Steve,

Thanks for the clarification. You are right that there are AC Stark shifts in the setting and pylcp is capable of capturing that.

My comment about large detunings came from a different application. I tried to implement rectified dipole forces, which essentially are the spatial variations of AC Stark shifts. I don't remember finding a good computationally efficient solution in pylcp.

Cheers,
Sheng-wey

Eckel, Stephen P. (Fed)

unread,
Feb 9, 2022, 8:55:43 AM2/9/22
to Sheng-wey Chiow, Karthick Ramanathan, pylcp

Hi Sheng-way,

 

If you are still interested in solving your problem with pylcp, let me know.  With more details about your particular problem, I might be able to suggest a setup that is computationally efficient.

 

Feel free to email me directly.

 

-Steve

Karthick Ramanathan

unread,
Feb 11, 2022, 12:20:32 PM2/11/22
to pylcp

Hi Stephen and Sheng-wey,

Thanks a lot for your detailed explanation and comments. It was quite illuminative.

The system I am dealing with is a bit more complicated and I discovered a few flaws in my code which I am trying to fix now. Nonetheless, I could benefit from your comments on my approach to the problem in general.

I am working on an experiment to perform stimulated Raman transitions between the hyperfine ground states through the D2 line of Rb87, and would like to understand the process theoretically. I initially have an ensemble in F=2 with populations equally distributed across the sublevels. To stimulate the transition, I have two lasers from F=1 and F=2 which are far detuned from the excited state manifold. To model this in pylcp, I constructed the Hamiltonian from scratch with the hamiltonians.singleF() for 4 F states – F=1,2 in 5S1/2 and F’=1,2 in 5D3/2. Since I am interested in only stimulated processes, I decided not to use hyperfine_coupled() as it would also contain F’=0,3, each of which is forbidden transitions from one of F=1 and F=2, and I didn’t want to unnecessarily increase the size of the matrices. For the dipole coupling, I used dqij_two_bare_hyperfine() for all 4 combinations from F=1,2 to F’=1,2. I gave the appropriate detunings in the rotating frame similar to the three-level-susceptibility example. The construction of the laserBeams also follows the example, only I now have ‘four’ laser beams in the model. I then evolve density using the OBEs and plot my populations.

I then tried comparing my results with the theoretical rabi frequency, which raised several doubts, both regarding the soundness of my construction, as well as pylcp’s computations. I went down the rabbit hole of figuring out how pylcp computes the rabi frequency, from the source code. From my understanding, the rabi frequency connecting any two magnetic sublevels between the ground (g) and excited (e) state, is given by (from Eq. (35) and (36) of Rubidium 87 Line D Data - https://steck.us/alkalidata/):

omega.png

I believe the last two terms (along with the J=1/2->3/2 transition dipole matrix element) in the above expression should constitute ‘d_q’. However, I only find the last term, in source code for dqij_two_bare_hyperfine(). I suppose this was intended - the F states are not associated with any particular J in this function. But, even in the dqij_two_hyperfine_manifolds(), I find a term missing - \sqrt(2J+1):

def matrix_element(J, F, m_F, Jp, Fp, m_Fp, I, q):

        return (-1)**(F-m_F+J+I+Fp+1)*np.sqrt((2*F+1)*(2*Fp+1))*wig3j(F, 1, Fp, -m_F, q, m_Fp)*wig6j(J, F, I, Fp, Jp, 1)

Further, pylcp computes the amplitude of ‘E_q’ from the saturation parameter ‘s’, as \sqrt{2*s}. So, the first part of my expression also appears in 0.5*\sum_q{d_qE_q}.

So, with all this in mind, a few questions/comments:

1.       Am I right about the missing \sqrt{2J+1} factor? I suppose it would get cancelled during normalization, but is the absolute value of the rabi frequency correct as well?

2.       I feel it would be helpful to include functionality where I can define dqij only select states in the hyperfine manifold. In my opinion, this would help both in saving computational cost, as well as understanding effects caused by coupling between select levels in isolation. I coded everything with dqij_two_bare_hyperfine(), and later realized the values would be incorrect considering there is no J information given. Again, am I missing something in the pylcp toolbox which does things like I require?

3.       In my construction, I would have several sets of magnetic sublevels which would form 3-level systems and would rabi-flop at their own frequencies. The behavior of F=1 and F=2 as a whole, would be a superposition of all these oscillations. But again, each of the sub-levels would have their own stark shifts. In our experiment, we wish to calculate and set the ratio of intensities of the 2 Raman beams at value which cancels out the stark shift, which would be possible for only one set of sub-levels. We are dealing with effective rabi frequencies of the order of kHz, and at one intensity ratio, theoretically, I expect the other sets of sublevels to be considerably farther detuned from Raman resonance and hence execute oscillations which doesn’t transfer populations as effectively. I am not really able to see this effect in the results. Do you have any thoughts on this? Experimentally, I have intensity ratio fluctuations in my system, and I want to understand how my oscillations get affected theoretically, and how much fluctuations can be tolerated.

On the whole, I would really appreciate getting some inputs/suggestions on the problem, and my approach towards modelling it in pylcp.

Thanks a lot for your time!!

Best regards,
Karthick

Sheng-wey Chiow

unread,
Feb 11, 2022, 1:43:31 PM2/11/22
to Karthick Ramanathan, pylcp
Hi Karthick,

1. The definition of Rabi frequency has several conventions, so I didn't check which one is used in pylcp or if it's self-consistent. Even I_sat has several values in Steck. I wouldn't worry about the absolute value
2,3. I think once the hamiltonian and the laser fields are set up, pylcp will calculate all couplings based on the Rabi frequency expression. You don't have to set up each 3-level system by hand. Physics-wise, I'd worry more about differential Zeeman shifts between different couplings than differential AC Stark shifts, unless you are considering a very quiet magnetic environment. Even if all transitions are on resonance, different 2-photon Rabi frequencies will mess up the detection if you only detect F=1 vs F=2, not to mention potential couplings between all 3-level systems due to polarization imperfections.

Sheng-wey

Karthick Ramanathan

unread,
Feb 11, 2022, 7:45:48 PM2/11/22
to Sheng-wey Chiow, pylcp
Hi Sheng-wey,

Yes, the zeeman shift would definitely be more of a concern, but we have compensated for the stray magnetic fields upto a few tenths of mG. The reason I am considering stark shifts, is because at the laser intensities we are dealing with, the stark shifts are of the same order of magnitude as the 2-photon rabi frequency. An intensity ratio not tuned to cancel this out would seriously impact the resonance. I want to study this effect properly in isolation. 

Experiment wise,  as the Raman beams' pulse duration is varied, I do think some semblance of oscillation should be visible with F=1 vs 2 detection itself, before they get washed away due to incoherence. One obvious way to simplify life is to perform magnetic state selection and apply a bias field to shift the undesired sublevels far out of 3-level resonance. But we aren't there yet. 

Thanks a lot for your reply!

Best regards,
Karthick.

Eckel, Stephen P. (Fed)

unread,
Feb 15, 2022, 11:10:03 PM2/15/22
to Karthick Ramanathan, Sheng-wey Chiow, pylcp

Hi Karthick,

 

Thank you for your questions.  Let me try to answer them in turn:

  1. You are strictly speaking correct about the missing sqrt(2 J+1) term.  But note that the d_q matrices you are using are normalized: the sum of the squares of all the d_q elements along a column and over all three polarizations sum to one.  This ensures that the decay rate of all the states in the same excited manifold is equal.  Thus common factors, like the sqrt(2 J+1), are normalized away.

    I suspect where your code is going wrong is that if you blindly construct a three-level system out of two normalized d_q matrices, you may unwittingly double your decay rate and, I think, quadrupole your saturation intensity.

    In this case, it is better to calculate the full D2 line first, then manually divide up your ground state into two different manifolds (dropping the Zeeman terms connecting them), manually divide the d_q, then put it together into your Hamiltonian class:

 

H_g, mu_g = pylcp.hamiltonians.hyperfine_coupled(1/2, 3/2, 2, 0, 100)

H_e, mu_e = pylcp.hamiltonians.hyperfine_coupled(3/2, 3/2, 3., 0, 10)

d_q = pylcp.hamiltonians.dqij_two_hyperfine_manifolds(1/2, 3/2, 3/2, normalize=True)

 

hamiltonian = pylcp.hamiltonian(mass=mass)

hamiltonian.add_H_0_block('g', 0.*np.eye(3))

hamiltonian.add_mu_q_block('g', muq_g[:, 0:3, 0:3])

hamiltonian.add_H_0_block('r', delta*np.eye(5))

hamiltonian.add_mu_q_block('r', muq_g[:, 3:, 3:])

   

hamiltonian.add_H_0_block('e', H_e) # Add two-level detuning properly

hamiltonian.add_mu_q_block('e', muq_e)

hamiltonian.add_d_q_block('g','e', d_q[:, :3, :])

hamiltonian.add_d_q_block('r','e', d_q[:, 3:, :])

 

  1. You can always restrict the space once you’ve created the matrices by simply dropping the appropriate rows and columns of the matrices before joining them together into the Hamiltonian class.  For example, in a 87-Rb-like D2 manifold, you can drop the F’=0 state by doing the following:

 

hamiltonian.add_H_0_block('e', H_e[1:, 1:]) # Add two-level detuning properly

hamiltonian.add_mu_q_block('e', muq_e[:, 1:, 1:])

hamiltonian.add_d_q_block('g','e', d_q[:, :3, 1:])

hamiltonian.add_d_q_block('r','e', d_q[:, 3:, 1:])


In my experience though, if a completely unconnected state is left in the Hamiltonian, you don’t take a big performance hit by leaving it in.

  1. Unfortunately, I don’t have any additional insight into your particular problem then Sheng-wey has already given you.  But it makes me think the best way to approach the problem is to build it up slowly: you can restrict your space to a problem you can understand, then slowly work your way up to the problem you don’t.

Reply all
Reply to author
Forward
0 new messages