Steady state within Bloch-Redfield framework?

[javierde...@gmail.com]

Hi all. Firstly I would like to appreciate all the work done for this quantum dynamics computing libraries and the impressive quality of it as well as the documentation available.

It has become a tool I daily use during my PhD.

   I had a question; perhaps too trivial. I'm using the Bloch-Redfield Master equation solver because I have some non-flat-in frequency correlation functions for a bath that describes a non markovian evolution; 

I need to compute the absorption spectrum of a system with a given Hamiltonian and correlation functions, so I need the steady state density matrix of this system (with maybe an external driving term),

 or at least, some of its components (which I'm able to obtain as the expected value of some projection operator). 

Since I see that the Bloch-Redfield equation has no associated routine to extract a steady state value of a given operator, I guess if the steady routine that gives you the steady density matrix for a given Linblad superoperator, 'steady(L)'

can be adapted for that purpose.

Another option I see is to compute the time evolution for a long interval in terms of non stationary dynamics, and then take a sort of mean value of the last-time values to construct a steady solution, but it seems to be very unstable numerically, as the 

precision for the ODE solver is not as good as I'd want for that thing...

Sorry for the long message; I'd like to collaborate to the development in case of the routines are not finished.

Thanks in advance

Javier

[jrjoh...@gmail.com]

Hi Javier

Once you have calculated the Bloch-Redfield tensor R, for example using

R, ekets = bloch_redfield_tensor(H, a_ops, [S_w])

then you could use the same steady state solver as for the Lindblad master equation, e.g. 

rhoss = steady(R)

or

rhoss = steadystate(R) 

in the development version of qutip. However, if your system is non-markovian this might not work (if you have an integral equation for the density matrix).

With time dependent systems (external drive) calculating the steadystate also becomes trickier. In that case it is probably best to transform to a frame where the drive is time-independent under a rotating wave approximation (if you can) or perhaps look into calculating the Floquet states for the time-dependent problem (if the drive is periodic). 

Rob

[javierde...@gmail.com]

Hi again; thank you very much for the answer; there is something I don't understand, sorry

in the development version of qutip. However, if your system is non-markovian this might not work (if you have an integral equation for the density matrix).

Once you have computed the correlation functions spectra S given at eq (4)  of the documentation you write the equation as (5), using the expression (6) for the Bloch-Redfield tensor so I think I don't get how you could still have an integral equation...

(Indeed understanding this would be very interesting because observables at long times using the "generation of tensor + time solver" and "generation of tensor + steady(R)" doesn't seem to coincide at all)

Concerning the other point, thank you very much for it; if i put an external driving I'm able to go to a rotating frame and make Rotating Wave Approximation so that constitutes no problem.

Best

Jaier

 

 

 

[Javier Del Pino Gutiérrez]

Sorry I was meaning using "steadystate(R)", no steady(R).  

Even if I compute the stationary density matrix using this, in the case of a flat density of states (so Markovian and Redfield are the same, like in this very well explained example made by yours: https://github.com/qutip/qutip-notebooks/blob/master/examples/example-bloch-redfield.ipynb ) the results of long-time expectation values and steadystate(R) don't agree (indeed, in this particular case steady(liouvilian) != steady(bloch redfield tensor)

Best

Javier 

[Javier Del Pino Gutiérrez]

up!

[jrjoh...@gmail.com]

Hi Javier

Sorry for being late with getting back to you on this. I'm not sure what your problem is, I haven't been able to reproduce it. I've created a notebook with some test cases and it seems to behave as I would expect:

http://nbviewer.ipython.org/github/qutip/qutip-notebooks/blob/master/examples/example-brmesolve-steadystate.ipynb

Have you taken into account that the Bloch-Redfield tensor is written in the eigenbasis? Meaning that after using 'rhoss_eb = steadystate(R)' you have to tranform rhoss_eb from the eigenbasis to the original basis. See the linked notebook for examples of what I mean.

If you still are having problems with this please provide some examples that demonstrates the issues.

Rob

[Javier Del Pino Gutiérrez]

Hi Rob, thank you very much for your response, and thank you for your effort.

I've realized that the problem arises when I change the system operators which I couple to the bath; I'm still using hermitian system operators so the the computation of the Bloch-Redfield tensor using the formula from the documentation should be correct.

I've prepared a short version of the code which give me problems and put as attached.

Thanks in advance; I hope the code is clear.

Best

Javier