Simulate non-hermitian Hamiltonian with QuTip

[Michele Cotrufo]

Hi, 

I was wondering if and up to which extent I can use Qutip to simulate non-Hermitian Hamiltonians. I'm particularly interested in the Hamiltonian described in eq. 4 of this paper https://www.nature.com/articles/natrevmats201648?WT.feed_name=subjects_quantum-mechanics or, similarly, in eqs. (1-4) of this paper http://www.jetp.ac.ru/cgi-bin/e/index/e/44/3/p477?a=list (note that in this second paper they work in a "first quantization language")

I will briefly summarize the problem. Consider the case of two cavities, with frequencies Omega1 and Omega2, interacting with each other with a coupling rate J. The two cavities are also interacting with a continuum of modes of a third system (e.g., a waveguide). We assume that the coupling between the two cavities and the modes of the waveguide does not depend on the frequency of the waveguide mode. However, the coupling between the two cavities and the waveguide can be different (Here are called Gamma1 and Gamma 2). The Hamiltonian is shown here https://pasteboard.co/GEsKrLs.png , where b_\omega is the annihilation operator of the modes in the waveguide. This Hamiltonian is quite nasty because of the presence of the infinite continuum of modes. 

However, after a bit of algebra you can actually show that the degree of freedom of the waveguide can be "eliminated", and the system can be equivalently described by the non-Hermitian Hamiltonian shown here https://pasteboard.co/GEsN1vj.png (I omitted some numerical constant and some 2 pi, just for clarity). As you see, the two cavities now become "lossy", and also their mutual interaction acquires an imaginary part.

Is Qutip able to handle this kind of Hamiltonian? I did a first test by just considering a single-cavity Hamiltonian and adding an imaginary part to the cavity frequency, but there is no decay in the cavity population.

Do you have any suggestion on how I could tackle this problem? 

thanks

[Michele Cotrufo]

(small typos in the last equation, this is the correct one https://pasteboard.co/GEtffh4.png , nothing fundamental is changed)

[Kevin Fischer]

Hi Michele,

QuTiP by default renormalizes the density matrix or state vector after each time step. If you want to remove this, all you need to do is add 

options=qutip.Options(normalize_output=False)

as an argument to the sesolver or mesolver, which will allow for simulation of non-Hermitian Hamiltonians.

Regards,

Kevin

[Michele Cotrufo]

Thanks for the suggestion!! I will let you know the results! :)

Best,

Michele

[Michele Cotrufo]

Hi Kevin, if I use that option in the mesolve function it seems to work! great!

I was wondering if the same could be done for the steadystate solver. It does not seem to accept that option and it generates an error message.

Maybe is there a fundamental reason for which I cannot use non-hermitian Hamiltonians when looking for the steady state?

Thanks,

Michele

[Michele Cotrufo]

Eventually, I'm interested in calculating the spectrum of one of the two cavities under an external pumping. The pumping is also implemented directly in the Hamiltonian. I tried using the first of the two methods shown here http://qutip.org/docs/3.1.0/guide/guide-correlation.html#emission-spectrum , in particular the lines 

tlist = np.linspace(0, 100, 5000)
corr = correlation_ss(H, tlist, c_ops, a.dag(), a)
wlist1, spec1 = spectrum_correlation_fft(tlist, corr)

where I additionally included the option that you suggested before when calling the correlation_ss function.

However, I get an error from correlation_ss saying that I didn't specify any collapse operator, and therefore it cannot calculate the correlation function. But in my case the dissipation is included directly in the Hamiltonian, and therefore I dont use any Lindblad operator.

Any suggestion?

Thanks,

Michele

[Kevin Fischer]

Hi Michele,

In the eq. 4 of the paper "Bound states in the continuum," they write that H_eff governs the dynamics of the system according to a first-order ODE that looks like a Schrodinger equation. My understanding is that the equation they write is for c-numbers (coherent state amplitudes) rather than for state vectors. Despite the visual similarity to the Schrodinger equation, it represents the equations of motions for the annihilation operators of each mode <a_1(t)> and <a_2(t)>. If you wrote down the Heisenberg equations of motion for the operator expectation values with the effective Hamiltonian, that would give you the right equations of motion, but this leaves out some quantum dynamical issues compared to Schrodinger evolution.

Non-Hermitian Hamiltonians have eigenvalues that either decay (to zero) or grow (non-physical) - thus the effective Hamiltonian evolution actually leaks probability at each time step. If you unravel the master equation into terms that go into the effective Hamiltonian compared to those that must be written in super-operator form, you find that the effective Hamiltonian is missing re-normalization of the state at each time-step. So as far as I understand, steady-state calculations can't work with the Schrodinger equation and a non-Hermitian Hamiltonian.

When you say you want to calculate the spectrum, do you mean watching one cavity's occupation vary as you sweep the laser frequency? There are a few ways to do this with QuTiP - the easiest being to incoherently pump the cavity (at a rate slow compared to the cavity decays) and use the correlation function you noted from the documentation. Although I work on cavity-dot systems, much of the ideas are similar, and I discussed a few of these points on cavity spectra in a book chapter. Please see section 2.4 of https://arxiv.org/pdf/1701.03039.pdf if you're interested.

Hope this is helpful.

Regards,

Kevin

[Michele Cotrufo]

Dear Kevin,

sorry for the huge delay of this answer. In the past months I have been trying to study this system with different methods, and I wanted to be sure that I explored all my options before coming back to annoy you :)

The Hamiltonian that I showed in the first post can indeed be handled by QuTip by setting the option that you mentioned. In that Hamiltonian, the presence of the waveguide is described by adding an imaginary part to the frequencies of the resonators and to their mutual interaction rate.

However, I became interested in simulating numerically the full dynamic of the system, i.e. including the waveguide which acts as an excitation channel for the system. Note that I cannot use the Lindblad operator to incoherently pump the system, because I am interested in the correlations that arise between the two cavities due to the fact that they interact with the same waveguide.

Do you think that it would be possible (and/or correct) to simulate the presence of the waveguide by adding a third cavity with a very large linewidth and which couples to the first two cavities? Or do you know of any other trick that could be used for simulating a waveguide in QuTip? I was thinking about describing the waveguide by a chain of coupled cavities, but that would increase dramatically the dimension of the Hilbert space...

(I have spent quite some time addressing the problem with the analytical approach of the single-photon transport in 1D waveguide that has been described in many papers by the group of Shanhui Fan and others. While being quite an elegant method, the calculations become rather complicated and lengthy when going from a 1-excitation sector to a 2-excitation sector, so having some numerical backup could be extremely useful)

Thanks in advance for any help!

p.S. I gave a look to that book chapter - very interesting, thanks! 

[Kevin Fischer]

Hi Michele,

I think I might be able to help, but I'm a bit confused about your problem setup. Perhaps it would be easiest if we setup a video chat sometime next week?

Regards,

Kevin

[Michele Cotrufo]

Great! I sent you an email!

Best,

Michele 

[Vishvendra Singh Poonia]

Hi Kevin,

This method doesn't seem to work for us in the latest version of QuTiP. We are simulating a non-Hermitian Hamiltonian with mesolve() in QuTiP 4.7. Even If we put:

options=qutip.Options(normalize_output=False)

Still the trace of the density matrix is preserved (equal to 1).

On the other hand, if we use a Hermitian Hamiltonian, both diagonal and off-diagonal elements are oscillating without decay. With non-Hermitian Hamiltonian, off-diagonal elements decay but the trace is still preserved. We expect the population to decay as well with non-Hermitian Hamiltonian.

Could you please help us out with this.