Simulating several manifolds

[Roberto Haubold]

Hi everyone,

I was just wondering if anyone has an example of setting up rate equations for more than 2 manifolds where the different transitions have different decay rates. It seems all the examples work in units that set Gamma = 1, so things get confusing when there is more than one Gamma. 

Best,

Roberto

[Barker, Daniel S. (Fed)]

Hi Roberto,

 

I don’t have an example that does this (though I’m working on a similar problem currently).

 

I suggest starting with the EIT example, which shows how to construct the Hamiltonian “by hand” with the add block functions: Three-level susceptibility & EIT — pylcp 1.0.0 documentation. I’ll see if I can dig up some old code that does EIT with hyperfine structure (as a more advanced starting place).

 

The add_d_q_block function has a keyword arguments gamma and k (both equal to 1 by default): Hamiltonian Class — pylcp 1.0.0 documentation. Start by picking a “main” manifold that will have gamma = 1, k = 1. Add “secondary” manifolds by specifying gamma != 1, k !=1 when you add the d_q block for that manifold. I believe the correct gamma and k will be the ratio of (e.g.) gamma for the secondary manifold to gamma for the main manifold.

 

Let us know if you have additional questions,

 

Daniel

 

From: py...@googlegroups.com <py...@googlegroups.com> On Behalf Of Roberto Haubold
Sent: Friday, February 6, 2026 8:21 PM
To: pylcp <py...@googlegroups.com>
Subject: [EXTERNAL] Simulating several manifolds

 

You don't often get email from roberto....@gmail.com. Learn why this is important

--
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 visit https://groups.google.com/d/msgid/pylcp/b3af0677-f719-4468-a616-b2adae68f7ebn%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Roberto Haubold]

Hi Daniel,

Thank you for your message. I will try again to use the add block functions and see if I can get things to work.

Best,

Roberto

[Roberto Haubold]

Hi Daniel,

I think I've understood what's going on. In rateeq.py, line 264, the pumping rates are calculated for each laser:

self.Rijl[key][ll] = gamma*intensity/2*\

fijq/(1 + 4*(-(E2 - E1) + delta - np.dot(kvec, v))**2/gamma**2)

This calculation follows equation (2) in  http://dx.doi.org/10.1088/1367-2630/17/1/015007:

But pylcp is missing the 2*pi/lambda that appears in the Tarbutt paper -- that is, k is assumed to be 1. This will give incorrect results when adding manifolds with k !=1.

I propose the following modifications:

Add a line below line 240 of rateeq.py:

from 

gamma = self.hamiltonian.blocks[ind].parameters['gamma']

to

gamma = self.hamiltonian.blocks[ind].parameters['gamma']

k = self.hamiltonian.blocks[ind].parameters['k']

then change line 265:

from

self.Rijl[key][ll] = gamma*intensity/2*\

fijq/(1 + 4*(-(E2 - E1) + delta - np.dot(kvec, v))**2/gamma**2)

to

self.Rijl[key][ll] = gamma*intensity/2*\

    fijq/(1 + 4*(-(E2 - E1) + delta - np.dot(kvec, v)*k)**2/gamma**2)

I've attached an example where we can explicitly change the unit system by selecting characteristic scales for the energy and wavevector. The equilibrium populations printed at the end should not depend on the unit system we choose. This works with the modifications I suggest above, but not with standard pylcp.

I should also mention that I am also using the modifications to the source code suggested in this previous post:

https://groups.google.com/g/pylcp/c/67nvnt4klfE/m/3rYcjH44AAAJ

I have not tried this without them.

Apologies if I have any of this wrong.

Best wishes,

Roberto

[Barker, Daniel S. (Fed)]

Hi Roberto,

 

This is very timely as we’re trying to get version 1.0.3 completed. I think that I agree with your changes and I’ll add them to the 1.0.3 release once we can get them tested.

 

Best,

 

Daniel

 

From: py...@googlegroups.com <py...@googlegroups.com> On Behalf Of Roberto Haubold
Sent: Wednesday, February 11, 2026 1:42 PM
To: pylcp <py...@googlegroups.com>
Subject: Re: [EXTERNAL] Simulating several manifolds

 

You don't often get email from roberto....@gmail.com. Learn why this is important

Hi Daniel,

 

I think I've understood what's going on. In rateeq.py, line 264, the pumping rates are calculated for each laser:

 

self.Rijl[key][ll]

=

gamma*intensity/2*\

 

