Question regarding generate_force_profile in OBE

[cc c]

Hi,

Thank you for this amazing package.

I am trying to generating a force profile using OBE, but I am confused with how pylcp handles this process . In obe.py, function generate_force_profile use find_equilibrium_force to calculate the force, in line 1120 of the code, why does pylcp need to use  set_initial_position_and_velocity ?  The equilibrium force should be related to a specific position and velocity, so these values should remain unchanged. However, in line 1120, only the velocity is kept constant (with acceleration set to 0), while the position is updated by v*t. After multiple iterations, how is the force related to the initial position?  Would it make more sense to comment out lines 1120 and 1121?

Best,

Chen

[Lajos Palanki]

Hi,

As an extension of this I thought I'd share my own observations and testing relevant to this topic.

First I want to highlight a (possible?) bug in the determination of deltat. The nested function default_deltat (obe.py:1166), copied here, does not seem to honor deltat_tmax if neither deltat_v nor deltat_r are supplied. I highlight this as testing might be affected by it.

        def default_deltat(r, v, deltat_v, deltat_r, deltat_tmax):

            deltat = None

            if deltat_v is not None:

                vabs = np.sqrt(np.sum(v**2))

                if vabs==0.:

                    deltat = deltat_tmax

                else:

                    deltat = np.min([2*np.pi*deltat_v/vabs, deltat_tmax])

            if deltat_r is not None:

                rabs = np.sqrt(np.sum(r**2))

                if rabs==0.:

                    deltat = deltat_tmax

                else:

                    deltat = np.min([2*np.pi*deltat_r/rabs, deltat_tmax])

            return deltat

Secondly, I have implemented the suggested change, by introducing reset_pos and replacing 1120 with:

                self.set_initial_position_and_velocity(self.sol.r[:, 0 if reset_pos else -1],

                                                       self.sol.v[:, 0 if reset_pos else -1])

The tests were performed without absolute or relative tolerances set to 0 so no early termination was possible. I attach the notebook I used, but for a toy system (F=2->F'=2, Gamma/2pi=6MHz,k=2pi/800nm, code attached) these are the results:

It seems to converge to something at high enough deltat. The resetting of the position seems to eliminate some sort of deviation that the original suffers from at long runtime. This issue is not solved however by resetting the position but is merely eliminated if deltat*itermax is too high. There is perhaps a way to use this in a way that auto-finds deltat and itermax based on the absolute and relative tolerances alone. It is somewhat late here, so I'll leave this here, but will revisit it next week.

Best,

Lajos

[Eckel, Stephen P. (Fed)]

Hi Lajos,

Sorry for the (super) delayed reply.  I haven’t had much time recently to devote to pylcp code. Before I address your observations, a little bit of detail about how find_equilibrium_force works: it evolves the OBE for some deltat and compute the average force over that deltat. It then continues the evolution further for another deltat, recomputing the force over that new time block, and checking for convergence with the last time block. Physically, the particle is being “dragged” through the light field to compute its force, and once the force stabilizes over a time block deltat, the force is declared to be the “equilibrium” force. There are probably a nearly infinite number of ways to do this better.

But, given the way it is coded, the selection of deltat is paramount.  The best way to pick your deltat is to provide a functio through the delta_func keyword argument that takes arguments r_initial and v_initial, and computes the relevant deltat.  By default, there is a function to do that, the nested default_deltat() function to which you refer. This function assigns deltat as either deltat_r/|r| or deltat_v/|v|, depending on whether delta_r or delta_v is provided, or the maximum deltat_max.  Use of deltat_max is motivated by the fact that you might get a v or r magnitude really close to zero, making deltat obsurdly large, larger than what is necessary for convergence.  Another possibility is that you can bypass this entire thing by setting deltat_func=None and providing deltat.

So, as to your first observation: yes, delta_tmax is only honored if delta_v or delta_r is provided.  If this is the use case though, you will get an error later in the code because deltat returned by default_deltat will be None.

None of this is in the documentation and clearly it should be; I will open an appropriate issue on GitHub.

Second, I think you are introducing a bug by reseting the position as you do. If you reset the position, then you effectively reset the spatial phase of the light field the atom last experienced in the time block to the spatial phase the atom first experienced when it started in that time block. This will wash out the sub-Doppler forces, as can be seen by your force profiles with reset_pos=True. Physically, it is akin to some sort of phase noise that is resetting the phase of the light field after each time block.   As you observe, it will probably converge faster though, because now you only have to have the density matrix itself evolve toward some sort equilibrium.

Hope that helps,

-Steve

--

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

On Oct 25, 2024, at 4:22 PM, Lajos Palanki <lal...@gmail.com> wrote:

Hi,

As an extension of this I thought I'd share my own observations and testing relevant to this topic.

First I want to highlight a (possible?) bug in the determination of deltat. The nested function default_deltat (obe.py:1166), copied here, does not seem to honor deltat_tmax if neither deltat_v nor deltat_r are supplied. I highlight this as testing might be affected by it.

        def default_deltat(r, v, deltat_v, deltat_r, deltat_tmax):

            deltat = None

            if deltat_v is not None:

                vabs = np.sqrt(np.sum(v**2))

                if vabs==0.:

                    deltat = deltat_tmax

                else:

                    deltat = np.min([2*np.pi*deltat_v/vabs, deltat_tmax])

            if deltat_r is not None:

                rabs = np.sqrt(np.sum(r**2))

                if rabs==0.:

                    deltat = deltat_tmax

                else:

                    deltat = np.min([2*np.pi*deltat_r/rabs, deltat_tmax])

            return deltat

Secondly, I have implemented the suggested change, by introducing reset_pos and replacing 1120 with:

                self.set_initial_position_and_velocity(self.sol.r[:, 0 if reset_pos else -1],

                                                       self.sol.v[:, 0 if reset_pos else -1])

The tests were performed without absolute or relative tolerances set to 0 so no early termination was possible. I attach the notebook I used, but for a toy system (F=2->F'=2, Gamma/2pi=6MHz,k=2pi/800nm, code attached) these are the results:

<params.png>

It seems to converge to something at high enough deltat. The resetting of the position seems to eliminate some sort of deviation that the original suffers from at long runtime. This issue is not solved however by resetting the position but is merely eliminated if deltat*itermax is too high. There is perhaps a way to use this in a way that auto-finds deltat and itermax based on the absolute and relative tolerances alone. It is somewhat late here, so I'll leave this here, but will revisit it next week.

Best,

Lajos

On Saturday 5 October 2024 at 13:37:26 UTC+1 cc c wrote:

Hi,

Thank you for this amazing package.

I am trying to generating a force profile using OBE, but I am confused with how pylcp handles this process . In obe.py, function generate_force_profile use find_equilibrium_force to calculate the force, in line 1120 of the code, why does pylcp need to use  set_initial_position_and_velocity ?  The equilibrium force should be related to a specific position and velocity, so these values should remain unchanged. However, in line 1120, only the velocity is kept constant (with acceleration set to 0), while the position is updated by v*t. After multiple iterations, how is the force related to the initial position?  Would it make more sense to comment out lines 1120 and 1121?

Best,

Chen

--
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/f37a08e5-8807-4e71-a22f-587f02c50d11n%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
<obe_test.ipynb><params.png>