Assistance needed for defining excited correlated initial state for HEOM emission spectrum

[AMIT UPADHYAY]

Dear QuTiP community,

I hope you are doing well.

I am currently working on a QuTiP implementation to simulate the linear emission spectrum of a single molecule modeled using a Frenkel-exciton Hamiltonian (as shown in the attached figure). However, I am obtaining a negligible Stokes shift, which suggests that the excited system–bath correlated equilibrium state used as the initial condition for the second-stage HEOM propagation may not be constructed correctly.

In particular, I would be grateful for guidance on how to properly define the hierarchy state (including ADOs) required as the initial condition for the second HEOMSolver run, which is needed to compute the dipole–dipole autocorrelation function.

For clarity, my current procedure is as follows:

1. I start from an initial uncorrelated system–bath state and propagate the HEOMSolver until it reaches the steady state, which I use as the equilibrium state for emission calculations.

2. I then take the full set of steady-state ADOs and pre-multiply each ADO by the dipole moment operator.

3. Finally, I use this modified set of ADOs as the initial hierarchy state for a second HEOMSolver propagation to obtain the reduced density matrix dynamics and compute the dipole correlation function.

Despite this, the resulting emission spectrum shows an unexpectedly small Stokes shift.

I would greatly appreciate any suggestions or corrections regarding the proper construction of the excited correlated initial hierarchy state in this context. I have attached my QuTiP code for reference.

Thank you very much for your time and support.

Sincerely,
Amit Upadhyay

[AMIT UPADHYAY]

Dear Dr. J.R. Johansson,

This is the reminder to my previous email. I request you to suggest me something on this problem, if possible.

Thank you for your time and support.

Best regards,

Amit

--
You received this message because you are subscribed to the Google Groups "QuTiP: Quantum Toolbox in Python" group.
To unsubscribe from this group and stop receiving emails from it, send an email to qutip+un...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/qutip/ff4b134a-6a2a-4daa-a7d8-56497de3bf32n%40googlegroups.com.

[Neill Lambert]

Dear Amit,

Your procedure sounds correct to me, and I had a quick look at the notebook and what you do with the ADOs looks sensible (there's some annoying ADO transpose step which is a bit user unfriendly, but if i remember i  think you do it correctly...).  I will play a little with it and see if I can spot any issue.

when you say the stokes shift is smaller than you expected, what did you expect? are you trying to reproduce a known result from a paper?

thanks

neill

--
You received this message because you are subscribed to the Google Groups "QuTiP: Quantum Toolbox in Python" group.
To unsubscribe from this group and stop receiving emails from it, send an email to qutip+un...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/qutip/ff4b134a-6a2a-4daa-a7d8-56497de3bf32n%40googlegroups.com.

--

-----------------------------------
Senior Research Scientist
RIKEN, Japan
Research Homepage

[AMIT UPADHYAY]

Dear Dr. Neill,

Thank you for your kind response. I greatly appreciate your time and consideration, and I look forward to hearing from you further.

I expect the Stokes shift (i.e., the difference between the absorption and emission peak positions) to be approximately twice the reorganization energy. This is also illustrated in the attached figure, where the emission lineshape is shown in green and the absorption lineshape in magenta. The figure was generated using a Fortran code written by one of my colleagues.

For your reference, I am also attaching my absorption spectra code.

Best regards,
Amit

To view this discussion visit https://groups.google.com/d/msgid/qutip/CAKKp-tP26Zko-BPHOJsxK0-4UjjQEPKyUovu-sjP0am-uF57mw%40mail.gmail.com.

[AMIT UPADHYAY]

Dear Dr. Neill,

I am writing to remind you regarding my previous email and waiting for any update on the issues in my emission spectra code.

Thank you very much for your time and support.

Best regards,

Amit Upadhyay

[Neill Lambert]

Dear Amit,

Sorry for the slow reply, I was digging a little bit.  My first impression was that perhaps the initial evolution time wasn't long enough;  if you inspect the ADOs themselves, they are not really in a steadystate at 0.04 ps, you have to evolve out to about 0.8 to see the dominant ADO reach a steadystate. However, then using that as the initial condition in the subsequent system correlation function calculation leads to very weird and sometimes unstable results. I have not quite understood why yet, I will keep playing around a bit, see if there's any hidden bugs causing this.

also, probably you can reduce max_depth and Nk quite a lot for your parameters, the coupling is quite weak and the temperature is quite high. you can always check the rough accuracy by plotting the baths[0].power_spectrum() and compare it  to the original bath power spectrum.

all the best

neill

To view this discussion visit https://groups.google.com/d/msgid/qutip/b14d1e27-3537-46d1-b29d-443ae30ad105n%40googlegroups.com.

[Neill Lambert]

Dear amit,

Just fyi, I can fix the instability if I rescale your energies up by a factor of 100, and all your timescales down by a factor of 100. I guess the large units are causing some issues with the ODE solvers + the HEOM (I think the HEOM can be more sensitive to unit issues sometimes due to how the ADO matrix elements scale.  Previously we used a renormalized version which was a bit more stable, but we removed it at some point. I will put it as an issue to revert it if it makes sense to do so).

So in summary, if in your first time evolution you set:

# HEOM dynamics for the same system

tmax = 8 * ps2au      # ps to a.u.

the renormalize your units

# Unit conversion factors

cm2au = 4.556335e-6  *100       # cm-1 to atomic unit

fs2au = 41.3413745758  /100     # femtosecond to atomic unit

ps2au = 41341.3746    /100      # picosecond to atomic unit

Debye2au = 0.393430307  *100    # Debye to atomic unit

sec2au = 2.418884e-17  *100     # 1/second to atomic unit

I see a plot like this (where the shift looks more like your fortran code picture?)

thanks

neill