fijq/(1

+

4*(-(E2

-

E1) +

delta

-

np.dot(kvec,

v))**2/gamma**2)

 

 

This calculation follows equation (2) in  http://dx.doi.org/10.1088/1367-2630/17/1/015007:

To view this discussion visit https://groups.google.com/d/msgid/pylcp/4ef1356e-d961-4d28-889f-41bd1416a63bn%40googlegroups.com.

[Lajos Palanki]

Hi All,

I would strongly advise against implementing this solution. The way to have k!=1 in pylcp is to have the norm of the kvec of the laser to not equal 1, but rather k in the simulation frame. This will make the rate calculation include the "missing" 2pi/lambda term.

If in the current state this change would be added to pylcp it would break all existing code, that already works by scaling the kvec of the lasers (even if that way is not intended). Apart from not breaking old code, there is a bigger reason not to make this change. The OBE's don't explicitly calculate the Doppler shift, but the Doppler shift enters through the phase evolution of the states and the laser field. If the "solution" proposed in the e-mail chain is implemented, then this would mean that for the same Hamiltonian and laser fields the OBEs and the rate equations will be simulating different systems. This is hardly the expected behaviour.

I think the better solution to the problem is to expand the documentation to discuss how to set the kvec for when k != 1 and include this in examples.

Best,

Lajos

[Barker, Daniel S. (Fed)]

Hi Lajos,

 

Thank you for pointing this out. I was wondering if scaling kvec was the correct solution overnight.

 

Rest assured, we would not push a change like this without making sure that it works.

 

Daniel

 

From: py...@googlegroups.com <py...@googlegroups.com> On Behalf Of Lajos Palanki
Sent: Thursday, February 12, 2026 4:28 AM
To: pylcp <py...@googlegroups.com>
Subject: Re: [EXTERNAL] Simulating several manifolds

 

You don't often get email from lal...@gmail.com. Learn why this is important

To view this discussion visit https://groups.google.com/d/msgid/pylcp/66565b17-66cf-4598-b8ca-7207100a2fc2n%40googlegroups.com.

[Roberto Haubold]

Hi Lajos and Daniel,

Thank you for the help! I agree that scaling kvec is the correct solution and I've used this to get the simulations working. I think some additional documentation explaining this might help other users in future.

Best,

Roberto

[Lajos Palanki]

Hi Roberto,

Sorry for not addressing the original question previously, but I wanted to double-check some things first.

If you want to add a new manifold that doesn't share gamma with the main system you need to tune dijq rather than gamma. This is most easily seen by examining the behaviour of a system with a large hyperfine splitting in the ground state.

This system we can set up either with the coupled hyperfne or by treating the two hyperfine manifolds of the ground state separately as separate pylcp manifolds (especially useful if the separation is large). Let's call the manifolds "g1", "g2" and "g" for the combined ground state. For the combined case g["gamma"] = 1 (k we will ignore, because in the hamiltonian it is used only for recoil, and in the lasers' case it was discussed above). In this case dijq will be normalised such that for all excited states summing over all decay paths dijq[e->q]**2 will evaluate to 1 (i.e. 

np.sum(np.abs(dq)**2,axis=(0,1))=[1]*dq.shape[2]

).

Considering the case of the split manifolds. g1["gamma"] and g2["gamma"] must equal g["gamma"] as otherwise the shape of the resonance would not match the case of the joint system. Splitting dijq however does mean that summing over all decay paths still evaluates to 1, but this is not true for all decay paths within a manifold. Specifically this will mean that the sum within a manifold will equal the branching ratio to the said manifold.

Applying this to a generic lambda system (and borrowing some code from the EIT example) one can implement it this way:

Delta = -2; delta = 0.

Gamma1 = 0.9

Gamma2 = 0.1

H0 = np.array([[1.]])

mu_q = np.zeros((3, 1, 1))

d_q = np.zeros((3, 1, 1))

d_q[1, 0, 0,] = 1

def return_three_level_hamiltonian(Delta, delta):

    hamiltonian = pylcp.hamiltonian()

    hamiltonian.add_H_0_block('g', 0.*H0)

    hamiltonian.add_H_0_block('r', delta*H0)

    hamiltonian.add_H_0_block('e', -Delta*H0)

    hamiltonian.add_d_q_block('g','e', d_q*np.sqrt(Gamma1), gamma = Gamma1 + Gamma2)

    hamiltonian.add_d_q_block('r','e', d_q*np.sqrt(Gamma2), gamma = Gamma1 + Gamma2)

Hopefully this helps.

Best,

Lajos